2016-11-17 209 views
2

我有numpy數組Z形狀(k,N)和第二個陣列X形狀(N,n)使用廣播與稀疏scipy矩陣

使用numpy廣播,我可以很容易地獲得新的陣列H與形狀(n,k,N)其切片陣列Z其行已經乘的X列:

H = Z.reshape((1, k, N)) * X.T.reshape((n, 1, N)) 

這工作得很好,是出奇的快。 現在,X是非常稀疏的,我想用稀疏矩陣操作進一步加速這個操作。

但是如果我執行以下操作:

import scipy.sparse as sprs 
spX = sprs.csr_matrix(X) 
H = (Z.reshape((1,k,N))*spX.T.reshape((n,1,N))).dot(Z.T) 

我得到以下錯誤:

Traceback (most recent call last): 
    File "<input>", line 1, in <module> 
    File "C:\Python27\lib\site-packages\scipy\sparse\base.py", line 126, in reshape 
    self.__class__.__name__) 
NotImplementedError: Reshaping not implemented for csc_matrix. 

是否有使用稀疏矩陣scipy廣播的方式?

+0

除了2d限制之外,稀疏數學只適用於其他稀疏矩陣。在處理密集數組時,稀疏矩陣轉換爲密集矩陣。用'spX'來玩耍,看看各種數學運算會發生什麼。 – hpaulj

回答

1

Scipy稀疏矩陣限於2D形狀。但是你可以在「疏」的方式來使用numpy的:

H = np.zeros((n,k,N), np.result_type(Z, X)) 
I, J = np.nonzero(X) 
Z_ = np.broadcast_to(Z, H.shape) 
H[J,:,I] = Z_[J,:,I] * X[I,J,None] 

不幸的是,結果H仍然是一個密集排列。

N.b.使用None進行索引是在所需軸上添加單位長度尺寸的簡便方法。將高級索引與切片結合時的結果順序爲explained in the docs