# Update: Calculate Tensor products and traces with Mathematica

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  SimId[SimTP[ TensorProduct[a2 + b1, n1].(id + TensorProduct[a1, a1 + b2]).TensorProduct[b1, b2] // TensorExpand]]  However, once I enter more complex expressions, it seems to go wrong and give me  SimId[SimTP[ 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 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 tp[id,v_]:=v tp[v_,id]:=v Tra[tp[a,b]]:=ip[a,b] Trans[tp[a,b]]:=tp[b,a] ip[tp[a,b],tp[c,d]]:=ip[b,c] tp[a,d]  and tp[a_+b_,c_]:=tp[a,c]+tp[b,c]  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 myVectorsQ[exp_]:=MatchQ[exp,a|b|c]  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). • 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. Commented Feb 8, 2016 at 19:53 • 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. – Berg Commented Feb 8, 2016 at 20:08 • Hm - the simplification I need is multiplying out products and taking the trace. (Expand) Commented Feb 8, 2016 at 20:31 • 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. – Berg Commented Feb 8, 2016 at 20:59 • 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] 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$.