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