SymGEigsShiftSolver available soon?

[ysonline]

For sparse symmetric 10Kx10K A and B (B is diagonal), generalized problem is solved for the 15 smallest eigenvalues via

• SymGEigsSolver in 20 secs
• converting A v = lambda B v into the standard problem B^{-1} A v = lambda v (B diagonal), and using GenEigsSolver (not SymEigsSolver 'cos B^{-1} not symmetric anymore), this takes 7 secs

I've a feeling that it'll run faster via SymGEigsShiftSolver but I couldn't see it in the current release spectra-0.8.1.tar.gz. So, am I right that SymGEigsShiftSolver (with sigma = -1e-6) will be faster/better (for the small eigenvalues that I need)? And when will it be released so I can use it (saw test/SymGEigsShift.cpp but can't adapt it as SymGEigsShiftSolver not in the current release)?

[yixuan]

SymGEigsShiftSolver will be introduced in the 1.0.0 release, but the API is different from pre-0.9.0 versions. You can find the example code of SymGEigsShiftSolver here and the migration guide for the upgrade.

[ysonline]

Can I use SymGEigsShiftSolver now without waiting 1.0.0? I saw SymGEigsShiftSolver.cpp but this is just testing by include'ing SymGEigsShiftSolver.h. I can't find that file to include in my project. Can you help me with that please?

[yixuan]

SymGEigsShiftSolver was newly added not long ago, and it made use of the refactored code. You can just download the master branch and follow the guide to update your code. The 1.0.0 release will be pretty much like what the library is right now.

[ysonline]

Yeah I did it. Properly set all the paths etc. However, it doesn't compile; here is the first of many errors coming from SymGEigsShiftInvertOp.h:

line31: using Scalar = typename OpType::Scalar; //error: 'Scalar' : symbol cannot be used in a using-declaration

Is it a C++11 issue? I'm using the old C++ in my Visual Studio 2010. Any idea how I can get it run there, or is this a different problem?

[yixuan]

Yeah, you need a modern compiler such as GCC and Clang, or Visual Studio 2015, according to this document.

[yixuan]

You can use this self-contained example to test the compiler:

#include <Eigen/Core>
#include <iostream>

#include <Spectra/SymGEigsShiftSolver.h>
#include <Spectra/MatOp/SymShiftInvert.h>
#include <Spectra/MatOp/DenseSymMatProd.h>

using namespace Spectra;

// Generate data for testing
void gen_dense_data(int n, Eigen::MatrixXd& A, Eigen::MatrixXd& B)
{
    // Eigen solver only uses the lower triangle of A,
    // so we don't need to make A symmetric here
    A = Eigen::MatrixXd::Random(n, n);
    B = A.transpose() * A;
    // To make sure B is positive definite
    B.diagonal() += Eigen::VectorXd::Random(n).cwiseAbs();
}

int main()
{
    std::srand(123);

    // Eigen solver only uses the lower triangle
    Eigen::MatrixXd A, B;
    gen_dense_data(100, A, B);
    int k = 10;
    int m = 20;
    double sigma = 1.2345;

    using OpType = SymShiftInvert<double, Eigen::Dense, Eigen::Dense>;
    using BOpType = DenseSymMatProd<double>;

    OpType op(A, B);
    BOpType Bop(B);
    SymGEigsShiftSolver<OpType, BOpType, GEigsMode::ShiftInvert> eigs(op, Bop, k, m, sigma);
    eigs.init();
    eigs.compute(SortRule::LargestMagn, 100);
    std::cout << eigs.eigenvalues() << std::endl;

    return 0;
}

[ysonline]

Makes sense but there's a problem: I tried to use clang++ -std=c++11 which is good for these new C++11 Spectra code but it fails to compile my whole project written in the old C++, e.g., it errors about Mesh::computeHeatKernelSignature() syntax and many others. I expect some back-compatibility but apparently it's not the case. So, do you recommend a solution, or do you have another SymGEigsShiftSolver compatible with the old C++? That's all I need.

[yixuan]

C++11 should be mostly back-compatible, and will eventually be the standard for future code. If the code does not compile with Clang then very likely there is something not following the standard. I suggest looking at those errors more carefully to see what happens.

[yixuan]

Also, is it simply because you changed the compiler?

[ysonline]

Yes, simply because I changed the compiler. My code that is including your old Spectra release (without SymGEigsShiftSolver) compiles perfectly with: g++ -O3 *.cpp -o output -w

When I include your new Spectra (C++11), I had to use (on the same machine): clang++ -std=c++11 -stdlib=libc++ -Weverything main.cpp -o output -w -O3
And it gives this error referring to my old C++ code (no error from Spectra-side):
Undefined symbols for architecture x86_64
"Mesh::computeHetKernelSignature()", referenced from _main in main-f0daa5.o
.. and many other Mesh::* functions are listed (not printing them here).

Do you have any comment on this compiler error? Thanks.

[yixuan]

It looks like a linking problem, see the discussion here.

Some things you can check:

1. Did you completely clean all .o files before compiling?
2. Is it really necessary to use -stdlib=libc++?

[ysonline]

1. Yes rm *.o can't find any files so no .o files around.
2. Removed -stdlib=libc++ and still the same error.

Thanks for the link; it is very related. As suggested there, I've replaced clang++ with g++, CC, clang, and gcc but still get the same error. I'm on macOS BigSur. Do you know another compiler or any other trick?

[yixuan]

To better debug, I suggest doing the following:

1. Compile the single-file example above and see if the error occurs. If no, then the compiler is good, and there must be problems with other files in the project.
2. Try to create a minimal example to reproduce the linking error (reduce as much code dependency as possible), and even better, paste the code here.
3. Report as much information as possible, for example the version of the compilers, detailed error logs, etc.

[ysonline]

OK will do that. I also want to ask an algebra question that still concerns SymGEigsShiftSolver:

Eigenvectors of a real symmetric matrix, like Laplacian, are already orthonormal: standard eigenvalue problem works: L v = lambda v.
It's advised to make these eigenvectors orthonormal w.r.t. a matrix, namely area matrix A, hence I end up with the generalized eigenvalue problem: L v = lambda A v, which can be converted to the standard version (A^{-1} L v = lambda v) but I want to stick with the generalized version and use SymGEigsShiftSolver.

My question is, what does it mean for eigenvectors to be orthonormal w.r.t. a matrix? What's the benefit of that? Any geometric intuition is highly appreciated.

[yixuan]

My understanding is the following: x and y are orthonormal w.r.t. matrix B if x'By=0, x'Bx=1, y'By=1. It is like redefining the inner product from <x, y>=x'y to <x, y>=x'By.

[Spammed]

what does it mean for eigenvectors to be orthonormal w.r.t. a matrix?

Here's how I see it (Attention: Mathematicians please look away):
To bring matrices on diagonal form means to change its coordinate system in which the system describe so skillfully that the degrees of freedom described individually in it become independent from each other.
Independency can be understood as 'has no common parts' and can be expressed by the condition x'y = 0.
This is exactly the trick of the eigenvector decomposition: Transform the problem into a representation in which all degrees of freedom run independent/separated from each other.
The eigenvector matrix E is exactly the transformation matrix we are looking for.
In the inverse this means that E' B E is a diagonal matrix, thus e_i' B e_j = 0 for i<>j.

In short: It is so that it is useful. How lucky we are! :-)

[stigersh]

Hey,
Thank you so much for developing this tool!
I'm trying to extract a null space of a matrix for geometry processing application, by finding the smallest eigenvalues and their eigenvectors.

I tried using SymGEigsShiftSolver on a positive semidefinite matrix A (Av=Bv) with a B which is diagonal and positive, with sigma offset -1e-6 to make it invertible. Unfortunately it didn't converge.
Any advice about extracting the null space of a matrix? If indeed finding the smallest eigenvalues is the way to go, what would you recommend? eigs of matlab should have worked by extracting it, but I want to stay in cpp.

Thanks a lot!

[yixuan]

@stigersh Can you provide an example matrix to reproduce the problem?

[stigersh]

Sorry didn't see the response.. Thank you for replying!
Apparently my matrix didn't have a null space, and when I took a matrix which has a null space it worked.
Can you please tell me where am I wrong - I expected to get the eigenvectors which belong to the smallest eigenvalues, even if they are not close to 0, I didn't think that the system will diverge in this case.