Logical indexing two dimensions. How do I avoid a nested for loop?

I have a matrix X with doubles of size = 450 8156 and another matrix C with logical values with size = 64 8156 where 400 of the 8156 elements of each of the 64 rows are true. And I also have yet another matrix T of logical values with size = 9 450 where 400 of the 450 elements of each row are true.

I want to extract all the true values of X so that the resulting matrix becomes of size = 400 400 9 64

I can do this easily with a double for loop like this:

for n = 1 : 64
    for z =  1 : 9
    A(:,:,z,n) = X(T(z,:)',C(n,:));
    end
end

but it is very slow. Is there a faster, more reasonable, vectorized way of doing thing? By using repmat, reshape, bsxfun, arrayfun etc.? Anyone good with logical indexing that can help me out?

Thanks.

[Jos (10584)]

Did you pre-allocate A? Doing so, should make the for-loops pretty fast.

[Peta]

Yes I did. It’s not incredibly slow, but it takes almost 2 seconds for me to run which is 1.99 seconds more than I would prefer. And since it’s based on a nested for loop my intuition tells me that there should be a much faster way of doing it in one call. But I’m too stupid to see it.

[Stephen23]

Your output array has 92160000 elements (92 million is nearly 1GB of double class) and it takes two seconds to process...

Generating intermediate variables is likely to take just as much time as the loops with a preallocated array, if not longer. Consider your time investment: are these loops really a bottleneck? Will the more complex code maintenance be worth it?

[Peta]

Yes, it is really a bottleneck. And the data isn’t always 1:64 and 1:9, those dimensions change over time so the problem is usually smaller than what I wrote in my example and it’s frustrating that it should slow down my code by several seconds per iteration since I’m not even calculating anything here, just rearranging existing data.

[Stephen23]

"just rearranging existing data"

Actually it is copying existing data... almost 1 GB of it.

Answer 1

[Jan]

You can try a logical indexing:

T = T.' ;
A = zeros(400,400,9,64);
Cn = false(1, size(X, 2));
for n = 1 : 64
  Cn(:)       = false;
  Cn(C(n, :)) = true;
  for z =  1 : 9
    A(:,:,z,n) = X(T(:,z), Cn);
  end
end

I have no idea if this is faster. Please provide some meaningful test data, such that the readers can perform some experiments before posting.

Do you have the parallel programming toolbox?

[Peta]

Thank you that did make it faster and now it only takes about 0.5 seconds for me to run.

And yes I have the parallel programming (computing??) toolbox, is there any functionality in there than could be useful here? I tried making the inner loop in the code you wrote into a parfor instead of for, but that makes it much much slower.

Code for generating the sample data variables can be created like this:

X = (rand(450,8156));
C = false(64,8156);
for i = 1 : 64
    while sum(C(i,:),2) ~= 400
    C(i,randi([1 8156])) = true;
    end
end
T = false(9,450);
for i = 1 : 9
    while sum(T(i,:),2) ~= 400
    T(i,randi([1 450])) = true;
    end
end

[Jan]

@Peta: You can use parfor instead of for.

Answer 2

[Azzi Abdelmalek]

You can improve your for loop by pre-allocating

A=zeros(400,400,9,64)

Answer 3

[Jos (10584)]

Pre-allocate but also put the transpose out of the loop!

T = T.' ;
A = zeros(400,400,9,64)
for n = 1 : 64
    for z =  1 : 9
    A(:,:,z,n) = X(T(:,z), C(n,:));
    end
end

Damn, this is fast :-)

[Peta]

Of course doing one isn’t a problem and appears fast enough. But I need to do this tens of thousands of times and then pretty soon the slowness of this code adds up to seconds, minutes, hours.. Which is why I’m searching for a more sophisticated method that doesn’t involve nested for loops. And I see from everyone’s answer here that I should have made it clear when formulating my question that I am familiar with the basic concept of pre-allocation :)

[Jos (10584)]

tens of thousands of times x 0.1s per time (which is fast) = thousands of seconds = a few hours. Maybe Matlab is not the right tool for the job?