Example for solving asymmetric generalized eigenvalues with shift-invert method

[pan3rock]

Could you please give an example about it in Eigen or Armadillo format? It's hard for me to understand how to add shifts_f and access eigenvalues. Thanks a lot.

[krivenko]

Hi there,

FYI, I am now working on a proper set of examples and documentation in general. It should be merged into master by the end of the month.

Meanwhile, you could have a look at the last two sections of test/eigen/asymmetric.cpp.

I can also come up with a simple example if you provide some extra information.

• Do you know what computation mode you are interested in? There are two versions of the shift-invert mode in the non-symmetric case, ShiftAndInvertReal and ShiftAndInvertImag. The second one is usable only when complex shift \sigma has a non-zero imaginary part (see section 3.7 of ARPACK Users' Guide).

• Are you sure you want to use shifts_f? shifts_f has nothing to do with the complex shift \sigma used in the spectral shift-invert transformation. Instead, it provides a way to customize the implicit restart procedure when something different from ARPACK's default "Exact Shift Strategy" is needed (section 4.4 of the Guide).

[pan3rock]

I'm interested in ShiftAndInvertReal. Maybe I misunderstood shifts_f, and I saw that shift-invert mode will determine a few minimum eigenvalues if sigma assigned in some books. I'm major in geophysics, and my recent work involves the quadratic eigenvalue problem that can be converted to the generalized eigenvalue problem. Therefore, I search for a solver for it and I'm stuck in it.

I used Eigen:: GeneralizedEigenSolver but the performance is not well. The matrix is sparse and I have known the eigenvalues scope, maybe there is an efficient solver.

[krivenko]

Thanks for the explanation.

First of all, you do not need shifts_f as it really is an expert level feature unrelated to the spectral shift \sigma. Just omit it from all calls to arpack_worker::operator().

Now the problem is that your matrix on the RHS does not seem to be symmetric, and ARPACK wants it to be symmetric and positive semi-definite. Would it be mathematically possible and numerically feasible to invert it and multiply both sides of the eigenproblem by the inverse? You would get a standard eigenproblem supported by ARPACK then.

[krivenko]

Looks like one has to do a transformation to the standard eigenproblem in the most general case...

[krivenko]

Alright, let me think a little bit. I guess there is a fairly straightforward way to implement this using tools from Eigen and the standard eigenproblem mode of arpack_worker.

[krivenko]

@pan3rock You should get a general idea of what to do from the following example. I have successfully compiled it but never ran as testing it would require some valid input matrices.

#include <iostream>

#include <ezarpack/arpack_worker.hpp>
#include <ezarpack/storages/eigen.hpp>

 #include<Eigen/SparseLU>

using namespace ezarpack;
using namespace Eigen;

int main() {

  // Solve a generalized eigenproblem A x = \lambda M x with real matrices A and M.
  // Neither A nor M is assumed to be symmetric positive semi-definite.
  // In this case one cannot rely on arpack_worker's ShiftAndInvertReal mode.
  // The spectral shift-and-invert transformation has to be done by the user.

  // Input matrices
  const int N = 1000;
  SparseMatrix<double> A;
  SparseMatrix<double> M;

  //
  // Fill A and M as needed ...
  //

  // Spectral shift for the shift-and-invert transformation.
  // It should be chosen to be close to the eigenvalues of interest.
  double sigma = 1.0;

  // The goal here is to reduce the original generalized eigenproblem to a
  // standard one, O x = \mu x, where O = (A - \sigma M)^{-1} M, and
  // \lambda = \sigma + 1 / \mu.

  // Use Eigen's sparse LU solver to decompose A - sigma * M.
  // Here, I assume that the LU decomposition exists, which is not always the
  // case. If it does not, one has to resort to computing (A - \sigma M)^{-1} M
  // at once, which is much slower.
  SparseLU<SparseMatrix<double>> lu;
  lu.compute(A - sigma * M);

  // Asymmetric arpack_worker
  using worker_t = arpack_worker<ezarpack::Asymmetric, eigen_storage>;
  using params_t = worker_t::params_t;
  using vector_view_t = worker_t::vector_view_t;
  using vector_const_view_t = worker_t::vector_const_view_t;

  // Implementation of linear operator O
  auto o = [&](vector_const_view_t x, vector_view_t y) {
    // Act with (A - \sigma M)^{-1} M on x in two steps.

    // I. Act on x with M
    y = M * x;

    // II. Solve (A -\sigma M) z = M x for z and put the result into y.
    y = lu.solve(y).eval();
  };

  worker_t worker(N);

  // Compute 10 eigenvalues \mu with the largest magnitudes.
  params_t params(10, params_t::LargestMagnitude, params_t::Ritz);
  worker(o, params);

  auto const& mu = worker.eigenvalues();

  // Back-transformation of the eigenvalues
  auto lambda = mu.cwiseInverse().eval();
  lambda.array() += sigma;

  std::cout << lambda << std::endl;
}

[pan3rock]

I wrote a simple shift-invert function that only one eigenvalue is obtained. I compared the version you offered, my simple version and Eigen::GeneralizedEigenSolver. The matrices used are dense matrix. It indicates that the solution you offered is the best. The codes is here. Then I will try it in my problem with sparse matrix. In addition, I replace LU decomposition with householderQr in codes above.

Thank you very very much.

20x20, 4 smallest eigenvalues | Function Name | Total | Be Called | Percentage | | Total Run Time | 1.69116 ms | 1 times | 100 % | | shift-invert | 0.0949707 ms | 1 times | 5.615 % | | ezsolver | 0.665039 ms | 1 times | 39.324 % | | GeneralizedEigenSolver | 0.864014 ms | 1 times | 51.089 % |

100x100, 4 | Function Name | Total | Be Called | Percentage | | Total Run Time | 31.3301 ms | 1 times | 100 % | | shift-invert | 3.42603 ms | 1 times | 10.935 % | | ezsolver | 3.59009 ms | 1 times | 11.458 % | | GeneralizedEigenSolver | 23.7961 ms | 1 times | 75.953 % |

1000x1000, 4 | Function Name | Total | Be Called | Percentage | | Total Run Time | 25213.1 ms | 1 times | 100 % | | shift-invert | 352.012 ms | 1 times | 1.396 % | | ezsolver | 158.247 ms | 1 times | 0.627 % | | GeneralizedEigenSolver | 24698.6 ms | 1 times | 97.959 % |

[krivenko]

@pan3rock You are very welcome!

The matrices used are dense matrix. It indicates that the solution you offered is the best.

ARPACK should normally be even more beneficial when you switch to the sparse matrices.

Btw, I checked your benchmark code and saw the following lines in ezsolver.cc.

  auto o = [&](vector_const_view_t x, vector_view_t y) {
    // y = matB_ * x;
    // y = lu.solve(y).eval();
    auto vecR = matB_ * x;
    y = lhh.solve(vecR).eval();
  };

I believe that moving declaration of vecR outside of the lambda function o will eliminate a memory allocation and potentially improve performance.

https://github.com/pan3rock/shift-invert/blob/3138fc86eba9be73b46e1888dfc7e265e71708a5/src/ezsolver.cc#L24-L29

[krivenko]

@pan3rock Please, let me know if you are satisfied so that this issue can be closed.

[pan3rock]

Applied to my problem, I met a strange situation:

when the eigenvalues are: 1.37200e+00 1.73669e+00 1.99408e+00 2.00003e+00 2.00010e+00 2.00023e+00 2.00043e+00 2.00075e+00 2.00115e+00 2.00158e+00 2.00224e+00 2.00358e+00 2.00571e+00 2.00878e+00 2.01116e+00 2.01544e+00 2.02443e+00 2.03765e+00 2.04701e+00 2.06769e+00 2.13337e+00 2.21554e+00 2.33741e+00 2.43624e+00 2.53733e+00 2.76864e+00 2.86346e+00 3.00000e+00 3.00000e+00 3.00000e+00 3.00000e+00

If I get the first 3 eigenvalues, the elapsed time is reasonable. However, if I get more than 3 eigenvalues, it is stuck and much slower than Eigen::GeneralizedEigenSolver. Maybe It is relate to the fact that the remaining eigenvalues are close to each other.

[krivenko]

Well, the run time of the Arnoldi iterations is very sensitive to the distribution of the eigenvalues. Also, it can grow very quickly with the amount of requested eigenvalues -- so quickly that pretty much any other method will become faster at some point. In general, ARPACK is designed to compute a few eigenvalues of large sparse matrices.

If the number of the eigenvalues you need is indeed small compared to the matrix size, you might want to experiment with different values of \sigma.

Btw, your eigenvalues are real. Are you still using the asymmetric version of arpack_worker or the symmetric one?

[pan3rock]

As you said, using different sigma with very small number of eigenvalues may be a good solution for my problem.

My eigenvalues are complex, but only eigenvalues that ev.real() > 0 and ev.imag() == 0 are accepted.

Thank you again.