Reworked upper Hessenberg decomposition
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 itHessenberg? I suppose havingDecompositionin 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_intowhensubdiagonal == 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 forsubdiagonal == 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 ofleft_multiply_into, as well asfirst_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.