Eigenvalues via rootfinding: Variational method with Lagrange multipliers

Introduction

One can find an eigenvalue and the corresponding eigenvector of a real symmetric matrix A by searching for the stationary point of the following constrained form

E(v)=vTAv , where vTv=1 .

This can be done using the method of Lagrange multipliers: the original problem is reformulated into unconstrained search for the stationary point (the point where all partial derivatives are zero) of the Lagrangian function defined as

L(v,λ) = vTAv - λ(vTv-1) .

in the space {v,λ}. At the stationary point: λ is the eigenvalue and v is the corresponding eigenvector normalized to unity. Indeed, at the stationary point

∂L(v,λ)/∂v = 2(Av - λv) = 0 ,
∂L(v,λ)/∂λ = vTv-1 = 0 .

Of course this is equivalent to finding the root in the n+1 dimensional space x={v,λ} of the n+1 dimensional vector-function

f(x) = {Av-λv, vTv-1}

where n is the size of v.

Note that the Jacobian of this function is analytic,

∂fi(x)/∂vj = Aij-λδij , where i,j≤n ,
∂fi(x)/∂λ = -vi , where i≤n ,
∂fn+1(x)/∂vj = 2vj , where j≤n 
∂fn+1(x)/∂λ = 0 .

The analytic Jacobian greatly facilitates the search for the root.

Project

Implement a subroutine that calculates the eigenvalue (close to a given value λ_start) and the corresponding eigenvector of a given symmetric matrix A by searching for the root of the function

f(v,λ) = {Av-λv, vTv-1}

using the components of the vector v and the (Lagrangian multiplier) λ as free parameters. Use your own Newton's rootfinding routine and customize it to use analytic Jacobian.

Extra

Investigate the scaling of this method (time to find the eigenvalue as function of the size of the marix) and find out how far in the size of the matrix you can realistically go (that is, within few seconds of calculation time).

Calculate several lowest state of the Hydrogen atom.

It seems that the line-search here can be done analytically, you might want to try to investigate this possibility.