Help to make a code faster/memory efficient

Question

Fernando 2017-9-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/356087-help-to-make-a-code-faster-memory-efficient

编辑： Cedric 2017-9-12

Dear all, I am trying to improve the computation speed/memory usage of a code I've done for a electromagnetism problem. I have converted the arrays to single values which has helped a bit and also tried converting to gpuArrays but this seems slower. I've copied the most time consuming part of my code, other parts are coded similarly so any suggestion in this piece of code would be helpful. It seems the cross product calculation consumes a lot of time, I think maybe using something like pagefun, arrayfun? but I have no idea how to convert this to that. Any ideas would be deeply appreciated.

 magnetic=zeros(19128960,3,'single');
 for u = 1:2592 
    for n = 1:7380 
    q =1+(n-1)+(7380*(u-1));
    magne1=cross((PosS(n,:)-PosSeg1(u,:)),(ExtPosSeg1(u,:)-PosSeg1(u,:)));
    magnetic(q,1)=magne1(:,1);
    magnetic(q,2)=magne1(:,2);
    magnetic(q,3)=magne1(:,3);
    end
 end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Cedric 2017-9-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/356087-help-to-make-a-code-faster-memory-efficient#answer_281104

编辑：Cedric 2017-9-11

在 MATLAB Online 中打开

Check this out:

 n1         = 7380 ;
 n2         = 2592 ;
 PosS       = rand( n1, 3 ) ;
 PosSeg1    = rand( n2, 3 ) ;
 ExtPosSeg1 = rand( n2, 3 ) ;

Loop-based:

 tic ;
 magnetic=zeros(n1*n2, 3,'single');
 for u = 1:n2 
    for n = 1:n1
        q =1+(n-1)+(n1*(u-1));
        magne1=cross((PosS(n,:)-PosSeg1(u,:)),(ExtPosSeg1(u,:)-PosSeg1(u,:)));
        magnetic(q,1)=magne1(:,1);
        magnetic(q,2)=magne1(:,2);
        magnetic(q,3)=magne1(:,3);
    end
 end
 toc

outputs

Elapsed time is 86.585071 seconds.

and here is a vector version

 tic ; 
 magnetic_CW = reshape( permute( cross( ...
    bsxfun( @minus, permute( PosS, [3,2,1] ), PosSeg1 ), ...
    repmat( ExtPosSeg1 - PosSeg1, 1, 1, n1 )), [2,3,1] ), 3, [] ).' ;
 toc

outputs

Elapsed time is 2.267693 seconds.

Both outputs are equal:

 isequal( magnetic, single( magnetic_CW ))
 ans =
  logical
   1

No need to convert to single, however, I did it for the comparison.

Cheers,

Cedric

NOTES

It looks awfully complicated because I wanted to avoid creating intermediary variables, but what we do is just a cross product on two 3D arrays, one containing PosS-PosSeg1 for all columns of both (hence the call to BSXFUN), and the second being the corresponding 3D array of ExtPosSeg1-PosSeg1. All the rest is about permuting/reshaping/repeating so dimensions are appropriate.
You may not need to work with 2D arrays actually, and you can probably eliminate the external RESHAPE and PERMUTE operations.

5 个评论
显示 3更早的评论隐藏 3更早的评论

Cedric 2017-9-11

在 MATLAB Online 中打开

If you split the one liner and create a few intermediary variables, to understand:

 A = bsxfun( @minus, permute( PosS, [3,2,1] ), PosSeg1 ) ;
 B = repmat( ExtPosSeg1 - PosSeg1, 1, 1, n1 ) ;
 C = cross( A, B ) ;
 magnetic_CW = reshape( permute( C, [2,3,1] ), 3, [] ).' ;

you see that we use the ability of CROSS to work on 3D arrays. What comes before is for building two 3D arrays with appropriate content, and what comes after is for reshaping the output so it matches the shape of your original output.

These operations take time, not as much as compute the cross product but still. If you look inside the CROSS function ( open cross ) you will see more calls to (I)PERMUTE and RESHAPE.

What happens is quite a few conversions 2D->3D->2D->3D... Often, the fastest computation that you can get is by working on the largest possible arrays (so MALTAB and underlying libs can optimize). The point of note 2 is just that the 3D array C may be already shaped in a way that is advantageous for the next computation, so you may want to evaluate if it is worth reshaping is back to your original shape at this stage, or if it is better to perform further computations first and reshaping later.

Cedric 2017-9-11

编辑：Cedric 2017-9-12

在 MATLAB Online 中打开

The most complicated part to understand is the call to BSXFUN. To better understand, look at this:

 >> bsxfun( @times, 1:5, [10;20;30] )
 ans =
    10    20    30    40    50
    20    40    60    80   100
    30    60    90   120   150

where the first array is 1x5 and the second 3x1. You can see what the implicit expansion does (here when the operation is a product, but you can pass handles to any appropriate functions).

In your case (think in 3D and xyz axes) I rotate PosS that is 7380x3x1 (in the xy plan) around the y axis (by permuting dims 1 and 3), so it becomes vertical (in the yz plan). So now I have the rotated/permuted PosS that is 1x3x7380 and PosSeg1 that is 2592x3x1 (in the xy plan) and I can perform an expansion (here the operation is minus). The expansion "fills the volume defined by the two arrays" in a sense.

The outcome is a 2592x3x7830 array with all relevant differences.

The next step is to create the corresponding 3D array of ExtPosSeg1-PosSeg1 ..

Fernando 2017-9-12

Oh I understand now.

Thanks a lot, I have managed to implement this vectorization in the different calculations of that for loop in my code, now it works 140 times faster in that part. Thanks a lot.