My goal is to algebraically simplify expressions involving tensor products and taking the trace.

That is, I would like to compute $\operatorname{Tr}\left( \prod_i(id+a_i \otimes b_i) \right)$,

using the rules of tensor algebra. E.g. $$\operatorname{Tr}(a_i \otimes b_i)=a_i.b_i$$ and $$a_1 \otimes b_1 a_2 \otimes b_2= b_1.a_2 a_1 \otimes b_2.$$

The final expression should be of the form

$$ \sum_i \prod_j a_{i_j}.b_{i_j}.$$

Any implementation to let Mathematica do this for me would be much appreciated.

So far I tried to use the built-in tensor algebra package

$Assumptions = {Element[a1, Vectors[3, Reals]], 
   Element[a2, Vectors[3, Reals]], Element[b1, Vectors[3, Reals]], 
   Element[b2, Vectors[3, Reals]], Element[n, Vectors[3, Reals]], 
   Element[m, Vectors[3, Reals]], Element[id, Matrices[{3, 3}]]};

and define simplification rules such as

SimId[expr_] := 
 expr //. {Dot[id, tensor___] :> tensor, Dot[tensor___, id] :> tensor,
    MatrixPower[id, 2] :> id, MatrixPower[id, 3] :> id, 
   MatrixPower[id, 4] :> id}
SimTP[expr_] := 
 expr //. {Dot[TensorProduct[before1_, before2_], 
     TensorProduct[after1_, after2_]] :> 
    Dot[before2, after1] TensorProduct[before1, after2]}

This works fine for simple expressions such as

  TensorProduct[a2 + b1, 
     n1].(id + TensorProduct[a1, a1 + b2]).TensorProduct[b1, b2] // 

However, once I enter more complex expressions, it seems to go wrong and give me

  TensorProduct[a2, n1].TensorProduct[b2, n1].TensorProduct[a2, 
     b1].TensorProduct[a2, b2] // TensorExpand]]

I am still getting invalid expression of the form

a2.b1 (n1.b2 a2\[TensorProduct]n1).a2\[TensorProduct]b2

and they don't really make sense.


2 Answers 2


You cannot use @ as it is used for assignments etc. So just set up your own tensor product tp and inner product ip and define

ip[tp[a,b],tp[c,d]]:=ip[b,c] tp[a,d]



It is unclear as to what you mean by "respect scalar multiplication". If you only want this with respect to explicit numbers you can set, for example,

tp[(s : _Integer | _Rational | _Real|_Complex) a_,b_]:=s tp[a,b]

If, however, you want to use variables that are scalar, you first need to tell Mathematica which symbols represent vectors or scalars. For example, you can define


which is supposed to mean that the symbols a, b and c are vectors and everything else a scalar. Now, you can set

tp[s_*a_?myVectorsQ,b_]:=s tp[a,b]
tp[a_,s_*b_?myVectorsQ]:=s tp[a,b]

Note that this approach is independent of the dimension and that one can define other functions on those tensors (like derivatives or integrals).

  • $\begingroup$ Hi - thank you for the suggestion but unfortunately doesn't solve my problem. I would like to do some algebra on it. Particularly, simplify and expand etc. So maybe a combination with built-in methods is better?! I will update the question. $\endgroup$ Commented Feb 8, 2016 at 19:53
  • $\begingroup$ The simplification rules above would either be contained in the definitions I suggested or could be added using UpValues. Also, keep in mind that "simplification" is subjective and context dependent. $\endgroup$
    – Berg
    Commented Feb 8, 2016 at 20:08
  • $\begingroup$ Hm - the simplification I need is multiplying out products and taking the trace. (Expand) $\endgroup$ Commented Feb 8, 2016 at 20:31
  • $\begingroup$ If the arguments of tp are not in the right form you could use the rule exp/.tp[a_,b_]:>tp[Expand@a,Expand@b] to expand the arguments. Then the linearity in the definition of tp would apply. $\endgroup$
    – Berg
    Commented Feb 8, 2016 at 20:59
  • $\begingroup$ Hello Berg, I am trying your suggestions but it is a little difficult to enter everything. Is there a way to use the built-in multiplication instead of ip. So that for instance tp[a,b].tp[c,d]=b.c tp[a,d] $\endgroup$ Commented Feb 8, 2016 at 22:36

I managed to to it relatively simply. This introduces the objects:

$Assumptions = {Element[x, Matrices[{3, 1}, Reals]], Element[a, Matrices[{3, 1}, Reals]], Element[n, Matrices[{3, 1}, Reals]], Element[b, Matrices[{3, 1}, Reals]], Element[m, Matrices[{3, 1}, Reals]], Element[id, Matrices[{3, 3}, Reals, Symmetric]]};

and this set of simplification rules does what I wanted it to do.

Simp[expr_] := expr //. {x__.Dot[Transpose[tensor1__], tensor2__].y__ :> tensor1.tensor2*x.y} SimId[expr_] := expr //. {Dot[id, tensor_] :> tensor, Dot[tensor_, id] :> tensor, MatrixPower[id, 2] :> id, MatrixPower[id, 3] :> id, MatrixPower[id, 4] :> id, MatrixPower[id, 5] :> id, MatrixPower[id, 6] :> id, MatrixPower[id, 7] :> id, MatrixPower[id, 8] :> id, MatrixPower[id, 9] :> id} SimTranspose[expr_] := expr //. {Transpose[n_, {2, 1}] :> Transpose[n], Dot[x_, {2, 1}] :> x} SimTr[expr_] := expr //. {Dot[tensor2__, Transpose[tensor1__]] :> tensor2.tensor1} Sim[expr_] := expr //. {Dot[a, n] :> 0, Dot[n, a] :> 0, Dot[b, m] :> 0, Dot[m, b] :> 0, Dot[n, n] :> 1, Dot[m, m] :> 1} SimAll[expr_] := Sim[SimTr[SimTranspose[Simp[SimId[TensorExpand[expr]]]]]]

Then one can just enter any expression of the described from and get the Trace. The function Sim is optional and can contain certain things things that cancel. E.g. $a \perp n \Rightarrow a.n=0$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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