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?