Julia言語 Nyström method による行列の低ランク近似

ようやく実装できた．参考にした論文は以下の２つ．

proceedings.neurips.cc

ojs.aaai.org

それにしても，こんなに著名な手法のjulia実装でどこにもあがってないのはちょっと驚き． サンプリングの方法はいろいろあるんだけど，最も古典的（？）と思われるランダムサンプリングにしている．

using LinearAlgebra
using StatsBase

"""
low-rank approximaion based on random sampling.
See the following reference,
[Petros Drineas et al.](https://www.jmlr.org/papers/v6/drineas05a.html)
# Aruguments
- 'K' : input positive definite matrix.
- 'p' : target rank

example:
    A = randn(100,100)
    K = A*A'
    for p in [1,5,10,20,30,50,100]
        Kr = nystrom(K, p, reconst=true)
        @show p, norm(K-Kr)
    end
"""

function nystrom(K,p;reconst=false)
    n, n = size(K)
    idxs = sample(1:n, p, replace=false, ordered=false)
    C = K[:,idxs]
    W = K[idxs,idxs]

    if reconst == true
        return C*( W \ C')
    else
        return C, inv(W)
    end
end

テストがてら，適当な正定値行列Kをつくって，ランクp近似をしてみる． pを大きくしていくにつれて，再構成誤差がさがるはず．．．

A = randn(100,100)
K = A*A'
for p in [1,5,10,20,30,50,100]
    Kr = nystrom(K, p, reconst=true)
    @show p, norm(K-Kr)
end

#(p, norm(K - Kr)) = (1, 1383.25284972177)
#(p, norm(K - Kr)) = (5, 1296.5715363623078)
#(p, norm(K - Kr)) = (10, 1211.8022048540815)
#(p, norm(K - Kr)) = (20, 1005.0471242087337)
#(p, norm(K - Kr)) = (30, 812.8200435737806)
#(p, norm(K - Kr)) = (50, 506.67283549630315)
#(p, norm(K - Kr)) = (100, 1.4370668157705365e-8)

pを大きくすると，どんどん再構成誤差が小さくなっている．問題なく実装できてそう．