Randomized eigenvalue decomposition

Question

I need the implementation of randomized eigenvalue decomposition. This algorithm can find the largest r-eigenvalues of a symmetric real n×n matrix A rapidly.

Unfortunately, RandomizedLinAlg.jl is out of work, so I had to implement it by myself.

I used the pseudocode in the paper Randomized algorithms for Generalized Hermitian Eigenvalue Problems with application to computing Karhunen-Loève expansion. This is the link and I also pasted the pseudocode as an image.

My implementation

This is my implementation.

using LinearAlgebra

function rEigen_singlepass(A,Ω)
    n, m = size(Ω)
    Y = A*Ω
    Q, R = qr(Y)
    Q = Q * Matrix(I, n, m)
    T = (Q' * Y) * ( (Q' * Ω)\I)
    λ, S = eigen(T)
    Λ = Diagonal(λ)
    U = Q*S
    return U, λ
end

function rEigen_twopass(A,Ω)
    n, m = size(Ω)
    Y = A*Ω
    Q, R = qr(Y)
    Q = Q * Matrix(I, n, m)
    T = Q' * A * Q
    λ, S = eigen(T)
    Λ = Diagonal(λ)
    U = Q*S'
    return U, λ
end

It works for r=n.

The outputs of my functions rEigen_singlepass and rEigen_twopass are the same as the outputs of eigen(A). It is good!

# example 
A = [1 2 3 4 8;
     2 5 6 2 4;
     3 6 4 9 2;
     4 2 9 0 6;
     8 4 2 6 3
]
Ω = randn(size(A)[1], r);

U, λ = rEigen_singlepass(A,Ω);
@show λ
# -9.281615829881577
# -5.421965974315312
#  1.7288438279829894
#  4.767834782574994
# 21.20690319363888

U, λ = rEigen_twopass(A,Ω);
@show λ
# -9.28161582988157
# -5.421965974315315
#  1.7288438279829903
#  4.767834782575
# 21.20690319363886

@show eigen(A).values[end-r+1:end]
#  -9.281615829881577
# -5.421965974315318
#  1.7288438279829883
#  4.767834782574983
# 21.20690319363888

It does not work for r≠n.

However, in the case of r<n, it does not work. The result of rEigen_singlepass and rEigen_twopass are not the same as the outputs of eigen(A). See the below example for r=2.

#example
A = [1 2 3 4 8;
     2 5 6 2 4;
     3 6 4 9 2;
     4 2 9 0 6;
     8 4 2 6 3
]
r = 2
Ω = randn(size(A)[1], r);

U, λ = rEigen_singlepass(A,Ω);
@show λ
# 15.217678410755674
# 94.74613399123132

U, λ = rEigen_twopass(A,Ω);
@show λ
# -2.260970145245041
# 20.64784624148779

@show eigen(A).values[end-r+1:end]
#  4.767834782574983
# 21.20690319363888

It means that my code has some wrong. Where did I make mistake?

I also refer to the original paper on randomized eigenvalue decomposition. This is the link.

EDIT: Example for large n. (response to a comment by BadZen)

This is an example of a large n. I define r=150 but let's see only the 10-largest eigenvalues.

The results of rEigen_singlepass and rEigen_twopass are not the same as the outputs of eigen(A) even for large n and small r. I feel my implemantion has problmes.

#example for large random matrix.
n = 2000
A = Symmetric(rand(n,n))
r = 150
Ω = randn(size(A)[1], r)

U, λ = rEigen_singlepass(A,Ω);
@show λ[end-10+1:end]
#  207.66110043406184
#  252.37543117950105
#  286.9021092350813
#  316.9180089435316
#  480.174283312442
#  570.4399205164796
#  670.1301653437732
# 1030.1556622604826
# 1684.7577297400524
# 2070.861149510907

U, λ = rEigen_twopass(A,Ω);
@show λ[end-10+1:end]
#  11.502794808482616
#  11.584253591602952
#  12.012179355201381
#  12.035950001524203
#  12.349115667274983
#  12.659437107692673
#  12.86745333781016
#  13.103010410849544
#  13.422670031703078
# 998.2723066471225

@show eigen(A).values[end-10+1:end]
#   24.86757230301933
#   24.901746966804254
#   24.940305083695115
#   24.98139716896731
#   25.15440459011352
#   25.269162557644382
#   25.312282048911026
#   25.374796985742325
#   25.572648646590533
# 1000.1370202221857

Answer 1

Apparently, my code has no problem. According to the original N.Halko paper, randomized eigenvalue decomposition is effective only when the eigenvalues of the input matrix are decrying. Let us define A as such a matrix. The following experiment shows the outputs of my functions are approximately the same as the output of lapack.

n = 30
r = 6
A = Symmetric(rand(n,n))
B = zeros(n,n)
e, v = eigen(A)
for i in 1:n
    B += e[i] * v[:,i] * v[:,i]' * (i>n-r ? 1 : 1e-3)
end
A = Symmetric(B);
? = randn(size(A)[1], r)

U, ? = rEigen_singlepass(A,?);
?[end-r+1:end]
#  1.6972269946076446
#  2.0678582596060338
#  2.114020877070649
#  2.645103836439442
#  2.80298502359695
# 14.78544974888867

U, ? = rEigen_twopass(A,?);
?[end-r+1:end]
#  1.6883681640086117
#  2.0536447837869174
#  2.1103655428699253
#  2.6438950659892435
#  2.797125374675521
# 14.78174207245349

eigen(A).values[end-r+1:end]
#  1.6883689201446117
#  2.053646020688538
#  2.1103671910241286
#  2.643897161145393
#  2.7971278616415454
# 14.781742117602374