Help us understand the problem. What is going on with this article?

Juliaで『前処理による特異値分解の高速化』

はじめに

以前『前処理による特異値分解の高速化』を書きました。他の言語で実装したことはあったのですが、最近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で実装してみました。
  • 分解の対象となる行列が低ランクであれば、かなり高速になりますね!
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした