Is it possible to vectorize this?

A = [starts ends];
b = zeros(size(A,1),1);
c = zeros(size(A,1),1);

for i = 1:1:size(A,1)
    b(i) = max(x(A(i,1):A(i,2)));
    c(i) = min(x(A(i,1):A(i,2)));
end

starts, ends and x are column vectors. The intervals does not overlap.

[madhan ravi]

whats your desired output of b and c? give an example

[Jan]

By the way, length(A) will fail, if starts and ends are scalars and A is a [1 x 2] matrix. Prefer size(A,1) instead to avoid such ambiguities in general.

[Joakim Hansen]

Yeah, it's a mistake I did when I switched out all the real variable names.

[Joakim Hansen]

x(A(i,1):A(i,2)) should output all the points of x with indices between A(i,1) and A(i,2). For example A(1,1)=1 and A(1,2)=5 will be [1 2 3 4 5]. So x([1 2 3 4 5]) = [x(1) x(2) x(3) x(4) x(5)]. Then I wanna find the largest and smallest number of those. So b will contain all the max values for each i, and c all the min values for i.

[Jan]

I guess that you want to vectorize this to increase the speed. Then it matters, how large A and x are and maybe if the intervals overlap. It is important if you run this piece of code on GB of data once only, of if you call it millions of times with small input arrays. So please provide some relevant input data, e.g. created by rand.

Do you have a C-compiler? Then a C-mex function might be useful.

[Joakim Hansen]

In my data sets A is 79x2 and x is 775000x1 and no overlap. I do not have a C-compiler.

[Jan]

If this is an important and time-critical part of your code, I recommend to install a C-compiler and use a C-Mex function. Did you try parfor?

[Joakim Hansen]

It's just a part of some coursework where vectorization is preferred over loops, and I wasn't able to vectorize that loop. If the only solution is some relatively advanced stuff like C-Mex, I will just go with the loop.

Answer 1

[Christine Tobler]

There is no general way to vectorize this. If the starts and ends vector define intervals that are not overlapping, and that cover every index in x, you can do the following:

ind = zeros(length(x), 1);
for i=1:size(A, 1)
    ind(starts(i):ends(i)) = i;
end
% ind(j) gives the number of the interval that j (and x(j)) are in.

b = accumarray(ind, x(:), [], @max);
c = accumarray(ind, x(:), [], @min);

Depending on where starts and ends are coming from, maybe you can construct ind in a vectorized format from there?

Answer 2

[Bruno Luong]

For generic x, A, no there is no direct way.

What you might do is decomposing the set of intervals in bigger set of non-overlapping intervals, find MIN and MAX in subintervals, which is faster then group them together depending how the original interval is composed of subintervals.

Not sure it's more efficient, but this conclusion might depends on how the overlapping/non-overlapping of the intervals are.

[Joakim Hansen]

None of the intervals are overlapping.

[Jan]

@Bruno: What do you think about sorting x at first? While this makes it trivial to find minimal and maximal values, locating the corresponding indices gets much more expensive.

[Bruno Luong]

No, I have not say anything about sorting x.

Sorting A break points yes.

[Bruno Luong]

'None of the intervals are overlapping.'

So I think the bottleneck is that MATLAB make a duplication of subdata for every iteration. In that case you might have to program MEX file to avoid such copying.

If you are running R2018, One intermediate alternative is using my on-table share data mex file as I show in this thread

Answer 3

[Guillaume]

First, note that you should never use length on a matrix. You're using length as synonymous to the number of rows but if A has more columns than rows then it's going to be the number of columns. So use size with the proper dimension index.

There is no straightforward way to vectorise your code. It can be done at the expense of a large temporary array and a fair bit of juggling, so I'm not convinced it's worthwile. Here is one way:

mask = zeros(numel(starts), numel(x) + 1);  %one more column than elements in x as we will be using ends+1 as index
mask(sub2ind(size(mask), 1:numel(starts), starts')) = 1;
mask(sub2ind(size(mask), 1:numel(ends), ends' + 1)) = -1;
mask = cumsum(mask(:, 1:end-1), 2);
mask(mask == 0) = NaN;
[c, b] = bounds(mask .* x', 2);  %requires R2017a. In previous versions use separate calls to min/max

This may be slightly faster but will output row vectors for c and b instead of column vectors:

mask = zeros(numel(x) + 1, numel(starts));
mask(sub2ind(size(mask), starts', 1:numel(starts))) = 1;
mask(sub2ind(size(mask), ends' + 1, 1:numel(ends))) = -1;
mask = cumsum(mask(:, 1:end-1));
mask(mask == 0) = NaN;
[c, b] = bounds(mask .* x);  %requires R2017a. In previous versions use separate calls to min/max

Answer 4

[Bruno Luong]

The orginal method is quite fast, for random intervals (they cover 50-60% of x) on my computer it takes 1.7 ms in average. Guillaume's method takes forever. I add my own method called Scanning here.

Orginal algo  : t=0.001669 [s]
Scanning algo : t=0.002096 [s]
Guillaume algo: t=1.077173 [s]
Christine algo: t=0.026156 [s]

Test code

ntest = 100;
t = zeros(3,ntest);
for j = 1:ntest

      x = rand(775000,1);
      n = 79;

      b = randperm(length(x),2*n);
      A = reshape(sort(b),2,[])';

      % Orginal method
      tic
      b = zeros(size(A,1),1);
      c = zeros(size(A,1),1);
      for i = 1:1:size(A,1)
          b(i) = max(x(A(i,1):A(i,2)));
          c(i) = min(x(A(i,1):A(i,2)));
      end
      t(1,j) = toc;

      % Scanning method, assuming Interval are sorted
      tic
      b = zeros(n,1);
      c = zeros(n,2);
      for k=1:n
          is = A(k,1);
          minx = x(is);
          maxx = minx;
          for i=is+1:A(k,2)
              if x(i) < minx
                  minx = x(i);
              elseif x(i) > maxx
                  maxx = x(i);
              end
          end
          b(k) = maxx;
          c(k) = minx;
      end
      t(2,j) = toc;

      % Guillaume method
      tic;
      m = size(x,1);
      mask = zeros(n, m + 1);  %one more column than elements in x as we will be using ends+1 as index
      mask(sub2ind(size(mask), 1:n, A(:,1)')) = 1;
      mask(sub2ind(size(mask), 1:n, A(:,2)' + 1)) = -1;
      mask = cumsum(mask(:, 1:end-1), 2);
      mask(mask == 0) = NaN;
      [c, b] = bounds(mask .* x', 2);
      t(3,j) = toc;
      clear mask

      % Modification of Christine Method
      tic
      ind = zeros(length(x), 1);
      for i=1:size(A, 1)
          ind(A(i,1):A(i,2)) = i;
      end
      % ind(j) gives the number of the interval that j (and x(j)) are in.
      tf = ind > 0;
      x = x(:);
      b = accumarray(ind(tf), x(tf), [], @max);
      c = accumarray(ind(tf), x(tf), [], @min);
      t(4,j) = toc;

end
t = mean(t,2);
fprintf('Orginal algo  : t=%f [s]\n', t(1));
fprintf('Scanning algo : t=%f [s]\n', t(2));
fprintf('Guillaume algo: t=%f [s]\n', t(3));
fprintf('Christine algo: t=%f [s]\n', t(4));

[Guillaume]

Guillaume's method takes forever

I'm not surprised by that. It's got to fill a 775001 x 79 matrix and then do a cumsum on that.