はじめに

以前『前処理による特異値分解の高速化』を書きました。他の言語で実装したことはあったのですが、最近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で実装してみました。
分解の対象となる行列が低ランクであれば、かなり高速になりますね！