In Matlab, there is no builtin function for the multiplication of multidimensional arrays. There is a function called TPROD on the file exchange, which does the job, but it depends on mex files, and therefore on a compiler, and it also does not support sparse matrices.

I wrote a pure matlab script, called tensor_mult, which does the same thing, sometimes faster, sometimes slower. I havenâ€™t figured out why, yet. It simply uses `reshape`, `permute` and basic matrix-matrix multiplication to do the job.

Below you can see the whole script. An example, using the (semi-)Einstein convention, is `C_{bde} = A_{abc} B_{dcea}`, i.e. summation is to be done over the indices `a` and `c`. The function call is `tensor_mult(A, B, [1 3], [4 2])`, which means: compute the tensor product between `A` and `B`, while summing over the 1st index of A and the 4th index of B, and over the 3rd index of A and the 2nd index of B.

``````function C = tensor_mult(A, B, sum_idx_A, sum_idx_B)

sum_idx_A = reshape(sum_idx_A, 1, numel(sum_idx_A));
sum_idx_B = reshape(sum_idx_B, 1, numel(sum_idx_B));

size_A = size(A);
size_B = size(B);

num_size_A = numel(size_A);
num_size_B = numel(size_B);

perm_A = 1:num_size_A;
perm_B = 1:num_size_B;

num_idx = size(sum_idx_A, 2);

assert(isequal(size(sum_idx_A) , size(sum_idx_B) ));
assert(isequal(size_A(sum_idx_A), size_B(sum_idx_B)));

sum_dim = prod(size_A(sum_idx_A));

for i = 1:num_idx
perm_A = perm_A(perm_A~=sum_idx_A(i));
perm_B = perm_B(perm_B~=sum_idx_B(i));
end

perm_A = [perm_A, sum_idx_A];
perm_B = [sum_idx_B, perm_B];

size_A(sum_idx_A) = [];
size_B(sum_idx_B) = [];

if isempty(size_A)
size_A = 1;
end
if isempty(size_B)
size_B = 1;
end

C = squeeze(reshape(...
reshape(permute(A, perm_A), [prod(size_A), sum_dim]) * ...
reshape(permute(B, perm_B), [sum_dim, prod(size_B)]), ...
[size_A,size_B]));

end
``````