まぃふぇいばりっと

機械学習やってます．Julia大好きです．勉強したことの殴り書きです．

Julia言語 Randomized eigenvalue decomposition via Nystrom Method を実装した．

Julia言語 行列分解 論文実装

n×nの正定値対称行列Aの固有値を大きい順にr個求めたいです．これを高速に計算するアルゴリズムとしてRandomized eigenvalue decomposition via Nystrom Methodを実装しました．以下の論文のFig.6.を実装しました．

arxiv.org

前に実装したREDは固有値を大きい順にr個だけ高速に計算できましたが，固有値が減衰していない行列については，（Fノルムで評価する低ランク近似はうまくできても）固有値の近似がよくないアルゴリズムでした．今回のRED via Nystrom Methodは入力が正定値対称行列に縛られますが，高速にかつより正確に固有値を大きい順に求めることができます．

とくに困ったところはなく，半日くらいで順調に実装できました．

using LinearAlgebra
function evdnys(A::Symmetric, k::Int;p=10)
    n, n = size(A)
    G = randn(n, k+p)
    Y = A * G
    Q = qr(Y).Q*Matrix(I,n,k+p)
    B1 = A*Q
    B2 = Q'*B1
    C = cholesky(Symmetric(B2))
    F = B1 * inv(C.U)
    U, Σ, _ = svd(F)
   #reconstract is U*diagm(Σ.^2)*U'
    return U, Σ.^2
end

サンプリングパラメータのpは適当に10にしていますが，p=2rにするのがよい，といっている文献もどこかにありました．

計算精度について

実験用にとりあえず，適当な正定値対称行列Aをつくります．

using Statistics
m = 50
n = 100
x = randn(n,m)
demean_x = x .- mean(x, dims=1)
A = demean_x'*demean_x ./(n-1)
isposdef(A) #true

比較のために前回実装した，EVDを用意しておきます．

function reigen(A,r;p=10, twopass=false)
    n = size(A)[1]
    m = r + p
    Ω = randn(n, m)
    n, m = size(Ω)
    Y = A*Ω
    if twopass
        Q = qr!(Y).Q * Matrix(I,n,m)
        T = Q' * A * Q
    else
         Q = qr(Y).Q * Matrix(I,n,m)
         T = (Q' * Y) * ( (Q' * Ω)\I)
    end
    λ, S = eigen(T)
    Λ = Diagonal(λ)
    U = Q*S'
    return U, λ
end

まず，lapackで行列Aの固有値を大きい順に求めて，10個表示させてみます．

result_lapack = eigen(A).values
reverse(result_lapack)[1:10]
# 2.499093654863234
# 2.449174582223532
# 2.3513196905628972
# 2.2851514992705284
# 2.1091186508409545
# 2.069503593052392
# 1.9434150310107634
# 1.804034622747044
# 1.7192090439026784
# 1.6285715863497319

適当にr=20にして，EVDで行列Aの固有値を大きい順に10個求めてみます． 結構激しくずれてます．

r=20
U, values = reigen(A,r)
reverse(values)[1:10]

# 2.3648795631492248
# 2.2652936494283664
# 2.065460162647662
# 2.053243332161027
# 2.0015998322530044
# 1.9269963048738135
# 1.8328635753401925
# 1.6750657308299814
# 1.611119258198799
# 1.542245495794801

では，今回実装したRED via Nystrom Methodで行列Aの固有値を大きい順に10個求めてみます．

r = 20
_, eigenvalues = evdnys(Symmetric(A),r)
display(eigenvalues[1:10])

# 2.499093654863234
# 2.449174582223532
# 2.3513196905628972
# 2.2851514992705284
# 2.1091186508409545
# 2.069503593052392
# 1.9434150310107634
# 1.804034622747044
# 1.7192090439026784
# 1.6285715863497319

lapackの結果をよく近似できているように思います．

計算速度について

既存の手法との比較をしてみました．

function evdnys(A::Symmetric, k::Int;p=10)
    n, n = size(A)
    G = randn(n, k+p)
    Y = A * G
    Q = qr(Y).Q*Matrix(I,n,k+p)
    B1 = A*Q
    B2 = Q'*B1
    C = cholesky(Symmetric(B2))
    F = B1 * inv(C.U)
    U, Σ, _ = svd(F)
    Λ = diagm(Σ.^2)
    return U, Σ.^2
end

function reigen(A,r;p=10, twopass=false)
    n = size(A)[1]
    m = r + p
    Ω = randn(n, m)
    n, m = size(Ω)
    Y = A*Ω
    if twopass
        Q = qr!(Y).Q * Matrix(I,n,m)
        T = Q' * A * Q
    else
         Q = qr(Y).Q * Matrix(I,n,m)
         T = (Q' * Y) * ( (Q' * Ω)\I)
    end
    λ, S = eigen(T)
    Λ = Diagonal(λ)
    U = Q*S'
    return U, λ
end

using Statistics
m = 2000
n = 10000
x = randn(n,m)
demean_x = x .- mean(x, dims=1)
A = demean_x'*demean_x ./(n-1)
isposdef(A) #Aは正定値対称行列

p0 = 5
rs = [ p0*2^(n-1) for n = 1:7];
for r in rs
    result_lapack, runtime_lapack = @timed eigen(Symmetric(A), m-r+1:m);
    result_arpack, runtime_arpack = @timed eigs(A, nev=r, which=:LR);
    result_krylov, runtime_krylov = @timed eigsolve(A, r, :LM, krylovdim=r+10);
    result_lobpcg, runtime_lobpcg = @timed lobpcg(A, true, r);
    result_reigen1, runtime_reigen2 = @timed reigen(A, r, twopass=false);
    result_reigen2, runtime_reigen2 = @timed reigen(A, r, twopass=true);
    result_evdnys, runtime_evdnys = @timed evdnys(Symmetric(A), r);
    println(
       "r $r
        lapack  $runtime_lapack \t arpack $runtime_arpack \t 
        krylov $runtime_krylov \t lobpcg $runtime_lobpcg
        reigen1 $runtime_reigen1 \t reigen2 $runtime_reigen2 \t evdnys $runtime_evdnys")
    println("******")
end

結果は以下の通り．やっぱりEVDviaNystromをはじめ，Randomizedな方法は速いですね．

r 5
        lapack  0.8177911    arpack 0.4960746    krylov 0.3070053    lobpcg 0.3421401
        reigen1 0.0058875    reigen2 0.0104825   **evdnys 0.011699**
******
r 10
        lapack  0.8837821    arpack 0.576648     krylov 0.3857683    lobpcg 0.4529197
        reigen1 0.0101869    reigen2 0.0106938   **evdnys 0.0127273**
******
r 20
        lapack  0.8577545    arpack 0.590889     krylov 0.4786105    lobpcg 0.5592837
        reigen1 0.0080609    reigen2 0.0101375   **evdnys 0.0151632**
******
r 40
        lapack  1.1575799    arpack 0.7349361    krylov 0.722302     lobpcg 0.8956097
        reigen1 0.0125024    reigen2 0.0159989   **evdnys 0.028382**
******
r 80
        lapack  0.9327568    arpack 1.1262472    krylov 1.2329233    lobpcg 2.6357616
        reigen1 0.0242507    reigen2 0.028978    **evdnys 0.051772**
******
r 160
        lapack  1.1156746    arpack 2.0878638    krylov 2.3862248    lobpcg 4.6012995
        reigen1 0.0991418    reigen2 0.1222233   **evdnys 0.1450164**
******
r 320
        lapack  1.4092244    arpack 5.3995158    krylov 4.1578276    lobpcg 15.1888958
        reigen1 0.3209115    reigen2 0.7082315   **evdnys 0.4270121**
******

reigen1,reigen2,evdnysでは，最初にランダム行列と入力行列の積を計算しているんですけど，これには，O(nmr)くらいの計算量が必要です．しかし，N.halko論文やこの論文によると，O(nmlogr)に減らせるらしい．fast Johnson-Lindenstrauss transforms とか，SRFT とかよばれる技術を使うらしいがよくわからない．具体的にどうしたら計算量がへるのか分からないので，とりあえず StackOverflowに投げた．

stackoverflow.com

これをちゃんと実装できれば，reigen1,reigen2,evdnysの計算速度はさらに半分くらいになりそう．．誰か教えてくれ．．．涙