I understand this question has been asked before, or something very similar to what I'm trying to accomplish at the very least, however I didn't understand much of what was happening in the other solutions. I'm trying to find the double contraction of A, a rank 4 tensor, and B, a rank 2 tensor (a matrix), or A:B, and I'm not sure how to go about doing this on Mathematica. I know that the double contraction (or double dot product) is meant to yield a rank 2 tensor. The equation that I was given and need to solve for does not indicate any indices so I'm not sure how to go about it.

    A = KroneckerDelta[9,9]
    B = 0.5*IdentityMatrix[3]

I have tried:

    Dot[A,B] (*which makes no sense*)
    A*B (*makes absolutely no sense, but I am desperate at this point*)
    Tr[A,Transpose[B]] (*But I think this only works for rank 2 tensors*)

Edit: Here is a bit more context. I'm trying to find the double dot product of the projection tensor P and a matrix which are denoted by the following:

    I = Array[KroneckerDelta, {3,3}]; 
    J = Array[KroneckerDelta,{9,9}];
    A = J -((1/3)*TensorProduct[I,I]);
and B = 0.5*IdentityMatrix[3];

Yeah. Just need help going about this, any would be great. Have a great day

  • $\begingroup$ As said before in this former post, KroneckerProduct will never ever return a 4-tensor. And KroneckerProduct[9,9] is absolutely meaningless---which would have been told to you by Mathematica if you had bothered to test your own code beforehand... $\endgroup$ Nov 21, 2019 at 12:43

1 Answer 1


What $A:B$ is supposed to mean is hard to tell without context (I consider it bad notation anyways because I was told in German primary school that "$:$" means division). But my guess is the following:

n = 20;
A = RandomReal[{-1, 1}, {n, n, n, n}];
B = RandomReal[{-1, 1}, {n, n}];

a = TensorContract[TensorProduct[A, B], {{3, 5}, {4, 6}}]; // AbsoluteTiming // First 


More performant variants are however these:

b = Activate@TensorContract[Inactive[TensorProduct][A, B], {{3, 5}, {4, 6}}]; // AbsoluteTiming // First 
c = ArrayReshape[A, {n, n, n^2}].Flatten[B]; // AbsoluteTiming // First 



All of them lead to eseentially the same result:

Max[Abs[a - b]]
Max[Abs[b - c]]




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.