Julia motivation: why weren't Numpy, Scipy, Numba, good enough?

Community
1

I recently gave a Why Julia? talk to a roomful of computational fluid dynamicists and astrophysicists. (John Gibson, UNH & KITP, Julia - a high-performance dynamic programming language for technical computing, Kavli Institute, UCSB 2017-02-02).

I got two main skeptical responses. One Python expert in the audience pointed out that much of what I was touting in Julia is available in Python: Numba for fast array math, %jit macros for just-in-time computation. I read up on these things afterwards and now have some questions about why the Julia developers thought starting afresh was the right move. I lean that way myself, based on experience(in the long run, clean, total solutions beat patches on imperfect ones). And from what I’ve read there are compatibility issues between the multiple Python extensions (Cython, PyPy, Numpy, and Numba), and maturity issues for Numba. But I imagine the Julia development team has very specific and clear ideas of the shortcomings of Python’ that warranted starting afresh with Julia. If anyone has time to spell these out or point me to a previous description, I’d be grateful.

The other skeptical response was about parallelism. I’ll post about this elsewhere. Suffice it to say that people in my community want to see a working example of a parallel PDE solver, and for this end I am thinking of implementing a Julia code for isotropic turbulence in a triply-periodic box. Should be easy if I can get Distributed Arrays and parallel FFTW to work together.

best regards,

John Gibson

26 Likes
2

I’m not a julia developer, but to answer the question about ‘not good enough’ i think you need put the questions and answers on a timeline. The ‘why we created julia’ blog entry dates back to feb2012 and at that time julia was already existing ~2 years. Numba was not seen.

2 Likes
3

This is pretty good reading: Personal thoughts about Pyston’s outcome – kmod’s blog

"The thing I wish people understood about Python performance is that the difficulties come from Python’s extremely rich object model, not from anything about its dynamic scopes or dynamic types. The problem is that every operation in Python will typically have multiple points at which the user can override the behavior, and these features are used, often very extensively. Some examples are inspecting the locals of a frame after the frame has exited, mutating functions in-place, or even something as banal as overriding isinstance. "

Sure you can make Python fast in some cases where you use a small subset of the language. However, Python as a whole is extremely difficult to make fast because of the incredible freedom you have. The “compiler” can assume very little, so it is difficult to generate efficient code. Just consider the incredible amount of time and money that has been spent on trying to make it faster and the relatively mediocre success it has had (Pyston, Numba, PyPy etc).

This is also a good watch: Jeff Bezanson - Why is Julia fast? - YouTube

19 Likes
4

Not a developer, but I have these kind of discussions with colleagues (economists).

My take on this is that if a language/environment is “good enough” for someone, no amount of persuasion will convince them to switch. In my experience, people who gravitate to Julia do so because it really solves the problems they had with other languages (in terms of speed, elegance, abstraction), or expect that it will, given time, so it is worth investing in.

The ideal Julia user at this stage should be someone who, when reading the manual for the first time, will just scream in delight every few minutes, because they find a solution (or an architecture for possible solutions) to language problems they have encountered. Otherwise, the language will just look like an unnecessarily complex construct with syntax/semantics that changes every few months, and lots of work in progress code in the package ecosystem.

23 Likes
5

@kristoffer.carlsson agreed; that article about pyston was the first thing that came to my mind.

Another piece of the “why Julia” puzzle is the fact that Julia does not treat a specific set of built-in types differently. That is, I can make my own type or immutable that can do anything that Julia’s built-in types can do, and my type can be just as fast as the one provided with Julia (if I were as skilled at numerical computing as the devs, which I am not). This isn’t just a nice feature; it’s essential to making things like automatic differentiation and other numeric types fast and easy to use. It’s also why we have really nice abstractions like the idea that an image is not a matrix of doubles or bytes but instead a matrix of colors. That works because Julia lets us define a whole family of color-related datatypes which provide useful behavior and abstraction at little or no run-time cost. That’s been the biggest benefit of Julia, in my experience.

Oh, and also JuMP. JuMP is great.

16 Likes
6

I think Julia’s advantage is some synergy between its type system, fast code on custom types and new abstractions, syntax, multiple dispatch, macros etc

Thus, its not just fast on specific types of code, but allows for unparalleled flexibility, code reuse, expression of ideas and abstractions, higher order programming etc in a way that makes numba and python seem stone age.

Here is a concrete elucidating example of something that wouldn’t be possible in python: 6 Months of DifferentialEquations.jl: Where We Are and Where We Are Going - Stochastic Lifestyle

1 Like
7

This is how I understand the “why Julia?” question. This was my round 2+ after I didn’t have an adequate answer written down the first time I did this workshop.

http://ucidatascienceinitiative.github.io/IntroToJulia/Html/WhyJulia

Julia essentially answers the Pyston question: how should a language be restricted and thought about so that way it still looks and acts dynamic, but in the end is very easy to extend and optimize? When you understand how Julia’s type system allows for JIT compilation to not just work, but work very well by specializing on exactly the functions its seeing, then you see what makes Julia special and how to apply this design to your own codes to get similar results.

I was recently giving a talk about DifferentialEquations.jl, at a part where I was highlighting the use of odd-number types (SymPy expressions, ArbFloats, ApproxFun Funs) in the differential equation solvers. I got an immediate question from the audience along the lines of

I know from Python that if you objects inside of numerical methods / loops, it will immediately be very slow. And it’s difficult/impossible to get user defined objects in something like Numba. So while all of this stuff is really great for expanding what we can solve… won’t it be very slow?"

Of course I had to answer (along the lines of)

That’s a limitation of Python. Julia’s type system is specifically designed so that way user-defined types are just as fast as those which are “built-in”. In fact, Float64 is just a type which is defined inside Julia. You can, the same way Base was built, build up your own crazy numbers and arrays, have them be fast, and they will “just work” via dispatch. ApproxFun is probably coolest example: it’s type Fun is essentially a “Function in a function space”, but from Functional Analysis we know that those functions are just points (numbers) in this new space. So you can define +, -, /, *, etc. between them… and then just stick that “number” in functions or libraries written for Numbers and get out the solution to a PDE.

The key is understanding why the abstraction works. Because Julia recompiles what each of those operators mean, when you put in some complex type it essentially compiles the code you would have wanted to write. When you start to get the hang of this, your Julia code can become bizarrely abstract, and you can compose interesting type + generic function to get fast implementations with no work. A good example of this is SymPy.jl or SymEngine.jl. These implement symbolic mathematical expressions in Julia. When you first start using them, it seems odd that the documentation is so little. But then you understand that “oh, +, -, *, and / is defined on a symbolic expression, and therefore to take the inverse of a matrix of SymEngine expressions, I can call Julia’s Base inv function. Nobody ever needed to implement such a function because the way inv is written only requires those dispatches.” And when you benchmark it, it’s super fast (it can sometimes be faster than built in SymPy expressions…). All of the “linear algebra fallbacks” will make matrices of SymPy/SymEngine expressions “just work”. You can string together fast algorithms using methods nobody ever wrote but the compiler makes for you from generic code. That’s something you don’t get with Python.

All of this “zero-cost abstraction” allows large software ecosystems like JuliaOpt, JuliaDiffEq, JuliaStats, JuliaML, etc. to be developed very fast, but have amazing flexibility and performance. At this point, I think we recommend Julia for its packages, since there are quite a few things which Julia packages can readily do that I wouldn’t know simple ways of doing in other languages.

20 Likes
8

In addition to performance reasons indicated by someone, I’d point also to maintainability, mostly thanks to a well devised type system. I’m not a Python user, so please do correct me if I’m wrong anywhere, but I’ve got the feeling that it’s often tricky to make two classes from different packages play nicely together. At least, some glue code is sometimes required.

In this regard, I find the example of my package Measurements.jl instructive. There is a similar Python package: uncertainties. They’re both uncertainty propagation library, but they have different support for some features:

  • complex numbers. uncertainties: not supported. Measurements.jl: work out-of-the-box.
  • arbitrary precision. uncertainties: not supported. Measurements.jl: work out-of-the-box.
  • arrays of numbers with uncertainty and operation with them. uncertainties: support to NumPy thanks to a submodule written for the purpose. Measurements.jl: work out-of-the-box.
  • linear algebra operations with matrices: uncertainties: limited support: a submodule written for the purpose provide methods to compute matrix inverse and pseudo-inverse. Measurements.jl: those operations defined for any AbstractFloat work out-of-the-box, things will improve if/when issue #5429 will be fixed. Note this is not a limitation of the language, it’s simply that we need someone who writes the methods for generic types.
  • operations with numbers with uncertainty and physical units. uncertainties: you have to write glue code, Pint package does this. Measurements.jl: works out-of-the-box with any type providing unit feature, see Examples — Measurements.jl 0.4.1-dev documentation.
  • numerical integration. uncertainties: not supported. Measurements.jl: quadgk from QuadGK.jl already works out-of-the-box with functions returning arbitrary types, including Measurement. In Measurements.jl I only taught quadgk the derivatives of the integral with respect to the interval bounds, so that one can use also Measurement objects as endpoints of the integral.

When I say that Measurements.jl works out-of-the-box I mean that I didn’t add any single line of code to get that feature. They all came for free during development. This is a huge leap for the developer, and users as well (see second last point: Measurements.jl, nor the physical units packages, doesn’t impose a choice on the user, they can use whatever package they like).

15 Likes
9

I think this is a really key point. I would even go as far to say that most people likely “don’t need” Julia. Most users in the general realm of scientific programming likely are stringing together results from various packages. And for those people, sure if there’s already a package in Python, why switch?

The problem comes when there isn’t a package. The next group of users are those who wrote some simple libraries in Python/MATLAB/R. This is the group that gets buffed when trying to make Cython/MEX work, finds out that PyPy doesn’t work with some libraries they want to use, and tries to use Numba but runs into difficulties the moment they have “non-trivial codes” that aren’t just a benchmark loop and involve some objects. These are the people who find Julia and go “yes, this is what I needed”.

Overtime, these converts make(/made) some really interesting packages in Julia, and that gives the incentive for the “package user” group to switch. But I think that it’s really important to understand that the vast majority of users of the language will probably not know what “broadcast fusion” is or that it occurs, that in-place operations / mutating functions are the “fast way” to write array functions, etc.

And that’s okay! It just means that Julia, because it solves the two language problem, also has a two audience problem, and one should evangelize Julia by appropriately discussing Julia to the different audiences. Some people want to hear about compiler details and type systems for writing libraries. Some people want to see cool examples of solving difficult optimization problems in JuMP, complex simulations using DifferentialEquations, and easily train large neural nets using MXNet. Both of these groups are Julia’s (potential and current) audience.

(I have now written too much so I will stop :slight_smile:)

15 Likes
10

I guess I can give some historical motivation here. When we started working on Julia in 2009, the Python ecosystem was nowhere near where it is today:

  • PyPy was pretty new: 1.0 was released in 2007 – it was the first version to be able to run any substantial Python apps and there were no plans at the time for supporting any of NumPy;
  • Cython was even newer: initial release was 2007;
  • Numba didn’t exist: initial release was in 2012.

Even NumPy itself was a bit rough back then. It’s entirely possible that if the SciPy ecosystem had been as well developed in 2009 as it is today, we never would have started Julia. But then again, there are a lot of aspects of Python that can still be improved upon for scientific computing – and for general computing as well. Note, however, that Julia wasn’t a reaction to Python – @viralbshah has used Python and NumPy extensively, but I had not, and as far as I know, neither had @jeff.bezanson. Julia was much more influenced by Matlab, C, Ruby and Scheme, although we certainly have looked to Python for inspiration on some designs (I straight up copied the path and file APIs, for example). Python is often a good reference for basic sanity in language and API design choices. Here are some of the areas where we felt Julia could do better than existing languages, Python included.

Real metaprogramming. There’s no good reason for a dynamic language not to do this, but aside from Lisps, they don’t. I’m not sure why. The greatness of JuMP is due to Julia’s metaprogramming facilities – and, of course, the brilliant API design work of the JuMP folks to make effective use of metaprogramming to allow targeting different optimization backends efficiently.

Make numerical programming a priority. Technical computing deserves a really first-class modern language that’s focused on it. After a couple of decades, Python finally threw NumPy a bone and added the @ operator for infix matrix multiplication. Numerical computing is just not a priority for the Python core devs – and reasonably so, since that’s not the language’s primary focus. Julia, on the other hand, has been about numerical and scientific work from the start, so we care about usability and convenience in those areas. This is not a diss on Python in any way, it’s just a fact that numerical programming will never be a top priority for Python the way it is the top priority for Julia.

Have an expressive built-in language for talking about types. Every numerical computing library ends up needing to talk about types. At a minimum, you can’t do anything practical without typed arrays and being able to pass them to C and Fortran libraries like FFTW and BLAS – both of which require talking about types. As a result, there seems to be a Greenspun’s Tenth Rule analogue for numerical libraries, which I’ll posit as follows:

Any sufficiently complicated numerical library in a dynamic language, contains an ad-hoc, informally-specified, bug-ridden, slow implementation of a type system.

I just made this up, but it’s true – the SciPy ecosystem is littered with dozens of incompatible DSLs for talking about types with varying degrees of complexity and expressiveness. Python 3.5 has standardized on syntax for type hints, but from what I can tell, it still isn’t sufficient for what most numerical libraries need and probably never will be (priorities again) – it’s more designed for type checking in the presence of duck typing. If it quacks like a Float64 and walks like a Float64, is it a Float64? Not necessarily, if you really need your code to be efficient. Julia circumvents Greenspun’s Eleventh Rule (the one I just made up) by fully committing to having a type system and making the most of it for dispatch, types, error checking, code generation, and documentation.

Don’t make numeric types or operators special. What does it take to make “primitives” like Int and Float64 normal types and operators like + normal functions? This requirement leads you seemingly inexorably to dynamic multiple dispatch – which happily turns out to also be incredibly useful for many other things. Not least, neatly addressing the expression problem, which allows Julia libraries to share types and generic algorithms remarkably effectively. Multiple dispatch was a strange experiment at first, but it has worked out surprisingly well. I’m still regularly shocked that more languages haven’t used this as a paradigm before Julia. Dylan is the only language that’s taken multiple dispatch nearly as seriously as Julia. (And apparently Perl 6, from what I hear.)

Use modern technologies and design for performance from the start. These seem separate, but they’re not really. Language performance dies the death of a thousand cuts. No single design choice destroys performance – it’s an accumulation of little choices that kills it, even though each would individually be survivable. These choices are made in the context of technologies. Designs that are easy to implement with one technology may be hard to implement with another. So we should design for modern, high-performance technologies from the beginning.

Most dynamic languages were designed to be interpreted, and interpreters are relatively slow. So no one thought that much about speed when these languages were designed. Later, when you try to apply modern technologies like JIT and LLVM, choices that were made on the assumption that a language would be interpreted are often a big problem. For example, representing objects and scopes as dictionaries is a reasonable choice in an interpreter. From there, it’s perfectly natural to let people manipulate objects and scopes just like dictionaries – after all, that’s what they are. But making code fast requires making local variables and objects completely vanish – the fastest computation is the one you don’t do. But once you’ve committed to exposing dictionary semantics and people have written tons of code using those features, it’s too late – you won’t ever be able to effectively optimize that code.

Julia, on the other hand, has never had a full interpreter. When we started the project, we looked at compiler frameworks, picked LLVM, and started using it right away. If a feature wasn’t easy to generate code for, it was never added to the language in the first place. Often there are slight variations on features that make them vastly easier to generate code for but don’t really eliminate any expressiveness. eval that operates in local scopes? Impossible. eval that operates only at global scope but allows you to define new functions (which have local scopes)? No problem. A lot of Julia’s design fell directly out of making it easy to generate LLVM IR – in some ways, Julia is just a really great DSL for LLVM’s JIT. Of course, code generation frameworks aren’t that different, so Julia is really just a nice DSL for any code generator that can be invoked at runtime.

103 Likes
12

@StefanKarpinski, very enlightening writing.

I really like the fact that Scientific Computing is a first priority.
Going forward, do you think it is achievable, in the future, to write code in its Vectorized Form + Decorations in Julia and have performance of highly optimized (SIMD + Threading) C Code?

I liked the guide line that what sacrifices performance won’t get in.

Python, in my opinion, went too fast with its flexibility and types of variables.
Sometimes its code looks like Black Magic.

I prefer language with less built in types which sometimes means you code more but also the code is simple to grasp.
I hope Julia will remain simple in its Base.

@John_Gibson,
I think Julia’s charm is how deep it lets you go.
For instance, @inbounds macro, Choice whether to have SIMD or not, use Multi Threading or not, etc…
It means that on the surface you have the usefulness of MATLAB like language yet when you want to speed things up and get your hands dirty you can do it within the same language, Julia (As opposed going C in MATLAB).

Yet still, the challenge of Julia will be to close the gap of performance when writing Vectorized “MATLAB Like” code to what can a user get when going down.
Hopefully it is doable and if it is doable, Julia seems to be the best candidate to achieve it.

Thank You.

2 Likes
21

A key point around the question of “will julia get various advanced optimizations?” is our philosophy and priorities. If a language+compiler can optimize more programs more effectively, we consider it better, and we want it. Somewhat surprisingly, not everybody agrees with this. The other camp is the “scripting language” camp, which says “programs should be written in two languages” (this is a John Ousterhout quote). You have a performance language, and a scripting language. Any thought to optimizing the “scripting” language is almost by definition a waste of time, since that’s what the performance language is for. This is where Python’s occasional “fast enough for me” attitude comes from.

Unfortunately, with the hardware changes happening now, I see a potential new “two-language problem” around accelerators like Halide, TensorFlow, ArrayFire, ParallelAccelerator, etc. We are starting to need optimizations that don’t naturally fall out of what we know about functional and object-oriented programming. Maybe we don’t know how good we had it — when all you needed for performance was C, well at least that’s a powerful, widely-adopted, general-purpose language. Now, we are willing to cram whatever hacks are needed into julia to get the best performance, but I have not yet seen a satisfying, disciplined approach to this problem.

23 Likes
27

Yes. C++ overloading is multi-argument but is also all static — it only works on types known to the compiler. If you want dynamic dispatch, you need to use virtual methods, which dispatch on a single argument and use different syntax. Dynamic dispatch appears to be important. People want to just throw whatever they have at some functions, and have the right thing happen, regardless of whether the code can be statically typed. But in C++ static type errors are not even the worst of it: the compiler will pick an overload to call based on whatever type it can determine, which can be different from the run-time type. That doesn’t seem to be a huge problem in practice, but it’s very unsatisfying. Only the ground truth of what an object actually is should matter, not compiler limitations.

22 Likes
29

I’ve noticed this too, and I think it’s a concern and an important problem. The most optimistic scenario I’ve come up with is that these kinds of optimizations are perhaps best incorporated into user-level julia code via a set of sophisticated iterators. For example, in ImageFiltering I’ve been able to achieve some of the same kinds of “fusion” as Halide with surprisingly few lines of code. The core element of this approach is itself a small package, TiledIteration, which I think does constitute a reusable nugget of ideas for moving forward in this problem space.

From the standpoint of implementing this kind of thing more broadly and for “micro” scale computations, a crucial optimization will be the whole stack vs heap for structs that contain “pointers” to heap-allocated memory: these iterators need to create a crazy number of wrappers for operations working on small chunks of arrays.

15 Likes
34

This is wildly off-topic. Please respect that a thread has a subject and stick somewhat to it.

4 Likes
35

Different for me. I saw a post saying look at our great new numerical language “Julia”. I thought, ah another one, best of luck to you. A couple of years later, Andy Ferris, who worked down the hall from me, said he had abandoned other languages for Julia. I know his background so I took this very seriously; that was the key factor. I then spent a couple of hours coding and running a small test problem in C++ and Julia (including downloading, building, learning). I did not believe the results. After I convinced myself that I hadn’t made a mistake and what I saw was real, I immediately dropped everything else and never looked back.

28 Likes
47

http://numba.pydata.org/numba-doc/dev/user/jitclass.html

Yes, generality is the key. Numba can get you going for Float64. And it can now, in very limited cases, get you JIT compiled objects, but that’s not even getting close to Julia’s types. I believe you cannot write your JIT functions to auto-specialize according to input types, so generic programming is out of the question. This only works on a small class of objects with limitations.

But more importantly: Python does not have a good language for talking about types. Numba is trying to bolt on features that get it closer to Julia, but it’s missing the language which makes it easy to actually use and discuss. For example, Julia has parametric types. Parametric types can be “infinitely many” different JIT compiled classes, and many times you want to dispatch differently depending on these type parameters. With what exists in Numba, you technically “can” do it, but it’s all manual and at that point you might as well be writing C++ (or… Julia).

So all of the examples that people give for “but Numba can do it” tend to be “here’s an example looping on Float64”. Yes, bravo, Numba got that. But what I have been finding out in my journey with Julia is that, that’s the simplest case (and you might as well just use C/Fortran if that’s all you want).

What’s interesting about Julia is that same exact code is efficient for arbitrary arithmetic, or AbstractArrays (and gets you auto-compilation of your functions to GPU variants using GPUArrays for example). Numba can keep bolting on what’s needed. It has bolted on stuff for GPUs. If they see that people are using Julia’s generic algorithms with DistributedArrays, they can make a distributed array and change the compiler so that case will work. But Julia is designed correctly so that way these kinds of specializations aren’t “top-down”: no compiler changes are needed. You can do all of this by adding a package.

In the end, I am sure that with enough work Numba can keep trying to keep up with “the standard use cases” of Julia that are beyond Float64, but it’ll still be in a language that has no way to discuss what it’s actually doing, and it’ll still be “given to you” by compiler changes in Numba itself, and making changes to add stuff like this won’t be accessible to standard Python developers without compiler knowledge.

7 Likes
59

I think you misunderstood what I said. I did not mean to imply that programmers who are currently not able to appreciate Julia never will, and that they should be ignored. I merely said people look for solutions when they learn more deeply about the problems, which takes time.

Also, I don’t think we should disparage Python or Python programmers. Julia has a lot of promise, but at the same time the core language and the library ecosystem need time to mature. Waiting a bit more with adoption is a reasonable decision in many cases.

8 Likes
62

I will always have a soft spot for Python, but its high-performance story remains messy. My flowchart proves it :joy:

4266×3200 406 KB

22 Likes