矩阵¶
>>> 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 imes 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 imes m\) 的 \(0\) 矩阵。
>>> zeros(2, 3)
⎡0 0 0⎤
⎢ ⎥
⎣0 0 0⎦
同样地，ones 创建一个全一矩阵。
>>> ones(3, 2)
⎡1 1⎤
⎢ ⎥
⎢1 1⎥
⎢ ⎥
⎣1 1⎦
要创建对角矩阵，请使用 diag。 diag 的参数可以是数字或矩阵。数字被解释为 \(1 imes 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 返回一个 特征值: 代数重数 对的字典（类似于 根 的输出）。
>>> 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 返回一个元组列表，格式为 (特征值, 代数重数, [特征向量])。
>>> 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)
可能的问题¶
零测试¶
如果你的矩阵运算失败或返回错误答案，常见的原因可能是零测试问题。如果表达式没有正确进行零测试，可能会在寻找高斯消元法的枢轴，或决定矩阵是否可逆，或任何依赖于先前过程的高级函数中带来问题。
目前，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] 以获得社区的详细帮助。
脚注