python - 使用 `scipy.sparse` 构建一个简单的线性系统

Question

想象一个简单的 3x3 矩阵，我在其上施加以下“边界条件”。以下是一个数组，显示扁平化后的相应索引以及边界值。

I =     B = 
0 3 6   3 0 1
1 4 7   3 0 1
2 5 8   3 0 1

很容易看出，展平后的边界值数组将是

b = 
3 3 3 0 0 0 1 1 1

我还有一个连接结构 C，以一组对的形式给出，

C = 
0 0 1 1 2 3 3 4 4 5 6 7
1 3 2 4 5 4 6 5 7 8 6 8

我想构建一个代表该线性系统的矩阵 A，以使用 spsolve(A, b) 并获取

A =
 1  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0
 0  0  1  0  0  0  0  0  0
 1  0  0 -3  1  0  1  0  0
 0  1  0  1 -4  1  0  1  0
 0  0  1  0  1 -3  0  0  1
 0  0  0  0  0  0  1  0  0
 0  0  0  0  0  0  0  1  0
 0  0  0  0  0  0  0  0  1

x = 
3 3 3 2 2 2 1 1 1

使用 numpy 并操作密集矩阵、将列归零等操作非常简单。然而，随着矩阵变大，我开始耗尽内存，并且求解器非常慢。

我想我应该用这个逻辑构建稀疏矩阵:

1. 使用连接数组和np.ones_like向量之一初始化矩阵
2. 添加矩阵及其转置以填充 LD 区域
3. 其中 b 非零:删除行并在对角线上放置 1(这表示初始条件)
4. 其中 b 为零:沿 axis=1 对矩阵求和，并将总和的负值放入每个对角线上(此平衡汇和源)

这在numpy中完美运行，但我发现在初始化稀疏矩阵后绝对不可能对它们做任何有用的事情。它们不处理项目删除、转置等。我可以了解这些操作吗？

Answer 1

我想到目前为止我都会这样做

看起来有点晦涩，好像我做错了什么

b = np.array([3, 3, 3, 0, 0, 0, 1, 1, 1])
i = (0, 0, 1, 1, 2, 3, 3, 4, 4, 5, 6, 7)
j = (1, 3, 2, 4, 5, 4, 6, 5, 7, 8, 7, 8)

f = b!=0
A = sparse.lil_matrix((b.size, b.size))
A[f,f] = 1
C = sparse.coo_matrix((np.ones_like(i+j), (i+j,j+i)), shape=(b.size, b.size)).tolil()
D = sparse.diags(np.asarray(C.sum(axis=1).T)[0], 0).tocsr()
A[~f,:] = C[~f, :] - D[~f]

print A.toarray()