4 ms·

As others have pointed out improvements could be made to improve performance of the Cython code. However it seems to me the take away from this post should be i

by hartror 13y ago

As others have pointed out improvements could be made to improve performance of the Cython code. However it seems to me the take away from this post should be in this application Julia has a better performance profile while being as simple and expressive as Python [1]. As a Python programmer that often works with these sorts of problems this reinforces my perception that Julia, while not suited to all programming problems offers a productivity gain when it comes to mathematical/algorithmic computing. Sure I could sit and tweak the Cython code or rewrite the hot loops in C++, but it seems to me Julia hits a sweet spot where I am less likely to have to compromise programmer performance for computational performance. [1] I will note however that the code in the OP doesn't seem like the clearest implementation of the algorithms in either language, some of that heavy nesting should be factored for a start.

ajtulloch 13y ago
In terms of the heavy nesting, that's born out of performance measurements - factoring out the inner loop into a `cdef inline` function in Cython appeared to be a performance hit, so I manually inlined it. When it came for the Julia implementation, I wanted to keep it as similar as possible to the tuned Cython implementation for pedagogical purposes. In general, if anyone has performance improvements or code cleanups to the Cython PAVA code [1] (the implementation used by scikit-learn), please benchmark and send a pull request over to the scikit-learn team. I'm sure they'll be very happy to accept improvements. [1]: https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/_isotonic.pyx https://github.com/scikit-learn/scikit-learn/blob/master/skl...
- sesqu 13y ago
  This isn't github, but I'll point out one cleanup: # solution[i:k + 1] is a decreasing subsequence, so should read # solution[i:k] is a decreasing subsequence, so
StefanKarpinski 13y ago
Another subtle but important point is that the Julia code can be very easily made completely generic without sacrificing any performance or clarity. All you need to do is remove those Float64 annotations and replace the 0.0 values with `zero(y[1]*weights[1])`. With that small change, the same exact code will work for Float32, Rational, BigFloat, Complex (if that even makes sense), etc. – even new, user-defined numeric types that didn't exist when you wrote this code and that you know nothing about, as long as they can be compared, multiplied, added, and divided. This kind of generality isn't necessary for a lot of user-level code, where there's a specific application, but it is quite important when writing libraries that need to be as reusable as possible. And this is where the Cython solution has some problems – sure, you can make it pretty fast, but in doing so, you lose the generality that the original Python code had, and your code becomes as monomorphic as C.
- btown 13y ago
  Cython does have limited support for generics in the form of its fused types [1]; it can be made to generate C code for different primitive types, for instance. It's nowhere near as flexible as Julia in that regard, since you can't redefine operators for working on Complex numbers for instance, but if you want the interoperability and maturity (read: available programmers and existing libraries) that Python offers, you don't need to trade away generics completely. I'd love to live in a completely-Julia ecosystem, but unfortunately that's not quite reality yet, and luckily Cython is more than capable for making do. [1] http://docs.cython.org/src/userguide/fusedtypes.html http://docs.cython.org/src/userguide/fusedtypes.html
heydenberk 13y ago
This Python code is "as expressive" as the Julia code because it's not written in idiomatic Python. Por ejemplo: for i in range(n): solution[i] = y[i] when one could simply solution[:n] = y[:n] Don't mean to bikeshed, but this is one of many cases.
- ajtulloch 13y ago
  FWIW, benchmarks indicated that the former was faster than the latter. Remember, this is Cython, so the goal is to write code that is easily compiled into efficient C code, not to write idiomatic Python.

by hartror 13y ago

ajtulloch 13y ago

In terms of the heavy nesting, that's born out of performance measurements - factoring out the inner loop into a `cdef inline` function in Cython appeared to be a performance hit, so I manually inlined it. When it came for the Julia implementation, I wanted to keep it as similar as possible to the tuned Cython implementation for pedagogical purposes. In general, if anyone has performance improvements or code cleanups to the Cython PAVA code [1] (the implementation used by scikit-learn), please benchmark and send a pull request over to the scikit-learn team. I'm sure they'll be very happy to accept improvements. [1]: https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/_isotonic.pyx https://github.com/scikit-learn/scikit-learn/blob/master/skl...

sesqu 13y ago
This isn't github, but I'll point out one cleanup: # solution[i:k + 1] is a decreasing subsequence, so should read # solution[i:k] is a decreasing subsequence, so

StefanKarpinski 13y ago

Another subtle but important point is that the Julia code can be very easily made completely generic without sacrificing any performance or clarity. All you need to do is remove those Float64 annotations and replace the 0.0 values with `zero(y[1]*weights[1])`. With that small change, the same exact code will work for Float32, Rational, BigFloat, Complex (if that even makes sense), etc. – even new, user-defined numeric types that didn't exist when you wrote this code and that you know nothing about, as long as they can be compared, multiplied, added, and divided. This kind of generality isn't necessary for a lot of user-level code, where there's a specific application, but it is quite important when writing libraries that need to be as reusable as possible. And this is where the Cython solution has some problems – sure, you can make it pretty fast, but in doing so, you lose the generality that the original Python code had, and your code becomes as monomorphic as C.

btown 13y ago
Cython does have limited support for generics in the form of its fused types [1]; it can be made to generate C code for different primitive types, for instance. It's nowhere near as flexible as Julia in that regard, since you can't redefine operators for working on Complex numbers for instance, but if you want the interoperability and maturity (read: available programmers and existing libraries) that Python offers, you don't need to trade away generics completely. I'd love to live in a completely-Julia ecosystem, but unfortunately that's not quite reality yet, and luckily Cython is more than capable for making do. [1] http://docs.cython.org/src/userguide/fusedtypes.html http://docs.cython.org/src/userguide/fusedtypes.html

heydenberk 13y ago

This Python code is "as expressive" as the Julia code because it's not written in idiomatic Python. Por ejemplo: for i in range(n): solution[i] = y[i] when one could simply solution[:n] = y[:n] Don't mean to bikeshed, but this is one of many cases.

ajtulloch 13y ago
FWIW, benchmarks indicated that the former was faster than the latter. Remember, this is Cython, so the goal is to write code that is easily compiled into efficient C code, not to write idiomatic Python.