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.

7 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:

10 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.

6 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 after the last 3 examples, but we’re actually iterating c across rows 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 indices (1:3, 1:2). If we want to emulate the first 3 examples where r represented rows, 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 adhere to conventional matrix dimensions where rows are to the left of columns. Column-major array dimensions thus expand from left to right (rows → columns → beyond…) in the dimensions’ order, and the most efficient for-loop iterates indices from right to left (innermost loop over rows) in the reverse order. That’s the root cause of the “inconsistency” between for-loops and comprehensions. Breaking changes would just move that “inconsistency” somewhere else; I’d be fine with reversing the order of the comma-delimited for-loop and maybe the comma-less comprehension to innermost-first, but someone else would probably consider that inconsistent with the irreversible nested for-loop.

On the other hand, row-major array dimensions expand from right to left (beyond… ← rows ← columns) in the reverse order, so the for-loop would iterate indices from 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. I doubt even Julia v2 would make this big a change to Array; there’s just too much tradition in column-major matrices, even preceding Julia e.g. BLAS buffers, though it’s worth mentioning that transpose flags handle row-major matrices there (I believe at zero cost but I’m not sure).

6 Likes
14

I like Benny’s recap, which is quite clear.

I want to emphasize one thing:

  • This topic is not a duplicate of the past discussions.

My idea is much simpler, which cares only the front-end behavior. To understand my claim, you don’t need any knowledge about “column-major” or the efficiency of the backend implementation.

My claim in the #1 post was only suggesting that the final outcome should obey

image
image505×360 17.2 KB

, as a user I don’t need to look into how julia actually process that comprehension under the hood.

The pain with the current design is essentially caused by this redundant condition.

This mindset was essentially needless, and is incompatible with Julia’s choice on Array mechanism.

What I’m thinking currently:

Some constructs in julia are just not that perfect and 100% consistent. What practitioners can do, from the safest point of view, is just to bypass them. e.g. if we think carefully, it appears that not all of them are used frequently, e.g. let, local etc.

I think the basic for loop (not those in comprehensions) works good. And we have functions, that’s basically enough. Just like to communicate thoughts I only need a vocabulary about 3000.

15

You’re making a bad assumption here that the comprehension syntax is only about column-major Arrays. We can easily replace the brackets [] with () and make a generator expression, which is not even an AbstractArray. There are no eager or lazy elements to index, but it must be materialized into the same abstract array regardless of how the concrete array is implemented. Matching the comprehension indices to that abstract array’s indices (axes and size also work on generators) is the only sensible neutral option. It’d be difficult to convince people to reverse the comprehension indices just because it matches the innermost-row-iterating for-loops over a specific Array type; it wouldn’t match the innermost-column-iterating for-loops over row-major array types or for-loops over more exotic types, so why bother reversing it from its only meaningful property if it’d be generally “inconsistent” anyway?

1 Like
16

I didn’t get your idea.

In this topic I proposed only one concept: self-consistency.
I didn’t find where generators can violate this self-consistency, since generators cannot be indexed.

Maybe you can come up with an instance to illustrate?

I think the current Array design—the elements in the first axis are organized contiguously in memory—is reasonable. The unfortunate point is that this style stands at the opposite of how we instinctively write normal loops

for i = 1:I, j = 1:J, k = 1:K
    A[i, j, k]
end

Although this is unfortunate, I think it is reasonable, not “inconsistent”.

The real place where I find inconsistent is that Julia already goes with the “unfortunate” choice, yet still resort to the “instinctive” style of Array comprehension syntax

[(i, j, k) for i = 1:I, j = 1:J, k = 1:K]

, which is perceived as an object that belongs to, e.g. Python (numpy).

Now my point is

  • Can we make our syntax consistent, despite it being “unfortunate” already?

I think it is possible.

If you disagree, I’d like to see a code example.

17

You’re contradicting yourself here. The “front-end behavior” of multidimensional comprehension is already independent of the “backend implementation” because its for syntax is based on the order of the dimensions/axes, but you’re insisting on a consistency with for-loops that do depend on Array’s layout. A nested for-loop for a column-major array is not written the same way as a for-loop for a row-major array, and a for-loop for an array with contiguous elements along a middle axis won’t be written like either of those loops. It’s impossible to make comprehension’s for syntax consistent with efficient nested for-loops over arrays in general (eachindex is a simple sequence that is not always as efficient as possible for a given concrete type), you can only be consistent with 1 particular layout’s. My earlier point was that the loops over row-major arrays just happen to align with the dimensions’ order, not that the loops are a better basis.

You mentioned NumPy, so I’ll note something funny here. 1-dimensional list comprehension is the only idiomatic option, so multidimensional NumPy arrays are converted from nested lists, the innermost list being a row. As a consequence of comprehension syntax being led by the element expression, the iterations are written in the reverse order of efficient for-loops over row-major arrays:

>>> for r in range(1, 3): 
...   for c in range(1, 4): # iterate columns directly in row-major
...     print(r, c)
... 
1 1
1 2
1 3
2 1
2 2
2 3

>>> np.array([[10*r+c for c in range(1,4)] for r in range(1,3)])
array([[11, 12, 13],
       [21, 22, 23]])

So even a language whose comprehension’s for syntax is based directly on nested for-loops doesn’t have your idea of “consistency”. Earlier when I said I could accept reversing orders in nested for-loops, I was following this lead.

The absence of indexing is irrelevant. A generator represents a lazily iterable version of a comprehension by definition, they share the same syntax, and comprehensions lower to collecting a corresponding generator. Since a generator can be used to fill an array of any layout, the syntax has a fundamental reason not to be based on any particular array layout.

1 Like
18

Considerations are already abundant.

To me, the simplest way out is to write customized code using the ordinary for loop instead.

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

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

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

julia> function my_procedure(x)
           x_new_2 = Matrix{Int}(undef,2,2)
           for c=1:2, r=1:2 if x[r, c] > 0
               x_new_2[r, c] = x[r, c]
           end end
           x_new_2
       end;

julia> x_new_2 = my_procedure(x)
2×2 Matrix{Int64}:
 1  2
 3  4

I have no additional comments. Different new views from others may still be appended subsequently…

19

FWIW, I think the relation to how some array implementation stores its data in memory is not so relevant for the discussion about the syntax.

To me the surprise here is, that the following things are different:

julia> [println((i,j)) for i in 1:2 for j in 3:4]
(1, 3)
(1, 4)
(2, 3)
(2, 4)

julia> [println((i,j)) for i in 1:2, j in 3:4]
(1, 3)
(2, 3)
(1, 4)
(2, 4)

whereas the same things written as for-loops are the same:

julia> for i in 1:2
           for j in 3:4
               println((i,j))
           end
       end
(1, 3)
(1, 4)
(2, 3)
(2, 4)

julia> for i in 1:2, j in 3:4
           println((i,j))
       end
(1, 3)
(1, 4)
(2, 3)
(2, 4)

I agree that this is a footgun. However, I also think that each of the constructs by itself makes sense:

  1. Nested for-loops: trivially correct
  2. For loop with comma: short-hand for the nested for-loop
  3. comprehension with 2 for: simple rewrite of the nested for-loops
  4. comprehension with comma: constructs an array/iterator with dimensions matching the iterators (oops different order than the other 3)

For me the fundamental reason was not stated clearly so far: The comma in the comprehension means something fundamentally different, then the comma in the for loops! The comprehension with comma is a completely different iteration construct conceptually. Using a comma for it is an abuse of notation, if you wish. To make this more explicit consider dependent loop bounds:

julia> for i in 1:2, j in i:3
           println((i,j))
       end
(1, 1)
(1, 2)
(1, 3)
(2, 2)
(2, 3)

julia> [println((i,j)) for i in 1:2, j in i:3]
ERROR: UndefVarError: `i` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
Stacktrace:
 [1] top-level scope
   @ REPL[13]:1

julia> [println((i,j)) for j in i:3, i in 1:2]
ERROR: UndefVarError: `i` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
Stacktrace:
 [1] top-level scope
   @ REPL[14]:1

This shows directly that the comprehension with , is a fundamentally different construct. So much so, that I would have liked a different syntax for it just to avoid exactly the confusion that we discussed here in this thread. To emphasize again: The true origin of the confusion is not there is an inconsistency in Julia’s iteration/comprehension construct but rather that there are actually 2 conceptually different constructs that use very similar notation such that they can be confused for another.

Personally, I think it would have been a better choice to use ; instead , in the comprehension because ; is also used to separate dimension in array literals. This would resolve the issue discussed here because it much easier to accept that

for i in 1:2, j in 3:4
    #...
end

and

[#= =# for i in 1:2; j in 3:4]

are different constructs following different rules.

11 Likes
20

Thanks. Your account is very well organized, I’m fully convinced.

The newly suggested [#= =# for i=1:2; j=3:4] reads much clearer so I can perceive that i=1:2 is the inner layer and i alters more frequent than j does (and this corresponds to julia’s Array memory layout, which is the correct thing).

One more concern: how will the if cond(i, j) be inserted so that I can still perceive the hierarchy correctly from the syntax?

e.g. how to explain

sum(x[i, j] if i>j for i=1:3; j=2:4)

or

sum(x[i, j] for i=1:3; j=2:4 if i>j)

? (I think the latter style reads awkward—my eye needs to jump left first and then jump right)

Or, say, the front-end syntax should let people know its direct programming counterpart: is it true that the counterpart of the former sum(...) style corresponds to

for j=2:4, i=1:3
    if i>j
       ....

? (as I currently understand).

Moreover, can we write

sum(x[i, j] for i=1:3 if i>2 for j=2:4)

, which corresponds to

for i=1:3
    if i> 2
        for j=2:4
            ...

?

Even going further to write something more crazy like

maximum(x[i, j, k] for i=1:4; j=1:3 if i>j for k=5:6)

We need to stipulate the rule carefully.