Code covered by the BSD License  

Highlights from
Multiple matrix multiplications, with array expansion enabled

4.95

5.0 | 42 ratings Rate this file 40 Downloads (last 30 days) File Size: 505 KB File ID: #8773
image thumbnail

Multiple matrix multiplications, with array expansion enabled

by

 

20 Oct 2005 (Updated )

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

| Watch this File

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

This file inspired Vector Algebra For Arrays Of Any Size, With Array Expansion Enabled, Frontal, Inversion Every 2 D Slice For Arbitrary Multi Dimension Array., Superkron, and Dot Product Of Two N Dimensional Arrays..

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)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (52)
01 Aug 2014 Mikhail Katliar

Very helpful, exactly what I needed.

23 Jul 2014 Adolfo Herbster  
29 Jun 2014 Manuel Diaz  
27 May 2014 Denis Anikiev  
14 Mar 2013 Romain W

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

Hi Sir

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

14 Nov 2012 Oliver Woodford

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

31 Oct 2012 YIXIAO  
26 May 2012 Mark Flynn

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.

25 May 2012 Paolo de Leva

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.

24 May 2012 Mark Flynn

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

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

11 Apr 2012 Paolo de Leva

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 :-)

10 Apr 2012 Michael Völker

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

10 Apr 2012 A

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

10 Apr 2012 Michael Völker

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.

10 Apr 2012 A

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?

19 Feb 2012 Dinakar

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

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).

03 Feb 2012 Michael Völker

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().

27 Jan 2012 Shatrughan

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

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

09 Jun 2011 Nicolas  
28 Apr 2011 Xiaodong

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

27 Apr 2011 Daniel Armyr

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

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  
08 Apr 2011 Kirankumar  
15 Jan 2011 Michael Völker

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  
05 Aug 2010 Tanyer Alan  
04 Mar 2010 Aleksandra  
03 Nov 2009 Deepak

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

12 Aug 2009 Denis Kozlov  
31 Jul 2009 Paolo de Leva

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

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  
21 May 2009 Noah Snavely

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!

14 Oct 2006 Juan D'Adamo  
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.

Updates
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

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

27 Feb 2009

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

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

26 Jul 2010

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 (.*)

Contact us