# Preconditioning ARPACK eigenvalue solver

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?

• 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. Commented Feb 28, 2017 at 13:20
• @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) Commented Feb 28, 2017 at 14:08
• 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. Commented Feb 28, 2017 at 14:37
• @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$). Commented Feb 28, 2017 at 14:40
• @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... Commented Mar 1, 2017 at 10:58

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.`

(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).