Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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 \[Element] Reals, x, Transpose[x]]

then

Simplify[g2@Simplify[Transpose[a*A1], TransformationFunctions -> {Automatic, g1}, 
   ComplexityFunction -> cpfunc], 
   Assumptions -> a \[Element] Reals && A1 \[NotElement] 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\[Element] Reals&&A1\[NotElement]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 \[Element] Reals, a2 \[Element] Reals, 
 a3 \[Element] Reals, a4 \[Element] Reals, A1 \[NotElement] Reals, 
 A2 \[NotElement] Reals, A3 \[NotElement] Reals, 
 A4 \[NotElement] Reals, A5 \[NotElement] Reals, 
 A6 \[NotElement] 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.

Thank you very much for the help.

share|improve this question
1  
Just checking : the second identity holds in your specific problem, not in general, right ? –  b.gatessucks Feb 22 '13 at 9:27
    
Hi @b.gatessucks, yes it's a special case, I'll make it clear in my question, thanks for point it out. –  xslittlegrass Feb 22 '13 at 13:49
2  
You might find useful our Quantum Mathematica add-on (it is free) for doing these kind of operations, see for example homepage.cem.itesm.mx/jose.luis.gomez/quantum/v7algebra.nb and homepage.cem.itesm.mx/jose.luis.gomez/quantum/v7triangles.pdf I hope this helps. –  MaTECmatica Feb 26 '13 at 23:54

1 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[_] | ...
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.