numpy矩陣乘法到三角形/稀疏存儲？

Question

我正在處理一個非常大的稀疏矩陣乘法（matmul）問題。舉個例子，讓我們假設：

• A是一個二元（75 x 200,000）矩陣。它很稀疏，所以我使用csc進行存儲。我需要做以下MATMUL操作：

• B = A.transpose（）* A

• 輸出將是尺寸200Kx200K的稀疏和對稱矩陣。

不幸的是，B將是方式到大的RAM來存儲（或「核心」）我的筆記本電腦。另一方面，我很幸運，因爲B有一些屬性可以解決這個問題。由於B將沿着對角線和稀疏對稱對稱，因此我可以使用三角矩陣（上/下）來存儲matmul操作的結果，並且稀疏矩陣存儲格式可以進一步減小大小。

我的問題是......可以提前告訴numpy或scipy，輸出存儲要求會是什麼樣子，以便我可以使用numpy選擇存儲解決方案，並避免「矩陣太大」計算幾分鐘（小時）後運行時錯誤？

換句話說，通過使用近似計數算法分析兩個輸入矩陣的內容，矩陣乘法的存儲需求是否可以近似？

• https://en.wikipedia.org/wiki/Approximate_counting_algorithm

如果不是這樣，我尋找到一個蠻力解決方案。一些涉及到的map/reduce，外的核心存儲，或從下面的網頁鏈接MATMUL細分的解決方案（施特拉森算法）：

一對夫婦的Map/Reduce問題細分解決方案

• http://www.norstad.org/matrix-multiply/index.html
• http://bpgergo.blogspot.com/2011/08/matrix-multiplication-in-python.html

甲外的芯（PyTables）存儲溶液

• Very large matrices using Python and NumPy

一個MATMUL細分的解決方案：

• https://en.wikipedia.org/wiki/Strassen_algorithm
• http://facultyfp.salisbury.edu/taanastasio/COSC490/Fall03/Lectures/FoxMM/example.pdf
• http://eli.thegreenplace.net/2012/01/16/python-parallelizing-cpu-bound-tasks-with-multiprocessing/

預先感謝任何建議，意見或指導！

Answer 1

由於您使用的是矩陣與其轉置的乘積，因此[m, n]的值基本上將會是原始矩陣中的列m和n的點積。

我將使用下面的矩陣作爲一個玩具示例

a = np.array([[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1], 
       [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0], 
       [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1]]) 
>>> np.dot(a.T, a) 
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0], 
     [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0], 
     [0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 2]]) 

它是形狀(3, 12)，並具有7個非零項。其轉置產品當然是(12, 12)的形狀，並有16個非零項，其中6個在對角線，因此它只需要存儲11個元素。

你能得到你的輸出矩陣的大小將是在以下兩種方式之一的一個好主意：

CSR FORMAT

如果你原來的矩陣具有C非零列，您的新矩陣將至多有C**2非零條目，其中C在對角線中，並且確保不爲零，其餘條目中您只需保留一半，這樣至多是(C**2 + C)/2非零值，零元素。當然，其中很多也是零，所以這可能是一個高估。

如果矩陣被存儲在csr格式，則相應的scipy對象的indices屬性具有與所有非零元素的列索引的陣列，所以可以很容易地計算上述估計爲：

>>> a_csr = scipy.sparse.csr_matrix(a) 
>>> a_csr.indices 
array([ 2, 11, 1, 7, 10, 4, 11]) 
>>> np.unique(a_csr.indices).shape[0] 
6 

因此，有6列非零項，因此估計將是至多36個非零項，比真正的16

CSC數據的方式更

如果不是非零元素的列索引我們有行索引，我們實際上可以做一個更好的估計。對於兩列的點積是非零的，它們必須在同一行中有一個非零元素。如果給定行中有非零元素，則它們將爲產品貢獻R**2非零元素。當你對所有行進行總結時，你必然會多次計算一些元素，所以這也是一個上限。

你的矩陣的非零元素的行索引是稀疏CSC矩陣indices屬性，所以這個估計可以計算如下：

>>> a_csc = scipy.sparse.csc_matrix(a) 
>>> a_csc.indices 
array([1, 0, 2, 1, 1, 0, 2]) 
>>> rows, where = np.unique(a_csc.indices, return_inverse=True) 
>>> where = np.bincount(where) 
>>> rows 
array([0, 1, 2]) 
>>> where 
array([2, 3, 2]) 
>>> np.sum(where**2) 
17 

這是混賬接近真實16！它其實並不是一個巧合，這個估計實際上是一樣的：

>>> np.sum(np.dot(a.T,a),axis=None) 
17 

在任何情況下，下面的代碼應該讓你看到，估計是相當不錯的：

def estimate(a) : 
    a_csc = scipy.sparse.csc_matrix(a) 
    _, where = np.unique(a_csc.indices, return_inverse=True) 
    where = np.bincount(where) 
    return np.sum(where**2) 

def test(shape=(10,1000), count=100) : 
    a = np.zeros(np.prod(shape), dtype=int) 
    a[np.random.randint(np.prod(shape), size=count)] = 1 
    print 'a non-zero = {0}'.format(np.sum(a)) 
    a = a.reshape(shape) 
    print 'a.T * a non-zero = {0}'.format(np.flatnonzero(np.dot(a.T, 
                   a)).shape[0]) 
    print 'csc estimate = {0}'.format(estimate(a)) 

>>> test(count=100) 
a non-zero = 100 
a.T * a non-zero = 1065 
csc estimate = 1072 
>>> test(count=200) 
a non-zero = 199 
a.T * a non-zero = 4056 
csc estimate = 4079 
>>> test(count=50) 
a non-zero = 50 
a.T * a non-zero = 293 
csc estimate = 294