Improvement the speed of array filling

I am performing some intensive array work for a simulation . Esentially , I have a 40X27000 (X) array which at first is filled with initial values and after doing some operations is updated , all this procedure is done inside a for loop . After doing profiling , it can be seen that the lines of code used for updating the values are considerably slow.

(As I am not providing all the code , please consider that not described variables are auxilliary variables of compatible size that come from intermediate computing steps )

Inside the for loop I have considered the following options I have not noticed any improvement :

Option A

A = (MM-X(8,:))./VV;
B =  (LL-X(3,:))./FF;
C = G.*(-I)./(1.0+(1.5./X(9,:)).^8.0); %This is slow , but expected

%Continues for C D E ...%
X = X + W .* [A; B; C; D; E] %This last array [A; B; C; D; E] is made up of 40 elements , just showing some for simplicity purposes.

Option B

A = (MM-X(8,:))./VV;
B =  (LL-X(3,:))./FF;
C = G.*(-I)./(1.0+(1.5./X(9,:)).^8.0); 

%Continues for C D E ...%

X(1,:) = X(1,:) + W .* A 
X(2,:) = X(2,:) + W .* B 
X(3,:) = X(3,:) + W .* C 
X(4,:) = X(4,:) + W .* D 

% ... Continues for all variables %

Option C

X(1,:) = X(1,:) + W .* (MM-X(8,:))./VV;
X(2,:) = X(2,:) + W .* (LL-X(3,:))./FF;
X(3,:) = X(3,:) + W .* G .*(-I)./(1.0+(1.5./X(9,:)).^8.0); 

% ... Continues for all variables%

Furthermore , I have tried initializing X as a gpuArray() without any improvement . Thanks for your help

[dpb]

"...I am not providing all the code ,..."

Too little to be able to do much, methinks...as is, nothing on the RHS shows any dependence on the loop variable...it could by what we see all be moved outside the loop. Oftentimes, something as simple as moving invariant calculations outside the loop is all it takes...

Couple other observations that also probably won't make a huge difference, but...

C = G.*(-I)./(1.0+(1.5./X(9,:)).^8.0); %This is slow , but expected

replace the floating point exponent with integer...I don't know if MATLAB optimizer is good enough to recognize it or not, but fractional exponents are more expensive to calculate than integer ones.

X = X + W .* [A; B; C; D; E]

I'd've thunk this would be expensive with the [] operation -- but since it isn't stored and operated on, maybe the optimizer can create the temporary by simple addressing if removing it really didn't make any appreciable difference.

[Jan]

Please post some code, which we can run by copy&paste. Without the dimensions of the used variables it is hard to compare the results or to improve the speed.

[Enrique]

Thanks for your help , I definetively understand your points of view . I tried to give as much context I can because the original file is 800+ lines with hundreds of auxiliary variables . The key of my question is how to optimize this kind of array updating inside a loop.

I can provide you a toy example with more context regarding the dimensionality of the problem .

size = 15000 % On the real problem this should be 30000

X = rand(45,size);
MM = rand(1,size);
VV = rand(1,size);
LL = rand(1,size);
FF = rand(1,size);
I = rand(1,size);
G = 0.823;
W = 5e-3;
t = 0 ;


while t < 10000 % On the real problem this should be 1e6

A = (MM-X(1,:))./VV;
B =  (LL-X(3,:))./FF;
C = G.*(-I)./(1.0+(1.5./X(2,:)).^8.0);

% Problematic Line %
X = X + W .* [A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C;A; B; C;A; B; C;] ; 
% ----------------- %
t = t +1 ;

end

Regarding the line of X update , please keep in mind that I am repeating [A B C ... A B C ] for explanatory purposes , but on the orignal problem each of the 40+ variables is different [A B C D E ... Z ... A1 ... Z2] (I am not looking for a repmat() style solution)

Thanks again

[dpb]

AHA! My suposition was, in fact, correct -- the building of the array from the individual vectors IS the bottleneck -- 65% of the code time is spent there.

See what the profiler does with

tmp=[A;B;C;A;B;C;A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C;A; B; C;A; B; C;] ; 
X = X + W .* tmp;

instead. I'm betting the time is now going to be spent in the line that builds tmp, not in the computation using it.

I've gotta' run now, but the key is then to NOT do that -- instead create/preallocate a 2D array and index into it when you create A, B, C, ... then you've already built the array instead of copying it as the present code does.

[Enrique]

Yes , if it is used tmp as you explained the bottleneck will be there .

With your last recommendation I understand that Option C will be the best (as X is preallocated ) or there is something better ..?

[Walter Roberson]

Also, you would improve performance by working with the transpose of your array, so that you are working on columns instead of rows. Memory for any individual column is consecutive and so faster.

[dpb]

It's tmp that's the bottleneck, not X, it is the one that needs preallocation and direct addressing into when each row/column(*) is defined.

If you were to also follow Walter's point regarding internal array storage order and reorient as column vectors (thus writing contiguous memory as much as possible), it would seem most probable, yes. I was thinking about A with the large array preallocated first. Not sure there whether there would/would not be difference between the two; would have to test to see; the larger memory footprint might override the internal array operations.

@Matt J's idea of the cell array is interesting to do similar with the "divide and conquer" idea although I'd think the cell addressing would add a little over the direct array addressing mode if the vectors were, again, column oriented. Again, probably only profiling/timing will really prove it one way or t'other; you can't predict the outcome of the optimizer's output.

Whether this would end up being a candidate for mex'ing I dunno, either...possibly.

Answer 1

[Matt J]

Furthermore , I have tried initializing X as a gpuArray() without any improvement.

There should be some improvement. If there isn't, the bottleneck is somewhere else.

% Problematic Line %

X = X + W .* [A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C; A; B; C;A; B; C;A; B; C;] ;

Do you really have such substantial repetition in the update matrix A,B,C,A,B,C... If so, that can be exploited.

[m,n]=size(X);

X=reshape(X,3,n,[]) + W.*[A;B;C]

If not, and if you really do have 40 different expressions for all 40 rows, (I would find it strange/surprising to have so many), the fastest thing would then be to keep the rows of X as a cell array:

X{1} = X{1} + W .* (MM-X{8})./VV;
X{2} = X{2} + W .* (LL-X{3})./FF;
X{3} = X{3} + W .* G .*(-I)./(1.0+(1.5./X{9}).^8.0); 

[Enrique]

Thanks for your answer .Yes! in fact there are 45 different variables . Regarding your second solution , What are the differences in terms of performance between Cell Array and the standard array ?

[Matt J]

It depends both on the number of cells and how big their contents are. Compare:

X=rand(40,270000);
dX=rand(40,270000);

tic;
for i=1:40
 X(i,:)=X(i,:)+dX(i,:);
end
toc;

Elapsed time is 0.193364 seconds.

X=num2cell(X,2);
dX=num2cell(dX,2);

tic;
for i=1:40
 X{i}=X{i}+dX{i};
end
toc;

Elapsed time is 0.011586 seconds.

[dpb]

That seems an issue with code generator...what is the reason for there being any other overhead in indexing through the full array as opposed to indexing to a set of arrays sequentially?

Something there justa ain't right.

But, certainly looks as though the way code is currently generated should fix OP's problem...

Thought I'd look at the memory order thingie, too...

X=rand(40,270000);
dX=rand(40,270000);

tic;
for i=1:40
 X(i,:)=X(i,:)+dX(i,:);
end
toc;
Elapsed time is 0.357787 seconds.
>> 

>> X=X.';
dX=dX.';

tic;
for i=1:40
 X(:,i)=X(:,i)+dX(:,i);
end
toc;
Elapsed time is 0.077559 seconds.
>>

so, indeed, the difference in memory access order also is significant it seems...

[Matt J]

That seems an issue with code generator...what is the reason for there being any other overhead in indexing through the full array as opposed to indexing to a set of arrays sequentially?

Because memory is allocated by the operation X(i,:) on the right hand side of,

 X(i,:)=X(i,:)+dX(i,:);

In the cell array version, there is no memory allocation. Similarly, in the original code, memory is allocated by the concatenation operation [A;B;C;D;...] whereas in the cell array version, the overhead of concatenation is avoided.

[dpb]

Ah....but that's an optimization missed, the target memory is already allocated.

[Bruno Luong]

I agree with @dpb, if X is not shared (meaning no Y = X somewhere upstream), there is no memory allocation due to (inplace optimization) during the statements

X(i,:) = somerow;

or

X(:,j) = somecolumn;

The speed difference is from contiguous or non-contiguous access of memory.

[Matt J]

It's definitely something to ponder, however a side-by-side timing test of all 3 implementations still shows cell array access to be significantly faster than column-wise access.

X=rand(40,27e5);
dX=rand(40,27e5);

tic;
for i=1:40
 X(i,:)=X(i,:)+dX(i,:);
end
toc;

Elapsed time is 3.178377 seconds.

Xt=X';dXt=dX';
tic;
for i=1:40
 Xt(:,i)=Xt(:,i)+dXt(:,i);
end
toc;

Elapsed time is 0.628279 seconds.

X=num2cell(X,2);
dX=num2cell(dX,2);

tic;
for i=1:40
 X{i}=X{i}+dX{i};
end
toc;

Elapsed time is 0.237045 seconds.

[Matt J]

OK. Need to put the code inside a function to see full optimization:

runTest()

Elapsed time is 0.706534 seconds.
Elapsed time is 0.070628 seconds.
Elapsed time is 0.089183 seconds.

Still, though, it's important to realize there will be scenarios where in-place updates are not optimized and cell array access turns out to be better, like in this second test:

runTest2()

Elapsed time is 2.356626 seconds.
Elapsed time is 0.898723 seconds.
Elapsed time is 0.276424 seconds.

function runTest()

    X=rand(40,27e5);
    dX=rand(40,27e5);
    
    tic;
    for i=1:40
     X(i,:)=X(i,:)+dX(i,:);
    end
    toc;
    
    Xt=X';dXt=dX';
    tic;
    for i=1:40
     Xt(:,i)=Xt(:,i)+dXt(:,i);
    end
    toc;
    
    X=num2cell(X,2);
    dX=num2cell(dX,2);
    
    tic;
    for i=1:40
     X{i}=X{i}+dX{i};
    end
    toc;
    
end

function runTest2()

    X=rand(40,27e5);
    dX=rand(40,27e5);
    
    I=speye(27e5);
    
    tic;
    for i=1:40
     X(i,:)=X(i,:)*I;
    end
    toc;
    
    Xt=X';dXt=dX';
    tic;
    for i=1:40
     Xt(:,i)=I*Xt(:,i);
    end
    toc;
    
    X=num2cell(X,2);
    dX=num2cell(dX,2);
    
    tic;
    for i=1:40
     X{i}=X{i}*I;
    end
    toc;
    
end

[Bruno Luong]

"it's important to realize there will be scenarios where in-place updates are not optimized "

I don't think your second test show the penalty is due to inplace. IMO the speec reduction is because that the Xt(:,i) on the RHS has to be extracted from the original vector to form a new vector for matrix x vector multiplication. Matlab JIT does not go yet into inplace-pointer in operation/function. Where as the cell array loop, no such extraction is needed.