音楽プログラミングの超入門(仮)

Python / 音楽情報処理 初心者が、初心者にも分かるような記事を書きたい。

【Python】SciPyの特異値分解ともうちょっと速い特異値分解【SVD】

svd

photo by Jared Tarbell

特異値分解とは?
特異値分解 - Wikipedia

SciPyの特異値分解

SciPyには特異値分解のための関数が2種類入っています。

  1. scipy.linalg.svd
  2. scipy.sparse.linalg.svds

前者(svd)が一般的な特異値分解、後者(svds)がスパース行列に対する特異値分解の関数となっています。しかし、おそらく svds は別の場面で使う場合が多いのではないでしょうか。

svd では全ての特異値を計算するのに対し、svds では上位 k 個の特異値のみを計算させるオプションが存在します。当然、k が小さいほど計算は速く終わるため、スパース行列でなくても、上位の特異値のみが必要な場合は svds を使ったほうが圧倒的に計算が速いことが多いです。

PROPACK

SVD のスピードについて調べているときに下のような記事を見つけました。

Sparse SVDs in Python

svds よりも計算の速い PROPACK なるライブラリが存在するそうです。
PROPACK は pypropack としてPythonライブラリ化されており、以下のリンクからダウンロードし、ビルドすることでインストール可能です。

jakevdp/pypropack · GitHub

比較

pypropack のほうが本当に速いの?みたいなことを調べるため適当にテストしました。本当に適当なので、状況によって結果はかなり変わると思います。
とりあえず、結果的にはちょっと速そうです。

テストコード

# -*- coding: utf-8 -*-
from scipy import rand
from time import time

from scipy.linalg import svd
from scipy.sparse.linalg import svds
from pypropack import svdp

def test_svd():
    shape = [1000, 1000]
    svd_time = []
    svds_time = []
    svdp_time = []

    k = 100 ## 計算したい特異値の数
    for ii in xrange(10):
        M = rand(shape[0], shape[1])
        t1 = time()
        svd(M)
        t2 = time()
        svds(M, k = k)
        t3 = time()
        svdp(M, k = k)
        t4 = time()

        svd_time.append(t2 - t1)
        svds_time.append(t3 - t2)
        svdp_time.append(t4 - t3)

    print "scipy.linalg.svd (full):\t\t", min(svd_time)
    print "scipy.sparse.linalg.svds (k = 100):\t", min(svds_time)
    print "pypropack.svdp (k = 100):\t\t", min(svdp_time)


if __name__ == "__main__":
    test_svd()

出力

scipy.linalg.svd (full):		2.58191680908
scipy.sparse.linalg.svds (k = 100):	1.10373687744
pypropack.svdp (k = 100):		0.83057808876