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.

Sign up for GitHub

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

Jump to bottom

Incorporate GSoC 2015 Eigensolver into Develop branch #68

Open
wants to merge 21 commits into
base: develop
Choose a base branch
from
Open

Incorporate GSoC 2015 Eigensolver into Develop branch #68

wants to merge 21 commits into from

Conversation

thomasyang1207
Copy link

After reviewing the code and tests produced from Rajaditya Mukherjee's GSoC 2015 project, I wish to work towards merging this with the main development branch.

The current develop branch does not contain methods for producing matrix inverses. In other words, there is no solver without external dependencies on existing routines like LAPACK, etc. Having one will allow many advanced algorithms involving matrices and tensors to be implemented.

I'm sure there is a reason why it wasn't merged, but at the moment, including it in the development branch would make it much more accessible to other developers

rajaditya-m and others added 17 commits June 3, 2015 20:33
@rajaditya-m
Hessenberg
cbdfd8f 
The Hessenberg decomposition module for square matrix is completed...
@rajaditya-m
Francis QR Step
060f437
@rajaditya-m
FRancis Double Step QR
6c8137a 
This finished the recovery of eigen values
@rajaditya-m
Transform Accumulations
b3dbc13 
This is needed for the computation of Eigen Vectors
@rajaditya-m
Eigen Solver Wrapper
6ac85ee 
eigen_solver.hpp which can compute real eigenvalyes of real matrices..
@rajaditya-m
Givens Rotation and Modified Schur Decomposition
62872a8 
Implemented Givens rotation and modified schur decomposition for better
convergence in complex case
@rajaditya-m
Complex Eigen Values
ae04672 
Complex Eigen Values are supported
@rajaditya-m
eigenvalue new routine
3d48bb2
@rajaditya-m
Complex Eigenvector
58e0bb8
@rajaditya-m
Small Numerical Fix
2e73352
@rajaditya-m
Matrix Balancing
803d24d
@rajaditya-m
Doxygen Comments
147766e 
Documentations
@rajaditya-m
SingleTest and sample
24531ff
@rajaditya-m
Test Cases
929ec70
@rajaditya-m
HTML Docs
ac89993
@rajaditya-m
Matlab's Magic Shift Function
1cbafef
@thomasyang1207
Merge GSOC 2015 eigensolver into Develop branch
44543f1
@bassoy
Copy link
Collaborator

bassoy commented Apr 21, 2019 via email

@yimyom
Copy link
Collaborator

yimyom commented Apr 23, 2019

Really, the only reason was lack of time. I've kind of lost contact with the student too. I think he finished his PhD and was busy with his new work life, which I understand (every morning :-D ).
We can work together on adding comparison with R, and Octave with Cem for example.

Should I merge the branch now or we wait until you've got more benchmarks?

@NAThompson
Copy link

I would recommend performance comparison with Eigen. I used Eigen in this PR:

https://github.com/boostorg/math/blob/b14975fa7f8eb49aec89aa8f51fb1615afe126b5/example/daubechies_scaling_integer_grid.cpp

But of course I would've liked to remove the external dependency.

@thomasyang1207
Copy link
Author

At the very least, should we include a reference (let's say to a paper in arxiv) for each of the algorithms?

I also think more benchmarks are important since the code is almost 4 years old

@bassoy
Copy link
Collaborator

bassoy commented Apr 23, 2019 via email

@yimyom
Copy link
Collaborator

yimyom commented Apr 27, 2019 via email

@thomasyang1207
Copy link
Author

I added 5 more 10x10 test matrices with numbers at different orders of magnitude. I also changed the test suite so that it now covers the double and float data types. For some reason, my changes aren't showing up when I click on my new commit - any reason why this is?

Right now, I'm still a bit uncertain whether the eigensolver can handle complex matrices... AFAIK Eigen has a separate eigensolver specifically for complex and hermitian matrices.

I just realized that the eigensolver test suite was never added to the jamfile, so it was never run during build. For the eigensolver test suite and future test suites, do we add it manually to the Jamfile? Until then, we cannot merge anything

For the new test cases I added, I also computed the same result in Eigen. Inside benchmarks I could potentially just copy the "Eigen" directory found inside the Eigen package, and then make the comparison there with the ublas eigensolver.

There are also some other minor formatting changes I could do regarding some of the function signatures - adding const specifiers, etc. Some of the comments that rajaditya left also give us some unfinished tasks

@bassoy
Copy link
Collaborator

bassoy commented Apr 28, 2019 via email

@thomasyang1207
Copy link
Author

Currently, all the tests for the eigensolver (both old ones and new ones) pass. I fixed some minor bugs in the test code and the code for the eigensolver as well.

Right now the test cases are not written using the full extent of the Boost Test Suite. Could I get some guidance on that?

The next step for me is adding benchmarks

@bassoy bassoy self-assigned this May 12, 2019
@yimyom
Copy link
Collaborator

yimyom commented May 12, 2019 via email

@penguian
Copy link

penguian commented May 12, 2019 via email

@bassoy
Copy link
Collaborator

bassoy commented May 12, 2019 via email •
edited
Loading

@thomasyang1207
Copy link
Author

I'll get going on reproducing this issue with ublas's test suite.

The line about the complex eigenvalues is in the file eigen_solver.hpp, throughout the eigen_solver class definition. I can fix that as well

@penguian I'd love to learn more about the specifics of implementing an eigensolver, since I'm actually quite new to all this. I have added you as a collaborator to my branch

penguian
penguian suggested changes May 19, 2019
View reviewed changes
Copy link

@penguian penguian left a comment

Choose a reason for hiding this comment

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

The proposed code is unfinished and contains unused functions which need to be incorporated and tested. See also these pull requests and branches of which two are in progress.

For a description of the LAPACK implementation of the multi-shift QR algorithm, as used in (e.g.) DGEES, see Lapack 3.1 xHSEQR: Tuning and Implementation Notes. This uses the double shift algorithm for small matrices.

See also the equivalent Eigen code.

include/boost/numeric/ublas/schur_decomposition.hpp
* @tparam[in] m matrix type (like matrix<double>) - input hessenberg form output - schur form
*/
template<class M>
void schur_decomposition(M &h) {
Copy link

Choose a reason for hiding this comment

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

This function is incomplete as a practical implementation of the Francis Double Shift QR Algorithm. [1, p. 173]: "In practice, the QR iteration with the Francis shift can fail to converge ... So the practical algorithm in use for decades had an "exceptional shift" every 10 shifts if convergence had not occurred. ... another exceptional shift was recently added to the algorithm."

[1] James Demmel, Applied Numerical Linear Algebra, SIAM, 1997.

include/boost/numeric/ublas/schur_decomposition.hpp
}

template<class M>
void find_small_diag_entry(M &h, typename M::size_type end, typename M::value_type l1_norm, typename M::size_type &small_index)
Copy link

Choose a reason for hiding this comment

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

Unused function. See the equivalent Eigen method RealSchur<MatrixType>::findSmallSubdiagEntry which is used in RealSchur<MatrixType>::computeFromHessenberg. This function should be used in schur_decomposition above, but is not.

include/boost/numeric/ublas/schur_decomposition.hpp
}

template<class M>
void row_split(M &h, M &qv, typename M::size_type end, typename M::value_type exceptional_shift_sum)
Copy link

Choose a reason for hiding this comment

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

Unused function. See the equivalent Eigen method RealSchur<MatrixType>::splitOffTwoRows which is used in RealSchur<MatrixType>::computeFromHessenberg. This function should be used in schur_decomposition above, but is not.

include/boost/numeric/ublas/schur_decomposition.hpp
}

template<class M>
void infer_shifts(M &h, typename M::size_type end, typename M::size_type iter_nos, typename M::value_type exceptional_shift_sum, vector<typename M::value_type> &shift_vector)
Copy link

Choose a reason for hiding this comment

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

Unused function. See the equivalent Eigen method RealSchur<MatrixType>::computeShift which is used in RealSchur<MatrixType>::computeFromHessenberg. This function should be used in schur_decomposition above, but is not.

This is the function that contains the exceptional shifts that I mentioned above at line 36. See below.

include/boost/numeric/ublas/schur_decomposition.hpp

//Original Shift
if (iter_nos == size_type(10))
{
Copy link

Choose a reason for hiding this comment

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

This is "Wilkinson's original ad hoc shift", as per the equivalent Eigen method RealSchur::computeShift.

include/boost/numeric/ublas/eigen_solver.hpp
template <class M>
class eigen_solver {
private:
M matrix;
Copy link

Choose a reason for hiding this comment

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

Most of the source files retain tab characters, which automatically translate to 8 spaces when displayed here. These tabs should be replaced by 4 spaces.

Copy link

Choose a reason for hiding this comment

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

See branch clean-indentation (in progress).

test/test_eigensolver.cpp
MATRIX Lambda_Real(3, 3);
MATRIX Lambda_Complex(3, 3);

for (int i = 0; i < 3; i++){
Copy link

Choose a reason for hiding this comment

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

This tests only the first 3 eigenvalues, instead of all eigenvalues. See pull request Test all eigs.

include/boost/numeric/ublas/eigen_solver.hpp
bool has_eigenvalues;
bool has_eigenvectors;
M eigenvalues_real;
M eigenvalues_complex;
Copy link

Choose a reason for hiding this comment

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

The term complex as used here and elsewhere is misleading and should be replaced by imag or `imaginary. See pull request API: Change misleading instances of complex to imag or imaginary.

include/boost/numeric/ublas/eigen_solver.hpp
}
}

for (size_type j = n; j--!=size_type(0);)
Copy link

Choose a reason for hiding this comment

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

It would be better to have

for (size_type j = n - 1; j != 0; j--)

if that is what is meant.

test/test_eigensolver.cpp
"[10,10]((4, 2, 3, 4, 3, 6, 7, 8, 9, 10),(0, 2, 3, 3, 5, 0, 7, 8, 9, 10),(3, 2, 0, 4.2, 5, 6.1, 0, 8, 9, 0),(1, 2, 3, 4, 5, 6, 0, 8, 9, 1),(3, 2, 3, 4.5, 5, 6, 0, 8, 9, 0),(5, 2, 0, 4, 5, 6, 1, 8, 9, 1),(0, 2, 3, 0, 5, 6, 7, 0, 9, 4),(1, 13, 0, 4, 5, 6, 2, 8, 9.3, 0),(100, 2, 3, 0, 0, 6, 7, 8, 9, 10),(0, 2, 3, 4, 5, 6, 0, 8, 9, 0),)\0",
"[10,10]((0, 2e+10, 3e+10, 4e+10, 3e+10, 6e+10, 7e+10, 8e+10, 9e+10, 1e+11),(1e+10, 2e+10, 3e+10, 3e+10, 5e+10, 0, 7e+10, 8e+10, 9e+10, 1e+11),(3e+10, 2e+10, 0, 4.2e+10, 5e+10, 6.1e+10, 0, 8e+10, 9e+10, 0),(1e+10, 2e+10, 3e+10, 4e+10, 5e+10, 6e+10, 0, 8e+10, 9e+10, 1e+10),(1e+10, 2e+10, 3e+10, 4.5e+10, 5e+10, 6e+10, 7e+10, 8e+10, 9e+10, 0),(5e+10, 6e+10, 0, 4e+10, 5e+10, 6e+10, 7e+10, 8e+10, 9e+10, 1e+10),(0, 2e+10, 3e+10, 0, 5e+10, 6e+10, 7e+10, 1e+10, 9e+10, 4e+10),(1e+10, 2e+10, 0, 4e+10, 5e+10, 6e+10, 7e+10, 8e+10, 9.3e+10, 0),(0, 2e+10, 3e+10, 4e+10, 0, 6e+10, 7e+10, 8e+10, 9e+10, 1e+11),(0, 2e+10, 3e+10, 4e+10, 5e+10, 6e+10, 7e+10, 8e+10, 9e+10, 0))\0"
};

Copy link

Choose a reason for hiding this comment

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

More tests are needed, including tests which cause the current implementation to fail to converge.

Copy link

Choose a reason for hiding this comment

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

See branch more-eigensolver-tests (in progess).

penguian
penguian reviewed May 19, 2019
View reviewed changes
include/boost/numeric/ublas/eigen_solver.hpp
size_type n = real_schur_form.size1();
size_type i = size_type(0);

eigenvalues_real = zero_matrix<value_type>(n,n);
Copy link

Choose a reason for hiding this comment

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

The real and imaginary components of the eigenvalues should be returned as vectors, especially when eigenvectors are not required. Compare with Eigen.

penguian
penguian reviewed May 19, 2019
View reviewed changes
test/test_eigensolver.cpp
}

if (!_fail_counter) {
std::cout << "\nEigen Solver regression suite passed.\n";
Copy link

Choose a reason for hiding this comment

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

Usually in C++, end of line is indicated by std::endl. The use of "\n" here is C usage.

penguian
penguian reviewed May 19, 2019
View reviewed changes
test/test_eigensolver.cpp
}
}

assertTrue("Real portion of verifier is zero:", compare_on_tolerance(Lambda_Real, Zero_Matrix));
Copy link

Choose a reason for hiding this comment

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

From line 121 to here: why not just compare (e.g.) norm _2(Lambda)/norm_2(prod(V,D)) with a small tolerance?

penguian
penguian reviewed May 19, 2019
View reviewed changes
test/test_eigensolver.cpp
M(i, j) = std::complex<TYPE>(A(i, j), 0.0);
}
}
Lambda = prod(M, V) - prod(V, D);
Copy link

Choose a reason for hiding this comment

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

The name Lambda here is misleading, because in this context, the letter lambda is used to denote a generic eigenvalue of a matrix. Maybe use Verifier, since that is what the message below calls it.

penguian
penguian reviewed May 26, 2019
View reviewed changes
include/boost/numeric/ublas/eigen_solver.hpp

value_type t = (std::max)((std::abs)(T(i, k - size_type(1))), (std::abs)(T(i, k)));
if ((std::numeric_limits<value_type>::epsilon() * t)*t > value_type(1)) {
for (size_type j = i; j < k + i - size_type(1); j++){
Copy link

Choose a reason for hiding this comment

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

The bound of k + i - size_type(1) here is incorrect and will cause j to index past the end of the matrix. Compare with EigenSolver.h lines 596-598.
See thomasyang1207 pull request #5 for the fix.

@penguian
Copy link

Hi @thomasyang1207 . I have commented on this merge request and have made suggested changes of my own. See the branches and merge requests at https://github.com/thomasyang1207/ublas
I have done about as much as I can without diving in to the eigensolver and Schur code to properly complete it. I now need to spend some extra time at work and will have little time to work on this pull request for at least 4 weeks.

In the meantime, you can study the Eigenvalues module of Eigen and read

[1] Demmel, Applied Numerical Linear Algebra, SIAM, 1997. Chapter 4 Nonsymmetric Eigenvalue Problems. (Section 4.4.8 covers QR iteration with implicit shifts, but it would be best to read from the start of the chapter.)

[2] Trefethen and Bau. Numerical Linear Algebra, SIAM 1997, Part V Eigenvalues. (Lecture 29 covers the QR algorithm with shifts, but it would be best to read from Lecture 24.)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers

@penguian penguian penguian requested changes

At least 2 approving reviews are required to merge this pull request.

Assignees

@bassoy bassoy

Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.
6 participants
@thomasyang1207 @bassoy @yimyom @NAThompson @penguian @rajaditya-m