如何在這個程序中創建一個稀疏矩陣而不是密集矩陣？

Question

我有這個三角功能有3種情況。 mask1，mask2如果他們沒有滿足delta = 0，因爲res = np.zeros

def delta(r, dr): 
    res = np.zeros(r.shape) 
    mask1 = (r >= 0.5*dr) & (r <= 1.5*dr) 
    res[mask1] = (5-3*np.abs(r[mask1])/dr \ 
     - np.sqrt(-3*(1-np.abs(r[mask1])/dr)**2+1)) \ 
     /(6*dr) 
    mask2 = np.logical_not(mask1) & (r <= 0.5*dr) 
    res[mask2] = (1+np.sqrt(-3*(r[mask2]/dr)**2+1))/(3*dr) 
    return res 

然後，我有這等功能，我打電話前，我建立一個數組，E

def matrix_E(nk,X,Y,xhi,eta,dx,dy): 
    rx = abs(X[np.newaxis,:] - xhi[:,np.newaxis]) 
    ry = abs(Y[np.newaxis,:] - eta[:,np.newaxis]) 
    deltx = delta(rx,dx) 
    delty = delta(ry,dy) 
    E = deltx*delty 
    return E 

的事情是， E的大部分元素屬於三角洲的第三種情況，大多數意味着大約99％。 所以，我想有一個稀疏矩陣，而不是一個密集的矩陣，而不是爲了節省內存而存儲0個元素。

任何想法，我怎麼能做到這一點？

Answer 1

創建稀疏矩陣的正常方法是構造三個具有非零值的1d數組，以及它們的i和j索引。然後將它們傳遞給coo_matrix函數。

座標不必按順序排列，因此您可以爲2個非零蒙版案例構建陣列並將它們連接起來。

下面是使用2個掩模

In [107]: x=np.arange(5) 

In [108]: i,j,data=[],[],[] 

In [110]: mask1=x%2==0 
In [111]: mask2=x%2!=0 

In [112]: i.append(x[mask1]) 
In [113]: j.append((x*2)[mask1]) 

In [114]: i.append(x[mask2]) 
In [115]: j.append(x[mask2]) 

In [116]: i=np.concatenate(i) 
In [117]: j=np.concatenate(j) 

In [118]: i 
Out[118]: array([0, 2, 4, 1, 3]) 

In [119]: j 
Out[119]: array([0, 4, 8, 1, 3]) 

In [120]: M=sparse.coo_matrix((x,(i,j))) 

In [121]: print(M) 
    (0, 0) 0 
    (2, 4) 1 
    (4, 8) 2 
    (1, 1) 3 
    (3, 3) 4 

In [122]: M.A 
Out[122]: 
array([[0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 3, 0, 0, 0, 0, 0, 0, 0], 
     [0, 0, 0, 0, 1, 0, 0, 0, 0], 
     [0, 0, 0, 4, 0, 0, 0, 0, 0], 
     [0, 0, 0, 0, 0, 0, 0, 0, 2]]) 

甲coo格式存儲這些3個陣列的樣品構造是，但他們得到排序，並且當轉換爲其它格式和印刷清理。

我可以將其適用於您的案例，但這可能足以讓您開始。

---

看起來像X,Y,xhi,eta是1d陣列。然後是rx和ry。 delta返回結果與其輸入相同的形狀。 E = deltx*delty建議deltax和deltay是相同的形狀（或至少可廣播）。

由於稀疏矩陣有一個.multiply方法來做元素乘法，我們可以專注於生成稀疏的delta矩陣。

如果你負擔的內存使rx，和幾個口罩，那麼你也可以承擔使deltax（所有相同的大小）。即使通過deltax有很多零，它可能是最快的使它密集。

但讓我們嘗試將delta計算爲稀疏構建。

這看起來像你正在做delta什麼，至少有一個屏蔽的essense：

開始用二維數組：

In [138]: r = np.arange(24).reshape(4,6)  
In [139]: mask1 = (r>=8) & (r<=16) 
In [140]: res1 = r[mask1]*0.2 
In [141]: I,J = np.where(mask1) 

得到的載體是：

In [142]: I 
Out[142]: array([1, 1, 1, 1, 2, 2, 2, 2, 2], dtype=int32)  
In [143]: J 
Out[143]: array([2, 3, 4, 5, 0, 1, 2, 3, 4], dtype=int32) 
In [144]: res1 
Out[144]: array([ 1.6, 1.8, 2. , 2.2, 2.4, 2.6, 2.8, 3. , 3.2]) 

製作稀疏矩陣：

In [145]: M=sparse.coo_matrix((res1,(I,J)), r.shape)  
In [146]: M.A 
Out[146]: 
array([[ 0. , 0. , 0. , 0. , 0. , 0. ], 
     [ 0. , 0. , 1.6, 1.8, 2. , 2.2], 
     [ 2.4, 2.6, 2.8, 3. , 3.2, 0. ], 
     [ 0. , 0. , 0. , 0. , 0. , 0. ]]) 

我可以製作另一個稀疏矩陣mask2，並添加這兩個。

In [147]: mask2 = (r>=17) & (r<=22)  
In [148]: res2 = r[mask2]*-0.4 
In [149]: I,J = np.where(mask2) 
In [150]: M2=sparse.coo_matrix((res2,(I,J)), r.shape) 
In [151]: M2.A 
Out[151]: 
array([[ 0. , 0. , 0. , 0. , 0. , 0. ], 
     [ 0. , 0. , 0. , 0. , 0. , 0. ], 
     [ 0. , 0. , 0. , 0. , 0. , -6.8], 
     [-7.2, -7.6, -8. , -8.4, -8.8, 0. ]]) 

... 
In [153]: (M1+M2).A 
Out[153]: 
array([[ 0. , 0. , 0. , 0. , 0. , 0. ], 
     [ 0. , 0. , 1.6, 1.8, 2. , 2.2], 
     [ 2.4, 2.6, 2.8, 3. , 3.2, -6.8], 
     [-7.2, -7.6, -8. , -8.4, -8.8, 0. ]]) 

或者我可以串聯的res1和res2等，使一個稀疏矩陣：

In [156]: I1,J1 = np.where(mask1) 
In [157]: I2,J2 = np.where(mask2) 
In [158]: res12=np.concatenate((res1,res2)) 
In [159]: I12=np.concatenate((I1,I2)) 
In [160]: J12=np.concatenate((J1,J2)) 
In [161]: M12=sparse.coo_matrix((res12,(I12,J12)), r.shape) 

In [162]: M12.A 
Out[162]: 
array([[ 0. , 0. , 0. , 0. , 0. , 0. ], 
     [ 0. , 0. , 1.6, 1.8, 2. , 2.2], 
     [ 2.4, 2.6, 2.8, 3. , 3.2, -6.8], 
     [-7.2, -7.6, -8. , -8.4, -8.8, 0. ]]) 

在這裏，我選擇了面具所以非零值不重疊，但如果這兩種方法工作，他們沒有。這是coo格式的重新設計功能，用於重複索引的值相加。在爲有限元素問題創建稀疏矩陣時非常方便。

---

我也可以從面具創建一個稀疏矩陣得到數組索引：

In [179]: rmask1=sparse.coo_matrix(mask1) 

In [180]: rmask1.row 
Out[180]: array([1, 1, 1, 1, 2, 2, 2, 2, 2], dtype=int32) 

In [181]: rmask1.col 
Out[181]: array([2, 3, 4, 5, 0, 1, 2, 3, 4], dtype=int32) 
In [184]: sparse.coo_matrix((res1, (rmask1.row, rmask1.col)),rmask1.shape).A 
Out[184]: 
array([[ 0. , 0. , 0. , 0. , 0. , 0. ], 
     [ 0. , 0. , 1.6, 1.8, 2. , 2.2], 
     [ 2.4, 2.6, 2.8, 3. , 3.2, 0. ], 
     [ 0. , 0. , 0. , 0. , 0. , 0. ]]) 

我不能，不過，從稀疏版本的r創建模板。 (r>=8) & (r<=16)。稀疏矩陣尚未實現這種不等式測試。但這可能並不重要，因爲r可能並不稀疏。