矩阵¶
>>> from sympy import *
>>> init_printing(use_unicode=True)
要在 SymPy 中创建矩阵，请使用 Matrix 对象。矩阵的构造方法是提供构成矩阵的行向量的列表。例如，要构造矩阵
使用
>>> Matrix([[1, -1], [3, 4], [0, 2]])
⎡1 -1⎤
⎢ ⎥
⎢3 4 ⎥
⎢ ⎥
⎣0 2 ⎦
为了便于创建列向量，元素列表被视为列向量。
>>> Matrix([1, 2, 3])
⎡1⎤
⎢ ⎥
⎢2⎥
⎢ ⎥
⎣3⎦
矩阵的操作与 SymPy 或 Python 中的任何其他对象相同。
>>> M = Matrix([[1, 2, 3], [3, 2, 1]])
>>> N = Matrix([0, 1, 1])
>>> M*N
⎡5⎤
⎢ ⎥
⎣3⎦
关于 SymPy 矩阵需要注意的一点是，与 SymPy 中的所有其他对象不同，它们是可变的。这意味着它们可以在原地修改，我们将在下面看到。这样做的缺点是 Matrix 不能在需要不可变性的地方使用，例如在其他 SymPy 表达式内部或作为字典的键。如果你需要 Matrix 的不可变版本，请使用 ImmutableMatrix。
基本操作¶
以下是一些 Matrix 的基本操作。
形状¶
要获取矩阵的形状，请使用 shape() 函数。
>>> from sympy import shape
>>> M = Matrix([[1, 2, 3], [-2, 0, 4]])
>>> M
⎡1 2 3⎤
⎢ ⎥
⎣-2 0 4⎦
>>> shape(M)
(2, 3)
访问行和列¶
要获取矩阵的单个行或列，请使用 row 或 col。例如，M.row(0) 将获取第一行。 M.col(-1) 将获取最后一列。
>>> M.row(0)
[1 2 3]
>>> M.col(-1)
⎡3⎤
⎢ ⎥
⎣4⎦
删除和插入行和列¶
要删除行或列，请使用 row_del 或 col_del。这些操作将就地修改矩阵。
>>> M.col_del(0)
>>> M
⎡2 3⎤
⎢ ⎥
⎣0 4⎦
>>> M.row_del(1)
>>> M
[2 3]
要插入行或列，请使用 row_insert 或 col_insert。这些操作不就地执行。
>>> M
[2 3]
>>> M = M.row_insert(1, Matrix([[0, 4]]))
>>> M
⎡2 3⎤
⎢ ⎥
⎣0 4⎦
>>> M = M.col_insert(0, Matrix([1, -2]))
>>> M
⎡1 2 3⎤
⎢ ⎥
⎣-2 0 4⎦
除非明确说明，否则下面提到的方法不就地执行。通常，不就地执行的方法将返回一个新的 Matrix，而就地执行的方法将返回 None。
基本方法¶
如上所述，简单的操作，如加法、乘法和幂运算，只需使用 +、 * 和 ** 即可完成。要查找矩阵的逆矩阵，只需将其提升到 -1 次方。
>>> M = Matrix([[1, 3], [-2, 3]])
>>> N = Matrix([[0, 3], [0, 7]])
>>> M + N
⎡1 6 ⎤
⎢ ⎥
⎣-2 10⎦
>>> M*N
⎡0 24⎤
⎢ ⎥
⎣0 15⎦
>>> 3*M
⎡3 9⎤
⎢ ⎥
⎣-6 9⎦
>>> M**2
⎡-5 12⎤
⎢ ⎥
⎣-8 3 ⎦
>>> M**-1
⎡1/3 -1/3⎤
⎢ ⎥
⎣2/9 1/9 ⎦
>>> N**-1
Traceback (most recent call last):
...
NonInvertibleMatrixError: Matrix det == 0; not invertible.
要对矩阵进行转置，请使用 T。
>>> M = Matrix([[1, 2, 3], [4, 5, 6]])
>>> M
⎡1 2 3⎤
⎢ ⎥
⎣4 5 6⎦
>>> M.T
⎡1 4⎤
⎢ ⎥
⎢2 5⎥
⎢ ⎥
⎣3 6⎦
矩阵构造函数¶
存在几个构造函数用于创建常用矩阵。要创建单位矩阵，请使用 eye。 eye(n) 将创建一个 \(n\times n\) 单位矩阵。
>>> eye(3)
⎡1 0 0⎤
⎢ ⎥
⎢0 1 0⎥
⎢ ⎥
⎣0 0 1⎦
>>> eye(4)
⎡1 0 0 0⎤
⎢ ⎥
⎢0 1 0 0⎥
⎢ ⎥
⎢0 0 1 0⎥
⎢ ⎥
⎣0 0 0 1⎦
要创建全零矩阵，请使用 zeros。 zeros(n, m) 创建一个 \(n\times m\) 的 \(0\) 矩阵。
>>> zeros(2, 3)
⎡0 0 0⎤
⎢ ⎥
⎣0 0 0⎦
类似地，ones 创建一个全一的矩阵。
>>> ones(3, 2)
⎡1 1⎤
⎢ ⎥
⎢1 1⎥
⎢ ⎥
⎣1 1⎦
要创建对角矩阵，请使用 diag。 diag 的参数可以是数字或矩阵。数字被解释为 \(1\times 1\) 矩阵。矩阵沿对角线堆叠。其余元素用 \(0\) 填充。
>>> diag(1, 2, 3)
⎡1 0 0⎤
⎢ ⎥
⎢0 2 0⎥
⎢ ⎥
⎣0 0 3⎦
>>> diag(-1, ones(2, 2), Matrix([5, 7, 5]))
⎡-1 0 0 0⎤
⎢ ⎥
⎢0 1 1 0⎥
⎢ ⎥
⎢0 1 1 0⎥
⎢ ⎥
⎢0 0 0 5⎥
⎢ ⎥
⎢0 0 0 7⎥
⎢ ⎥
⎣0 0 0 5⎦
高级方法¶
行列式¶
要计算矩阵的行列式，请使用 det。
>>> M = Matrix([[1, 0, 1], [2, -1, 3], [4, 3, 2]])
>>> M
⎡1 0 1⎤
⎢ ⎥
⎢2 -1 3⎥
⎢ ⎥
⎣4 3 2⎦
>>> M.det()
-1
RREF¶
要将矩阵转换为行阶梯形矩阵，请使用 rref。 rref 返回一个包含两个元素的元组。第一个是行阶梯形矩阵，第二个是主元列的索引元组。
>>> M = Matrix([[1, 0, 1, 3], [2, 3, 4, 7], [-1, -3, -3, -4]])
>>> M
⎡1 0 1 3 ⎤
⎢ ⎥
⎢2 3 4 7 ⎥
⎢ ⎥
⎣-1 -3 -3 -4⎦
>>> M.rref()
⎛⎡1 0 1 3 ⎤ ⎞
⎜⎢ ⎥ ⎟
⎜⎢0 1 2/3 1/3⎥, (0, 1)⎟
⎜⎢ ⎥ ⎟
⎝⎣0 0 0 0 ⎦ ⎠
注意
rref 返回的元组的第一个元素是 Matrix 类型。第二个是 tuple 类型。
零空间¶
要查找矩阵的零空间，请使用 nullspace。 nullspace 返回一个 list，其中包含跨越矩阵零空间的列向量。
>>> M = Matrix([[1, 2, 3, 0, 0], [4, 10, 0, 0, 1]])
>>> M
⎡1 2 3 0 0⎤
⎢ ⎥
⎣4 10 0 0 1⎦
>>> M.nullspace()
⎡⎡-15⎤ ⎡0⎤ ⎡ 1 ⎤⎤
⎢⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎥
⎢⎢ 6 ⎥ ⎢0⎥ ⎢-1/2⎥⎥
⎢⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎥
⎢⎢ 1 ⎥, ⎢0⎥, ⎢ 0 ⎥⎥
⎢⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎥
⎢⎢ 0 ⎥ ⎢1⎥ ⎢ 0 ⎥⎥
⎢⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎥
⎣⎣ 0 ⎦ ⎣0⎦ ⎣ 1 ⎦⎦
列空间¶
要查找矩阵的列空间，请使用 columnspace。 columnspace 返回一个 list，其中包含跨越矩阵列空间的列向量。
>>> M = Matrix([[1, 1, 2], [2 ,1 , 3], [3 , 1, 4]])
>>> M
⎡1 1 2⎤
⎢ ⎥
⎢2 1 3⎥
⎢ ⎥
⎣3 1 4⎦
>>> M.columnspace()
⎡⎡1⎤ ⎡1⎤⎤
⎢⎢ ⎥ ⎢ ⎥⎥
⎢⎢2⎥, ⎢1⎥⎥
⎢⎢ ⎥ ⎢ ⎥⎥
⎣⎣3⎦ ⎣1⎦⎦
特征值、特征向量和对角化¶
要查找矩阵的特征值，请使用 eigenvals。 eigenvals 返回一个字典，其中包含 eigenvalue: algebraic_multiplicity 对（类似于 roots 的输出）。
>>> M = Matrix([[3, -2, 4, -2], [5, 3, -3, -2], [5, -2, 2, -2], [5, -2, -3, 3]])
>>> M
⎡3 -2 4 -2⎤
⎢ ⎥
⎢5 3 -3 -2⎥
⎢ ⎥
⎢5 -2 2 -2⎥
⎢ ⎥
⎣5 -2 -3 3 ⎦
>>> M.eigenvals()
{-2: 1, 3: 1, 5: 2}
这意味着 M 的特征值为 -2、3 和 5，特征值 -2 和 3 的代数重数为 1，特征值 5 的代数重数为 2。
要查找矩阵的特征向量，请使用 eigenvects。 eigenvects 返回一个包含以下形式元组的列表： (eigenvalue, algebraic_multiplicity, [eigenvectors])。
>>> M.eigenvects()
⎡⎛ ⎡⎡0⎤⎤⎞ ⎛ ⎡⎡1⎤⎤⎞ ⎛ ⎡⎡1⎤ ⎡0 ⎤⎤⎞⎤
⎢⎜ ⎢⎢ ⎥⎥⎟ ⎜ ⎢⎢ ⎥⎥⎟ ⎜ ⎢⎢ ⎥ ⎢ ⎥⎥⎟⎥
⎢⎜ ⎢⎢1⎥⎥⎟ ⎜ ⎢⎢1⎥⎥⎟ ⎜ ⎢⎢1⎥ ⎢-1⎥⎥⎟⎥
⎢⎜-2, 1, ⎢⎢ ⎥⎥⎟, ⎜3, 1, ⎢⎢ ⎥⎥⎟, ⎜5, 2, ⎢⎢ ⎥, ⎢ ⎥⎥⎟⎥
⎢⎜ ⎢⎢1⎥⎥⎟ ⎜ ⎢⎢1⎥⎥⎟ ⎜ ⎢⎢1⎥ ⎢0 ⎥⎥⎟⎥
⎢⎜ ⎢⎢ ⎥⎥⎟ ⎜ ⎢⎢ ⎥⎥⎟ ⎜ ⎢⎢ ⎥ ⎢ ⎥⎥⎟⎥
⎣⎝ ⎣⎣1⎦⎦⎠ ⎝ ⎣⎣1⎦⎦⎠ ⎝ ⎣⎣0⎦ ⎣1 ⎦⎦⎠⎦
这表明，例如，特征值 5 也具有 2 的几何重数，因为它有两个特征向量。由于所有特征值的代数重数和几何重数相同，因此 M 是可对角化的。
要对角化矩阵，请使用 diagonalize。 diagonalize 返回一个元组 \((P, D)\)，其中 \(D\) 是对角矩阵，并且 \(M = PDP^{-1}\)。
>>> P, D = M.diagonalize()
>>> P
⎡0 1 1 0 ⎤
⎢ ⎥
⎢1 1 1 -1⎥
⎢ ⎥
⎢1 1 1 0 ⎥
⎢ ⎥
⎣1 1 0 1 ⎦
>>> D
⎡-2 0 0 0⎤
⎢ ⎥
⎢0 3 0 0⎥
⎢ ⎥
⎢0 0 5 0⎥
⎢ ⎥
⎣0 0 0 5⎦
>>> P*D*P**-1
⎡3 -2 4 -2⎤
⎢ ⎥
⎢5 3 -3 -2⎥
⎢ ⎥
⎢5 -2 2 -2⎥
⎢ ⎥
⎣5 -2 -3 3 ⎦
>>> P*D*P**-1 == M
True
请注意，由于 eigenvects 还包括特征值，因此如果您也想要特征向量，则应使用它而不是 eigenvals。但是，由于计算特征向量通常会很昂贵，因此如果只希望查找特征值，则应优先使用 eigenvals。
如果您只想要特征多项式，请使用 charpoly。这比 eigenvals 更有效，因为有时符号根的计算成本很高。
>>> lamda = symbols('lamda')
>>> p = M.charpoly(lamda)
>>> factor(p.as_expr())
2
(λ - 5) ⋅(λ - 3)⋅(λ + 2)
可能出现的问题¶
零测试¶
如果您的矩阵运算失败或返回错误答案，常见原因可能是零测试。如果某个表达式没有被正确地零测试，它可能会在寻找高斯消元的 pivots、判断矩阵是否可逆或依赖于先前过程的任何高级函数时带来问题。
目前，SymPy 的默认零测试方法 _iszero 仅保证在某些有限的数值和符号域中是准确的，而超出其可判定性的任何复杂表达式都被视为 None，其行为类似于逻辑 False。
使用零测试过程的方法列表如下
echelon_form , is_echelon , rank , rref , nullspace , eigenvects , inverse_ADJ , inverse_GE , inverse_LU , LUdecomposition , LUdecomposition_Simple , LUsolve
他们具有 iszerofunc 属性，允许用户指定零测试方法，该方法可以接受任何具有单个输入和布尔输出的函数，同时默认为 _iszero。
以下是如何解决由不足的零测试导致的问题的示例。虽然该特定矩阵的输出已经得到了改进，但下面的技术仍然值得关注。 [1] [2] [3]
>>> from sympy import *
>>> q = Symbol("q", positive = True)
>>> m = Matrix([
... [-2*cosh(q/3), exp(-q), 1],
... [ exp(q), -2*cosh(q/3), 1],
... [ 1, 1, -2*cosh(q/3)]])
>>> m.nullspace()
[]
您可以通过在启用警告的情况下注入自定义零测试来跟踪哪个表达式被低估了。
>>> import warnings
>>>
>>> def my_iszero(x):
... result = x.is_zero
...
... # Warnings if evaluated into None
... if result is None:
... warnings.warn("Zero testing of {} evaluated into None".format(x))
... return result
...
>>> m.nullspace(iszerofunc=my_iszero)
__main__:9: UserWarning: Zero testing of 4*cosh(q/3)**2 - 1 evaluated into None
__main__:9: UserWarning: Zero testing of (-exp(q) - 2*cosh(q/3))*(-2*cosh(q/3) - exp(-q)) - (4*cosh(q/3)**2 - 1)**2 evaluated into None
__main__:9: UserWarning: Zero testing of 2*exp(q)*cosh(q/3) - 16*cosh(q/3)**4 + 12*cosh(q/3)**2 + 2*exp(-q)*cosh(q/3) evaluated into None
__main__:9: UserWarning: Zero testing of -(4*cosh(q/3)**2 - 1)*exp(-q) - 2*cosh(q/3) - exp(-q) evaluated into None
[]
在这种情况下，(-exp(q) - 2*cosh(q/3))*(-2*cosh(q/3) - exp(-q)) - (4*cosh(q/3)**2 - 1)**2 应该为零，但零测试未能捕获。这可能意味着应该引入更强的零测试。对于此特定示例，将表达式重写为指数形式并应用简化可以使双曲线的零测试更强，同时对其他多项式或超越函数无害。
>>> def my_iszero(x):
... result = x.rewrite(exp).simplify().is_zero
...
... # Warnings if evaluated into None
... if result is None:
... warnings.warn("Zero testing of {} evaluated into None".format(x))
... return result
...
>>> m.nullspace(iszerofunc=my_iszero)
__main__:9: UserWarning: Zero testing of -2*cosh(q/3) - exp(-q) evaluated into None
⎡⎡ ⎛ q ⎛q⎞⎞ -q 2⎛q⎞ ⎤⎤
⎢⎢- ⎜- ℯ - 2⋅cosh⎜─⎟⎟⋅ℯ + 4⋅cosh ⎜─⎟ - 1⎥⎥
⎢⎢ ⎝ ⎝3⎠⎠ ⎝3⎠ ⎥⎥
⎢⎢─────────────────────────────────────────⎥⎥
⎢⎢ ⎛ 2⎛q⎞ ⎞ ⎛q⎞ ⎥⎥
⎢⎢ 2⋅⎜4⋅cosh ⎜─⎟ - 1⎟⋅cosh⎜─⎟ ⎥⎥
⎢⎢ ⎝ ⎝3⎠ ⎠ ⎝3⎠ ⎥⎥
⎢⎢ ⎥⎥
⎢⎢ ⎛ q ⎛q⎞⎞ ⎥⎥
⎢⎢ -⎜- ℯ - 2⋅cosh⎜─⎟⎟ ⎥⎥
⎢⎢ ⎝ ⎝3⎠⎠ ⎥⎥
⎢⎢ ──────────────────── ⎥⎥
⎢⎢ 2⎛q⎞ ⎥⎥
⎢⎢ 4⋅cosh ⎜─⎟ - 1 ⎥⎥
⎢⎢ ⎝3⎠ ⎥⎥
⎢⎢ ⎥⎥
⎣⎣ 1 ⎦⎦
您可以清楚地看到 nullspace 在注入替代零测试后返回了正确的结果。
请注意，这种方法仅适用于某些有限的矩阵情况，这些矩阵仅包含数值、双曲函数和指数函数。对于其他矩阵，您应该使用针对其域选择的不同方法。
可能的建议是，要么利用重写和简化，以速度为代价 [4]，要么使用随机数值测试，以准确性为代价 [5] 。
如果您想知道为什么没有通用的零测试算法可以适用于任何符号实体，那是因为零测试的判定性问题一直存在 [6] ，不仅是 SymPy，其他计算机代数系统 [7] [8] 也会面临同样的根本问题。
然而，发现任何零测试的缺陷都可以提供一些很好的例子来改进 SymPy，所以如果您遇到过这种情况，可以向 SymPy 问题跟踪器报告问题 [9] ，以获得社区的详细帮助。
脚注