In Defense of Matlab Code

> I don't think Julia really solves any problems that aren't already solved by Python.

I don't really need proper furniture, the cardboard boxes and books setup I had previously "solved" the same problems, but I feel less worried about random parts of it suddenly buckling, and it is much more ergonomic in practice too.

This is a huge understatement. At the hedge fund I work at, I learned Julia by porting a heavily optimized Python pipeline. Hundreds of hours had gone into the Python version – it was essentially entirely glue code over C.

In about two weeks of learning Julia, I ported the pipeline and got it 14x faster. This was worth multiple senior FTE salaries. With the same amount of effort, my coworkers – who are much better engineers than I am – had not managed to get any significant part of the pipeline onto Numba.

> And if something is truly performance critical, it should be written or rewritten in C++ anyway.

Part of our interview process is a take-home where we ask candidates to build the fastest version of a pipeline they possibly can. People usually use C++ or Julia. All of the fastest answers are in Julia.

But isn't the whole point of this article that Matlab is more readable than Python (i.e. solves the readability problem)? The Matlab and Julia code for the provided example are equivalent[1]: which means Julia has more readable math than Python.

[1]: Technically, the article's code will not work in Julia because Julia gives semantic meaning to commas in brackets, while Matlab does not. It is perfectly valid to use spaces as separators in Matlab, meaning that the following Julia code is also valid Matlab which is equivalent to the Matlab code block provided in the article.

    X = [ 1 2 3 ];
    Y = [ 1 2 3;
          4 5 6;
          7 8 9 ];
    Z = Y * X';
    W = [ Z Z ];

It's wild what people get used to. Rustaceans adapt to excruciating compile times and borrowchecker nonsense, and apparently Pythonistas think it's a great argument in favor of Python that all performance sensitive Python libraries must be rewritten in another language.

In fairness, we Julians have to adapt to a script having a 10 second JIT latency before even starting...

Our goal is to make a runtime that lets people stay at the math layer as much as possible, and run the math as fast as possible.

Doesn't just (Y @ X)[None] work? None adding an extra dimension works in practice but I don't know if you're "supposed" to do that

I know what all of these do, but I just can’t shake the feeling that I’m constantly fighting with an actual python. Very aptly named.

I also think it’s more to do with the libraries than with the language, which I actually like the syntax of. But numpy tries to be this highly unopinionated tool that can do everything, and ends up annoying to use anywhere. Matplotlib is even worse, possibly the worst API I’ve ever used.

When I've had to write similar code in Python, it's a massive pain to "prove" that my conversation code is correct. Often I've resorted to using MATLAB's trusted functions to generate "truth" data and then feeding that to Python to verify it gets the same results.

Obviously this is more work than just using the premade stuff that comes with the toolbox.

Any MATLAB alternative faces the same trust issue. Until it reaches enough mindshare that people assume that it's too popular to have incorrect math (which might not be a good assumption but it is one that people make about MATLAB) then it doesn't actually mimic the main benefit of MATLAB which is that I don't need to check its work.

We like Octave a lot, but the reason we started fresh is architectural: RunMat is a new runtime written in Rust with a design centered on aggressive fusion and CPU/GPU execution. That’s not a small feature you bolt onto an older interpreter; it changes the core execution model, dataflow, and how you represent/optimize array programs.

Could you add a JIT to Octave? Maybe in theory, but in practice you’d still be fighting the existing stack and end up with a very long, risky rewrite inside a mature codebase. Starting clean let us move fast (first release in August, Fusion landed last month, ~250 built-ins already) and build toward things that depend on the new engine.

This isn’t a knock on Octave, it’s just a different goal: Octave prioritizes broad compatibility and maturity; we’re prioritizing a modern, high-performance runtime for math workloads.

Julia has OffsetArrays.jl implementing arbitrary-base indexing: https://juliaarrays.github.io/OffsetArrays.jl/stable/

The experience with this has been quite mixed, creating a new surface for bugs to appear. Used well, it can be very convenient for the reasons you state.

  julia> A = collect(1:5)
  5-element Vector{Int64}:
   1
   2
   3
   4
   5

  julia> B = OffsetArray(A, -1)
  5-element OffsetArray(::Vector{Int64}, 0:4) with eltype Int64 with indices 0:4:
   1
   2
   3
   4
   5

  julia> A[1]
  1

  julia> B[0]
  1

But okay -- as I mentioned, you're entitled to your views!

Cheers for the comment!

(There is unhelpful subtext here that I can't possibly have done serious math, but putting that aside...) On the contrary, most actual linear algebra is easier when you have real 1D arrays. Compare an inner product form in Matlab:

   x' * A * y

   x @ A @ y

> math texts ... are all extremely clear about column vs row vectors and notation too, and all make it clear whether column vs. row vector is the default notation, and use superscript transpose accordingly.

That's because they use the blunt tool of matrix multiplication for composing their tensors. If they had an equivalent of the @ operator then there would be no need, as in the above formula. (It does mean that, conversely, numpy needs a special notation for the outer product, whereas if you only ever use matrix multiplication and column vectors then you can do x * y', but I don't think that's a big deal.)

> This is also a constant issue working with scikit-learn, and if you regularly read through the source there, you see why.

I don't often use scikit-learn but I tried to look for 1D/2D agreement issues in the source as you suggested. I found a couple, and maybe they weren't representative, but they were for functions that could operate on a single 1D vector or could be passed as a 2D numpy array but, philosophically, with a meaning more like "list of vectors to operate on in parallel" rather than an actual matrix. So if you only care about 1d arrays then you can just pass it in (there's a np.newaxis in the implementation, but you as the user don't need to care). If you do want to take advantage of passing multiple vectors then, yes, you would need to care about whether those are treated column-wise or row-wise but that's no different from having to check the same thing in Matlab.

Notably, this fuss is precisely not because you're doing "real linear algebra" - again, those formulae are (usually) easiest with real 1D arrays. It when you want to do software-ish things, like vectorise operations as part of a library function, that you might start to worry about axes.

> unless you ingrain certain habits to always call e.g. .ravel or .flatten or [:, :, None] arcana

You shouldn't have to call .ravel or .flatten if you want a 1D array - you should already have one! Unless you needlessly went to the extra effort of turning it into a 2D row/column vector. (Or unless you want to flatten an actual multidimensional array to 1D, which does happen; but that's the same as doing A(:) in Matlab.)

Writing foo[:, None] vs foo[None, :] is no different from deciding whether to make a column or row vector (respectively) in MATLAB. I will admit it's a bit harder to remember - I can never remember which index is which (but I also couldn't remember without checking back when I used Matlab either). But that's because it's more general in that it works for higher dimensions too. Plus, as I've said, you should rarely need it in practice.