# Multidimensional array comprehension
**URL:** https://discourse.julialang.org/t/multidimensional-array-comprehension/134004
**Category:** General Usage
**Tags:** matrix, comprehension
**Created:** [November 20, 2025, 9:25am UTC](https://discourse.julialang.org/t/multidimensional-array-comprehension/134004 "2025-11-20T09:25:18Z")
**Posts on this page:** 14
**Page:** 1
### Author:  [@fps](https://discourse.julialang.org/u/fps)
#### Post date: [November 20, 2025, 9:25am UTC](https://discourse.julialang.org/t/multidimensional-array-comprehension/134004/1 "2025-11-20T09:25:18Z")
I wonder it is possible to get the following matrix with a single array comprehension that’s also efficient to evaluate:
```julia-auto
(1, 1, 1) (1, 2, 1) (1, 1, 2) (1, 2, 2) (1, 1, 3) (1, 2, 3)....
(2, 1, 1) (2, 2, 1) (2, 1, 2) (2, 2, 2) (2, 1, 3) (2, 2, 3)
(3, 1, 1) (3, 2, 1) (3, 1, 2) (3, 2, 2) (3, 1, 3) (3, 2, 3)
.
.
.
```
Note the alternating second entry and the third entry increment every second column. I sometimes need this kind of comprehension to quickly build up a fourier basis. Something like this works:
```julia-auto
julia> [(a, mod(b-1, 2)+1, div(b+1, 2)) for a = 1:3, b = 1:6]
3×6 Matrix{Tuple{Int64, Int64, Int64}}:
(1, 1, 1) (1, 2, 1) (1, 1, 2) (1, 2, 2) (1, 1, 3) (1, 2, 3)
(2, 1, 1) (2, 2, 1) (2, 1, 2) (2, 2, 2) (2, 1, 3) (2, 2, 3)
(3, 1, 1) (3, 2, 1) (3, 1, 2) (3, 2, 2) (3, 1, 3) (3, 2, 3)
```
Or more explicitly for this use-case:
```julia-auto
function fourier_basis(Ts, order)
[sin(2f0 * pi .* (div(o+1, 2) .* t) .+ (mod(o, 2) * pi/2f0)) for t in Ts, o in 1:(2*order)]
end
julia> fourier_basis(0:0.01:1, 2)
101×4 Matrix{Float64}:
1.0 0.0 1.0 0.0
0.998027 0.0627905 0.992115 0.125333
0.992115 0.125333 0.968583 0.24869
0.982287 0.187381 0.929776 0.368125
0.968583 0.24869 0.876307 0.481754
0.951057 0.309017 0.809017 0.587785
0.929776 0.368125 0.728969 0.684547
0.904827 0.425779 0.637424 0.770513
0.876307 0.481754 0.535827 0.844328
0.844328 0.535827 0.425779 0.904827
0.809017 0.587785 0.309017 0.951057
0.770513 0.637424 0.187381 0.982287
0.728969 0.684547 0.0627905 0.998027
0.684547 0.728969 -0.0627906 0.998027
0.637424 0.770513 -0.187381 0.982287
⋮
0.637424 -0.770513 -0.187381 -0.982287
0.684547 -0.728969 -0.0627902 -0.998027
0.728969 -0.684547 0.0627908 -0.998027
0.770513 -0.637424 0.187382 -0.982287
0.809017 -0.587785 0.309017 -0.951056
0.844328 -0.535827 0.42578 -0.904827
0.876307 -0.481754 0.535827 -0.844328
0.904827 -0.425779 0.637424 -0.770513
0.929777 -0.368124 0.728969 -0.684547
0.951057 -0.309017 0.809017 -0.587785
0.968583 -0.24869 0.876307 -0.481753
0.982287 -0.187381 0.929777 -0.368124
0.992115 -0.125333 0.968583 -0.24869
0.998027 -0.0627903 0.992115 -0.125333
1.0 1.74846e-7 1.0 3.49691e-7
```
but it’s a little unwieldy. Of course I could go the complex number route, but then I would need to destructure the complex entries in another step which might be unnecessary 🙂
I tried different combinations of the two iteration primitives in array comprehensions `for a in A, b in B` and `for a in A for b in B` and I couldn’t find one that did the above in a single step.
---
### Author:  [@fps](https://discourse.julialang.org/u/fps)
#### Post date: [November 20, 2025, 9:31am UTC](https://discourse.julialang.org/t/multidimensional-array-comprehension/134004/2 "2025-11-20T09:31:55Z")
OK, sometimes you just need to talk about it to then find a solution. In this case I can just do:
```julia-auto
julia> reshape([(a, b, c) for a = 1:3, b = 1:2, c = 1:3], 3, 6)
3×6 Matrix{Tuple{Int64, Int64, Int64}}:
(1, 1, 1) (1, 2, 1) (1, 1, 2) (1, 2, 2) (1, 1, 3) (1, 2, 3)
(2, 1, 1) (2, 2, 1) (2, 1, 2) (2, 2, 2) (2, 1, 3) (2, 2, 3)
(3, 1, 1) (3, 2, 1) (3, 1, 2) (3, 2, 2) (3, 1, 3) (3, 2, 3)
```
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [November 20, 2025, 12:31pm UTC](https://discourse.julialang.org/t/multidimensional-array-comprehension/134004/3 "2025-11-20T12:31:01Z")
I would just define something like
```julia
basis_index(a, b) = ((d, r) = divrem(b + 1, 2); (a, r + 1, d))
```
and then use comprehension; but are you sure that you actually need to instantiate this matrix?
---
### Author:  [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)
#### Post date: [November 20, 2025, 12:58pm UTC](https://discourse.julialang.org/t/multidimensional-array-comprehension/134004/5 "2025-11-20T12:58:58Z")
Minor side comment:
> [@fps](#):
>
> `mod(b-1, 2)+1, div(b+1, 2)`
This is the same as `fldmod1(b, 2)` although in reverse order (assuming b isn’t negative).
---
### Author:  [@fps](https://discourse.julialang.org/u/fps)
#### Post date: [November 20, 2025, 1:35pm UTC](https://discourse.julialang.org/t/multidimensional-array-comprehension/134004/6 "2025-11-20T13:35:06Z")
Yes, it is used as an input to a neural network. So sadly it has to “exist” 😉
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [November 20, 2025, 1:42pm UTC](https://discourse.julialang.org/t/multidimensional-array-comprehension/134004/7 "2025-11-20T13:42:24Z")
Is the calculation in Julia though? If yes, you can code up an `<:AbstractArray` that is never instantiated, just calculated on demand when `Base.getindex` is called.
---
### Author:  [@jishnub](https://discourse.julialang.org/u/jishnub)
#### Post date: [November 20, 2025, 3:37pm UTC](https://discourse.julialang.org/t/multidimensional-array-comprehension/134004/8 "2025-11-20T15:37:02Z")
Alternately,
```julia-auto
julia> Tuple.(reshape(CartesianIndices((3,2,3)), 3, 6))
3×6 Matrix{Tuple{Int64, Int64, Int64}}:
(1, 1, 1) (1, 2, 1) (1, 1, 2) (1, 2, 2) (1, 1, 3) (1, 2, 3)
(2, 1, 1) (2, 2, 1) (2, 1, 2) (2, 2, 2) (2, 1, 3) (2, 2, 3)
(3, 1, 1) (3, 2, 1) (3, 1, 2) (3, 2, 2) (3, 1, 3) (3, 2, 3)
```
---
### Author:  [@mikmoore](https://discourse.julialang.org/u/mikmoore)
#### Post date: [November 20, 2025, 3:45pm UTC](https://discourse.julialang.org/t/multidimensional-array-comprehension/134004/9 "2025-11-20T15:45:40Z")
There are other cute games you can play with this particular use case:
```julia-repl
julia> reinterpret(Float64, [cispi(2*Ts*o) for o in 1:2, Ts in 0:0.01:1])'
101×4 adjoint(reinterpret(Float64, ::Matrix{ComplexF64})) with eltype Float64:
1.0 0.0 1.0 0.0
0.998027 0.0627905 0.992115 0.125333
0.992115 0.125333 0.968583 0.24869
0.982287 0.187381 0.929776 0.368125
0.968583 0.24869 0.876307 0.481754
0.951057 0.309017 0.809017 0.587785
⋮
0.968583 -0.24869 0.876307 -0.481754
0.982287 -0.187381 0.929776 -0.368125
0.992115 -0.125333 0.968583 -0.24869
0.998027 -0.0627905 0.992115 -0.125333
1.0 0.0 1.0 0.0
```
(The axis reordering and transposing was necessary here so that the real/imag values would go into alternating columns after the `reinterpret`)
---
### Author:  [@fps](https://discourse.julialang.org/u/fps)
#### Post date: [November 20, 2025, 4:16pm UTC](https://discourse.julialang.org/t/multidimensional-array-comprehension/134004/10 "2025-11-20T16:16:30Z")
That is a good question! For bigger problems it might be worth instantiating them “on the fly” for each batch. Right now they fit nicely into memory (even on the GPU).
---
### Author:  [@fps](https://discourse.julialang.org/u/fps)
#### Post date: [November 20, 2025, 4:18pm UTC](https://discourse.julialang.org/t/multidimensional-array-comprehension/134004/11 "2025-11-20T16:18:26Z")
Good to know - thanks!
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [November 20, 2025, 6:10pm UTC](https://discourse.julialang.org/t/multidimensional-array-comprehension/134004/12 "2025-11-20T18:10:02Z")
A comprehension solution (but not efficient):
```julia
[(a, d...) for a=1:3, d in Tuple((b, c) for b=1:2, c=1:3)]
```
This other comprehension is efficient but not as nice looking:
```julia
[(a, b, c) for a=1:3, (b,c) in vec(collect(Iterators.product(1:2, 1:3)))]
```
---
### Author:  [@fps](https://discourse.julialang.org/u/fps)
#### Post date: [November 21, 2025, 6:45am UTC](https://discourse.julialang.org/t/multidimensional-array-comprehension/134004/13 "2025-11-21T06:45:01Z")
Great answers! Thanks to all of you. Hard to choose any particular one as “the” solution 🙂
---
### Author:  [@fps](https://discourse.julialang.org/u/fps)
#### Post date: [November 21, 2025, 8:00am UTC](https://discourse.julialang.org/t/multidimensional-array-comprehension/134004/14 "2025-11-21T08:00:24Z")
Sounds interesting. Do you happen to have a pointer or two to relevant documentation?
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [November 21, 2025, 8:31am UTC](https://discourse.julialang.org/t/multidimensional-array-comprehension/134004/15 "2025-11-21T08:31:59Z")
The [array interface](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array) is super-simple to implement, eg
```julia
struct BasisMatrix <: AbstractMatrix{NTuple{3,Int}}
nrow::Int
ncol::Int
end
basis_index(a, b) = ((d, r) = divrem(b + 1, 2); (a, r + 1, d))
Base.size(bm::BasisMatrix) = (bm.nrow, bm.ncol)
function Base.getindex(bm::BasisMatrix, r, c)
(; nrow, ncol) = bm
(1 ≤ r ≤ nrow && 1 ≤ c ≤ ncol) || throw(BoundsError(bm, (r, c)))
basis_index(r, c)
end
```
gives you
```julia
julia> BasisMatrix(3,4)
3×4 BasisMatrix:
(1, 1, 1) (1, 2, 1) (1, 1, 2) (1, 2, 2)
(2, 1, 1) (2, 2, 1) (2, 1, 2) (2, 2, 2)
(3, 1, 1) (3, 2, 1) (3, 1, 2) (3, 2, 2)
```