Simplify matrix algebra

Question

I'm trying to simplify some matrix linear algebra, for example, simplify

$$\big(a1\times(A1\cdot A2)\big)\cdot\Big(a2\times A3\cdot A4+(a3\times A5)\cdot(a4\times A6)\Big)^{T}$$

where lower case variables $(a1-a4)$ are numbers and upper case variables $(A1-A6)$ are matrixes. "$\times$" is the Times function and "$\cdot$" is the Dot function, $\text{T}$ is the transpose operation.

I have these math identities for my matrixes in my case:

$$ (A\cdot B)^{T}=B^{T}\cdot A^{T}\\ (A+B)^{T}=A^{T}+B^{T}\\ c\times(A+B)=c\times A+c\times B\\ c^{T}=c $$

so the expresion will simplify into $$ a1~a2~A1\cdot A2\cdot A4^{T}\cdot A3^{T}+a1~a3~a4(A1\cdot A2\cdot A6^{T} \cdot A5^{T}) $$

My major question is how can I use Simplify or FullSimplify to do this simplification in MMA in an elegant and extensible way? For extensible I mean that one can easily add more operation rules into the simplification identities, for example, we could add inverse operation and its identities such as $(A\cdot B)^{-1}=B^{-1}\cdot A^{-1}$ and $(A^{-1})^T=(A^T)^{-1}$ etc.

---

Here is my try (I'm a novice and apparently my code is not elegant in any sense so please don't hesitate to advice suggestions or to point out any mistakes I made, I'm very welling to learn from you guys. And please pardon my poor English :P):

1.First simplify $(A\cdot B)^{T}$ to $B^{T}\cdot A^{T}$ and $(A+B)^{T}$ to $A^{T}+B^{T}$, I define a transformation function

g1[expr_] := expr /. Transpose[A_] :> If[Head[A] == Dot || Head[A] == Times || Head[A] == Plus,
 Head[A] @@ Reverse[Transpose /@ A], Transpose[A]]

which changes Transpose[Dot[A,B]] to Dot[Transpose[B],Transpose[A]], and Transpose[Plus[A,B]] to Plus[Transpose[B],Transpose[A]]. And then define a ComplexFunction that prefers the form that Transpose is in the inside, for example, prefer Dot[Transpose[B],Transpose[A]] to Transpose[Dot[A,B]] :

cpfunc[expr_] := Module[{funcls},
 If[Length[expr] > 1, 0, 
  If[Depth[expr] <= 2, 0,
   funcls = Reverse[Part[#, 0] & /@ (Level[#, {0, Depth[#] - 2}])] &@expr;
   If[Max[Position[funcls, Transpose]] - Max[Position[funcls, Dot]] < 0 ||
      Max[Position[funcls, Transpose]] - Max[Position[funcls, Times]] < 0 || 
      Max[Position[funcls, Transpose]] - Max[Position[funcls, Plus]] < 0,
   10^3, 0] (*gives punishment when found Dot inside Transpose, etc.*)
]]]

This works for Transpose[A1.B1.C1]:

Simplify[Transpose[A1.B1.C1], TransformationFunctions -> {Automatic, g1},
  ComplexityFunction -> cpfunc]
(*Transpose[C1].Transpose[B1].Transpose[A1]*)

Simplify[Transpose[a*A1], TransformationFunctions -> {Automatic, g1}, 
  ComplexityFunction -> cpfunc]
(*Transpose[a] Transpose[A1]*)

2.Second simplify $(a\times A)^{T}$ to $a~A^{T}$ : I define a transformation function:

g2[expr_] := 
 expr /. Transpose[x_] :> If[x ∈ Reals, x, Transpose[x]]

then

Simplify[g2@Simplify[Transpose[a*A1], TransformationFunctions -> {Automatic, g1}, 
   ComplexityFunction -> cpfunc], 
   Assumptions -> a ∈ Reals && A1 ∉ Reals]
(*a Transpose[A1]*)

Initially I was trying to do something like this but it doesn't work(could you point out the mistakes?)

g2[expr_] := expr /. (Transpose[x_];/NumberQ[x]) :> x
Simplify[Transpose[a*A],TransformationFunctions->{Automatic,g1,g2},
   ComplexityFunction->cpfunc,Assumptions->a∈ Reals&&A1∉Reals]
(*Transpose[a] Transpose[A]*)

3.Simplify the whole expression:

Simplify[(a1*A1.A2).Transpose[a2*A3.A4 + (a3*A5).(a4*A6)], 
 TransformationFunctions -> {Automatic, g1},
 ComplexityFunction -> cpfunc]
(*
 (a1 A1.A2).((Transpose[a4] Transpose[A6]).(Transpose[a3] Transpose[A5]) 
 + Transpose[A4].Transpose[A3] Transpose[a2])
*)

Simplify[g2@%, Assumptions -> {a1 ∈ Reals, a2 ∈ Reals, 
 a3 ∈ Reals, a4 ∈ Reals, A1 ∉ Reals, 
 A2 ∉ Reals, A3 ∉ Reals, 
 A4 ∉ Reals, A5 ∉ Reals, 
 A6 ∉ Reals}]
(*
 (a1 A1.A2).(a2 Transpose[A4].Transpose[A3] 
 + (a4 Transpose[A6]).(a3 Transpose[A5]))
*)

the output is $$ (\text{a1} \text{A1}\cdot\text{A2})\cdot(\text{a2} \text{A4}^{\mathsf{T}}\cdot\text{A3}^{\mathsf{T}}+(\text{a4} \text{A6}^{\mathsf{T}})\cdot(\text{a3} \text{A5}^{\mathsf{T}})) $$ This final answer is close to the result I want to achieve. I tried to go further but without any luck. Here is my thinking(questions) about this:

1. Is there some more elegant way to do this? For example, it is possible to set some kind of attributes to these different operations so that Dot, Times, Pulse etc have different priority in the calculation, and one can sort them directly by their attributes in order to simplify or expand an expression?

2. I'm very curious how MMA Expand an algebra expression like (a + b)*(c + d*e).

3. How to make a transformation function work on different assumptions? For example, when only x is a number, that the transformation Transpose[x_]:>x will apply.

Answer 1

I was really intrigued by this problem and I tried multiple solutions before posting one here. The most elegant approach I found is to use ReplaceRepeated with a custom list of scalar variables.

In order for it to work we need the proper rule set

rules := {
  T[Dot[A_, B__]] :> T[Dot[B]].T[A],
  T[Plus[A_, B__]] :> T[A] + T[Plus[B]],
  T[a_ A_] /; MemberQ[scalars, a] :> a T[A],
  T[a_] /; MemberQ[scalars, a] :> a,

  x__.(y_ + z__).w___ :> (x.y.w) + (x.z.w),
  (x___).(a_ y_).(z___) /; MemberQ[scalars, a] :> a (x.y.z)
  }

Let's see what they are doing. The first two,

  T[Dot[A_, B__]] :> T[Dot[B]].T[A],
  T[Plus[A_, B__]] :> T[A] + T[Plus[B]],

correspond exactly to your first two math identities. Using a sequence (B__) instead of a single expression (B_) allows to apply this rules one argument at the time. The third and the forth,

  T[a_ A_] /; MemberQ[scalars, a] :> a T[A],
  T[a_] /; MemberQ[scalars, a] :> a,

provide for some simplification whether a scalar is involved, and the last two,

  x__.(y_ + z__).w___ :> (x.y.w) + (x.z.w),
  (x___).(a_ y_).(z___) /; MemberQ[scalars, a] :> a (x.y.z)

implement some properties of the dot product.

This is the expression which has to be simplified:

expr = (a1*A1.A2).T[a2*A3.A4 + (a3*A5).(a4*A6)]

This the result without any assumptions:

In[484]:= expr //. rules // Expand

Out[484]= (a1 A1.A2).T[a2 A3.A4] + (a1 A1.A2).T[a4 A6].T[a3 A5]

And this providing an explicit list of scalars:

In[485]:= Block[
 {scalars = {a1, a2, a3, a4}},
 expr //. rules // Expand
 ]

Out[485]= a1 a2 A1.A2.T[A4].T[A3] + a1 a3 a4 A1.A2.T[A6].T[A5]

where I use Expand only to avoid the automatic collection of a1.

---

All of that can also be packed into a function

MatrixSimplify[expr_, s : _List : {}] :=
 Module[{},
  If[OwnValues[$scalars] === {},
       $scalars = {}
   ];
  Block[
   {$scalars = Join[$scalars, Flatten[{s}]]},

   With[{rules = {
       T[Dot[A_, B__]] :> T[Dot[B]].T[A],
       T[Plus[A_, B__]] :> T[A] + T[Plus[B]],
       T[a_ A_] /; MemberQ[$scalars, a] :> a T[A],
           T[a_] /; MemberQ[$scalars, a] :> a,

       x__.(y_ + z__).w___ :> (x.y.w) + (x.z.w),
       (x___).(a_ y_).(z___) /; MemberQ[$scalars, a] :> a (x.y.z)
       }},

    expr //. rules // Expand
    ]
   ]
  ]

MatrixSimplify must be provided with an expression and an optional list of scalars.

In[562]:= expr // MatrixSimplify

Out[562]= (a1 A1.A2).T[a2 A3.A4] + (a1 A1.A2).T[a4 A6].T[a3 A5]

In[563]:= MatrixSimplify[expr, {a4, a1}]

Out[563]= a1 A1.A2.T[a2 A3.A4] + a1 a4 A1.A2.T[A6].T[a3 A5]

It also relies on the global variable $scalars (the bit of code within the If statement is to check whether $scalar is initialized and, if not, it does set it to an empty list), so you can do something like

In[564]:= $scalars = {a1, a2, a3, a4};

In[565]:= expr // MatrixSimplify

Out[565]= a1 a2 A1.A2.T[A4].T[A3] + a1 a3 a4 A1.A2.T[A6].T[A5]

---

P.S.

An alternative way to do this is by using MatchQ instead of MemberQ. This makes it possible to specify the scalar variables with patterns as in

$scalars = a[_] | b[_] | ...