n×nの実対称行列Aの固有値を大きい順にr<n個得たいです．lapackを使うと
val, vec = eigen(Symmetric(A), n-r+1:n)
でOKです．この時，eigenの第一引数はSymmetric型じゃないといけないので注意．（ところで，ここでeigenではどんなアルゴリズムが動いてるんですかね？全部の固有値を求める方法は調べると出てきますが，大きい順に，っていう時は何が裏で動いているんでしょう．．．？知っている人がいたら教えてほしいです．）
またArpackを使って，
using Arpack val, vec=eigs(A, nev=r, which=:LR);
でもOKでした．ArpackはKrylov部分空間法の固有値解法ライブラリで，疎行列でも動くライブラリみたいです．このブログではArpackのことを何度か「疎行列向けのライブラリ」と表現しましたがそれはミスリーディングっぽいです．ごめんなさい．ここで，eigs(A)のAは実対称行列じゃなくても動きます．
自分のいま取り組んでいる研究で，Aの固有値全部はいらないから，大きい順にr個だけを爆速で取り出したい，という要請があります．こないだちょっと試してみたところ，eigen(Symmetric(A), n-r+1:n)とeigs(A, nev=r, which=:LR)は固有値を全部求めるeigen(A)と比べても，そんなに速くないみたいでした．そこで，噂に聞く結構速いアルゴリズムであるらしい Randomized eigenvalue decomposition(RED) を調べてみました．REDの特異値分解版のRSVDに関する実装や説明は沢山あるみたいですけど，REDに関しては，少なくともJulia実装は見当たりませんでした．一応，RandomizedLinAlg.jlというJuliaのライブラリがありますが，試してみても現在は少なくとも密行列についてはreigenは動かないです．（疎行列については分からない）
最適化済みLAPACKはSIMD、スレッド並列化が両方効いて相当速いが，ARPACKは並列化がうまく効かないから遅いらしいです．一方で，REDを使うと、ARPACKと計算量は同じですが、r≠1ではBLAS3を使える分、SIMDが効いて早くなることが期待されるらしいです．これは強い人に教わりました．
そこで次の2つの論文を参考に，REDを自分で実装してみました．固有値分解の周辺のトピックが全然分からないので論文読んだり調べたりしながら実装するのは，大変でした．１０数行のコードを書くのに４日以上かかりました．．．
ちなみに，１つ目のN.Halko論文に関連して，この方のD論はこちらです．
REDのJulia実装
見よう見まねで書いたRandomized eigenvalue decompositionのjulia実装がこちら．．
function reigen(A,r;p=10, twopass=false) n = size(A)[1] m = r + p Ω = randn(n, m) 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
twopass=falseの方が速いが，より不正確です．twopass=trueの方がより正確です．pはover-sampling parameterで，この文献に従ってp=10としています．よい決め方知っている人がいたら教えてほしいです．
ちょっと詰まったのがQR分解のところですね．qr(Y).Qはそのままでは正方行列なのでn×mの矩形行列として取り出したいときはqr(y).Q*Matrix(I,n,m)とします．
論文によって，QがQR分解のQだったり，特異値分解の主要な右特異ベクトルだったりしますが，これは，Johnson-Lindenstraussの補題で正当化できるって，強い人に教わりました．
なお，N Halko 論文によると，Y=A*Ωの積は，FFTを使うと，もっとはやく計算できるっぽいです．FFTとかあんまり分からないので，とりあえずこのままで進めます．
全体の計算量については，強い方に教えてもらいました．
精度について
さて，じゃぁ，計算時間の比較の前にeigenと同じ結果を得られるか実験してみます．
適当な30×30の対称行列Aの固有値をr=6個取り出してみます．
n = 30 r = 6 A = Symmetric(rand(n,n)) U, λ = reigen(A,r,twopass=false); λ[end-r+1:end] # 1.7535736114874219 # 2.2731278141844804 # 2.778424916550649 # 2.9919084354863856 # 8.475088453581433 # 14.989196447947041 U, λ = reigen(A,r,twopass=true); λ[end-r+1:end] # 1.2197855614132485 # 1.55690645146256 # 2.0387086285543066 # 2.111767689505717 # 2.3878222860931304 # 15.152842264762363 eigen(A).values[end-r+1:end] # 1.7131238291530306 # 2.0291194695459063 # 2.2267398731934835 # 2.3956403884643858 # 2.6719356178806133 # 15.209793778533955
eigen(A)の出力結果と比べると，セルフ実装の結果が全然違うやん！！！ どこに誤りがあるのか分からなくて，StackOverflowで質問してしまったが，どこにも誤りはないです． というのも，上であげたどちらの論文でも説明されているように，REDは，固有値が減衰する行列に対して有用な手法のようです．知らなかった！！． ちなみに，固有値が減衰している場合もしていない場合も，RSVDによる行列の再構成に関しては，フロベニウスノルムの意味ではよく近似できるらしいが，特異値についてはあまりよく近似できていないっぽい．
If they decay fast, then the spectral norm error and the Frobenius norm error are similar,and the RSVD works well. If they decay slowly, then the RSVD performs fine when errors are measuredin the Frobenius norm, but not very well when the spectral norm is the one of interest. ([1607.01649] Randomized methods for matrix computations, P.G. Martinsson)
N halko論文の実験パートもよく読むと，固有値が減衰する手法についてRSVDを実験していますね．ということで，行列Aをいい感じに固有値が減衰する行列に書き直して再度実験すると，
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); U, λ = reigen(A,r,twopass=false); λ[end-r+1:end] # 1.6972269946076446 # 2.0678582596060338 # 2.114020877070649 # 2.645103836439442 # 2.80298502359695 # 14.78544974888867 U, λ = reigen(A,r,twopass=true); λ[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=1500にして，固有値が100個目以降は小さくなるような行列を用意して，
ターゲットランクrを徐々に大きくして計算時間の変化を観測しました．
using LinearAlgebra using Arpack using KrylovKit using IterativeSolvers n = 1500 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-100 ? 1 : 1e-3) end A = Symmetric(B); p0 = 5 rs = [ p0*2^(n-1) for n = 1:7]; for r in rs result_lapack, runtime_lapack = @timed eigen(A, n-r+1:n); 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_reigen1 = @timed reigen(A, r, twopass=false); result_reigen2, runtime_reigen2 = @timed reigen(A, r, twopass=true); 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") println("******") end
結果はこんな感じ．
r 5
lapack 0.478838 arpack 0.0602754 krylov 0.0251191 lobpcg 0.1881867
reigen1 0.0050092 reigen2 0.1158657
******
r 10
lapack 0.5068946 arpack 0.0926673 krylov 0.0395786 lobpcg 0.3958909
reigen1 0.0096395 reigen2 0.1367981
******
r 20
lapack 0.5186003 arpack 0.0878039 krylov 0.0501138 lobpcg 0.4257776
reigen1 0.0059776 reigen2 0.1848728
******
r 40
lapack 0.5065934 arpack 0.1264882 krylov 0.0722774 lobpcg 0.8301309
reigen1 0.0111938 reigen2 0.4549637
******
r 80
lapack 0.794966 arpack 0.1150471 krylov 0.1482746 lobpcg 0.2263155
reigen1 0.0181737 reigen2 0.7018366
******
r 160
lapack 0.8155617 arpack 0.9227155 krylov 1.2557783 lobpcg 2.1098337
reigen1 0.078549 reigen2 1.6126655
******
r 320
lapack 0.8963038 arpack 3.7780156 krylov 3.7544088 lobpcg 6.2717872
reigen1 1.1731766 reigen2 3.3344803
r<<nだとREDの圧勝ですね．めでたしめでたし．
さて，私の研究では，分解したい行列Aの固有値は別に減衰していなさそうなので，REDは使えなそう．代わりにNystrom手法による固有値分解の手法とかでも実装しようかなと思います．
