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: