How can I improve the speed of the following matrix multiplication?

My problem is doing the following calculation. Suppose the size of the matrix is as follows:

Mat A: [100 X 500]
Mat B: [500 X 1]
Mat C: [500 X 200]

All matrices are not sparse. What I want to calculate is

for i = 1:n
  D = A*(B.*C)
end

Here, (B.*C) creates a matrix having a size of [500 X 200]. The problem is elements of the Mat B are changed in every loop whereas Mat A and C have constant elements. I think the cpu time for this calculation can be drastically reduced by calculating Mat A and C outside of the for-loop first. (Now, a matrix multiplication between Mat A, B, C is done in every loop again and again even though Mat A and C are matrices having constant elements.) But I don't know what I should do.. Can you please help me?

The Matlab code is

A = rand(100, 500);
C = rand(500, 200);
n = 1000;
for i = 1:n
B = rand(500, 1);
D(:,:,i) = A*(B.*C);
end

[Stephen23]

"How can I improve the speed of the following matrix multiplication?"

Probably just write the code inside a function and let the MATLAB JIT engine do any optimization that it can.

Have you profiled the entire code? Is this really the major bottleneck?

[HYDYN]

Thank you for your comment. The following is the code. (not exact but similar)

A = rand(100, 500);
C = rand(500, 200);
n = 1000;
for i = 1:n
B = rand(500, 1);
D(:,:,i) = A*(B.*C);
end

[HYDYN]

Actually, this code is used for the nonlinear ODE solver. This calculation is iterated more than 1,000,000 times for 1 simulation. This is why I'm trying to reduce the time taken for this calculation.

[HYDYN]

I tried to use this kind of additional toolbox.

https://kr.mathworks.com/matlabcentral/fileexchange/25977-mtimesx-fast-matrix-multiply-with-multi-dimensional-support

but, it gave me many error messages to me. I think it's because of version problem or anything like that.

[dpb]

Anything from James is bound to be good; I've not use the mtimesx package, however, so I can't comment on it. I suspect whatever issues you have are those owing to installing and building the necessary libraries which may be somewhat challenging since the mex setup evolves with releases of Matlab and one has to rebuild the package to run it.

[dpb]

" elements of the Mat B are changed in every loop whereas Mat A and C have constant elements. I think the cpu time for this calculation can be drastically reduced by calculating Mat A and C outside of the for-loop first. (Now, a matrix multiplication between Mat A, B, C is done in every loop again and again even though Mat A and C are matrices having constant elements.)"

I don't see why you think it should matter that A, C are constant; when you change B then the products must be recomputed.

Answer 1

[Bruno Luong]

AC = reshape(reshape(A,size(A,1),1,[]).*reshape(C.',1,size(C,2)[]),[],size(B,1));
for i = 1:n
  reshape(AC*B,size(A,1),[]); % = A*(B.*C)
end

[HYDYN]

Thank you for your comment. But, have you run the code you just attached? I just tried to run it but it failed. Could you please upload the entire code? Thank you.

[Bruno Luong]

What MATLAB version you have, if old version you have to replace

the first operation auto-expansion operation ... .* ...

by

bsxfun(@times, reshape(A,size(A,1),1,[]), reshape(C.',1,size(C,2),[]))

[HYDYN]

I'm using the latest version.

AC is successfully calculated but AC*B failed. Because AC is a 3-D matrix whereas C is a vector.

Could you upload the entire matlab code you used? It will be the most straightforward way I think.

[Bruno Luong]

I just correctly edit my code above

[HYDYN]

Thank you so much for your effort. I've checked the code you made gives me the exact solution. But unfortunately, it takes about 3~4 times longer time to complete the calculation.

I think it's because "reshape" function spends quite a long time.

But thank you so much for your effort again.

[Bruno Luong]

You are wrong, reshape doesn't take any significant time.

[dpb]

I'm confused, Bruno, ...

for i = 1:n
  reshape(AC*B,size(A,1),[]); % = A*(B.*C)
end

has no dependency on i in the loop????

[Bruno Luong]

Read OP's question: B changes within the loop, meaning depends on i

[dpb]

Not in that loop it doesn't... :)

If OP added something to it, ok, that could make it be different...

[dpb]

", it takes about 3~4 times longer time"

It would be interesting to see what the profiler says about who's the culprit in that case.

Is your code written as a function or a script? While improving, JIT compilation/optimization works "more better" for functions over scripts so at least, as Stephen C suggested first off, give it all the help you can by wrapping in a function if it isn't.

[HYDYN]

I think you r right in many cases.

But in this case, I've run both version of codes (original one and the one proposed by Bruno) in a script and also function.

Unfortunately, the original one is much faster than Bruno's.

Because "reshape" function spent a much longer time than MATLAB (*) or (.*) functions.

[Walter Roberson]

You are mistaken about the time required for reshape. reshape() is the fastest function in MATLAB other than just a plain assignment. It just checks to be sure that the number of elements of the new shape matches the number in the original, and then it copies the data pointer and puts in the new dimensions and increments the reference count in the data pointer.

[Bruno Luong]

Reshape virtually takes no time, as Walter and I told you. The extra time is I suspect memory accessing of much larger array, which cannot be cache. So even my code requires less arithmetic operations (about 1/3 less) it is slower. You can profile by separate RESHAPE and MTIMES to check what we told you. But I also found the extra CPU time is quite a bad surprise.

[Bruno Luong]

To deal with sparse matrix, see here.