lobpcg.jl

using Random
using LinearAlgebra

function ortho(X)
    Array(qr(X).Q)
end

function align(X,Y)
    # returns an orthogonal combination X * c with c unitary such that X*c is closest to Y
    proj = X'Y
    U,S,V = svd(proj)
    any(x -> x <= .1, S) && error("Not enough overlap")
    U*V'
end

function LOBPCG(A, X, tol=1e-16, maxiter=100, do_align=true)
    N = size(X,1)
    M = size(X,2)
    X = ortho(X)
    Y = zeros(N,3M)
    P = zeros(N,M)
    R = zeros(N,M)
    niter = 1
    while true
        Y[:,1:M] = X
        R .= A*X - X*(X'*A*X)
        Y[:,M+1:2M] = R
        Y[:,2M+1:3M] = P
        Y .= ortho(Y)
        c = eigen(Symmetric(Y'*A*Y)).vectors[:,1:M]
        new_X = Y*c
        e = Array{Float64}(I,3M,M)
        do_align && (e .= e*align(e,c))
        P .= Y*(c - e)
        println(niter, " ", norm(R))
        X .= new_X
        norm(R) < tol && break
        niter = niter + 1
        niter <= maxiter || break
    end

    X
end

Random.seed!(0)
N = 1000
M = 10
A = Diagonal(randn(N))
X = ortho(randn(N,M))
LOBPCG(A,X,1e-12,1000,true)