4
$\begingroup$

I'm trying to solve a generalized eigenvalue problem

$$Ax = \lambda Bx, \quad A = A^\top > 0,\; B = B^\top > 0$$

with $\lambda \approx \sigma$ using Rayleigh Quotient Iteration (RQI) (RQI is applied to $B^{-1/2}AB^{-1/2}$, but the $B^{-1/2}$ cancels from the final formulas): $$ (A - \mu_k B) \tilde x_{k+1} = B x_k\\ x_{k+1} = \frac{\tilde x_{k+1}}{||\tilde x_{k+1}||}\\ \mu_k = \frac{x_m^\top A x_m}{x_m^\top B x_m}\\ \mu_0 = \sigma, \quad x_0 = \text{random vector} $$ I've decided to use some iterative method to solve $$ (A - \mu_k B) \tilde x_{k+1} = B x_k $$ for $\tilde x_{k+1}$, since I do not want to store $A$ and $B$ explicitly. I only wish to code computing $Ax$ and $Bx$ products.

The matrix is symmetric, but not positive definite, otherwise I would have used conjugate gradients method. Which method would you suggest? Are there any pitfalls when using Krylov subspace methods with RQI?

$\endgroup$
  • 1
    $\begingroup$ @nicoguaro $A > 0$ is a standard notation to denote a positive definite matrix en.wikipedia.org/wiki/… . $\endgroup$ – uranix Aug 18 '15 at 6:07
  • 1
    $\begingroup$ $A \geq 0$ is psd, $A > 0$ is pd. Yes, both matrices are positive definite $\endgroup$ – uranix Aug 18 '15 at 6:12
4
+50
$\begingroup$

If your matrices are large, why not use a library like ARPACK? The shift-and-invert mode of ARPACK will help you calculate the eigenvalues close to $\sigma$.

There are interfaces to ARPACK for most high-level programming languages used for numerical computations (Fortran, C, Python, MATLAB, etc.).

For example, a quick tutorial for Python can be found here (reasonably small, dense matrices), or here (large, sparse matrices, available only in operator form).

Later edit: I read your post more carefully and noticed you are looking for a way to solve the shift-invert linear system... so why not try MINRES or SYMMLQ? Both algorithms are based on a Lanczos (3-term Arnoldi) iteration specifically designed for symmetric (but not necessarily positive definite) matrices. SYMMLQ should perform better if $A$ is non-singular.

$\endgroup$
  • $\begingroup$ Since I'm trying to solve near-to-singular system $(A - \mu B)$ when $\mu$ is close to an eigenvalue, I'd better stick with MINRES instead of SYMMLQ. But I have troubles with finding some canonical pseudocode for MINRES, $\endgroup$ – uranix Aug 20 '15 at 20:58
  • $\begingroup$ Look here. It is a rather concise MATLAB implementation of MINRES with lots of helpful comments (and a "tutorial" PDF). $\endgroup$ – GoHokies Aug 21 '15 at 7:31
  • $\begingroup$ That implementation is mostly unreadable, already tried. I found one good here $\endgroup$ – uranix Aug 21 '15 at 15:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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