2
$\begingroup$

I am working on a generalized eigenvalue problem of the form

$$ \boldsymbol{A}\cdot\boldsymbol{x}=\lambda\boldsymbol{B}\cdot\boldsymbol{x} $$

where $\boldsymbol{B}$ is not symmetric positive. Therefore I recast the problem to $$ \boldsymbol{B}^{-1}\cdot(\boldsymbol{A}\cdot\boldsymbol{x})=\lambda\boldsymbol{x}\ . $$ I am lucky: the matrix $\boldsymbol{B}$ is block diagonal. Therefore I can compute the inverse by an LU-factorization of the diagonal blocks.

This eigenvalue problem is then solved with the ARPACK solver. ARPACK is used because of its feature to extract specific parts of the spectrum like eigenvalues with the largest real part. The method works well for small problems.

However, if I increase the problem size, the performance of the ARPACK algorithmen in terms of iterations required increases. The condition number of $\boldsymbol{B}^{-1}\cdot\boldsymbol{A}$ increases heavily up to the order $10^6$.

Is there a way of preconditioning the ARPACK algorithmen to accelerate the convergence?

Improve this question
$\endgroup$
10
  • 1
    $\begingroup$ Do I understand correctly that you are looking only for the largest eigenvalues? How many iteration vectors (ncv) are you using? You might try increasing this number to see how that affects the computation cost. I suggest you provide the details on exactly how you are using arpack. $\endgroup$
    – Bill Greene
    Feb 28, 2017 at 13:20
  • $\begingroup$ @BillGreene I am not always looking for the same part of the spectrum. Sometimes I need the say 12 eigenvalues with largest real part, sometimes the 12 eigenvalues of smallest magnitude. I am not specifying ncv at all. But I will try. I am using arpack through scipy. My call looks like this `l, v = scipy.sparse.linalg(LinOp, k = 12, which = "SM", tol = 1e-12, maxiter = 200) $\endgroup$
    – sebastian_g
    Feb 28, 2017 at 14:08
  • 1
    $\begingroup$ That complicates things considerably. I don't know the scipy interface to arpack. Does your LinOp simply multiply a vector with a matrix? If so, the convergence of the smallest eigenvalues can be very slow. As I suggested previously, you can try making ncv much larger. The only efficient way to calculate the smallest eigenvalues is to essentially "invert" A. $\endgroup$
    – Bill Greene
    Feb 28, 2017 at 14:37
  • 1
    $\begingroup$ @sebastian_g How are you implementing your LinOp applied to a vector $v$? It it something along the lines of (1) $ z \leftarrow A v$ (2) Solve $M w = z$ (for $w$, using the pre-computed factorization of $M$). $\endgroup$
    – GoHokies
    Feb 28, 2017 at 14:40
  • 1
    $\begingroup$ @GoHokies Thank you for your comment. Is the shift-invert approach also suitable with $\sigma=0$? I think that is the best guess for my problem. Do you have experience regarding Krylov-Schur methods? If I understood the SLEPc documentation right, these methods also allow to search for certain parts of the spectrum. I could incorporate PETSc and SLEPc using their python wrappers. But that would be a lot of effort... $\endgroup$
    – sebastian_g
    Mar 1, 2017 at 10:58

1 Answer 1

Reset to default
3
$\begingroup$

I've summarized the comments thread of your original question into an answer.

Here are a few things that you can try:

  1. Increase the number of your Arnoldi vectors (NCV) generated at each iteration. Here's what the ARPACK documentation for DSAUPD says about this:

    At present there is no a-priori analysis to guide the selection of NCV relative to NEV (the number of eigenvalues you're looking for). The only formal requirement is that NCV > NEV. However, it is recommended that NCV >= 2*NEV.

    If many problems of the same type are to be solved, one should experiment with increasing NCV while keeping NEV fixed for a given test problem. This will usually decrease the required number of OP*x operations but it also increases the work and storage required to maintain the orthogonal basis vectors. The optimal "cross-over" with respect to CPU time is problem dependent and must be determined empirically.

  2. Revisit your implementation of the linear operator. In your case, it should read something like:

    (a) $ z \leftarrow A v$

    (b) Solve $M w = z$ for $w$, using the pre-computed factorization of $M$.

  3. If you're looking for the smallest eigenvalues - or eigenvalues close to a given real value, it is advisable to use the shift-and-invert mode of ARPACK. The shift $\sigma$ can be zero (see the example they give here) or any other real/complex value. Your linear operator becomes $(A - \sigma M)^{-1} M $ (for $M$ symmetric and indefinite).

  4. If SciPy-ARPACK works well enough for your problem, then PETSc/SLEPc are probably overkill.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.