如何在这个程序中创建一个稀疏矩阵而不是密集矩阵？

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可能并不稀疏。