A Matlab programmer's take on Julia

M. Bazdresch

[I am running Julia Version 0.0.0-prerelease, Commit 6f25d4d84a (2012-03-01 07:50:31)]

Julia is a new language for numerical computing. It is fast (comparable to C), its syntax is easy to pick up if you already know Matlab, supports parallelism and distributed computing, has a neat and powerful typing system, can call C and Fortran code, and includes a pretty web interface. It also has excellent online documentation. Crucially, and contrary to SciPy, it indexes from 1 instead of 0.

For a Matlab programmer, it sounds like a dream come true, right? (I mean Matlab the language, not the commercial software; I use Octave exclusively). The Matlab language is slow, it is crufty, and has many idiosyncracies. A few examples:

  • A line of code not ending in a semicolon will print its result to stdout.

  • One function per .m file, both with the same name. The file may contain other functions but they are private to the function whose name corresponds to the file name.

  • No namespaces.

  • You can do things such as

    sin=3;
    

    which shadows the sin function; doing clear sin will clear the variable and restore the function.

  • Many functions do very different things depending on how they are called and even how many output arguments are specified.

And the list goes on and on.

I strongly disagree, however, with the opinion, common among some circles, that Matlab is to be dismissed just because it is crufty or "not well designed". It is actually a very productive language that is very well suited to numerical computing and algorithm exploration. Cruftiness and slowness are the price we pay for its convenience and flexibility.

However, Julia's allure is strong, so strong that I have spent some hours exploring it, reading the manual, and perusing their GitHub repository. My conclusion is that I will not switch to it in the short term, but I will keep an eye on it and periodically re-evaluate it.

I believe some of the concerns I have about Julia will be solved eventually. Some relate to missing infrastructure. Others are small things that would make me less productive if I switched to Julia today, compared to Octave, regardless of its strengths in speed, parallelism, etc. Some things that come to mind are:

  • No plotting, except in the web REPL, and only very simple (if pretty) 2D plots supported.

  • No debugging facilities. No need for fancy IDE-style setting of breakpoints; Octave-style keyboard and db* functions would be good enough.

  • No test suite.

  • No code autoloading. When exploring an algorithm or ironing out bugs, I usually have an Octave session in a terminal and my code in a gvim window. I change something in the code, alt-tab to the terminal and call my code. Octave will figure out by itself that the file containing the code has changed from the last time it was called, so it will reload it automatically. I can be sure Octave is always running the latest version of my code.

    In Julia, no autoloading is done. This means that I have to issue load("mycode.jl") before running my modified code. This is a minor but very real nuisance.

  • No package management or global configuration. Julia will read a custom.jl file in the directory it's run from, if present, but there seems to be no way to do this globally. I can invoke Octave from anywhere in my system, and I will have access to all my packages without further action. Readline, the prompt, etcetera will all be set up the way I want.

  • Sandboxing and security in the web REPL. How will I be able to read and write data files, while preserving filesystem privacy and security? I don't have the answer, and nothing about this is mentioned in the manual. Julia has the ability to run any function on the system's libc, which in principle means it will be able to wipe my home directory or upload my files to the cloud.

  • The way vectors are handled is, frankly, weird. Let's build a column vector:

    julia> x = [1,2,3]
    [1, 2, 3]
    

    Note that, even though it's printed "horizontally", x is actually a column vector and it has only one dimension.

    julia> size(x)
    (3,)
    julia> ndims(x)
    1
    

    Now let's build a row vector:

    julia> y = [1 2 3]
    1x3 Int64 Array:
     1  2  3
    

    Note that, for some reason, Julia now prints the type of y along with its value. y is a row vector and, contrary to x, it is two-dimensional.

    julia> size(y)
    (1,3)
    julia> ndims(y)
    2
    

    Under multiplication, x does indeed behave like a column vector:

    julia> y*x
    [14]
    

    Curiously, if we transpose y, we can build a two-dimensional column vector:

    julia> yt = permute(y,[2,1])
    3x1 Int64 Array:
     1
     2
     3
    julia> ndims(yt)
    2
    

    x and y may be multiplied, as shown above, but yt and x cannot be added, even though both are column vectors, since they have different number of dimensions:

    julia> x+yt
    argument dimensions must match
    

    The purpose of this way of handling vectors escape me, and I fail to see how it will make numerical programming easier or more powerful. Maybe a few years of Matlab programming will do this to you, but I find Matlab's conventions much more natural and unobtrusive.

These and other similar issues are probably to be expected in a language as young as Julia. Some of these are the kind of convenience features you want to build once the rest of the infrastructure is solid. I am not too worried about them; Julia will probably solve these issues in time, or I will just decide to put up with them the way I put up with Matlab today. Ohter issues, however, I am more concerned about since they seem to be fundamental to the way the language is being designed, and are not so easy to work around:

  • About that web REPL, it seems that the authors' intention is that it should be the main interface to Julia, relegating the terminal to second class status (see this discussion). Even if the web REPL is "very slick", I am not sure I will ever want it as my primary development environment. For instance, I frequently have four or five figures open, and I want to be able to see them at all times. In the web repl, plots will scroll back as more commands are input. This is a serious inconvenience. A web interface is awesome for use in tablets or phones, but when I'm sitting on my desk I want a terminal and multiple windows.

  • Julia's sqrt(-1) returns NaN and according to the manual, a math function with real arguments will return reals, while the same function with complex arguments will return complex numbers. I can't quite put my finger on it, but at first sight this behavior bothers me a lot, and I don't think it makes mathematical sense.

    I can't think of a single instance where I have wished sqrt(-1) returned NaN just because -1 is real. If I wanted that behavior, I would wrap sqrt in a realsqrt function that checked the argument's type. What happens most of the time is that I don't know in advance whether I will be taking square roots of negative numbers, and I would hate getting NaNs polluting my results if that happened.

    Furthermore, this behavior is not consistent. If you calculate the fast fourier transform of a real vector, you get a complex result:

    julia> fft([1.0, 2.0, 1.0])
    [6.0+0.0im, -1.5+0.866025im, -1.5-0.866025im]
    

    Speaking of the fft command, for some reason it won't take integer arguments:

    julia> fft([1, 2, 1])
    no method fft(Array{Int64,1},(),Int64)
    

    A type system such as Julia's may give a language a lot of power and speed; it may also straightjacket the programmer. I don't want to have to make sure all types are correct before calling fft. I just want to transform the vector and go on with my calculations.

  • Optional arguments in functions, versus multiple dispatch. It is a fact of life that, in numerical computing, you will want to create functions with lots of arguments, some of which are optional. I can appreciate the elegance and cleanliness of multiple dispatch versus the "crufty" Matlab approach (using introspection to figure out just what arguments the function was provided with). For an example, see my spectrum front-end function for the fast fourier transform. It has seven input arguments, and only one is mandatory. This provides for a function that is very convenient (I can just type spectrum(x) and instantly get an idea of x's frequency content, or I can fine-tune every aspect of the plot if I choose to).

    Under Julia's multiple dispatch philosophy, I would have to write seven versions of spectrum. Or, I could ignore the benefits of multple dispatch and instead write spectrum as a vararg function. Doing so, however, results in a less readable function; compare Matlab's

    function [ssf afxs pfxs] = spectrum(x, Ts, N, win, dp, clip, pctrl)
    

    to Julia's vararg

    function spectrum(x, ...)
    

    In Matlab's case, the function argument names self-document their purpose to anyone familiar with the knowledge domain: win is for windowing, Ts is the sampling interval, etcetera. Furthermore, in Julia the varargs tuple x will have to be processed without the convenience of Matlab's nargin function.

  • Single return values from functions. Julia's functions always return a single value; multiple values are returned as elements in a tuple. This is quite inconvenient. In the case of my example function spectrum, if I want just one of its outputs, then I call it with just one assignment:

    s = spectrum(x);
    

    This means I don't have to extract the result I'm interested in from a tuple. In Julia, I could do

    (s,z,zz) = spectrum(x);
    

but I end up with variables z and zz that clutter the workspace and take up memory. Furthermore, using Matlab's introspective function nargout, I can avoid calculating values that the user didn't ask for, speeding up computation.

As I said, I am not sure these problems will ever be solved, since they seem to be inherent to the way Julia has been designed. If they are worked around, one risks losing Julia's elegance and simplicity. For instance, see this discussion on named function arguments. This is a great feature (Python has it, I wish Octave had it too). Implementing this might involve sidestepping multiple dispatch and an "out of band" structure that is passed to the functions. The danger of this is, of course, that the language might end up being just as crufty as Matlab's.

As a conclusion, I trust the authors will stay true to their vision and turn Julia into a language that preserves most of Matlab's convenience and at the same time meets all of their promises. It'd be a dream come true for those of us who spend a good deal of time doing numerical computation. This post may sound too negative, but I complain because I care. For the time being, alas, it's back to Octave for me.