18

私は次の方程式を実装するために働いています：

X =(Y.T * Y + Y.T * C * Y) ^ -1

Yは（nxf）行列で、Cは（nxn）対角行列です。nは約300kで、fは100から200の間で変化します。最適化プロセスの一部として、この方程式はほぼ1億回使用されるため、非常に高速に処理する必要があります。

Yはランダムに初期化され、Cは非常にスパースな行列であり、対角線上の300kのうち数個だけが0とは異なります。Numpyの対角関数は密な行列を作成するため、Cをスパースなcsr行列として作成しました。しかし、方程式の最初の部分を解こうとすると、次のようになります。

r = dot(C, Y)

メモリ制限が原因でコンピュータがクラッシュします。次に、Yをcsr_matrixに変換して、同じ操作を行うことにしました。

r = dot(C, Ysparse)

このアプローチには1.38ミリ秒かかりました。しかし、スパース行列を使用して密な行列を格納しているため、このソリューションはやや「トリッキー」です。これが実際にどれほど効率的か疑問に思います。

だから私の質問は、Yをスパースに変えてパフォーマンスを向上させることなく、スパースCとデンスYを乗算する方法があるかどうかです。どういうわけか、大量のメモリを消費せずにCを対角線密度で表すことができれば、これは非常に効率的なパフォーマンスにつながる可能性がありますが、これが可能かどうかはわかりません。

私はあなたの助けに感謝します！

4

4 に答える 4

29

r = dot（C、Y）を計算するときにドット積がメモリの問題に遭遇する理由は、numpyのドット関数がスパース行列を処理するためのネイティブサポートを持っていないためです。何が起こっているのかというと、numpyはスパース行列Cをnumpy配列ではなく、pythonオブジェクトと見なしています。小規模で検査すると、問題を直接確認できます。

>>> from numpy import dot, array
>>> from scipy import sparse
>>> Y = array([[1,2],[3,4]])
>>> C = sparse.csr_matrix(array([[1,0], [0,2]]))
>>> dot(C,Y)
array([[  (0, 0)    1
  (1, 1)    2,   (0, 0) 2
  (1, 1)    4],
  [  (0, 0) 3
  (1, 1)    6,   (0, 0) 4
  (1, 1)    8]], dtype=object)

明らかに、上記はあなたが興味を持っている結果ではありません。代わりに、あなたがしたいのは、scipyのsparse.csr_matrix.dot関数を使用して計算することです。

r = sparse.csr_matrix.dot(C, Y)

以上コンパクトに

r = C.dot(Y)
于 2013-05-25T22:26:38.127 に答える
9

試す：

import numpy as np
from scipy import sparse

f = 100
n = 300000

Y = np.random.rand(n, f)
Cdiag = np.random.rand(n) # diagonal of C
Cdiag[np.random.rand(n) < 0.99] = 0

# Compute Y.T * C * Y, skipping zero elements
mask = np.flatnonzero(Cdiag)
Cskip = Cdiag[mask]

def ytcy_fast(Y):
    Yskip = Y[mask,:]
    CY = Cskip[:,None] * Yskip  # broadcasting
    return Yskip.T.dot(CY)

%timeit ytcy_fast(Y)

# For comparison: all-sparse matrices
C_sparse = sparse.spdiags([Cdiag], [0], n, n)
Y_sparse = sparse.csr_matrix(Y)
%timeit Y_sparse.T.dot(C_sparse * Y_sparse)

私のタイミング：

In [59]: %timeit ytcy_fast(Y)
100 loops, best of 3: 16.1 ms per loop

In [18]: %timeit Y_sparse.T.dot(C_sparse * Y_sparse)
1 loops, best of 3: 282 ms per loop
于 2012-11-07T16:48:43.837 に答える
2

まず、問題で完全な行列反転を実行する必要があると本当に確信していますか？ほとんどの場合、実際に計算する必要があるのはx = A ^ -1 yだけであり、これははるかに簡単に解決できる問題です。

これが本当にそうであれば、完全な行列の逆行列ではなく、逆行列の近似を計算することを検討します。行列の反転は本当にコストがかかるためです。逆行列の効率的な近似については、たとえばランチョスアルゴリズムを参照してください。近似値は、ボーナスとしてまばらに保存できます。さらに、行列とベクトルの演算のみが必要なため、逆行列に完全な行列を格納する必要もありません。

別の方法として、pyoperatorsを使用して、.todenseメソッドを使用して、効率的な行列ベクトル演算を使用して逆行列を計算することもできます。対角行列用の特別なスパースコンテナがあります。

Lanczosアルゴリズムの実装については、pyoperatorsを見ることができます（免責事項：私はこのソフトウェアの共著者の1人です）。

于 2012-11-07T15:46:38.380 に答える
0

質問があったときにそれが可能だったかどうかはわかりません。しかし、今日では、放送はあなたの友達です。n * n対角行列は、行列積で使用される対角要素の配列である必要があります。

>>> n, f = 5, 3
>>> Y = np.random.randint(0, 10, (n, f))
>>> C = np.random.randint(0, 10, (n,))
>>> Y.shape
(5, 3)
>>> C.shape
(5,)
>>> np.all(Y.T @ np.diag(C) @ Y == Y.T*C @ Y)
True

Y.T*C @ Yこれは非連想的であることに注意してください。

>>> Y.T*(C @ Y)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: operands could not be broadcast together with shapes (3,5) (3,)

しかしY.T @ (C[:, np.newaxis]*Y)、期待される結果が得られます。

>>> np.all(Y.T*C @ Y == Y.T@(C[:, np.newaxis]*Y))
True
于 2021-01-19T04:20:19.013 に答える