Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Reworked upper Hessenberg decomposition #179

Open
wants to merge 10 commits into
base: master
Choose a base branch
from

Conversation

Andlon
Copy link
Collaborator

@Andlon Andlon commented May 7, 2017

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 added 10 commits May 2, 2017 20:17
Necessary for testing the validity of the output
from Hessenberg decomposition.
We don't need the buffer to be a specific size,
it's sufficient that it's large enough to hold
the intermediate data. This simplifies usage
a little, since we don't need to truncate or
resize the buffer when it is reused for different
instances of HouseholderReflection.
The composition of Householder matrices for
QR decomposition and Upper Hessenberg decomposition
is very similar. We can exploit this to reuse virtually
all code we have written for HouseholderComposition in the
QR decomposition to implement the Hessenberg decomposition.
We do this by adding a "subdiagonal" parameter to
HouseholderComposition which controls the
subdiagonal to pull Householder reflectors from.

for element in y.iter_mut() {
*element = T::zero();
}
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This zeroing isn't actually necessary! Should remove

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

1 participant