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

# はじめに
以前[『前処理による特異値分解の高速化』](https://qiita.com/msekino/items/3ec058e8b363e3f396e7)を書きました。他の言語で実装したことはあったのですが、最近Juliaの勉強を始めたので試しにJuliaで実装してみました。

# 通常の特異値分解
LinearAlgebraパッケージにSVDがあります。randn(2,3)で2×3の正規乱数行列を生成してSVDしてみます。

```Julia
> 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を再構成できるか確認してみましょう。

```Julia
> 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パッケージを用います。

```Julia
> 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への実装
## 『前処理による特異値分解の高速化』の復習
[『前処理による特異値分解の高速化』](https://qiita.com/msekino/items/3ec058e8b363e3f396e7)を復習すると、

- 行列$\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への実装
次のように実装してみました。

```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が求まる（！）

の部分です。

```Julia
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が横長だったら転置して縦長にしています。


```Julia
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本の基底を求めに行くのではなく、インクリメンタルに基底を求めていきます。

```Julia
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 を計算しています。グラムシュミットの直行化：

```Julia
normv = gramschmidt!(q, qs)
```
で、(j-1)ループ目までに求まっている基底qsで説明される分を差し引いています。このノルムが十分小さかったとき「(j-1)ループ目までで説明しきることができた」＝「Aのランクに達した」と判定してループを抜けています。

### gramschmidt!(q, qs)
（懐かしの）グラムシュミットの正規直行化です。

```Julia
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できていることの確認

```Julia
> 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できていることを確認するため、低ランク行列を生成する次の関数を書きました。

```Julia
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
```
これを用いてもう少し大きい行列で試しても

```Julia
> 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がフルランクの場合には低次元化の恩恵は受けられない・・どころか無駄に正規直行基底を算出している分、より遅くなってしまうはずです。確認してみましょう。

```Julia
> 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を生成してベンチマークしてみます。

```Julia
> 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で実装してみました。
- 分解の対象となる行列が低ランクであれば、かなり高速になりますね！