### Highlights from Multiple matrix multiplications, with array expansion enabled

Join the 15-year community celebration.

Play games and win prizes!

4.93617
4.9 | 49 ratings Rate this file 30 Downloads (last 30 days) File Size: 505 KB File ID: #8773 Version: 1.37

# Multiple matrix multiplications, with array expansion enabled

### Paolo de Leva (view profile)

20 Oct 2005 (Updated )

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

File Information
Description

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.

Acknowledgements
MATLAB release MATLAB 7 (R14)
Other requirements MULTIPROD calls the builtin function BSXFUN, available in MATLAB R2007a. For earlier MATLAB releases, use Schwarz's BSXFUN substitute (MATLAB Central, file #23005, http://www.mathworks.com/matlabcentral/fileexchange/23005)
07 Sep 2016 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.

Comment only
07 Sep 2016 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).

Comment only
07 Sep 2016 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.

27 May 2016 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)

27 Apr 2016 Ilan

### Ilan (view profile)

Works great with cuda

15 Feb 2016 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?

24 Oct 2015 Yuan Zhou

### Yuan Zhou (view profile)

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

27 Jul 2015 Em Rickinson

### Em Rickinson (view profile)

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

01 Aug 2014 Mikhail Katliar

### Mikhail Katliar (view profile)

Very helpful, exactly what I needed.

29 Jun 2014 Manuel A. Diaz

### Manuel A. Diaz (view profile)

27 May 2014 Denis Anikiev

### Denis Anikiev (view profile)

14 Mar 2013 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

18 Feb 2013 Rajab Legnain

### Rajab Legnain (view profile)

Hi Sir

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

14 Nov 2012 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.

31 Oct 2012 YIXIAO

### YIXIAO (view profile)

26 May 2012 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.

Comment only
25 May 2012 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.

Comment only
24 May 2012 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?

21 Apr 2012 uniquebird

### uniquebird (view profile)

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

11 Apr 2012 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 :-)

Comment only
10 Apr 2012 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.

Comment only
10 Apr 2012 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

Comment only
10 Apr 2012 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.

Comment only
10 Apr 2012 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?

Comment only
19 Feb 2012 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.

03 Feb 2012 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).

Comment only
03 Feb 2012 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().

Comment only
27 Jan 2012 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.!

22 Jun 2011 Chris

### Chris (view profile)

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

09 Jun 2011 Nicolas

### Nicolas (view profile)

28 Apr 2011 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

Comment only
27 Apr 2011 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.

14 Apr 2011 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!

14 Apr 2011 Paul

### Paul (view profile)

08 Apr 2011 Kirankumar

### Kirankumar (view profile)

15 Jan 2011 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.

10 Aug 2010 Soravit

### Soravit (view profile)

05 Aug 2010 Tanyer Alan

### Tanyer Alan (view profile)

04 Mar 2010 Aleksandra

### Aleksandra (view profile)

03 Nov 2009 Deepak

### Deepak (view profile)

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

12 Aug 2009 Denis Kozlov

### Denis Kozlov (view profile)

31 Jul 2009 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.

31 Jul 2009 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?

29 Jun 2009 Erdal Bizkevelci

### Erdal Bizkevelci (view profile)

21 May 2009 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.

10 Oct 2007 Friedemann Groh

Cool. Thank You!

28 Sep 2007 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.

04 Mar 2007 Lucas Finn

Well documented, well optimized, and works.

05 Dec 2006 Giuseppe Vannozzi

Awesome! Thank you for sharing this useful package

04 Dec 2006 Pietro Picerno

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

04 Aug 2006 Arun Mano

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

01 Aug 2006 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.

27 Apr 2006 Ryan Eustice

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

25 Apr 2006 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!

25 Apr 2006 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,[1]), it callibarted my file in this ONE step!

27 Jan 2006 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.

23 Nov 2005 Bulent Ozer

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

08 Nov 2005 Randy Wells

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

31 Oct 2005

I discovered an error in the discussion, which I hadn't found in my first editing. I corrected it and also decided to further refine the whole text to improve its readability.

22 Feb 2006

The satellite function REPMAT2 was substituted with MATEXP (matrix expansion). The manual was refined.

22 Feb 2006

The satellite function REPMAT2 was substituted with MATEXP (matrix expansion). The manual was refined.

01 Mar 2006

1. Solved character conversion problem in "File Description" text.
2. Included a "MULTIPROD toolbox manual" in Acrobat (PDF) format.

30 Mar 2006

The code was optimized according to suggestions automatically generated by function M-Lint.

30 Mar 2006

The code was optimized according to suggestions automatically generated by function M-Lint.

21 Nov 2006

Version 1.2.1 (2006b). Improved management of syntax errors by the user. Revised comments and manual.

03 Oct 2007

Fixed a bug: when input contained complex numbers, in some circumstances MULTIPROD returned an incorrect result. MULTIPROD is defined as a multiple matrix algebra multiplication over both the real and complex fields.

03 Oct 2007

Fixed a bug: when input contained complex numbers, in some circumstances MULTIPROD returned an incorrect result. MULTIPROD is defined as a multiple matrix algebra multiplication over both the real and complex fields.

27 Feb 2009 1.30

Warning. New version available as soon as technical problems are solved by The Mathworks.

27 Feb 2009 1.33

Version 2.1 (2009) exploits a quicker, more efficient, and more powerful double-kernel engine. It also introduces virtual array expansion, which allows you, for instance, to multiply a single matrix by an array of matrices.

07 Mar 2009 1.35

Version 2.1 (2009), with array expansion enabled, and memory efficient double-kernel engine. Updating manual and adding screenshot.

26 Jul 2010 1.37

Version 2.1 (2009) - Updating manual and description. MULTIPROD operates "block-by-block" matrix multiplication, on block arrays. It is a generalization of both matrix multiplication (*) and element-by-element multiplication (.*)