どの程度sparseだと疎行列ライブラリで計算が速くなるのか

気になったので試してみました。

TL; DR

Objective

「疎行列クラスの計算速度向上効果は、データがどの程度疎（sparse）になったら発揮されるのか」について調べた。

Methods & Data

0か1の要素からなるデータを作り、その計算時間を

1. 和の演算
2. 積の演算
3. XGBoostの学習
4. LightGBMの学習

の4つのタスクについて計測した。

Findings

1. 行列の和・積の演算は、非ゼロ要素の割合がどの程度であっても（1%でも100%でも）np.arrayよりcsr_matrixやcsc_matrixのほうが高速であった。
2. XGBoostは非ゼロ要素の割合が約90％以下であれば疎行列クラスの使用による学習時間向上がみられた。
3. LightGBMは非ゼロ要素の割合がどの程度であっても、疎行列クラスによる学習時間の向上はみられなかった。
4. データがsparseでない状況で疎行列クラスを使用することで計算速度が低下するようなことはなかった。

なお、今回使用したコードです。

github.com

はじめに

疎行列とは

疎行列（sparse matrix）とは、成分のほとんどがゼロである行列のこと。

この疎行列のようなゼロの多いデータを効率よく格納するためのクラスがscipyなどのライブラリに実装されています。

今回試すこと

疎行列クラスでデータを保持すると、非ゼロ要素のデータしか持たないため、メモリ使用量の観点では極めて効率がよくなります。

また、計算速度も高くなると言われています。

今回は「どの程度非ゼロ要素が少なくなると疎行列ライブラリによる計算速度向上の効果が出るのか」について試してみます。

方法

環境

• Google Colaboratory
• Intel(R) Xeon(R) CPU @ 2.30GHz
• Python 3.6.9
• numpy 1.17.4
• scipy 1.3.3

データ

0.0か1.0が確率pで現れる二項分布を使って任意の密度の行列を作ります。

def gen_dataset(nrow, ncol, p):
    x = np.random.binomial(n=1, p=p, size=nrow*ncol)
    X = x.reshape(nrow, ncol).astype(np.float16)
    return X

行列クラス

以下の4つを比較します。

1. numpy.array()

疎行列ではない行列クラス

2. scipy.sparse.lil_matrix

疎行列。リストの中に非ゼロ要素の列インデックスと値を保持したリストを入れて保持する（List of lists: LIL）方式。

3. scipy.sparse.csr_matrix

圧縮行格納方式（Compressed Sparse Row: CSR）といって、非ゼロ要素の行インデックス情報を圧縮して保持する。

これについてはWikipediaの説明がわかりやすい（疎行列 - Wikipedia）

4. scipy.sparse.csc_matrix

圧縮列格納方式（Compressed Sparse Column）。CSRの列バージョン。

検証１：和の演算

方法

• 同じ大きさの正方行列を作ってA + Bします。
• 演算を3回繰り返して平均処理時間で評価します。
• 非ゼロ要素の割合を[0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]に変化させつつ演算の処理時間を計測していきます。
• 行列のサイズに依存して結果が変わるかもしれないので、行数・列数を[100, 1000, 10000]に変えて試します。

import pandas as pd
import timeit
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix, csc_matrix

def gen_dataset(nrow=100000, ncol=100, p=0.5):
    x = np.random.binomial(n=1, p=p, size=nrow*ncol)
    X = x.reshape(nrow, ncol).astype(np.float16)
    return X

results = []
for nrow in [100, 1000, 10000]:
    ncol = nrow
    for p in [0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]:
        print(f"# nrow={nrow}, ncol={ncol}, p={p} "+"-"*50)
        A_np = gen_dataset(nrow=nrow, ncol=ncol, p=p)
        B_np = gen_dataset(nrow=nrow, ncol=ncol, p=p)
        for matrix in ["np", "lil", "csr", "csc"]:
            if matrix == "lil":
                A = lil_matrix(A_np)
                B = lil_matrix(B_np)
            elif matrix == "csr":
                A = csr_matrix(A_np)
                B = csr_matrix(B_np)
            elif matrix == "csc":
                A = csc_matrix(A_np)
                B = csc_matrix(B_np)
            else:
                A = A_np
                B = B_np
            t = timeit.Timer('A + B', globals=globals())
            fit_times = t.repeat(repeat=3, number=1)
            result = {"nrow": nrow, "ncol": ncol, "p": p, "matrix": matrix, 
                            "mean": np.mean(fit_times), "std": np.std(fit_times)}
            results.append(result)
results = pd.DataFrame(results)

結果

lil_matrixは公式ドキュメントでも「lil + lilの演算は遅い」と書かれているだけあって、遅いです。

lilを除いてグラフを再描画するとこうなります

そこそこの大きさの行列になってくると、CSC・CSRがnp.arrayと大差をつけて非常に速いです。

そしてこの関係は密な（非ゼロ要素が多い）行列でも変わりませんでした。

疎行列でなくてもCSR・CSCのほうがnp.arrayよりも高速に演算できるみたいです。（なんで？）

なお、CSRとCSCで比べても大差はありませんでした。

検証２：積の演算

方法

• 同じ大きさの正方行列を作ってA.dot(B)します。
• 演算を3回繰り返して平均処理時間で評価します。
• 非ゼロ要素の割合を[0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]に変化させつつ演算の処理時間を計測していきます。
• 行列のサイズに依存して結果が変わるかもしれないので、行数・列数を[100, 1000]に変えて試します。

結果

CSC = CSR < LIL << np.arrayという感じの順番になりました。

こちらもCSRとCSCの結果は同程度でした。

検証３：XGBoostの学習

方法

• fit()を10回繰り返して平均処理時間で評価します。
• 非ゼロ要素の割合を[0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]に変化させつつ演算の処理時間を計測していきます。
• 行列のサイズに依存して結果が変わるかもしれないので、行数を[100, 1000, 10000, 100000]と変化させ、列数は[10, 100]と変化させてすべての組み合わせについて試してみます。
• lil_matrixには対応していないので、np.array, csr_matrix, csc_matrixの３種で比較します。

import pandas as pd
import timeit
import xgboost as xgb
import lightgbm as lgb
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix, csc_matrix


def gen_dataset(nrow=100000, ncol=100, p=0.5):
    x = np.random.binomial(n=1, p=p, size=nrow*ncol)
    X = x.reshape(nrow, ncol).astype(np.float16)
    y = np.apply_along_axis(_func, 1, X)
    return X, y


def _func(row: np.array):
    return np.sum(row) + np.random.rand()

results = []
for nrow in [100, 1000, 10000, 100000]:
    for ncol in [10, 100]:
        for p in [0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]:
            print(f"# nrow={nrow}, ncol={ncol}, p={p} "+"-"*50)
            X_np, y = gen_dataset(nrow=nrow, ncol=ncol, p=p)
            for algorithm in ["xgb", "lgb"]:
                model = xgb.XGBRegressor(objective="reg:squarederror") if algorithm == "xgb" else lgb.LGBMRegressor()
                for matrix in ["np", "csr", "csc"]:
                    if matrix == "csr":
                        X = csr_matrix(X_np)
                    elif matrix == "csc":
                        X = csc_matrix(X_np)
                    else:
                        X = X_np
                    t = timeit.Timer('model.fit(X, y)', globals=globals())
                    fit_times = t.repeat(repeat=3, number=1)
                    result = {"nrow": nrow, "ncol": ncol, "p": p, "matrix": matrix, "algorithm": algorithm,
                              "mean": np.mean(fit_times), "std": np.std(fit_times)}
                    results.append(result)
results = pd.DataFrame(results)

結果

非ゼロ要素の割合が少なくなる（よりsparseになる）ほど、CSR・CSCのほうがnp.arrayよりもずっと高速に計算できるようです。

そして、非ゼロ要素の割合が90%くらいであっても、CSC・CSRを使ったほうが高速であるようです。

検証４：LightGBMの学習

方法

XGBoostと同じです。

結果

行数・列数、非ゼロ要素の割合を変えてみても、np.array、csr、cscは同程度の計算時間であり、疎行列クラスの計算速度向上はみられませんでした。

なぜなのかはわかりません・・・

おわりに

まとめると

1. CSR・CSCはnp.arrayより和・積の演算が高速
2. XGBoostは非ゼロ要素が約90%以下であればCSC・CSRを使ったほうがいい
3. LightGBMは疎行列クラスを使っても計算速度の向上は期待できなさそう¹
4. ゼロ要素の少ないデータ（密なデータ）であっても、疎行列クラスを使って計算が遅くなることはない

という感じでしょうか。

CSR・CSCのデメリットは、scipyのドキュメントを見ると、圧縮している行/列方向へのスライス操作が遅くなることや構造を変更することであるようなので、計算上はメリットが多くてデメリットはほぼ無いんじゃないか、という印象です。

前処理では使わないほうがいいでしょうが、機械学習アルゴリズムに投入する段階になったらCSR・CSCに変換してから突っ込めば疎行列クラスの恩恵をうまく受けられるかもしれません。

1. データに10％以上ゼロ要素がある
2. 使おうとしている機械学習アルゴリズムが疎行列クラスを受け入れている（疎行列用の高効率な演算を実装している）

という状況であれば疎行列クラスを試してみるのはアリだと思います。

---

1. とはいえ、sparseなデータであれば疎行列クラスのほうがメモリ効率は良いはずなので「どんなに疎なデータであってもnp.arrayを使えばいい」というわけではない↩