推薦や検索の前処理に使う行列分解は、理論上の計算量だけでは速くならない。疎行列を何回読むか、どちら側に基底を持つか、直交化をどうするか、CPUのメモリ帯域をどう使うかが効いてくる。
そこで、速度優先の近似行列分解ライブラリmymfを実装した。公開クラスは MyMF。呼び出し方は scikit-learn の TruncatedSVD に近く、重い疎行列処理は Cython/C にコンパイルする。GPUは不要である。
v0.2.0では、小規模な実データのMovieLens 100Kでも、単一スレッド・同じ反復回数で TruncatedSVD より高速になった。反復2回で1.64倍、5回で1.88倍である。各手法を300回測定し、推薦指標は10個のseedで評価した。
前の版はMovieLens 100Kで遅かった。大規模疎行列で速いという結果だけでは、小規模でも速いことにはならない。今回は実データを用意し、固定費とCython/Cカーネルを修正した。世界最速、全CPUでの勝利、全推薦モデルより高精度という主張ではない。
何を分解するライブラリなのか
入力をユーザー数 × アイテム数の行列 $X$ とすると、返すものは次の形になる。
$$
Z = XC^{\mathsf T},\qquad \hat X = ZC
$$
components_ が $C$、fit_transform(X) が $Z$ である。$C$ の行はおおむね直交しており、$\hat X$ のランクは指定した n_components 以下になる。
目的は、非中心化行列の低ランク近似である。TruncatedSVD も入力を中心化しないため、疎行列を扱える。[1]
重要なのは、疎行列の未格納要素はゼロとして再構成誤差に含めるという点だ。「観測した評価だけを当てる」「クリックに信頼度の重みを付ける」「ペアのランキング損失を最小化する」モデルではない。観測値だけを扱う行列分解、implicit ALS、BPR、NMFとは区別する必要がある。
推薦にはユーザー・アイテムの埋め込みやスコア計算として利用できるが、再構成誤差が小さいことと、推薦のRecallやNDCGが高いことは別である。そのため、今回は両方を測定した。
手法選び。新しい名前より、対象条件を見る
検討した主な候補は次のとおり。
| 候補 | 今回の位置づけ |
|---|---|
| 通常のrandomized SVD | 比較の基準。実際にインストールされているsklearn 1.8.0で測定する。[1][2] |
| frPCA | 疎行列向けの高速化を扱う先行研究。独立した実装比較は今回は行っていない。[3] |
| dashSVD | 動的シフト付き反復を採用。今回のアルゴリズムの土台。[4] |
| Randomized block Krylov | 収束と保持する基底サイズのトレードオフが異なる候補。未実装・未測定。[5] |
| PRIMME_SVDS、PCAone | 高性能SVD/PCAの比較候補として確認。今回の速度表には含めない。[6][7] |
| Row-aware Randomized SVD | 2024年のプレプリントと2026年の論文掲載を確認。公開年が新しいことだけで採用は決めない。[8] |
選んだのは、2024年の Algorithm 1043: Faster Randomized SVD with Dynamic Shifts、通称dashSVDである。大規模疎行列の、低〜中程度の近似精度を狙う計算に適した有力な研究として採用した。[4]
mymfは論文の固定反復版、Algorithm 4を土台にした独立実装である。著者らのC/MKL実装のラッパーではないし、動的シフト自体を新しい発明として主張するものでもない。追加したのは、Cython/OpenMPの疎行列カーネル、数値的な保護、Python API、推薦補助機能である。
また、著者らのdashSVD/MKL実装と速度を直接比較してはいない。「既存の有力アルゴリズムを採用し、sklearnとの比較で効果を確認した」という範囲の成果である。
まず使ってみる
ソースを展開し、pyproject.toml があるディレクトリでインストールする。PyPIに公開したという意味ではないので、コマンドは pip install mymf ではなく pip install . になる。
Debian系では、CコンパイラとPythonの開発用ヘッダが必要になる。
sudo apt-get install build-essential python3-dev python3-venv
cd mymf
python3 -m venv .venv
. .venv/bin/activate
python -m pip install --upgrade pip
python -m pip install .
利用例は、ランダム行列を作ってshapeを見るだけでは何をしているのか分かりにくい。そこで alice_in_wonder_land.txt をTF-IDF行列にし、MyMF で低次元化して、指定した文に近い文を検索してみる。
ここでは各文が行、単語が列になる。つまり推薦システムでいう「ユーザー × アイテム」の代わりに、「文 × 単語」の疎行列を分解する。同じ低ランク近似を検索用途へ流用した例である。
import re
import numpy as np
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.metrics.pairwise import cosine_similarity
from mymf import MyMF
# 1. Alice in Wonderland を読み込む
with open("alice_in_wonder_land.txt", "r", encoding="utf-8") as f:
text = f.read()
# 改行を空白にしてから、. ! ? を文末として大まかに文分割する。
# 厳密な英文分割器ではないが、サンプル用途ならこれで十分。
text = re.sub(r"\s+", " ", text).strip()
sentences = [
s.strip()
for s in re.split(r"(?<=[.!?])\s+", text)
if len(s.strip()) >= 20
]
# 2. 文 × 単語 のTF-IDF疎行列を作る
vectorizer = TfidfVectorizer(
lowercase=True,
stop_words="english",
min_df=2,
dtype=np.float32,
)
X = vectorizer.fit_transform(sentences)
print("TF-IDF:", X.shape, "nnz=", X.nnz)
# 3. MyMFで低ランク化する
# n_componentsは入力サイズを超えないように調整する。
n_components = min(128, X.shape[0] - 1, X.shape[1] - 1)
if n_components < 1:
raise ValueError("文または語彙が少なすぎて行列分解できません")
mf = MyMF(
n_components=n_components,
n_iter=2,
n_oversamples=10,
random_state=42,
n_jobs=-1,
blas_threads=1,
)
# 各文を n_components 次元の潜在ベクトルへ変換
Z = mf.fit_transform(X)
print("latent:", Z.shape)
# 4. index=200 の文に近い文 Top-k を探す
index = 200
k = 10
if index >= len(sentences):
raise IndexError(
f"index={index} だが、分割後の文数は {len(sentences)} しかありません"
)
# 潜在空間上のcosine similarityを計算
similarities = cosine_similarity(
Z[index:index + 1],
Z,
).ravel()
# 自分自身は除外
similarities[index] = -np.inf
# 全件sortせず、まずTop-k候補だけ拾う
k = min(k, len(sentences) - 1)
topk = np.argpartition(similarities, -k)[-k:]
topk = topk[np.argsort(similarities[topk])[::-1]]
print(f"\nquery [{index}]")
print(sentences[index])
print(f"\nTop {k} similar sentences")
for rank, i in enumerate(topk, 1):
print(f"{rank:2d}. score={similarities[i]:.4f} index={i}")
print(f" {sentences[i]}")
処理の流れは単純で、TfidfVectorizer が作った高次元疎行列 X を MyMF で Z に圧縮し、その潜在空間でcosine similarityを取っている。元のTF-IDF空間で直接類似度を計算することもできるが、文書数や語彙数が大きくなるほど、低次元表現を後段で何度も使い回せる利点が出てくる。
なお、この例は「TF-IDF + SVD系低ランク近似」によるLSA的な検索例であり、Transformer系埋め込みのように語順や深い意味まで理解するものではない。行列分解そのものの使い方を確認するサンプルとして割り切っている。ここで測っているのは低ランク空間上の近さである。
fit、fit_transform、transform、inverse_transform のほか、get_params、set_params、sklearnのクローン、Pipeline、特徴量名、set_output に対応する。
ただし完全な引数互換ではない。algorithm="arpack" や power_iteration_normalizer は用意していない。n_components は両方の入力次元以下に制限する。tol もsklearnと同じ意味ではない。
学習済みモデルで新しい文を処理するときは、同じ TfidfVectorizer で transform してから mf.transform すればよい。語彙の意味と列順は学習時と一致している必要があり、vectorizer.fit_transform を新規文ごとにやり直してはいけない。
動的シフトで何をしているのか
縦長行列の場合は、列側に小さな基底 $Q$ を持たせる。出力ランクを $k$、余分に確保する次元を $s$ とすると、基底の幅は $\ell=k+s$ である。実装では行列の小さいほうの次元までに制限する。
初期化はガウス乱数 $\Omega$ を使って行う。
$$
Q_0 = \operatorname{orth}(X^{\mathsf T}\Omega)
$$
続く反復の中心は、次の計算である。
$$
Y_j = X^{\mathsf T}(XQ_j)-\alpha_jQ_j
$$
大きな $X^{\mathsf T}X$ は作らない。疎行列と細長い密行列の積を2回実行する。
$Y_j$ の小規模SVDから次の基底と特異値を求め、最小の保持特異値に応じてシフト $\alpha$ を更新する。シフトの決め方はdashSVDのものを使っている。[4]
実際の処理を簡略化すると、以下になる。これは説明用の擬似コードで、完全な処理は src/mymf/estimator.py にある。
Q = orth(X.T @ gaussian_matrix)
alpha = 0.0
for _ in range(n_iter):
Y = X.T @ (X @ Q) - alpha * Q
Q, shifted_s, _ = thin_svd(Y)
if shifted_s[-1] > alpha:
alpha = (shifted_s[-1] + alpha) / 2
B = X @ Q
_, singular_values, rotation_t = thin_svd(B)
components = (Q @ rotation_t[:k].T).T
embedding = B @ rotation_t[:k].T
初回のシフトはゼロで、更新されたシフトが使われるのは2回目からである。n_iter=0 はもちろん、n_iter=1 でも動的シフトの効果はまだ出ない。
algorithm="randomized" を指定すると、mymfの他の処理は変えずにシフトだけを無効化できる。アルゴリズムの効果と実装上の効果を切り分けるための比較モードであり、sklearnの実装をそのまま再現するモードではない。
速度のために変更した部分
1. 基底を持つ向きを選ぶ
縦長なら右側、横長なら左側に基底を持つ。小さい側に反復用の基底を置くことで、その直交化の負担を抑える。
ただし、最後に返す値はどちらの向きでも X @ components_.T にそろえる。近似計算の途中で得た U * S を無条件に返すと、fit_transform(X) と transform(X) が違う意味になりうる。
右側の基底を使う場合は、最後の B = X @ Q を再利用して埋め込みを作れる。このため、出力のためだけの追加の疎行列積を省ける。左側の場合は、公開APIとの整合性を優先して最後に投影を計算する。
結果として、学習中の全幅の行列積は右側で 2 * n_iter + 2 回、左側で 2 * n_iter + 3 回になる。ただし、入力検証、転置の準備、分散の計算でも別にデータを走査する。これらまで省かれたという意味ではない。
2. 疎行列積をCython/OpenMPで行う
最も頻繁に使うのは、CSR行列と細長い密行列の積である。各行を別のスレッドに割り当てれば、出力の同じ行を複数のスレッドが更新しないため、アトミック加算が不要になる。
Cython部分では、型付きmemoryview、nogil、prange を使う。[9] 次は実装から抜き出した並列化の中心部分である。
for i in prange(
m,
nogil=True,
num_threads=n_threads,
schedule='guided',
use_threads_if=(data.shape[0] * width >= 200000),
):
_row(data, indices, indptr, b, out, i, width)
非ゼロ数が偏った行にも対応しやすいようにguidedスケジューリングを使う。小さい行列積では並列化の準備のほうが高くつくため、仕事量が少ないときは逐次実行にする。
データはfloat32/float64、インデックスはint32/int64に特殊化する。読み取り専用入力にも対応し、CSRインデックスの検証はPython側の入口で行ってから、内部の境界チェックを省略したカーネルへ渡す。
3. SIMDは実行時に選ぶ
行内の計算はCヘッダ _simd.h と、レジスタ保持用の _reg.h に分離した。対応するx86-64のGCC/Clang環境では、ベースライン、AVX2/FMA、AVX-512/FMAを用意し、CPUの機能を調べて実行時に選ぶ。
切り替え先の関数は非ゼロ要素ごとではなく、行ごとに呼ぶ。係数方向に連続した配列へ加算するため、SIMD化しやすい。
実装を調べたところ、使用したGCCでは外側ループのunroll-and-jam最適化が内側SIMDループをスカラー化するケースがあった。この変換だけをGCCの関数属性で無効化し、他の -O3 最適化は残した。単に「最適化オプションを増やせば速くなる」わけではなかった。
既定ビルドに -march=native や -ffast-math は付けない。前者は配布先CPUの互換性、後者は浮動小数点の扱いを壊しやすいためである。FMAなどによる丸めの差まで完全になくなるわけではない。
AVX-512が常に最速とも限らないので、実行時の切り替えも用意した。環境変数はPythonでimportする前に設定する。
MYMF_SIMD=avx2 python examples/basic.py
MYMF_SIMD=baseline python examples/basic.py
4. 小さな分解は高速経路と安全側を分ける
細長い行列のSVDでは、float64で小さなGram行列を計算して固有値分解する経路を用意した。計算するのはスケッチ幅に応じた小さな行列であって、入力全体の巨大なGram行列ではない。
ただし、この方法だけではランク不足や悪条件の入力に弱い。そこで、最小・最大Gram固有値の比が小さすぎる場合は、Householder QRと小さな行列のSVDに切り替える。ゼロの特異値で割る処理はしない。
normalizer="qr" を指定すれば、最初から安全側の経路を使える。ただし、その前の行列積で失われた精度まで復元できるわけではない。高精度の残差保証が必要な問題を、この速度優先の近似法だけで済ませる意図はない。
また、密行列の積までCythonで手書きの三重ループにしない。その部分はNumPy/SciPyのBLAS/LAPACKを使う。
v0.2で削った固定費と行内の無駄
毎回のライブラリ探索をやめる
旧版では threadpool_limits(...) を呼ぶたびに、読み込まれた共有ライブラリを探索していた。ThreadpoolController を初回だけ作り、各呼び出しでは制限の設定と復元だけを行う。これはthreadpoolctl自身が用意している使い方である。[11]
キャッシュするのはライブラリの検出情報だけだ。行列、転置、乱数、因子、分解結果はキャッシュしない。初回探索の時間は、後述する「新しいPythonプロセスで最初のfit」の測定には含まれる。初回のあとに別のBLASライブラリを読み込んだ場合は、並行実行の外側で mymf.refresh_threadpools() を呼ぶ。
また、n_jobs=1 が明示されたときはCPU数・cgroupの探索をせず、CythonもOpenMP領域に入らない逐次ループを使う。公開API入口でのCSR検証は残し、検証済み行列に対するスカラー倍とSciPyの転置結果を重複して検証することだけを省いた。
行ごとの途中和をレジスタに置く
従来は各非ゼロ要素の処理で、出力行を何度も読み書きしていた。新しい短幅カーネルではSIMDの累積変数に途中和を保持し、CSRの1行を処理し終わってから出力へ書く。端数の幅にはマスク付きのload/storeを使う。例えば幅42でも末尾を越えて読まない。
幅ごとのケースは tools/generate_simd.py で生成する。MovieLensの行列サイズ、データ内容、seedを判定する分岐はない。大きな幅は従来の汎用カーネルへ進む。これは実装の高速化であって、ランクやoversamplingの削減ではない。
乱数と統計量も入力サイズに合わせる
整数seedはNumPyのPCG64へ渡し、float32の作業行列にはfloat32の標準正規乱数を直接作る。[12] 生成する乱数の分布と行列サイズは変えていないが、同じ整数seedでもv0.1とは因子が変わる。旧乱数列を使うには random_state=np.random.RandomState(42) とする。None は従来どおり np.random.seed に従い、v0.2では Generator オブジェクトも受け付ける。
分散はCythonで列ごとのWelford統計量を計算し、未格納のゼロの集団と結合する。非ゼロ値全体をfloat64へコピーする必要がなくなった。explained_variance_ や explained_variance_ratio_ を省略したわけではない。小さなGram行列はLAPACKのDSYEVDを直接呼び、悪条件時のQR/SVDフォールバックも残す。
MovieLens 100Kの測定条件
入力は実際の ml-100k/u.data。[10] 評価4以上を正例とし、各ユーザーの新しい20%をテストへ回す。正例が5件以上あるユーザーを評価する。この分割処理は旧版と同じである。
行列は 943 × 1,682、学習の非ゼロ数は 43,933、テスト正例は 11,442、評価ユーザーは 938。学習済み正例を除いた全アイテムを推薦候補とし、負例サンプリングはしない。両モデルとも共通の transform(query) @ components_ と同じ順位評価コードを使う。
ユーザーごとの時系列分割なので、他ユーザーも含めた全体で未来情報の混入がない分割ではない。また、未観測をゼロとして近似するSVDであって、評価値の欠損を無視する行列分解やimplicit ALSではない。
| 項目 | 条件 |
|---|---|
| Python / NumPy / SciPy | 3.13.5 / 2.3.5 / 1.17.0 |
| scikit-learn | 1.8.0。オンラインの最新版を測ったという意味ではない |
| Cython / OpenBLAS | 3.2.4 / 0.3.30 |
| 型 / ランク / oversampling | float32 / 32 / 10 |
| スレッド | 両者のBLASは1、MyMFの疎行列積も1 |
| ビルド | OpenMPあり、native指定なし、実行時AVX-512/FMA |
| 反復 | q=0, 2, 4, 5。早期停止なし |
| 回数 | seed 0〜9、それぞれ30回。各手法300回の新規fit |
毎回新しいモデルで fit_transform 全体を測る。入力検証、CSR準備、転置、乱数、分解、投影、分散、内部のBLAS制限を含める。import、データ読込、モデルの構築、ウォームアップ、精度評価は除く。各測定ラウンドで実行順を変え、遅い回を削除せず全記録を保存した。
時間は300回をまとめた中央値、推薦指標は10個のseedの平均である。同じseedの30回を、独立した精度評価30件として数えてはいない。
MovieLens 100Kの結果
| 手法 | 学習時間中央値［ms］ | Recall@20 | NDCG@20 | 相対Frobenius誤差 |
|---|---|---|---|---|
| TruncatedSVD_q2 | 5.647 | 0.195720 | 0.158065 | 0.734004 |
| MyMF_q2 | 3.451 | 0.195770 | 0.158434 | 0.731105 |
| TruncatedSVD_q5 | 9.448 | 0.197715 | 0.160791 | 0.728483 |
| MyMF_q5 | 5.018 | 0.198344 | 0.161064 | 0.727808 |
| MyMF_q4 | 4.512 | 0.197540 | 0.160583 | 0.728152 |
| MyMF_q0 | 2.375 | 0.177086 | 0.143387 | 0.759070 |
同じq=2では 1.64倍、同じq=5では 1.88倍高速だった。今回の平均Recall/NDCGはsklearnを下回らなかったが、この小さな差から推薦精度の統計的優越までは主張しない。
q=0 は品質を落としている。速度比較の主張はこれではなく、同じ反復回数のq=2同士・q=5同士に基づく。
旧乱数列を使う RandomState 指定でも、q=2は 3.955ms、q=5は 5.535msで、対応するsklearnより速かった。こちらはRandomStateオブジェクトの作成をfit時間に含めない補助測定であり、旧版の実行時間そのものではない。
AVX2を指定した追加測定でも同じ反復数で速度優位を確認した。これは同じCPUで命令セットを制限したテストであり、古い物理CPUでも同じ倍率になることの証明ではない。
初回だけ遅い実装になっていないか
各測定で新しいPythonプロセスを起動し、MyMFのウォームアップなしで最初の fit_transform を測った。各手法5プロセスの中央値は次のとおり。初回の内部コントローラ探索も含むが、importとデータ読込は両者とも含まない。
| 手法 | 初回fitの中央値［ms］ |
|---|---|
| sklearn_q2 | 6.673 |
| sklearn_q5 | 10.903 |
| mymf_q2 | 4.922 |
| mymf_q5 | 6.370 |
今回の環境では初回もsklearnより速かった。起動時間や全ハードウェアでの性能まで保証する結果ではない。
大規模疎行列でも再測定する
小規模だけを速くして大規模を壊していないかも確認した。旧版と同じ生成規則の 50,000 × 20,000、非ゼロ2,444,696個の合成疎行列で、float32、ランク64、oversampling 10、3個のseed、1スレッドで測った。
| 手法 | 時間の中央値［秒］ | 相対再構成誤差 |
|---|---|---|
| sklearn_q2 | 0.447861 | 0.97535604 |
| mymf_q2 | 0.219736 | 0.97534519 |
| sklearn_q5 | 0.767076 | 0.97519049 |
| mymf_q4 | 0.319305 | 0.97517225 |
| mymf_q5 | 0.349339 | 0.97513339 |
ここでも速度優位は残った。誤差が0.975というのは「精度97.5%」ではない。ノイズの多い行列に対する低ランク近似の残差を比較している。
旧版の密行列、4スレッド、ピークRSS、合成推薦の測定は docs/benchmarks_v0.1.md と元のJSONへ履歴として残した。それらの値はv0.2の測定ではない。今回、密行列全般やメモリ使用量の改善は主張しない。
推薦時も、全スコア行列は作らない
実際の呼び出しには recommend を用意した。
item_ids, scores = model.recommend(
X[:100],
n_recommendations=20,
filter_seen=True,
batch_size=32,
item_batch_size=4096,
)
ユーザー側とアイテム側を両方ブロックに分け、各ブロックの上位候補を更新する。全ユーザー × 全アイテムの密スコア行列を一度に確保しない。ただし全候補の内積は計算するので、近似近傍探索のように候補を飛ばす方式ではない。
既に入力に非ゼロ値があるアイテムを除外し、候補が足りなければIDを -1、スコアを -inf にする。IDは元の行列の0始まりの列番号である。同点の境界では選ばれるアイテムが変わりうる。
入力が全ゼロの新規ユーザーには、個人の好みを推定する情報がない。この場合の同点ランキングを、意味のあるコールドスタート対策と考えるべきではない。必要ならアプリケーション側で人気順などの方針を用意する。
パラメータの選び方
速度優先なら、まずfloat32のCSR、n_iter=2、n_oversamples=10 から始める。n_iter=0、2、4を検証用データで比較し、精度が必要なら反復数や幅を増やす。q=0を無条件で使うと、今回のように大きく品質を落とすことがある。
難しい数値条件には dtype="float64", normalizer="qr" を使える。ただし処理時間は増える。tol を正にすると早期停止できるが、保証されたSVD誤差のしきい値ではなく、変化量によるヒューリスティックである。掲載した測定はすべて tol=0 で行った。
n_jobs は疎行列積、blas_threads はBLASを制御する。小さな行列を大量に分解するときに両方を大きくすると、スレッドの奪い合いで遅くなることがある。疎行列ではBLASを1に固定し、OpenMPを1、2、4と試す。密行列ではBLASのスレッド数を別に調整する。
また、threadpoolctlによるBLAS制限はプロセス全体へ一時的に作用する。同一プロセスのPythonスレッドから複数モデルを同時実行する場合は、外側でBLASの制限を設定し、blas_threads=None を使うほうが扱いやすい。
コンパイルと検証
この実装は .pyx を置いただけではない。CythonからCを生成してコンパイルし、Pythonから拡張モジュールを読み込むことを実際に確認した。
Linux CPython 3.13.5で 291件のテストが通った。scikit-learnの check_estimator もテストに含めている。疎・密行列、float32/float64、32/64ビットの疎行列インデックス、ゼロ行列、ランク不足、極端なスケール、読取専用配列、不正インデックス、保存と再読込、Pipeline、DataFrameの列名検証などを確認した。
さらに、各反復でLAPACKの通常SVDを使う別の小規模参照実装をテスト内に置き、特異値、シフト履歴、部分空間を比較した。成分の符号が一致しない場合でも比較できるよう、部分空間は射影行列で検証している。
OpenMPを無効にした別ビルドでも291件が通った。実行時SIMDの自動選択、AVX2指定、ベースライン指定でもテストした。 ソース配布物からwheelを作り、別の仮想環境にインストールした状態でも291件が通った。後者は既存の実行時依存を共有したオフライン検証であり、依存関係を新規ダウンロードする検証ではない。
python -m pip install -e '.[dev]'
python -m pytest -q
python examples/basic.py
python examples/recommendation.py
python examples/pipeline.py
テストのpandas項目はpandas未導入ならスキップされる。CIではpandasも入れる設定にしている。Python 3.10〜3.13のLinuxと、macOS/Windows向けCI定義も用意したが、これらのリモートジョブをすべて実行済みという意味ではない。
OpenMPの使えない環境では、逐次版を明示的にビルドできる。
MYMF_OPENMP=0 python -m pip install .
ビルド依存が導入済みで、一時領域やネットワークに制約がある場合は、次のようにする。
mkdir -p .tmp
TMPDIR="$PWD/.tmp" python -m pip install \
--no-cache-dir --no-build-isolation --no-deps .
この指定は足りない依存関係を自動でそろえない。必要なNumPy、SciPy、scikit-learn、Cythonなどはあらかじめ導入する。通常のisolated build用の設定は pyproject.toml にあるが、この環境ではネットワーク経由の依存取得を含めて検証したわけではない。
同じ測定を再実行する
python benchmarks/recommendation.py \
--movielens /path/to/ml-100k/u.data \
--threads 1 --rank 32 --oversamples 10 \
--repeat 10 --timing-repeats 30 \
--include-legacy --assert-faster \
--output benchmarks/results/local_movielens.json
python benchmarks/cold_start.py \
--movielens /path/to/ml-100k/u.data --repeat 5 \
--output benchmarks/results/local_cold.json
python benchmarks/run.py \
--case sparse_large --threads 1 --repeat 3 \
--output benchmarks/results/local_large_t1.json
--assert-faster は、q=2またはq=5でsklearnより遅ければ、JSONを保存してから非ゼロ終了する。成功するまで測り直したり、失敗した測定値を捨てたりはしない。CPU、BLAS、コンパイラ、データ、ランクを結果と一緒に確認できる。
MovieLens本体はZIPへ再配布せず、手元の公式 u.data を使う。記事のAliceサンプルは、手元の alice_in_wonder_land.txt を使う。
結論
行列分解の高速化は、反復回数を削るだけでも、全部をCythonに書き換えるだけでも不十分だった。今回効いたのは、動的シフトを使う既存アルゴリズムを選び、疎行列積とデータの流れをCPU向けに実装し、品質と一緒に測ることである。
mymfは、非中心化の低ランク近似が欲しい疎行列ワークロードで、速度と精度の調整幅を持たせた実装になった。一方、MovieLens以外への一般化、高精度のSVD、全SOTA実装との比較は別の検証課題として残る。何を確認できたかと、何をまだ確認していないかを分けて利用してほしい。
参考資料
-
scikit-learn: TruncatedSVD。APIの比較資料。実測バージョンは1.8.0。
-
Halko, Martinsson, Tropp. Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions, 2011.
-
Feng et al. Fast Randomized PCA for Sparse Data, 2018.
-
Feng, Yu, Xie, Tang. Algorithm 1043: Faster Randomized SVD with Dynamic Shifts, ACM TOMS, 2024. プレプリント、著者実装。
-
Musco and Musco. Randomized Block Krylov Methods for Stronger and Faster Approximate Singular Value Decomposition, 2015.
-
GroupLens: MovieLens 100K Dataset。v0.2の実データ評価に使用。
-
threadpoolctl公式README。ThreadpoolControllerの再利用と、作成後にロードしたライブラリの扱い。
-
NumPy Generator.standard_normal。float32/float64を指定した標準正規乱数生成。