Matrix times array

Fast multiplication of multi-dimensional array x by matrix A along dim. n (or: mode-n product of tensor x and matrix A)

out = MatrixTimesArray(A,x,i);

Multiplication of multi-dimensional array x by matrix A along dimension i. This is equivalent to mode-i product of tensor x by matrix A in tensor language, see for example De Lathauwer et al. (SIAM 2000). No further mex-files needed.

Elements of output are determined by the formula:

The function reshapes x into a matrix (performing a permute, if necessary) and then either pre- or post-multiplies it by A, whichever minimizes permute operations). If x is a matrix, the following equivalence holds to conventional matrix algebra:

A*x = MatrixTimesArray(A,x,1);   % dim. 1: matrix pre-multiplication
x*B'= MatrixTimesArray(B,x,2);   % dim. 2: matrix post-multiplication

As for speed: For large arrays, note that the permute operations cost time time, not the reshape. Thus, put most computationally intense dimensions on first (preferred) and last dimension.

Inputs:

A: matrix, full or sparse.

x: k-dimensional array of size

i: Scalar, optional. Along which dimension to multiply. If omitted, the function multiplies along the first compatible dimension

Output:

: k-dimensional array of size . Note that the size of dimension i has changed from to m with respect to the input x.

Example:

A = rand(60,30);                  % Matrix A: 60-by-30.
x = rand(50,30,20);               % Array/tensor x: (50,30,20)
y = MatrixTimesArray(A,x,2);      % Multiply along dimension 2:
size(y)                           % y is (50,60,20) array.
yy= MatrixTimesArray(A,x);        % Omitting last input (dim.) gives same result.
norm(yy(:)-y(:));  

Version: 20-Dec-2021

Author: Matthias Kredler (Universidad Carlos III de Madrid), [e-mail]

Acknowledgement: This code follows the same idea as my function 'kronm.m', which carries out a series of MatrixTimesArray / n-mode products on a tensor.