File Exchange

## Multiple matrix multiplications, with array expansion enabled

version 1.37.0.0 (505 KB) by Paolo de Leva

### Paolo de Leva (view profile)

Multiplying matrices, vectors, or scalars contained in two N-D arrays, with array expansion enabled.

Updated 26 Jul 2010

MULTIPROD is a powerful, quick and memory efficient generalization for N-D arrays of the MATLAB matrix multiplication operator (*). While the latter works only with 2-D arrays, MULTIPROD works also with multidimensional arrays.
MULTIPROD performs multiple multiplications between matrices, vectors, or scalars contained in two multidimensional arrays, with automatic virtual array expansion (AX) enabled. AX allows you, for instance, to multiply a single matrix A by an array of matrices B, by virtually replicating A to obtain an array compatible with B.

Multidimensional arrays may contain matrices or vectors or even scalars along one or two of their dimensions. For instance, a 4×5×3 array A contains three 4×5 matrices along its first and second dimension (fig. 1). Thus, array A can be described as a block array the elements of which are matrices, and its size can be denoted by (4×5)×3.

MULTIPROD can be also described as a generalization of the built-in function TIMES. While TIMES operates element-by-element multiplications (e.g. A .* B), MULTIPROD operates block-by-block matrix multiplications.

EXAMPLES

Let's say that

A is (2×5)×6, and
B is (5×3)×6.

With MULTIPROD the six matrices in A can be multiplied by those in B in a single intuitively appealing step:

C = MULTIPROD(A, B).

where C is (2×3)×6.

By automatically applying AX, MULTIPROD can multiply a single matrix by all the blocks of a block array. So, if

A is 2×5 (single matrix), and
B is (5×3)×1000×10,

then C = MULTIPROD(A, B) yields a (2×3)×1000×10 array. A is virtually expanded to a (2×5)×1000×10 size, then multi-multiplied by B. This is done without using loops, and without actually replicating the matrix (see Appendix A). We refer to this particular application of AX as virtual matrix expansion. In a system running MATLAB R2008a, MULTIPROD performs it about 380 times faster than the following equivalent loop (see Appendix B):

for i = 1:1000
for j = 1:10
C(:,:,i,j) = A * B(:,:,i,j);
end
end

AX generalizes matrix expansion to multidimensional arrays of any size. For instance, if

A is (2×5)×10, and
B is (5×3)×1×6,

then C = MULTIPROD(A, B) multiplies each of the 10 matrices in A by each of the 6 matrices in B, obtaining 60 matrices stored in a (2×3)×10×6 array C. It does that by virtually expanding A to (2×5)×10×6, and B to (5×3)×10×6. A detailed definition of AX is provided in the manual.

APPLICATIONS

MULTIPROD has a broad field of potential applications. By calling MULTIPROD, multiple geometrical transformations such as rotations or roto-translations can be performed on large arrays of vectors in a single step and with no loops. Multiple operations such as normalizing an array of vectors, or finding their projection along the axes indicated by another array of vectors can be performed easily, with no loops and with two or three rows of code.

Sample functions performing some of these tasks by calling MULTIPROD are included in the separate toolbox "Vector algebra for multidimensional arrays of vectors" (MATLAB Central, file #8782).

OPTIMIZATION AND TESTING

Since I wanted to be of service to as many people as possible, MULTIPROD was designed, debugged, and optimized for speed and memory efficiency with extreme care. Precious advices by Jinhui Bai (Georgetown University) helped me to make it even faster, more efficient and more versatile. Suggestions to improve it further will be welcome. The code ("testMULTIPROD.m") I used to systematically test the function output is included in this package.

THE ARRAYLAB TOOLBOX

In sum, MULTIPROD is a generalization for N-D arrays of the matrix multiplication function MTIMES, with AX enabled.

Vector inner, outer, and cross products generalized for N-D arrays and with AX enabled are performed by DOT2, OUTER, and CROSS2 (MATLAB Central, file #8782, http://www.mathworks.com/matlabcentral/fileexchange/8782).

Element-by-element multiplications (see TIMES) and other elementwise binary operations (such as PLUS and EQ) with AX enabled are performed by BAXFUN (MATLAB Central, file #23084, http://www.mathworks.com/matlabcentral/fileexchange/23084).

Together, these functions make up the “ARRAYLAB toolbox”. I hope that The MathWorks will include it in the next version of MATLAB.

MULTITRANSP

This package includes the function MULTITRANSP, performing multiple matrix transpositions. B = MULTITRANSP(A, DIM) transposes all the matrices contained along dimensions DIM and DIM+1 of A.

### Cite As

Paolo de Leva (2019). Multiple matrix multiplications, with array expansion enabled (https://www.mathworks.com/matlabcentral/fileexchange/8773-multiple-matrix-multiplications-with-array-expansion-enabled), MATLAB Central File Exchange. Retrieved .

Louis Vallance

### Louis Vallance (view profile)

Excellent code which does exactly what I was hoping for. Bravo!

Sebastian Kraemer

### Sebastian Kraemer (view profile)

Dear Paolo,

I believe one example you give is wrong:

% A single matrix A pre-multiplies each matrix in B
% If A is ........................... a (1�3) single matrix,
% and B is ........................... a 10�(3�4) 3-D array,
% C = MULTIPROD(A, B, [1 2], [3 4]) is a 10�(1�4) 3-D array.

It should be (A,B,[1 2], [2 3]), correct?

Xiaodong

### Xiaodong (view profile)

On Matlab 2016a, when I use multiprod, I got a warning: "Warning: NARGCHK will be removed in a future release. Use NARGINCHK or NARGOUTCHK instead. " I guess this can be easily fixed by replacing a function the program calls. Would you like to make an update?

Timothy Sinclair

### Timothy Sinclair (view profile)

I tried to use bsxfun() to multiply multiple sheets of matrices with multiple sheets of vectors but I didn't realy bsxfun() couldn't do that ... Stackexchange suggestions didn't work either and then I finally found this! Thanks!

SK

### SK (view profile)

In general, from the point of view of speed and memory efficiency, it is not a good ides to convert matrix multiplications into element-wise multiplications (even via bsxfun).

I see from your examples that it looks like your intention is to use this for large numbers of multiplications of small matrices - in which case the comments below may not be so relevant.

SK

### SK (view profile)

In general, from the point of view of speed and memory efficiency, it is not a good ides to convert matrix multiplications into element-wise multiplications (even via bsxfun).

SK

### SK (view profile)

Nice Idea, but I would not call it memory efficient. For example trying to multiply the slices of two 1000x1000x10 double arrays results in multiprod calling bsxfun with arguments of size 1000x1000x1x10 and 1x1000x1000x10 which (at least naively) requires 160 Gb of memory. This results in Out of Memory error from Matlab.

Also speed does not appear to be great. Try:
A = rand(250, 250);
B = rand(250, 250);
E = repmat(A, [1,1,10]);
F = repmat(B, [1,1,10]);

tic;
C = zeros(250,250,10);
for i = 1 : 10
C(:,:,i) = E(:,:,i)*F(:,:,i);
end
toc;

Elapsed time is 0.044976 seconds.

But:
tic;
X = multiprod(E, F, [1,2], [1,2]);
toc;

Elapsed time is 0.744825 seconds.
This is about 16 time slower.

Nathan Thern

### Nathan Thern (view profile)

@cszczyglowski: if M2 is square, you can perform M1/M2 with multiprod and multinv (http://www.mathworks.com/matlabcentral/fileexchange/31222)

Ilan

### Ilan (view profile)

Works great with cuda

cszczyglowski

### cszczyglowski (view profile)

Fantastic bit of code, thank you very much! Is there an equivalent function for performing multiple right-divide operations in a single step instead of using a loop?

Yuan Zhou

### Yuan Zhou (view profile)

This is a great piece of work! It saves a lot of time.

Em Rickinson

### Em Rickinson (view profile)

Great function! I was wondering if there is a similar function for matrix power. Thanks.

Mikhail Katliar

### Mikhail Katliar (view profile)

Very helpful, exactly what I needed.

Manuel A. Diaz

Denis Anikiev

Romain W

### Romain W (view profile)

Hi Paolo,
I was wondering if it is possible without using a for-loop to multiply each row of a matrix by another matrix using multiprod in a vectorised way.

In my case: A [55000,3] and B [55000,3] and I want each row A [1,3] * B'[3,55000]? And I want something like:
t=1:55000
multiprod(A(t,:),B)

Very good job for this function! Cheers

Rajab Legnain

### Rajab Legnain (view profile)

Hi Sir

It is really a great function.
Thanks for sharing.
Regards

Oliver Woodford

### Oliver Woodford (view profile)

Very useful in general, but there is an overhead which becomes significant when applying this function to smaller arrays.

YIXIAO

Mark Flynn

### Mark Flynn (view profile)

Hi Paolo,

Thank you so much for your help! I've been trying to figure out how to vectorize this bit of my code for a long time. Dot2 really did the trick. This was by far the most valuable tool I've gotten from the Matlab Central.

Paolo de Leva

### Paolo de Leva (view profile)

MULTIPROD does multiple matrix products. For inner products, you should use function DOT2, which is included in my vector algebra toolbox (File ID: #8782).

Notice that, if you use complex numbers, there is a difference between an inner product (MATLAB builtin function DOT) between two row vectors A and B and the matrix product between A and B' (the transpose of B).
Namely,

DOT(A, B) = CONJ(A) * B'.

When you use complex numbers, this makes a difference. When you use real numbers, CONJ(A) is the same as A and as a consequence

DOT(A, B) = A * B'.

Since my function DOT2 has array expansion enabled, it can perform the operation you need with a single command, provided that you define the two block arrays as a 24x1x10 array A and a 1x64x10 array B. Both arrays contain 10-element vectors along dimension 3.

For instance:
A = rand(24,1,10)
B = rand(1,64,10)
PRODUCTS = DOT2(A, B, 3)

PRODUCTS will be a 24 x 64 matrix.

See MULTIPROD manual for details about the Array expansion (AX) feature implemented in MULTIPROD's engine.

Enjoy my Arraylab toolboxes. Thank you for your rating.

Mark Flynn

### Mark Flynn (view profile)

This looks like something I've been searching for for a while. I have 2 sets of vectors where I want to get the inner product of each vector in the first set with each vector in the second. E.g. 24 10x1 vectors with 64 10x1 vectors. Is this something multiprod can do and could you give me some suggestions for how I could do it?

uniquebird

### uniquebird (view profile)

I want to multiply an matrix of nxnxn by a nxnxn matrix with the same place.

Paolo de Leva

### Paolo de Leva (view profile)

Hi, A.

You don't need to reshape. MULTIPROD's array expansion (AX) engine does everything for you automatically, provided you use the correct syntax.

Since you need to multiply vectors (1-D blocks) by a matrix M (2-D block), you need to use syntax 3:

vectors = randn(15,15,15,5);
matrix = randn(5,5);
C = multiprod(vectors, matrix, 4, [1 2]);

It is assumed, since you want vector-by-matrix products, that you use the row-vector convention (each vector is considered to be a 1x5 matrix). If you use the column vector convention (each vector is considered to be a 5x1 matrix), then you need matrix-by-vector products:

C = multiprod(matrix, vectors, [1 2], 4);

In both cases, the size of C is [15,15,15,5], as expected.

I suggest you to read the PDF manual. If some part of the manual is not clear for you, please let me know by e-mail.

If you like multiprod, please rate it. I greatly enjoy receiving five-star ratings :-)

Michael Völker

### Michael Völker (view profile)

Hi A, then I think that
C = reshape(multiprod( A, B, 4, 1 ), size(A));
will work.

A

### A (view profile)

Michael,

that worked; except sometimes I might have singleton dimensions in my actual data that I don't want touched by squeeze function; e.g. 15x1x15x5 times 5x5 = 15x1x15x5

Michael Völker

### Michael Völker (view profile)

Sure you can:
A = randn(15,15,15,5);
B = randn(5,5);
C = squeeze(multiprod( A, B, 4, 1 ));

I think this is what you expect.

A

### A (view profile)

I would like to multiply an array of 5-element vectors by a 5x5 matrix. In other words, multiply 15x15x15x5 by 5x5 to get 15x15x15x5

can I use multiprod?

Dinakar

### Dinakar (view profile)

A fantastic piece of code. Was able to get rid of the loop and significantly speed up my code.

Mathworks should consider including this in their standard distribution.

Paolo de Leva

### Paolo de Leva (view profile)

Shatrughan, I do not know outerm and PLS, and I cannot completely understand your request, but I can try and explain what you can do with my functions. Allow me to start with the simplest case:

A = rand(7,1);
B = rand(9,1);
X = outer(A, B);

Here, X is a matrix, and its size is 7x9. It would make sense to write a function to generalize this to 3-D:

A = rand(7,1);
B = rand(9,1);
C = rand(10,1);
X = simple_outer3D(A, B, C);

where X is a 3-D array and its size is 7x9x10, which is the size of your expected result. This can be done using this code (where BAXFUN is a function that I posted on MATLAB Central File Exchange; File ID: #23084):

function X = simple_outer3D(A, B, C)
% A, B, and C must be column vectors
X = baxfun(@times, A, B, 0, 1); % Size: 7x9
X = baxfun(@times, X, C, 0, 2); % Size: 7x9x10

Warning: this code is completely valid as a generalization of an outer product only if A, B, and C do not contain complex numbers. See the help of function OUTER for details.

Notice that I got the result you expected, but the input matrices have not the same size as in your example. Conversely, input matrices with the same size as in your example would give a different result. Simple case:

A = rand(7,2);
B = rand(9,2);
X = outer(A, B);

Here, A and B are not two vectors, but two “arrays of vectors” (each containing two vectors). Hence, X is not a single outer product, but a 7x9x2 array containing two outer products. The first outer product is X(:,:,1) and the second is X(:,:,2). If we generalize to 3-D, we have:

A = rand(7,2);
B = rand(9,2);
C = rand(10,2);
X = outer3D(A, B, C);

In this case, X is supposed to contain two 3-D arrays. Hence, its size should be 7x9x10x2 (rather than 7x9x10). This can be done by defining outer3D as follows:

function X = outer3D(A, B, C)
A = reshape(A, [size(A,1),1,1,size(A,2)]); % Size: 7x1x1x2
B = reshape(B, [1,size(B,1),1,size(B,2)]); % Size: 1x9x1x2
C = reshape(C, [1,1,size(C)]); % Size: 1x1x10x2
X = baxfun(@times, A, B); % Size: 7x9x1x2
X = baxfun(@times, X, C); % Size: 7x9x10x2

Warning: this code is completely valid as a generalization of an outer product only if A, B, and C do not contain complex numbers. See the help of function OUTER for details.
Note: In this code, BAXFUN (binary array expansion) is equivalent to BSXFUN (binary singleton expansion).

Michael Völker

### Michael Völker (view profile)

Shatrughan, what is the code you would write without having multiprod()?

As far as I can see, multiprod() accepts 2 input arrays while you want to do something with 3 arrays.

multiprod() does nothing which you cannot do with regular code, it just does it more elegantly and faster.

So, please write down an example code which you expect to get "beautified" or accelerated with the help of multiprod().

Shatrughan

### Shatrughan (view profile)

Hey Paolo,

Very interesting toolbox, appreciate your hard work towards its development...however, I could not implement it on my need may be I am not applying it correctly, your help is needed..

I have three matices (say, a=7x2; b=9x2; c=10x2)

Output I want, X=7x9x10; (similar as you can get in outerm using PLS but I do not have access to the same);

I tried using using outer script in your Vector toolbox also, but to no avail..Any help will be much appreciated
Thanks.!

Chris

### Chris (view profile)

This is great -- it has been a tremendous improvement for my project.

Nicolas

Xiaodong

### Xiaodong (view profile)

Hi Paolo de Leva

Great job! I like your toolbox very much. As inspired by your code, I made a MULTINV function to do the inversion for N-D array in a similar way: http://www.mathworks.com/matlabcentral/fileexchange/31222-inversion-every-2d-slice-for-arbitrary-multi-dimension-array

Now we have both multiplication and inversion functions for array operations!

If you are interested, welcome to send comments on my code! Hopefully, we can jointly make a fully functional toolbox for array operations.

With best,
Xiaodong

Daniel Armyr

### Daniel Armyr (view profile)

Ambitious documentation, fast implementation, and best of all, it worked out-of-the box. The lastes resident in my private toolbox.

Paul

### Paul (view profile)

Incredible! I have to do a series of structured matrix-vector multiplications where the matrices have identical structure. There was no way to do this with built-in matlab functionality, outside of using a FOR loop. A MULTIPROD-based implementation turned out to be 100-1000 times faster on problems of practical size!

Paul

Kirankumar

Michael Völker

### Michael Völker (view profile)

Excellent. Not only does it what it does but you can learn a lot by playing with the code.
By inlining the essential part I could dramatically increase speed, since less copying of data is required.

Soravit

Tanyer Alan

Aleksandra

Deepak

### Deepak (view profile)

Nice program to change all the lines of code into one command. really useful in simulations etc.

Denis Kozlov

Paolo de Leva

### Paolo de Leva (view profile)

Chris, thak you for your rating. What you describe is indeed a useful operation which is described in geometry as a transformation matrix composition (a special kind of function composition), but it is different from what MULTIPROD does. MULTIPROD is designed to operate on 2 block arrays of matrices or vectors or scalars, not on a single block array. You can compare it with DOT and CROSS, which also operate on 2 block arrays. In other words, MULTIPROD, DOT and CROSS can be described as binary operations, while you need a unary operation (operating on a single operand). Possibly, you or some other author will be willing to design and publish an optimized function to do what you need. Such a function would be consistent with what I called the "ARRAYLAB" phylosophy (an evolution of the original MATLAB approach, see MULTIPROD manual).
Sorry for rating myself. I did that by mistake and did not find a way to undo it.

Chris Sarantos

### Chris Sarantos (view profile)

Question: Can multiprod be used to multiply a single array of matrices together along one dimension? For example if you have a 2x2xN array, and want to multiply the 2x2 matrices together along the third dimension to get a 2x2x1 output?

Erdal Bizkevelci

Noah Snavely

### Noah Snavely (view profile)

This is an amazingly useful tool. Anyone who uses matrices will---at some point in their life---be stymied by a problem for which this software is the perfect solution. It just works.

Friedemann Groh

Cool. Thank You!

Paul Macey

Brilliant - this does exactly what I want, performing vector re-orientation on a 3D array (i.e., 3D x 3 x 3 rotation matrix multiplied by 3D x 3 vector matrix); very fast.

Lucas Finn

Well documented, well optimized, and works.

Giuseppe Vannozzi

Awesome! Thank you for sharing this useful package

Pietro Picerno

just fantastic, many compliments, it's a great and well-done job!

Arun Mano

awesome stuff! just got what i was looking for!
Pretty efficient too!

Matthew G

This is an extremely useful library, I was trying to figure out how to do this on my own, but I'm lucky this page came up in the google search.

Ryan Eustice

Just what I was looking for. This capability really should be a builtin for Matlab.

paul naude

Now, where do I find a file that can subtract my DC offset (vector) from my data matrix? I.e B = A - V where A may have many rows and the columns in A and number of elements in V agree. Obviously I don't want to find the size of A first and then run a FOR loop!

paul naude

I needed to multiply a matrix M (with variable number of rows) with a calibration vector V and by using multiprod(M,V,), it callibarted my file in this ONE step!

Magnus Karlsson

Super!
Have I longed for this functionality in Matlab! Farewell for-looops!

This should be an obvious part of Matlab, and it is amazing it has not been done until now.
It extends the functionality of Matlab in amazing ways, as simulations requiring for-loops now can take advantage of Matlab's inherent fast matrix multiplication algorithms.

Mr deLeva deserves Mathwork's medal in gold for this.

Bulent Ozer

Good job. Thanks to your program got rid of the loop. My code is much faster now. Thanks

Randy Wells

Good job. Very useful when manipulating time dependant transformation matrices (vector of 3x3 matrices). No longer need a for loop.