# Generalized n-fold inner product for tensors, iteration and indexing or built in functions?

I'm very new to Mathematica, and struggling to convert one of the codes I'd written in MATLAB. I am trying to program a function for a generalised inner product between two tensors, that is for example:

In this particular case intuitively I would have three "for" loops on a,b,c, and another three "for" loops on variables l,r and t to place them accordingly in the resulting tensor. Obviously I cannot do it like that, I do not know beforehand the order of input tensors, nor the "fold" of my inner product.

function t_out = inner_product(A, B, nfold)

%Outputs array 1xN indicating the number of elements in each dimension.
%The length of this vector corresponds to the number of dimensions N of
%input arrays
dimA = size(A); dimB = size(B);

%Collect the free indicies of inner product.  length(out_dim) determines the
%dimension, its contents represent number of elements in each dimension.
out_dim = [dimA(nfold+1:end), dimB(nfold+1:end)];

%Dimensions of free indicies associated to input matricies A and B
eblock_A = length(dimA(nfold+1:end)); eblock_B = length(dimB(nfold+1:end));

%Memory pre-allocation
t_out = zeros(out_dim);

if isempty(out_dim)
t_out = sum(sum(A.*B));
else
for i = 1:numel(t_out)

%Identifies which element is being worked on in an output matrix
out_access = cell(1,numel(out_dim));
[out_access{:}] = ind2sub(size(t_out),i);

%generates accessors for matricies A and B
sum_accg(1:nfold) = {':'};
sum_accA = [sum_accg, out_access{1:eblock_A}];
sum_accB = [sum_accg, out_access{eblock_A+1:end}];

%Element by element multiplications on "n-fold" dimensions and
%their sum
tmp = A(sum_accA{:}).*B(sum_accB{:});
t_out(out_access{:}) = sum(tmp(:));

end
end
end


For those who do not want to read into it, I work around the summation on a,b and c by choosing a portion of tensors A and B that fall under the summation, multiply them element-by-element, and add all the elements of the resulting tensor. I work around the "for" loops needed to place the elements in the resulting tensor by converting between MATLAB's linear and subscript indexing, in one "for" loop.

The issue I have is that Mathematica does not support linear indexing and conversion between the two (none that I know of). Flattening the table is obviously not an option. I had thought I can introduce a variable number of loops by doing something to the effect of:

InnerProduct[tensor1_, tensor2_, nfold_] := (
dim1free = Dimensions[tensor1][[nfold + 1 ;;]];
dim2free = Dimensions[tensor2][[nfold + 1 ;;]];
t2offset = Dimensions[dim1free][[1]];
iterList = {Table[{x[a], dim1free[[a]]}, {a, t2offset}],
Table[{x[a + t2offset], dim2free[[a]]}, {a,
Dimensions[dim2free][[1]]}]};

Table[Total[Times[Extract[tensor1, pos1], Extract[tensor2, pos2]]], iterList]
)


but that didn't end up working out for me either.

Maybe you know a better, more efficient way of solving this problem using internal mathematica functions?

Many thanks!

Does this function do what you are looking for?

InnerProduct[T1_, T2_, k_] := Dot[
Transpose[Flatten[T1, k - 1],RotateRight[Range[TensorRank[T1] + 1 - k]]],
Flatten[T2, k - 1]
]


It flattens the first k slots together, rotates the flattened slot of the first tensor to the right and uses Dot to compute the summation efficiently.

Here is a test example for demonstration, where the formula

$R_{lrt} = (A \cdot^3 B)_{lrt}=\sum_{a}\sum_{b}\sum_{c}A_{abcl}B_{abcrt}$

is computed explicitely, but less efficiently:

SeedRandom[1]
A = RandomReal[{-1, 1}, {3, 5, 7, 11}];
B = RandomReal[{-1, 1}, {3, 5, 7, 13, 17}];

n1 = Dimensions[A][[1]];
n2 = Dimensions[A][[2]];
n3 = Dimensions[A][[3]];
nfold = 3;
R = ConstantArray[0.,Join[Dimensions[A][[nfold + 1 ;;]], Dimensions[B][[nfold + 1 ;;]]]];

Do[
Do[
Do[
R[[l, r, t]] = Sum[
A[[a, b, c, l]] B[[a, b, c, r, t]],
{a, 1, n1},
{b, 1, n2},
{c, 1, n3}
],
{l, 1, Dimensions[A][[4]]}],
{r, 1, Dimensions[B][[4]]}],
{t, 1, Dimensions[B][[5]]}];


Compared to InnerProduct[A,B,nfold], we have only negligible error:

Max[Abs[R - InnerProduct[A, B, nfold]]]

(* 1.06581*10^-14 *)

• Thank you for your answer, but unfortunately not. While this worked on a two matricies for k=1 and k=2, it gives me an incorrect answer for tensors of higher orders. For example: AAA = {{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}}; InnerProduct[AAA,AAA,2]. Due to symmetry of the tensor in the 3rd dimension, my output matrix should have all the same values (the double summation comes out the same). Matlab output: {{285,285,285},{285,285,285},{285,285,285}} Mathematica output: {{198, 234, 270}, {234, 279, 324}, {270, 324, 378}} – Slava K. Jul 28 '17 at 9:50
• @SlavaK." Are you sure that the texed formula in your question displays what you look for? Because that one reduces to KroneckerProduct if A and B are vectors and nfold equals 0. And that is in general not a constant matrix... – Henrik Schumacher Jul 28 '17 at 10:03
• I'm not sure what you mean. The result ususally is not a constant matrix, it just should happen to be so in the case of InnerProduct[AAA,AAA,2]. InnerProduct is a generalisation, so you correct in saying that InnerProduct[T1,T2,0] = TensorProduct[T1,T2] (or KroneckerProduct for vectors), InnerProduct[T1,T2,1] = Dot[Transpose[T1],T2]. I just need to generalise for any higher nfold (so more summations and lower order output tensor). Unfortunately Flatten[] cannot be used on T1 and T2 because it reduces the dimensions but does not insure correct subsequent multiplication – Slava K. Jul 28 '17 at 11:01
• You are correct, sorry. After manually tesing with for loops in MATLAB as well and still getting different answers turns out I was inputing different matricies into the programs due to my lacking knowledge in Mathematica syntax. Thank you! – Slava K. Jul 28 '17 at 12:45
• You're welcome! – Henrik Schumacher Jul 28 '17 at 13:57

Conceptually, I think the following is the clearest:

ip[a_, b_, k_] := With[{rank = TensorRank[a]},
TensorContract[
TensorProduct[a, b],
Table[{i, rank+i}, {i, k}]
]
]


For example:

A = RandomReal[1, {2, 3, 4, 5}];
B = RandomReal[1, {2, 3, 4, 6, 7}];

ip[A, B, 3] //Dimensions


{5, 6, 7}

For large tensors, this approach is inefficient because it constructs a really large intermediate tensor that it contracts over. In this case, it makes sense to use Inactive/Activate in the following way:

ip2[a_, b_, k_] := With[{rank = TensorRank[a]},
Activate @ TensorContract[
Inactive[TensorProduct][a, b],
Table[{i, rank+i}, {i, k}]
]
]


A speed comparison:

r1 = ip[A, B, 3]; //RepeatedTiming
r2 = ip2[A, B, 3]; //RepeatedTiming

MinMax[r1-r2]


{0.00032, Null}

{0.000091, Null}

{-3.55271*10^-15, 2.66454*10^-15}

Note that the approach by @Henrik is even faster.

• Hi, thank you for your answer. As I had said to Henrik, the dimensions of the output tensors are correct, but the numerical values are wrong – Slava K. Jul 28 '17 at 9:28