# Trying to find the double dot product of a Rank 4 and Rank 2 tensor

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];

• 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... – Henrik Schumacher Nov 21 '19 at 12:43

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

1.63814

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

0.06622

0.000292

All of them lead to eseentially the same result:

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

8.88178*10^-15

0.