[julia-users] Re: smallest eigenvalue of sparse matrix

[Steven G. Johnson]

For Hermitian problems, you should also try the LOBPCG algorithm, which is 
not implemented in Julia yet, but which is a nice alternative to Lanczos 
for Hermitian problems.  Unfortunately, there is no Julia implementation of 
LOBPCG yet that I know of, but you can probably translate one from Matlab.

Here is a toy implementation of LOBPCG for the smallest eigenvalue that I 
put together for class (warning: if you run it for too many iterations it 
goes a little wonky because the X matrix becomes singular, though it is 
straightforward to modify it to desingularize .... also, a serious 
implementation should update X'*X etcetera in-place):

# run n iterations of conjugate-gradient Rayleigh-quotient iteration on A, 
starting from x0,
# returning the list of n eigenvalue estimates
function cgrq(A, x0, n)
    λ = Array(typeof(real(zero(eltype(A)))), n)
    xprev = x0
    x = x0
    Ax = A*x
    Axprev = Ax
    for i = 1:n
        Ad = A*Ax
        X = [xprev x Ax]
        AX = [Axprev Ax Ad]
        ν, Z = eig(X' * AX, X' * X)
        iν = sortperm(ν, by=real)[1]
        xprev = x
        Axprev = Ax
        x = X * Z[:,iν]
        Ax = AX * Z[:,iν]
        λ[i] = real(ν[iν])
    end
    return λ
end