Julia言語 実対称行列の固有値を大きい順にr個抽出する．

Julia標準搭載のeigenを使うと，N×N行列Aの固有値分解がO(N^3)になっちゃうと思ってたからArpack使ってたんだけど，Arpackはそもそも疎行列向きっぽい． どうやら入力行列Aが実対称行列ならば

val, vec = eigen(A, 1:r)

で小さい順に固有値をr個とりだせるみたい．大きい順に取り出したい場合は

val, vec = eigen(A, N-r+1:N)

でいけるっぽい．ついでにArpackでは，

val, vec=eigs(A, nev=r, which=:LR);

で大きい順にrの固有値を計算する．

ということで，N=3500にして，rをうごかしてみて，計算時間がどんな感じになるかちょっと試してみた．

using LinearAlgebra
using Arpack

N = 3500
A = Symmetric(rand(N,N));

@time _,_ = eigen(A, N-5+1:N);
@time _,_ = eigen(A, N-10+1:N);
@time _,_ = eigen(A, N-20+1:N);
@time _,_ = eigen(A, N-40+1:N);
@time _,_ = eigen(A, N-80+1:N);
@time _,_ = eigen(A, N-160+1:N);
@time _,_ = eigen(A, N-320+1:N);
@time _,_ = eigen(A, N-640+1:N);
@time _,_ = eigen(A, N-1280+1:N);
@time _,_ = eigen(A, N-2560+1:N);

すると，こんな感じ．やっぱり，rを大きくすると徐々に計算量が増えていく．でも理論通りなら，計算量はO(N² r)なので，倍々で計算時間が増えていくと思うんだけどなぁ，おかしいな．

  2.236443 seconds (21 allocations: 138.497 MiB, 1.99% gc time)
  2.287247 seconds (21 allocations: 138.612 MiB)
  2.423010 seconds (21 allocations: 138.841 MiB)
  2.569268 seconds (21 allocations: 139.299 MiB)
  2.423246 seconds (21 allocations: 140.215 MiB)
  2.690961 seconds (21 allocations: 142.046 MiB, 0.28% gc time)
  3.180988 seconds (21 allocations: 145.710 MiB)
  4.213903 seconds (21 allocations: 153.036 MiB, 0.15% gc time)
  9.066842 seconds (21 allocations: 167.690 MiB)
 26.820623 seconds (20 allocations: 196.996 MiB, 0.03% gc time)

ちなみに

@time _,_ = eigen(A);

だと， 5.097789 seconds (19 allocations: 207.070 MiB, 0.12% gc time) でした．つまり，rがある程度大きい場合は，全部求めちゃった方が速かったりするってことですね．

次に，Arpackの結果はこんな感じ．

@time _,_ = eigs(A,nev=5, which=:LR);
@time _,_ = eigs(A,nev=10, which=:LR);
@time _,_ = eigs(A,nev=20, which=:LR);
@time _,_ = eigs(A,nev=40, which=:LR);
@time _,_ = eigs(A,nev=80, which=:LR);
@time _,_ = eigs(A,nev=160, which=:LR);
@time _,_ = eigs(A,nev=320, which=:LR);
@time _,_ = eigs(A,nev=640, which=:LR);
@time _,_ = eigs(A,nev=1280, which=:LR);
@time _,_ = eigs(A,nev=2560, which=:LR);

  0.548882 seconds (2.33 k allocations: 775.391 KiB)
  1.044984 seconds (4.04 k allocations: 983.375 KiB)
  0.806430 seconds (3.35 k allocations: 1.632 MiB)
  1.210613 seconds (4.00 k allocations: 3.070 MiB)
  1.765235 seconds (5.56 k allocations: 6.031 MiB)
  3.085386 seconds (8.10 k allocations: 12.222 MiB)
  8.542838 seconds (12.02 k allocations: 25.732 MiB)
 24.201966 seconds (16.58 k allocations: 57.314 MiB)
 101.972415 seconds (23.18 k allocations: 139.133 MiB)
  97.998941 seconds (18.07 k allocations: 196.968 MiB, 0.03% gc time)

いまAは密行列ですが，rが小さいときは，疎行列用のArpackの方が速いという結論ぽいですね．ちなみにeigen(A)とすると，デフォルトで6個だけ固有値を抽出する仕様っぽいです．

ちなみに，KrylovKit.jlというライブラリもあって，こちらは

using KrylovKit
@time _,_,_ = eigsolve(A, 5, :LR);
@time _,_,_ = eigsolve(A, 10, :LR);
@time _,_,_ = eigsolve(A, 20, :LR);

 0.322175 seconds (2.65 k allocations: 15.157 MiB)
 0.371782 seconds (3.62 k allocations: 20.360 MiB)
 0.899626 seconds (12.60 k allocations: 71.841 MiB)

という風にeigsolve(A, r, :LR)でr個の固有値を取り出します．rを大きくすると，krylovdimが小さいとエラーがでるのですが，krylovdimをどれくらい大きくするのが最適なのかよくわからなかったです．