Mathematica中具有未知维度的符号矩阵

Question

在Mathematica中，是否有一种方法可以对维数未知的矩阵进行符号矩阵代数运算？例如，如果我有一个MxL矩阵A和一个LxN矩阵B，我想输入

A.B

并且让它给我一个矩阵，其元素ab[i,j]由以下内容给出

Sum[a[i,l]*b[l,j],{l,1,L}]

我正在处理的问题与这个问题类似，但涉及到 12 个矩阵的乘积，其中包括相同的矩阵（及其转置）被多次重复使用。可能可以简化结果矩阵的值，但在进行代数运算之前，不确定是否有这种可能性。这可能是一个必须手动解决的问题，但如果Mathematica能够提供一些帮助来简化代数运算，那将更容易些。

- jack

---

@belisarius：我认为@rcollyer指的是在手动计算时使用较短的符号更方便，而不是在 MMA 中。例如， a_1b_1+...+a_nb_n 可以方便地写成 a_ib_i，其中隐含着对所有 n 的求和。点积、双重积、迹等都可以用这种方式简洁地表示。当涉及到指数时，事情会变得更混乱，因为符号会变得令人困惑（我不确定在这种情况下他们是否使用这种符号）。 - user616736

@yoda 我错过了“手动”部分。所以问题是：有没有办法在Mma中使用索引缩并？ - Dr. belisarius

@belisarius：一个快速的谷歌搜索让我找到了这个包，它允许你使用爱因斯坦求和约定执行符号张量计算。我没有测试过，但如果它不错，就应该将其添加到关于Mathematica工具的CW帖子中。 - user616736

@yoda 不错的发现！然而，由于它自2006年以来没有更新，我们应该预料到兼容性问题... - Dr. belisarius

目前这并不能帮助我们，但是Stephen Wolfram承诺将很快推出一个功能，Real Soon Now，它将使得处理张量分析的感觉几乎像处理普通代数一样简单。我非常期待看到他的意思是什么，并且这个功能是否能解决像这个问题所提出的那样的问题。 - WReach

Answer 1

这是浪费我早上时间的代码 [失效链接已移除]... 它不完整，但基本上可以工作。你可以从之前的链接 [失效] 获取笔记本或复制下面的代码。

请注意，不久前在ask.sagemath上出现了一个类似的问题。

与Sasha的解决方案几乎相同，使用符号矩阵定义

A = SymbolicMatrix["A", {n, k}]

对于一些字符串"A"，它不必与符号A相同。好的，这是代码：

ClearAll[SymbolicMatrix]
Options[SymbolicMatrix] = {Transpose -> False, Conjugate -> False, MatrixPower -> 1};

输入方阵的简写形式（也可以适用于不同类型的矩阵...）

SymbolicMatrix[name_String, n:_Symbol|_Integer, opts : OptionsPattern[]] := SymbolicMatrix[name, {n, n}, opts]

转置、共轭、共轭转置和逆运算下的行为

SymbolicMatrix/:Transpose[SymbolicMatrix[name_String,{m_,n_},opts:OptionsPattern[]]]:=SymbolicMatrix[name,{n,m},
  Transpose->!OptionValue[SymbolicMatrix,Transpose],Sequence@@FilterRules[{opts},Except[Transpose]]]
SymbolicMatrix/:Conjugate[SymbolicMatrix[name_String,{m_,n_},opts:OptionsPattern[]]]:=SymbolicMatrix[name,{m,n},
  Conjugate->!OptionValue[SymbolicMatrix,Conjugate],Sequence@@FilterRules[{opts},Except[Conjugate]]]
SymbolicMatrix/:ConjugateTranspose[A:SymbolicMatrix[name_String,{m_,n_},opts:OptionsPattern[]]]:=Conjugate[Transpose[A]]
SymbolicMatrix/:Inverse[SymbolicMatrix[name_String,{n_,n_},opts:OptionsPattern[]]]:=SymbolicMatrix[name,{n,n},
  MatrixPower->-OptionValue[SymbolicMatrix,MatrixPower],Sequence@@FilterRules[{opts},Except[MatrixPower]]]

SymbolicMatrix/:(Transpose|Conjugate|ConjugateTranspose|Inverse)[eye:SymbolicMatrix[IdentityMatrix,{n_,n_}]]:=eye

组合矩阵幂（包括单位矩阵）

SymbolicMatrix/:SymbolicMatrix[a_String,{n_,n_},opt1:OptionsPattern[]].SymbolicMatrix[a_,{n_,n_},opt2:OptionsPattern[]]:=SymbolicMatrix[a,{n,n},Sequence@@FilterRules[{opt1},Except[MatrixPower]],MatrixPower->Total[OptionValue[SymbolicMatrix,#,MatrixPower]&/@{{opt1},{opt2}}]]/;FilterRules[{opt1},Except[MatrixPower]]==FilterRules[{opt2},Except[MatrixPower]]

SymbolicMatrix[a_String,{n_,n_},opts:OptionsPattern[]]:=SymbolicMatrix[IdentityMatrix,{n,n}]/;OptionValue[SymbolicMatrix,{opts},MatrixPower]===0

SymbolicMatrix/:(A:SymbolicMatrix[a_String,{n_,m_},OptionsPattern[]]).SymbolicMatrix[IdentityMatrix,{m_,m_}]:=A
SymbolicMatrix/:SymbolicMatrix[IdentityMatrix,{n_,n_}].(A:SymbolicMatrix[a_String,{n_,m_},OptionsPattern[]]):=A

使用尺寸作为工具提示的漂亮打印。

Format[SymbolicMatrix[name_String,{m_,n_},opts:OptionsPattern[]]]:=With[{
  base=If[OptionValue[SymbolicMatrix,MatrixPower]===1,
    StyleBox[name,FontWeight->Bold,FontColor->Darker@Brown],
    SuperscriptBox[StyleBox[name,FontWeight->Bold,FontColor->Darker@Brown],OptionValue[SymbolicMatrix,MatrixPower]]],
  c=Which[
    OptionValue[SymbolicMatrix,Transpose]&&OptionValue[SymbolicMatrix,Conjugate],"\[ConjugateTranspose]",
    OptionValue[SymbolicMatrix,Transpose],"\[Transpose]",
    OptionValue[SymbolicMatrix,Conjugate],"\[Conjugate]",
  True,Null]},
  Interpretation[Tooltip[DisplayForm@RowBox[{base,c}/.Null->Sequence[]],{m,n}],SymbolicMatrix[name,{m,n},opts]]]

Format[SymbolicMatrix[IdentityMatrix,{n_,n_}]]:=Interpretation[Tooltip[Style[\[ScriptCapitalI],Bold,Darker@Brown],n],SymbolicMatrix[IdentityMatrix,{n,n}]]

定义一些规则来处理 Dot。需要扩展它，使其能够处理标量等。同时，即使 A 或 B 不是方阵，如果 A.B 是方阵，则可以取 A.B 的逆。

SymbolicMatrix::dotdims = "The dimensions of `1` and `2` are not compatible";
Unprotect[Dot]; (*Clear[Dot];*)
Dot/:(a:SymbolicMatrix[_,{_,n_},___]).(b:SymbolicMatrix[_,{m_,_},___]):=(Message[SymbolicMatrix::dotdims,HoldForm[a],HoldForm[b]];Hold[a.b])/;Not[m===n]
Dot/:Conjugate[d:Dot[A_SymbolicMatrix,B__SymbolicMatrix]]:=Map[Conjugate,d]
Dot/:(t:Transpose|ConjugateTranspose)[d:Dot[A_SymbolicMatrix,B__SymbolicMatrix]]:=Dot@@Map[t,Reverse[List@@d]]
Dot/:Inverse[HoldPattern[d:Dot[SymbolicMatrix[_,{n_,n_},___]...]]]:=Reverse@Map[Inverse,d]
A_ .(B_+C__):=A.B+A.Plus[C]
(B_+C__).A_:=B.A+Plus[C].A
Protect[Dot];

使转置、共轭和共轭转置在加法运算中具有分配律。

Unprotect[Transpose, Conjugate, ConjugateTranspose];
Clear[Transpose, Conjugate, ConjugateTranspose];
Do[With[{c = c}, c[p : Plus[a_, b__]] := c /@ p], {c, {Transpose, Conjugate, ConjugateTranspose}}]
Protect[Transpose, Conjugate, ConjugateTranspose];

这里有一些简单的测试/示例

现在讨论处理组件扩展的代码。像Sasha的解决方案一样，我会重载Part函数。

Clear[SymbolicMatrixComponent]
Options[SymbolicMatrixComponent]={Conjugate->False,MatrixPower->1};

一些符号

Format[SymbolicMatrixComponent[A_String,{i_,j_},opts:OptionsPattern[]]]:=Interpretation[DisplayForm[SubsuperscriptBox[StyleBox[A,Darker@Brown],RowBox[{i,",",j}],
RowBox[{If[OptionValue[SymbolicMatrixComponent,{opts},MatrixPower]===1,Null,OptionValue[SymbolicMatrixComponent,{opts},MatrixPower]],If[OptionValue[SymbolicMatrixComponent,{opts},Conjugate],"*",Null]}/.Null->Sequence[]]]],
SymbolicMatrixComponent[A,{i,j},opts]]

代码提取矩阵的部分和矩阵的Dot乘积

需要添加检查以确保显式求和范围都是合理的。

SymbolicMatrix/:SymbolicMatrix[A_String,{m_,n_},opts:OptionsPattern[]][[i_,j_]]:=SymbolicMatrixComponent[A,If[OptionValue[SymbolicMatrix,{opts},Transpose],Reverse,Identity]@{i,j},Sequence@@FilterRules[{opts},Options[SymbolicMatrixComponent]]]

SymbolicMatrix/:SymbolicMatrix[IdentityMatrix,{m_,n_}][[i_,j_]]:=KroneckerDelta[i,j]

Unprotect[Part]; (*Clear[Part]*)
Part/:((c___.b:SymbolicMatrix[_,{o_,n_},OptionsPattern[]]).SymbolicMatrix[A_String,{n_,m_},opts:OptionsPattern[]])[[i_,j_]]:=With[{s=Unique["i",Temporary]},Sum[(c.b)[[i,s]]SymbolicMatrixComponent[A,If[OptionValue[SymbolicMatrix,{opts},Transpose],Reverse,Identity]@{s,j},Sequence @@ FilterRules[{opts}, Options[SymbolicMatrixComponent]]],{s,n}]]
Part/:(a_+b_)[[i_,j_]]:=a[[i,j]]+b[[i,j]]/;!And@@(FreeQ[#,SymbolicMatrix]&/@{a,b})
Part/:Hold[a_][[i_,j_]]:=Hold[a[[i,j]]]/;!FreeQ[a,SymbolicMatrix]
Protect[Part];

一些例子：

- Simon

---

这对我来说有些难度，但看起来很高级，所以：+1 - Mr.Wizard

@Mr.Wizard：我怀疑这些内容是否超出了你的理解范围。它只是一堆Option*命令的混乱。 - Simon

哇，你应该把这个发布到Mathematica工具包社区维基https://dev59.com/nuo6XIcBkEYKwwoYTS1D，非常棒的工作。 - Timo

@jack：目前我还没有决定实现实际矩阵替换的最佳方法。我想要一个简单的接口，你只需要说SymbolicMatrix -> Matrix，它就会处理所有的转置、共轭、矩阵幂和组成表达式。我也没有太多时间去玩它。暂时你可以使用显式替换规则。（抱歉...我最终会写出代码！） - Simon

这是非常有用的功能。我想知道您是否可以添加以下可能性：1）创建具有符号子矩阵的实际矩阵，例如 Mi = {{1,Transpose[vi]},{vi,Ai}};（其中 vi 和 Ai 是符号），使得取 M1.M2 将执行 2x2 矩阵乘法并将 . 运算符放在其中的符号元素之间。2）创建隐含依赖于一个或多个变量的符号矩阵，例如 A[u]，以便我们在进行导数计算时获得正确的结果，而不是 D[A[u],u]=0。这将非常棒！ - Kagaratsch

Answer 2

我不确定这是否真的有帮助，但它可能是一个开始：

ClearAll[SymbolicMatrix]
SymbolicMatrix /: Transpose[SymbolicMatrix[a_, {m_, n_}]] := 
     SymbolicMatrix[Evaluate[a[#2, #1]] & , {n, m}]
SymbolicMatrix /: 
 SymbolicMatrix[a_, {m_, n_}] . SymbolicMatrix[b_, {n_, p_}] := 
     With[{s = Unique[\[FormalI], Temporary]}, 
  SymbolicMatrix[Function[{\[FormalN], \[FormalM]}, 

    Evaluate[Sum[a[\[FormalN], s]*b[s, \[FormalM]], {s, 1, n}]]], {m, 
    p}]]
SymbolicMatrix /: SymbolicMatrix[a_, {m_, n_}][[i_, j_]] := a[i, j]

---

现在，定义一些符号矩阵并进行点乘：

In[109]:= amat = SymbolicMatrix[a, {n, k}]; 
bmat = SymbolicMatrix[b, {k, k}]; 

评估矩阵元素：

In[111]:= (amat . bmat . Transpose[amat])[[i, j]]

Out[111]= Sum[
 a[j, \[FormalI]$1485]*
  Sum[a[i, \[FormalI]$1464]*
    b[\[FormalI]$1464, \[FormalI]$1485], {\[FormalI]$1464, 1, k}], 
   {\[FormalI]$1485, 1, k}]

- Sasha

---

我看到使用 In 和 Out 的一个劣势，它使复制变得困难得多。你介意将定义放在一行输入中，以便我们更轻松地访问它吗？ - rcollyer

+1 我喜欢这种方法，但你需要定义大量的矩阵函数，例如转置。此外，对于每个符号矩阵的定义，你还需要提供一个函数名称（第109行中的a、b）。最好不要让用户受到干扰。 - Sjoerd C. de Vries

ToMatrix原本是用来根据给定的显式维度实例化矩阵的，但它并不好用，所以我在编辑中将其删除了。 - Sasha

非常感谢您的帮助，@Sasha。我不是一个特别擅长Mathematica编程的人，所以在理解您的代码之前，我需要花费一些时间仔细阅读文档（我可能会阅读有关模块、作用域等方面的内容）。与此同时，我可能会遵循@rcollyer的建议，手动解决问题。 - jack

Answer 3

如果您愿意从Mathematica转换到Python，您所需的功能在SymPy的开发分支中。它应该在0.72版本中。

In [1]: from sympy import *
In [2]: X = MatrixSymbol('X', 2,3)
In [3]: Y = MatrixSymbol('Y', 3,3)
In [4]: X*Y*X.T
Out[4]: X⋅Y⋅X'

In [5]: (X*Y*X.T)[0,1]
Out[5]: 
X₀₀⋅(X₁₀⋅Y₀₀ + X₁₁⋅Y₀₁ + X₁₂⋅Y₀₂) + X₀₁⋅(X₁₀⋅Y₁₀ + X₁₁⋅Y₁₁ + X₁₂⋅Y₁₂) + X₀₂⋅(X₁₀⋅Y₂₀ + X₁₁⋅Y₂₁ + X₁₂⋅Y₂₂)

这些纯符号对象也可以利用所有标准矩阵算法来明确定义矩阵

In [14]: X = MatrixSymbol('X', 2,2)
In [14]: X.as_explicit().det()
Out[14]: X₀₀⋅X₁₁ - X₀₁⋅X₁₀

符号形状的矩阵也是可以实现的。

In [7]: n,m,k = symbols('n,m,k')
In [8]: X = MatrixSymbol('X', n,m)
In [9]: Y = MatrixSymbol('Y', m,k)
In [10]: (X*Y)[3,4]
Out[10]: 
m - 1                
 ___                 
 ╲                   
  ╲   X(3, k)⋅Y(k, 4)
  ╱                  
 ╱                   
 ‾‾‾                 
k = 0                

- MRocklin

---

不幸的是，sympy 似乎无法通过符号矩阵进行通用代数操作。例如，您无法使用 solve 重新排列表达式：https://github.com/sympy/sympy/issues/11124 - mpeac

Answer 4

你可以使用 NCAlgebra 进行操作。例如：

MM = 3
LL = 2
NN = 3
AA = Table[Subscript[a, i, j], {i, 1, MM}, {j, 1, LL}]
BB = Table[Subscript[b, i, j], {i, 1, LL}, {j, 1, NN}]

将定义具有非交换条目 $a_{i,j}$ 和 $b_{i,j}$ 的符号矩阵 AA 和 BB。这些可以被操作并产生您寻找的结果。例如：

NCDot[AA, BB]

将使用**乘以这两个矩阵。

AA ** BB // NCMatrixExpand

将会做同样的事情，而且。

tp[AA ** BB] // NCMatrixExpand

将使用 tp 进行转置。

- Mauricio de Oliveira

Answer 5

我使用这种方法：

SymbolicMatrix[symbol_String, m_Integer, n_Integer] := Table[
  ToExpression[symbol <> ToString[i] <> ToString[j]],
  {i, 1, m}, {j, 1, n}
];
SymbolicMatrix[symbol_Symbol, m_Integer, n_Integer] := SymbolicMatrix[ToString[symbol], m, n];
SymbolicMatrix[symbol_, m_Integer] := SymbolicMatrix[symbol, m, 1];

当以A = SymbolicMatrix["a", 2, 3];或A = SymbolicMatrix[a, 2, 3];调用时，它会创建一个矩阵。

{{a11, a12, a13}, {a21, a22, a23}}

所以它创建了mxn个符号，但我发现它们很具描述性，整个东西很容易使用（至少对于我的目的而言）。

- Vedran Šego