TRANSPOSE(.') / RESHAPE too much TIME CONSUMING--> Looking for ALTERNATIVE

HI guys,

I wrote a function that uses th sum of a raw and the transpose of a colomn in a triple loop.

As my data is quite large, this line alone consume 904s (I did a run and time) or 75% of the total time.

I am looking for an alternative. Here is the code:

for t = 1:24
    for i = 1:1858
        for j = 1:1858
            if i == j
                D(i,j) = 0;
            else
                C = A(i,:) + permute(B(:,j),[2,1]); % 904 s or 75% of the total (I removed part of the code)  <---------
            end  
        end
    end
end

I tried the transpose function .' and the permute function. Both gave high running time which is a real problem in my study.

I really don't understand how a sum and a transpose consume so much time.

Please help me :)

[madhan ravi]

Please provide a simple example and the desired output , you seem to be mixing loop and permute (havn't examined thoroughly).

[Adam]

Can you not permute your entire DistPT array just once before the loops?

Answer 1

[Guillaume]

I fund it really weird that, for a computer, selecting a simple raw from a matrix is that time consuming

It's always a lot slower than accessing columns in a loop, particularly on modern processors.

Elements of a column are stored together in memory, so accessing a single column is just a straight contiguous memory access. With modern processors, it will be even more optimised as the processor will guess that you'll want to access the next column and prefetch the data before you even request it (read up on processor caches to learn more). With a row access, the elements are scattered across memory. That's not cacheable and each element has to be fetched one by one. Repeat that on a lot of rows, and yes, you're wasting time.

[Habib]

You are so right! Thanks, I didn't know that Matlab has some confusions if you talk raws instead of columns.

Answer 2

[Jan]

Is D pre-allocated? This is an important detail if runtime matters.

permute(DistPT(:,j),[2,1]) ?! What about: DistPT(:,j).' ?

nDist    = size(DistRoad, 1);       % Nicer to do this once
D        = zeros(nDist, nDist);     % Pre-allocate!!!
E_roadPT = zeros(nDist, timeSteps); % Pre-allocate!!!
for t = 1:timeSteps
    for i = 1:nDist
        Ai = A(i, :);
        for j = 1:nDist
            if i ~= j
                D(i,j) = min(Ai + DistPT(:,j).');
            end
        end
    end
    add           = DistFF ./ D;
    E_roadPT(:,t) = sum(add, 1) / (size(DistFF,1) - 1);
end

If you provide some input values, e.g. created by rand, we could test the suggested improvements.

[EDITED] You have modified your question. The new code is not as meaningful as the original one was. The pre-allocation problem is still existing. C = A(i,:) + B(:,j).' consumes some time, because it creates a (I guess!) 1858x1858 matrix in each iteration, which needs 28 MB. Doing this 3.4 millions of times must need some time.

Prefer to post a minimal working example, which includes input data (e.g. created by rand) and which creates the wanted output data. The current example is too much simplified to see, what you want to achieve.

Do you have the parallel processing toolbox? Then parfor might be an elegant option.

Do you have a C compiler installed? Then a tiny C function might save some time.

[John D'Errico]

+1 of course. But there are too many unknown factors in here to really know how to resolve this issue. I would wonder if the transpose cannot just be done on the entire array in advance. DistPT does not seem to change in the code we see. So just transpose the entire mess once. That would be a significant time reduction alone.

So, is it a preallocation issue? Very very possibly so. But without seeing the entire code, there is no way to know.

Finally, it just looks like this is just computing the minimal difference between the rows of A, and the columns of DistPT. That begs other questions, like, should there be an absolute value in there? And why is it that when i==j, D(i,j) was assigned as zero in the original version? That would only be true if A was just DistPT.' already.

So there are many things I find confusing in this question. And also, why I personally would have avoided answering it in the first place. So I admire Jan's willingness to try.

[Habib]

Here is the code (most of it) of the function to really understand what I did.

I tried to transpose DistPT before beggining the loop, but it is still low.

@Jan, I used permute rather than .' because it is faster according to what I read in the internet.

I don't think it is because of the number of iteration because the line 'add(i,j) = DistFF(i,j)/D(i,j);' which is in the same loop takes only 1% of the running time, and the min function only 21%.

function [E_roadPT,DistRoadPT,perc_road,centr2,E_roadPTLinks,BikeCarCentr] = merging_layers_centrality_bikes(DistFF,DistRoad,DistPT,timeSteps,NearestStation,NearestStationX,StationsOnRange,Link,BikeMetroCentr,BikeStations,NearestBikeStation)
%merging_layers_centrality: Calculates the efficiency of the 2 layers that
%we want to merge, and the centrality of the metro stations

%let us combine the road and the metro to compute the station centrality
DistRoadPT = zeros(size(DistRoad,1),size(DistRoad,1),timeSteps);
centr = zeros(timeSteps,size(DistRoad,1));
perc_road = zeros(size(DistRoad,1),size(DistRoad,1),timeSteps);
BikeCarCentr = zeros(length(BikeStations),timeSteps); % centrality of bike stations
E_roadPT = zeros(size(DistFF,1),timeSteps);
DistPTpermuted = permute(DistPT,[2,1]);
for t = 1:timeSteps
    A = DistRoad(:,:,t);
    D = zeros(size(DistRoadPT,1),size(DistRoadPT,2));
    C = zeros(size(DistRoad,1),size(DistRoad,1));
    add = zeros(size(DistFF,1),size(DistFF,2));
    BCC = zeros(length(BikeStations),1); %bike car centrality
    for i = 1:size(DistRoad,1)
        for j = 1:size(DistRoad,1)
            if i == j
                D(i,j) = 0;
            else
                add2 = A(i,:) + DistPTpermuted(j,:); %the minimum index is the changing node + calculate the pourcentage distance road
                [v,p] = min(add2); %v is the value, p is the changing node
                D(i,j) = v;
                add(i,j) = DistFF(i,j)/D(i,j);
                if p == i || p == j
                    centr(t,p) = centr(t,p);
                else
                    centr(t,p) = centr(t,p)+1; %don't add +1 if it changes in the 1st or last node if p=i or j
                    if BikeMetroCentr(p,j) ~= 0 %bike metro centrality (which bike station to use if going from i to j by walking biking and metro)
                        BCC(NearestBikeStation(p)) = BCC(NearestBikeStation(p)) + 1;
                    end
                    if BikeMetroCentr(j,p) ~= 0
                        BCC(NearestBikeStation(j)) = BCC(NearestBikeStation(j)) + 1;
                    end
                end
                C(i,j) = A(i,p);
            end
        end
        E_roadPT(i,t) = (1/(size(DistFF,1)-1))*sum(add(i,:));
    end
    DistRoadPT(:,:,t) = D;
    perc_road(:,:,t) = C;
    BikeCarCentr(:,t) = BCC;
end
end

[Guillaume]

"I used permute rather than .' because it is faster according to what I read in the internet"

If there is a difference, it will be insignificant. Both do exactly the same rearranging of elements.

In your question title, you mention reshape. reshape is a lot faster than permute (and .') since it does not change the order of the elements (just the matrix header). However, you can't emulate .' with reshape.

[Habib]

I isolated the problem and it seems it is because of the selection of the raw of the DistPTpermuted matrix. Do you have an idea to avoid that time consumption?

[Adam]

I don't really have time to look in any detail, but can you not replace this:

A(i,:) + DistPTpermuted(j,:);

with just

A + DistPTpermuted;

outside that loop since, from my very quick scanning of that code they are the same size? Then just deal with the minority of cases where i == j separately.

May not be feasible and maybe not faster, but at a glance it does seem as though you are mostly just adding two arrays together in many small parts.

As an aside: Good to see someone using the profiler to focus in on the true sources of slow speed in their algorithm though. A lot of speed-related questions in Matlab Answers come from someone who is just guessing at what they think is the slow part!

[Jan]

Do not overestimate the results of the profiler. Matlab's JIT acceleration reorders and re-interpretes code to accelerate it. During the profiler runs, the JIT is disabled. This means, that the highlighted line must not be the the bottleneck of the code running without profiling. Compare the runtime of the code with enabled profiling and debugging with a run without these braking features.

This line is bewildering:

add2 = plus(Ai, + Dj);

The only effect of such tricks is that it can confuse the JIT. I hope Matlab ignores the unnecessary unitary +, but if not, this is a waste of time only. Prefer the simple:

add2 = Ai + Dj;

Obtaining a rowvector is time-consuming. Prefer to process arrays columnwise in Matlab, because this let it use the processor cache much more efficiently. So in your case using DistPTpermuted instead of DistPT is counter-productive. Try:

add2 = Ai + DistPT(:,j).';

as I have suggested in my answer already.

By the way, reshape can replace the tranpose operation in this case, when the argument is a vector. For vectors transpoasing and reshaping is the same operation.

If you really find the opinion in the internet, that the powerful permute is faster than the simpler transposition, this demonstrates, that you cannot believe everything you find there.

I asked already, if you have a C-compiler installed and are able to write small C-Mex functions. I could post the needed code also, but as long as you do not clarify if you can use it, I hestitate to spend the time for writing a C-mex funtion.

This line is just a waste of time:

centr(t,p) = centr(t,p);

[Habib]

Actually I tried all this stuff and the problem is stil in selecting this raw in the DistPTpermuted matrix.

@Adam: I can't add both matrices because I need to add only specific raws which are not essencially in the same position.

@Jan: Thanks for the:

centr(t,p) = centr(t,p);

which is clearly useless.

But the:

add2 = Ai + DistPT(:,j).';

Isn't efficient as it does a transpose at each loop, and if I do it at the beginning, I don't need to transpose anymore, juste call the concerned raw (which takes too much time).

I fund it really weird that, for a computer, selecting a simple raw from a matrix is that time consuming, even if it has to do it 24x1800x1800 times (3 loops). I really think there is a solution to that.

Please help me :)

[Jan]

DistPT(:,j).' is efficient, because the input of the transposition is a vector. Then transposing is equivalent to reshaping, which does not modify the contents of the RAM. But of course tranposing the array before the loops is a good idea also, as long as you do not access row-vectors from the marix afterwards.