c语言解三元一次方程组_python求解含有稀疏矩阵的线性方程组

python的数据科学计算包numpy几乎可以解决所有的计算问题，使用np.linalg.solve即可解线性方程组。但是最近遇到一个问题，求解Ax=b方程组，如果A是一个很大的稀疏矩阵，numpy初始化矩阵会导致内存不足。例如：

import numpy as np
a = np.zeros([4000000,3000000], dtype=float)

运行代码会报如下错误：

MemoryError: Unable to allocate array with shape (4000000, 3000000) and data type 
float64

c++中Eigen库提供了稀疏矩阵创建方法以及相应的线性方程组求解接口。查询资料后了解到python也有scipy包支持创建稀疏矩阵，并且提供了求解线性方程组的接口函数，下文做一整理记录。

• 核心方法：scipy.sparse.linalg.spsolve
• 稀疏矩阵结构体：scipy.sparse.csc_matrix

首先，使用三元组构建稀疏矩阵：

import numpy as np
from scipy.sparse import csc_matrix


row = np.array([0, 2, 2, 0, 1, 2, 0])
col = np.array([0, 0, 1, 2, 2, 2, 0])
data = np.array([1, 2, 3, 4, 5, 6, 5])
sp_A = csc_matrix((data, (row, col)), shape=(3, 3), dtype=float)

使用toarray函数，将稀疏矩阵转换为np.array格式：

sp_A.toarray()

输出：

array([[6, 0, 4],
       [0, 0, 5],
       [2, 3, 6]])

这里的data只支持一维数据，row、col可以重复，同一位置的数值相加作为矩阵该位置的最后数值。当稀疏矩阵的shape很大的时候，使用toarray函数会出现内存不足。

若需要解方程Ax=b，此时已有稀疏矩阵sp_A， 构造b向量，解方程即可：

from scipy.sparse.linalg import spsolve


b = np.array([4, 5, 6])
x = spsolve(sp_A, b)
x

解得x：

array([0., 0., 1.])

验证方程解正确性：

>>> sp_A*x
array([4., 5., 6.])

与b向量一致。

更多scipy库中的矩阵操作、解方程函数，请参考：

scipy中解方程相关接口

scipy矩阵初始化及矩阵相关操作

---

完。