Matrix multiplication-optimal speed and memory
显示 更早的评论
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 个评论
James Tursa
2011-4-28
Are your matrices real or complex?
James Tursa
2011-4-28
Also, what are the actual sizes of the matrices you are using?
Ronni
2011-4-28
bym
2011-4-29
is sparse an option? what type of problem is it?
Ronni
2011-4-29
Oleg Komarov
2011-4-29
What are you trying to do A for? Can you post more code so that we can see the final point you want to get to.
采纳的回答
James Tursa
2011-4-29
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 个评论
James Tursa
2011-4-29
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.
更多回答（6 个）
Oleg Komarov
2011-4-28
sum(D.^2.*diag(W))
or
(D.*diag(W)).'*D
diag(W) can be replaced by the vector containing the diagonal elements.
James Tursa
2011-4-28
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.
Ronni
2011-4-28
0 个投票
2 个评论
James Tursa
2011-4-29
Your comment doesn't make sense. How can you do A=D' if you first clear D? If you start with a 9GB D and then do A=D' you will get 18GB total since a data copy takes place. If you then clear D the memory should drop back to 9GB. Are you saying this doesn't happen?
Ronni
2011-4-29
James Tursa
2011-4-29
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 个评论
James Tursa
2011-4-29
A thought just occurred to me that you might be able to coax MATLAB into doing the D = sqrt(W) * D operation in-place if you use the bsxfun formulation shown above and you put it inside of a function. I will have to check on this.
Ronni
2011-4-29
James Tursa
2011-4-29
Update: I just checked on R2011a and the JIT is not yet smart enough to do the sqrt or the bsxfun(@times,etc) operations in-place. So it looks like mex is the only option for this at present.
类别
在 帮助中心 和 File Exchange 中查找有关 Operators and Elementary Operations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
