[ENH] Improve performance of isotonic regression

[ajtulloch]

This supercedes #2940.

PAVA runs in linear time, but the existing algorithm was scaling
approximately quadratically, due to the array pop in the inner loop
being O(N).

I implemented the O(N) version of PAVA using the decreasing subsequences
trick, and the performance is significantly faster in benchmarking.

The benchmarks in this diff (adapted from one of the existing unit
tests) show significant performance improvements

Problem Size	Relative Speedup of PAVA
10	2x
100	3x
1,000	14x
10,000	57x
100,000	561x
1,000,000	4645x

On correctness - unit tests cover the isotonic regression code fairly
well, and all pass before and after the change. It's a fairly well known
algorithm with a bunch of implementations, so I think this is correct.
In coding up this algorithm I made some mistakes and the unit tests
caught the failures, which makes me more confident in the correctness
now. Still, the performance improvements seem suspiciously large.

[coveralls]

Coverage remained the same when pulling 8e25dbc on ajtulloch:isotonic-regression-in-linear-time into 4437d56 on scikit-learn:master.

[ogrisel]

Thanks very much for the fix @ajtulloch. It looks good to me. One cosmetics note: please replace k+1 by k + 1 to be consistent with the project style. Also please add an entry to the doc/whats_new.rst changelog.

@fabianp @NelleV you might be interested in having a look at this.

[ogrisel]

Please make it work under Python 3 by adding parens to print statements.

[fabianp]

I didn't run the code but it looks awesome :-)

[ajtulloch]

To replicate/test for yourself, run

cython sklearn/_isotonic.pyx
make inplace
nosetests sklearn.tests.test_isotonic
python benchmarks/bench_isotonic.py 2 8 1 | gnuplot -e "set terminal dumb; set logscale xy; plot '<cat' title 'performance'"

[NelleV]

It's funny, because I originnally implemented pava, and it was slower than @fabianp's active set's implementation. Can you show us how you generated the benchmarks ?

[NelleV]

My mistake: the benchmarks are in the PR…

[ogrisel]

@ajtulloch could you please commit the result of the cythonization of the edited pyx file in this PR? Please use the latest stable version of cython (run pip install -U cython to be sure).

[ogrisel]

If you want to link to your home page you must give the URL at the end of the file. Otherwise, just remove the link markup.

[ajtulloch]

Hi @ogrisel,

The updated sklearn/_isotonic.c is attached (using version 0.20.1), see https://github.com/ajtulloch/scikit-learn/blob/isotonic-regression-in-linear-time/sklearn/_isotonic.c. GitHub doesn't show it in the diff due the line /sklearn/_isotonic.c -diff in .gitattributes.

I'll update the link in whats_new.rst, thanks!

@NelleV, I suspect the performance differences are down to a few things:

• this is a somewhat tuned implementation of PAVA for Cython - entirely inlined, all type information given, so it's very easy for GCC/Clang to generate efficient code.
• the performance of the algorithms depends on the given dataset, so if you have different example problems we should benchmark those as well.

[ogrisel]

My bad I somehow did not see the C file even though I thought I had clicked on the list of updated files.

[ogrisel]

+1 for merging.

[ajtulloch]

I reran some benchmarks - with the symmetrically perturbed log1p from the unit tests, and a simulated logistic regression dataset, results below.

Does anyone have other suggestions for example problems to test perf on?

[larsmans]

The benchmark script is undocumented and uncommented. It takes three numbers and produces two other numbers, without a hint as to how to interpret them.

[ajtulloch]

@larsmans - I've updated the benchmark script with more usage details, a proper parameter parser, and added comments where appropriate. Thanks for your comment.

[coveralls]

Coverage remained the same when pulling a3dba48 on ajtulloch:isotonic-regression-in-linear-time into 4437d56 on scikit-learn:master.

[larsmans]

Ok, merged as 3753563. Thanks!

[ajtulloch]

Great, thanks again for your comments.