Julia: Faster than Fortran, cleaner than Numpy | Hacker News Reader

Julia: Faster than Fortran, cleaner than Numpy

matecdev.com

281 pointscs702232 comments

This an unfair example of trying to make the numpy readable:

    n = len(a)
    M = np.exp(1j*k*(np.tile(a,(n,1))**2 + np.tile(a.reshape(n,1),(1,n))**2))

This has a common issue in functional programming and there's an easy trick to fix it. Change the big one liner by naming things:

    n = len(a)
    A = np.tile(a,(n,1))
    A_T = np.tile(a.reshape(n,1),(1,n))
    M = np.exp(1j*k*(A**2 + A_T**2))

Much easier to read now, and it even exposes an error (the sqrt is missing!). Before we fix that, let's make the obvious simplification of using the .T for transpose:

    n = len(a)
    A = np.tile(a,(n,1))
    M = np.exp(1j*k*(A**2 + A.T**2))

Then let's use numpy's broadcasting to be more idiomatic:

    A = a[:, np.newaxis]
    M = np.exp(1j*k*(A**2 + A.T**2))

Now let's fix the bug:

    A = a[:, np.newaxis]
    M = np.exp(1j*k*np.sqrt(A**2 + A.T**2))

or even

    A = a[:, np.newaxis]
    root_sum_squares = np.sqrt(A**2 + A.T**2)
    M = np.exp(1j*k*root_sum_squares)

I think the last two of these are fairer comparisons.

--------

And to pitch the hype train for jax: if you want to fuse the loops and/or get gpu/tpu for free, just swap np for jnp:

    A = a[:, jnp.newaxis]
    M = jnp.exp(1j*k*jnp.sqrt(A**2 + A.T**2))

The julia implementation was also suboptimal to be fair. I sped it up considerably here:

https://discourse.julialang.org/t/i-just-decided-to-migrate-...

and then even further here:

https://discourse.julialang.org/t/i-just-decided-to-migrate-...

I'd actually be very interested if there's anything other than handwritten assembly that's faster than the second one I posted.

______________

Regarding your more readable version of the Numpy syntax though, I think you have to admit it's still not as readable as

    for i in 1:N
        A[i,j] = exp((100+im)*im*sqrt(a[i]^2 + a[j]^2))
    end

right?

I'll also note that the author of the post didn't come from the Julia community. He was a heavy user of Python and Fortran, calling Fortran from Python to speed up performance hotspots. He wrote this blogpost as his first introduction to Julia.

So I think if it does cast numpy in a poor light, it was not done so intentionally to make it look bad.

Regarding your performance work, you've nerd sniped me into looking for analytical tricks to speed it up ;) We'll see...

Regarding the readability, I suppose it's a matter of taste at this point, but if I swapped `import numpy as np` for `from numpy import newaxis, exp, sqrt`, then IMO, the numpy is more readable:

    A = a[newaxis]
    M = exp(1j * k * sqrt(A**2 + A.T**2))

But then my tastes also run towards preferring the handwritten definitions to be in terms of vectors and transposes instead of elements and indices, at least until you get to many indexed tensors.

Side note: any idea why the python and fortran are exp(i k... while the julia is exp((100 +i) i...)? Is it something I overlooked?

> Regarding your performance work, you've nerd sniped me into looking for analytical tricks to speed it up ;) We'll see...

Looking forward to it! These microbenchmarks are very fun to explore.

> Side note: any idea why the python and fortran are exp(i k... while the julia is exp((100 +i) i...)? Is it something I overlooked?

Oh, that is a partially applied edit I guess. When the author posted on the julia forum, people pointed out that since the exponent was pure imaginary, it could be speeded up even more with cis(...) instead of exp(im * ...) but then the author claimed that exp((100 + im) * ...) was more representative of his actual workflow and I guess changed the julia version in his blogpost but not the Python or Fortran versions.

I got bored with trying to find an analytical boost, but I benchmarked a couple IMO super readable python versions (basically what's in my original comment after making the (100+i)i change):

https://colab.research.google.com/drive/1ABrZJlm8pwB6_Sd6ayO...

On my macbook, using XLA's jit in python gave about a 12-15x speedup on CPU over OP's solution, which was pretty cool, but I'm too lazy to figure out how to install and benchmark Julia on my machine. Applying a 12-15x speedup would at least beat the Julia MT solution in OP, and you've got to admit `exp(CONST * sqrt(A**2 + A.T**2))` is a pretty clean way to do it.

Then I ran on whatever GPU colab decided to give me (a P100), and for just adding a decorator, it's a 1000x-1900x speedup (better as n goes up). Hence my current honeymoon period with jax. I love the speed vs readability tradeoff.

> I'm too lazy to figure out how to install and benchmark Julia on my machine.

Assuming you're eager to try once you've found out:

You should be able to simply download and unpack a binary from: https://julialang.org/downloads/ I'd strongly recommend going with the current stable release (1.6.1).

To start the Julia REPL:

  bin/julia

To install packages:

  using Pkg
  Pkg.add("BenchmarkTools")
  Pkg.add("LoopVectorization")

then you can copy/paste code from eigenspaces link to discourse to run Julia benchmarks. E.g., you can copy/paste from https://discourse.julialang.org/t/i-just-decided-to-migrate-... if you define

  using LoopVectorization
  const var"@tvectorize" = var"@tturbo"

(because `@tvectorize` has since been renamed on the latest releases)

Also, note that Julia starts with only a single thread by default. You'd either need to add `-t4` when starting julia, e.g. `julia -t4`, or set the environmental variable `JULIA_NUM_THREADS=4` for 4 threads, for example.

Looks good, but I actually love that Julia gets at least the same performance right out of the box, with full interoperability and the great multiple dispatch logic.

Would be interesting to run your optimized code against the Julia optimized code that was linked to in this thread. Or run a Julia GPU benchmark as well.

As for loops versus vectorization - I personally am used to vectorized code as it has always been a requirement. Sometimes it makes sense to use (e.g. Matrix Algebra stuff).

On the other hand, I have also been in many situations where I could not vectorize my code. More code than I'm happy to admit includes list or dict comprehensions or even loops. Sure, not all performance critical. Nevertheless, the prospect of performance no matter what is pretty exciting.

Finally, threading and multiprocessing in python is a headache compared to Julia imo. Sometimes I just miss the good ol' "parfor" from Matlab.

All in all, Julia is a really appealing value proposition to people doing numeric computing.

> the great multiple dispatch logic

Curiously, as much as I love the Julia language, I cannot stand multiple dispatch. I find it extremely confusing, ugly, and downright disturbing. Why would you ever want two different functions with the same name? And if this is good, why stop here? Why not two different variables with the same name but different types? Depending on the context where they are used, the language could pick one or the other.

But this is what most languages do, except they have _singular_ dispatch. Every OOP language has different behaviour for different input types, except it only works for the _first_ input argument. Julia simply extends this to every argument.

How would you feel about having to write `sqrt_int32(x)`, `sqrt_int64(x)`, `sqrt_float64(x)`, `sqrt_rational_complex_float16(x)`, etc, etc, etc, etc, instead of just `sqrt`, and let dispatch or overloading take care of the different implementations?

> Every OOP language has different behaviour for different input types

Yeah. And I find OOP horrific.

> How would you feel about having to write `sqrt_int32(x)`, `sqrt_int64(x)`, `sqrt_float64(x)`,

What kind of savage wants the square root of an integer? I'd prefer if the language only had a single number type (maybe configurable at once by an external option) and a single sqrt function.

But if you talk about the general problem of naming functions differently according to their types, I feel that it's perfectly OK. A tiny price that I'm eager to pay for the large benefit of being able to identify which function is called just by looking at its name (and not at the---possibly yet undetermined---types of its arguments).

But this is precisely what Julia allows you to do.

Your function may be written for numeric, and you use your favorite kind exclusively.

Now person B from this thread comes along, using both ints and floats. Well guess what, your package still works seamlessly. And where it doesn't, it must be due to a particularity that you yourself don't care about. In any case, the function can then be extended without namespace issues. In other languages, one would have to rewrite your package since you are using an obscure numeric type which others don't agree with.

Consider for example how people have used solvers for differential equations with entirely new datatypes, simply because the compatibility comes for free.

I think it's really neat!

I think it's great!

> I'd prefer if the language only had a single number type (maybe configurable at once by an external option) and a single sqrt function.

If you meant a numeric single type class, like haskell's Num, I think it's a great option. If you mean a single numeric type, like javascript, it unfortunately leads to a bunch of issues. Integers and floats really are both necessary very frequently. For example floats to represent something like speed, and ints to represent your bank balance. And you need different sizes of them (like float32 vs float64) are still necessary in a ton of applications, so you'll need them eventually if you want the language/numerics library to be truly general purpose.

It's also a huge convenience to be able to apply e.g. `+` to arrays as well as numbers and other polymorphism niceties.

> What kind of savage wants the square root of an integer?

It comes up moderately frequently (e.g., in prime sieves). It's usually defined as the largest integer X whose square X^2 is no greater than N. Search any sufficiently large codebase (postgres or something) and you'll probably find one or more functions called something like "isqrt".

My favorite implementation just takes Newton's method and blindly applies it to integers. It happens to converge for at least one starting value.

good point! but notice that this is an entirely different function than flotaing-point sqrt; and it has, appropriately, a different name.

The integer square root is mild a generalization of the floating point square root (in the sense that if you swap ints for floats you get the ordinary sqrt back out). With that in mind, do you think that the different name is chosen because it's a fundamentally different operation or just because of limitations in the dispatch mechanisms of common languages?

it would be very confusing to have both functions go by the same name. Very often you want the floating-point valued square root of an integer. Having the two functions with different names is much clearer, and better than forcing a conversion just to dispatch another function that performs a different computation.

> I'd prefer if the language only had a single number type

> But if you talk about the general problem of naming functions differently according to their types, I feel that it's perfectly OK.

You and I have nothing in common in this entire world.

Why on earth would you even want to touch Julia with 12 lightyear stick? The whole point of Julia is diametrically opposite to your preferences.

Being opposed to multiple dispatch, polymorphism and generic programming, while 'liking' Julia, is similar to loving bycicles, but abhorring wheels.

But I really love Julia! I was seduced by the idea of "like octave, but with fast loops", which is exactly what I need and the interpreter pretty much lives to that ambition. Also I love the examples with single-letter greek variables, the "f.(x)" notation, everything. And the fact that it feels like a well thought out language for numerical computing, cleaner and more beautiful than matlab/octave, and not an ugly kludge like numpy.

The multiple dispatch and oo features are some unfortunate warts, but I can live with those. Following your analogy, I love bikes and Julia is an electrically assisted bike. Sure, I would prefer if it was lighter and without the stupid motor, but it's still a bike. Not like numpy which is a horse carriage.

> You and I have nothing in common in this entire world.

for one, I 100% agree with your views on variable naming elsewhere on this thread ;)

> "like octave, but with fast loops"

But the performance relies on the aggressive specialization which depends on multiple dispatch. And the adaptation to numerical computing, the cleanness and beauty is all about multiple dispatch.

To stay with the bike analogy, multiple dispatch is definitely the wheels, not the motor.

> for one, I 100% agree with your views on variable naming elsewhere on this thread ;)

Waaah! (pulls hair.) And yet, so far apart on the function naming ;)

But seriously, though. Without the incredible polymorphism and genericity, there's really nothing at all left of Julia. Multiple dispatch isn't a feature bolted onto Julia. It is the core philosophy, and the central organizing principle.

> But the performance relies on the aggressive specialization which depends on multiple dispatch.

This is not a necessity, but an implementation choice. "Specialization" is an unnecessary step when everything is explicit from the beginning. I'm not talking on philosophical grounds, but thinking on concrete examples of jit systems which are really fast but have nothing to do with multiple dispatch (for example luajit).

> Without the incredible polymorphism and genericity, there's really nothing at all left of Julia.

If "matlab with fast loops" is "nothing" to you, sure. As a user of Julia, this is the killer feature for me. The rest I see as unnecessary complexity and mumbo-jumbo. But there's nothing wrong that each user has different favorite parts of the language!

Actually, come to think of it, Matlab is now fully jit-compiled. As long as you stick with double precision floats, there is already a "Matlab with fast loops": it's called Matlab!

Just stay away from classes and other complex data structures, which you don't like anyway, and performance is very good indeed.

Cannot reply to deeply nested posts, so I put the answer up here:

> This is not a necessity, but an implementation choice.

The performance of Julia relies on algorithmic specialization, and since you do not know what input types your code will receive, you cannot take advantage of clever specializations. I understand that you reject the concept of allowing user types that can be used by generic functions, but then you throw all the performance out the window. I want to do fast linear algebra, but your library says "no special static arrays!", "no special handling of diagonal matrices", and in general no exploitation of special structure in user types.

I want to calculate fast derivatives, but your library doesn't accept dual number types.

Your code will never be fast if it's only fast Float64, because it's when my special type elides all computation that the real speedups arrive.

> "Specialization" is an unnecessary step when everything is explicit from the beginning.

But you cannot know what it should be explicit about. I give your `sum` function a special Range vector, or a BitVector, or a OneHot vector. What should it do? It has to fall back to the same thing it does with all arrays: add-add-add-add-add. What is the norm of my UnitVector type? Your library has to ploddingly square, sum and sqrt. In the meantime, multiple dispatch just replaces the call with the number 1.0.

> If "matlab with fast loops" is "nothing" to you, sure.

As a daily user of Matlab for 25 years, several tens of thousands of hours of grudging use, I assure you it would be less than nothing. Matlab is a tragic, horrible mess of a language, and it's not even that slow! I rue every hour I've spent on it. If Julia were slow, and Matlab faster, I would still prefer Julia.

Also, Matlab has OOP and sigular dispatch.

> But there's nothing wrong that each user has different favorite parts of the language!

Well, as I am trying to say, multiple dispatch isn't a part of the language. It is the language. Everything revolves around that concept.

You can print out the assembly it generated with code_native [0] and check if it is seems there are superfluous instructions.

To really tune things you'd have to check the instruction timings among other things[1] for your processor and to be fair apply the equivalent setting for your compiler back end as well if it supports tuning for your processor and anything else applicable.

While it is likely there is a possibility for improvement, like most optimization work, it is usually for diminishing returns for your time and may be small in nature, but sometimes compilers, even those like llvm, do something strange, so you never know for sure unless you check.

[0] https://docs.julialang.org/en/v1/stdlib/InteractiveUtils/#In...

[1] https://www.agner.org/optimize/

I only glanced at what you did but if this operating matrices or vectors, there are often crazy optimizations that can be made for CPU cache

>if you want to fuse the loops and get gpu/tpu for free

the work being done at JAX is great and Julia shares many of the goals (Julia would even be the pioneer in some cases). The best part is your library doesn't even need numpy/JAX as dependency. Yet your function that works on AbstractArray will work on arrays that live on GPU/TPU, for free. thanks to multi-dispatch, you don't need to write a ton of boiler plate code for interface and inherent some classes from JAX/numpy.

It seems to break down a bunch in real use cases though. For example, we found a differential equation benchmarks the Jax devs were optimizing, where it was still 6x slower than Julia (https://gist.github.com/ChrisRackauckas/62a063f23cccf3a55a4a...). Then the gradient calculations don't mix well with GPUs, at least for ODEs (https://github.com/google/jax/issues/5006). That lack of composability compared to Julia is the main issue that then creeps up in the bigger problems.

For fun, I ported your 3rd solution to J

  f =: 3 : 0
  'k a' =: y
  n =: #a
  A =: (n, n) $ a
  ^ 0j1 * k * 0.5 ^~ (*: A) + (*: |: A)
  )

I'm curious how the speed compares to Numpy, but I don't have a python environment installed. It ran in 6.6 seconds on my computer for n=10,000, but it used a ridiculous 10 gigabytes of memory. When you consider a 10,000x10,000 matrix has 100M elements, that you need two of them since you add two together and that a short has 4 bytes, you technically should only need 800 MB for the whole thing.

Hi, I'm the post author. As others have pointed out, I'm not a Julia evangelist of any sort, I am just starting to use Julia and I'm very enthusiastic so far.

As for more readable Python alternatives, indeed, the ‘np.tile’ isn’t strictly required, as Python will add a (1,n) matrix and a (n,1) matrix into an (n,n) matrix by default (I would be happier with an error message here, though).

There are many alternative ways to take advantage of this fact, for example what you have suggested

  A = a[:, np.newaxis]
  M = np.exp(1j*k*np.sqrt(A**2 + A.T**2))

Note that A is a (n,1) matrix here, which is not quite obvious at first sight. Equivalently, my new preferred alternative is resorting to two reshapes

  n = len(a)
  M = np.exp(1j*k*sqrt(a.reshape(n,1)**2 + a.reshape(1,n)**2)

which is more explicit about the shapes of the arrays involved. I have updated the post to better reflect this discussion.

In some more complex manual-vectorization cases I have encountered, the np.tile cannot be dropped, and the code looks a lot like my original posting.

Being able to resort to loops, and not even having to think about this manual-vectorization issues is a big plus, if you need to do this very often in your code.

Regards

Hi! Thanks for engaging with the comments.

I believe you've misunderstood the main point of my comment. I don't have too strong opinions about reshape vs newaxis.

The real point was about excessive inlining, which your update does not to fix. A fairer comparison regarding readability would be pulling out the matrix version of a into a named variable:

  n = len(a)
  A = a.reshape(n,1)
  M = exp(1j * k * sqrt(A**2 + A.T**2))

I believe these are much much more readable.

When you need loops, you need loops, of course, and python tends to suck here (though with jax's functional programming constructs the boundary is shifting). But the readability/cleanliness comparison in your post is still an unfair comparison.

If you're going to go through the effort to split things out into separate lines and name them, why the single letter names?

In mathematically oriented code, long variable names make code less readable, not more. Mathematical expressions become long and difficult to parse, and the variables often have no intrinsic meaning, so trying to make up 'meaningful names' is pointless.

Some things are just very common in math notation, enough that people recognize it better than fully named variables. Things like uppercase letter for a matrix, _T for a transpose, etc.

For someone with maths training (be it years in undergraduate engineering or whatever), one does get quite used to one letter names. Note that they always should have limited scope.

It's funny, but I find your example way more illustrative and deserving to be shared on HN, than the original post. Honestly, I don't find microbenchmarks that compelling, because there are too many caveats to take them seriously without really understanding what's happening under the hood. And I wouldn't participate in arguing about manual vs auto-vectorization, because these are just 2 different use-cases, so there's no question for me that auto-vectorization being possible in the language is just better than it being impossible. (I mean, disregarding the programming, writing kind-of-imperative formula when discussing math can serve a purpose, as it is a different perspective, than writing the same thing declaratively. I would prefer "vectorized" declarative form in the most cases, but surely one can find exceptions, where M[i,j] = ... is just much clearer.)

But your comment just vividly shows, how much proficiency in writing clear code actually matters.

I think most of the criticisms against Julia come from the fact that people see it as a Python competitor. In its actual state (even if the last version improved compile time a bit), it is not really usable as an exact python replacement (short scripts, opening a new interpreter each time).

However it is a revolution for the HPC crowd. If your program is supposed to run for several days on multiple nodes, 10 seconds of compilation is a negligible price to pay for a much cleaner and shorter code than Fortran or C/C++.

I've been using Julia for a few months for physics related simulations and it's just so much better than Fortran to work with.

While I definitely agree, I'd say Julia is totally fine for short scripts too. Just spin up https://github.com/dmolina/DaemonMode.jl and you can launch lots of short scripts without excessively restarting julia.

There's also PackageCompiler.jl https://github.com/JuliaLang/PackageCompiler.jl for AOT compilation, but that's heavy enough that I wouldn't bother with it for short scripts.

Of course, these packages are fantastic, they do help a lot for more experienced users. But for a pure beginner, there is a high risk he will give up Julia after 10s of Plots.jl compilation, way before we would hear about these.

(and tbh, even leaving a REPL open all the time works perfectly fine, since include("myscript.jl") really does reload the script, unlike the import keyword in Python)

Yeah, that's fair. I guess I read your comment as saying that julia isn't a good replacement for Python, rather than it takes a tiny bit of learning when you're coming from other languages.

I think the way we talk about julia sometimes makes Python and Matlab programmers think they can just take some of their old code and copy paste it into julia, switch around some keywords and have everything be faster. It's important to emphasize to these people that it's a different language with its own idioms that need to be learned.

I strongly agree with that, coming from Matlab to both Python and Julia. While I enjoy some of the similarities between Julia and Matlab, it does take some learning and effort to write proper code in another language. Array indexing with square brackets in Julia tripped me up initially for instance as Matlab uses parentheses, and it still take non-zero mental effort to adjust when writing Julia (despite having done C in a distant past...).

Writing proper and idiomatic code that gets the best out of each language takes even longer as it requires more knowledge and experience. Having spent the last ten years or so predominantly in Matlab, and the corporate environment moving increasingly to Python, and my own interest in Julia, I am getting a double dose of this.

Pragmatically I like Matlab best, in large part because I am so comfortable in my workflow and the large existing codebase, but also the IDE, debugger, etc. I most fascinated by Julia but find that exploiting its potential has its own learning curve, and wrangling with type stability has its own challenges. I am least enthused by Python, which I am learning mostly from necessity, but this may be colored by my extreme aversion to its use of indentation.

I also started in Matlab + C, and didn't personally like Python very much, but wound up liking Julia after a bit of adjustment.

Still a work in progress, but I've written up a few of my lessons learned in the process [1] in case they're useful to anyone else. Properly figuring out type stability and dispatch-oriented programming took me way too long, but things started making a lot more sense after that.

I haven't used it a ton, but [2] also seems potentially useful for anyone else interested in making the same switch.

[1] https://github.com/brenhinkeller/JuliaAdviceForMatlabProgram...

[2] https://lakras.github.io/matlab-to-julia/

Python has importlib.reload(), though I haven't found much need to use it

Indeed. Given the fast startup time of Python, I've always just restarted the session (which is most of the time faster than typing "import importlib; importlib.reload()" haha)

I really want to pick up Julia as my general purpose scripting language, but neither of those solutions meet the standard for “totally fine” imo. For reference, my current use case for Python is 95% scripting, 5% ad hoc exploratory data science.

I need a way to package and distribute scripts, so that they can run in a container for example. Running a server to execute scripts only works well for local workflows. PackageCompiler.jl has terrible ergonomics. I don’t think Julia is quite there for scripting yet, unless you’re already deeply invested into it for other reasons.

I recently (last week) started using numba, for similar reasons to why the author seems to like Julia. I tested translating his example to numba:

  @numba.njit(parallel=True, fastmath=True)
  def w(M, a):
      n = len(a)
      for i in numba.prange(n):
          for j in range(n):
              M[i,j] = np.exp(1j*k * np.sqrt(a[i]**2 + a[j]**2))

and timed it like this:

  %%timeit
  n = len(a)
  M = np.zeros((n,n), dtype=complex)
  w(M, a)

On my 8-core system, this ends up more than 10x as fast as the numpy version he listed (which seems to lack the sqrt, though), which would place it close to the multithreaded Julia, even considering that ran it on a 4-core system. As an added bonus, it can also pretty much automatically translate to GPU as well.

the raising of numba shows us why numpy and "just write vectorize-styled code with a C++ backend" is not enough.

yet Numba basically makes your python code not python. It doesn't support so many things: pandas dataframe, or even as simple as a dict(), which means you often have to manually feed your numba function separate arguments.

To separate a complicated calculation into numba-infer-able parts and the not ones is not fun and sometimes just impossible.

Yep, completely agree. For the project I'm currently doing, it seems like a fairly good fit though. Lots of prototyping different approximations, and needs to be faster than plain numpy.

Also, the jitclass things help somewhat. I use them as plain data containers, to work around the hideously long argument lists that otherwise would be required, but with no methods. jitclass breaks the GPU option, though.

I have personally experimented quite a lot with numba. When it works, it's great. However, numba can have very cryptic error messages making it difficult to debug.

Which is why I switched to dask, which even though slower integrates better with numpy.

Actually, Numba does support dicts now.(You can't have a mix of types in the dicts unless that's changed, but that isn't an actual problem for most ML work.) I have used numba very effectively to make my machine learning research projects run very quickly. I don't use pandas; I do use a lot of numpy and scipy. I understand that pandas can use numpy arrays for at least some things. Since numba works great with numpy, it seems like that might be an approach for using it with pandas, in at least some cases.

This example is not really testing Fortran vs. Julia, but rather testing the speed of the math library you linked to.

The most expensive part of the fortran code is evaluating sqrt and exp.

For example on a 5600X the current Fortran code (with n=10 000) takes about 2.7 seconds. If I change the definition of M to remove the special functions (M(i,j) = a(i)+ a(j) then it takes .3 seconds. (Using gfortran with -Ofast)

When the special functions are the issue, using the Intel Fortran compiler and Intel MKL will make the code about 2x faster (this is what Matlab or Julia is probably using)

`vsqrtpd` is an assembly instruction. `exp` and `sincos` are implemented in Julia.

If you want to use a library to speed up the Julia code, checkout https://discourse.julialang.org/t/i-just-decided-to-migrate-... (but replace `@tvectorize` with `@tturbo`), which is code by eigenspace.

My point is that for codes like this, the issue is not vectorization.

Just evaluating the special functions exp and sqrt is taking 90% of the time. You can evaluate these special functions in parallel, but that is about it.

Julia ships with OpenBLAS, in some cases there are pure-Julia "blas-like" routine that can be as fast:

https://github.com/mcabbott/Tullio.jl

Or much faster (at small sizes): https://github.com/JuliaLinearAlgebra/Octavian.jl

Also much faster at very large sizes :)

Just another data point. I've just tested pythran (a numeric python to C++ compiler) with the code from https://discourse.julialang.org/t/i-just-decided-to-migrate-...

On my machine the julia code on my PC is:

  N = 2000
    88.273 ms (8 allocations: 61.04 MiB)
    49.445 ms (3 allocations: 61.05 MiB)
    14.599 ms (2 allocations: 61.04 MiB)

If I compile the following code with pythran:

   #pythran export texp2(int)
   def texp2(N):
       a = np.linspace(0,2*np.pi,N)
       return np.exp(1j*100*(a[:,np.newaxis]**2+a[np.newaxis,:]**2))

I get:

  timeit texp2(2000)
  12.3 ms ± 27.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

So actually slightly faster than the optimized julia code.

How does Pythran fare against the code in this link:

https://discourse.julialang.org/t/i-just-decided-to-migrate-...

The author made some julia performance mistakes in his benchmarks.

The code in this version:

https://discourse.julialang.org/t/i-just-decided-to-migrate-...

should be very tricky to beat, and I suspect that the only way is doing a lot of customized assembly.

That is the code I'm comparing against. So the fastest Julia version is actually the one using @tturbo in your link. I'm not surprised, pythran gets me the best performance across numba, customized cython and (for some comparisons) julia in most of my tests.

Ah very interesting, this is great to know! Sorry I didn't properly read your initial comment. I'm now quite interested to see what we can do to match or beat Pythran here.

What CPU are you using? I know my CPU, a Zen+ Ryzen 5 2600, is one of the worst modern CPUs for this sort of thing and there's much better scaling on Intel CPUs or the Zen 3 CPUs, but those things are probably also effecting the Pythran code as well.

Is Pythran any good at producing BLAS microkernels?

Would be good to confirm the number of threads being used by both, but it could also be that pythran is using a better `sincos` implementation. On very old hardware (which it would have to be if multithreaded), the exp implementation might not be great either. That'd especially be the case if it has no AVX, as the exp implementation relies on a table (which requires gather instructions).

For 2000x2000, their time was 14ms, vs 6ms for the 2600.

Yeah, I tried to get Pythran working on my computer to compare against, but it wouldn't build for me.

Just an additional comment, using sincos in the comparison is not really fair, because it is a change to the algorithm (it assumes that a is real), because in general

   exp(x) = exp(x.real())+ sincos(x.imag())

(for simplification sincos(x) returns cos(x) + jsin(x), but you get the point).

The other code can't easily optimise for that (although it would be a nice optimisation). Interestingly (and I can't explain why, maybe someone knowing more about julia), if I change the code of that third method using tturbo to use the cis function (or exp) it becomes a factor of 20 slower and is actually the slowest of all the different algorithms)

  exp(x) = exp(x.real()) * sincos(x.imag())

Anyway, it is the same algorithm as long as `x` is real. The reason for manually inlining exp(::Complex) and manually decomposing the complex number into its real and imaginary parts is that `@tturbo` currently only supports real inputs.

Given types that it does not support, it will silently run the original loop single-threaded and optimized simply with `@inbounds @fastmath`. I imagine in changing to cis or exp, you also switched to using complex numbers? If so, that would explain the dramatic slowdown.

If not, mind sharing your example, here or in an issue at LoopVectorization.jl?

It's planned that "some day" LoopVectorization will support complex numbers through the AbstractInterpreter interface, but this is probably still a ways off.

I agree, that it's the same algorithm as long as it x is real. My point was that with this optimisation we give more information to this function than the comparable python function (or the other julia implementations), so it gains an advantage for optimising.

Yes I was using complex inputs when switching to to cis/exp. LoopVectorization not supporting complex does explain this.

As a side note, I generally dislike silent fall-backs, for optimisation libraries/modules. That's why I always do python=False when using numba (and I wish it was default). I'm using the library/module because I want speed up, the problem with the fall-backs is they are often slower than the default code (e.g. unwrapped loops in numba, and in this example), so I would like to know if it doesn't speed up the code, because rather not use it then. A failure would save me benchmarking the slow-down.

Here’s the code for the fully complex case: https://discourse.julialang.org/t/i-just-decided-to-migrate-...

Hopefully after the LoopVectoruzation.jl rewrite we can make it understand complex numbers automatically. We could probably put in special-case complex support today, but we’d really like to support arbitrary structs of bits.

Did you start Julia with multiple threads? Julia unfortunately starts with only a single thread by default, and the slow time you reported (14 ms) makes it appear that this may have been the case.

You can start `julia -t4` for 4 threads, for example. Or you can set the environmental variable `JULIA_NUM_THREADS=4`.

If you were using the same number of threads for both, mind letting me know which CPU you're using?

You were correct, I just tested with julia -t4 and that improved times by about a factor of 3.

  N = 2000
  27.137 ms (24 allocations: 61.04 MiB)
  49.496 ms (3 allocations: 61.05 MiB)
  4.597 ms (2 allocations: 61.04 MiB)

The CPU is a ryzen 3600 btw.

-t6 would give 6 threads. FWIW, the Ryzen 2600 with 6 threads took 6 ms, so the 3600's floating point and vector improvement is quite evident from your timing. (The 2600 emulates 256-bit operations with 2x 128-bit, while the 3600 has full support.)

Yes, I've also run with 6 threads and get an additional minor improvement on the loopvectorized code (3.9 vs 4.5ms).

I have to say LoopVectorization.jl is quite an impressive piece of work. I have to try out some more julia at some point, unfortunately pretty much everything I do is using complex values.

Can you post the C++ that was generated somewhere?

It's heavy templated code, which I don't think is very useful without the pythran source, but here it is:

https://gist.github.com/cycomanic/a0b16ab0c45b3d888a818d17a0...

I just wish the Julia community wasn't so invested in using UTF8 symbols in the code. I understand that for Physicists, Engineers, Mathematicians, etc it look cleaner and more familiar, but when writing code that is meant to be maintained and read over and over agains it makes a difference if that µ looks too much like a u sometimes. What's so bad about just writing mu?

Anyway, this is not a complaint about the language (which I like very much), just a dislike of the popular usage

The right name for a variable is related to what you think the way people will check it will be. In most code the attestation to the correctness of the code is the programmer's understanding of the code. Clearly then names should be descriptive. In contrast, when the attestation to the correctness of a formula is that it is the same as a formula printed in a published paper, then it should look as much like the formula in the paper as possible. A programmer who does not know that rho means density will not be able to catch scientific misconceptions by virtue of the fact that they see the word "density." On the other hand, the will be able to catch differences between the symbols in the source material and the source code.

In any code you write, you should be asking, "how will people check this?" If the answer is "by comparing it to something else," then it should resemble the original as much as possible. If the answer is "by thinking about it," then it should be formatted for self-contained comprehensibility. The distinct use cases lead to different best practices.

I wish non-Julia community would adopt this practice, because unicode symbols exist for a reason, and using proper symbols instead of making look every language like Perl with unix-style abbreviations that look bogus to pretty much everyone but the most seasoned C-programmers. And complaining u and μ are too similar is like complaining 1 and l or 0 and O are similar. I mean, yes, they are, unfortunately, but if you are noticing it's a real problem for you when writing code, you should probably check out better fonts for your text editor.

For the record, I'm not a blackboard addict, I never had a career in academia and I'm not even really a member of Julia community any more than any other person who wrote maybe a couple thousands of LoC in Julia. But I wish people would realize it's not 1970's anymore, and ASCII is culturally outdated. For that matter, I would love if I could write formulas in my code the same way I do on paper, like in that MathJax sample from the article. That's why math notation exists at all, imperfect as it may be (and I hate many things about established math notation, but it's the best we have). I just don't see how it's possible without abandoning plain-text, which I surely wouldn't like because we (currently) don't have tools that would make handling it as effortless, as editing text in Vim.

For that to happens common keyboards need to include these esoteric symbols. Yes, esoteric. Who's in the market for a keyboard with Greek alts? Or do you prefer a space cadet keyboard? What if you already have a US-Arabic keyboard which is fairly common?

The fact is it's easier to type "mu" instead of "u" (like right now, I can't be bothered to look up "mu unicode", or open the special characters picker, or exhaustively brute force alt-gr key combos to find the mu).

Take it a step further, I wish physicists and mathematicians would stop naming their variables so poorly so that only the most seasoned of their field can follow along.

Sometimes, like v (velocity) it's quite obvious, at other times I find myself staring at scripts with insane variables like a,o,t_o, etc. And this is a habbit that seems to carry over to non-academic code as well. It's trivial to spot code written by a mathematician because it tends to be incredibly obfuscated.

Do you really need to name your constants K_(single letter)? Or can you just write "AIRDENSITY" once and let your editor autocomplete it in the future?

The reason I complain about these symbols is because they don't have clearly defined meaning. Maths/Physics had to pull out the Greek alphabet after it exhausted the Latin alphabet (and overloaded most letters with 2 or more meanings) and then promptly molested the Greek one in the same way. Just look at Omega. It's useless. It has been overloaded with so many possible meanings that it has none. Stop naming your variables "omega". Unlike Maths, programming languages support variable names of more than 1 character long.

Long variable names make mathematically oriented code unreadable. A short half-liner becomes a multi-line monstrosity. What you need to understand is that in mathematical code, the variables themselves are of secondary importance, while the structure of the expression is primary. With long names, you cannot see what's going on.

Julia allows you to input Unicode characters via latex names though. I think that's a sweet spot of ease of typing and expressivity of characters.

Honestly I really hate seeing mu_a, delta_k and company all over the place. If we are gonna name our variables like this might as well give us the actual characters

> making look every language like Perl with unix-style abbreviations

your criticism of Perl is unfair. Base perl has supported unicode variable names since more than 20 years ago, way before all the other languages that you speak about.

Physicists, mathematicians, and engineers have been writing µ and u on chalkboards for years with no issues. Why should it be harder on a screen?

The difference is that form is important in mathematics and not very important in programming. Therefore symbols are used because it makes it easier for the eye to recognize patterns.

For example these equations look similar,

x^2+4x+5=0

y"+4y'+5=0

Even though one is a polynomial equation and the other is a differential equation, the common visual pattern suggests that the techniques needed to solve them may be similar.

On the other hand in programming there's no need to do this, all the math has already been worked out on paper, so it's better to use clear, distinct, easy to type variable names.

> form is important in mathematics and not very important in programming

How important is form when I'm programming mathematics?

> How important is form when I'm programming mathematics?

Not important.

As a physicist who does a lot of mathematical programming for my every day job, I'm going to disagree pretty hard.

With how I've defined form, I'm not sure how you could.

Who in their right mind does math on a computer, math is implemented on a computer, but it's done by hand.

> Who in their right mind does math on a computer, math is implemented on a computer, but it's done by hand.

Maybe you think this becuase your used to writing math on computers that make writing that math ugly and obtuse?

Besides, having the math you've written by hand closely match the syntax you use on the computer can greatly reduce the cognitive burden of switching between code and paper.

I'll just make one brief comment in regard to wanting the math you've written by hand to closely match the syntax on the computer, and say that this is rarely the correct thing to do when writing high quality numerical software.

Making your code performant, and doing a thorough numerical analysis, usually requires a significant rearrangement of your formulae.

I'm sure there is so much room for improvement here. The closest thing for me is using TeXmacs: https://m.youtube.com/watch?v=H46ON2FB30U

Usually, the formulas I derive on paper end up quite complicated and can not be simplified away. When I am numerically implementing the formula, I want that it matches closely what I wrote in the paper, and long variable names would make my code less clear mathematically.

Why would I use the variable name "probability_distribution_on_m" instead of "rho_m", when I know all along what rho_m means in the contexte of my code (the symbol use throughout the paper). Usually, I comment at the beginning and specify what the variables mean and to what they correspond in the paper. And if somebody needs to read my code, that person will need to understand the paper first. More descriptive variable names won't make the person understand the paper better.

Of course, when I am writing some non-scientific software, I will use descriptive variable names, because there is no complicated formula, no paper that sets the context of the code, and it makes sense in general to have descriptive variables for logic elements for a clear code.

I didn't say to make your variable names verbose. rho_m is fine in my opinion, although I find it strange that you would name a probability distribution rho unless it's a density matrix. I would probably name it prob_dist_m, or even pdist_m if you want to be more terse.

What I do think is bad form would be naming it ρ_m.

How exactly do you type mu? On a chalkboard writing the mu sign is trivial, on a computer... well I don't have a mu button.

Edit:

There are like 5 replies to this mentioning different ways of doing it (including memorizing 3 digit unicode values) none of which seem more intuitive than writing "mu".

I just do \mu (or ;mu) in Emacs.

We're in the 21st century, all the tools we use should have reasonable Unicode support. The fact that we collectively keep on talking about typing non-ASCII characters as a real problem is a pretty depressing reflection of our shared computer infrastructure :(.

For me: alt+enter (switches to Greek keyboard layout), m, alt+enter (back to usual layout). I use Greek symbols in Python and in Javascript. I do it almost without thinking now, as lots of my coding is tightly linked to maths.

you can map a useless key like capslock to do that in a single step. Thus to type μ, for example, you simply press CAPS+m, and similarly for all greek letters

What is a capslock key, I only have ctrl and hyper ;)

You probably have what is called a GUI key, hyper died out with the space cadet keyboard. Caps lock is above the left shift.

man, it's the only key that has an assigned light in my keyboard! Funnily, the capslock light in my laptop seems to be detached to the configuration of that key, so that the light shows the parity of the total number of greek characters that I have typed.

> none of which seem more intuitive than writing "mu"

Write for the convenience of the reader, not of the writer!

I use on Linux the XCompose key (which is assigned the menu key right from the space bar, or the "Print Screen" key), and then "\mu" like in LaTeX. I don't need to think about it because I have used LaTeX for years. I also use the XCompose to generate € and ß which are special symbols in my native language. This has the advantage that I can use the UK International layout for everything, which makes specifically writing code much more comfortable. Also, the xcompose keymap is just a text file in my home directory, I can extend it at any time and take it with me. The method is also independent from any application.

Also, Emacs has an input method which supports "\mu", and this convention is also used in Racket for things like λ.

Vim/Neovim makes it very easy to input all kinds of interesting Unicode, using the "digraph" system. A lot of other code editors have similar functionality, either through keystrokes or some kind of "palette" interface. And Julia REPL itself has Tex-like escape sequences.

I agree that we need better system-wide tools for arbitrary Unicode input. Font support and confusable glyphs are also issues. But I think these are solvable problems, and they will be good problems to have solved.

I type <Greek>-m. One extra keypress, just like when I type é, ü, č, ç, π, etc. If you type in more than one language, or math, your keyboard is already set up to handle this.

in Julia:

\mu<tab>

On Windows you can use Autohotkey and map it to "Right Windows Key" + u. Or use "Right Alt Key"+ "Numpad 230" ALT Code, https://usefulshortcuts.com/alt-codes/greek-alt-codes.php)

In the Julia REPL and most Julia editor extensions you can just use \mu <tab> and call it a day. This is a non-issue for Julia users.

It's not the screen so much as the keyboard that is the problem

They’re very easy to type in Julia’s repl, or any text editor with Julia support.

You just type \mu and then hit the tab button and you get μ. I actually switch over to my Julia REPL all the time when I’m emailing someone and want to type a greek character.

> with Julia support.

So you can see how this becomes a bit of a barrier to entry.

I'm not sure I'd want to do any serious coding in an editor that doesn't have the most basic support for the language I'm using, but to each their own I guess.

For new users, if they're not used to writing unicode characters and don't have an editor with nice support, then they don't have to write them.

The julia community has a pretty strong convention of giving both unicode and ascii versions of almost all functions. Some packages don't have this convention, but I'm not aware of any big important ones that don't follow it.

Anyone that doesn't use a language aware editor is making themselves a disservice.

that's not really acceptable. I'm a big julia fanboy, have been since 0.4 (proof in bio) but dismissing real pain felt by newcomers is kind of bad for the ecosystem.

What editor do you use, if I may ask, and have you asked them to add unicode support? Part of the goal in permitting unicode variable names is that more and more niche editors might be forced to become fluent in languages other than english, benefits everyone (and might be a bit overdue at this point).

I use vim, and have a completions plugin (part of julia-vim mode), so that I can write \mu<tab> for any unicode name, which makes it as easy to type as "mu", but easier for me to read (as it matches my textbook expectation of the formula).

It's not a display issue. Do I want to type ctrl-u <whatever hex> or copy-paste if I want to write some code that calls a function when I'm playing around with a new library? Well those are the two strategies I've used with experience... But some person coming new to Julia will grab some package, and then see some function call and get deer in the headlights "how the hell do i do this".

And LaTeX codes are not always obvious. For example it's \triangleleft but \leftarrow

that's two extra strokes than just "mu"

If Emacs doesn't have enough input methods, you can add one reasonably easily... A good deal of early-ish work on multilingual support in Emacs was done out of interest in editing technical text, which is basically the same problem as more-or-less complex natural language scripts (which deserve support anyway, of course).

Code is read more often than it's written, even scientific computing code.

agree, that doesn't make it not a problem.

It is a non-issue when you can use just use autocompletion for mathematical symbols on almost all Julia-supported editors

Besides the good points about the difficulty typing characters, not all fonts have good glyphs for non-alphanumeric (by that I mean low ASCII) characters. Sometimes they don't have those glyphs at all and the system has to go to a fallback font. Terminals can also shit the bed with Unicode characters.

The solutions to those problems is to use better fonts, terminals, and editors. Using Unicode characters is fine but some people do have legitimate problems with them.

On that topic, I'll throw in a shoutout to Cormullion's excellent JuliaMono font [1,2]!

Typing certainly does add more friction between u and µ than a chalkboard would, though it is perhaps also notable that mathematicians seem to have felt that it's worth the effort, and have come up with a quite nice system for it in LaTeX -- and I would certainly not mind seeing more languages and editors add support for LaTeX completions the way the Julia ecosystem has.

[1] https://juliamono.netlify.app/

[2] https://github.com/cormullion/juliamono

> Besides the good points about the difficulty typing characters, not all fonts have good glyphs for non-alphanumeric (by that I mean low ASCII) characters

Not all fonts are good fonts for an evey use use; heck, not all fonts make 0 and O or 1 and l and I clearly distinguishable.

> The solutions to those problems is to use better fonts, terminals, and editors

Exactly.

Because of keyboards?

Because they dont do it at all on a screen. All papers written in Latex use \mu .

The thing is, symbols like µ and λ often have well-understood meanings in context, especially when using a particular package. It would be like a programmer complaining about using `i` to iterate through an array or `row` and `col` to iterate through a 2D array.

I'm personally very empathetic to this. Fortunately, I think most "real" developers are pretty restraint regarding this in popular packages' code.

Some proper use casees IMHO are: 1. in "terminal" code: scripts, notebooks. 2. function internal variables 3. for making the code look like their counterpart of a paper.

idk what to think of 3. if you look at any paper, non of them only uses ASCII, which raises the question, if we're happy with reading papers (even CS ones) with symbols, why not in our code?

I definitely agree with this. It's just a pain to type out.

Maybe I'm behind on this, but needing to use the mouse to find symbols in a huge menu is very painful.

You can copy and paste too but again, huge break of flow.

Do some editors support discord-style emoji syntax, where typing :fo would bring up a menu of emojis that might match foo. Then hitting enter inserts the emoji over the :foo: representation. You can also not use the auto complete menu.

Ex. :pow2: might turn into ²

In the Julia REPL you can type LaTeX commands for Greek letters and math symbols, and a <TAB> converts them to the actual character; so \mu<TAB> gets you a μ. Or you can just input the Unicode directly. I’m constantly doing such things everywhere, so the extra few keystrokes seem worth it to get code where the math looks a more like math. Using the mouse and a menu or character palette is the worst possible solution.

EDIT: In general it makes sense to set up your keyboard with a compose key and a “dead Greek” key, and set up your Compose table so the shortcuts make sense to you. Then you can use an expanded set of symbols everywhere, including comment boxes like this one. You can even put things like your email address in your Compose table.

Julia has the best REPL for this. You can \<type><tab> to type unicode symbol and emoji.

Better yet, you can reverse look up how to type a thing:

  help?> χ²ᵢ
  "χ²ᵢ" can be typed by \chi<tab>\^2<tab>\_i<tab>

Yeah. Now go edit your code.

I actually agree with this sentiment, using non-keyboard utf8 characters in code is just asking for massive problems down the line. A fun party trick perhaps, but such a bad idea for anything that isn't a hobby Project.

I disagree, in that it's the opposite. You might want to put some internal calls in Unicode as a way to gatekeep contributors and users to people who are serious about using Julia. But if you think your project could be used by casual explorers, your public API better not have Unicode calls, or, at least provide sensible non-unicode aliases for unicode calls.

I can edit it in Vim with no problem.

You can easily search for omega and lambda via vim slash(/) search command?

Yes, just like searching for any other character.

How do you type it?

On my keyboard, if I’m in Vim and want to search for λ, I type /λ. I type the λ with <right-CTRL>l.

Any special setup for that? I'd be hard pressed to touch type capital/lowercase sigma, ro, pi and quite a few others, I think.

I fully think this is a problem with archaic input - but that doesn't mean it would not be inconvenient when navigating through a code base :/

I just use what X-Windows supplies for customizing the keyboard. The way to tell the machine what modifier keys to use for what is with xmodmap. My .Xmodmap file contains:

    clear control
    clear mod2
    clear mod4
    keycode 134 = dead_greek dead_greek dead_greek dead_greek
    keycode 108 = Multi_key
    remove Lock = Caps_Lock
    keycode 66 = Control_L
    add control = Control_L

But you have to figure out the keycodes that your keyboard is sending to make it work for your particular machine. The program xev is good for this.

This mapping also makes the capslock key into another control key, which is good if you use Vim.

The Greek letters will be on keys that make sense, and the other Unicode characters are pretty rational too: to make é, type <mod>'e. You can make any key combination produce anything, even long strings, by putting the mappings in the Compose file. My Compose file is in /usr/share/X11/locale/en_US.UTF-8. I have snippets in there, such as my email address. This way the mappings work anywhere, not just in a particular editor.

Thank you. Looks like custom input is still a little up in the air on Wayland. I guess one benefit is that some solutions will also work for the console (though not for the bootloader, so no "pi" boot options... but i can live with that...).

https://unix.stackexchange.com/questions/292868/how-to-custo...

In particular: https://gitlab.com/interception/linux/tools

https://gitlab.com/interception/linux/plugins/caps2esc

(for remapping caps lock to escape - in the terminal and under Wayland).

Oh, sorry, I don’t know anything about Wayland. It took me long enough to learn how to get X11 to do what I want, so I’m not about to learn a new system unless there is a big benefit. I don’t use multiple monitors or do anything that X11 can’t handle, so I'll stay here in 1987.

Sure, because you have a julia plugin. Not every casual user trying out the language will bother, or want to bother installing someone else's plugin. They'll just hit a function call they need to make with a Unicode character, go WTF man, and give up.

> They'll just hit a function call they need to make with a Unicode character

This is optional, the language adapts to the capability of the user. He can write ∑ or sum for a function, µ or mu for an identifier.

https://docs.julialang.org/en/v1/manual/functions/

You get people riled up in this thread with your false assumption, next time simply ask questions.

that's for stuff in the stdlib. You can't control library writers, and julia does not force a function identifier that is mapped to unicode to have a "non-unicode equivalent", though, that might be interesting.

If libraries use unicode in their internals, it doesn't affect you or me. And I have yet to see in the wild any examples of unicode identifiers that are required for use. The closest thing i can think of is the composition operator in Base, though it has a somewhat awkward ascii equivalent.

Do you have any example? As far as I know this is strongly discouraged in the community.

I don’t have a Julia plugin. I don’t even understand what you think the problem is.

Oh dang I didn't know this. Thanks, TIL!

Yeah. Many editors (VS-code, julia repl, atom), support something very similar where I can type in `\ome` tab complete to `\omega` tab complete to ω. (similar completions exist for many other common unicode symbols).

That might be a problem with your tooling, but most likely not. I know how to do that in Linux (XCompose) and I've heard there are similar tools in both Windows and MacOS as well. So more likely just a problem of not being aware/proficient with using these tools properly. Typing λ² is 6 keystrokes for me, no touching mouse required. And while maybe not very illustrative in this post, I'd argue it's way more readable than lambda_squared when being a part of long and complex formula.

Absolutely it's my lack of experience in tooling. Never had to type unicode symbols consistently.

If your editor doesn’t help, on macOS ^⌘Space brings up the character picker with the search bar focused for the keyboard, so you can just search the name of the symbol and enter.

I think restraint and common sense are key here (once you're comfortable typing the non-ASCII characters). In parts of my code it's much more intuitive to use \Omega or \Sigma and \mu than giving them ASCII names, and others looking at my code will immediately understand what the matrices and vectors do.

On the other hand, if you're going to use \mathfrak{i} to index a simple loop you're doing it wrong.

In certain domains, people will immediately recognize the symbols for delta, lambda, mu...etc. That's how they're written in the textbooks and academic white papers and so on. Seeing the English word is actually more confusing.

> but when writing code that is meant to be maintained and read over and over agains it makes a difference if that µ looks too much like a u sometimes.

So, if you are maintaining code where that is the idiom, use a font where that isn’t an issue.

Mathematica has many non-latin symbols as well, and it works great. The symbols in my paper are the same symbols in my code. No confusion at all.

Have you had to use Symbolize to make compound symbols?

I have not had the need to do it. I have written two large programs, with variables as a[i,j] or as a_{ij} and in both cases Mathematica treats the whole thing as a single symbol.

You can also put Greek symbols in superscripts or subscripts and that works just fine.

Do you think the Symbolize (or the Notation package) are useful? I don't see any examples in the docs, so hard to tell.

I'm a mechanical engineer by training, and a software and control systems engineer by trade. IMO, neither alternative is a good one. I prefer to write out the un-abbreviated form wherever practicable. ex: "microradians", "microseconds", "permeability", etc.

And slower (to start) than both of them. I tried julia to write some scripts that were previously in numpy, but the slow start threw me off. I hope they solve this because it's a great language.

I have also recently been playing around with Julia coming from numpy/pandas. I found that 1.6.1 felt way smoother on startup and import than 1.5.x. Might be worth trying to upgrade if you haven't already. At the very least it's encouraging that it seems like they're making progress in that direction.

Unfortunately Julia is a “repl-based language”. You are expected to open a repl and work in it for as long as possible so that your compilations all remain cached. Packages like Revise can help reevaluate code that changes so that you don’t have to leave the repl.

I do wish that Julia would start up in an interpreted mode and compile in the background so that it would be fast enough when first opened and then attain maximum speed later on. (I think this is how JavaScript engines work?)

You can have what you want here by using DaemonMode.jl

https://github.com/dmolina/DaemonMode.jl

It would be nice if their was an "offline version" of this, e.g. automagically cache all compiled code to disk. Daemon mode isn't restart-friendly.

That's basically what PackageCompiler does. That said, StaticCompiler.jl is expected to improve the state of this in the future.

Not quite. PackageCompiler isn't automagic. The sysimage doesn't get automatically updated in the REPL. You have to create one manually (and this requires tracing instead of static type inference which Julia can do). So the workflow isn't zero-friction. And StaticCompiler hasn't been updated in over a year.

The problem with Julia tooling in general is that they feel 90% done. And tooling is, in my opinion, more important than the language itself.

> So the workflow isn't zero-friction. And StaticCompiler hasn't been updated in over a year.

There's much more recent work here: https://github.com/tshort/StaticCompiler.jl/pull/46 and apparently some more is still ongoing privately.

That said, yeah. I would not recommend Julia currently for people who truly believe they need AOT compilation, and that they need to trigger AOT compilation very often and with low friction. That definitely needs more work, but it's happening.

That said, a lot of people overestimate how much they actually need AOT compilation.

Julia has some fantastic tooling in other areas though. Especially the package management and interactive analysis tools.

I really wish there was an option for better non-REPL usage - basically Cargo for package installation and ahead-of-time compilation, but with the Julia syntax and ethos.

> (I think this is how JavaScript engines work?)

Yes, modern V8 does this, as does the JVM. Common Lisp runtimes also often support this, though I think they usually leave it up to the programmer to choose when to compile a function, they don't always do it automatically. The .NET CLR also behaves like Julia - JIT compile on first execution of every function.

The startup times have improved a lot in version 1.6, the latest release, and are I Proving again in the next upcoming release candidate.

Are you hitting the startup time very often? If so, you might want to try some things to keep one julia session open and sending code to it like a daemon instead of constantly closing and opening sessions.

This package makes that workflow really easy: https://github.com/dmolina/DaemonMode.jl

I most often need to restart the REPL when redefining structs.

Same for me. Ported my calculations to Python from Julia. Even before Julia has started python was done. And Julia also did jit compilation in runtime. A function call could take seconds on first call.

if it's that fast (<O(1s)) and your use-case is to call it many many times from ground-up, you sure should use Python or even bash since there's literally no performance to speak of. (i.e. nothing really changes the usability anyways)

Julia startup can easily be minutes with big libraries

`precompile` is a one-time thing.

Actually loading them when running scripts is O(s) even for some of the biggest library (Plotting, Differential equations etc.)

I thought Julia is a compiled language. Surprising that the startup is then so slow.

JIT compiled. Has a warmup time.

It’s just-ahead-of-time compiled and no «interpreted mode» so it’s always slow to start, if I understand correctly. Really annoying, especially the first time you use it.

If you need to run small scripts and can't switch to a persistent-REPL-based workflow, you might consider starting Julia with the `--compile=min` option. If you know which packages you're going to need, you can also reduce startup times dramatically by building a sysimg with PackageCompiler.jl

There is also technically an interpreter if you want to go that way [1], so in principle it might be possible to do the same trick javascript does, but someone would have to implement that.

[1] https://github.com/JuliaDebug/JuliaInterpreter.jl

Ah, I see. That explains a lot. Could also be the reason why it beats gfortran in the test. A JIT knows more about the data the program sees/runtime behavior than a normal compiler. Still surprising. Or is gfortran just so much behind state of the art fortran compilers?

Not this test. Julia’s JIT compiler does not know about runtime values, it’s not a tracing JIT like the ones you may be used to. Instead, it’s ‘just ahead of time’, compiling a specialized executable for each new function signature the first time that signature gets hit. I.e. if I do

    f(x) = 2x + 1

The first time I call

    f(1)

it’ll compile a function specialization for f( ::Int), and the first time you call f on that integer it’ll be slow and all the subsequent calls will be equally fast.

Next if you do

    f(1.0 + 2im)

it’ll compile a new specialization for f(::Complex{Float64}) which will be slow the first time while it compiles and then fast on all the subsequent runs.

The genius of Julia’s design is that the JIT compiler is designed around the semantics of multiple dispatch, and the multiple dispatch semantics are designed around having a JIT

That’s strictly accurate but kinda misleading. It doesn’t have multiple optimization passes it does just for hot code paths, it fully compiled per function. They call it Just Ahead Of Time compiled.

Think about compiling a large C++ program and its dependencies every time you run the generated executable. This is sort of like what is Julia is doing.

Wow, I may have to give Julia a try. I've been holding back because of that classic language issue "it doesn't have the libraries I like" but these performance numbers are impressive and hard to ignore.

The need for libraries is a lot different, too. In Python when you need a library what you specifically might need is a high performance Python wrapper around a C library. Just writing your own short thing isn’t an option because pure Python is so slow. In Julia that’s not the case.

In addition to everything else, Julia makes it easy to call functions from Python, C, and Fortran libraries.

The Fortran doesn't look quite equivalent, and we don't know how it was run, though that's irrelevant to the analysis except perhaps for OpenMP. As usual, it's a single number from a benchmark without understanding it, e.g. with -fopt-info. The general conclusion about Fortran invalid.

It happens to be an example of a pathological vectorization case for GCC with an ancient issue open, sigh [1]. Much as it pains me to say so, apart from NAG's, I think, all the other compilers I know will vectorize such loops: ifort, xlf, flang, PGI, nvfortran (but the last three may be essentially the same here).

Also, last time I tried it, -floop-nest-optimize (marked "experimental") was quite broken, again with an issue open either in the gcc or redhat tracker. -O3 does unroll-and-jam, at least.

1. https://gcc.gnu.org/bugzilla/show_bug.cgi?id=40770

You can speed up the numpy code by 2x by replacing the line:

    M_old = np.exp(1j*k*(np.tile(a,(N,1))**2 + np.tile(a.reshape(N,1),(1,N))**2))

with

    a_squared = a ** 2
    M = np.exp(1j * k * (a_squared[:, None] + a_squared[None, :]))

This gets you to the same performance as single-threaded julia.

Whenever people run benchmarks on numerical code, especially involving square root, I immediately wonder if the generated output is either ´correct’ or just approximate.

There is no ‘correct’ in floating point, just approximate. But I would guess that such a talented group as Julia’s would pay extreme attention to numerical stability/best implementation of algorithms, or fall back to a lib everyone uses already.

i did notice in the benchmark game they were using fast math with a float and an iteration to get the inverse sqrt. but the c++ and rust examples were doing the same

Is the python code missing a square root?

https://github.com/mdmaas/julia-numpy-fortran-test/blob/main...

For fun I ran in Matlab with a 2.9 GHz i7-7820HQ and get about 1.83s for N=10,000 single threaded.

    A = exp((k*1i)*sqrt(a.^2 + (a.^2)'))

I find the argument that vectorized code is more difficult to read weird. Having to write out the loops becomes tedious very quickly and does often significantly decrease readibility. Take for example the element wise multiplication of two 2d square arrays and compare that to the multiplication of an array and a transpose. In the loop the only difference is that the indices are flipped, easy to overlook in a quick read.

My biggest annoyance with Julia is however that they decided to follow matlab and do matrix operations by default and even worse use the '.' as the element wise operation modifier. From my own experience and from teaching many students, one of the primary source of issues is getting that wrong. I don't know how many times when debugging some weird matlab issue it was a matter of find the missing '.' they could have chosen any other symbol but instead the opted for the easiest one to overlook.

IMO, doing matrix operations by default is the mathematically correct thing to do. E.g. having exp(M::Matrix) give anything other than the matrix exponential is insanity. Same with A::Matrix * B::Matrix giving anything other than matrix multiplication.

As far as the ‘.’ for elementwise operations, I don’t think it’s that hard to miss, but I guess my eyes are trained for it by now.

There’s also an @. macro to make an entire expression element wise though if needed, e.g.

    @. h(f(x) + g(y))

is equivalent to

    h.(f.(x) .+ g.(y))

I think this broadcasting machinery is one of the closet parts of julia because you get to choose at the call site if you want the function to apply element-wise or to see the whole array, and it works on any container, with user customizable behaviour for user defined containers, and has all sorts of goodies like automatic loop fusion.

I personally think a vector library that allowed some form of Einstein notation [1] would be a killer feature for a scientific computing environment.

[1] https://en.wikipedia.org/wiki/Einstein_notation

Ask and ye shall receive (more than ye bargained for):

https://github.com/mcabbott/Tullio.jl

https://github.com/Jutho/TensorOperations.jl

https://github.com/under-Peter/OMEinsum.jl

Tullio.jl is the most powerful and interesting to me of the above libraries, but they all have cool strengths.

Having zero experience with this I could see myself getting used to the notation. The dot before the addition sign is strange looking to unfamiliar eyes, is it calling + on the result of f.(x)?

Yeah that’s eighth, it’s element wise addition between the results of f.(x) and g.(y) with loop fusion.

the syntax for broadcasting a regular function is

    f.(x)

but the syntax for broadcasting an infix function like +, is

    x .+ y

so the dot goes at the start. I think there were some convincing reasons it ended up this way, but I forget what those reasons are. It’s just muscle memory for me at this point. You can also treat any infix function as a regular prefix one by wrapping it in parens, e.g.

    (+).(x, y)

I would argue that an array is something else than a matrix. It's a design decision that matlab (and julia) made, regard every multi-element variable as a matrix. I disagree with that decision, because for most of the code I'm exposed to I don't want matrix operations and it causes issues, e.g. when I move from a scalar implementation to an array implementation of a function.

Ah sure, yes we did not need to conflate vectors and matrices with arrays and I sometimes wonder if we made the wrong choice there.

I guess I'd say that then someone would get to own the syntax

   [a, b, c]

it'd either be a vector of a non-mathematical array and whichever choice was made, someone would be unhappy. I think not too much was lost by having them be the same thing, but the ship has sailed anyways.

I think this depends upon the language intrinsic support for vectors/vector operations. I would (generally) prefer to specify operations using a simpler, more concise syntax, and have the compiler map that into the hardware. The problem with this is, most compilers don't do such a good job with that.

I like the Matrix/vector operations by default for the reason that it becomes quite a bit simpler to understand the written code. Expressing things the way Numpy forces you to, means that you have to reason about the interface of your code to Numpy, as well as how Numpy will operate.

For matrix mult, I agree, writing out loops is tedious. This is why things like M = A * B, where M, A, B are all matrices, is far, far better than writing out the loop representation. I agree as well, messing up the index ordering is annoying (and I still do this 30+ years onward with Fortran, C, etc.)

Julia makes this stuff easy to trivial. It makes the interface to using this capability easy to trivial. It makes for good overall performance (I rarely run codes only once).

1. in Julia you can't write "vecotrize-styled" code by using `@.` without manually adding the dots everywhere and the result will be as fast in numpy "classic use case"

2. Have not-slow for-loops is invaluable precisely because sometimes your workload is not (either too hard, or impossible/bad in terms of RAM usage) suitable for vectorzie-styled code.

1. I know and pretty much every line of code I have would need to be preceded with a @. I should also add that julia did not follow the matlab insanity of regarding 1D vectors as arrays, but instead throwing errors helps catch errors before they cause weird issues.

2. I agree that this is one of the nice things about Julia.

In general, what you want to do instead of having `@.` in front of every line is write scalar functions and `@.` their application. For example, instead of writing ``` function f1(A,B,C) tmp = A.+B return tmp1 .* C end

f1(::Vector,::Vector,::Vector) ``` You can write

``` function f2(a,b,c) tmp = a + b return tmp1 *c end

f2.(::Vector,::Vector,::Vector) ```

Can someone explain how Julia MT got a super-linear improvement in performance by using 4 cores + HT vs 1 core?

Fortran’s MT performance is more along the lines of what I would expect.

Genuine question: Is Julia used for anything other than mathematics yet? Is there any sign of the Julia community expanding beyond science, maths, machine learning etc?

Yes, people are using it for databases, web servers, video games, and a bunch of other things. These domains are less mature than julia's traditional wheelhouse of numerical / technical computing, but I think eventually it can be very competent at almost anything people put the time into building.

Julia is an incredibly flexible language that's all about enabling powerful code transformations, so it can be adapted to very different niches.

which library do you use for databases? when i checked a while ago, the situation was not very good and the best option was the JDBC interface

is this still the case

If you are interested in datasets that fit in memory, Dataframes.jl is really good now. For interfacing with SQL, there still isn't anything much better than jdbc.

Does the Julia runtime have good properties for something like a Web server?

Slow startup makes it questionable for CLI tools, although that is supposedly improving a lot.

But what about high-concurrency "server-like" workloads? The language itself I think would be great for such an application, but I have no idea if the runtime itself would be good.

Currently Julia will probably lose to C++ here due to some GC quirks and less well optimized strings. That said, both of these issues are being worked on, and I would expect fairly major progress by 1.8. (some of this work is already done. Specifically, as of 2 days ago memory alignment of strings has changed to save 8 bytes for most strings). Future improvements will likely include some form of small string optimization to prevent them from being heap allocated, and making julia's GC more efficient for highly threaded code.

Yes, but a lot of it isn't as highly optimized as some of the more mathy stuff. Julia has a ton of really good readers of various file formats (CSV, JSON, Arrow, etc). Many of these are as fast as anything out there, but some of them have some performance issues still. https://www.genieframework.com/ is a web framework entirely in Julia, (but probably not the fastest).

If you have specific things you are interested in, feel free to follow up.

Does it have to? IMO it's good to have a language focused on that domain as the tradeoffs tend to be different there, hence why Fortran is still actively used.

I use it for general application development. Some of the libraries still need quite a bit of work, but it's getting there. For practical web-service development, you need to hook in Revise.jl to get updates to your request handler without having to restart the server.

There’s a handful of things in a bunch of different domains, although mathy things are the most well developed. The language is definitely general purpose and suitable in a ton of different domains. It’s just going to take time for the ecosystem to build up.

I have to ask to ask the obvious question: why?

What problems would they be trying to solve?

Gradually-typed, syntactic macros, tidy syntax, good Unicode support in the standard library. All positive traits for app and/or server development in my opinion.

Maybe a nitpick but why would someone put time on Y axis? I remember it's DRYMIX for graphs - dependent on Y axis (N), independent on X axis (time)

Because what's measured is how how the runtime (on the Y axis) increases as the matrix size N (on the X axis) increases (i.e., you set N and measure t, rather than setting t and measuring N)

The article is surprising as it illustrates the reason for switching to Julia with only a very simplistic single example. Surely there is more to that than this example especially considering that the Python implementation could be improved for readability and speed.

Another thing to add, the python example code is a written in quite a weird (and in my opinion difficult to read way).

n = len(a) M = np.exp(1jk(np.tile(a,(n,1))*2 + np.tile(a.reshape(n,1),(1,n))*2))

is equivalent to

M = np.exp(1jk(a[:,None]*2 + a[None,:]*2))

which I find very easy to grasp.

A big improvement, but i'm not sure if indexing with None is very intuitive.

np.newaxis instead of None is the idiomatic way

Interesting to see.. for so many years julia was so much in a blur.

Is gfortran considered a good fortran compiler? Back in my day, I think it still generated C code that then was compiled with gcc. I might be wrong, this was in the 90s...

It is. Maybe you are confusing it with f2c, which was a different project that translated Fortran 77 to C.

Probably. I think there was g77, but gfortran and g95 apparently have very limited relations to this.

The i used in the index Mij, is not the same i as in the multiplication ik?

How does one make that out in the formula?

It's too hard to read on mobile device.

What makes Julia so fast?

Being designed to be fast. Lots of other high level languages were designed with the assumption that performance critical code would be written in low-level languages so they didn't bother designing the language to make various optimzations possible.

Julia was designed for people who don't want to use a low level language for speed. It took a lot of great ideas from languages like Common Lisp, Fortran, etc. as well as a few good ideas of its own and packaged them together in a rather novel way.

The language is dynamic, but it's dynamism is limited in specific ways that makes it very easy for the compiler to optimize code.

Functions compile just-ahead-of-time to native code using LLVM. Functions specialize based on the input types used and use type inference, meaning that in idiomatic Julia code, all types are known at compile time. Finally, Julia uses struct layout similar to static languages like C or Rust.

The trifecta of complete type information, compilation to native code, and efficient data layout, is what makes the difference.

Notably, Julia is not any faster than most static languages, nearly all of whom also use the same trio.

Numba exists, which is a vectorizing, optimizing JIT compiler for a subset of Python that uses NumPy. You simply install it and add a decorator.

https://numba.pydata.org

The problem is that Numba's subset is really limited and often means you have to spend a significant time investment on making your code compatible, often while decreasing readability or doing unnecessary type conversions.

Numba is very cool, but it's inherently limited.

Which interpreted language has the fastest startup time.

Question is which Fortran vs Which Version of Julia

Reply on news.ycombinator.com