LOBPCG not works well

[Meisenheimer]

Hello yixuan,

Thank you for developing this wonderful package. I'm working on a large-scale generalized eigenvalue problem, where the inverse matrix or decomposition may be computationally prohibitive. Thus, I want to use the contrib/LOBPCGSolver.h.

However, I found that the example yields incorrect results (as shown below), where the eigenvalues are exactly what I want, and the residuals are also reasonable; however, I have no idea what the eigenvectors are. Did I make any mistakes? Could you please provide some suggestions to help me figure it out?

P.S. The generalized eigenvalue problem is given by the spectral element method with the domain of freedom (or the dimension of matrices) more than $10^5$, where the matrices $A$ and $B$ are symmetric sparse matrices and $B$ is guaranteed to be positive definite. The required eigenpairs are those with the smallest algebraic value. Please let me know if there is a way to apply the SymGEigsSolver or the SymGEigsShiftSolver.

The code is

#include <iostream>
#include <Spectra/contrib/LOBPCGSolver.h>

using Matrix = Eigen::Matrix<long double, Eigen::Dynamic, Eigen::Dynamic>;
using Vector = Eigen::Matrix<long double, Eigen::Dynamic, 1>;

int main()
{
    // random A
    Matrix a;
    a = (Matrix::Random(10, 10).array() > 0.6).cast<long double>() * Matrix::Random(10, 10).array() * 5;
    a = Matrix((a).triangularView<Eigen::Lower>()) + Matrix((a).triangularView<Eigen::Lower>()).transpose();
    for (int i = 0; i < 10; i++)
        a(i, i) = i + 0.5;
    std::cout << a << "\n";
    Eigen::SparseMatrix<long double> A(a.sparseView());
    // random X
    Eigen::Matrix<long double, 10, 2> x;
    x = Matrix::Random(10, 2).array();
    Eigen::SparseMatrix<long double> X(x.sparseView());
    // solve Ax = lambda*x
    Spectra::LOBPCGSolver<long double> solver(A, X);
    solver.compute(10, 1e-4); // 10 iterations, L2_tolerance = 1e-4*N
    std::cout << "info\n"
              << solver.info() << std::endl;
    std::cout << "eigenvalues\n"
              << solver.eigenvalues() << std::endl;
    std::cout << "eigenvectors\n"
              << solver.eigenvectors() << std::endl;
    std::cout << "residuals\n"
              << solver.residuals() << std::endl;
    return 0;
}

and the result gives

      0.5  -3.06696         0  -1.49709         0  -3.25892         0         0         0         0
 -3.06696       1.5         0         0         0         0         0         0   2.83319         0
        0         0       2.5         0         0   2.79656   4.96796         0         0         0
 -1.49709         0         0       3.5         0         0         0         0 -0.420392         0
        0         0         0         0       4.5  0.171056         0         0         0         0
 -3.25892         0   2.79656         0  0.171056       5.5         0         0         0         0
        0         0   4.96796         0         0         0       6.5         0   4.49339         0
        0         0         0         0         0         0         0       7.5  -4.72533         0
        0   2.83319         0 -0.420392         0         0   4.49339  -4.72533       8.5         0
        0         0         0         0         0         0         0         0         0       9.5
info
0
eigenvalues
-3.87057
-1.58281
eigenvectors
          -1 -1.42913e-06
-1.15977e-07    -0.999996
-1.35075e-05  0.000345232
 1.69391e-06 -5.50784e-05
residuals
-4.37446e-05  0.000267876
-1.01832e-05 -0.000291125
 0.000170494 -0.000435531
-5.20398e-05  0.000400166
 0.000192459  5.48295e-05
 0.000130473 -0.000370708
 8.85456e-05 -0.000270534
-4.85901e-05 -0.000103941
-0.000115008  0.000165789
 0.000157712   0.00029169