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$ Sep 14 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$ Sep 15 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$ Sep 15 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$ Sep 16 at 3:51
  • 1
    $\begingroup$ DSYEV is a good option, you will need to change the routine depending the working precision $\endgroup$ Sep 16 at 7:30

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, privacy policy and cookie policy

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