Why column major?

Why Julia arrays are stored in column-major order?

Column-major vs row-major are not entirely symmetric regarding code style, because when we access the i,j'th element of an array, we write A[i,j] with the row index first, which can break the symmetry.

For example, I find it surprising that in:

[(println(i, " ", j); (i,j)) for i=1:5, j=1:5]

i is traversed faster than j, even though it looks like j is an inner loop.
This is completely different from a true nested for loop:

for i=1:5, j=1:5 
  println(i, " ", j)
end 

Here j is indeed varying faster. If Arrays were row-major, these two code snippets would vary j faster, which looks more consistent to me given that in both cases the iteration is specified by for i=1:5, j=1:5.

Another example is:

[(println(i, " ", j); (i,j)) for i=1:5 for j=1:5]

Here also j varies faster, which makes [... for i=1:5, j=1:5] different from for i=1:5 for j=1:5. The difference is not just that the later is linearized, but also the order in which elements are traversed. Again, there would be no difference in order if arrays were row major.

I’m just wondering about the reasoning behind the column-major choice, and if there is a more “natural” way to think about this that I am missing.

Probably this has been discussed before on Github, though I can’t find it now.

Row major is a really strange choice if you think about what happens when you add dimensions… How would you iterate A[i,j,k,...]? In row major languages, the fastest way is typically to iterate A (from inner loop to outer) j, i, k,... which is just horrible. In julia, you always iterate the dimensions in order A[i,j,k,...].

The order in the loop comprehension might be a bit confusing though but it still makes sense. The first one produces a vector and has the last loop as the tightest. The second one produces a matrix and has the first dimension as the tightest.

1 Like

Row-major layout does not require any odd iteration order. In numpy, A[i,j,k,l] is most efficiently iterated in the the order i,j,k,l from outer to inner loop (and each A[i,j], which is the same as A[i,j,:,:], is a row-major matrix). Ditto for the typical way to emulate multi-dimensional arrays in C.
The thing I miss most about numpy after moving to Julia is its indexing behavior where A[i]==A[i,...], A[i,j] == A[i,j,...] etc, and A[i,j]==A[i][j]. Seems much more natural and useful to me than Julia’s punning of linear indexing with multi-dimensional indexing, which I only found to get in the way of finding bugs.

As I understand it, it was mainly historical precedent, driven by the fact that most linear-algebra libraries (LAPACK etc.) are column-major (because Fortran is column-major, and hence Matlab is column-major), though it would have been possible to work around this with judicious use of transposition.

It’s been discussed before, e.g. https://groups.google.com/d/msg/julia-users/GkvDunqE0ww/2gIDmA5cAAAJ

9 Likes

It also aligns quite nicely with how we (and the typical mathematical notations) treat 1-dimensional vectors: as columns.

7 Likes

Row vs column major is one of those topics where good arguments can be made for either choice, but at the end of the day it is best to just pick something and proceed. Given that BLAS/LAPACK mostly use colum-major, I think it is a reasonable choice to follow.

For people who really like row-major, it should be feasible to implement it using a custom type without a performance penalty — with a lot of work, of course.

I did this when implementing a Common Lisp (which is row-major) wrapper for BLAS/LAPACK. It works in about 90% of the cases (not all primitives handle transposes), and was always a reliable source of bugs.

3 Likes