[Math] Coding the Power Method in Python

linear algebranumerical linear algebrapython

I am a very beginner to coding and am trying to write a code for the power method/von Mises iteration for finding the largest eigenvalue (in magnitude) of a given matrix in Python.

I'm referencing Numerical Linear Algebra by Holger Wendland, where he gives his definition of the power method in the following way:

First, we pick a start vector $y_0 = \sum_{j=1}^n c_jw_j$ with $c_1 \neq 0$ and set $x_0 = y_0/||y_0||_2$. Then, for $m = 1, 2, \dots$ we define

1) $y_m = Ax_{m-1}$

2) $x_m = \sigma_m y_m / ||y_m||_2$, where the sign $\sigma_m \in \{-1, 1\}$ is chosen such that $x_m^T x_{m-1} \geq 0$.

He then states

Since the algorithm can directly be derived from this definition, we refrain from stating it explicitly.

To a beginner like myself, it is not so obvious how it can be derived from that definition, and I am struggling to code with this definition. I've referenced several other articles on it to try and write my code, to no avail. Can anyone help me understand how to code this into Python using NumPy?

Thanks!

Best Answer

I'm not sure what is your difficulty... The steps you have posted are already a pseudo-code for your python algorithm.

  1. You pick an in initial guess, $y_0$, for an eigenvector associated to $\lambda$, the largest eigenvalue, and normalize it, obtaining $x_0$. Normally, using a vector of random numbers is ok. If $y_0$ (or some $y_k$, for that matter) is an eigenvector associated to other eigenvalue, the method will fail to converge to the expected eigenvalue.

  2. For $m \ge 1$ you compute $y_m = A x_{m-1}$, and then normalize it computing $x_m=\sigma_m \frac{y_m}{\|y_m\|_2}$. The introduction of $\sigma_m$ is not really necessary in terms of convergence to the dominant eigenvalue, but it ensures that the sequence of eigenvectors is actually convergent. The approximation of the eigenvalue is then given by $\lambda_m = \langle A x_m, x_m\rangle$.

Related Question