Solution (source code)

= Solution

Since $A$ is diagonalizable over the complex numbers, write $x^{(0)}=\sum_i a_iv_i$. Provided iteration does not encounter zero, its direction is that of
$$
A^kx^{(0)}=\sum_i a_i\lambda_i^kv_i.
$$
With one strictly dominant <eigenvalue> and $a_n\ne0$, all other contributions are smaller by factors $(|\lambda_i|/|\lambda_n|)^k$. The dominant <eigenvalue> is real, since nonreal <eigenvalues> of a real matrix occur in conjugate pairs of equal modulus. The normalized <power method> therefore converges to a signed unit <eigenvector> if $\lambda_n>0$, and alternates between opposite unit <eigenvectors> if $\lambda_n<0$. A quotient such as $(x^{(k)})^TAx^{(k)}$ approaches $\lambda_n$ in either case. If $a_n=0$, replace it by the largest active modulus; if $Ax^{(k)}=0$, the prescribed iteration is undefined.

With two dominant <eigenvalues> of equal modulus and both components active, the iterates approach their real invariant plane rather than one direction. Distinct real equal-modulus <eigenvalues> are $\rho,-\rho$ and produce a limiting two-cycle. A conjugate pair $\rho e^{\pm i\theta}$ produces a rotation in suitable coordinates on that plane: the normalized directions are periodic when the angle is rational relative to $2\pi$ and otherwise generally keep rotating without convergence. If the two <eigenvalues> are equal and real, the initial projection is itself one <eigenvector> in their repeated eigenspace, and converges or alternates as in the single-value case; one orbit cannot recover every <independent> <eigenvector> of that eigenspace.

For distinct dominant values, retain unnormalized consecutive vectors $y,Ay,A^2y$ in the nearly invariant plane and solve
$$
A^2y=sAy-py.
$$
The roots of $\lambda^2-s\lambda+p$ estimate the two <eigenvalues>. Corresponding <eigenvectors> are $Ay-\lambda_2y$ and $Ay-\lambda_1y$, since applying $A$ and using the recurrence multiplies them by $\lambda_1$ and $\lambda_2$. This is <dominant eigenpair extraction from a two-step Krylov recurrence>. Equivalent extraction uses the $2\times2$ matrix of $A$ on the computed plane. It needs two <independent> iterates and nonzero projections onto both distinct <eigenvectors>.