# Matrix multiplication-optimal speed and memory

31 views (last 30 days)

Show older comments

I have large matrices that I'm multiplying, not only is it very slow, it also runs out of memory. I have to do this for many iterations so I really need to find a better way to implement it.

Basically I have A=D'*W*D W is a diagnoal matrix with dimensions MxM and D has dimensions M*N where M>>N.

If it's unclear, someone else had the same problem and he has a nice diagram to illustrate what the problem is. Could someone provide suggestions on how to improve this.

http://stackoverflow.com/questions/4420235/efficient-multiplication-of-very-large-matrices-in-matlab

Thanks.

##### 6 Comments

Oleg Komarov
on 29 Apr 2011

### Accepted Answer

James Tursa
on 29 Apr 2011

Here is the mex routine. Comments should explain its intended use. You will need to install a C compiler as stated earlier. Then do this:

mex -setup

(then press Enter)

(then enter the number of the C compiler)

(then press Enter again)

mex sqrttimes.c

Now you are ready to use it. I didn't bother to parallelize this with OpenMP, but if you happen to install an OpenMP capable compiler then it will be very easy to put a #pragma omp parallel for in front of the loops (with appropriate private clauses) to potentially get more speed.

/*****************************************************************************/

/* sqrttimes.c */

/* Syntax: sqrttimes(A,B) */

/* Does the following: */

/* A = sqrt(A) in-place first, then */

/* B = diag(A) * B in-place second. */

/* Where A = a double vector of size 1 x M or M x 1 */

/* B = a double matrix of size M x N */

/* */

/* The idea is to call this function as follows: */

/* */

/* >> sqrttimes(A,B); */

/* >> C = B' * B; */

/* */

/* The resulting C matrix will be the same as if you had done this: */

/* */

/* >> C = B' * diag(A) * B */

/* */

/* But using sqrttimes will result in practically 0 extra temporary memory */

/* usage so it is very memory efficient. This can be important if the A and */

/* B matrices are very large. Note that the original A and B matrices have */

/* been changed in-place, which is necessary to avoid extra memory usage. */

/* */

/* Programmer: James Tursa */

/* Date: 2011-Apr-28 */

/*****************************************************************************/

#include <math.h>

#include "mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])

{

double *Apr, *Bpr, *Apr0;

mwSize i, j, k, l, m, n;

/* Check inputs */

if( nrhs != 2 ) {

mexErrMsgTxt("Need exactly two inputs.");

}

if( nlhs != 0 ) {

mexErrMsgTxt("Too many outputs.");

}

if( !mxIsDouble(prhs[0]) || !mxIsDouble(prhs[1]) || mxIsSparse(prhs[0]) || mxIsSparse(prhs[1]) ) {

mexErrMsgTxt("Inputs must be double non-sparse.");

}

if( mxGetNumberOfDimensions(prhs[0]) != 2 || mxGetNumberOfDimensions(prhs[1]) != 2 ) {

mexErrMsgTxt("Inputs must be 2D.");

}

m = mxGetM(prhs[0]);

k = mxGetN(prhs[0]);

l = mxGetM(prhs[1]);

n = mxGetN(prhs[1]);

if( m != 1 && k != 1 ) {

mexErrMsgTxt("1st Input must be a vector.");

}

m = m * k;

if( m != l ) {

mexErrMsgTxt("1st vector input must have same number of elements as 1st dim of 2nd input.");

}

/* Get pointers to the data */

Apr = Apr0 = mxGetPr(prhs[0]);

Bpr = mxGetPr(prhs[1]);

/* Do A = sqrt(A) in-place */

for( i=0; i<m; i++ ) {

*Apr = sqrt(*Apr);

Apr++;

}

/* Do B = diag(A)*B in-place */

for( j=0; j<n; j++ ) {

Apr = Apr0;

for( i=0; i<m; i++ ) {

*Bpr++ *= *Apr++;

}

}

}

##### 1 Comment

James Tursa
on 29 Apr 2011

P.S. I need to caution you that sqrttimes modifies its input arguments in-place, meaning that if they are sharing data with another variable then the other variable will be changed also. You may not want this behavior. If you don't know what this means then do this test:

>> A = 1:3

A =

1 2 3

>> B = reshape(1:15,3,5)

B =

1 4 7 10 13

2 5 8 11 14

3 6 9 12 15

>> a = A % a is now sharing data memory with A

a =

1 2 3

>> b = B % b is now sharing data memory with B

b =

1 4 7 10 13

2 5 8 11 14

3 6 9 12 15

>> sqrttimes(A,B)

>> A

A =

1.0000 1.4142 1.7321

>> a

a =

1.0000 1.4142 1.7321

>> B

B =

1.0000 4.0000 7.0000 10.0000 13.0000

2.8284 7.0711 11.3137 15.5563 19.7990

5.1962 10.3923 15.5885 20.7846 25.9808

>> b

b =

1.0000 4.0000 7.0000 10.0000 13.0000

2.8284 7.0711 11.3137 15.5563 19.7990

5.1962 10.3923 15.5885 20.7846 25.9808

You can see that by changing A and B in-place, we have also changed a and b in-place. This is the side-effect of doing in-place operations. So if you don't want this behavior you need to be careful to not share A and B with any other variables.

### More Answers (6)

Oleg Komarov
on 28 Apr 2011

sum(D.^2.*diag(W))

or

(D.*diag(W)).'*D

diag(W) can be replaced by the vector containing the diagonal elements.

##### 0 Comments

James Tursa
on 28 Apr 2011

I don't understand why the bsxfun-based proposed solutions listed in the link do not work for you. I would not expect this to run out of memory if you are implementing it correctly unless you are up against the limit with just W and D. e.g.,

M = large number

N = small number

D = M x N matrix

W = 1 x M vector (i.e., diagonal stored as a vector)

bsxfun(@times,D',W) * D

or

D' * bsxfun(@times,W',D)

The latter will likely result in D' not being explicitly calculated. Rather, a transpose flag will simply be passed to the dgemm routine in the BLAS call. And the W' in the latter will not result in a data copy, it will be just be a simple reshaped shared data copy of the W vector.

##### 0 Comments

James Tursa
on 29 Apr 2011

OK, so how about his approach (assumes W elements are non-negative):

D' * W * D

= D' * sqrt(W) * sqrt(W) * D

= D' * sqrt(W') * sqrt(W) * D

= (sqrt(W) * D) ' * (sqrt(W) * D)

Then do the sqrt(W) * D operation in-place on D in a mex routine. i.e., do D = sqrt(W) * D in-place inside the mex routine. Should be pretty fast and no extra memory requirement (could even parallelize it easily with OpenMP if you had a compatible compiler). Then do this operation:

A = D' * D

That will call a BLAS routine with a transpose flag to do the multiply, so no explicit transpose D' is formed and MATLAB will recognize the symmetry and call a specialized BLAS routine that is faster than a generic matrix multiply.

The big IF in all of this is: Can you tolerate changing D in-place? i.e., do you need the original D explicitly after this operation? If so, you could recover it by undoing the sqrt(W) * D operation in-place with another mex routine.

Let me know how this sounds to you and if it looks like it will work for you I can write up the mex routine and post it here. You will need a 64-bit C compiler to do this. Since one does not ship with MATLAB you will have to install it yourself. e.g., see this thread:

##### 3 Comments

James Tursa
on 29 Apr 2011

### See Also

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!