はじめに
以前『前処理による特異値分解の高速化』を書きました。他の言語で実装したことはあったのですが、最近Juliaの勉強を始めたので試しにJuliaで実装してみました。
通常の特異値分解
LinearAlgebraパッケージにSVDがあります。randn(2,3)で2×3の正規乱数行列を生成してSVDしてみます。
> using LinearAlgebra
> A = randn(2,3)
> U, s, V = svd(A)
SVD{Float64,Float64,Array{Float64,2}}
U factor:
2×2 Array{Float64,2}:
-0.997385 0.0722677
0.0722677 0.997385
singular values:
2-element Array{Float64,1}:
2.179336847718715
1.7295114827398714
Vt factor:
2×3 Array{Float64,2}:
-0.477078 0.846436 0.236522
0.0305668 0.284941 -0.958058
と結果が出ます。Aを再構成できるか確認してみましょう。
> A .-= U * Diagonal(s) * V'
2×3 Array{Float64,2}:
-2.22045e-16 4.44089e-16 2.22045e-16
7.63278e-17 0.0 2.22045e-16
Juliaで .（ピリオド）は要素に対する演算を全要素に一括で適用する記法です。Aを再構成できていることがわかります。もっと大きな行列で速度を測定してみましょう。速度測定にBenchmarkToolsパッケージを用います。
> using BenchmarkTools
> A = randn(1000,500)
> @benchmark U, s, V = svd(A)
BenchmarkTools.Trial:
memory estimate: 17.23 MiB
allocs estimate: 17
--------------
minimum time: 213.673 ms (0.00% GC)
median time: 221.038 ms (0.00% GC)
mean time: 228.602 ms (2.98% GC)
maximum time: 315.837 ms (28.63% GC)
--------------
samples: 22
evals/sample: 1
となり、私の環境では1000×500の正規乱数行列のSVDに228ms程度掛かっています。
『前処理による特異値分解の高速化』のJuliaへの実装
『前処理による特異値分解の高速化』の復習
『前処理による特異値分解の高速化』を復習すると、
- 行列$\boldsymbol{A} \in \mathbb{R}^{M \times N}$（$M > N$）の正規直行基底$\boldsymbol{Q} \in \mathbb{R}^{M \times R}$を用いて圧縮した行列：$\boldsymbol{Q}^T \boldsymbol{A} \in \mathbb{R}^{R \times N}$についてSVDする。
- 求まった$\tilde{\boldsymbol{U}}$, $\boldsymbol{\sigma}$, $\boldsymbol{V}$のうち、$\boldsymbol{\sigma}$, $\boldsymbol{V}$はそのままでよく、$\tilde{\boldsymbol{U}}$だけ左から$\boldsymbol{Q}$を掛ければ元の行列のSVDが求まる（！）
- $R$は行列$\boldsymbol{A}$のランクで、$R \ll N$のとき大幅に高速化できる。
- そして、行列$\boldsymbol{A}$の正規直行基底$\boldsymbol{Q}$は、乱数ベクトル$\boldsymbol{w}_r$を用いて$\boldsymbol{y}_r = \boldsymbol{A} \boldsymbol{w}_r$（$r=1,\cdots,R$）を正規直行化すれば得られる（！）
というものでした。
Juliaへの実装
次のように実装してみました。
module RandomApprox
using LinearAlgebra
export svd
function svd(A::AbstractArray{Float64,2}, eta::Float64=1.0E-8)
nrow = size(A,1)
ncol = size(A,2)
transposed = (nrow < ncol)
if transposed
A = transpose(A)
end
qs = orthonormalbases(A, eta)
if length(qs) > 0
Q = hcat(qs...)
QT = transpose(Q)
Ut, s, V = LinearAlgebra.svd(QT * A)
U = Q * Ut
else
U = zeros(size(A,1))
s = zeros(1)
V = zeros(size(A,2))
end
if !transposed
return U, s, V
else
return V, s, U
end
end
function orthonormalbases(A::AbstractArray{Float64,2}, eta::Float64)
qs = Vector{Vector{Float64}}(undef,0)
N = size(A,2)
for j in 1:N
w = randn(N)
q = A * w
normv = gramschmidt!(q, qs)
if normv < eta
break
end
push!(qs, q)
end
return qs
end
function gramschmidt!(q::AbstractVector{Float64}, qs::AbstractVector{Vector{Float64}})
#Orthogonalization
for q2 in qs
q .-= dot(q,q2) * q2
end
#Normalization
normv = norm(q)
q ./= normv
return normv
end
end # module
以下、関数ごとに見ていきます。
svd(A, eta)
ここは
- 行列$\boldsymbol{A} \in \mathbb{R}^{M \times N}$（$M > N$）の正規直行基底$\boldsymbol{Q} \in \mathbb{R}^{M \times R}$を用いて圧縮した行列：$\boldsymbol{Q}^T \boldsymbol{A} \in \mathbb{R}^{R \times N}$についてSVDする。
- 求まった$\tilde{\boldsymbol{U}}$, $\boldsymbol{\sigma}$, $\boldsymbol{V}$のうち、$\boldsymbol{\sigma}$, $\boldsymbol{V}$はそのままでよく、$\tilde{\boldsymbol{U}}$だけ左から$\boldsymbol{Q}$を掛ければ元の行列のSVDが求まる（！）
の部分です。
function svd(A::AbstractArray{Float64,2}, eta::Float64=1.0E-8)
nrow = size(A,1)
ncol = size(A,2)
transposed = (nrow < ncol)
if transposed
A = transpose(A)
end
qs = orthonormalbases(A, eta)
if length(qs) > 0
Q = hcat(qs...)
QT = transpose(Q)
Ut, s, V = LinearAlgebra.svd(QT * A)
U = Q * Ut
else
U = zeros(size(A,1))
s = zeros(1)
V = zeros(size(A,2))
end
if !transposed
return U, s, V
else
return V, s, U
end
end
最初に行列Aが横長だったら転置して縦長にしています。
qs = orthonormalbases(A, eta)
で正規直行基底を求めて、基底が1つ以上求まったら QT * A についてのSVDしています。このSVDで求まったUtについては左からQを掛けて、sとVはそのまま返しています。基底が1つも求まらなかったらゼロ行列ですので、とりあえず空のU, s, Vを返しています。
orthonormalbases(A, eta)
ここは
- そして、行列$\boldsymbol{A}$の正規直行基底$\boldsymbol{Q}$は、乱数ベクトル$\boldsymbol{w}_r$を用いて$\boldsymbol{y}_r = \boldsymbol{A} \boldsymbol{w}_r$（$r=1,\cdots,R$）を正規直行化すれば得られる（！）
の部分です。ランクRが小さいほど高速化の恩恵を受けられるので、フルランクを想定してN本の基底を求めに行くのではなく、インクリメンタルに基底を求めていきます。
function orthonormalbases(A::Array{Float64,2}, eta::Float64)
qs = Vector{Vector{Float64}}(undef,0)
N = size(A,2)
for j in 1:N
w = randn(N)
q = A * w
normv = gramschmidt!(q, qs)
if normv < eta
break
end
push!(qs, q)
end
return qs
end
まず正規乱数ベクトルwを生成し、q = A * w を計算しています。グラムシュミットの直行化：
normv = gramschmidt!(q, qs)
で、(j-1)ループ目までに求まっている基底qsで説明される分を差し引いています。このノルムが十分小さかったとき「(j-1)ループ目までで説明しきることができた」＝「Aのランクに達した」と判定してループを抜けています。
gramschmidt!(q, qs)
（懐かしの）グラムシュミットの正規直行化です。
function gramschmidt!(q::AbstractVector{Float64}, qs::AbstractVector{Vector{Float64}})
#Orthogonalization
for q2 in qs
q .-= dot(q,q2) * q2
end
#Normalization
normv = norm(q)
q ./= normv
return normv
end
forループで、すでに求まっている基底qsで説明される分を差し引き、ノルムが1になるように正規化しています。ランクの判定のために正規化のために求めた（正規化前の）ノルムを返しています。
SVDできていることの確認
> using LinearAlgebra, RandomApprox
> A = randn(2,3)
> U, s, V = RandomApprox.svd(A)
> A .-= U * Diagonal(s) * V'
2×3 Array{Float64,2}:
-2.88658e-15 3.55271e-15 5.55112e-15
0.0 3.10862e-15 3.05311e-15
となり、行列Aを再構成できています。 Aが低ランクな場合にもSVDできていることを確認するため、低ランク行列を生成する次の関数を書きました。
function generatelowrankmat(nrow, ncol, rank)
A = zeros(nrow, ncol)
for i in 1:rank
u = randn(nrow)
v = randn(ncol)
s = rand()
A .+= s * (u * v')
end
return A
end
これを用いてもう少し大きい行列で試しても
> A = generatelowrankmat(100, 50, 10)
> U, s, V = RandomApprox.svd(A)
> A .-= U * Diagonal(s) * V'
> maximum(abs.(A))
6.838973831690964e-14
と、Aを再構成できていることが確認できました。
速度比較
フルランクの場合
RandomApprox.svdはAのランクまで圧縮してLinearAlgebra.svdを呼び出しているので、Aがフルランクの場合には低次元化の恩恵は受けられない・・どころか無駄に正規直行基底を算出している分、より遅くなってしまうはずです。確認してみましょう。
> A = randn(1000,500)
> @benchmark U, s, V = LinearAlgebra.svd(A)
211.089 ms (17 allocations: 17.23 MiB)
> @benchmark U, s, V = RandomApprox.svd(A)
499.089 ms (125788 allocations: 993.91 MiB)
1000×500のAがフルランクのとき、通常のLinearAlgebra.svdで211msのところ、RandomApprox.svdでは499ms掛かってしまっています。
低ランクの場合
先ほど定義したgeneratelowrankmatで低ランクのAを生成してベンチマークしてみます。
> A = generatelowrankmat(1000, 500, 10)
> @btime U, s, V = LinearAlgebra.svd(A)
182.693 ms (17 allocations: 17.23 MiB)
> @btime U, s, V = RandomApprox.svd(A)
2.229 ms (108 allocations: 851.08 KiB)
この場合、RandomApprox.svdは2.20msとなり、LinearAlgebra.svdの182msの83倍速になりました!
おわりに
- 『前処理による特異値分解の高速化』をJuliaで実装してみました。
- 分解の対象となる行列が低ランクであれば、かなり高速になりますね！