CODEWITHSUNDEEP
matlabloopsfor-loop
Fred
Fred

Reputation: 233

How to get only the diagonal of the matrix from these loops?

I have the following code for a 8 dimensional empirical copula that creates a 8d matrix but I only need the diagonal of the matrix which is named EC in this code. Since this code is very slow, is there anyway that I can get "EC" without computing all of "ecop"?

function EC = ecopula8d(x)

[m n] = size(x);

y = sort(x);

for r=1:m
        for q=1:m
            for p=1:m
                for o=1:m
                    for l=1:m
                        for k=1:m
                            for j=1:m
                                for i=1:m
             ecop(i,j,k,l,o,p,q,r) = sum( (x(:,1)<=y(i,1)).*(x(:,2)<=y(j,2)).*(x(:,3)<=y(k,3)).*(x(:,4)<=y(l,4))...
             .*(x(:,5)<=y(o,5)).*(x(:,6)<=y(p,6)).*(x(:,7)<=y(q,7)).*(x(:,8)<=y(r,8)) )/(m+1);
                                end
                            end
                        end
                    end
                end
            end
        end
end

for i=1:m
EC(i,1)=ecop(i,i,i,i,i,i,i,i);
end

Upvotes: 1

Views: 74

Answers (3)

Divakar
Divakar

Reputation: 221514

You can employ a completely vectorized solution with bsxfun -

EC = squeeze(sum(all(bsxfun(@le,x,permute(y,[3 2 1])),2),1))/(m+1)

The magic here happens with the use of permute enabling us to go full throttle on vectorization.

Here's a benchmarking test to compare this approach and the other partially vectorized approach with bsxfun on runtime efficiency -

x = rand(2000,2000);
y = sort(x);
m = size(x,1);

%// Warm up tic/toc.
for k = 1:100000
    tic(); elapsed = toc();
end

disp('----------- With completely vectorized solution')
tic
EC1 = squeeze(sum(all(bsxfun(@le,x,permute(y,[3 2 1])),2),1))/(m+1);
toc, clear EC1

disp('----------- With partial vectorized solution')
tic
for i = 1:m
    EC2(i,1) = sum(all(bsxfun(@le, x, y(i,:)), 2), 1)/(m+1);
end
toc

The runtimes thus obtained were -

----------- With completely vectorized solution
Elapsed time is 2.883594 seconds.
----------- With partial vectorized solution
Elapsed time is 4.752508 seconds.

One can pre-allocate for the other partially vectorized approach -

EC2 = zeros(m,1);
for i = 1:m
    EC2(i,1) = sum(all(bsxfun(@le, x, y(i,:)), 2), 1)/(m+1);
end

The runtimes thus obtained weren't that different though -

----------- With completely vectorized solution
Elapsed time is 2.963835 seconds.
----------- With partial vectorized solution
Elapsed time is 4.620455 seconds.

Upvotes: 1

knedlsepp
knedlsepp

Reputation: 6084

I haven't checked if your initial computation is correct (as in compared your implementation with the wikipedia article's formula), but your code should be equivalent to:

[m n] = size(x);
y = sort(x);
for i = 1:m
    EC(i,1) = sum(all(bsxfun(@le, x, y(i,:)), 2), 1)/(m+1);
end

Upvotes: 3

madbitloman
madbitloman

Reputation: 826

Once of the approaches I would use is to convert N-D array into square 2-D(if possible) and then simply extract diagonal term as they should be equal in both cases:

EC=diag(reshape(ecop,size1,size2));

I would suggest to try use Python because numpy has really nice and efficient linear algebra package to deal with N-D arrays. Matlab is pretty slow in adding and updating its libraries.

Upvotes: 0

Related Questions