Symbolic Matrices in Mathematica with unknown dime

2019-02-01 08:30发布

站内文章 / 前端开发

30 0

爷、活的狠高调

女 | 书童

私信

可以将文章内容翻译成中文,广告屏蔽插件可能会导致该功能失效(如失效，请关闭广告屏蔽插件后再试):

问题:

Is there a way to do symbolic matrix algebra in Mathematica for matrices where the dimensions are unknown? For example, if I have an MxL matrix A and an LxN matrix B, I would like to be able to enter

A.B

And have it give me a matrix whose element ab[i,j] is given by

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

The problem I'm working on is like this one, but involves the product of 12 matrices, including the same matrix (and its transpose) repeated several times. It will probably be possible to simplify the values of the resulting matrix, but it's not obvious whether this is possible until after I do the algebra. This may be a problem that I have to solve by hand, but it would be much easier if Mathematica could provide some help in simplifying the algebra.

回答1:

Here's the code [dead-link removed] that wasted my morning... It's not complete, but it basically works. You can either get the notebook from the previous link [dead] or copy the code below.

Note that a similar question turned up on ask.sagemath not so long ago.

Almost like Sasha's solution, you define a symbolic matrix using

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

for some string "A" that does not have to be the same as the symbol A. Ok, here's the code:

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

Short hand for entering square matrices (could make it work for different heads...)

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

Behavior under Transpose, Conjugate, ConjugateTranspose and 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

Combining matrix powers (including the identity matrix)

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

Pretty printing with the dimension as a tooltip.

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

Define some rules for Dot. Need to extend then so that it can handle scalar quantities etc... Also so that inverses of A.B can be taken if A.B is square, even if neither A nor B are square.

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

Make Transpose, Conjugate and ConjugateTranspose distribute over 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];

Here's some simple tests/examples

Now for code that deals with the component expansion. Like Sasha's solution, I'll overload Part.

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

Some notation

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

Code to extract parts of matrices and Dot products of matrices Need to add checks to ensure that explicit summation ranges are all sensible.

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

Some examples:

回答2:

I am not sure if this is really very helpful, but it could be a start:

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]

Now, define some symbolic matrices and do dot product:

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

Evaluate matrix elements:

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

回答3:

If you're willing to switch from Mathematica to Python the functionality you need is in the development branch of SymPy. It should be in the 0.72 release.

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₂₂)

These purely symbolic objects can also make use of all of the standard matrix algorithms for explicitly defined matrices

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

Symbolic shaped matrices are doable too

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:

I use this approach:

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

When invoked as A = SymbolicMatrix["a", 2, 3]; or A = SymbolicMatrix[a, 2, 3];, it creates a matrix

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

So, it creates mxn symbols, but I find them descriptive, and the whole thing is easy to use (at least for my purposes).

回答5:

You can use NCAlgebra for that. For example:

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

will define symbolic matrices AA and BB with noncommutative entries $a_{i,j}$ and $b_{i,j}$. These can get manipulated and produce the results you're looking for. For example: