Technical Articles and Newsletters

# An Adventure of Sorts–Behind the Scenes of a MATLAB Upgrade

By Bill McKeeman, MathWorks and Loren Shure, MathWorks

While choosing how to upgrade sort for MATLAB 7, we used MATLAB to prototype `sort` algorithms and to add the newsort direction qualifiers `'ascend'` and `'descend'`. We translated the prototype from MATLAB into C++ to add this latest version of sort to Release 14—and found a few surprises in the process.

Algorithms in MathWorks products and user applications use the MATLAB `sort`

There are 12 types of sortable matrices in MATLAB: the eight integer types `([u]int8/16/32/64);` plus single, double, char, and logical. Single and double may be real or complex; double and logical can be sparse. Full n-dimensional (ND) arrays can be sorted along any selected dimension, in either ascending or descending order. `sort` can also return the permutation vector relating the input values to the outputs, with call syntax ranging from the most simple

```r = sort(v);
```

to the most complete

```[r,i] = sort(v, dim, dirflag);
```

The most readily available upgrade would use qsort from the C library, but qsort was too specific for the sort tasks that MATLAB requires. The C library qsort has a hard-coded stride of 1 (data is assumed to be arranged in contiguous memory). MATLAB also needs to sort data containing NaN values, which qsort does not accommodate.

The available C qsort implementations share a common characteristic: to achieve type independence, data access is indirect (via void* pointers) and the user supplies a comparison function. As a result, general-purpose sort algorithms written in C are inelegant and somewhat slow.

### Selecting an Algorithm

Sorting experts recommend a three-level algorithm, insert/quick/bin , where each level calls the one to its left if the size of the array being sorted passes below some experimentally discovered cutoff. When we updated `sort` , we expected to require four levels ( insert/quick/bin/quick ). To achieve speed proportional to the size of the data (order N ) bin sort requires as much auxiliary storage for the bins as it does for the input data. Because quicksort manipulates the data in place, it can sort larger arrays, albeit somewhat more slowly. (Hoare’s 1962 quick sort algorithm has undergone a generation of changes, particularly to improve performance when the data is not random.)

Resources
Bentley, J. L. and M. D. McIlroy, 1993. Engineering a sort function. Software—Practice and Experience 23, 11, 1249--1265.
Knuth, D. E. 1998. The Art of Computer Programming, Volume 3: Sorting and Searching, 2nd ed. Boston: Addison-Wesley.
Cormen et al. 2001. Introduction to Algorithms, 2nd ed. Columbus: McGraw-Hill.

### Prototyping Alternative sort Algorithms

We explored speed and memory tradeoffs by prototyping the basic algorithms in M. The prototypes work only on vectors. (Full ND arrays are discussed in Extending Vector Sorts to N-Dimensions later in this article.)

`````` function v = viSort(v)    % vector insert
for i=1:numel(v)-1
[m, k] = min(v(i:end));
v(i+1:i+k-1) = v(i:i+k-2)  % make room
v(i) = m;                  % place min
end```
```
`````` function v = vqSort(v)    % vector quick sort
if numel(v) > 1
piv = median(v);           % divide
lss = v(v < piv);          % and
eql = v(v == piv);         % conquer
gtr = v(v > piv);
v = [vqSort(lss), eql, vqSort(gtr)];
end```
```
`````` function v = vbSort(v)          % vector bin sort
n = numel(v);
if n > 1
mn = min(v);
mx = max(v);
binAvg = 200;                    % experimental
binCt = ceil(n/binAvg);
slope = binCt/(mx-mn);
binNo = floor(slope.*(v-mn))+1;  % assign bin
bins = cell(1, binCt+1);         % allocate bins
for i=1:binCt+1
bins{i} = vqSort(v(binNo==i));
end
v = [bins{:}];                  % reassemble
end```
```

As we planned to have `vqSort` call `viSort` when there were few enough data points, we modified the quick sort. design.

`````` function v = vqSort(v)  % quick sort
if numel(v) < 32
v = viSort(v);           % experimental
else
piv = median(v);         % need approximation
lss = v(v < piv);
eql = v(v == piv);
gtr = v(v > piv);
v = [vqSort(lss), eql, vqSort(gtr)];
end```
```

We translated the M-code vector sorts into C in-place sorts so that we could test our four-level design hypothesis. The results surprised us. For vectors of random data on Intel x86 hardware, the most modern quick sorts were slower than the already-implemented quick sort in MATLAB! Algorithmic changes made to the modern sorts to ensure good behavior on less-than-random input were costly, while the existing MATLAB sort cleverly uses its pivot-picking approximation to do some sorting work.

We finally decided that best performance on random data was more useful than adequate performance on all data, and kept the existing MATLAB quick sort. (Performance results depend somewhat on the underlying hardware.)

We took into account that the C implementation of bin sort needs three passes: one pass to decide on the number of bins, one to bin the data, and one to call qsort on each bin. The results are left in place, conveniently ordered for output. We can justify this considerable overhead because it results in an order N, instead of NlogN, run-time algorithm.

Another surprise! We experimented to find a cutover from quick sort to bin sort, but never found it. Our best effort at a C-implemented bin sort was slower than our C-implemented quick sort, even on datasets containing up to 10 million elements. As a result, we abbreviated our development plan to a two-stage sort: insert/quick. Both these algorithms already have satisfactory implementations in MATLAB.

### Extending Double Sorts to New Types: Plus `NaN` and `Inf`

Next, we extended the double sorts to the sortable MATLAB types. We recoded the `insert/quick` sort pair as C++ template functions parameterized by the underlying MATLAB types.

MATLAB sorts `NaN` values high, beyond all other single/double values. Because `NaN` and `Inf` need to be dealt with only for the real types, we implemented a function (called a jacket) to handle the cases for single and double. The jacket calls the usual sort routine on the finite values and handles the nonfinite values correctly for MATLAB. This jacket’s sorting prepass counts the `NaN`s and replaces them with `+Inf`. After the sort, the jacket overwrites an appropriate number of `+Inf` values with `NaN`s at the top end of the sorted vector. (Had we still been doing bin sort, this jacket would also have had to count and remove `Inf` values, as the bin selection algorithm depends on finite values.)

`````` function v = vSort(v)             % sort NaN jacket
if isfloat(v)
nans = isnan(v);
v(nans) = inf;
v = vqSort(v);
v(end-sum(nans)+1:end) = nan;
else
v = vqSort(v);
end```
```

### Extending Vector Sorts to N Dimensions

We extended the vector sorts to sort along a user-selected dimension of an ND array. Array data is laid out in memory in column-major order: elements in a column are adjacent in memory. Elements along another dimension are separated by a constant distance (the so-called stride, usually larger than 1). We prototyped in M first

`````` function v = strideSort(v, d)       % stride jacket
s = size(v);
stride = prod(s(1:d-1));
for i = 1 : stride*s(d) : numel(v)   % #of elements
for j = i : i+stride-1
idx = j : stride : j+stride*s(d)-1;
v(idx) = vSort(v(idx));
end
end```
```

We recoded the C++ in-place insert/quick sort pair to step by a stride, and then wrote a second level of C++ jacket to call them.

### Extending Array Sorts to Permutation Vectors

We extended the array sorts to provide the permutation vectors, which required that the sort be stable, (i.e., sorting may not reorder equal values). The algorithms `viSort`, `vqSort`, and `vbSort` are stable, but our C++ implementations were not. Without requesting the permutation vector, you cannot tell whether the sort is unstable.

Rather than recode the insert/quick sorts, we chose to make a post-sort pass over the sorted data, detect runs of equal values, and then sort the corresponding subvector in the parallel index vector. Any `NaN`s presented a special problem because they had been interspersed with the `Inf` values. We found the original indices of the `NaN`s and `Inf`s by making a pass over the original unsorted data.

### Extending Array Sorts to Complex Values

The MATLAB algorithm for complex values first sorts by absolute value, then sorts by phase angle, and finally sorts to ensure stability for duplicate values. We decided to leave the MATLAB C code for sorting complex double values in place and call it from the jackets described above.

### Implementing Ascend and Descend Flags

We passed the direction flags '`ascend`' and '`descend`' to the jackets. If the flag is '`descend`', the output is reversed just before being returned. After the reversal, we ran the index vector sorts to restore stability.

### Testing the Sorts

We tested for the following criteria:

• The sorts should accept all valid inputs (12 types of data, shapes from empty arrays to large ND ensembles, and valid dimension and direction flags) and otherwise fail with an error message.
• If requested, the sorts should produce the permutation vector.
• The sorts should “run as fast as C.” Suppose T means “type” and S means “shape.” Then we guarantee the following.
Statement Guarantees
r = sort(v, dim, flags) T(r)=T(v) ^ S(r)=S(v)
[r,i] = sort(v, dim, flags) T(r)=T(v) ^ T(i)=double ^ S(r)=S(i)=S(v)

We checked that in each case the outputs were sorted and the permutation vectors applied to the inputs, given the outputs. We coded all the test routines in M, using the MATLAB functions `class` and `size`.

To test performance, we raced our algorithm against the C standard library `qsort` and the C++ Standard Template Library `sort` for simple cases. The MATLAB `sort` function is about five times faster than C `qsort` and a little faster than the C++ STL `sort`.

MATLAB provided a productive environment to explore enhancements to the `sort` algorithm. Prototyping the possibilities in MATLAB before recoding the algorithms in C++ helped us quickly understand our programming choices. We can still use the prototype M versions of `sort` as test components when we verify the MATLAB `sort` function’s behavior. Using MATLAB to upgrade MATLAB enabled us to deliver a higher performance `sort`, and we can use M-code to explain how we did it.

````find`

Similar to the MATLAB function sort with options 'ascend' and 'descend’, the function find has added
functionality in MATLAB 7 that lets you find the 'first' or 'last' elements that meet some criterion. For
example, if you want to find the first element in array X less than 42, use

`find(X<42,1,'first')`

This operation stops as soon as the first element meeting the condition is encountered, which can
result in significant time savings if the array is large and the value you find happens to be far from the
end.

If, instead, you want to find the final 10 elements > 17, then use

`find(X>17,10,'last')`

Again, the operation returns as soon as 10 values are found. If fewer than the requested number is found
after the entire array is scanned, the return vector will have fewer than the requested number of elements.

Before `sort` had the option to return items in descending order, you had to perform the  `sort `and then
index from the end. With MATLAB 7, `sort` has the new '`descend'` option, and `find` has the new `'last'` option to
combine the find with the indexing into a single function call.
```

Published 2004