Julia言語 逆行列の高速な計算のために，右除算演算子と左除算演算子を利用しよう．

数値計算で逆行列の計算を避けましょうというのはもはや常識と思います．例えば，givenの行列AとBについて，C = inv(A) * B を計算したいときは，左除算演算子\を使って

using LinearAlgebra
C = A \ B

とした方が速いです．実際に

using BenchmarkTools

n = 6000
A = rand(n,n)
B = rand(n,n)
@btime inv($A) * $B
# 3.458 s (8 allocations: 552.29 MiB)
@btime $A \ $B 
# 2.206 s (6 allocations: 552.29 MiB)

念のため，inv(A) * B と A \ B がほとんど同じ出力であることを確認しておきます．

norm( inv(A) * B - A \ B)
#2.2633599357010303e-9

一方で，C = B * inv(A)を高速に計算したいときはどうしたらよいのでしょうか．こういう時は，実は右除算演算子/を使えばよいらしいです．知らなかったので，早速比較実験してみると

@btime $B * inv($A)
# 3.474 s (8 allocations:980.53 MiB)
@btime $B / $A 
# 2.717 s (8 allocations: 976.62 MiB)

有意差があるように思います．念のため B * inv(A) と B/A が同じことを確認します．

norm( B * inv(A) - B/A)
#4.5818497409971097e-8

一緒ですね！

ということで，右除算演算子と左除算演算子を駆使して，今後はinvは使わないようにしましょう～．

追記

最近書いてた，Randomizedな固有値分解の実装にinvが含まれていたので，改めてinvが含まれていない形に書き直しておいた．

using LinearAlgebra

"""
Randomized eigenvalue decomposition via Nystrom Method.
See Fig 6 in the [original paper](https://arxiv.org/abs/1607.01649)
or See Algorithm 5.5 in the [original paper](https://epubs.siam.org/doi/10.1137/090771806)
This method is available for only a positive definite matrix.

# Aruguments
- 'A' : input positive definite matrix
- 'k' : target rank, number of wanted eigenvalues
- 'p' : over-sampling parameter, it should be 10 or 2k.
- 'reconst' : if ture, the output is reconstracted matrix which approximates 'A'

example:
    # make a positive definite matrix `A`.
    n = 5000; m = 25000; x = randn(m,n)
    demean_x = x .- mean(x, dims=1)
    A = demean_x'*demean_x ./(m-1)

    # run evdnys
    lowrank_mtx = evdnys(A,10, reconst=true)
"""
function evdnys(A::Symmetric, k::Int;p=10, reconst=false)
    n, n = size(A)
    G = randn(n, k+p)
    Y = A * G
    Q = qr(Y).Q*Matrix(I,n,k+p)
    B1 = A*Q
    B2 = Q'*B1
    C = cholesky(Symmetric(B2))
    #F = B1 * inv(C.U)
    F = B1 / C.U
    U, Σ, _ = svd(F)

    if reconst
        return U*diagm(Σ.^2)*U'
    else
        return U, Σ.^2
    end
end

"""
Randomized eigenvalue decomposition proposed by N.halko et al.
See Algorithm 5.6. in the
[original paper](https://epubs.siam.org/doi/10.1137/090771806).
This method is effective matrices with decaying eigenvalues.

# Aruguments
- 'A' : input positive definite matrix
- 'r' : target rank, number of wanted eigenvalues
- 'p' : over-sampling parameter, it should be 10 or 2r.
- 'reconst' : if ture, the output is reconstracted matrix which approximates 'A'

example:
    # make a matrix with decaying eigenvalues
    n = 150; A = rand(n,n); B = zeros(n,n)
    e, v = eigen(A)
    for i in 1:n
        B += e[i] * v[:,i] * v[:,i]' * (i>n-100 ? 1 : 1e-3)
    end

    # run reigen
    lowrank_mtx = reigen(A, 10, reconst=true)
"""
function reigen(A, r ;p=10, twopass=false, reconst=false)
    n = size(A)[1]
    m = r + p

    Ω = randn(n, m)
    Y = A*Ω

    
    if twopass
        Q = qr!(Y).Q * Matrix(I,n,m)
        T = Q' * A * Q
    else
        Q = qr(Y).Q * Matrix(I,n,m)
        T = (Q' * Y) / ( Q' * Ω)
    end
    λ, S = eigen(T)
    U = Q*S'
    if reconst
        Λ = Diagonal(λ)
        return U*Λ*U'
    else
        return U, λ
    end
end