I am considering implementing the accelerated inverse power (AIP) method with Rayleigh quotient to speed up eigendecomposition of a real square symmetric matrix. Halton (1996) gives an example algorithm which is listed below

Step 1. $\quad \mathbf{M}= \mathbf{H} - \mu \mathbf{I} $

Step 2. $\quad \text{Solve } \mathbf{M}\mathbf{y}= \mathbf{z}$

Step 3. $\quad \mathbf{z}= (\mathbf{y}^\top \mathbf{y})^{-1/2} \mathbf{y} $

Step 4. $\quad \rho = z_{k_{h}} / y_{k_{h}} $

Step 5. $\quad \mu = \mu + \rho $

Repeat Steps 1-5 until the delta between the updated eigenvector $\mathbf{z}$ and previous $\mathbf{z}$ are below some value.

Question is, why does he have subscripts for $k$ and $h$ on the ratio elements in the Rayleigh quotient? Which vector elements of $\mathbf{z}$ and $\mathbf{y}$ is he referring to?

I have also seen a matrix form of the Rayleigh quotient used in AIP which looked similar to a conjugate gradient step involving a numerator and denominator with vectors sandwiched around a matrix. What would that be, as an example?

Last, I also need all eigenvalues and their eigenvectors, and AIP is only geared toward finding the closest eigenvalue to $\mu$. Since my matrix $\mathbf{H}$ is a $p \times p$ correlation matrix, I know the sum of eigenvalues is equal to the $p$-dimensions of $\mathbf{H}$. Thus, I believe that if I first set $\mu=p$, then I will obtain the major eigenvalue/eigenvector. After $\lambda_1$ and $\mathbf{e}_1$ are found, dimension both $\mathbf{z}$ and $\mathbf{y}$ to be $p \times 2$, and still set $\mu=p$, but don't update $\lambda_1$ inside the $2 \times 2$ $\boldsymbol{\Lambda}$ or $\mathbf{e}_1$ inside the $p \times 2$ $\mathbf{E}$. The results for $\lambda_2$ and $\mathbf{e}_2$ will be inside $\boldsymbol{\Lambda}$ and $\mathbf{E}$.

Has anyone developed steps for finding all eigenvalue/eigenvector pairs iteratively using AIP?

  • 1
    $\begingroup$ do you want to implement the method to use in some project or simply to understand the mathematics better? $\endgroup$ Commented Sep 14, 2021 at 7:33
  • $\begingroup$ Yes, implement. However, it is looking like any variant of AIP won't be able to outperform QL and QR iteration for determining all pairs at once, which I already use (QL). But they're both $6n^3 + O(n^2)$. Dhillon's MMRR is $O(n^2)$, but I am not sure if it's implemented in LAPACK yet. $\endgroup$
    – user39981
    Commented Sep 15, 2021 at 16:06
  • 1
    $\begingroup$ In that case, just use LAPACK. In fact, the QR algorithm can be seen as a more complicated variant of the AIP method. And the reduction step that precedes MMRR is still $O(n^3)$ $\endgroup$ Commented Sep 15, 2021 at 19:19
  • $\begingroup$ Thanks! What's the specific LAPACK package for QR algorithm for a real square symmetric matrix [which is always the correlation matrix having no pathologies, and is always positive-semidefinite (could have a few zero eigenvalues)]? I ask because there's dozens of similar LAPACK algorithms - I know MATLAB studies a matrix, then decides which LAPACK algorithm to use. $\endgroup$
    – user39981
    Commented Sep 16, 2021 at 3:51
  • 1
    $\begingroup$ DSYEV is a good option, you will need to change the routine depending the working precision $\endgroup$ Commented Sep 16, 2021 at 7:30

1 Answer 1


Per comments from @Thijs Steel, this has been answered via suggestion to use the LAPACK DSYEV algorithm.


Your Answer

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