Fastest selection of a coordinate in 3d array that fulfills a condition?

I have a relatively simple problem: I need to generate random coordinates of an array, but they must fulfill a condition. I thought find+ind2sub would make this easy to sample from condition-positive sets - and it does - but i was surprised to find that these are relatively slow, and on the arrays i'm using (>1e7 elements) it can take up to seconds. Needing to run this thousands of times makes that prohibitive.

Turns out that a loop is ~10-100 times faster merely generating completely random coordinates/indices and testing them for the condition. It's even more disparately faster with ever larger arrays, and approaches similarity only under extreme conditions (<<.001% fulfilling condition) in smaller arrays.

I'm guessing that this is the overhead of find+ind2sub exhaustively searching and indexing the entire array. Is there possibly a vectorized solution that avoids the overhead, or at least overcomes the loop through internal optimizations?

[Matt J]

Is there possibly a vectorized solution that avoids the overhead, or at least overcomes the loop through internal optimizations?

Yes, possibly, but the task hasn't been described fully enough to say more than that.

[Jan]

@Carson Purnell: I agree to Matt J. Please provide a small data set for the input and the wanted output. Then mention the typical dimensions of the problem. Afterwards we can serach an efficient method.

[Carson Purnell]

A sufficient test case is using randi(10,500,500,200) to generate the array. The only output is a single randomly selected voxel within the array (as coordinates or a linear index) that satisfies a basic comparison, for this test case ==1 works.

Again, to make it clear: I do not need many random points from this single array. ind2sub would be fast at that. I need exactly one point, but quickly because i'm doing this for thousands of non-identical arrays.

Answer 1

[Jan]

If the condition is not rare in the array, a gun-shot method can be very fast:

T = randi(10, 500,500,200);
tic
doSearch = true;
n = numel(T);
index = randi([1, n]);
while T(index) ~= 1
    index = randi([1, n]);
end
toc

Elapsed time is 0.009403 seconds.

A comparison with an exhaustive search:

tic;  % Matt J's suggestion:
 logIndex=(T(:)==1); %logical index
 n=sum(logIndex);
 k=randi(n,1);
 z=find(logIndex,k,'first');
 result=z(end); %linear index
toc

Elapsed time is 0.159166 seconds.

If the data have a small range, this is working, but of course the less frequent the condition is true, the more likely is a huge run-time.

Answer 2

[Matt J]

Why not as follows? I don't know if it's as fast as you were hoping, but it definitely doesn't take seconds per iteration.

T=randi(10, 500,500,200); %hypothetical input

tic;
 logIndex=(T(:)==1); %logical index

 n=sum(logIndex);

 k=randi(n,1);

 z=find(logIndex,k,'first');

 result=z(end); %linear index
toc

Elapsed time is 0.055249 seconds.

[Matt J]

If you have the Parallel Computing Toolbox, there is substantial speed-up to be had from working on the GPU. This comparison used the GTX 1080 Ti:

T=randi(10, 500,500,200);
Tgpu=gpuArray(T);

timeit(@()doIt(T)) %on CPU

ans =

0.0828

gputimeit(@()doIt(Tgpu))  %on GPU

ans =

0.0126

function result=doIt(T)

     logIndex=(T(:)==1);
    
     n=sum(logIndex);
    
     k=randi(n,1);
    
     z=find(logIndex,k,'first');
    
     result=z(end);
     
end