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

[dotnotlock]

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

[dotnotlock]

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]

@yixuan is the one with the keys. :)

[yixuan]

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]

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

[dariomangoni]

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 https://github.com/projectchrono/chrono/issues/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.

[dotnotlock]

@dariomangoni thanks for your comment. I also had problems with the existing Arnoldi procedures in Spectra. Thatswhy i decided to implement this Krylov-Schur algorithm. It works very stable for symmetric FE-models which can have eigengaps near zero (equal eigenvalues). In my current test FE-models the LU-Solver works very fast and stable. Didnt realized any problem for this. But i have realized that iterative solvers, like CG, sometimes not converging or are very slow even for simple FE-Models. Thatswhy I have integrated all available solver from Eigen into Spectra::SparseRegularInverse.

Short test explaination: The model matrices used in the test file (tests/KrylovSchur.cpp) are equivalent to the following: B-Matrix = Stiffness Matrix (constraint) A-Matrix = Mass Matrix The current test solve the inverse eigenproblem: M-1/eigenvalue*K=0 Normally a constraint stiffness matrix is always invertable, while the mass matrix can be singular. The inverse eigenproblem bypasses the inversion of the mass matrix and inverts the stiffness matrix instead. I have build everything with Eigen 3.4.0.

Are you refering to a damped modal analysis when you say "damped FE systems"? ANSYS for example uses a QR-decomposition for damped eigenvalue and -vector extraction.

[dariomangoni]

I see. I misread that you had problems with SparseLU. Now I read better.

Actually for your setup I have quite consistent results already with SymGEigsShiftSolver; the Shift did a lot of difference to me. But it's only for symmetric matrices and it's not (always) my case.

What I'm trying to solve is a problem with potentially unsymmetric stiffness and damping matrices problem, together with constraints (and possibly without computing the null space of the constraints). Like this (taken from Yang - A Direct Eigenanalysis of Multibody System in Equilibrium)

Matlab solves this quite robustly, while Arpack is having a lot of problems though. I'll surely give a try to your solver!

[yixuan]

The code did not pass the clang-format check, so I further edited this PR in the develop branch. I feel it is already in good shape, and I will try to merge it soon, possibly on the weekend.

[yixuan]

Just noticed that I can edit the PR, so we can directly work on your branch. I just pushed the changes.

[yixuan]

I found a compatibility issue here. The template Eigen::Vector<> in KrylovSchurGEigsBase.h was introduced by Eigen 3.4.0, so the code does not compile with Eigen 3.x.y. Could you make a simple change to Eigen::Matrix<>? Thanks.

[dotnotlock]

@yixuan thanks for the check run! yes, it didnt passed the clang format and there is also an compilation error with this line in KrylovSchurSolverBase.h: reorder.indices() = Eigen::Map<Eigen::Vector<Index, Eigen::Dynamic>>(ind.data(), ind.size()).cast(); i was using Eigen 3.4.0 and a local VS compiler. i will try to recheck again by myself using your workflow checks and give another update to pass all checkups. Best Regards :)

[JensWehner]

@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]

finally passed all workflow checks (https://github.com/dotnotlock/spectra/pull/1)

[dotnotlock]

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]

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]

Codecov Report

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

@@           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%)	:arrow_down:
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%)	:arrow_down:

---

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]

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

[dotnotlock]

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

[yixuan]

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]

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]

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]

@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]

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

[tasora]

Dear David:

@tasorahttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Ftasora&data=04%7C01%7Calessandro.tasora%40unipr.it%7C23abca8baf7e4b42ee8a08d9ac007381%7Cbb064bc5b7a841ecbabed7beb3faeb1c%7C0%7C0%7C637729940213152892%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=h6bt4Y2TrC84CtCoEcQb8nmoVZF0c6jIoohHjgvyelk%3D&reserved=0 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.

Thank you!

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.

In fact, the current Krylov Schur solver is inherited as

class KrylovSchurGEigsBase : public SymEigsBase<OpType, BOpType>

so I suppose that your future version for general matrices could be inherited as

class ComplexKrylovSchurGEigsBase : public GenEigsBase<OpType, BOpType>

(by the way, maybe it is better to rename the current KrylovSchurGEigsBase as SymKrylovSchurGEigsBase and the new one as GenKrylovSchurGEigsBase just to be coherent with the pre-existing name schemes of all other solvers, that always start with Gen.... Sym... ?)

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 ]

Yes this is a classical eigenvalue problem in structural dynamics and aeroelasticity. It is one of the various ways to solve the quadratic eigenvalue problem (s^2 M + s R + K)x = 0. Differently from the "classical" eigenvalue problem with just M and K, here the eigenvalue is a complex number (its modulus is the undamped natural frequency), and one also gets the damping factor (as a ratio between the real part of the eigenvalue and its modulus; if negative it means instability).

A small added complication is that here I also introduce a jacobian Cq that expresses the presence of constraints on the possible displacements, as when using joints.

If you want, I can send you an example of M,K,R, Cq matrices for benchmarking and comparison with Matlab "eigs".

Thanks

Alessandro Tasora

Thank you :)

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fyixuan%2Fspectra%2Fpull%2F132%23issuecomment-974616024&data=04%7C01%7Calessandro.tasora%40unipr.it%7C23abca8baf7e4b42ee8a08d9ac007381%7Cbb064bc5b7a841ecbabed7beb3faeb1c%7C0%7C0%7C637729940213152892%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=AECA2TV1g4f8WOfs3hUs0gakFKHP%2FBmnq2hb8xmGvjc%3D&reserved=0, or unsubscribehttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FABHTOBAQPJ6TQJDEBAE4NT3UM5MODANCNFSM5FTPX2EQ&data=04%7C01%7Calessandro.tasora%40unipr.it%7C23abca8baf7e4b42ee8a08d9ac007381%7Cbb064bc5b7a841ecbabed7beb3faeb1c%7C0%7C0%7C637729940213162849%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=1EZG2gqNzXzdYWdIYj%2F%2F3a52a5pvvtp3vPylO4nRyNw%3D&reserved=0. Triage notifications on the go with GitHub Mobile for iOShttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fapps.apple.com%2Fapp%2Fapple-store%2Fid1477376905%3Fct%3Dnotification-email%26mt%3D8%26pt%3D524675&data=04%7C01%7Calessandro.tasora%40unipr.it%7C23abca8baf7e4b42ee8a08d9ac007381%7Cbb064bc5b7a841ecbabed7beb3faeb1c%7C0%7C0%7C637729940213162849%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=PW8aUX4TtGGZtPhnZfL0rDCYSLX2Lv%2F%2BPVldZ9thfNo%3D&reserved=0 or Androidhttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fplay.google.com%2Fstore%2Fapps%2Fdetails%3Fid%3Dcom.github.android%26referrer%3Dutm_campaign%253Dnotification-email%2526utm_medium%253Demail%2526utm_source%253Dgithub&data=04%7C01%7Calessandro.tasora%40unipr.it%7C23abca8baf7e4b42ee8a08d9ac007381%7Cbb064bc5b7a841ecbabed7beb3faeb1c%7C0%7C0%7C637729940213172813%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=RrsPH48RbrrmkqozRxTS8vhYue2TE1S1Pu3Ka1PGmU0%3D&reserved=0.

Firma il tuo 5x1000 all’Università di Parma, aiutaci a essere sempre più accoglienti e inclusivi verso le nostre studentesse e i nostri studenti - Indica 00308780345 nella tua denuncia dei redditi.

[dotnotlock]

Hello,

yes you can send a sample for the systemmatrices. You can send me this also in Matlab MAT file format, so i can create some input for c++ test file. Thank you :)

Best Regards David

[yixuan]

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]

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]

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]

Dear Yixuan Qiu

i am interested in the case of the complex-value eigenvalues/eigenvectors in fact. My problems involve sparse real-valued but /not symmetric/ A and B,  with shift-and-invert, and that's the case that I cannot solve efficiently with the current version of Spectra.

I tried to modify the experimental Krylov Schur solver some weeks ago but I gave up because it took too many hacks from my side and it did not work as expected, so I preferred to wait until someone from the official Spectra team can do an "official" implementation of Krylov Schur for non-symmetric matrices.

I hope that the Krylov Schur solver could support this case because there are not so many c++ alternatives for doing this...

best regards

Alessandro Tasora

On 29/03/2022 02:55, Yixuan Qiu wrote:

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.

— Reply to this email directly, view it on GitHub https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fyixuan%2Fspectra%2Fpull%2F132%23issuecomment-1081291676&data=04%7C01%7Calessandro.tasora%40unipr.it%7C9158009108e54574d61608da111ecf1a%7Cbb064bc5b7a841ecbabed7beb3faeb1c%7C0%7C0%7C637841121279273440%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=ymoTEuSkXFVYs8penvKx98H93x6KiA4GfyyFURWVH8A%3D&reserved=0, or unsubscribe https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FABHTOBFRYKLGAKBNY3WXDELVCJIHXANCNFSM5FTPX2EQ&data=04%7C01%7Calessandro.tasora%40unipr.it%7C9158009108e54574d61608da111ecf1a%7Cbb064bc5b7a841ecbabed7beb3faeb1c%7C0%7C0%7C637841121279273440%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=otft4CLx4dboDCVU9d8LFi5DzZMjXnhqsJOnx4Llmok%3D&reserved=0. You are receiving this because you were mentioned.Message ID: @.***>

[dotnotlock]

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]

Any news? This PR can merge as it is?