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

Question

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.

I had already written MATLAB code which works around these problems:

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!

Answer 1

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 *)

Answer 2

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.