The power method

Here is a matrix, a starting vector, and the first three rounds of one repeated operation, then the eighth.

!Eight rounds of the power method on $\left(\begin{array}{cc}2& 0\\ 0& 1\end{array}\right)$, with the length divided out at every step. The iterates slide along the unit circle toward $(1,0)$, the dominant eigenvector: the angle falls from $45^{\circ }$ at $t=0$ to $0.22^{\circ }$ at $t=8$. Status: re-simulated.

The table skips rounds four to seven because each one repeats the same arithmetic. The vector is sliding toward $(1,0)$. Nobody solved an equation to make that happen. Each row comes from multiplying the previous vector by the matrix $A=\left(\begin{array}{cc}2& 0\\ 0& 1\end{array}\right)$ and then dividing by its length. The length of $(x,y)$ is $\sqrt{x^2+y^2}$, and it is thrown away on purpose, because only the direction is wanted. !The distance to the dominant direction shrinks by $|{\lambda }_2/{\lambda }_1|$ per step. Reaching one part in a thousand costs about ten steps at a ratio of 0.5 and 135 at 0.95, so a narrow gap between the two largest eigenvalues is expensive. Status: re-simulated.

That procedure is the power method. Start with any nonzero vector, multiply by the matrix, rescale to length one, and repeat. The method finds the dominant eigenvalue and its eigenvector, using nothing but matrix times vector [S:28:136].

Each step costs one matrix-vector product, so the method scales to matrices with millions of rows. It never forms the characteristic polynomial, and it never needs the matrix to be symmetric. When one direction is enough, no step more expensive than a matrix-vector product is needed.

Why does it work? Sort the eigenvalues ${\lambda }_1,{\lambda }_2,\dots$ by absolute value so that $|{\lambda }_1|>|{\lambda }_2|\ge \cdots$. The first inequality is strict: that is what makes every ratio inside the bracket below have modulus below one. Write the start in the eigenbasis, as $v=c_1{v}_1+c_2{v}_2+\cdots$. After $t$ multiplications every coefficient has been multiplied by its own eigenvalue, so

The bracket tends to $c_1{v}_1$, because every ratio inside it has absolute value below one. The overall factor ${\lambda }_1^t$ does not change the direction, and the rescaling removes it anyway. So the direction converges to $v_1$, and the error shrinks by the factor $|{\lambda }_2/{\lambda }_1|$ per step [S:28:135-136]. That factor is the convergence rate.

Worked example. Take $A=\left(\begin{array}{cc}2& 0\\ 0& 1\end{array}\right)$ and the start $(1,1)$. The eigenvectors are $(1,0)$ with eigenvalue $2$ and $(0,1)$ with eigenvalue $1$. After $t$ raw multiplications the vector is $(2^t,1)$, whose length is $\sqrt{4^t+1}$. The normalised second coordinate is therefore

At $t=0$ this gives $1/\sqrt{2}=0.7071$, the first row of the table. At $t=8$ it is $1/\sqrt{65537}$, which is $0.0039$ to four decimal places, matching the table. At $t=10$ it is $1/1024$, about $0.00098$. The ratio $|{\lambda }_2/{\lambda }_1|$ is $1/2$, so each step halves the contamination from the second direction. Ten steps divide it by about one thousand.

One practical bonus: if the rescaling divides by the largest entry instead of by the length, that entry tends to the absolute value of the dominant eigenvalue [S:28:136]. The same iterates then give both the direction and the size.

Rescaling is not cosmetic. Without it, the length grows like $|{\lambda }_1{|}^t$ and passes the largest number a computer can hold once $t{\log }_{10}|{\lambda }_1|$ passes the exponent of that number, which is about three hundred. For $|{\lambda }_1|=2$ that takes about a thousand steps, for $|{\lambda }_1|=10$ about three hundred, and for $|{\lambda }_1|=10^{10}$ about thirty. Rescaling is necessary for large $|{\lambda }_1|$ and for long runs, not within a few dozen steps for every eigenvalue above one.

The error is everything in the iterate that is not a multiple of $v_1$. In the bracket, that part is everything after $c_1{v}_1$, and every term in it carries a factor of the form $({\lambda }_i/{\lambda }_1{)}^t$. The largest such factor belongs to the second eigenvalue, so it sets the pace.

The same count tells you when to stop. If you need the direction to one part in a thousand, you need about ten steps here. If the two eigenvalues were $2$ and $1.9$, the ratio would be $0.95$, and the number of steps for that accuracy is the logarithm of one thousandth divided by the logarithm of $0.95$, which is about $135$. A small gap is expensive.

The method needs two hypotheses for the direction to converge, and it fails without either of them. The largest modulus must be attained by exactly one eigenvalue, and the start must have a nonzero component along its eigenvector. A positive real dominant eigenvalue is not needed for the direction: with ${\lambda }_1<0$ the iterates still become parallel to $v_1$, with the sign alternating. It is needed only if the unnormalised iterates themselves, rather than the line, are to settle.

If $\lambda$ and $-\lambda$ are the two largest eigenvalues, the iterates alternate between two directions forever [S:28:136-137]. A complex pair makes a real start rotate. The method also returns one eigenvector, not the whole spectrum.

There is a cheaper failure to notice: if the start is exactly orthogonal to the dominant eigenvector, the coefficient $c_1$ is zero and the method converges to the wrong direction. In exact arithmetic that is a real possibility; with floating-point rounding it is unlikely, but it is the reason the hypothesis is stated rather than assumed [S:28:135-137].

The formal statement. Let $A$ be an $n\times n$ matrix with eigenvalues ordered $|{\lambda }_1|>|{\lambda }_2|\ge \cdots \ge |{\lambda }_n|$. Let $v_1$ be a unit eigenvector for ${\lambda }_1$, and let the start $v_0$ be nonzero with coefficient $c_1\ne 0$ in the expansion $v_0={\sum }_i{c}_i{v}_i$. Define $w_k=A{v}_{k-1}$ and $v_k=w_k/\Vert w_k\Vert$. Then the line spanned by $v_k$ converges to the line spanned by $v_1$, and every other component shrinks by the factor $|{\lambda }_i/{\lambda }_1|$ per step. When ${\lambda }_1>0$ the vectors themselves converge. When ${\lambda }_1<0$ they alternate in sign, and only the line converges. The method needs two hypotheses for the line to converge: the largest modulus is attained by exactly one eigenvalue, and the start has a nonzero component along that eigenvector. A basis of eigenvectors is what makes the expansion $v_0={\sum }_i{c}_i{v}_i$ and the clean ratio $|{\lambda }_i/{\lambda }_1|$ available; it is sufficient, not necessary. Solomon states that power iteration is still guaranteed to converge in the defective case, where the ratio governs the exponential part of the error and the Jordan blocks add polynomial factors [S:28:136].

The derivation. Expand the unnormalised iterate: $A^k{v}_0={\sum }_i{c}_i{\lambda }_i^k{v}_i={\lambda }_1^k\left(c_1{v}_1+{\sum }_{i\ge 2}{c}_i{\left({\lambda }_i/{\lambda }_1\right)}^k{v}_i\right)$. The overall factor ${\lambda }_1^k$ does not change direction. Every ratio inside the bracket has modulus below one, so the bracket tends to $c_1{v}_1$ and the line tends to the line through $v_1$. The omitted step is the passage from the bracket to the normalised iterates. It follows because the bracket stays away from zero once $k$ is large, and the normalisation divides by a positive number.

The edge case. Suppose ${\lambda }_2=-{\lambda }_1$. The ratio then has modulus one and never decays. The iterates alternate between two directions, and no amount of iteration settles them [S:28:136-137]. The core states this. What the core does not state is that the same failure appears whenever the two largest moduli are equal, including a complex conjugate pair. There the iterates rotate instead of alternating, the direction has no limit at all, and a stopping rule based on the change between consecutive iterates never fires. The usual repairs are a shift or inverse iteration [S:28:137].

What the core hid. The core presented $|{\lambda }_2/{\lambda }_1|$ as the convergence factor and used it to estimate a step count. That is the asymptotic rate. The first iterates can move away from $v_1$ before they approach it, when the coefficient $c_2$ is much larger than $c_1$. The rate describes the tail, not the whole path. The core also hid what a small gap costs. A shift moves the eigenvalues to ${\lambda }_i-\sigma$ and changes the ratio, which is why shifting is the first repair to try [S:28:137]. Rayleigh quotient iteration is the second, and it usually needs fewer steps [S:28:138].

Where the analogy breaks. "The direction converges" is true of the line, not of the arrow. With a negative dominant eigenvalue the arrow flips at every step. Reporting $v_k$ as "the eigenvector" without a sign convention is therefore ambiguous, and the ambiguity is not numerical noise.

Original sources. The algorithm appears in Solomon as two listings, one without normalisation and one with it [S:28:136]. The same section states the convergence factor and the failure when $\lambda$ and $-\lambda$ tie for largest magnitude [S:28:135-136]. The same loop computes the limiting distribution of a Markov chain, and that is how PageRank is computed [S:57:199].

An open question. The method returns one vector, and its cost is set by the gap to the second eigenvalue. The field has not settled how to update that vector when the matrix changes a little, as a graph gains and loses edges. Recomputing from scratch is correct and wasteful, and incremental updates carry no general error bound. The difficulty is that one edge can move the leading vector by an amount that depends on the whole spectrum. Evidence that would settle a given method: a stream of graph edits with the exact eigenvector recomputed at checkpoints, and a measured error that stays below a stated tolerance between checkpoints.

A frontier result. PageRank starts from a uniform vector and multiplies by the transition matrix until the vector stops changing. For the Web, 50 to 75 iterations suffice at double precision [S:57:199]. The same chapter explains why iteration replaced elimination at this scale. The matrix is sparse, and each round costs one matrix-vector product [S:57:200]. Eigenvector centrality uses the identical loop on an adjacency matrix [S:55:20]. In both cases the matrix is huge and only the leading direction is wanted.

What to read next. Read Solomon 6.3 for inverse iteration, shifting and the Rayleigh quotient, which are the three repairs for a small gap [S:28:137-138]. Read the PageRank chapter for the scale argument [S:57:199]. Then read the centrality section of Graph Representation Learning for the same iteration used as a definition rather than an algorithm [S:55:20].