Multiplying matrices, vectors, or scalars contained in two ND arrays, with array expansion enabled.
MULTIPROD is a powerful, quick and memory efficient generalization for ND arrays of the MATLAB matrix multiplication operator (*). While the latter works only with 2D 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 builtin function TIMES. While TIMES operates elementbyelement multiplications (e.g. A .* B), MULTIPROD operates blockbyblock 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 multimultiplied 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 rototranslations 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 ND arrays of the matrix multiplication function MTIMES, with AX enabled.
Vector inner, outer, and cross products generalized for ND arrays and with AX enabled are performed by DOT2, OUTER, and CROSS2 (MATLAB Central, file #8782, http://www.mathworks.com/matlabcentral/fileexchange/8782).
Elementbyelement 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.
1.37  Version 2.1 (2009)  Updating manual and description. MULTIPROD operates "blockbyblock" matrix multiplication, on block arrays. It is a generalization of both matrix multiplication (*) and elementbyelement multiplication (.*) 

1.35  Version 2.1 (2009), with array expansion enabled, and memory efficient doublekernel engine. Updating manual and adding screenshot. 

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

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

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. 

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. 

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

The code was optimized according to suggestions automatically generated by function MLint. 

The code was optimized according to suggestions automatically generated by function MLint. 

1. Solved character conversion problem in "File Description" text.


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

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

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. 
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 (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 elementwise multiplications (even via bsxfun).
Added later:
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 (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 elementwise multiplications (even via bsxfun).
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 (view profile)
@cszczyglowski: if M2 is square, you can perform M1/M2 with multiprod and multinv (http://www.mathworks.com/matlabcentral/fileexchange/31222)
Ilan (view profile)
Works great with cuda
cszczyglowski (view profile)
Fantastic bit of code, thank you very much! Is there an equivalent function for performing multiple rightdivide operations in a single step instead of using a loop?
Yuan Zhou (view profile)
This is a great piece of work! It saves a lot of time.
Em Rickinson (view profile)
Andrada (view profile)
Great function! I was wondering if there is a similar function for matrix power. Thanks.
Mikhail Katliar (view profile)
Very helpful, exactly what I needed.
Adolfo Herbster (view profile)
Manuel A. Diaz (view profile)
Denis Anikiev (view profile)
Romain W (view profile)
Hi Paolo,
I was wondering if it is possible without using a forloop 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 (view profile)
Hi Sir
It is really a great function.
Thanks for sharing.
Regards
Oliver Woodford (view profile)
Very useful in general, but there is an overhead which becomes significant when applying this function to smaller arrays.
YIXIAO (view profile)
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 (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 10element 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 (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 (view profile)
I want to multiply an matrix of nxnxn by a nxnxn matrix with the same place.
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 (1D blocks) by a matrix M (2D 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 vectorbymatrix products, that you use the rowvector 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 matrixbyvector 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 email.
If you like multiprod, please rate it. I greatly enjoy receiving fivestar ratings :)
Michael Völker (view profile)
Hi A, then I think that
C = reshape(multiprod( A, B, 4, 1 ), size(A));
will work.
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 (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 (view profile)
I would like to multiply an array of 5element vectors by a 5x5 matrix. In other words, multiply 15x15x15x5 by 5x5 to get 15x15x15x5
can I use multiprod?
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 (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 3D:
A = rand(7,1);
B = rand(9,1);
C = rand(10,1);
X = simple_outer3D(A, B, C);
where X is a 3D 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 3D, 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 3D 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 (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 (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 (view profile)
This is great  it has been a tremendous improvement for my project.
Nicolas (view profile)
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 ND array in a similar way: http://www.mathworks.com/matlabcentral/fileexchange/31222inversionevery2dsliceforarbitrarymultidimensionarray
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 (view profile)
Ambitious documentation, fast implementation, and best of all, it worked outofthe box. The lastes resident in my private toolbox.
Paul (view profile)
Incredible! I have to do a series of structured matrixvector multiplications where the matrices have identical structure. There was no way to do this with builtin matlab functionality, outside of using a FOR loop. A MULTIPRODbased implementation turned out to be 1001000 times faster on problems of practical size!
Paul (view profile)
Kirankumar (view profile)
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 (view profile)
Tanyer Alan (view profile)
Aleksandra (view profile)
Deepak (view profile)
Nice program to change all the lines of code into one command. really useful in simulations etc.
Denis Kozlov (view profile)
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 (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 (view profile)
Noah Snavely (view profile)
This is an amazingly useful tool. Anyone who uses matrices willat some point in their lifebe stymied by a problem for which this software is the perfect solution. It just works.
Cool. Thank You!
Brilliant  this does exactly what I want, performing vector reorientation on a 3D array (i.e., 3D x 3 x 3 rotation matrix multiplied by 3D x 3 vector matrix); very fast.
Well documented, well optimized, and works.
Awesome! Thank you for sharing this useful package
just fantastic, many compliments, it's a great and welldone job!
awesome stuff! just got what i was looking for!
Pretty efficient too!
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.
Just what I was looking for. This capability really should be a builtin for Matlab.
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!
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!
Super!
Have I longed for this functionality in Matlab! Farewell forlooops!
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 forloops now can take advantage of Matlab's inherent fast matrix multiplication algorithms.
Mr deLeva deserves Mathwork's medal in gold for this.
Good job. Thanks to your program got rid of the loop. My code is much faster now. Thanks
Good job. Very useful when manipulating time dependant transformation matrices (vector of 3x3 matrices). No longer need a for loop.