rulinalg icon indicating copy to clipboard operation
rulinalg copied to clipboard

Reworked upper Hessenberg decomposition

Open Andlon opened this issue 8 years ago • 0 comments

This PR continues the trend of reworking our decompositions one-by-one. Because the Hessenberg decomposition is very similar to the QR decomposition, we are able to reuse a large amount of functionality. In order to do this, we must extend the concept of HouseholderComposition to also allow Householder transformations that operate on the subdiagonal of the matrix rather than the main diagonal (which correspond to QR).

I recommend to look through the commit messages for an overview of how this PR is structured and what changes it makes to existing infrastructure.

Some discussion points

  • I've named the decomposition HessenbergDecomposition. Would it be sufficient to call it Hessenberg? I suppose having Decomposition in the name makes it explicitly clear that we are talking about the decomposition. It's strictly speaking an upper Hessenberg decomposition, but I've never even heard of anyone using a lower Hessenberg decomposition, so I think this is OK - the documentation makes it clear that we are talking about an upper Hessenberg decomposition.
  • I haven't actually written tests for HouseholderComposition::left_multiply_into when subdiagonal == 1 (corresponding to Hessenberg decomposition). These tests are rather painful to write due to the need for generating test data. The current test only tests for subdiagonal == 0 (corresponding to QR). However, I think this is OK because the Hessenberg tests verify the correctness of the Hessenberg decomposition, which implicitly depends on the correctness of left_multiply_into, as well as first_k_columns.

Performance

Benchmarks show quite considerably increased performance compared to the existing hessenberg routines:

running 8 tests
test linalg::hessenberg::hessenberg_decomposition_decompose_100x100             ... bench:     687,857 ns/iter (+/- 7,512)
test linalg::hessenberg::hessenberg_decomposition_decompose_200x200             ... bench:   5,333,774 ns/iter (+/- 46,163)
test linalg::hessenberg::hessenberg_decomposition_decompose_unpack_100x100      ... bench:     944,946 ns/iter (+/- 6,644)
test linalg::hessenberg::hessenberg_decomposition_decompose_unpack_200x200      ... bench:   7,405,509 ns/iter (+/- 89,003)
test linalg::hessenberg::upper_hess_decomp_100x100                              ... bench:   4,446,958 ns/iter (+/- 71,237)
test linalg::hessenberg::upper_hess_decomp_200x200                              ... bench:  36,832,089 ns/iter (+/- 137,874)
test linalg::hessenberg::upper_hessenberg_100x100                               ... bench:   3,085,405 ns/iter (+/- 45,021)
test linalg::hessenberg::upper_hessenberg_200x200                               ... bench:  25,728,909 ns/iter (+/- 85,632)

test result: ok. 0 passed; 0 failed; 0 ignored; 8 measured

In the above, hessenberg_decomposition_* corresponds to the new functionality. The new decompose_NxN is comparable in functionality to the old upper_hessenberg function, and decompose_unpack_NxN is comparable to the old upper_hess_decomp function.

Andlon avatar May 07 '17 10:05 Andlon