My mex file is slower than my original matlab equivalent

Hello friends,

I need to calculate some quantities of linear algeibra type, so they are merely matrix and vector products. The following is an example

EZ=[(1.0./Ds0.^2.*(Ds0.*(Dm0.*4.0+Dm0.*Ds1.^2.*2.0-Ds0.*(Ds1.*2.0+Dm0.*Ds2.*2.0+Dm1.*Ds1.*2.0-Dm2.*Ds0)+Dm0.*Dm1.*2.0)-Dm0.^2.*Ds1.*2.0))./4.0;
    (1.0./Ds0.^2.*(Ds0.*(Ds0.*(Dm1.*4.0+Ds1.^2-Ds0.*Ds2.*2.0+4.0)-Dm0.*Ds1.*8.0)+Dm0.^2.*4.0))./4.0;
    (Dm0.*6.0-Ds0.*Ds1.*3.0)./(Ds0.*2.0)];

where Ds0,Ds1,Ds2,Dm0,Dm1,Dm2 are 1*n vectors. When I do the calculations using matlabFunction (attached) it is fast. However, I am not satisfied since I really need to do such calculations thousands of time s(if not millions of times). To overcome this issue I decided to give mex a try. Unfortunately, the equivalent mex file (which I made by matlab coder) is slower 2-3 times (I could not upload it here, unforetunately).

Is there any hope to create a mex file out of this function which is much faster? I hope so!

Thanks for your help in advance,

Babak

[s pernot]

that's normal... for standard operations are highly optimized in built in matlab... and mex equivalent is therefore

not the solution to improve the cputime. in addition mex function add an overhead time call. so in your case you should explore the use of gpuArray if you have an Nvidia gpu card...

[Walter Roberson]

MATLAB Coder does not generate calls to the high performance multi-threaded library, so you should expect the code to be slower for linear algebra.

[Bruno Luong]

What linear algebra? All I see is element wise construction of array.

It looks awfully written as a longline copied from some unworking symbolic formula withouttany care about simplification or factorizing common terms.

[Mohammad Shojaei Arani]

Thanks for your comment about using "gpuArray". I did test it but got a surprisingly opposite results : with "gpuArray" my calculations become actually much slower!!! Why??? If you run the following lines you see that in this case using "gpuArray" computing time becomes more than 200 times bigger !!! Any explanation?

clc;
n=100;N=10000;
Dm0=rand(1,n,"gpuArray");Dm1=rand(1,n,"gpuArray");Dm2=rand(1,n,"gpuArray");Ds0=rand(1,n,"gpuArray");Ds1=rand(1,n,"gpuArray");Ds2=rand(1,n,"gpuArray");
tic;
for i=1:N
    EZ = EZ_QQ32(Dm0,Dm1,Dm2,Ds0,Ds1,Ds2);
end
t1=toc;
Dm0=rand(1,n);Dm1=rand(1,n);Dm2=rand(1,n);Ds0=rand(1,n);Ds1=rand(1,n);Ds2=rand(1,n);
tic;
for i=1:N
    EZ = EZ_QQ32(Dm0,Dm1,Dm2,Ds0,Ds1,Ds2);
end
t2=toc;
t1=t1/N;t2=t2/N;
vpa([t1 t2 t2/t1])

[0.00055874864999999997519641858900741, 0.0000022361800000000001622147034752164, 0.0040021215263786324989325748902047]

[Mohammad Shojaei Arani]

Hi Bruno,

I agree with you that I could make it (EZ_QQ32) much more simplified by using techniques of common subexpression elimination, factorization, using horner schemes, etc. But, this is not the point in my question. Please understand my point. My question still applies to any simplified expression you can get from my current longer expression.

Why "gpuArray" does not reduce the computing time for my expression or any simplified version of it?

[Bruno Luong]

So you don't think simplying expression matters? I do think the contrary. Everytime in the expression there is a gpu array involves there is a whole transfer data from cpu to gpu, you mighte have 200 such terms in your expression, I don't even try to count or understand your code as it is a so unreable and messy expression.

If you a raw unsimplified expression like yours, throw it in the computer and ask why it doesn't accelerate, you need to think a much more lower level how it works.

[Mohammad Shojaei Arani]

Bruno,

Of course, simplification matters a lot. My actual expressions are way longer than this. Matlab is not able to simplify them in an efficint way (and in many cases it simplifies a little). I have spent a lot of time on how to simplify my expressions. Unfortunately, using matlab I do not have any hope to simplify my expressions more than this (yes, you can perhaps simplify this expression more because it is not extremely long but can you do it for an expression which is 1KM long???) My expressions are in rational form. So, typically I perform 2 operations to simplify them: 1) first I apply [n,d] = numden(EZ), and then 2) EZ = horner(n,Ds0)./horner(d,Ds0). Unfortunately, matlab does support a multivariate horner scheme and I could only benifit the univariate horner scheme here (I apply horner scheme with respect to variable Ds0 as it is the most repeated variable. Typically, you should apply horner scheme with respect to such variables). So, at this point I convinced myself that I canot hope to simplify my expressions more using matlab. Therefore, I should find strategies to ask C or C++ to perform the calculations.

So, my question is not about how to come up with a better simplification (as it does not work with the current capacities of matlab). My question is "how can I use C/C++ or perhaps resort to stuff like gpuArray, etc to reduce the computational burden".

Answer 1

[Jan]

Just some experiments. You can gain some clarity, but hardly improve the speed with this simplifications. I've tried a loop version also.

n   = 1e4;
Ds0 = rand(1, n);
Ds1 = rand(1, n);
Ds2 = rand(1, n);
Dm0 = rand(1, n);
Dm1 = rand(1, n);
Dm2 = rand(1, n);

tic;
for rep = 1:1e4
EZ = [(1.0./Ds0.^2.*(Ds0.*(Dm0.*4.0+Dm0.*Ds1.^2.*2.0-Ds0.*(Ds1.*2.0+Dm0.*Ds2.*2.0+Dm1.*Ds1.*2.0-Dm2.*Ds0)+Dm0.*Dm1.*2.0)-Dm0.^2.*Ds1.*2.0))./4.0;
    (1.0./Ds0.^2.*(Ds0.*(Ds0.*(Dm1.*4.0+Ds1.^2-Ds0.*Ds2.*2.0+4.0)-Dm0.*Ds1.*8.0)+Dm0.^2.*4.0))./4.0;
    (Dm0.*6.0-Ds0.*Ds1.*3.0)./(Ds0.*2.0)];
end
toc

Elapsed time is 0.794406 seconds.

tic;
for rep = 1:1e4
   Ds0_2 = Ds0 .* Ds0;
   Dm0_2 = Dm0 .* Dm0;
   EZ2 = [(1 ./ Ds0_2 .* (Ds0 .* (Dm0 * 2 + Dm0 .* Ds1 .^ 2 - ...
          Ds0 .* (Ds1 + Dm0 .* Ds2 + Dm1 .* Ds1 - Dm2 .* Ds0 ./ 2) + ...
          Dm0 .* Dm1) - Dm0_2 .* Ds1)) / 2; ...
          1 ./ Ds0_2 .* (Ds0 .* (Ds0 .* (Dm1 + Ds1 .^ 2 / 4 - Ds0 .* Ds2 / 2 + 1) - ...
          Dm0 .* Ds1 * 2) + Dm0_2);
          (Dm0 * 3 - Ds0 .* Ds1 * 1.5) ./ Ds0];
end
toc

Elapsed time is 0.775081 seconds.

tic;
for rep = 1:1e4
   EZ3 = zeros(3, n);
   for k = 1:n
       a = Ds0(k);
       b = Dm0(k);
       c = Ds1(k);
       d = Dm1(k);
       e = Ds2(k);
     EZ3(1, k) = (1 / a^2 * (a * (b * 2 + b * c ^ 2 - ...
                  a * (c + b * e + d * c - Dm2(k) * a / 2) + b * d) - b^2 * c)) / 2;
     EZ3(2, k) = (a * (a * (d + c ^ 2 / 4 - a * e / 2 + 1) - b * c * 2) + b^2) / a^2;
     EZ3(3, k) = b * 3 / a - c * 1.5;
   end
end
toc

Elapsed time is 1.140882 seconds.

max(abs(EZ(:) - EZ2(:)))

ans = 0

max(abs(EZ(:) - EZ3(:)))

ans = 2.3283e-10

[Bruno Luong]

Well there is room for simplification, e.g.

(a * (a * (d + c ^ 2 / 4 - a * e / 2 + 1) - b * c * 2) + b^2) / a^2

can be replaced by

ia = 1/a;
bia = b * ia;
d + (c*c*0.5 - a*e)*0.5 + 1 + (-2*c + bia) * bia

Funny that people throw the raw stuff into gpu, mex, etc... whereas they can do better (I guess) in term of reducing operations by hand.

[Mohammad Shojaei Arani]

Hi Jan,

You say that "we can hardly improve the speed with these simplifications". But, if you show me a systematic way to do what you did then for me what you did is quite great. Why? Forget about this example. My actual expressions can be much much much longer. In your first experiment the elapsed time was 0.794406 and in your second experiment it was 0.775081. So, 0.0193 seconds performance gain!

If, instead, I would give you an expression which would be 1000 times longer do you know how to apply your idea to this much longer expression? If you teach me how to do this then it is quite a lot of performance gain for me.

So, my question is how can we ask a code to do what you did?

[Bruno Luong]

"If, instead, I would give you an expression which would be 1000 times longer do you know how to apply your idea to this much longer expression?"

You don't even have any warranty that the expresssion you compute using floating points is accurate at any degree with 1000 terms. People takling about Horner rule and it matters for 10 terms.

[Jan]

@Mohammad Shojaei Arani: The rules are straight:

1. Avoid repeated work. If a calculation appears repeatedly, compute it once and store it in a temporary variable.
2. Reduce the call to expensive functions: exp, power, trigonometric functions, faculty, ...
3. Combine operations, but keep in mind, that the result can be influenced by rounding effects. E.g. 1/a*b takes more time than b/a, but the result can be slightly different.

The clarity of the code improves the time needed for debugging:

1. Spaces around operators.
2. Compact names of variables.
3. Be careful with using parentheses, if they are not required.
4. Avoid elementwise operators, if the calculation does not need it. 3.0.*2.0 is harder to read then 3 * 2.

Bruno's point is important: The result of numerically instable functions can be influenced massively by simplifications. A basic example:

1e17 + 1 - 1e17

ans = 0

1e17 - 1e17 + 1

ans = 1

[Mohammad Shojaei Arani]

Thanks a lot Jan and Bruno!