Array comprehension is inconsistent, given the column-major property of julia

Internals & Design
,
1

Related doc: Performance Tips · The Julia Language

For a matrix

julia> x = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

The correct order of element is

julia> for e = x
           println(e)
       end
1
3
2
4

To exactly reproduce the above sequence in a double-index fashion, we have to write

julia> for c=1:2, r=1:2
           println(x[r, c])
       end
1
3
2
4

Now I’m expecting that

[x[r, c] for c=1:2, r=1:2] == x

returns true, which is deemed consistent. But the existing behavior is not.

(I found this yesterday when I was thinking about the JuMP container Specifying objective coefficient when creating a JuMP variable)

1 Like
2

In both cases the fastest changing iterator is closest to the loop body, the comprehension is just written with the body first which is why it has the opposite order of the for loop. I think this makes sense somewhat although I also have to think for a moment first sometimes to pick the right order for each.

4 Likes
3

[x[r, c] for r=1:2, c=1:2]: r is the rows, c is the columns. It will give a multi-dimension array, where the index does not rely on each other/another/any other when the array is generating. Totally it is not the syntactic sugar of for loop.

1 Like
4

@WalterMadelim the other thing to keep in mind is: the behaviour is what it is. We won’t be changing it because it would be breaking, both in JuMP and in Julia. Just learn to live with previous decisions :smile:

7 Likes
5

You’re suggesting that

[x[r, c] for r=1:2, c=1:2] == x

where == reads “recovers”. Okay I think this make sense at the beginning for everyone, as it is highly intuitive.

The main concern is that

  • It doesn’t co-exist well with julia’s column-major Arrays, which I pointed out in my #1 post.
6

I think the crux of this topic is that:

for a julia Matrix A_{rc}, whose element is called by A[r, c], the c index is considered at a higher rank.

In this expression for c=1:2, r=1:2, the c index is considered at a higher rank.

The word “higher” is also perceived as being “outer”, “independent”, e.g. the month index relative to the day index (31 days in the first month, 28 days in the second month).

Therefore I think in the julia world, letting

be true should make more sense.

Some more instances I come up with 
julia> Calendar2026 = [[1, 2, "...", 31], [1, 2, "...", 28]];

julia> for m = 1:2 # 🟦
           for d = Calendar2026[m]
               println("month $m, day $d")
           end
       end
month 1, day 1
month 1, day 2
month 1, day ...
month 1, day 31
month 2, day 1
month 2, day 2
month 2, day ...
month 2, day 28

julia> for m = 1:2, d = Calendar2026[m] # 🟧 readily derived from the above
           println("month $m, day $d")
       end
month 1, day 1
month 1, day 2
month 1, day ...
month 1, day 31
month 2, day 1
month 2, day 2
month 2, day ...
month 2, day 28

julia> @assert [[d for d = Calendar2026[m]] for m = 1:2] == Calendar2026 # corresponds to 🟦, this is okay

julia> Calendar2026Matrix = [[1, 2, "...", 31] [1, 2, "...", 28]] # if organized into Matrix, we should end up with this
4×2 Matrix{Any}:
  1        1
  2        2
   "..."    "..."
 31       28

julia> for m = 1:2, d_ind = 1:4
            println("month $m, day $(Calendar2026Matrix[d_ind, m])")
        end
month 1, day 1
month 1, day 2
month 1, day ...
month 1, day 31
month 2, day 1
month 2, day 2
month 2, day ...
month 2, day 28

julia> [Calendar2026Matrix[d_ind, m] for m = 1:2, d_ind = 1:4] # fails to recover the primal
2×4 Matrix{Any}:
 1  2  "..."  31
 1  2  "..."  28

Opened an issue

7

Identical issue has been discussed a number of times here. In particular, I recommend Why column major?, in which you can find a few useful explanations. Perhaps my take Why column major? - #36 by zdenek_hurak could be one of them.

1 Like
8

Also relevant is this post. :wink:

5 Likes
9

I agree with @jules that this convention makes sense when you switch from prefix for to postfix for. For me it’s natural to consider

[x[r, c] for r=1:2, c=1:2]

to be equivalent to something like

[(x[r, c] for r=1:2) for c=1:2]

(though this is already valid syntax creating a Vector of generators). I.e. here r changes the most rapidly.

On the other hand

for r=1:2, c=1:2
    something_with(x[r, c])
end

is the same as

for r=1:2
    for c=1:2
        something_with(x[r, c])
    end
end

where then c changes the most rapidly.

You could also consider [x[r, c] for r=1:2 for c=1:2] (without the parentheses resulting in generators), which for some reason behaves differently

julia> vec([x[r, c] for r=1:2, c=1:2])  # as you would expect if you agree with my reasoning above
4-element Vector{Int64}:
 1
 3
 2
 4

julia> vcat([[x[r, c] for r=1:2] for c=1:2]...)  # again, expected
4-element Vector{Int64}:
 1
 3
 2
 4

julia> [x[r, c] for r=1:2 for c=1:2]  # unexpected
4-element Vector{Int64}:
 1
 2
 3
 4

I’m not sure what’s going on here.

julia> ex = :([x[r, c] for r=1:2 for c=1:2])
:([x[r, c] for r = 1:2 for c = 1:2])

julia> ex.args[1].args[1]
:(((x[r, c] for c = 1:2) for r = 1:2))

I don’t immediately see the connection with the topic at hand, so you should probably expand on this a bit.

4 Likes
10

This is exactly where the current comprehension design is not consistent.

To follow the current comprehension grammar, one needs to remember the “folklore” convention

, which is a bit indirect by itself. And it clashes somewhat with

julia> [(i, j) for i=1:3 for j=1:i]
6-element Vector{Tuple{Int64, Int64}}:
 (1, 1)
 (2, 1)
 (2, 2)
 (3, 1)
 (3, 2)
 (3, 3)

It’s also inconsistent in JuMP when we use different type of containers:

julia> import JuMP

julia> model = JuMP.Model();

julia> JuMP.@variable(model, x[m = 1:2, d = 1:2]) # month, day
2×2 Matrix{JuMP.VariableRef}:
 x[1,1]  x[1,2]
 x[2,1]  x[2,2]

julia> for e = x # day-major (i.e. slowest changing)
           println(e)
       end
x[1,1]
x[2,1]
x[1,2]
x[2,2]

julia> JuMP.@variable(model, y[m = 1:2, d = ifelse(m>1, (1, 28), (1, 31))])
JuMP.Containers.SparseAxisArray{JuMP.VariableRef, 2, Tuple{Int64, Int64}} with 4 entries:
  [1, 1 ]  =  y[1,1]
  [1, 31]  =  y[1,31]
  [2, 1 ]  =  y[2,1]
  [2, 28]  =  y[2,28]

julia> for e = y # month-major
           println(e)
       end
y[1,1]
y[1,31]
y[2,1]
y[2,28]

I’m not sure if there would be any inconsistency if the (convention-free) idea in my #1 post were adopted. I’m not sure if that would be a more coherent design. Admittedly, I’ll have to get used to the existing behavior at the current stage.

11

By the way, this topic is essentially a duplicate of

In particular, check out okatsn’s summary at the end.

4 Likes
12

The discussion there is of high quality. But I think my depiction of this issue in my #1 post here conveys a sharper contradiction.

First of all, that post is quite lengthy and make use of if-else cases analysis. So the current design is intrinsically involved.

I think it doesn’t make sense to explain or understand it. What we can do in practice is just to use them warily—just like from

x, y = a, b
...

to

let x, y = a, b
    ...
end

, where the let stealthily modifies the coder’s intention, which is unfortunate.

13

Simpler example that intends 1-based indices of 2 rows and 3 columns. These are “consistent”:

julia> for c in 1:3
         for r in 1:2 # iterate rows directly in column-major
           println((r, c))
         end
       end
(1, 1)
(2, 1)
(1, 2)
(2, 2)
(1, 3)
(2, 3)

julia> for c in 1:3, r in 1:2
           println((r, c))
       end
(1, 1)
(2, 1)
(1, 2)
(2, 2)
(1, 3)
(2, 3)

julia> [(r, c) for c in 1:3 for r in 1:2]
6-element Vector{Tuple{Int64, Int64}}:
 (1, 1)
 (2, 1)
 (1, 2)
 (2, 2)
 (1, 3)
 (2, 3)

but this is not. It looks row-major, but we’re actually iterating c directly instead:

julia> [(r, c) for c in 1:3, r in 1:2]
3×2 Matrix{Tuple{Int64, Int64}}:
 (1, 1)  (2, 1)
 (1, 2)  (2, 2)
 (1, 3)  (2, 3)

The reason is actually pretty intuitive: the array’s dimensions 3x2 match the order of the iterables’ lengths.

If we want to emulate the first 3 examples, we have to reverse the order so the implied dimensions 2x3 also match:

julia> [(r, c) for r in 1:2, c in 1:3]
2×3 Matrix{Tuple{Int64, Int64}}:
 (1, 1)  (1, 2)  (1, 3)
 (2, 1)  (2, 2)  (2, 3)

Let’s say the order of the dimensions is left to right as written, and that rows must be to the left of columns for matrices. Column-major array dimensions expand from left to right (rows, columns, beyond) in the dimensions’ order, so the most efficient for-loop must iterate right to left (innermost loop over rows) in the reverse order. That’s the root cause of the for c in 1:3, r in 1:2 “inconsistency” between for-loops and comprehensions, and any reasonable switch would just move that “inconsistency” somewhere else.

On the other hand, row-major array dimensions expand from right to left (columns, rows, beyond) in the reverse order, so the for-loop would iterate left to right (innermost loop over columns) in the dimensions’ order. The reverse order isn’t written as often e.g. aligning dimensions for broadcasting, so this is one of the aesthetic arguments in favor of row-major.

Something to consider for a row-major array package or a different base language, I doubt even Julia v2 would make this big a change to Array instead of disambiguating for r in 1:2, c in 1:3 between for-loops and comprehensions. There’s just too much tradition in column-major matrices e.g. BLAS buffers, though it’s worth mentioning that transpose flags handle row-major matrices (I believe at zero cost but I’m not sure).

4 Likes