Just another Frontier

scipyにはscipy.sparseという疎行列の生成を扱うパッケージとscipy.sparse.linalgという疎行列の線形代数を扱うパッケージが存在する。しかし、この2つのリファレンスをいくら読んでも粗行列の行列積については触れていないので、公式ドキュメントと以下のスレッドを元に補足してみた。
[SciPy-User] Sparse matrices and dot product

sparse_matrixの扱いかた

scipy.sparseにはいくつかのフォーマットがあるが、基本的には生成時にlil_matrixを使用し、それをcsr_matrixに変換して計算するというのが効率的なようだ。

import scipy.sparse as sp

# generate 1000x1000 sparse matrix, whose components are initialized with 0
A_lil = sp.lil_matrix((1000,1000))

# do some thing to A_lil
do some thing to A_lil

# convert to csr_matrix
A_csr = A_lil.tocsr()

# do main calculation to A_csr
do some thing to A_csr

scipy.sparaseのsparsematrixはその存在自体が謎なnumpy.ndarrayのサブクラスの”matrix”ライクに操作できる。
すなわち”*”演算子やnp.dotを用いて行列積を計算する。

B_csr = A_csr.T * A_csr
B_csr2 = dot(A_csr.T,A_csr)

しかし、sparsematrixとndarrayに対してnp.dotを使うと変なことになる。

# this is invalid, C is 1D-ndarray of dtype=object
C = np.dot(A_csr.toarray(), A_csr)

この場合でも”*”演算子はちゃんと機能する。

# this is valid C2 is 2D-ndarray of dtype=np.float64
C2 = A_csr.toarray() * A_csr

そういったわけで、scikits.learn.util.extmathにはsafe_sparse_dotと言う関数が用意されているのでこれを使うのを推奨します。

def safe_sparse_dot(a, b, dense_output=False):
    """Dot product that handle the sparse matrix case correctly"""
    from scipy import sparse
    if sparse.issparse(a) or sparse.issparse(b):
        ret = a * b
        if dense_output and hasattr(ret, "toarray"):
            ret = ret.toarray()
        return ret
    else:
        return np.dot(a,b)

行列積の速度

以下のコードを走らせてみるとやはり行列積の速度は
csr_matrix > lil_matrix > ndarray
の順になった。というかcsr_matrixはかなり速い!もちろん行列が密になってくると速度の差はどんどん縮まっていきますのでご注意。

5000×5000成分のうち、90%が0の行列の行列積の結果

lil format needed 9.778 sec
csr format needed 5.387 sec
dense ndarray format needed 26.733 sec

error of lil = 1.34506066158e-22
error of csr = 1.34506066158e-22

テストコード

# -*- coding: utf-8 -*-
"""
Created on Tue Aug  2 11:04:57 2011

@author: lucidfrontier45
"""

from time import time
import numpy as np
import scipy.sparse as sp

def prepare_matrices(l=1000):
    # indexes and values of non-zero components
    print "generating indexes and values"
    ii = np.random.randint(0,l,(int(l*l/10),2))
    r = np.random.rand(int(l*l/10))

    # create a lil sparse matrix
    print "making lil matrix"
    A_lil = sp.lil_matrix((l,l))
    for n,i in enumerate(ii):
        A_lil[i] = r[n]
    A_lil.setdiag(np.ones(l))

    # convert to csr sparse matrix
    print "converting lil -> csr"
    A_csr = A_lil.tocsr()

    # convert sparse matrix to dense matrix
    print "converting lil -> dense ndarray"
    A_dense = A_lil.toarray()

    return A_lil, A_csr, A_dense

def tests_sparse_matmul(A_lil,A_csr,A_dense):
    """
    compare speed of matrix product
    """

    # mesure speed of matrix product for each format
    print "######################################"
    print "computing A_lil^T * A_lil"
    t1 = time()
    A2_lil = A_lil.T * A_lil
    t2 = time()
    print "lil format needed %8.3f sec" %(t2-t1)

    print "######################################"
    print "computing A_csr^T * A_csr"
    t1 = time()
    A2_csr = A_csr.T * A_csr
    t2 = time()
    print "csr format needed %8.3f sec" %(t2-t1)

    print "######################################"
    print "computing A_dense^T * A_dense"
    t1 = time()
    A2_dense = np.dot(A_dense.T,A_dense)
    t2 = time()
    print "dense ndarray format needed %8.3f sec" %(t2-t1)

    print "######################################"
    print "error of lil = ",((A2_lil.toarray()-A2_dense)**2).sum()
    print "error of csr = ",((A2_csr.toarray()-A2_dense)**2).sum()

if __name__ == "__main__":
    from sys import argv
    A_lil, A_csr, A_dense = prepare_matrices(int(argv[1]))
    tests_sparse_matmul(A_lil,A_csr,A_dense)
広告

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中

タグクラウド

%d人のブロガーが「いいね」をつけました。