How to effectively run loops and save time in computation? I have a matrix of size 'm' and run five loops from 1 to m. The logic is explained below. How to optimize the code and save time of calculation.

I need to find projection statistics for my matrix of size 'm' by initially calculating the sum of all possible observations (using logic as in code-Rousseuuw and Croux method). This needs a quite a large number of iterations. When my matrix size is 1000x1000, it takes close to 250seconds for this code to complete. Can someone help me in optimizing the code below thereby reducing the operation time.

    m=1000;
    H=rand(m);
    rng default;

    x=zeros(m,1);
    y=zeros(m,1);
    output1 = zeros(m,1);
    output = zeros(m,1);
    PS = zeros(m,1);

    for k=1:m
        for i=1:m
            for j=1:m
                if j~=i
                    x(j)=abs(H(i,k)+H(j,k));
                end
            end
            mask=x~=0; % mark the non-zero elements
            x=x(mask); % keep the non-zero elements
            y(i)=median(x);
            x=zeros(m,1); %Clearing x after calculating median and proceed to next iter
        end
        y(y==0)=NaN;
        output1(k)=1.1926*nanmedian(y);
    end
    for k=1:m
        for i=1:m
            output(i)=abs(H(k,i))/output1(i);
        end
        PS(k,1)=max(output);
    end

Even after making the suggestions in comments, the computation time is still high. Without pre-initialising, - 250s With pre-initialising - 250s After changing namedian to median - 200s (recent)

[Guillaume]

What is the point of

x(x==0)=NaN;
y(i)=nanmedian(x);

When the previous line has guaranteed that there are no 0 left in x and thus will never be any nan replacement.

[SanthoshKumar C]

@Guillaume, Yes. Thats the commented part of code.. Edited now.

[Guillaume]

Also I'm fairly certain that your code contains a bug because of the if j~=i because you never erase x. So when j == i you reuse the x(j) value calculated in a previous iteration. For example at the last step of the loop, the first m-1 elements of x are those calculated for j=1:m-1 and i = m, but the last element is the one calculated at the previous i step, hence j = m and i = m-1

%at step k = m, j=m, x is:
x = [abs(H(1:m-1, m) + H(m, m)); abs(H(m-1, m) + H(m, m))].'

If you had

 for k=1:m
    for i=1:m
       x = zeros(1, m);
       for j=1:m
           %...

You'd get a very different result

[SanthoshKumar C]

Thanks @Guillaume for noting it.. Now i cleared x.. Updated code in question.

Answer 1

[Andrei Bobrov]

a = abs(H + permute(H,[3,2,1])).*permute(diag(nan(m,1))+ones(m),[1,3,2]);
a(a == 0) = nan;
b = 1.1926*median(median(a,3,'omitnan'));
PS = max(abs(H./b(:)'),[],2);

other variant

m = 1000;
H = rand(m);
b = zeros(m);
for k = 1:m
    a = abs(H(:,k) + H(:,k).');
    % a = abs(bsxfun(@plus,H(:,k),H(:,k).')); - for MATLAB <= R2016a
    a(1:m+1:end) = nan;
    a(a == 0) = nan;
    b(:,k) = median(a,2,'omitnan');
end
c = 1.1926*median(b);
PS = max(abs(H./c(:)'),[],2);

[SanthoshKumar C]

Dear Andrei,

I tried your suggestion and got the following error.. Seems I am missing something.. Can you help..

    m=1000;
    H=rand(m);
    rng default;

    x=zeros(m,1);
    y=zeros(m,1);
    output1 = zeros(m,1);
    output = zeros(m,1);
    PS = zeros(m,1);

    y = abs(H + permute(H,[3,2,1])).*permute(diag(nan(m,1))+ones(m),[1,3,2]);
    y(y == 0) = nan;
    output1 = 1.1926*median(median(y,3,'omitnan'));
    PS = max(abs(H./b(:)'),[],2);

Error using + Matrix dimensions must agree.

>> Error in Test (line 33) y = abs(H + permute(H,[3,2,1])).*permute(diag(nan(m,1))+ones(m),[1,3,2]); >>

[Andrei Bobrov]

>> m = 3;
H = rand(m);

a = abs(H + permute(H,[3,2,1])).*permute(diag(nan(m,1))+ones(m),[1,3,2]);
a(a == 0) = nan;
b = 1.1926*median(median(a,3,'omitnan'));
PS = max(abs(H./b(:)'),[],2)
PS =
    1.0193
    0.6041
    0.5336
>>

[Andrei Bobrov]

@SanthoshKumar C: See my answer - "other variant"

[SanthoshKumar C]

Dear Andrei Bobrov,

In edit-1: (a = abs(H + permute(H,[3,2,1])).*permute(diag(nan(m,1))+ones(m),[1,3,2]); )

H is a matrix of size 1000x1000 whereas permute(H,[3,2,1])) is a mtrix of size 1x1000x1000. Therfore it gives error of * _ Error using + Matrix dimensions must agree._*

In other variant: H(:,k) + H(:,k).' are of different matrix dimensions (1000x1) and (1x1000) respectively... It also gives the same error as above.. Also , taking transpose of it wont work as it doesnt satisfy (abs(H(i,k)+H(j,k))

[Andrei Bobrov]

You use old version of MATLAB.

Please use "other variant".

[SanthoshKumar C]

Yes. I am using Matlab 2016a. Even 'other variant' is giving the error as follows:

Error using sum Dimension argument must be a positive integer scalar within indexing range.

Error in Untitled2 (line 14) a = abs(bsxfun(@sum,H(:,k),H(:,k).'));

[SanthoshKumar C]

@Andrei Bobrov,

I installed MATLAB 2017a today and tested the Other Variant solution given by you. But it still shows the following error:

Error using sum
Dimension argument must be a positive integer scalar within
indexing range.

Error in Untitled2 (line 76)
      a = abs(bsxfun(@sum,H(:,k),H(:,k).'));

Can you please help?

[Andrei Bobrov]

Excuse me! My typo! I'm fixed my "other variant". Please use his.

[SanthoshKumar C]

@Andrei Bobrov.. After changing @sum to @plus, it worked... Thanks a lot

Answer 2

[Jos (10584)]

A few observations:

1. you should pre-allocate x, y and output before the loops
2. H(i,k) + H(j,k) equals H(j,k)+H(i,k), so you can half the computation time by having j running from i+1 to m.
3. this will also make the "if i~=j" redundant!
4. after removing all zeros from x using mask, the statement x(x==0) = NaN not needed!
5. therefore you can use the faster median rather than nanmedian

[SanthoshKumar C]

In snapshot, Without pre-initialising, - 251s With pre-initialising - 250s After changing namedian to median - 200s

But i cannot use your suggestion 2 and 3, as for calulating median y, I need all the data points of x from 1:m.. Please let me know if you feel otherwise.

[SanthoshKumar C]

This is the modified code after the suggestions 1, 4 and 5.

    m=1000;
    H=rand(m);
    rng default;

    x=zeros(m,1);
    y=zeros(m,1);
    output1 = zeros(m,1);
    output = zeros(m,1);
    PS = zeros(m,1);

    for k=1:m
        for i=1:m
            for j=1:m
                if j~=i
                    x(j)=abs(H(i,k)+H(j,k));
                end
            end
            mask=x~=0; % mark the non-zero elements
            x=x(mask); % keep the non-zero elements
            y(i)=median(x);
            x=zeros(m,1); % Clearing x after calculating median and proceed to next iter
        end
        y(y==0)=NaN;
        output1(k)=1.1926*nanmedian(y);
    end
    for k=1:m
        for i=1:m
            output(i)=abs(H(k,i))/output1(i);
        end
        PS(k,1)=max(output);
    end