我怎样才能在这个程序中创建一个稀疏矩阵而不是一个密集矩阵？

Question

我有这个 delta 函数，它有 3 个案例。 mask1、mask2 并且如果 none 满足 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的大部分元素属于delta的第三种情况，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 returns 结果与其输入的形状相同。 E = deltx*delty 表明 deltax 和 deltay 是相同的形状（或至少可广播）。

由于稀疏矩阵有一个 .multiply 方法来进行元素乘法，我们可以专注于生成稀疏 delta 矩阵。

如果您负担得起制造 rx 的内存和几个面具，那么您也可以负担得起制造 deltax（所有尺寸相同）。即使 deltax 有很多零，使其变稠密也可能是最快的。

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

这看起来像你在 delta 中所做的事情的本质，至少有一个面具：

从二维数组开始：

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

我怎样才能在这个程序中创建一个稀疏矩阵而不是一个密集矩阵？

How can I create a sparse matrix instead of a dense one in this program?

python

numpy

sparse-matrix