随机特征值分解
Randomized eigenvalue decomposition
我需要实施随机特征值分解。该算法可以快速找到对称实数n×n矩阵A的最大r-特征值。
可惜RandomizedLinAlg.jl没工作了，只好自己实现了
我使用了论文Randomized algorithms for Generalized Hermitian Eigenvalue Problems with application to computing Karhunen-Loève expansion中的伪代码。这是 link 我还将伪代码粘贴为图像。
我的实现
这是我的实现。
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
适用于 r=n。
我的函数 rEigen_singlepass 和 rEigen_twopass 的输出与 eigen(A) 的输出相同。不错！
# 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
它不适用于 r≠n。
但是，在r<n的情况下，它不起作用。 rEigen_singlepass 和 rEigen_twopass 的结果与 eigen(A) 的输出不同。请参阅下面的 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
这意味着我的代码有问题。我哪里弄错了？
我也参考了关于随机特征值分解的原始论文。这是 link.
编辑：大 n 的示例。 （回应 BadZen 的评论）
这是一个很大的例子n。我定义了 r=150 但我们只看 10 个最大的特征值。
rEigen_singlepass 和 rEigen_twopass 的结果与 eigen(A) 的输出不同，即使对于大 n 和小 r。感觉自己实现有问题
#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
看来，我的代码没有问题。
根据原始 N.Halko 论文，仅当输入矩阵的特征值衰减时，随机特征值分解才有效。让我们将 A 定义为这样的矩阵。下面的实验表明我的函数的输出与 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
我需要实施随机特征值分解。该算法可以快速找到对称实数n×n矩阵A的最大r-特征值。
可惜RandomizedLinAlg.jl没工作了，只好自己实现了
我使用了论文Randomized algorithms for Generalized Hermitian Eigenvalue Problems with application to computing Karhunen-Loève expansion中的伪代码。这是 link 我还将伪代码粘贴为图像。
我的实现
这是我的实现。
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
适用于 r=n。
我的函数 rEigen_singlepass 和 rEigen_twopass 的输出与 eigen(A) 的输出相同。不错！
# 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
它不适用于 r≠n。
但是，在r<n的情况下，它不起作用。 rEigen_singlepass 和 rEigen_twopass 的结果与 eigen(A) 的输出不同。请参阅下面的 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
这意味着我的代码有问题。我哪里弄错了？
我也参考了关于随机特征值分解的原始论文。这是 link.
编辑：大 n 的示例。 （回应 BadZen 的评论）
这是一个很大的例子n。我定义了 r=150 但我们只看 10 个最大的特征值。
rEigen_singlepass 和 rEigen_twopass 的结果与 eigen(A) 的输出不同，即使对于大 n 和小 r。感觉自己实现有问题
#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
看来，我的代码没有问题。
根据原始 N.Halko 论文，仅当输入矩阵的特征值衰减时，随机特征值分解才有效。让我们将 A 定义为这样的矩阵。下面的实验表明我的函数的输出与 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