Symbolic Matrices in Mathematica with unknown dimensions

Question

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.

Answer 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:

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

Answer 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                

Answer 4

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:

NCDot[AA, BB]

will multiply the two matrices using **,

AA ** BB // NCMatrixExpand

will do the same, and

tp[AA ** BB] // NCMatrixExpand

will transpose using tp.

Answer 5

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