2010-03-29 108 views
40

所以,我正在做一些Kmeans分類,使用非常稀疏的numpy數組 - 很多很多的零。我想我會用scipy的'sparse'包來減少存儲開銷,但是我對如何創建數組而不是矩陣有些困惑。Scipy稀疏...數組?

我已經通過本教程了關於如何創建稀疏矩陣: http://www.scipy.org/SciPy_Tutorial#head-c60163f2fd2bab79edd94be43682414f18b90df7

要模仿一個數組,我只創建一個1×N個矩陣,但正如你可能已經猜到,Asp.dot(BSP)沒有按」因爲你不能乘以兩個1xN矩陣,所以很有效。我不得不將每個數組轉換爲Nx1,這是非常蹩腳的,因爲我會爲每個點積計算做這件事。

接下來,我試着創建一個NxN矩陣,其中第1行==第1行(這樣您可以乘以兩個矩陣,只需將左上角作爲點乘積),但事實證明這是真的效率低下。

我很樂意使用scipy的稀疏包作爲numpy的數組()的魔術替代品,但是至今我並不確定該怎麼做。

有什麼建議嗎?

+0

請參見下面的註釋,但我最終只是滾動了我自己的稀疏矢量實現,一個「dok」矩陣 – spitzanator 2010-03-30 18:55:41

+0

原始問題鏈接似乎已經死亡。@spitzanator。 – Mark 2016-07-26 13:14:06

回答

31

使用基於行或列的scipy.sparse格式:csc_matrixcsr_matrix

這些使用高效的C語言實現(包括乘法),並且移位是無操作的(特別是如果您調用transpose(copy=False)),就像numpy數組一樣。

編輯:通過ipython一些計時:

import numpy, scipy.sparse 
n = 100000 
x = (numpy.random.rand(n) * 2).astype(int).astype(float) # 50% sparse vector 
x_csr = scipy.sparse.csr_matrix(x) 
x_dok = scipy.sparse.dok_matrix(x.reshape(x_csr.shape)) 

現在x_csrx_dok 50%疏:

print repr(x_csr) 
<1x100000 sparse matrix of type '<type 'numpy.float64'>' 
     with 49757 stored elements in Compressed Sparse Row format> 

而且時機:

timeit numpy.dot(x, x) 
10000 loops, best of 3: 123 us per loop 

timeit x_dok * x_dok.T 
1 loops, best of 3: 1.73 s per loop 

timeit x_csr.multiply(x_csr).sum() 
1000 loops, best of 3: 1.64 ms per loop 

timeit x_csr * x_csr.T 
100 loops, best of 3: 3.62 ms per loop 

所以它看起來像我說謊。轉置非常便宜,但沒有有效的C實現csr * csc(在最新的scipy 0.9.0中)。新的CSR對象在:-(

每個呼叫作爲一個黑客構建的(儘管SciPy的相對穩定,這些天),您可以在稀疏的數據直接做點積:

timeit numpy.dot(x_csr.data, x_csr.data) 
10000 loops, best of 3: 62.9 us per loop 

注意這最後一種方法再次進行了一次numpy密集乘法運算,其稀疏性爲50%,所以它實際上比dot(x, x)快了2倍。

+5

+1 for plain numpy.dot。對於kmeans,您需要argmax(點(k x N箇中心,每個Nvec x));無論如何,中心都會變得密集,所以不妨保持密集。 (雖然爲新中心平均許多稀疏的xs是非常緩慢的。) – denis 2011-07-23 16:53:51

+0

好吧,如果我們把乘法速度放在一邊,OP可能會使用'scipy.cluster.kmeans' ... – Radim 2011-07-23 20:41:07

+3

合理。我更喜歡(advt)[this code](http://stackoverflow.com/questions/5529625/is-it-possible-to-specify-your-own-distance-function-using-scikits-learn-k-means) ,它可以使用scipy.spatial.distance中的任何20多個度量標準;度量對於高維kmeans比算法更重要。 – denis 2011-07-24 09:58:50

1

您可以創建現有的2D稀疏陣列中的一個子類

from scipy.sparse import dok_matrix 

class sparse1d(dok_matrix): 
    def __init__(self, v): 
     dok_matrix.__init__(self, (v,)) 
    def dot(self, other): 
     return dok_matrix.dot(self, other.transpose())[0,0] 

a=sparse1d((1,2,3)) 
b=sparse1d((4,5,6)) 
print a.dot(b) 
+0

不幸的是,這個問題是你必須在飛行中改變dang的東西,當你進行數百萬次比較時,這並沒有什麼意義。我嘗試緩存點產品,但不幸的是,我們不會經常做同樣的點產品,所以沒有多大幫助。 – spitzanator 2010-03-30 18:53:44

0

我不知道它是真的要好得多或更快,但你可以這樣做是爲了避免使用轉:

Asp.multiply(Bsp).sum() 

這隻需要兩個矩陣的元素和元素的乘積並對產品進行求和。你可以製作你使用的任何矩陣格式的子類,它具有上述語句作爲點積。

但是,它可能只是更容易TRANSPOSE(移調)他們:

Asp*Bsp.T 

似乎並不像這麼多的事,但你也可以做一個子類,並修改MUL()方法。

+0

我也嘗試,對於一個矢量[1,2,3],從而形成矩陣: [1,2,3] [2,0,0] [3,0,0] 以兩個這些和乘以(以任何順序)在結果矩陣的左上角給出所需的點積。不幸的是,這種速度嚴重受到負面影響。 – spitzanator 2010-03-30 18:55:12