Why are relational operators so slow in this case?

I have a 2305-by-4-by-4455 array, "P", from which I want to generate a logical vector, "isdis", for indexing as follows:

isdis = (P(:,2,:) >= P(:,3,:) & P(:,3,:) >= P(:,4,:) & P(:,4,:) >= 0) &... 
        (P(:,2,:) <= (sqrt(2)-1)*P(:,1,:)) &...
        (P(:,2,:)+P(:,3,:)+P(:,4,:) <= P(:,1,:));

Unfortunately, this takes a relatively long time to evaluate for some reason. On my machine (Windows 8, 4GB memory, Core i7 processor, Matlab R2012b), this statement takes about 0.8 sec to evaluate. Which may not seem like a long time, but it accounts for 55% (2 sec) of the larger function that it is a part of, which has to be called many times in my application (at 2 sec a pop, that's an expensive function call). I would have thought that relational operators like this (especially when vectorized as it is) would cost virtually nothing, can anyone help me understand why this is so slow?

If it helps, the one multiplication of:

 (sqrt(2)-1)*P(:,1,:))

takes 0.07 sec. And the one addition of:

 P(:,2,:)+P(:,3,:)+P(:,4,:)

takes 0.19 sec.

Any suggestions? This line needs to evaluate in at least an order of magnitude less time, and it seems like it ought to be possible to do it in at least 2 orders of magnitude less time.

[Teja Muppirala]

The suggestions from others below will help reduce the time a little from 0.8 to something like maybe 0.5~0.6 or somewhere in that ballpark, but I don't think you're going to do much better than that. Certainly not one or two orders of magnitute better. The matrices you are dealing with are large, 330MB in double precision. Sometimes a computation just inherently takes time.

Consider, even simply this:

tic
P > 0;
P & P;
P * sqrt(2);
toc;

already takes 1/3 of a second on my system (which also runs your code in 0.8 seconds).

So given that, I don't see how your calculation, which involves more work than this, could go much faster.

[Oliver]

Thanks for the comment, I guess I just assumed that relational operators would be cheaper than they really are.

Answer 1

[Jan]

Creating the temporary array P(:,i,:) is time consuming, because the memory must be copied from non-contiguous blocks. After applying permute(P, [1,3,2]) the comparison is much faster.

P(:,2,:)+P(:,3,:)+P(:,4,:) is 30% slower than sum(P, 2) - P(:,1,:) for the same reason.

[Oliver]

Thanks so much for your help, Jan. I had some clarifying questions:

(1) Are you saying that for each logical comparison Matlab creates a copy of the relevant portion of P that it will use for the comparison? So for instance if I wanted P(:,1:) < P(:,2,:), matlab would essentially create a = P(:,1,:) and b = P(:,2,:), and then perform the comparison between a < b? That seems really inefficient, I thought that indexing into an array would be faster than explicitly creating slices as suggested by Cedric. Using Cedric's suggestion I see a 2x speedup, so it must be helping, though.

(2) Could you explain what you were suggesting with the permute(P,[1,3,2]) command? I didn't quite understand what you meant. Are 2D matrices stored contiguously, while 3D arrays (and higher-D) are split up?

Thanks again for the help!

[James Tursa]

(1) Yes.

(2) P(:,2,1) is not contiguous with P(:,2,2) etc. In general, a bunch of non-contiguous memory has to be copied to form P(:,2,:) and other similar expressions. By permuting the array first, you get all of the blocks you are interested in into contiguous memory, potentially making subsequent indexing faster (less memory run through the high level cache for each copy).

[Matlab2010]

Can someone please clarify where the ordering of the permute function comes from?

[Jan]

The leasing dimensions of a multi-dim array have contiguous memory. Therefore the permute operation swaps the 2nd and 3rd dimension:

P = rand(2305, 4, 4455);
Q = permute(P, [1,3,2]);
size(Q)   % 2305, 4455, 4

Does this clear the ordering?

Answer 2

[Cedric]

The first thing that comes to my mind is that you could store slices in 2D variables to eliminate repetitive similar block-indexing operations in the expression for isdis:

 slice1 = P(:,1,:) ;
 slice2 = P(:,2,:) ;
 slice3 = P(:,3,:) ;
 slice4 = P(:,4,:) ;

You would then have something like (to check on your side):

 isdis = (slice2 >= slice3 & slice3 >= slice4 & slice4 >= 0) &... 
         (slice2 <= (sqrt(2)-1)*slice1) &...
         (slice2+slice3+slice4 <= slice1);

You won't get an order of magnitude improvement in computation time with only that though (maybe 20-30%).

[Jan]

This has the same advantage as using PERMUTE for a reordering of the dimensions. Using these copies avoids the need of creating temporary arrays Q(:,:,i) repeatedly, therefore I assume, this is faster. +1

[Cedric]

Actually your solution is the fastest; I tried to beat it using some hybrid approach, but I can't even get close to your permute().