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

Question

我有这个三角函数，有三个病例。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的大多数元素属于增量的第三种情况，即0。大部分意味着大约99%。因此，我希望有一个稀疏矩阵，而不是密集矩阵，而不是存储0元素，以节省内存。

我该怎么做呢？

Answer 1

创建稀疏矩阵的通常方法是构造三个具有非零值的一维数组，以及它们的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是一维数组。rx和ry分别为2d。delta返回与其输入形状相同的结果。E = deltx*delty认为deltax和deltay是相同的形状(或者至少是可广播的)。

由于稀疏矩阵有一种.multiply方法来进行逐元乘法，所以我们可以把重点放在产生稀疏delta矩阵上。

如果您负担得起制造rx的内存和几个掩码，那么您也可以制作deltax (大小相同)。即使通过deltax有大量的零，它可能是最快使其稠密。

但是，让我们尝试将delta计算作为一个稀疏构建来处理。

这看起来就像您在delta中所做的事情的本质一样，至少有一个掩码：

从2d数组开始：

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