Python SVDによる低ランク近似

PythonでSVDで入力行列Aのランクをrankに落とすプログラム． ちなみに，C++バージョンはこっち．

import numpy as np
import scipy as sp

def SVD(A, rank):
    u, s, v = np.linalg.svd(A)
    ur = u[:, :rank]
    sr = np.matrix(sp.linalg.diagsvd(s[:rank], rank, rank))
    vr = v[:rank, :]

    Rr = ur*sr*vr

    return Rr

ちなみにこのコードだと，O(n³)かかる．Arpack使うと，O(n² r)でできるはず．．． scipy.sparse.linalg.svds を参考．（スパース行列じゃなくてもできる）

噂では，sklearn.decomposition.truncatedSVD も使えるようだ． 計算量をさらにO(n² log r + r² n)まで落とすこともできるようで，これはこちらの論文に載っている（らしい．私は読んでない）

SVDでrランク近似するの、普通にnumpy.linalg.svd使うとO(n^3)やん。

でも噂では、O(n^2 r)でできるらしいじゃん。どうやるの？Juliaでも良い。既に実装されてるライブラリがあるはずだよね？#Julia言語

— ゆうだい.jl (@physics303) 2020年11月12日