【问题标题】:Symbolic Matrices in Mathematica with unknown dimensionsMathematica 中未知维度的符号矩阵
【发布时间】:2011-08-08 04:04:29
【问题描述】:

有没有办法在 Mathematica 中为维度未知的矩阵做符号矩阵代数?例如,如果我有一个 MxL 矩阵 A 和一个 LxN 矩阵 B,我希望能够输入

A.B

让它给我一个矩阵,其元素 ab[i,j] 由下式给出

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

我正在处理的问题与此类似,但涉及 12 个矩阵的乘积,包括重复多次的相同矩阵(及其转置)。可能有可能简化结果矩阵的值,但直到我做代数之后,这是否可能并不明显。这可能是一个我必须手动解决的问题,但如果 Mathematica 可以在简化代数方面提供一些帮助会容易得多。

【问题讨论】:

  • 作为 Sasha pointed out,这是可行的。但是,我会这样做,并依靠 Mathematica 来检查它。不过,我会使用Einstein summation convention,因为它可以节省大量时间,并且可以让大多数操作变得更容易。
  • @rcollyer 我想我在你的评论中遗漏了一些东西。使用索引收缩是什么意思?如果您的解释超出了评论的范围,我可以发布一个临时问题
  • @belisarius:我认为@rcollyer 指的是在手动计算时使用较短的符号,因为它更方便;不在mma。例如,a_1b_1+...+a_nb_n 可以方便地写成a_ib_i,其中暗示了所有n 的总和。点积、二元积、迹线等都可以这样紧凑地编写。当涉及到指数时,事情会变得更加混乱,因为符号会变得混乱(我不确定他们是否在这种情况下甚至使用这种符号)。
  • @belisarius:通过谷歌快速搜索,我找到了this package,它可以让您使用爱因斯坦求和约定执行符号张量计算。我还没有测试过,但如果它有任何好处,应该将它添加到 Mathematica 工具上的 CW 帖子中。
  • 目前它对我们没有帮助,但 Stephen Wolfram 已经承诺Real Soon Now 即将推出一项功能,即“让进行严肃的张量分析感觉就像进行普通代数一样”。我很想知道他的意思是什么,以及该功能是否会解决这个问题提出的问题。

标签: wolfram-mathematica


【解决方案1】:

我不确定这是否真的很有帮助,但这可能是一个开始:

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}]

【讨论】:

  • 我看到了使用InOut 的缺点,它使复制变得更加困难。您是否介意将定义放在输入的“一行”中,以便我们更轻松地访问它?
  • +1 我喜欢这种方法,但您需要定义大量矩阵函数,例如转置。此外,对于您定义的每个符号矩阵,您还必须提供一个函数名称(第 109 行中的 a,b)。用户最好不要为此烦恼。
  • ToMatrix 旨在为给定的显式维度实例化一个矩阵,但这并不好,所以我在编辑中删除了它。
  • 非常感谢您的帮助,@Sasha。我不是一个非常优秀的 Mathematica 程序员,所以在我理解你的代码之前,我需要花一些时间阅读文档(我可能会阅读模块、范围等)同时,我会可能只是遵循@rcollyer 的建议并手动解决。
【解决方案2】:

这是浪费我早上的代码 [dead-link removed]... 它不完整,但它基本上可以工作。您可以从上一个链接 [dead] 中获取 notebook,也可以复制下面的代码。

请注意,不久前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]

Transpose、Conjugate、ConjugateTranspose 和 Inverse 下的行为

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];

使转置、共轭和共轭转置分布在 Plus 上。

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];

一些例子:

【讨论】:

  • 这超出了我的想象,但看起来很花哨,所以:+1
  • @Mr.Wizard:我怀疑这是否超出了您的想象。这只是一堆乱七八糟的Option* 命令……
  • 哇,你应该在 Mathematica Tool Bag 社区 wiki stackoverflow.com/q/4198961/181759 中发布这个,非常棒的工作。
  • 谢谢,@西蒙!这非常有帮助。我不是一个非常好的 Mathematica 程序员,所以我不太明白发生了什么,但我几乎可以按原样使用你的代码。作为我试图评估的矩阵乘积的一部分,我需要反复乘以其条目全为 1 的方阵。我尝试(但失败)将其实现为 SymbolicMatrix 上的函数,因此每个条目都替换为该列元素的总和。您向我展示执行此操作的代码是否容易?再次感谢您的帮助 - 这启发了我学习更多 Mathematica!
  • @jack:目前我还没有决定实现实际矩阵替换的最佳方法。我想要一个简单的界面,你说 SymbolicMatrix -> Matrix,然后它将为你处理所有的转置、共轭、MatrixPowers 和组件表达式。我也没有太多时间玩它。目前,您可以使用显式替换规则。 (对不起...我最终会写代码!)
【解决方案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                

【讨论】:

【解决方案4】:

我使用这种方法:

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 符号,但我发现它们是描述性的,而且整个东西很容易使用(至少对我来说是这样)。

【讨论】:

    【解决方案5】:

    您可以为此使用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}$ 定义符号矩阵 AABB。这些可以被操纵并产生您正在寻找的结果。例如:

    NCDot[AA, BB]
    

    将使用** 将两个矩阵相乘,

    AA ** BB // NCMatrixExpand
    

    也会这样做,并且

    tp[AA ** BB] // NCMatrixExpand
    

    将使用tp进行转置。

    【讨论】:

      猜你喜欢
      • 2013-05-31
      • 2012-04-26
      • 1970-01-01
      • 2014-06-04
      • 1970-01-01
      • 2010-12-05
      • 2016-09-10
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多