Randomized eigenvalue decomposition

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

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