Julia言語 SVDによる低ランク近似

Julia でSVD をするときに，Lapackを使って

using LinearAlgebra

n = 10
A = rand(n, n)
svd(A)

としても良いのですが，これだとO(n³)必要ですね． ランクをrにするのに，完全なSVDをする必要はありません．Arpackのsvdsを使うべし．

using Arpack

n = 10
r = 5
A = rand(n, n)
svd = svds(A ; nsv=r, ritzvec=true)[1]
Ar = svd.U * diagm(svd.S) * svd.Vt

で10×10行列Aの5ランク近似がO(n² r)で手に入ります．

julialinearalgebra.github.io

実際に時間はかってみるとこんな感じ． Arpackのほうがだいぶはやいが，rが大きいときは，Lapackの方がはやいっぽい．

using LinearAlgebra
using Arpack

n = 2000
A = rand(n, n)
@time svd(A);
# 4.236471 seconds (12 allocations: 183.350 MiB, 3.53% gc time)

r = 5
@time svds(A ; nsv=r, ritzvec=true)[1];
# 0.479803 seconds (1.35 k allocations: 840.125 KiB)

r = 10
@time svds(A ; nsv=r, ritzvec=true)[1];
# 0.654874 seconds (1.78 k allocations: 1.237 MiB)

r = 100
@time svds(A ; nsv=r, ritzvec=true)[1];
# 1.538005 seconds (3.57 k allocations: 11.446 MiB)

r = 200
@time svds(A ; nsv=r, ritzvec=true)[1];
# 2.692379 seconds (5.21 k allocations: 23.644 MiB)

r = 400
@time svds(A ; nsv=r, ritzvec=true)[1];
# 5.119713 seconds (7.29 k allocations: 50.736 MiB)

関連記事

Python SVDによる低ランク近似

C++ SVDによる低ランク近似