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));
    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(:));


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, 

  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!


2 Answers 2


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:

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

    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 *)
  • $\begingroup$ 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}} $\endgroup$ Commented Jul 28, 2017 at 9:50
  • $\begingroup$ @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... $\endgroup$ Commented Jul 28, 2017 at 10:03
  • $\begingroup$ 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 $\endgroup$ Commented Jul 28, 2017 at 11:01
  • $\begingroup$ 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! $\endgroup$ Commented Jul 28, 2017 at 12:45
  • $\begingroup$ You're welcome! $\endgroup$ Commented Jul 28, 2017 at 13:57

Conceptually, I think the following is the clearest:

ip[a_, b_, k_] := With[{rank = TensorRank[a]},
        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


{0.00032, Null}

{0.000091, Null}

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

Note that the approach by @Henrik is even faster.

  • $\begingroup$ 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 $\endgroup$ Commented Jul 28, 2017 at 9:28

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.