I would assume scipy uses MKL, and Julia has wrappers for that too: https://github.com/JuliaSparse/MKLSparse.jl. For the GPU CUDA.jl has support for sparse matrices as well.

My only problem too is I just want to write high level mathematical code and have it work. With Python's anaconda distribution, I just download it and the sparse stuff works. I don't want the hassle of dealing with wrappers and code that is not very mature. I know it's a chicken and the egg thing.

Have you even read the readme in that link?

Start a julia repl, hit ] for the package manager, paste `add MKLSparse`, write `using SparseArrays, MKLSparse`, run `sprand(100,100, 0.1) * rand(100)`

I think the functionality of Julia's sparse arrays is mostly on par with scipy.sparse, it's just a bit rough around some edges still, and spread out into non-base packages (e.g. into IterativeSolvers.jl, ARPACK.jl, MUMPS(*).jl etc) of various maturity. Having the sparse arrays built-in might (hopefully) lead to a large ecosystem of interoperable sparse packages, but as with much of the Julia ecosystem this is still a work in progress.

Personally, working on sparse eigensolvers, my only complaint is that IterativeSolvers.jl is still somewhat lacking (though the ARPACK.jl bindings are more comparable to what scipy does), but otherwise cannot imagine going back to python for this kind of low-level numerical research. Having code run comparably to fortran speed (barring some memory annoyances) is huge for working on numerical algorithms.

The nice thing about julia is that its much simpler to roll your own things. With numpy, you might have to be writing C.

Just the opposite, because more people know C than Julia.

A Julia user knows more about Julia than a Python user knows about C.

So more people knowing C than Julia somehow makes it easier for me? Absolute nonsense.

They do, what's your point.

My point is that there are far more C programmers and C code in the world than there are Julia programmers and Julia code. This is the advantage of Python: it was able to leverage billions of lines of C code by simply acting as a glue language.

But I don't know C, and it's easier for me to be able to write fast code in a language I actually know.

Having to switch to a language you don't know is not a good situation, it doesn't help that a lot of other people know that other language.

Just fyi, Julia has a fantastic and very fast C foreign function interface, just like Python.

> With numpy, you might have to be writing Python

FTFY cf Numba

http://numba.pydata.org

Numba is great (as well as Jax, which is even better), but its limitations are drastic compared to Julia. So much of what is great about python (e.g., compossibility, introspection, simplicity) is lost with these libraries, while Julia still retains them.

The main thing missing from Numba is user defined structures, the jit classes are a performance penalty usually.

That said, in terms of composability you can jit over the closure to achieve a lot of what you might want, e.g.

    def make_loop(f):
        @jit
        def fn(x):
            for i in range(x.shape[0]):
                 x[i] = f(x[i])
        return fn

for any jit'd function f will be just as fast as if you had inlined the body of f.

For introspection and simplicity, I think, in high performance, you simply have to choose two of fast, simple and generic. Julia clearly chooses fast and generic.

Could you elaborate on why you are saying julia is not as simple? Certainly, there are performance tips one needs to be aware of, but that is the case with python too. I have used python for relatively high-performance numerics for research work for the last decade, and now that I am exploring Julia, it certainly addresses all of the concerns I mentioned above with code which is at least as legible/simple and much more introspective. It has its warts, but those warts are on the roadmap and I have seen significant improvement between v1.3 and v1.6 (e.g. compilation latency and debuggability).

Compiler errors can be on par with heavily templates C++, with similarly difficult to read code bases. It’s of course easier to work with than C++ but not the panacea that seems to be implied in a lot of the comments.

In terms of deployment, you essentially have to have Julia installed wherever you want to run, which frequently ok, and if not there’s always PackageCompiler except it’s not that easy to get a working shared lib.

There are other things I find complex but probably because I’ve used Python too long so Julia is different etc. I think Julia is a great choice for HPC generally except for challenging SIMD codes where it’s hard to vectorize without an explicit ILP model.

How does Jax lose composability or introspection?

E.g. jax does not autodifferentiate anything that is not jax (scipy ode solvers, special functions, image processing libraries, special number types (mpmath), domain-specific libraries). Compare that to Zygote.jl

It is true that Jax cannot differentiate through C code. But it can differentiate through python code that was written to accept Numpy.

Which is extremely limited compared to Zygote which can do custom types, dicts, custom arrays, complex type system and multiple dispatch uses etc

Try reading the docs before making sweeping negative comments about what a piece of software can and cannot do.

https://jax.readthedocs.io/en/latest/notebooks/autodiff_cook...

Are you talking about "Differentiating with respect to nested lists, tuples, and dicts" from that page? The comment to which you are responding covers quite a bit more. The jax documentation specifically says "standard Python containers". Zygote.jl and other less stable Julia auto-diff libraries go far beyond the built-ins and can work with structures defined by packages never designed to be used with automatic differentiation. Of course, there are limitations, but quite a bit less severe than the one in jax (and again, I am saying this while being a big jax fan).

As the document I linked to says, Jax autograd supports custom data types and custom gradients.

It’s honestly exhausting arguing with all you Julia boosters. You can down vote me to hell, I don’t care. I’m done engaging with this community.

You all are not winning over any market share from Python with your dismissive, arrogant, closed minded culture.

I understand you are frustrated, however, please remember

> Please don't comment on whether someone read an article. "Did you even read the article? It mentions that" can be shortened to "The article mentions that."

> Please don't comment about the voting on comments. It never does any good, and it makes boring reading.

> Please respond to the strongest plausible interpretation of what someone says, not a weaker one that's easier to criticize. Assume good faith.

https://news.ycombinator.com/newsguidelines.html

I am confused why you assume I am a "Julia booster" or use such combative language. I love Python and Jax and use it for much of my research work, I just also like learning of other approaches. Please try to honestly address the sibling comments. We have repeatedly claimed that tools like Zygote.jl can autodifferentiate efficiently things that Jax can not (without a lot of extra special code and hand-defined backprop methods), e.g., an array of structs with scalar and vector properties over which a scalar cost is defined. Just give examples, so that we can both learn something new about these wonderful tools instead of using such offensive language. It is hard to not take your own comments as the ones being dismissive.

Also, look from where this conversation started. My claim was that jax does not work with "(scipy ode solvers, special functions, image processing libraries, special number types (mpmath), domain-specific libraries)". A julia library does not need to know of Zygote.jl to be autodifferentiable. A python library needs to be pure-python numpy-based library to work with jax.

In order to try to contribute to the discussion: I think this paper describes relatively well what is so special about the Julia autodiff tools: https://arxiv.org/abs/1810.07951

For a separate approach, which is also very original, check out https://github.com/jrevels/Cassette.jl

I don't see anywhere about a better type system or multiple dispatch. Try being less salty

In what language are you defining these custom arrays or types? certainly not in python, or they'll be too slow to be worthwhile.

Indeed, but "python code written to accept numpy" is a pretty restrictive subset (comparatively; I do still enjoy using python). It does not even cover most of scipy, let alone the domain specific libraries, which frequently end up using cython or C for their tightest loops.

But then writing something in C is not really difficult isnt it? (The real difficult part is the algorithm not the implementation in a language)

I'm a pretty competent C programmer, but doing even fairly simple matrix programming in C is fiendishly hard not to screw up. When I was porting this code [1] to this [2], it took an absurdly long time to iron out all the bugs. And just look at the difference in verbosity and difficulty. Worst of all, C has no guard rails: it's so easy to think you have it right, only to realize much later (or not at all, frankly) that you're computing complete nonsense.

[1] https://github.com/JuliaLang/Microbenchmarks/blob/af3d18f7b3...

[2] https://github.com/JuliaLang/Microbenchmarks/blob/af3d18f7b3...

Thanks for the example. Could you also tell how much the speedup of this led to?

That C code is 39% faster than the Julia version, so faster but not massively so. You can optimize the Julia code to be as fast as the C code by using in-place operations, but it's probably not worth if for such a small speedup unless the code is a real bottleneck. By comparison, the Python version (using NumPy, of course) is 25x slower than C, which is a whole different situation.

Most the speed advantage of the C here is due to reusing the same memory over and over, which you can do pretty easily in Julia as well. Here's the Julia code using in-place operations and it's still quite a bit more readable than the C version: https://gist.github.com/StefanKarpinski/e57f5a36b7890b261a0d.... I timed this and it's the same speed as the C version. When developing this, it follows the same general outline as the C version, but you have several benefits: (1) You can use asserts to compare to the easy version; (2) There are niceties like bounds checks and array indexing. In C you can't do (1) because there is on easy version to compare with. And doing the array index computations in C is kind of a nightmare—it's so easy to accidentally screw them up.

I think you're overestimating how many Python developers are willing to write and compile C. Debugging gets more difficult across two languages. Shipping binaries with a Python package correctly can be a hassle too (what about -march=native performance).

That is not too hard if one uses Pybind11 C++ library, which allows creating native code Python modules (shared libraries) in a simple and declarative manner. Other library that can be used with Pybind11 is Eigen for linear algebra that can perform loop fusion in numerical computation, for instance, it can add multiple vectors of same size in a single loop.

Numba compiles Python to LLVM as a backend, why is this even a discussion?

Because numba (and jax), although being wonderful pieces of engineering, work on a fairly limited subset of python, do not permit the same level of introspection as python, and do not play nice with other python libraries. Julia does not make that sacrifice.

For Jax I believe this is false.

Jax is composable. In fact it’s a core design goal. Jax arrays implement the Numpy API. I routinely drop Jax arrays into other python libraries designed for Numpy. It works quite well. It’s not effortless 100% of the time but no library interop is (including Julia multiple dispatch).

I can introspect Jax. Until you wrap your function with jit(foo) it’s as introspectable as any other Python code, at least if I’m understanding what you mean by introspection.

Jax has implemented most of the Numpy functions, certainly most of the ones anyone needs to use on a regular basis. I rarely find anything missing. And if it is, I can write it myself, in python, and have it work seamlessly with the rest of Jax (autodiff, jit, etc)

Jax is awesome. But supporting most of numpy isn't enough. Because numpy isn't composable. You want to add banded-block-banded matrices to numpy? Then you need to fork numpy (or in this case fork jax); this is a package in julia and it works with everything. You want add names to your array dimentions like PyTorch recently did, then like PyTorch you need to fork numpy; again this is package in julia. You want to do both? You have to merge those two forks into each other. In julia this isn't even a package this is just using the aformentioned two packages.

You want to work with Units or track Measurment error (or both?). Basically same story. Except better in some ways worse in others. Better because you don't have to fork numpy, it is extensible enough to allow that. Packages exist that use that etendability for exactly that. Worse because those are scalar types, why are you even having to write code to deal with array support at all. Agian 2 julia packages and they don't even mention arrays internally.

The problem's not Jax. The problem is numpy. Or rather the problem is this level of composability is really hard most of the time in most languages (including the python + C combo. Especially so even).

Its true that this is not always trivial 100$% of the time with julia's multiple dispatch. but it is truer there than anywhere else i have seen.

> do not play nice with other python libraries.

Julia hasn't been around long enough to build ecosystem with multiple commerical giants building competing products.

> Julia does not make that sacrifice.

Sure, Julia sacrifices your sanity on the alter of stack traces, method resolution, time to first plot, among others.

> fairly limited subset of python

Fantastic comment, always brought up, except that this same subset of Python (really how to use CPU efficiently) is about the same that anyone writing about Julia performance is preaching.

I think a really good example of how Julian tools are composable and work with the whole of the language, while that is not the case with Pythonic tools, is autodifferentiation. The many different julian autodiff libraries have little trouble performing the autodiff through special functions, weird data structures, and 3rd party libraries. In python, none of the autodiff tools play nice with scipy or with mpmath or with sympy (or qutip, which is something I need). This is not about commercial competitors either - all of the libraries I mentioned are non-commercial.

Julia provides the interoperability between such libraries, including custom datatypes, frequently with fast zero-cost abstractions.

I agree but think that Julia had the benefit of watching Python and MATLAB suck at auto diff and do better as a community. That said PyTorch and Jax are great.

> fast zero-cost abstraction

Until you get a compiler error and the stack trace takes up a few screens of text. Who maintains that abstraction by the way? Nothing is free.

>Fantastic comment, always brought up, except that this same >subset of Python (really how to use CPU efficiently) is ?about the same that anyone writing about Julia performance is >preaching.

This is completely and demonstrably false. Julia allows fast programming with macros, multiple dispatch, abstract types and more.

All those features have a cost that, in a performance oriented setting, require restricted the code to a subset of Julia that matches what the hardware is fast at.

Where are you getting this from? None of those features have a performance cost, if used according to some simple guidelines.

In fact, multiple dispatch is essential to get good performance, abstract types have a cost if explicitly used inside struct definitions, but that's why you use parametric types, and abstract types do not have a cost in function signatures.

Macros are expanded at compile-time, and at least have no runtime cost, in fact they are often used for improved performance.

I think there's some fundamental misunderstanding going on here, but I'm not sure what it is. Or do you mean to say that if it is possible to write slow code in a language, then only a subset of that language has good performance? If that is the case, I don't know how to respond.

No, that's not true. I specifically mentioned features that can easily be made zero cost on both CPU and GPU targets.

In fact, those aren't special features and there is no other subset. Those are the core abstractions on which everything in Julia, down to primitive types and wrappers for llvm intrinsics are built. Without them you wouldn't have Julia.

Julia's GPU ecosystem wouldn't be where it is with just one or two people maintaining it without being able to reuse those features to plug into existing machinery.

I beg your pardon. What parts of Julia are not supported by its jit compiler?

Numba isn't as slam dunk of a solution as it seems, at least when your code needs to run on something besides x86. I recently witnessed someone getting numba based code running on an ARM based system, and they fought llvm package and performance issues for weeks. The code in question was generally working and seemed pretty good on x86, but took up multiple CPU cores at 100% on the ARM platform. Doing a straight forward port of the algorithm to C++ and wrapping it with pybind11 resulted in the same algorithm using 14% of one CPU core. Perhaps the numba implementation could have been made to work eventually, but it seemed pretty difficult to debug due to all the layers involved.

A subset of Python for one, and debugging stays awkward with gdb.

Makes me think of this little pun: https://github.com/FluxML/Flux.jl/blob/master/src/utils.jl#L...

The subset of Python compatible with Julia's whole performance story (mostly, structs are the missing piece), yes sure.

You mean the type system and multiple dispatch? I guess those are fairly minor points, yes.

One performance feature I find really nice in Julia, is that you can write your own function (cost function, function to integrate, whatever) that you can pass into a solver/integrator/etc. Then it will be inlined and compiled into that external library.

Does this work with numba? Are your jit-compiled functions compiled into external modules?

Even if structs worked, they wouldn't be part of a good type system, don't work with abstract interfaces/inheritance,traits, multiple dispatch, fast union types, zero allocation immutables or any myriad of other things

That really depends but on average it will be more difficult and there will be a lot of boiler plate.