「SymPy」符号运算(7) Matrix进阶运算与线性代数

在前几篇文章中，我们学习SymPy的各种用法，以及符号矩阵Matrix的基础用法，这一篇我们继续深入学习矩阵Matrix的相关内容，忘记前面内容的小伙伴可以复习一下—— 传送链接： 「SymPy」符号运算(6) 矩阵Matrix及基础运算 「SymPy」符号运算(5) Vector向量及运算 「SymPy」符号运算(4) 微积分与有限差分 「SymPy」符号运算(3) (非)线性方程(组)求解、数列求和、连乘、求极限 「SymPy」符号运算(2) 各种形式输出、表达式的化简合并与展开 「SymPy」符号运算(1) 简介/符号/变量/函数/表达式/等式/不等式/运算符

sympy.matrices官方文档1：https://docs.sympy.org/latest/tutorials/intro-tutorial/matrices.html?highlight=matrix

在线性代数中，一个矩阵 A \boldsymbol A A的列秩是 A \boldsymbol A A的线性独立的纵列的最大数目，方阵的列秩与行秩总是相等的。通过.rank()获得矩阵的秩

输出： [ 0 1 1 0 0 0 1 2 0 ] \left[\begin{matrix}0 & 1 & 1\\0 & 0 & 0\\1 & 2 & 0\end{matrix}\right] ​001​102​100​ ​ 查看A的秩：

用符号矩阵试试：

矩阵 A \boldsymbol A A的逆矩阵为 A − 1 \boldsymbol A^{-1} A−1，通过A.inv()求逆

输出： [ 1 2 x 1 − 1 x ] \left[\begin{matrix}1 & 2\\x & 1 - \frac{1}{x}\end{matrix}\right] [1x​21−x1​​]

输出： [ 1 − x 2 x 2 − x + 1 2 x 2 x 2 − x + 1 x 2 2 x 2 − x + 1 1 − 2 x + 1 − 1 x ] \left[\begin{matrix}\frac{1 - x}{2 x^{2} - x + 1} & \frac{2 x}{2 x^{2} - x + 1}\\\frac{x^{2}}{2 x^{2} - x + 1} & \frac{1}{- 2 x + 1 - \frac{1}{x}}\end{matrix}\right] [2x2−x+11−x​2x2−x+1x2​​2x2−x+12x​−2x+1−x1​1​​] .inv()函数有一个参数method，矩阵为稠密矩阵时默认为使用高斯消去法求逆、为系数矩阵时默认用LDL分解法求逆。所有方法包括：

求迹：.trace()

直接.T属性就可以：

输出： [ 1 x 2 1 − 1 x ] \left[\begin{matrix}1 & x\\2 & 1 - \frac{1}{x}\end{matrix}\right] [12​x1−x1​​]

.C属性或.conjugate()方法可以按元素取共轭：

输出： [ 1 − i 2 − i 2 + i 3 i − i 3 4 5 + i ] \left[\begin{matrix}1 - i & 2 & - i\\2 + i & 3 i & - i\\3 & 4 & 5 + i\end{matrix}\right] ​1−i2+i3​23i4​−i−i5+i​ ​ 获得Hermite共轭阵

[ 1 − i 2 + i 3 2 3 i 4 − i − i 5 + i ] \left[\begin{matrix}1 - i & 2 + i & 3\\2 & 3 i & 4\\- i & - i & 5 + i\end{matrix}\right] ​1−i2−i​2+i3i−i​345+i​ ​

求矩阵M的伴随矩阵，可以使用M.adjugate()方法：

输出： − 2 x + 1 − 1 x -2 x + 1 - \frac{1}{x} −2x+1−x1​ .det()函数也有一个可选的method参数，可以设为：'bareiss', 'domain-ge', berkowitz, 'lu'。

特征值：.eigenvals()

函数：eigenvals(error_when_incomplete=True, **flags)，error_when_incomplete如果为True，特征值未全部求解时会报错；其他参数还有，simplify = True时结果自动化简，rational = True时做计算前浮点数会被替换为有理数，multiple = True则结果返回列表、否则返回字典。

特征向量：.eigenvects()

函数：eigenvects(error_when_incomplete=True, iszerofunc=, **flags)，其中error_when_incomplete和simplify参数同上，其他还有chop=True时将矩阵含的浮点数转为有理数，但是计算的结果会经过eval函数变为小数，可以通过chop=n其中n为正整数来指定精度。

特征多项式： d e l ( λ I − A ) \mathrm{del}(\lambda I - A) del(λI−A)

函数：charpoly(x='lambda', simplify=) ，x是特征值的标识符，默认是lambda；simplify为对特征多项式使用的简化函数，默认为simplify。

输出： PurePoly ⁡ ( λ 2 − λ − 6 , λ , d o m a i n = Z ) \operatorname{PurePoly}{\left( \lambda^{2} - \lambda - 6, \lambda, domain=\mathbb{Z} \right)} PurePoly(λ2−λ−6,λ,domain=Z) 符号x取其他符号都是等价的：

计算向量或矩阵的范数：.norm(ord=None)

范数有如下可选的形式：

简单例子：

计算向量V的标准化向量：V.normalized() ，行、列向量均可，但是矩阵不能计算标准化向量。

[ 14 14 14 7 3 14 14 ] \left[\begin{matrix}\frac{\sqrt{14}}{14} & \frac{\sqrt{14}}{7} & \frac{3 \sqrt{14}}{14}\end{matrix}\right] [1414 ​​​714 ​​​14314 ​​​]

计算可逆矩阵的条件数：A.condition_number()

常用的矩阵分解包括LU分解、LDL分解、QR分解等。

对于LU分解，函数为LUdecomposition()：

则 L = [ 1 0 3 2 1 ] , U = [ 4 3 0 − 3 2 ] L = \left[\begin{matrix}1 & 0\\\frac{3}{2} & 1\end{matrix}\right] , \quad U = \left[\begin{matrix}4 & 3\\0 & - \frac{3}{2}\end{matrix}\right] L=[123​​01​],U=[40​3−23​​] 对于QR分解，函数为.QRdecomposition()：

则 Q = [ 6 7 − 69 175 − 58 175 3 7 158 175 6 175 − 2 7 6 35 − 33 35 ] , R = [ 14 21 − 14 0 175 − 70 0 0 35 ] Q = \left[\begin{matrix}\frac{6}{7} & - \frac{69}{175} & - \frac{58}{175}\\\frac{3}{7} & \frac{158}{175} & \frac{6}{175}\\- \frac{2}{7} & \frac{6}{35} & - \frac{33}{35}\end{matrix}\right] ,\quad R = \left[\begin{matrix}14 & 21 & -14\\0 & 175 & -70\\0 & 0 & 35\end{matrix}\right] Q= ​76​73​−72​​−17569​175158​356​​−17558​1756​−3533​​ ​,R= ​1400​211750​−14−7035​ ​

对于LDL分解，函数为.LDLdecomposition()：

则 A = [ 25 15 − 5 15 18 0 − 5 0 11 ] , D = [ 25 0 0 0 9 0 0 0 9 ] A = \left[\begin{matrix}25 & 15 & -5\\15 & 18 & 0\\-5 & 0 & 11\end{matrix}\right] , \quad D = \left[\begin{matrix}25 & 0 & 0\\0 & 9 & 0\\0 & 0 & 9\end{matrix}\right] A= ​2515−5​15180​−5011​ ​,D= ​2500​090​009​ ​ 此外还有cholesky分解，函数为.cholesky(hermitian=True)：

[ 5 0 0 3 3 0 − 1 1 3 ] \left[\begin{matrix}5 & 0 & 0\\3 & 3 & 0\\-1 & 1 & 3\end{matrix}\right] ​53−1​031​003​ ​

黑塞矩阵（Hessian Matrix），又译作海森矩阵、海瑟矩阵、海塞矩阵等，是一个多元函数的二阶偏导数构成的方阵，描述了函数的局部曲率。常用于牛顿法解决优化问题，利用黑塞矩阵可判定多元函数的极值问题。——百度百科

函数：hessian(f, varlist, constraints=()) ，计算函数f(varlist)的黑塞矩阵，varlist为变量列表，可以是一个序列数据/符号或者行/列向量。

例如，有一个函数f(x, y)，求其黑塞矩阵：

输出： [ ∂ 2 ∂ x 2 f ( x , y ) ∂ 2 ∂ y ∂ x f ( x , y ) ∂ 2 ∂ y ∂ x f ( x , y ) ∂ 2 ∂ y 2 f ( x , y ) ] \left[\begin{matrix}\frac{\partial^{2}}{\partial x^{2}} f{\left(x,y \right)} & \frac{\partial^{2}}{\partial y\partial x} f{\left(x,y \right)}\\\frac{\partial^{2}}{\partial y\partial x} f{\left(x,y \right)} & \frac{\partial^{2}}{\partial y^{2}} f{\left(x,y \right)}\end{matrix}\right] [∂x2∂2​f(x,y)∂y∂x∂2​f(x,y)​∂y∂x∂2​f(x,y)∂y2∂2​f(x,y)​] 再如：

[ 0 0 ∂ ∂ x g ( x , y ) ∂ ∂ y g ( x , y ) 0 0 ∂ ∂ x h ( x , y ) ∂ ∂ y h ( x , y ) ∂ ∂ x g ( x , y ) ∂ ∂ x h ( x , y ) ∂ 2 ∂ x 2 f ( x , y ) ∂ 2 ∂ y ∂ x f ( x , y ) ∂ ∂ y g ( x , y ) ∂ ∂ y h ( x , y ) ∂ 2 ∂ y ∂ x f ( x , y ) ∂ 2 ∂ y 2 f ( x , y ) ] \left[\begin{matrix}0 & 0 & \frac{\partial}{\partial x} g{\left(x,y \right)} & \frac{\partial}{\partial y} g{\left(x,y \right)}\\0 & 0 & \frac{\partial}{\partial x} h{\left(x,y \right)} & \frac{\partial}{\partial y} h{\left(x,y \right)}\\\frac{\partial}{\partial x} g{\left(x,y \right)} & \frac{\partial}{\partial x} h{\left(x,y \right)} & \frac{\partial^{2}}{\partial x^{2}} f{\left(x,y \right)} & \frac{\partial^{2}}{\partial y\partial x} f{\left(x,y \right)}\\\frac{\partial}{\partial y} g{\left(x,y \right)} & \frac{\partial}{\partial y} h{\left(x,y \right)} & \frac{\partial^{2}}{\partial y\partial x} f{\left(x,y \right)} & \frac{\partial^{2}}{\partial y^{2}} f{\left(x,y \right)}\end{matrix}\right] ​00∂x∂​g(x,y)∂y∂​g(x,y)​00∂x∂​h(x,y)∂y∂​h(x,y)​∂x∂​g(x,y)∂x∂​h(x,y)∂x2∂2​f(x,y)∂y∂x∂2​f(x,y)​∂y∂​g(x,y)∂y∂​h(x,y)∂y∂x∂2​f(x,y)∂y2∂2​f(x,y)​ ​

在向量微积分中，雅可比矩阵是一阶偏导数以一定方式排列成的矩阵，其行列式称为雅可比行列式。雅可比矩阵的重要性在于它体现了一个可微与给出点的最优线性逼近。因此，雅可比矩阵类似于多元函数的导数。——百度百科

求列向量 X \boldsymbol{X} X对列向量 Y \boldsymbol{Y} Y的雅克比矩阵，

[ cos ⁡ ( ϕ ) − ρ sin ⁡ ( ϕ ) sin ⁡ ( ϕ ) ρ cos ⁡ ( ϕ ) 2 ρ 0 ] \left[\begin{matrix}\cos{\left(\phi \right)} & - \rho \sin{\left(\phi \right)}\\\sin{\left(\phi \right)} & \rho \cos{\left(\phi \right)}\\2 \rho & 0\end{matrix}\right] ​cos(ϕ)sin(ϕ)2ρ​−ρsin(ϕ)ρcos(ϕ)0​ ​

函数：jordan_cell()，产生一个Jordan块

[ x 1 0 0 0 x 1 0 0 0 x 1 0 0 0 x ] \left[\begin{matrix}x & 1 & 0 & 0\\0 & x & 1 & 0\\0 & 0 & x & 1\\0 & 0 & 0 & x\end{matrix}\right] ​x000​1x00​01x0​001x​ ​

也就是求离散数据的最佳平方逼近，函数：solve_least_squares(rhs, method='CH')，rhs为右端向量，默认方法为CH (Cholesky分解) 算法，此外还有‘LDL’、‘QR’、PINV等算法。

来个例子，设待拟合的参数为 x x x, y y y，实验测得结果 f f f与 x x x, y y y的关系为： { x + 2 y → 8 2 x + 3 y → 14 3 x + 4 y → 18 \left\{\begin{matrix} x + 2y \rightarrow 8\\ 2x+3y\rightarrow 14\\ 3x+4y\rightarrow 18 \end{matrix}\right. ⎩ ⎨ ⎧​x+2y→82x+3y→143x+4y→18​

[ 5 3 10 3 ] \left[\begin{matrix}\dfrac{5}{3}\\\dfrac{10}{3}\end{matrix}\right] ​35​310​​ ​

线性方程组求解， A x = B \boldsymbol{Ax} = \boldsymbol{B} Ax=B

[ 3 7 5 ] \left[\begin{matrix}3\\7\\5\end{matrix}\right] ​375​ ​

此外还有.LDLsolve()、.QRsolve()、.diagonal_solve() (当A为对角阵时使用)、.cholesky_solve()、.lower_triangular_solve() (当A为下三角矩阵)、.upper_triangular_solve() (当A为上三角矩阵)、.gauss_jordan_solve() (消去法）等，这些求解方法基于不同的算法，根据实际需要进行选用和尝试即可。

等价地，还可以直接使用.solve(rhs, method=‘GJ’)函数，rhs为右端向量，方法包括：'GJ'(Gauss Jordan消去法）、‘LU’（LU分解）、‘QR’（QR分解）、‘PINV’、‘CH’（Cholesky算法）、'LDL’（LDL分解）等。

[ 3 7 5 ] \left[\begin{matrix}3\\7\\5\end{matrix}\right] ​375​ ​

格拉姆-施密特正交化：把线性无关向量系进行正交化。

方法：GramSchmidt(vlist, orthonormal=False) ，其中vlist是列向量组成的列表，orthonormal为是否进行标准化。

则 o u t 1 = [ [ 2 3 5 ] ,   [ 23 19 63 19 − 47 19 ] ,   [ 1692 353 − 1551 706 − 423 706 ] ] , o u t 2 = [ [ 38 19 3 38 38 5 38 38 ] ,   [ 23 6707 6707 63 6707 6707 − 47 6707 6707 ] ,   [ 12 706 353 − 11 706 706 − 3 706 706 ] ] \mathrm{out1} = \left[ \left[\begin{matrix}2\\3\\5\end{matrix}\right], \ \left[\begin{matrix}\frac{23}{19}\\\frac{63}{19}\\- \frac{47}{19}\end{matrix}\right], \ \left[\begin{matrix}\frac{1692}{353}\\- \frac{1551}{706}\\- \frac{423}{706}\end{matrix}\right]\right] ,\quad \mathrm{out2} = \left[ \left[\begin{matrix}\frac{\sqrt{38}}{19}\\\frac{3 \sqrt{38}}{38}\\\frac{5 \sqrt{38}}{38}\end{matrix}\right], \ \left[\begin{matrix}\frac{23 \sqrt{6707}}{6707}\\\frac{63 \sqrt{6707}}{6707}\\- \frac{47 \sqrt{6707}}{6707}\end{matrix}\right], \ \left[\begin{matrix}\frac{12 \sqrt{706}}{353}\\- \frac{11 \sqrt{706}}{706}\\- \frac{3 \sqrt{706}}{706}\end{matrix}\right]\right] out1= ​ ​235​ ​,  ​1923​1963​−1947​​ ​,  ​3531692​−7061551​−706423​​ ​ ​,out2= ​ ​1938 ​​38338 ​​38538 ​​​ ​,  ​6707236707 ​​6707636707 ​​−6707476707 ​​​ ​,  ​35312706 ​​−70611706 ​​−7063706 ​​​ ​ ​ 有

将矩阵内的符号或表达式替换为数值或者表达式：.subs(, )

输出： [ 3 0 0 0 3 0 0 0 3 ] \left[\begin{matrix}3 & 0 & 0\\0 & 3 & 0\\0 & 0 & 3\end{matrix}\right] ​300​030​003​ ​ 再替换回来：

[ x 0 0 0 x 0 0 0 x ] \left[\begin{matrix}x & 0 & 0\\0 & x & 0\\0 & 0 & x\end{matrix}\right] ​x00​0x0​00x​ ​

矩阵的运算结果可以通过.n()方法来限制输出数值的有效位数。创建一个随机 3 × 3 3\times3 3×3矩阵

[ 0.0173639460819103 0.304478747148788 0.0396927801622903 0.64611327050618 0.749989544306484 0.665449822109664 0.881072456803195 0.503832261345867 0.984496661648428 ] \left[\begin{matrix}0.0173639460819103 & 0.304478747148788 & 0.0396927801622903\\0.64611327050618 & 0.749989544306484 & 0.665449822109664\\0.881072456803195 & 0.503832261345867 & 0.984496661648428\end{matrix}\right] ​0.01736394608191030.646113270506180.881072456803195​0.3044787471487880.7499895443064840.503832261345867​0.03969278016229030.6654498221096640.984496661648428​ ​

控制输出两位有效数字：

[ 0.017 0.3 0.04 0.65 0.75 0.67 0.88 0.5 0.98 ] \left[\begin{matrix}0.017 & 0.3 & 0.04\\0.65 & 0.75 & 0.67\\0.88 & 0.5 & 0.98\end{matrix}\right] ​0.0170.650.88​0.30.750.5​0.040.670.98​ ​

SympPy矩阵之间直接用=符号赋值时，两个标识符指向的是同一个矩阵，修改其中一个、另一个随之变化。

可以通过.copy() 方法产生A矩阵的副本，对副本的修改不会影响的原来的矩阵：

通过.reshape()方法改变矩阵的结构，按行展开后再填充到指定的形状中，返回副本：

[ 0 1 2 3 4 5 6 7 8 ] \left[\begin{matrix}0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8\end{matrix}\right] [0​1​2​3​4​5​6​7​8​]

[ 0 1 2 3 4 5 6 7 8 ] \left[\begin{matrix}0 & 1 & 2\\3 & 4 & 5\\6 & 7 & 8\end{matrix}\right] ​036​147​258​ ​

[ 2 2 1 5 3 1 3 4 26 ] \left[\begin{matrix}\sqrt{2} & 2 & 1\\\sqrt{5} & 3 & 1\\3 & 4 & \sqrt{26}\end{matrix}\right] ​2 ​5 ​3​234​1126 ​​ ​

[ − 2 3 8 5 ] \left[\begin{matrix}-2 & 3\\8 & 5\end{matrix}\right] [−28​35​]

[ 7 3 11 7 ] \left[\begin{matrix}7 & 3\\11 & 7\end{matrix}\right] [711​37​]

[ − 3 2 15 4 ] \left[\begin{matrix}-3 & 2\\15 & 4\end{matrix}\right] [−315​24​]

按元素取微分：.diff(var)

[ 1 0 0 0 ] \left[\begin{matrix}1 & 0\\0 & 0\end{matrix}\right] [10​00​]

按元素取积分：.integrate(var)

[ x 2 2 x y x 0 ] \left[\begin{matrix}\frac{x^{2}}{2} & x y\\x & 0\end{matrix}\right] [2x2​x​xy0​]

以上基本就是SymPy中关于Matrices的大部分内容了，具体函数参数也可以通过help在命令行中详询。后续将讲解SymPy模块中的张量Tensor及其对应的数据类型Array，并且打算针对张量进行图文浅析，感谢关注🌹！

Meurer A, Smith CP, Paprocki M, Čertík O, Kirpichev SB, Rocklin M, Kumar A, Ivanov S, Moore JK, Singh S, Rathnayake T, Vig S, Granger BE, Muller RP, Bonazzi F, Gupta H, Vats S, Johansson F, Pedregosa F, Curry MJ, Terrel AR, Roučka Š, Saboo A, Fernando I, Kulal S, Cimrman R, Scopatz A. (2017) SymPy: symbolic computing in Python. PeerJ Computer Science 3:e103 https://doi.org/10.7717/peerj-cs.103 ↩︎