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

Krylov-Schur EigenSolver for (real) non-Hermitian matrices #132

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

Krylov-Schur EigenSolver for (real) non-Hermitian matrices #132

wants to merge 19 commits into from

Conversation

dotnotlock
Copy link

@dotnotlock dotnotlock commented Oct 8, 2021 •
edited
Loading

Dear authors and contributers of Spectra,

I have implemented a Krylov Schur algorithm for non-Hermitian matrices which computes eigenvalues and eigenvectors for the generalized Eigenproblem. I have created this type of solver, because the existing solvers in Spectra doesnt worked for my problems well. I'm dealing with system, in which some eigenvalue gaps can be near 0 (equal eigenvalues). The system matrices are generated in a typical Finite Element procedure. The SymGEigsSolver doesnt handled these system good. Thatswhy i decided to implement another algorithm. The structure of the code is based on SymGEigsBase and SymGEigsSolver.

The code itsself is based on MATLAB's "eigs" function. I have tested the solver with matrix size of A and B up to 300k×300k entries. The solver only needs some seconds to find correctly the first 6 eigenvalues (LargestMagn or for the inverse problem SmallestMagn).
It works stable and very fast with large sparse Eigenproblems. It uses all tweaks which MATLAB also using. For that reason i also replaced the original Arnoldi iterations with the one from MATLAB (class KrylovSchur [Spectra\LinAlg\KrylovSchur.h]). Its a little bit different from the existing Anroldi iterations, but very stable because of a robust reorthogonalization algorithm (see Matlabs "eigs" code).

Currently only real input Matrices A and B are allowed, but the algorithm whould also work on Matrices with complex scalars. The Schur decomposition is also available for complex matrices (Eigen::ComplexSchur), as well as the eigensolver "Eigen::ComplexEigenSolver". From this point of view it should be possible with some minor changes on the current code.

For some reasons i also have implemented the Regular Inverse Algorithm using "Eigen::SparseLU" alogirthm. Sometimes the CG-Solver doesnt converge or is much slower than a classical LU decomposition. Maybe you can also use it with the other solvers. It would be also worth to consider a new enum in GEigsMode for the Regular Inverse LU operator.

Tell me what are you thinking. I know the code is currently not perfect now, but its a good start.
I can also create some more significant test cases.

It would be a pleasure for me and a good reference to contribute something useful to the Spectra project.
Best Regards
David

jschueller and others added 2 commits September 15, 2021 21:18
@jschueller
CMake: include Eigen targets
819cc48 
The targets file links to the Eigen3::Eigen target so we must include here instead of the client side
@dotnotlock
Krylov-Schur EigenSolver for (real) non-Hermitian matrices
1ddd4f3 
Applies Krylov Schur algorithm for non-Hermitian matrices. Computes eigenvalues and eigenvectors for the generalized eigen problem.
JensWehner
JensWehner reviewed Oct 9, 2021
View reviewed changes
include/Spectra/LinAlg/KrylovSchur.h Outdated
Vector m_fac_v; // current v vector (residual)
Scalar m_beta; // ||f||, B-norm of f

bool robust_reorthogonalize(MapConstMat& Vjj, Vector& f, Scalar& fnorm, const Index jj, const Index seed, Vector& wIn)
Copy link
Contributor

Choose a reason for hiding this comment

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

could you add more documentation of what this does?

Copy link
Author

Choose a reason for hiding this comment

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

done in respect to the reference:
G. W. Stewart; Matrix algorithms. Volume 1, Basic decompositions; 1998; page 287, Algorithm 4.1.13

JensWehner
JensWehner reviewed Oct 9, 2021
View reviewed changes
include/Spectra/LinAlg/KrylovSchur.h Outdated
Comment on lines 101 to 118
for (int reorth = 1; reorth <= 5; reorth++)
{
// Check orthogonality
Vector VMf(jj + 1);
VMf.noalias() = Vjj.transpose() * f;
fMf = f.norm();

Scalar ortho_err = VMf.cwiseAbs().maxCoeff();
if (abs(fMf - 1) <= 1e-10 && ortho_err <= 1e-10)
{
stopAlgorithm = false;
break;
}

// Re-orthogonalize
f.noalias() -= Vjj * VMf;
f.normalize();
}
Copy link
Contributor

Choose a reason for hiding this comment

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

can you reuse the code from above?

Copy link
Author

Choose a reason for hiding this comment

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

in my opinion difficult because of the stopAlgorithm part

JensWehner
JensWehner reviewed Oct 9, 2021
View reviewed changes
test/KrylovSchur.cpp
using Matrix = Eigen::MatrixXd;
using Vector = Eigen::VectorXd;
using SpMatrix = Eigen::SparseMatrix<double>;

Copy link
Contributor

Choose a reason for hiding this comment

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

Can we add a test for non-symmetric matrices and also one where the cholesky factorisation is not used?

JensWehner
JensWehner reviewed Oct 9, 2021
View reviewed changes
Copy link
Contributor

@JensWehner JensWehner left a comment

Choose a reason for hiding this comment

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

Really nice contribution. I do not know the algorithm of matlab, just gave some thoughts.

Thanks, really nice.

@dotnotlock
Copy link
Author

Thanks for your comments and hints! i will correct and comment the code in the coming days.

Thank you.
It would be really nice for me to become a contributer to the project.

@JensWehner
Copy link
Contributor

@yixuan is the one with the keys. :)

@yixuan
Copy link
Owner

yixuan commented Oct 11, 2021

Thank you David. I'm sorry that I'm a bit busy recently, so it would be nice if you allow me some time reviewing the code. One quick impression is that much of the code is copied from existing classes, so you may consider if there is some better way reusing the code. Thanks.

@dotnotlock
review correction #1
e52e51d
@dotnotlock
correction 2
1ca85c3 
- derived KrylovSchurSolverBase from SymEigsBase to reuse existing code
- derived KrylovSchur from Arnoldi to reuse existing code
- deleted SparseRegularInverseLU
- extended SparseRegularInverse to make all solvers in Eigen available to use
- KrylovSchur test minimalized to use system matrices out of a regular Finite Element beam model (matrix_A and matrix_B file)
- added more documentation
@dotnotlock
Copy link
Author

dotnotlock commented Oct 13, 2021 •
edited
Loading

Thank you @yixuan and @JensWehner for the hints so far.
i have made some major correction to reuse the existing code:

  • derived KrylovSchurSolverBase from SymEigsBase to reuse existing code
  • derived KrylovSchur from Arnoldi to reuse existing code
  • deleted SparseRegularInverseLU
  • extended SparseRegularInverse to make all solvers in Eigen available to use
  • KrylovSchur test minimalized to use system matrices out of a regular Finite Element beam model (matrix_A and matrix_B file)
  • added more documentation

The extension of the SparseRegularInverse class can be very helpful if one wants to use other Eigen sparsematrix solvers. The helper class "Solver" here solves the polymorphism problem with the Eigen library solvers. They are indeed derived from Eigen::SparseSolverBase but the template parameter makes it difficult to implement any polymorphism in the "SparseRegularInverse" class. So i simply declared all available solvers in the helper class. With the enum SparseRegularInverse ::SolverType the solver type can be switched. This method shouldnt make any problems in respect of memory usage compared to the original SparseRegularInverse class. This extension is also fully compatible with existing syntax and examples/tests. The default solver is still the CG solver.

Best Regards
David

JensWehner
JensWehner reviewed Oct 13, 2021
View reviewed changes
include/Spectra/MatOp/SparseRegularInverse.h
@@ -64,7 +220,7 @@ class SparseRegularInverse
/// `Eigen::Map<Eigen::SparseMatrix<Scalar, ...> >`.
///
template <typename Derived>
SparseRegularInverse(const Eigen::SparseMatrixBase<Derived>& mat) :
SparseRegularInverse(const Eigen::SparseMatrixBase<Derived>& mat, SolverType type = SolverType::ConjugateGradient) :
Copy link
Contributor

Choose a reason for hiding this comment

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

this works, instead we could also just have a template argumement which takes any type which provides a solve interface. Sort of like the preconditioners in Eigen

Copy link
Contributor

Choose a reason for hiding this comment

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

so basically compile time polymorphism vs runtime. @yixuan what do you prefer?

JensWehner
JensWehner reviewed Oct 13, 2021
View reviewed changes
test/EigenIOFile.h Outdated Show resolved Hide resolved
JensWehner
JensWehner reviewed Oct 13, 2021
View reviewed changes
test/EigenIOFile.h Outdated Show resolved Hide resolved
dotnotlock and others added 3 commits October 13, 2021 22:35
@dotnotlock
correction 3
c4b5baa 
some minor review corrections
@yixuan
Merge pull request yixuan#128 from jschueller/patch-2
1690f6a 
CMake: include Eigen targets
@dotnotlock
test file correction
2caf6f7 
uses marketmatrix file import instead of binary sparse matrix import
@dariomangoni
Copy link

Hi @dotnotlock, I'm really really interested in this new branch; I was looking exactly for this kind of solvers (probably because we are working on the same topic/applications). Matlab can solve damped FE systems with no issues at all (cannot tell the performance though) while IR-Arnoldi is really really struggling in retrieving consistent results, even using ARPACK directly.

I see that you are having problems with Eigen::SparseLU; I had a lot of problems myself and I realized they were caused by AVX2 instruction sets (see projectchrono/chrono#290 and https://gitlab.com/libeigen/eigen/-/issues/2126 for reference). For example, the same issues didn't show up in Debug mode. And I think that the same problem affects also Eigen::SparseQR.
Might be the same for you? Please mind that I was on Eigen 3.3.9; I didn't try with 3.4 yet.
That's actually the reason for which I decided to always use either Pardiso MKL or MUMPS by wrapping them in Spectra-compliant classes and they work smoothly.

@JensWehner
Copy link
Contributor

@yixuan could you enable the CI for this PR?
https://docs.github.com/en/actions/managing-workflow-runs/approving-workflow-runs-from-public-forks

@dotnotlock
Update Basic.yml
75566fa
@dotnotlock
Update checkformat.yml
23d50fa
@dotnotlock
Update codecov.yml
87467a9
@dotnotlock
Eigen downward compatibility
01f7dbb 
compatible with Eigen < 3.4
@dotnotlock
Merge branch 'krylovschur' of https://github.com/dotnotlock/spectra i…
9bce82d 
…nto krylovschur
@dotnotlock
test validation
7448bc3 
- SparseRegularInverse::Solver::compute function corrected
- KrylovSchur test file path correction
@dotnotlock
Copy link
Author

finally passed all workflow checks (dotnotlock#1)

@dotnotlock
Copy link
Author

When everything went good with this PR, i can extend the Krylov-Schur Solver for Hermitian matrices.
MATLAB has 2 different implementations for Hermitian and non-Hermitian matrices. But the code just slightly differs, so it would be possible to merge both with the current implementation in Spectra.
Also the support for complex sparse matrices can be considered. This could be a great extension :D

@yixuan
Copy link
Owner

yixuan commented Oct 29, 2021

I hope I can finish this soon. Just a quick question: for Hermitian matrices, Krylov-Schur is equivalent to the implicitly restarted Lanczos algorithm, it that correct? Because Schur decomposition on a Hermitian matrix is basically an eigen decomposition.

@codecov-commenter
Copy link

Codecov Report

Merging #132 (7448bc3) into master (904bbbe) will decrease coverage by 3.1%.
The diff coverage is 66.6%.

Impacted file tree graph

@@           Coverage Diff            @@
##           master    #132     +/-   ##
========================================
- Coverage    92.1%   88.9%   -3.2%     
========================================
  Files          45      48      +3     
  Lines        2140    2436    +296     
========================================
+ Hits         1971    2167    +196     
- Misses        169     269    +100
Impacted Files Coverage Δ
include/Spectra/MatOp/SparseRegularInverse.h 46.1% <39.2%> (-42.8%) ⬇️
include/Spectra/LinAlg/KrylovSchur.h 67.0% <67.0%> (ø)
include/Spectra/KrylovSchurGEigsBase.h 80.8% <80.8%> (ø)
include/Spectra/KrylovSchurGEigsSolver.h 100.0% <100.0%> (ø)
include/Spectra/LinAlg/Arnoldi.h 60.3% <0.0%> (-0.4%) ⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 904bbbe...7448bc3. Read the comment docs.

@JensWehner
Copy link
Contributor

@yixuan Schur vs Eigendecomposition of hermitian matrices, that is also my understanding.

@dotnotlock
Copy link
Author

Yes you are right, the Schur decomposition is not needed for hermitian matrices. Its indeed equal with implicitly restarted Lanczos algorithm

@yixuan
Copy link
Owner

yixuan commented Nov 16, 2021

I just want to leave another comment to say sorry that I'm a bit tied up recently. I put this PR in a high priority, so once I get some free time I will come back here. Thanks for your understanding.

@tasora
Copy link

tasora commented Nov 18, 2021 •
edited
Loading

Dear @dotnotlock et al., I am testing the Krylov Schur solver and it works fine.
Now, I implemented a specialized version for the shift-and-invert problem, because I have the following matrices for constrained structural dynamics:

         A  =  [ -K   -Cq' ]
	       [ -Cq    0  ]
	
	 B  =  [  M     0  ]
	       [  0     0  ]

where K may be singular too (ex. when I need to extract the six rigid body modes at zero frequency in a free-free structure), so using shift-and-invert with a small sigma shift works perfectly in all cases.

What I did is to inherit a specialized Krylov-Schur class, inspired to SymGEigsShiftSolver , that is:

template <typename OpType, typename BOpType>
class KrylovSchurGEigsShiftInvert :
    public KrylovSchurGEigsBase<SymGEigsShiftInvertOp<OpType, BOpType>, BOpType>
{
private:
    using Scalar = typename OpType::Scalar;
    using Index = Eigen::Index;
    using Array = Eigen::Array<Scalar, Eigen::Dynamic, 1>;

    using ModeMatOp = SymGEigsShiftInvertOp<OpType, BOpType>;
    using Base = KrylovSchurGEigsBase<ModeMatOp, BOpType>;
    using Base::m_nev;
    using Base::m_ritz_val;

    const Scalar m_sigma;

    // Set shift and forward
    static ModeMatOp set_shift_and_move(ModeMatOp&& op, const Scalar& sigma)
    {
        op.set_shift(sigma);
        return std::move(op);
    }

    // First transform back the Ritz values, and then sort
    void sort_ritzpair(SortRule sort_rule) override
    {
        // The eigenvalues we get from the iteration is nu = 1 / (lambda - sigma)
        // So the eigenvalues of the original problem is lambda = 1 / nu + sigma
        m_ritz_val.head(m_nev).array() = Scalar(1) / m_ritz_val.head(m_nev).array() + m_sigma;
        Base::sort_ritzpair(sort_rule);
    }

public:
    KrylovSchurGEigsShiftInvert(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma) :
        Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv),
        m_sigma(sigma)
    {}
};

then I can use it as:

    using OpType = SymShiftInvert<double, Eigen::Sparse, Eigen::Sparse>;
    using BOpType = SparseSymMatProd<double>;
    OpType op(A, B);
    BOpType Bop(B);

	// The Krylov-Schur solver, using the shift and invert mode:
	KrylovSchurGEigsShiftInvert<OpType, BOpType> eigen_solver(op, Bop, n_modes, m, sigma);

It works perfectly but the issue here is that the set_shift_and_move is never called from your KrylovSchurGEigsBase, so after the compute(), I must postprocess the eigenvalues doing the shift by myself:

for (int i = 0; i < eigen_values.rows() ; ++i)
		eigen_values(i) = (1.0 / eigen_values(i)) + sigma;

Are you planning to support the "shift and invert" scheme later in a new push to develop, or did I miss something?
Please tell me if I am reinventing the wheel...

@tasora
Copy link

tasora commented Nov 18, 2021

A second comment: the Krylov-Schur in Matlab works also with non-symmetric matrices, for complex eigenpairs.
However, I was not able to make KrylovSchurGEigsBase working to extract complex eigenvalues, as I expect when passing the following matrices:

         A  =  [  0     I     0 ]
	       [ -K    -R  -Cq' ]
	       [ -Cq    0     0 ]
	
	 B  =  [  I     0     0 ]
	       [  0     M     0 ]
	       [  0     0     0 ]

Note that I need to solve this with a shift-and-invert, because alternative approaches requiring only the inversion will fail some special cases.
I tried to use the solver that I presented in my previous comment, but it does not return any complex eigenvalue, as I expect from a comparison with Matlab.
Is the complex eigenvalues case (with shift and invert) not yet supported? Or planned later?

Thanks and compliments for the great library.

@dotnotlock
Copy link
Author

@tasora thanks for the comment. Nice to hear that it works well so far.
I will take a look on the Shift and Invert scheme implementation in the coming days. Should be no problem.
For the return of complex eigenvalues i will try to add a ComplexKrylovSchurSolver or something like that. Currently only real input matrices are supported, because of the RealSchur function. But yes it is planned to get implemented. I will take your example as a test case

         A  =  [  0     I     0 ]
	       [ -K    -R  -Cq' ]
	       [ -Cq    0     0 ]
	
	 B  =  [  I     0     0 ]
	       [  0     M     0 ]
	       [  0     0     0 ]

Thank you :)

@JensWehner
Copy link
Contributor

@tasora can you just submit a second PR, so we can have a look?

@tasora
Copy link

tasora commented Nov 21, 2021 via email

@dotnotlock
Copy link
Author

dotnotlock commented Nov 21, 2021 via email •
edited
Loading

@yixuan
Copy link
Owner

yixuan commented Mar 28, 2022 •
edited
Loading

Today I spent some time reading the whole implementation, and got a quick question to ask: from the title of this PR and our previous discussions it seems that you designed the solver for non-Hermitian matrices. However, looking at the class member functions, I only see real-valued eigenvalues and eigenvectors. Did I miss anything?

@dotnotlock
Copy link
Author

Yes exactly, unfortunatly i had implemented and tested only for real values and RealSchur function. But it should be no problem to expand the algorithm for complex values.

@yixuan
Copy link
Owner

yixuan commented Mar 29, 2022

Oh, I mean, the A and B matrices can be both real-valued, but if A is not symmetric, then the eigenvalues and eigenvectors are supposed to be complex-valued, right? But I only see eigenvalues() with a return type of real vector.

@tasora
Copy link

tasora commented Mar 29, 2022 via email

@dotnotlock
Copy link
Author

I think i could expand this to complex values. Eigen offers the ComplexEigenSolver which return complex values. Also ComplexSchur is available in Eigen. So from this point its possible. I just have to find some time for this, then i will update the PR and offer some Test example for the complex case.

@talregev
Copy link

Any news? This PR can merge as it is?

@talregev talregev mentioned this pull request Jan 20, 2023
Closed
@dsl2501 dsl2501 mentioned this pull request Aug 19, 2023
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers

@JensWehner JensWehner JensWehner left review comments

Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.
8 participants
@dotnotlock @JensWehner @yixuan @dariomangoni @codecov-commenter @tasora @talregev @jschueller