MATLAB: Union of closed intervals

arrayintervalMATLAB

I have a matrix A of size 2 x N matrix representing N closed intervals. That is, the i-th interval is [ A(1, i) .. A(2, i)]. I would like to construct a matrix C of size 2 x M with the set of disjoint intervals of minimum size whose union is exactly the union of the intervals in A.
I was able to write the code using a loop. Specifically, I sorted the intervals in ascending order according to their lower bound, and built the matrix C with a greedy algorithm which scanned the sorted intervals one by one.
Is there a more efficient way of computing C, preferrably without any loops?
For completeness, here is my code:
function union = unify_intervals(intervals)
% Computes a set of disjoint intervals representing the union of the given
% set of intervals. Usage:
%   union = unify_intervals(intervals)
%

% Input arguments
%   intervals
%       2 by N matrix containing the input set of intervals. Each column
%       specifies the interval bounds.
%
% Output arguments
%   union
%       A 2 by N matrix containing a set of disjoint intervals of minimum
%       size, whose union is equal to the union of the given input 
%       intervals.
      %%sort intervals by starting point
      [~, order] = sort(intervals(1, :), 2);
      sorted_intervals = intervals(:, order);
      union = sorted_intervals(:, 1);
      for i = 2:size(sorted_intervals, 2)
          union_head = union(:, end);
          interval = sorted_intervals(:, i);
          overlap = max(interval(1), union_head(1)) <= min(interval(2), union_head(2));
          if overlap
              union(2, end) = max(union_head(2), interval(2)); % update union head
          else
              union = [union, interval]; % append new union head
          end
      end
  end

Best Answer

Start with a pre-allocation of the output instead of the iterative growing. Then update Head on demand only.
% sort intervals by starting point
[~, index] = sort(intervals(1, :), 2);
s          = intervals(:, index);
u    = zeros(size(s));
uk   = 1;
Head = s(:, 1);
for k = 2:size(s, 2)
   if s(1, k) < Head(2)    % Overlap: Update current interval
      if s(2,k) > Head(2)  % Faster than: Head(2)=max(Head(2), s(2,k))
         Head(2) = s(2, k);
      end
   else                    % New interval:
      u(:, uk) = Head;     % Append current interval
      uk       = uk + 1;
      Head     = s(:, k);  % Store new interval
   end
end
u(:, uk) = Head;           % Append last interval:
u        = u(:, 1:uk);     % Crop unused memory
This is 60% faster for the input containing 100'000 elements and 1'000 intervals.
Related Question