SymPy 矩阵表达式Matrix Expressions完全指南从 MatrixSymbol 符号运算到分块矩阵【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy导读本文围绕 SymPy 矩阵表达式Matrix Expressions模块展开讲解如何用MatrixSymbol等符号化对象以整个矩阵为单位进行抽象运算而非展开成标量元素内容包括表达式构造、常见表达式类型、矩阵微积分、索引与显式化以及基于BlockMatrix的分块矩阵构造、Schur 补、分块 LDU/LU 分解与block_collapse化简。读完本文你将掌握矩阵表达式的核心对象与运算符、求导与化简机制以及分块矩阵的构建与折叠技巧并能直接复现仓库中的实际示例。本文以 doc/src/modules/matrices/expressions.rst 为骨架辅以 sympy/matrices/expressions 源码与 tests 测试佐证。一、什么是矩阵表达式以矩阵为单位的符号代数普通 SymPy 符号Symbol代表标量而MatrixSymbol代表一个具有确定形状的抽象矩阵。矩阵表达式模块sympy.matrices.expressions允许把X.T*X、(X.T*X).I*Y这样的表达式保持为结构化的表达式对象直到需要时才求值或展开。这是 SymPy 矩阵计算中符号化与显式化的分水岭显式矩阵Matrix展开存储每个元素矩阵表达式则按运算结构存储。1.1 第一个示例转置、乘法与求逆的组合文档开篇给出了最典型的用法 from sympy import MatrixSymbol, Matrix X MatrixSymbol(X, 3, 3) Y MatrixSymbol(Y, 3, 3) (X.T*X).I*Y X**(-1)*X.T**(-1)*Y这里X、Y都是MatrixSymbol而非标量符号。SymPy 自动把(X.T*X).I改写为X**(-1)*X.T**(-1)——这正是乘积转置与乘积求逆的恒等式(A·B)ᵀBᵀ·Aᵀ、(A·B)⁻¹B⁻¹·A⁻¹的符号化体现其实现位于 inverse.py 的Inverse与 matmul.py 的MatMul中MatMul._eval_inverse()在所有因子都是方阵时逐因子反转并翻转顺序。1.2 索引与显式化表达式到稠密矩阵矩阵表达式支持两种落地方式元素索引(X*Y)[1, 2]返回求和式X[1, 0]*Y[0, 2] X[1, 1]*Y[1, 2] X[1, 2]*Y[2, 2]这是由 matmul.py 中MatMul._entry()通过引入哑指标并求和得到的。整体展开Matrix(X)返回 3×3 的稠密矩阵等价于MatrixSymbol.as_explicit()见 matexpr.py后者构造ImmutableDenseMatrixas_mutable()则返回可变的Matrix。MatrixExpr是所有矩阵表达式的公共基类定义于 matexpr.py。它在Expr之上叠加了矩阵语义shape、rows、cols属性is_square判定以及通过_add_handler/_mul_handler把、*、**等运算符重定向到MatAdd/MatMul/MatPow。值得注意的是_iterable False——矩阵表达式默认不被当作可迭代对象从而避免被标量表达式误处理。二、矩阵表达式核心对象速查本节列出文档 Matrix Expressions Core Reference 中收录的核心类并给出源码中的位置与关键行为。全部对象均可从sympy顶层导入参见 sympy/matrices/expressions/init.py。类 / 函数含义源码位置关键行为MatrixExpr矩阵表达式基类matexpr.py定义shape、T、I、.as_explicit()、.as_mutable()、.applyfunc()MatrixSymbol(name, n, m)符号矩阵matexpr.py形状可为符号非负整数支持索引为MatrixElementMatAdd矩阵加法matadd.py支持标量参与加法用于算法内部的非交换处理MatMul矩阵乘法matmul.py自动规范形合并同底幂、约掉X*X⁻¹、合并显式矩阵等MatPow矩阵幂matpow.pyX**-1即逆hadamard_product/HadamardProduct逐元素Hadamard乘积hadamard.py可结合、可交换、全一矩阵为单位元、零矩阵为零元hadamard_power/HadamardPower逐元素幂hadamard.py支持矩阵^标量、标量^矩阵、矩阵^矩阵Inverse矩阵逆inverse.pyA**(-1).doit()才真正求值Transpose转置transpose.pyA.Ttranspose(A*B) → B.T*A.TTrace/trace迹trace.py可重写为Sum支持循环置换规范化FunctionMatrix(rows, cols, lamda)用函数定义矩阵funcmatrix.py支持Lambda、SymPy 函数或 Pythonlambda字符串惰性求值PermutationMatrix/MatrixPermute置换矩阵 / 矩阵置换permutation.py逆即转置可重写为BlockDiagMatrixIdentity(n)单位矩阵special.py乘法单位元I*A → AZeroMatrix(m, n)零矩阵special.py加法单位元Z*A → 0不可逆OneMatrix(m, n)全一矩阵special.pyO.shape[0]**(exp-1) * O的幂规则MatrixUnit(rows, cols, i, j)单元素矩阵special.py只有一个 1 的矩阵转置即换位CompanionMatrix(poly)多项式友矩阵companion.py要求首一、单变量、次数 ≥ 1 的PolyMatrixSet(n, m, set)矩阵集合sets.pyMatrix([[1,2],[3,4]]) in MatrixSet(2,2,S.Reals)2.1 特殊矩阵的实际效果Identity与ZeroMatrix是代数运算中的单位元/零元它们的求值规则直接写在各类的_eval_*方法中。例如Identity._eval_inverse()返回自身ZeroMatrix._eval_inverse()抛出NonInvertibleMatrixErrorMatrix det 0; not invertible.见 special.py。 from sympy import Identity, ZeroMatrix I Identity(3) I*A A Z ZeroMatrix(3, 3) Z.A # 任何零矩阵与矩阵相乘都吸收为 0 0MatrixUnit是单元素矩阵MatrixUnit(3, 4, 1, 2)表示在 (1, 2) 处为 1 的 3×4 矩阵其元素索引借助 KroneckerDelta 表达special.py。它在矩阵求导中充当雅可比矩阵的单元角色——见下文第四节。2.2 FunctionMatrix用函数惰性描述超大矩阵FunctionMatrix(rows, cols, lamda)以Lambda或可接受两个参数的函数按坐标生成元素适合描述元素有规律可言的超大矩阵且不会真的展开存储 from sympy import FunctionMatrix, Lambda, symbols, MatPow i, j, n, m symbols(i,j,n,m) FunctionMatrix(n, m, Lambda((i, j), i j)) FunctionMatrix(n, m, Lambda((i, j), i j)) Y FunctionMatrix(1000, 1000, Lambda((i, j), i j)) isinstance(Y*Y, MatPow) # 仍是表达式对象 True (Y**2)[10, 10] # 惰性求值只算需要的元素 342923500从源码看funcmatrix.py构造时会校验维度为非负整数、函数能接受 2 个参数并把普通函数包装为Lambda_entry(i, j)直接调用self.lamda(i, j)。文档注释明确说明这是以最稀疏的方式表示元素成序列规律的极稠密矩阵的替代方案。2.3 MatrixSet矩阵的集合与成员判定MatrixSet(n, m, set)表示形状为 (n, m)、元素属于某集合的全体矩阵是 SymPy 集合体系Set与矩阵表达式的交汇点 from sympy.matrices import MatrixSet from sympy import S, I, Matrix M MatrixSet(2, 2, setS.Reals) Matrix([[1, 2], [3, 4]]) in M True Matrix([[1, 2], [I, 4]]) in M False其成员判定实现于 sets.py先校验形状匹配形状含符号时返回None表示未定再对每个元素执行set.contains(x)的模糊与运算。三、Hadamard逐元素乘积与逐元素幂hadamard_product(A, B, ...)返回逐元素乘积对象HadamardProduct其第 (i, j) 元素为A[i,j]*B[i,j]hadamard.py from sympy import hadamard_product, MatrixSymbol A MatrixSymbol(A, 2, 3) B MatrixSymbol(B, 2, 3) hadamard_product(A, B)[0, 1] A[0, 1]*B[0, 1]canonicalize()hadamard.py把 Hadamard 积规约到规范形依次应用以下代数规则结合律A.*(B.*C)展平为A.*B.*C单位元全一矩阵OneMatrix是单位元A.*1 → A吸收元任何零矩阵使整个积为零A.*0 → 0重复因子合并A.*A.*A改写为HadamardPower(A, 3)交换律因子按默认排序键排序B.*A → A.*B。HadamardPower(base, exp)覆盖四种情形见 hadamard.py 的文档公式矩阵^标量、标量^矩阵、矩阵^矩阵、标量^标量其导数实现为d(A^∘b)/dx (db/dx · log(A) b · dlog(A)/dx) ⊙ A^∘b对应源码 hadamard.py。测试见 test_hadamard.py例如test_hadamard_product_with_explicit_mat验证显式矩阵与符号矩阵混合时的逐元素合并。四、矩阵表达式的求导矩阵表达式支持对矩阵或矩阵元素求导。由于矩阵对矩阵求导一般是四维数组SymPy 会尽力把平凡或对角维度的结果压回矩阵表达式压不回去时返回数组表达式。文档示例一对矩阵求导得到矩阵表达式 a MatrixSymbol(a, 3, 1) b MatrixSymbol(b, 3, 1) (a.T*X**2*b).diff(X) a*b.T*X.T X.T*a*b.T文档示例二全维情形返回数组表达式 X.diff(X) PermuteDims(ArrayTensorProduct(I, I), (3)(1 2))第二个输出是四维的SymPy 用sympy.tensor.array的ArrayTensorProduct与PermuteDims表示PermuteDims的(3)(1 2)是置换的循环记号。其底层走的是矩阵表达式 → 数组表达式 → 求导 → 转回矩阵表达式的管线核心函数是 matexpr.py 的_matrix_derivative先convert_matrix_to_array转为数组表达式再array_derive求导最后convert_array_to_matrix尽量压回矩阵表达式若expr或x是显式MatrixBase则退化为逐元素求导。求导规则的更多细节MatrixSymbol._eval_derivative对不包含自身的变量求导返回零矩阵对自身求导返回MatrixUnitmatexpr.pyMatrixElement._eval_derivatived(A[i,j])/dA是 (i,j) 处的单元矩阵d(A[i,j])/d(A[k,l]) KroneckerDelta(i,k)·KroneckerDelta(j,l)对逆矩阵元素求导会展开为带Sum的形式matexpr.pyHadamardProduct._eval_derivative按乘积求导法则对每个因子分别求导再求和hadamard.pyTrace._eval_derivatived trace(X)/dX等标量函数对矩阵的导数通过rewrite(Sum)实现trace.py。相关测试集中在 test_derivatives.py如test_derivatives_of_hadamard_expressions、test_derivatives_matrix_norms、test_derivatives_of_scalar_functions_of_traces等可作为行为规范的参考。五、从表达式到稠密矩阵as_explicit 与 Matrix矩阵表达式可通过以下方式求值落地方法返回类型说明expr.as_explicit()ImmutableDenseMatrix逐元素展开为不可变稠密矩阵形状为符号时抛ValueErrormatexpr.pyexpr.as_mutable()Matrix展开为可变稠密矩阵Matrix(expr)Matrix顶层Matrix构造器接受表达式expr[i, j]/expr[i]元素单索引按行优先展开要求列数已知见 matexpr.pyexpr[i:j, k:l]MatrixSlice子矩阵切片视图slice.py X MatrixSymbol(X, 3, 3) Matrix(X) Matrix([ [X[0, 0], X[0, 1], X[0, 2]], [X[1, 0], X[1, 1], X[1, 2]], [X[2, 0], X[2, 1], X[2, 2]]])注意形状为符号时如MatrixSymbol(A, n, n)不能显式化as_explicit()会抛错这类表达式应保持符号形态参与后续推导。六、分块矩阵BlockMatrix 与 BlockDiagMatrix分块矩阵允许用较小的子块构造大矩阵子块既可以是MatrixExpr也可以是ImmutableMatrix。相关类与函数位于 blockmatrix.py。6.1 构造与形状约束BlockMatrix用一个矩阵的矩阵blocks存储子块 from sympy import (MatrixSymbol, BlockMatrix, symbols, ... Identity, ZeroMatrix, block_collapse) n, m, l symbols(n m l) X MatrixSymbol(X, n, n) Y MatrixSymbol(Y, m, m) Z MatrixSymbol(Z, n, m) B BlockMatrix([[X, Z], [ZeroMatrix(m, n), Y]]) print(B) Matrix([ [X, Z], [0, Y]]) C BlockMatrix([[Identity(n), Z]]) print(C) Matrix([[I, Z]])构造器会做严格的分块规则性校验blockmatrix.py每行块数相同、同一行内各块的行数相同、同一列内各块的列数相同若不满足则抛出ValueError。文档特别演示了反例行内块数相同但每行总列数不一致的伪分块矩阵必须交给Matrix构造BlockMatrix会明确报错提示 from sympy import ones, Matrix dat [[ones(3,2), ones(3,3)*2], [ones(2,3)*3, ones(2,2)*4]] BlockMatrix(dat) ValueError: Although this matrix is comprised of blocks, ... Matrix(dat) Matrix([ [1, 1, 2, 2, 2], ... [3, 3, 3, 4, 4]])BlockMatrix提供若干结构属性blockshape块的排列形状、rowblocksizes/colblocksizes各块行/列尺寸、is_structurally_symmetric等shape通过对块尺寸求和得到blockmatrix.py。6.2 block_collapse分块表达式的化简引擎block_collapse(expr)自底向上遍历表达式对每个含BlockMatrix的子表达式应用化简规则blockmatrix.py块与块相乘时逐块执行块乘法bc_matmul借助_blockmul块与块相加时逐块相加bc_mataddBlockMatrix若退化为 1×1 块则解包bc_unpack嵌套的块矩阵被展平deblock转置、求逆按块结构处理bc_transpose、bc_inverse求逆支持 1×1 与 2×2 分块求逆公式blockinverse_1x1/blockinverse_2x2后者依据 Lu–Shiou 的 2×2 分块求逆公式并借助_choose_2x2_inversion_formula挑选可逆子块标量系数分发进各块bc_dist。文档示例 print(block_collapse(C*B)) Matrix([[X, Z Z*Y]])C是 1×2 块行矩阵B是 2×2 块矩阵两者相乘后坍缩为单行块矩阵。相关测试见 test_blockmatrix.py如test_block_collapse_explicit_matrices显式稠密/稀疏矩阵进入块矩阵后能被block_collapse还原与test_issue_17624幂等块矩阵的幂坍缩为BlockMatrix([[a**2, z], [z, z]])。6.3 BlockDiagMatrix 与 blockcutBlockDiagMatrix(X, Y, ...)构造对角分块矩阵非对角位置自动填充ZeroMatrixblockmatrix.py BlockDiagMatrix(X, Y) Matrix([ [X, 0], [0, Y]])其求逆、转置、行列式都按逐块方式计算_eval_inverse()对所有块求逆后拼回对角块矩阵_eval_determinant()等于各块行列式之积仅当所有块为方阵否则整体秩亏、行列式为 0。对角线块可通过get_diag_blocks()取回。blockcut(expr, rowsizes, colsizes)则相反——按给定的行/列尺寸把一个矩阵表达式切割成BlockMatrix内部使用MatrixSliceblockmatrix.py from sympy import ImmutableMatrix, blockcut M ImmutableMatrix(4, 4, range(16)) B blockcut(M, (1, 3), (1, 3)) type(B).__name__ BlockMatrix ImmutableMatrix(B.blocks[0, 1]) Matrix([[1, 2, 3]])6.4 分块矩阵的 Schur 补与分块分解对 2×2 分块矩阵[[A, B], [C, D]]BlockMatrix提供成套的分块算法blockmatrix.pyschur(matA, generalizedFalse)Schur 补。默认D - C*A⁻¹*B指定matD得A - B*D⁻¹*C。子块不可逆时可用generalizedTrue走 Moore-Penrose 广义逆如X.schur(B, generalizedTrue)返回C - D*(B.T*B)**(-1)*B.T*A。非 2×2 块矩阵抛ShapeErrorLDUdecomposition()/UDLdecomposition()/LUdecomposition()分块 LDU、UDL、LU 分解要求对应子块可逆A或D奇异时抛NonInvertibleMatrixError并可用block_collapse(L*D*U)还原原矩阵验证。例如文档式测试test_blockmatrix.py验证在假设Q.invertible(A)成立的上下文中det(X) det(A) * det(X.schur(A))——即分块行列式的 Schur 补公式。七、BlockMatrix 内部的转置、求逆与迹分块结构上的运算无需逐元素展开转置先转置每个子块再转置块结构本身即BlockMatrix([[A,B],[C,D]]).T→[[A.T, 0], [Z.T, Y.T]]_eval_transposeblockmatrix.py迹仅当分块结构对称rowblocksizes colblocksizes时trace Σ trace(对角块)否则保持未求值_eval_trace行列式1×1 退化为子块行列式2×2 时若A或D已知可逆使用det(A)*det(D - C*A.I*B)或det(D)*det(A - B*D.I*C)否则保持Determinant(BlockMatrix)未求值_eval_determinantblockmatrix.py。这些规则与前面block_collapse的分块求逆逻辑相互配合构成尽可能在块级完成运算、避免展开到元素级的设计。八、配套工具与打印matrix_symbols(expr)提取表达式中的矩阵符号列表matexpr.pyMatrixExpr.from_index_summation(expr)把显式带指标求和Sum(A[i,j]*B[j,k], (j,0,N-1))反向压缩成矩阵表达式A*B还能识别转置A.T*B与迹Trace(A)matexpr.py底层借助sympy.tensor.array.expressions的convert_indexed_to_array与convert_array_to_matrix.applyfunc(func)逐元素应用标量函数返回ElementwiseApplyFunctionmatexpr.py.T、.I、.Hadjoint()、.det()、.inv()等便捷属性/方法定义于 matexpr.py数值互操作__array__支持转成 NumPy 数组dtypeobjectequals()支持与显式矩阵的逐元素比较matexpr.py。打印方面HadamardProduct在终端中显示为A.*BHadamardPower显示为A.^3样式见 hadamard.py 的 docstring 输出示例。九、最佳实践与注意事项符号形状不可显式化MatrixSymbol(A, n, n)n 为符号无法as_explicit()但可参与转置、乘法、求逆、block_collapse等全部符号运算——这正是符号推导的价值所在。求逆是惰性的Inverse(A)与A**(-1)只记录逆这一结构真正求数值逆需.doit()或对显式矩阵调用.inv()inverse.py。在assuming(Q.orthogonal(X))上下文中refine(X.I)能化简为X.Trefine_Inverseinverse.py。维度校验MatrixSymbol、ZeroMatrix、Identity等构造时通过_check_dim强制维度为非负整数matexpr.py乘法、加法与 Hadamard 积均有形状校验。不要混用标量与矩阵的歧义运算MatMul要求非交换系数标量系数会通过factor_in_front提取到最前matmul.py。大矩阵优先用表达式像FunctionMatrix(1000, 1000, ...)这类稠密但有规律的矩阵应保持表达式并依赖惰性求值(Y**2)[10,10]避免展开出百万级元素的稠密矩阵。十、延伸阅读矩阵表达式求导与数组表达式的衔接sympy/tensor/array/expressions 下的from_matrix_to_array.py、arrayexpr_derivatives.py、from_array_to_matrix.py分块矩阵的完整测试规范test_blockmatrix.py各表达式对象的逐模块测试test_matmul.py、test_inverse.py、test_hadamard.py、test_derivatives.py、test_trace.py、test_sets.py等均位于 sympy/matrices/expressions/tests更广义的矩阵模块文档与源码sympy.matrices下的matrixbase.py、dense.py、sparse.py提供显式矩阵的对应实现。【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考