The power method is one of the simplest methods to find the maximum eigenpair of an operator. We simply compute \(x_{k+1}=\frac{Ax_k}{\left\lVert Ax_k \right\rVert}\) until the corresponding Rayleigh quotient \(\frac{x^TAx}{x^Tx}\) stops changing.

This simple algorithm requires only the action of the operator. We start at iteration \(k=0\) using a random initial vector \(x_0\). Then, until some maximum number of iteration steps, we compute \(Ax\) and the current Rayleigh quotient. We then normalize \(x\). If the change from the previous Rayleigh quotient is sufficiently small, we break out of the loop. At the end of the iterations, our vector \(x\) should contain the maximum eigenvector. We return the corresponding maximum eigenvalue.

Implementation

Using linalgcpp, we implement the above steps.

double PowerMethod(const Operator& A,
                  int max_iter = 1000,
                  double tol = 1e-8)
{
    Vector x(A.Rows());
    Vector y(A.Rows());

    x.Randomize(-1.0, 1.0);
    x.Normalize();

    double rayleigh = 0.0;
    double old_rayleigh = 0.0;

    for (int i = 0; i < max_iter; ++i)
    {
        A.Mult(x, y);

        rayleigh = (x * y) / (x * x);

        std::swap(x, y);
        x /= x.L2Norm();

        double change_ratio = std::fabs(rayleigh - old_rayleigh) /
                              std::fabs(rayleigh);

        if (change_ratio < tol)
        {
            break;
        }

        old_rayleigh = rayleigh;
    }

    return rayleigh;
}

Distributed Environment

A very similar function will work for distributed operators. The vector dot products should be replaced with parallel versions.