Reductions, reshaping and the dims keyword

sum(M; dims = 2) of a 2×3 matrix does not give the two row sums as a vector: it gives a 2×1 matrix, which Julia prints as [6; 15;;]. The kept dimension is deliberate, because it lines the sums up with M for the next broadcast, and it is the first thing to check when code written for a vector rejects the result.

This lesson reduces along a chosen dimension with sum, maximum, count, mean and mapreduce, runs argmax, extrema, cumsum and diff along rows, joins arrays with hcat, vcat, cat and stack, tells reshape from permutedims and from the lazy adjoint A', and walks a station × day × variable array with eachslice. Every run below was recorded on Julia 1.13.1; one function needs Julia 1.12, and its run on the 1.10 LTS is shown too.

dims says which dimension to reduce, and keeps it as length 1

A reduction such as sum folds many values into one. Called on a whole matrix it adds everything, and with the keyword dimsdims: The keyword of sum, maximum, count, mean and other reductions that names the dimensions to reduce: sum(M; dims = 1) adds down each column of M and keeps that dimension with length 1. it adds along one dimension only:

# A reduction over the whole matrix, down its columns (dims = 1) and along its rows (dims = 2).
using Statistics                     # mean: a standard library that comes with Julia
M = [1 2 3; 4 5 6]
println("sum(M)           = ", sum(M))
println("sum(M; dims = 1) = ", sum(M; dims = 1), ", size ", size(sum(M; dims = 1)))
println("sum(M; dims = 2) = ", sum(M; dims = 2), ", size ", size(sum(M; dims = 2)))
println("vec(sum(M; dims = 2))                = ", vec(sum(M; dims = 2)))
println("dropdims(sum(M; dims = 2); dims = 2) = ", dropdims(sum(M; dims = 2); dims = 2))
println("maximum(M; dims = 1) = ", maximum(M; dims = 1), ", minimum(M; dims = 2) = ", minimum(M; dims = 2))
println("mean(M; dims = 2) = ", mean(M; dims = 2), ", sum.(eachrow(M)) = ", sum.(eachrow(M)))
println("each value as a share of its row total: ", round.(M ./ sum(M; dims = 2); digits = 3))

from session import files, run

# Sums, extremes and means of a 2×3 matrix, whole and along each dimension.
files("dims.jl")
run("""
julia dims.jl
""")

$ julia dims.jl
sum(M)           = 21
sum(M; dims = 1) = [5 7 9], size (1, 3)
sum(M; dims = 2) = [6; 15;;], size (2, 1)
vec(sum(M; dims = 2))                = [6, 15]
dropdims(sum(M; dims = 2); dims = 2) = [6, 15]
maximum(M; dims = 1) = [4 5 6], minimum(M; dims = 2) = [1; 4;;]
mean(M; dims = 2) = [2.0; 5.0;;], sum.(eachrow(M)) = [6, 15]
each value as a share of its row total: [0.167 0.333 0.5; 0.267 0.333 0.4]

Step through both reductions:

sum(M) gave 21. sum(M; dims = 1) added down each column and gave the 1×3 matrix [5 7 9], and sum(M; dims = 2) added along each row and gave a 2×1 matrix. The documentation of sum says it sums the elements of an array over the given dimensions, and its examples keep the reduced dimension with length 1. vec and dropdims(…; dims = 2) remove it when a plain vector is wanted: both gave [6, 15]. maximum, minimum and count take dims the same way, and so does mean from Statistics, the standard library for basic statistics: mean(M; dims = 2) gave the row means 2.0 and 5.0.

The kept dimension is what makes the next step a single broadcast. M ./ sum(M; dims = 2) divided each row by its own total, 0.167, 0.333 and 0.5 in the first row, because a 2×1 matrix stretches across the 3 columns of M.

argmax, extrema, cumsum and diff work along a dimension too

The other reductions and the running operations follow the same keyword. A vector first, then the rows of M:

# Other reductions, and running operations, along a dimension.
v = [3, 9, 2, 5]
println("argmax(v) = ", argmax(v), ", findmax(v) = ", findmax(v), ", extrema(v) = ", extrema(v))
println("cumsum(v) = ", cumsum(v), ", diff(v) = ", diff(v))
M = [1 2 3; 4 5 6]
println("extrema(M; dims = 2) = ", extrema(M; dims = 2))
println("findmax(M; dims = 2) = ", findmax(M; dims = 2))
println("cumsum(M; dims = 2)  = ", cumsum(M; dims = 2))
println("diff(M; dims = 2)    = ", diff(M; dims = 2))
println("count(>(2), M; dims = 1)        = ", count(>(2), M; dims = 1))
println("mapreduce(abs2, +, M; dims = 1) = ", mapreduce(abs2, +, M; dims = 1))
try
    cumsum(M)                       # a matrix has no single order to accumulate in
catch err
    println(sprint(showerror, err))
end

from session import files, run

# argmax, extrema, findmax, cumsum, diff, count and mapreduce, on a vector and along the rows of a matrix.
files("along.jl")
run("""
julia along.jl
""")

$ julia along.jl
argmax(v) = 2, findmax(v) = (9, 2), extrema(v) = (2, 9)
cumsum(v) = [3, 12, 14, 19], diff(v) = [6, -7, 3]
extrema(M; dims = 2) = [(1, 3); (4, 6);;]
findmax(M; dims = 2) = ([3; 6;;], CartesianIndex{2}[CartesianIndex(1, 3); CartesianIndex(2, 3);;])
cumsum(M; dims = 2)  = [1 3 6; 4 9 15]
diff(M; dims = 2)    = [1 1; 1 1]
count(>(2), M; dims = 1)        = [1 1 2]
mapreduce(abs2, +, M; dims = 1) = [17 29 45]
UndefKeywordError: keyword argument `dims` not assigned

argmax([3, 9, 2, 5]) gave 2, the position of the largest value; findmax gave the value and its position together, (9, 2), and extrema the smallest and the largest as one tuple, (2, 9). With dims = 2 each works row by row: extrema gave one tuple per row, and findmax returned the row maxima with their positions as CartesianIndex values, the positions of the previous lesson but one.

cumsum and diff are running operations rather than reductions: the running totals keep every position, and the differences between neighbours are one shorter than their input. On a matrix both need dims, because a matrix has no single order to accumulate in: cumsum(M) stopped with UndefKeywordError, as the array form in the documentation takes dims without a default. count(>(2), M; dims = 1) counted per column, and mapreduce(abs2, +, M; dims = 1) squared and added per column, [17 29 45], without an array of squares in between.

hcat, vcat, cat and stack join arrays along a dimension

Joining is the opposite of slicing, and it too is a question of which dimension:

# Joining vectors and matrices: hcat, vcat, cat with dims, stack and reduce.
a = [1, 2, 3]
b = [4, 5, 6]
println("hcat(a, b) = ", hcat(a, b), ", size ", size(hcat(a, b)))
println("vcat(a, b) = ", vcat(a, b), ", size ", size(vcat(a, b)))
println("[a b] == hcat(a, b): ", [a b] == hcat(a, b), ", [a; b] == vcat(a, b): ", [a; b] == vcat(a, b))
println("size(cat(ones(2, 2), zeros(2, 2); dims = 3)) = ", size(cat(ones(2, 2), zeros(2, 2); dims = 3)))
days = [[12.0, 0.0, 7.5], [0.0, 2.5, 11.0], [30.5, 18.0, 0.0], [4.0, 9.5, 1.5]]   # one vector per day
println("stack(days) = ", stack(days), ", size ", size(stack(days)))
println("size(stack(days; dims = 1)) = ", size(stack(days; dims = 1)))
println("reduce(hcat, days) == stack(days): ", reduce(hcat, days) == stack(days))
println("size(reduce(vcat, days)) = ", size(reduce(vcat, days)))

from session import files, run

# Two vectors joined side by side and end to end, two matrices stacked in a third dimension, and four days stacked.
files("concat_shapes.jl")
run("""
julia concat_shapes.jl
""")

$ julia concat_shapes.jl
hcat(a, b) = [1 4; 2 5; 3 6], size (3, 2)
vcat(a, b) = [1, 2, 3, 4, 5, 6], size (6,)
[a b] == hcat(a, b): true, [a; b] == vcat(a, b): true
size(cat(ones(2, 2), zeros(2, 2); dims = 3)) = (2, 2, 2)
stack(days) = [12.0 0.0 30.5 4.0; 0.0 2.5 18.0 9.5; 7.5 11.0 0.0 1.5], size (3, 4)
size(stack(days; dims = 1)) = (4, 3)
reduce(hcat, days) == stack(days): true
size(reduce(vcat, days)) = (12,)

hcat(a, b) put the two vectors side by side as the columns of a 3×2 matrix, and vcat(a, b) end to end as one vector of 6. [a b] and [a; b] are the same two calls, as the manual’s table of concatenation syntax shows, which is why the literals of the first lesson behave as they do. cat(…; dims = 3) joins along any dimension, so two 2×2 matrices became one 2×2×2 array.

To build a matrix from many vectors, stack puts each vector in a column of its own; its documentation says it needs Julia 1.9 or later. The four days of rainfall became the 3×4 matrix of the lesson on views, and stack(days; dims = 1) put them in rows instead, a 4×3 matrix. reduce(hcat, days) gave the same matrix as stack, and the documentation of vcat notes that reduce(vcat, A) calls an efficient method for a large vector of arrays.

reshape keeps the storage order, permutedims moves the elements

reshape and permutedims both change a shape, and they do different things to the values:

# reshape keeps the storage order; permutedims moves the elements to swap the dimensions.
M = [1 2 3; 4 5 6]
println("reshape(M, 3, 2) = ", reshape(M, 3, 2))
println("permutedims(M)   = ", permutedims(M))
v = [1, 2, 3]
println("reshape(v, 1, 3) = ", reshape(v, 1, 3), ", permutedims(v) = ", permutedims(v))
r = insertdims(v; dims = 1)          # needs Julia 1.12 or later
println("insertdims(v; dims = 1) = ", r, ", size ", size(r))
println("dropdims(r; dims = 1) = ", dropdims(r; dims = 1))

from session import files, run, use

# reshape against permutedims, then insertdims, which the 1.10 LTS does not have.
files("reshape_permute.jl")
run("""
julia --version
julia reshape_permute.jl
""")
use("1.10.12")
run("""
julia --version
!julia reshape_permute.jl
""")

$ julia --version
julia version 1.13.1
$ julia reshape_permute.jl
reshape(M, 3, 2) = [1 5; 4 3; 2 6]
permutedims(M)   = [1 4; 2 5; 3 6]
reshape(v, 1, 3) = [1 2 3], permutedims(v) = [1 2 3]
insertdims(v; dims = 1) = [1 2 3], size (1, 3)
dropdims(r; dims = 1) = [1, 2, 3]
$ julia --version
julia version 1.10.12
$ julia reshape_permute.jl
reshape(M, 3, 2) = [1 5; 4 3; 2 6]
permutedims(M)   = [1 4; 2 5; 3 6]
reshape(v, 1, 3) = [1 2 3], permutedims(v) = [1 2 3]
ERROR: LoadError: UndefVarError: `insertdims` not defined
Stacktrace:
 [1] top-level scope
   @ ~/demo/reshape_permute.jl:7
in expression starting at ~/demo/reshape_permute.jl:7
[exit status 1]

reshape(M, 3, 2) reread the same six slots as 3 rows and 2 columns and gave [1 5; 4 3; 2 6]: the storage order of the first lesson, not a transpose. permutedims(M) moved every element to the swapped position and gave [1 4; 2 5; 3 6], the rows of M as columns. For a vector the two agree: reshape(v, 1, 3) and permutedims(v) both gave the row [1 2 3].

insertdims(v; dims = 1) adds a dimension of length 1 at the position named, here turning the vector into a 1×3 matrix, and dropdims removes it again. The 1.12 release notes list insertdims as new, and on the 1.10.12 LTS the same file stopped at that line with UndefVarError: `insertdims` not defined; reshape(v, 1, 3) gives the shape on every version.

M’ is a lazy adjoint, permutedims a new matrix

The postfix ' gives the adjointadjoint: The conjugate transpose of a matrix or vector, written A'. Julia returns it as a lazy wrapper that reads A, swapping the indices and conjugating complex values as they are read. of a matrix, its conjugate transpose, and Julia does not build a new matrix for it. A complex matrix shows what that means:

# A' is a lazy wrapper that conjugates; transpose wraps without conjugating; permutedims copies.
A = [1+2im 3-1im; 0+1im 4+0im]
println("typeof(A') = ", typeof(A'))
println("A'             = ", A')
println("transpose(A)   = ", transpose(A))
println("permutedims(A) = ", permutedims(A))
At = A'
P = permutedims(A)
A[1, 2] = 100 + 0im
println("after A[1, 2] = 100: A' reads ", At, ", the copy reads ", P)
x = [1.0, 2.0, 3.0]
println("size(x') = ", size(x'), ", typeof(x') = ", typeof(x'))
println("typeof(collect(A')) = ", typeof(collect(A')))

from session import files, run

# A complex matrix through ', transpose and permutedims, then one write into the original.
files("adjoint_vs_permutedims.jl")
run("""
julia adjoint_vs_permutedims.jl
""")

$ julia adjoint_vs_permutedims.jl
typeof(A') = LinearAlgebra.Adjoint{Complex{Int64}, Matrix{Complex{Int64}}}
A'             = Complex{Int64}[1 - 2im 0 - 1im; 3 + 1im 4 + 0im]
transpose(A)   = Complex{Int64}[1 + 2im 0 + 1im; 3 - 1im 4 + 0im]
permutedims(A) = Complex{Int64}[1 + 2im 0 + 1im; 3 - 1im 4 + 0im]
after A[1, 2] = 100: A' reads Complex{Int64}[1 - 2im 0 - 1im; 100 + 0im 4 + 0im], the copy reads Complex{Int64}[1 + 2im 0 + 1im; 3 - 1im 4 + 0im]
size(x') = (1, 3), typeof(x') = LinearAlgebra.Adjoint{Float64, Vector{Float64}}
typeof(collect(A')) = Matrix{Complex{Int64}}

Step through it and watch which result sees the write into A:

The type of A' is LinearAlgebra.Adjoint{Complex{Int64}, Matrix{Complex{Int64}}}, which its documentation describes as a lazy wrapper: it stores no values of its own and reads A as its adjoint. It swaps the indices and conjugates every value as it is read, so 3 - 1im at A[1, 2] reads as 3 + 1im at position (2, 1). transpose(A) swaps without conjugating, and permutedims(A) built an ordinary Matrix with the same swap and no conjugation. The write A[1, 2] = 100 made the difference visible: the adjoint, still reading A, showed 100 + 0im, while the copy kept 3 - 1im.

The documentation says the adjoint type is intended for linear algebra and points to permutedims for general data manipulation, and copy materialises the wrapper: collect(A') returned an ordinary Matrix{Complex{Int64}}. For real numbers ' and transpose give the same values, and for a vector x' is a 1×3 Adjoint row, which is how x' * x became a dot product in the previous lesson.

eachslice walks a 3-D array in views, mapslices builds a new one

A third dimension adds a third index. Here it is the variable, rainfall or maximum temperature, beside station and day, in invented readings of three stations:

# Rainfall in mm and maximum temperature in °C of 3 stations over 6 days (invented), as one 3×6×2 array.
stations = ["Pune", "Chennai", "Delhi"]
rain = [12.0  0.0 30.5  4.0 18.0  0.0;
         0.0  2.5  8.0 41.0  9.5 22.0;
         7.5 11.0  0.0  1.5  0.0  3.0]
tmax = [29.5 31.0 27.0 30.5 28.0 31.5;
        33.0 34.5 32.0 29.5 31.0 30.0;
        35.5 34.0 36.5 37.0 36.0 35.0]
A = cat(rain, tmax; dims = 3)                       # station × day × variable
println("size(A) = ", size(A))
println("total rain per station: ", vec(sum(A[:, :, 1]; dims = 2)))
println("mean maximum temperature per station: ", round.(vec(sum(A[:, :, 2]; dims = 2)) ./ size(A, 2); digits = 2))
wettest, at = findmax(A[:, :, 1]; dims = 2)
for s in 1:3
    println(rpad(stations[s], 8), "wettest on day ", at[s][2], ", ", wettest[s], " mm")
end
anomaly = A[:, :, 1] .- sum(A[:, :, 1]; dims = 2) ./ size(A, 2)
println("Pune's rain minus its mean: ", round.(anomaly[1, :]; digits = 2))
for (k, layer) in enumerate(eachslice(A; dims = 3))  # one view per variable
    println("slice ", k, ": size ", size(layer), ", a view: ", layer isa SubArray, ", maximum ", maximum(layer))
end
println("size(mapslices(maximum, A; dims = (1, 2))) = ", size(mapslices(maximum, A; dims = (1, 2))))

from session import files, run

# A station × day × variable array reduced to station totals, wettest days and anomalies, then sliced.
files("station_cube.jl")
run("""
julia station_cube.jl
""")

$ julia station_cube.jl
size(A) = (3, 6, 2)
total rain per station: [64.5, 83.0, 23.0]
mean maximum temperature per station: [29.58, 31.67, 35.67]
Pune    wettest on day 3, 30.5 mm
Chennai wettest on day 4, 41.0 mm
Delhi   wettest on day 2, 11.0 mm
Pune's rain minus its mean: [1.25, -10.75, 19.75, -6.75, 7.25, -10.75]
slice 1: size (3, 6), a view: true, maximum 41.0
slice 2: size (3, 6), a view: true, maximum 37.0
size(mapslices(maximum, A; dims = (1, 2))) = (1, 1, 2)

Step down the reduction from the 3×6×2 array to the station totals:

cat(rain, tmax; dims = 3) stacked the two layers. A[:, :, 1] took the rainfall layer, sum(…; dims = 2) totalled each station’s six days, and findmax(…; dims = 2) found each station’s wettest day: day 3 for Pune with 30.5 mm, day 4 for Chennai with 41.0 mm and day 2 for Delhi with 11.0 mm. Subtracting the 3×1 matrix of station means from the 3×6 layer gave each day’s departure from its station’s mean, a broadcast that works because dims kept the dimension.

eachslice(A; dims = 3) handed out the two 3×6 layers, and its documentation says it returns views, as the run’s true confirmed. mapslices(maximum, A; dims = (1, 2)) applies a function to every slice and assembles the results into a new array, here of size (1, 1, 2). Its documentation explains that it passes each slice as a copy, which is why the function may change a slice without affecting A, and why eachslice is the cheaper choice inside a loop.

Practise: summaries of a station cube

The self-check works through the reductions and shapes of this lesson on matrices it has not shown yet, beginning with the idea you have answered least well; each answer explains itself, the wrong ones by the slip they come from:

Interview questions

Warm-up (fresher to 2–5 years): why does sum(M; dims = 2) return a matrix, and how do you get a vector?

sum(M; dims = 2) adds along the second dimension, the columns of each row, and keeps that dimension in the result with length 1, so a 2×3 M gives a 2×1 Matrix, printed [6; 15;;]. Keeping it is deliberate: the result lines up with M, so M ./ sum(M; dims = 2) divides each row by its own total in one broadcast. When a plain vector is needed, vec(sum(M; dims = 2)) or dropdims(sum(M; dims = 2); dims = 2) gives [6, 15], and sum.(eachrow(M)) builds the vector directly. The same holds for maximum, minimum, count, mean and mapreduce with dims.

Two questions take the shapes of this lesson into real data: numbers that come in from a row-major program, and a complex matrix whose phases change sign. Commit to an answer before you compare it with each first step:

Key takeaways

• dims names the dimension a reduction adds up, and the result keeps it with length 1: sum(M; dims = 2) of a 2×3 matrix is a 2×1 matrix, and vec or dropdims makes it a vector.
• argmax, findmax and extrema give positions and extremes, row by row with dims, and cumsum and diff need dims on a matrix.
• hcat and vcat join side by side and end to end, cat along any dimension, and stack or reduce(hcat, vs) turns many vectors into one matrix.
• reshape keeps the storage order, so it never transposes; permutedims moves the elements, and insertdims (Julia 1.12) and dropdims add and remove dimensions of length 1.
• A' is a lazy, conjugating view of A, transpose a view that only swaps, and permutedims or copy gives an independent matrix.
• eachslice hands out views of the slices of a 3-D array, while mapslices copies each slice and assembles the results into a new array.