はじめに
機械学習ではトレースノルム正則化やLatent Semantic AnalysisなどでSVDが用いられます。行列が大規模になってくるとフルのSVDを求めるのはとても時間が掛かってしまいます。
『前処理による特異値分解の高速化』、『Juliaで『前処理による特異値分解の高速化』』では、乱数ベクトルを用いてSVDしたい行列をその行列のランク（行列の近似誤差がεで抑えられるランク）に射影して、行列を小さくしてからSVDする手法（あくまで前処理）を紹介しました。
この手法の欠点は、行列のランクがそもそも大きい場合には行列を小さくできない、あるいは無理やり小さくしてSVDすると求まる特異値や特異ベクトルの精度が劣化してしまうことです。
今回は前処理ではなくSVD本体の部分で、フルのSVDでなく上位k本だけを取り出すTruncated SVDによる高速化を検討します。
『JuliaのKrylov部分空間法ライブラリ、KrylovKit.jlの紹介』で紹介されているKrylovKitのsvdsolveを用いてTruncated SVDしてみましたので、紹介させて頂きます。
JuliaのKrylovKitパッケージでTruncated SVD
KrylovKitでTruncated SVDする場合svdsolve(A, k; krylovdim)を用います。Krylov部分空間法も、問題を小さい次元に射影してから解く手法で、その次元をkrylovdimで指定します。krylovdim=30がデフォルトになっているようで、krylovdimはkより大きい必要があるので、30以上のkを指定する場合にはkrylovdimも明示的に指定する必要があります。
k=20, krylovdim=30の場合、次のように実行します。
using KrylovKit
A = randn(3000, 1000)
vals, lvecs, rvecs = svdsolve(A, 20, krylovdim=30)
krylovdimの指定の目安
krylovdim=k+10で指定した場合の、「フルのSVDで求まる特異値との誤差の最大値」は次のようになります。kを大きくしていっても、krylovdim=k+10程度に設定しておけばフルのSVDとほぼ同様の結果が得られています。
> using LinearAlgebra, KrylovKit
> A = randn(3000, 3000)
> P, q, R = svd(A)
>
> for k in [10, 50, 100, 150, 200, 250, 300, 350, 400, 450, 500]
> krylovdim=k+10
> ss, us, vs = svdsolve(A, k, krylovdim=krylovdim)
> println("k=", k, ", krylovdim=", krylovdim, "\t", maximum(abs.(ss[1:k]-q[1:k])))
> end
k=10, krylovdim=20 1.4352963262354024e-12
k=50, krylovdim=60 1.2647660696529783e-12
k=100, krylovdim=110 5.172751116333529e-12
k=150, krylovdim=160 7.261746759468224e-12
k=200, krylovdim=210 4.050093593832571e-12
k=250, krylovdim=260 6.0538241086760536e-12
k=300, krylovdim=310 3.254285729781259e-12
k=350, krylovdim=360 3.623767952376511e-12
k=400, krylovdim=410 2.7853275241795927e-12
k=450, krylovdim=460 2.3732127374387346e-12
k=500, krylovdim=510 2.9416469260468148e-12
フルのSVDとの速度比較
だいたいの目安として、SVDの対象の行列の次元に対してkが5%程度までであればsvdsolveを用いた方が速く、10%程度になってくるとフルのSVDをした方が速い、という結果になりそうです。
A=randn(100,100)
LinearAlgebra.svd(A) 5.609 ms (16 allocations: 482.02 KiB)
KrylovKit.svdsolve(A, k, krylovdim)
k=10 krylovdim=20 6.840 ms (1706 allocations: 564.78 KiB)
k=50 krylovdim=60 29.501 ms (5463 allocations: 2.06 MiB)
A=randn(500,500)
LinearAlgebra.svd(A) 121.581 ms (17 allocations: 11.51 MiB)
KrylovKit.svdsolve(A, k, krylovdim)
k=10 krylovdim=20 26.460 ms (3701 allocations: 4.62 MiB)
k=50 krylovdim=60 137.964 ms (16699 allocations: 22.45 MiB)
k=100 krylovdim=110 435.836 ms (32790 allocations: 45.79 MiB)
k=200 krylovdim=210 2.200 s (71405 allocations: 104.50 MiB)
k=300 krylovdim=310 3.272 s (70057 allocations: 108.71 MiB)
A=randn(1000,1000)
LinearAlgebra.svd(A) 662.412 ms (17 allocations: 45.90 MiB)
KrylovKit.svdsolve(A, k, krylovdim)
k=10 krylovdim=20 115.463 ms (6023 allocations: 14.17 MiB)
k=50 krylovdim=60 392.134 ms (22167 allocations: 56.18 MiB)
k=100 krylovdim=110 1.157 s (51502 allocations: 132.54 MiB)
k=200 krylovdim=210 5.059 s (104333 allocations: 273.74 MiB)
k=300 krylovdim=310 9.600 s (150208 allocations: 393.20 MiB)
A=randn(2000,2000)
LinearAlgebra.svd(A) 4.277 s (17 allocations: 183.35 MiB)
KrylovKit.svdsolve(A, k, krylovdim)
k=10 krylovdim=20 654.430 ms (6830 allocations: 31.58 MiB)
k=50 krylovdim=60 3.227 s (32212 allocations: 157.72 MiB)
k=100 krylovdim=110 6.690 s (84249 allocations: 415.87 MiB)
k=200 krylovdim=210 15.591 s (175471 allocations: 859.24 MiB)
k=300 krylovdim=310 24.699 s (208094 allocations: 1010.38 MiB)
A=randn(3000,3000)
LinearAlgebra.svd(A)(A) 14.424 s (18 allocations: 412.35 MiB)
KrylovKit.svdsolve(A, k, krylovdim)
k=10 krylovdim=20 1.937 s (10053 allocations: 53.07 MiB)
k=50 krylovdim=60 5.929 s (50573 allocations: 279.86 MiB)
k=100 krylovdim=110 10.858 s (91099 allocations: 511.03 MiB)
k=200 krylovdim=210 29.131 s (215364 allocations: 1.17 GiB)
k=300 krylovdim=310 53.876 s (394399 allocations: 2.08 GiB)