How to Vectorize Dependent For-Loops

Question

I'm working on a function that takes a 1xn vector x as input and returns a nxn matrix L.
I'd like to speed things up by vectorizing the loops, but there's a catch that puzzles me: loop index b depends on loop index a. Any help would be appreciated.

x = x(:); 
n = length(x);
L = zeros(n, n);
for a = 1 : n,
    for b = 1 : a-1,
        c = b+1 : a-1;
        if all(x(c)' < x(b) + (x(a) - x(b)) * ((b - c)/(b-a))),
            L(a,b) = 1;
        end
    end
end

Answer 1

From a quick test, it looks like you are doing something with the lower triangle only. You might be able to vectorize using ugly tricks like ind2sub and arrayfun similar to this

tril_lin_idx = find(tril(ones(n), -1));
[A, B] = ind2sub([n,n], tril_lin_idx);
C = arrayfun(@(a,b) b+1 : a-1, A, B, 'uniformoutput', false); %cell array
f = @(a,b,c) all(x(c{:})' < x(b) + (x(a) - x(b)) * ((b - c{:})/(b-a)));
L = zeros(n, n);
L(tril_lin_idx) = arrayfun(f, A, B, C);

I cannot test it, since I do not have x and I don't know the expected result. I normally like vectorized solutions, but this is maybe pushing it a bit too much :). I would stick to your explicit for-loop, which might be much clearer and which Matlab's JIT should be able to speed up easily. You could replace the if with L(a,b) = all(...).

Edit1

Updated version, to prevent wasting ~ n^3 space on C:

tril_lin_idx = find(tril(ones(n), -1));
[A, B] = ind2sub([n,n], tril_lin_idx);
c = @(a,b) b+1 : a-1;
f = @(a,b) all(x(c(a, b))' < x(b) + (x(a) - x(b)) * ((b - c(a, b))/(b-a)));
L = zeros(n, n);
L(tril_lin_idx) = arrayfun(f, A, B);

Edit2

Slight variant, which does not use ind2sub and which should be more easy to modify in case b would depend in a more complex way on a. I inlined c for speed, it seems that especially calling the function handles is expensive.

[A,B] = ndgrid(1:n);
v = B<A; % which elements to evaluate
f = @(a,b) all(x(b+1:a-1)' < x(b) + (x(a) - x(b)) * ((b - (b+1:a-1))/(b-a)));
L = false(n);
L(v) = arrayfun(f, A(v), B(v));

Answer 2

If I understand your problem correctly, L(a, b) == 1 if for any c with a < c < b, (c, x(c)) is “below” the line connecting (a, x(a)) and (b, x(b)), right?

It is not a vectorization, but I found the other approach. Rather than comparing all c with a < c < b for each new b, I saved the maximum slope from a to c in (a, b), and used it for (a, b + 1). (I tried with only one direction, but I think that using both directions is also possible.)

x = x(:);
n = length(x);
L = zeros(n);

for a = 1:(n - 1)
  L(a, a + 1) = 1;
  maxSlope = x(a + 1) - x(a);
  for b = (a + 2):n
    currSlope = (x(b) - x(a)) / (b - a);
    if currSlope > maxSlope
      maxSlope = currSlope;
      L(a, b) = 1;
    end
  end
end

I don't know your data, but with some random data, the result is the same with original code (with transpose).

Answer 3

An esoteric answer: You could do the calculations for every a,b,c from 1:n, exclude the don't cares, and then do the all along the c dimension.

[a, b, c] = ndgrid(1:n, 1:n, 1:n);

La = x(c)' < x(b) + (x(a) - x(b)) .* ((b - c)./(b-a));
La(b >= a | c <= b | c >= a) = true;

L = all(La, 3);

Though the jit would probably do just fine with the for loops since they do very little.

Edit: still uses all of the memory, but with less maths

[A, B, C] = ndgrid(1:n, 1:n, 1:n);

valid = B < A & C > B & C < A;
a = A(valid); b = B(valid); c = C(valid);

La = true(size(A));
La(valid) = x(c)' < x(b) + (x(a) - x(b)) .* ((b - c)./(b-a));
L = all(La, 3);

Edit2: alternate last line to add the clause that c of no elements is true

L = all(La,3) | ~any(valid,3);