Julia for MATLAB, NumPy and R users

A loop that adds up a 4000×4000 matrix row by row, the order that suits NumPy’s row-majorrow-major order: Storing a matrix row by row, so that the last index changes fastest from one memory slot to the next. NumPy arrays use it unless told otherwise; Julia, MATLAB and R store column by column. arrays, took several times as long in Julia as the same loop column by column, and no message said why. A port from MATLAB, NumPy or R can go wrong in this way without any error: the code runs, and the order, the copies or the arithmetic are quietly different.

This lesson collects those differences in one place, each shown by a run recorded on Julia 1.13.1: where indices start and how matrices are stored, which assignments copy and which share, what *, ^ and ' mean, when a loop beats vectorised code and when a library call beats both, how characters, ranges and missing values differ, and twenty everyday idioms translated and run.

Julia counts from 1 and stores by column, as MATLAB and R do

Two choices decide how array code is written, the first index and the order of storage, and Julia makes the same two as MATLAB and R:

The sources agree. Julia’s list of differences says that Julia indexes from 1 and keeps its arrays column-majorcolumn-major order: Storing a matrix column by column, so that the first index changes fastest from one memory slot to the next: [1 2 3; 4 5 6] is stored as 1, 4, 2, 5, 3, 6. Julia stores every Array this way. More in “Arrays: vectors, matrices and column-major storage”, in the order Fortran uses, while NumPy keeps them in the order of C, row by row, unless told otherwise. MathWorks says MATLAB stores elements with a column-major layout by default, and R’s manual says its arrays are filled in the order Fortran uses, with the first subscript moving fastest. Step through the same 2×3 matrix in both layouts:

The order decides which loop is fast. This program adds up a 250 × 250 matrix and a 4000 × 4000 one, each twice, with the row index innermost and then with the column index innermost, and keeps the best of several runs of each:

# The same sum over a small and a large matrix, column by column and then row by row.
function sum_by_columns(A)
    s = 0.0
    for j in axes(A, 2)
        for i in axes(A, 1)          # down a column: the order Julia stores A in
            s += A[i, j]
        end
    end
    return s
end
function sum_by_rows(A)
    s = 0.0
    for i in axes(A, 1)
        for j in axes(A, 2)          # along a row: each step jumps a whole column ahead
            s += A[i, j]
        end
    end
    return s
end
function best_time(f, A, runs)
    t = Inf
    for _ in 1:runs
        t = min(t, @elapsed(s = f(A)))
        s == -1 && println("never")  # uses the result, so the call cannot be skipped
    end
    return t
end
for n in (250, 4000)
    A = [Float64((i + 2j) % 7) for i in 1:n, j in 1:n]
    same = sum_by_columns(A) == sum_by_rows(A)
    runs = n < 1000 ? 50 : 5
    tc = best_time(sum_by_columns, A, runs)
    tr = best_time(sum_by_rows, A, runs)
    println(n, "×", n, " (", round(sizeof(A) / 2^20; digits = 1), " MiB): same sum ", same, "; by columns ",
            round(tc * 1000; digits = 2), " ms, by rows ", round(tr * 1000; digits = 2), " ms, ratio ", round(tr / tc; digits = 1))
end

from session import files, run

# One sum over a 250×250 and a 4000×4000 matrix in the two loop orders, each timed as the best of several runs.
files("row_vs_column.jl")
run("""
julia row_vs_column.jl
""")

$ julia row_vs_column.jl
250×250 (0.5 MiB): same sum true; by columns 0.06 ms, by rows 0.05 ms, ratio 1.0
4000×4000 (122.1 MiB): same sum true; by columns 16.23 ms, by rows 50.11 ms, ratio 3.1

Both orders gave the same sums. On the small matrix, half a megabyte, they took about the same time; on the large one, 122 MiB, the row-by-row loop, the natural order for NumPy, took several times as long, as the ratios in the run show. Each of its steps jumps 4,000 slots ahead in memory, while the column order reads the slots one after another, and the run shows the cost only once the matrix is large. Julia’s performance tips turn this into a rule of thumb: in a column-major array, the loop over the first index belongs innermost.

Assignment shares and slices copy, unlike MATLAB and NumPy

The second set of differences is about copies, and here each language made its own choice:

# Habits from MATLAB and NumPy that change other data in Julia, or leave it unchanged.
A = [1.0 2.0; 3.0 4.0]
B = A                    # MATLAB would give B values of its own; here B is a second name
B[1, 1] = 0.0
println("after B = A and B[1, 1] = 0: A = ", A)
A = [1.0 2.0; 3.0 4.0]
B = copy(A)              # an independent copy
B[1, 1] = 0.0
println("after B = copy(A) and B[1, 1] = 0: A = ", A)
col = A[:, 1]            # in NumPy a basic slice is a view; here it is a copy
col .= 0.0
println("after col = A[:, 1] and col .= 0: A = ", A)
col = @view A[:, 1]      # a view, for when writing through it is the point
col .= 0.0
println("after col = @view A[:, 1] and col .= 0: A = ", A)
x = [1, 2, 3]
y = x
x += [1, 1, 1]           # x = x + [1, 1, 1]: a new vector, now named x
println("after y = x and x += [1, 1, 1]: x = ", x, ", y = ", y)
y = x
x .+= 1                  # .+= writes into the vector that x and y both name
println("after y = x and x .+= 1: x = ", x, ", y = ", y)

from session import files, run

# Assignment, slices, views and updating operators, each followed by the data it did or did not change.
files("aliasing_bug.jl")
run("""
julia aliasing_bug.jl
""")

$ julia aliasing_bug.jl
after B = A and B[1, 1] = 0: A = [0.0 2.0; 3.0 4.0]
after B = copy(A) and B[1, 1] = 0: A = [1.0 2.0; 3.0 4.0]
after col = A[:, 1] and col .= 0: A = [1.0 2.0; 3.0 4.0]
after col = @view A[:, 1] and col .= 0: A = [0.0 2.0; 0.0 4.0]
after y = x and x += [1, 1, 1]: x = [2, 3, 4], y = [1, 2, 3]
after y = x and x .+= 1: x = [3, 4, 5], y = [3, 4, 5]

B = A followed by B[1, 1] = 0 changed A, because B is a second name for the same matrix; copy(A) gave an independent one. Julia’s list of differences states it plainly: giving an array a second name copies nothing. MATLAB behaves otherwise. MathWorks’ page on copies explains that MATLAB has no references to values: an array is passed by value, and a shared copy is separated the first time one side changes, which it calls copy-on-writecopy-on-write: Sharing one copy of data between names until one of them is changed, and copying it only then. MATLAB treats arrays this way, so after B = A a change to B leaves A as it was; Julia does not copy at all., so changing B there leaves A alone.

Slices go the other way. col = A[:, 1] followed by col .= 0 left A unchanged, because a Julia slice is a copy, while NumPy’s documentation says basic indexing always creates views, so the same line in NumPy would have cleared the column. In Julia, @view A[:, 1] is the view, and writing through it cleared the column. Updating operators differ as well: the list of differences notes that NumPy’s += works in place and Julia’s does not, and the run showed it, as x += [1, 1, 1] made a new vector and left y with the old values, while x .+= 1 wrote into the vector that both names reach.

A star multiplies matrices, and a dot works element by element

Operators are where a line-by-line port changes the arithmetic without an error. Julia keeps the meanings of linear algebra for *, ^ and ', as MATLAB does, and asks for a dot for element-wise work:

# The symbols that mean different things in Julia, MATLAB, NumPy and R, with what Julia computes.
A = [1 2; 3 4]
B = [5 6; 7 8]
println("A * B  = ", A * B)
println("A .* B = ", A .* B)
println("A^2    = ", A^2)
println("A .^ 2 = ", A .^ 2)
C = [1+2im 3; 0 1]
println("C'     = ", C')
S = [2.0 1.0; 1.0 3.0]
b = [3.0, 5.0]
println("S \\ b  = ", S \ b, ", and S * (S \\ b) = ", S * (S \ b))
println("7 / 2 = ", 7 / 2, ", 7 ÷ 2 = ", 7 ÷ 2)

from session import files, run

# Matrix and element-wise products and powers, the adjoint, a linear solve and the two divisions.
files("operators.jl")
run("""
julia operators.jl
""")

$ julia operators.jl
A * B  = [19 22; 43 50]
A .* B = [5 12; 21 32]
A^2    = [7 10; 15 22]
A .^ 2 = [1 4; 9 16]
C'     = Complex{Int64}[1 - 2im 0 + 0im; 3 + 0im 1 + 0im]
S \ b  = [0.8, 1.4], and S * (S \ b) = [3.0, 5.0]
7 / 2 = 3.5, 7 ÷ 2 = 3

A * B was the matrix product [19 22; 43 50] and A .* B the element-wise one, A^2 the matrix power [7 10; 15 22] and A .^ 2 each element squared. C' conjugated as well as transposed, and S \ b solved the system S x = b, which S * (S \ b) confirmed. 7 / 2 gave the float 3.5 and 7 ÷ 2 the integer 3, as the lesson on integer division showed.

The NumPy rows come from NumPy’s own guide for MATLAB users, and the R rows from its manual, which defines A * B as the element by element product and A %*% B as the matrix product. The trap is the *: in NumPy and R it multiplies element by element, in Julia it is the matrix product, so a port that keeps the * computes something else, or fails when the shapes do not fit a matrix product.

Loops are compiled, so vectorise for clarity or for a library

Julia’s list of differences says that in R, performance requires vectorization, and its Python section that in Julia loops are fast, with no need to write vectorised code for performance. The first lesson of this track measured a plain loop against sum. Here the question is the other one: when does a vectorised call still pay?

# A vectorised expression against a loop, and a library matrix product against loops.
using LinearAlgebra
BLAS.set_num_threads(1)                         # one thread, so the comparison is about the method
x = collect(range(0, 1; length = 1_000_000))
w = collect(range(1, 2; length = 1_000_000))
vectorised(x, w) = sum(x .^ 2 .* w) / sum(w)    # the MATLAB or NumPy habit: a temporary array
function looped(x, w)
    num = 0.0
    den = 0.0
    for i in eachindex(x, w)
        num += x[i]^2 * w[i]
        den += w[i]
    end
    return num / den
end
vectorised(x, w); looped(x, w)
println("weighted mean of squares: ", round(vectorised(x, w); digits = 6), " and ", round(looped(x, w); digits = 6))
println("bytes allocated: vectorised ", @allocated(vectorised(x, w)), ", loop ", @allocated(looped(x, w)))
function matmul_loops(A, B)
    C = zeros(size(A, 1), size(B, 2))
    for j in axes(B, 2), k in axes(A, 2), i in axes(A, 1)   # i innermost, down the columns
        C[i, j] += A[i, k] * B[k, j]
    end
    return C
end
n = 600
A = [sin(i + j) for i in 1:n, j in 1:n]
B = [cos(i - j) for i in 1:n, j in 1:n]
println("same product: ", matmul_loops(A, B) ≈ A * B)
best(f) = minimum(@elapsed(f(A, B)) for _ in 1:3)
println("A * B was ", round(best(matmul_loops) / best(*); digits = 1), " times as fast as the loops, on one thread")

from session import files, run

# A vectorised expression against a loop, and BLAS's matrix product against three nested loops.
files("vectorise_or_loop.jl")
run("""
julia vectorise_or_loop.jl
""")

$ julia vectorise_or_loop.jl
weighted mean of squares: 0.388889 and 0.388889
bytes allocated: vectorised 8011840, loop 0
same product: true
A * B was 4.6 times as fast as the loops, on one thread

The weighted mean of squares came out the same both ways, 0.388889, but the vectorised habit, x .^ 2 .* w, built a temporary array of a million values, 8,011,840 bytes, which the loop did not need. The matrix product went the other way: A * B on Float64 matrices hands the work to BLAS, the optimised library behind LinearAlgebra’s products, and even on one thread it was several times faster than three nested loops written in the right order. Vectorise when it names a library operation or reads more clearly, and loop, or fuse with dots, when the formula is your own. Answer three questions about your own computation:

In short: a library call for products, solves and other operations that a library provides, a fused broadcast for a formula applied element by element, and a plain loop for everything that combines elements in its own way.

Characters, ranges and missing values differ in small ways

A handful of smaller differences catch ported code as well, and this program runs into each of them:

# Strings, ranges and missing values, where the habits of other languages mislead.
println(typeof('a'), " and ", typeof("a"), "; 'a' * 'b' = ", 'a' * 'b')
try
    Meta.parse("x = 'ab'")
catch err
    println(strip(last(split(sprint(showerror, err), "── "))))
end
r = 1:1_000_000
println(typeof(r), ": ", sizeof(r), " bytes for a million values; collect(1:5) = ", collect(1:5))
data = [21.5, missing, 22.4]
println("sum(data) = ", sum(data), ", sum(skipmissing(data)) = ", sum(skipmissing(data)))
println("missing == 1 is ", missing == 1, ", NaN == NaN is ", NaN == NaN)
try
    NA
catch err
    println(first(split(sprint(showerror, err), '\n')))
end

from session import files, run

# Characters and strings, a lazy range, missing values, and a name that R users expect.
files("traps.jl")
run("""
julia traps.jl
""")

$ julia traps.jl
Char and String; 'a' * 'b' = ab
character literal contains multiple characters
UnitRange{Int64}: 16 bytes for a million values; collect(1:5) = [1, 2, 3, 4, 5]
sum(data) = missing, sum(skipmissing(data)) = 43.9
missing == 1 is missing, NaN == NaN is false
UndefVarError: `NA` not defined in `Main`

Single quotes make a Char and double quotes a String, so 'ab' does not even parse: the message says a character literal contains multiple characters. The list of differences tells R users the same: in Julia, single quotes are for characters only. 1:1_000_000 is a UnitRange of 16 bytes, which computes its values when asked, as the manual says of a:b for MATLAB users, and collect turns a range into a vector. R’s NA does not exist in Julia, where missing values are the missing object: sum of a vector with a missing is missing, as R’s manual says any operation on an NA becomes NA, and skipmissing leaves the gaps out, which gave 43.9. A missing compared with anything is missing, and a NaN is not equal even to itself.

R’s manual describes negative indices as the values to be excluded, which is why x[-1] means all but the first in R and the last element in NumPy, and is a BoundsError in Julia.

Twenty idioms, translated and run

Here are twenty idioms as a Julia programmer writes them, each run once with its result, and below the program the forms they translate from:

# Twenty everyday idioms as Julia writes them, each with what it returned.
x = [4, 1, 7, 2]
A = [2.0 1.0; 1.0 3.0]
b = [3.0, 5.0]
idioms = [
    "zeros(3)" => zeros(3),
    "zeros(2, 2)" => zeros(2, 2),
    "range(0, 1; length = 5)" => range(0, 1; length = 5),
    "collect(0:0.25:1)" => collect(0:0.25:1),
    "x[end]" => x[end],
    "x[2:end]" => x[2:end],
    "length(x)" => length(x),
    "findall(>(2), x)" => findall(>(2), x),
    "x[x .> 2]" => x[x .> 2],
    "cumsum(x)" => cumsum(x),
    "sort(x; rev = true)" => sort(x; rev = true),
    "repeat([1, 2], 3)" => repeat([1, 2], 3),
    "reshape(1:6, 2, 3)" => reshape(1:6, 2, 3),
    "vec([1 2; 3 4])" => vec([1 2; 3 4]),
    "A \\ b" => A \ b,
    "A .^ 2" => A .^ 2,
    "permutedims([1 2; 3 4])" => permutedims([1 2; 3 4]),
    "size(A, 1)" => size(A, 1),
    "sum(A; dims = 2)" => sum(A; dims = 2),
    "string(\"run-\", 7)" => string("run-", 7),
]
for (k, (code, value)) in enumerate(idioms)
    println(lpad(k, 2), "  ", rpad(code, 25), repr(value))
end

from session import files, run

# The Julia side of twenty translated idioms, with what each returned.
files("translation_table.jl")
run("""
julia translation_table.jl
""")

$ julia translation_table.jl
 1  zeros(3)                 [0.0, 0.0, 0.0]
 2  zeros(2, 2)              [0.0 0.0; 0.0 0.0]
 3  range(0, 1; length = 5)  0.0:0.25:1.0
 4  collect(0:0.25:1)        [0.0, 0.25, 0.5, 0.75, 1.0]
 5  x[end]                   2
 6  x[2:end]                 [1, 7, 2]
 7  length(x)                4
 8  findall(>(2), x)         [1, 3]
 9  x[x .> 2]                [4, 7]
10  cumsum(x)                [4, 5, 12, 14]
11  sort(x; rev = true)      [7, 4, 2, 1]
12  repeat([1, 2], 3)        [1, 2, 1, 2, 1, 2]
13  reshape(1:6, 2, 3)       [1 3 5; 2 4 6]
14  vec([1 2; 3 4])          [1, 3, 2, 4]
15  A \ b                    [0.8, 1.4]
16  A .^ 2                   [4.0 1.0; 1.0 9.0]
17  permutedims([1 2; 3 4])  [1 3; 2 4]
18  size(A, 1)               2
19  sum(A; dims = 2)         [3.0; 4.0;;]
20  string("run-", 7)        "run-7"

The numbers in the tables are the lines of the run. The MATLAB and NumPy forms follow NumPy’s guide for MATLAB users and MathWorks’ page on zeros, and the R forms R’s own manual:

Line 13 is the one to read twice: reshape(1:6, 2, 3) gave [1 3 5; 2 4 6], filled by columns as R’s matrix and MATLAB’s reshape fill it, while NumPy’s reshape fills by rows unless it is given order='F', as its guide notes. Line 19 differs in its numbering: NumPy counts axes from 0, so its axis=1, the columns of each row, is Julia’s dims = 2, and Julia keeps the reduced dimension, which NumPy’s sum removes unless it is given keepdims=True. Packages differ too: plotting and data frames come from packages in Julia rather than from the language, and the scientific-computing module of this track covers them.

Practise: port six snippets

The module’s last self-check mixes the translations of this lesson. Its first questions come from whichever idea you know least well so far, and each answer explains itself in the terms of the language the habit came from:

Interview questions

Warm-up (fresher to 2–5 years): name three differences that break a line-by-line port of NumPy array code to Julia.

Indices start at 1, so a[0] is an error and a[-1] must become a[end]. Storage is column-major, where NumPy’s default is row-major, so a loop with the column index innermost, fast in NumPy, jumps through memory in Julia: summing a 4000×4000 matrix row by row took several times as long as column by column. And the copy rules are reversed for slices: NumPy’s basic slices are views, so writing into a[:, 0] changes a, while Julia’s A[:, 1] is a copy and the write is lost; @view gives the view. Two more follow from the operators: NumPy’s * is element-wise and Julia’s is the matrix product, and NumPy’s += updates in place where Julia’s builds a new array unless written .+=.

Two porting problems follow on from here: a MATLAB function that corrupts its caller’s data, and a NumPy loop that runs several times slower in Julia. Note the first thing you would measure or print in each case before you open its first step:

Key takeaways

• Julia, MATLAB and R count from 1 and store matrices by column; NumPy counts from 0 and stores by row by default, so ported loops need the first index innermost.
• B = A shares the array in Julia, as in NumPy and unlike MATLAB, while A[:, 1] copies, as in MATLAB and unlike NumPy; copy and @view say which one you mean.
• *, ^ and ' keep their linear-algebra meanings, a dot makes any operation element-wise, and NumPy’s and R’s * is Julia’s .*.
• Loops and fused broadcasts are compiled, so vectorised code with temporaries gains no speed and costs memory, while a library call such as A * B still beats hand-written loops.
• 'a' is a Char, 1:5 is a lazy range, missing replaces NA, and skipmissing leaves gaps out of a sum.
• The twenty translated idioms run as shown; reshape fills by columns and dims = 2 is NumPy’s axis=1.