Dense data are stored contiguously in memory, addressed by a single index (the memory address). Array memory ordering schemes translate that single index into multiple indices corresponding to the array coordinates. For example, matrices have two indices: rows and columns. Three-d arrays have three, and so on.
Column-major order is used by Fortran, Matlab, R, and most underlying core linear algebra libraries (BLAS). Sequential address locations are translated into array coordinates i, j, k, … so that the first array coordinates vary most rapidly with address, the next array coordinates less rapidly, and so on. For instance, four address locations 1, 2, 3, 4 are translated into a two by two matrix coordinate system (i, j) as:
address i j
1 1 1
2 2 1
3 1 2
4 2 2The phrase column-major comes from the matrix example, where sequentially addressed data are laid out sequentially along columns of the matrix.
The concept of “row” and “column” don’t directly apply to n-d arrays, but the same idea holds. For instance the R language lays out sequential addresses from 1, 2, …, 8, into a 2x2x2 3-d array as:
address i j k
1 1 1 1
2 2 1 1
3 1 2 1
4 2 2 1
5 1 1 2
6 2 1 2
7 1 2 2
8 2 2 2Consider the 3-d case shown above. Given array dimensions d1=2, d2=2, d3=2, a formula that takes 1-based coordinates i, j, k and returns address location a is
a = i + (j - 1) * d1 + (k - 1) * d2 * d1.
Row-major ordering is a different translation between sequential address indices and array coordinates, instead laying sequential data in memory across rows in matrices. Row-major ordering is sometimes called “C” style ordering and column-major ordering “Fortran” style. Python/NumPy refers to the orderings in array flags as C_CONTIGUOUS and F_CONTIGUOUS, respectively. For instance address locations 1, 2, 3, 4 are translated into a 2x2 matrix coordinate system (i, j) as:
address i j
1 1 1
2 1 2
3 2 1
4 2 2Efficient wrappers to BLAS routines exist for row-major ordered arrays. For completeness, here is a 2x2x2 3-d example layout:
address i j k
1 1 1 1
2 1 1 2
3 1 2 1
4 1 2 2
5 2 1 1
6 2 1 2
7 2 2 1
8 2 2 2And similarly to above a formula for this example that converts these 1-based array coordinates to address indices is:
a = k + (j - 1) * d3 + (i - 1) * d3 * d2.
See the following notes for a general formula for row- and column-order coordinate to address mapping, but note the use of zero-based indexing.
The Python NumPy library is very general. It can use either row-major or column-major ordered arrays, but it defaults to row-major ordering. NumPy also supports sophisticated views of data with custom strides across non-contiguous regions of memory.
R displays array data with unambiguously-labeled coordinate indices. Python doesn’t show this and displays n-d array data in different order than R (making matters somewhat confusing for R users). Consider the following example that creates and displays identical 4x3x2 arrays in R and Python:
array(1:24, c(4,3,2))
## , , 1
## [,1] [,2] [,3]
## [1,] 1 5 9
## [2,] 2 6 10
## [3,] 3 7 11
## [4,] 4 8 12
##
## , , 2
## [,1] [,2] [,3]
## [1,] 13 17 21
## [2,] 14 18 22
## [3,] 15 19 23
## [4,] 16 20 24import numpy as np
np.reshape(np.arange(1,25), (4,3,2), "F")
## array([[[ 1, 13],
## [ 5, 17],
## [ 9, 21]],
##
## [[ 2, 14],
## [ 6, 18],
## [10, 22]],
##
## [[ 3, 15],
## [ 7, 19],
## [11, 23]],
##
## [[ 4, 16],
## [ 8, 20],
## [12, 24]]])It’s easier to know which coordinates go where in R because they are labeled. Python, using column-major ordering, displays the same thing but in a different order where the first indices are grouped together in order. To see that these arrays are, in fact, the same, let’s pick out values along just the first “row”, that is values with a first index of 1 (R) or 0 (Python):
array(1:24, c(4, 3, 2))[1,, ,drop=FALSE]
## , , 1
## [,1] [,2] [,3]
## [1,] 1 5 9
##
## , , 2
## [,1] [,2] [,3]
## [1,] 13 17 21np.reshape(np.arange(1,25), (4,3,2), "F")[0]
## array([[ 1, 13],
## [ 5, 17],
## [ 9, 21]])I specified R’s drop=FALSE argument to preserve array dimensionality above. If we use drop=TRUE (the default) then R returns a 3x2 array in column-major order–exactly the same result as Python above.
array(1:24, c(4, 3, 2))[1,, ,drop=TRUE]
## [,1] [,2]
## [1,] 1 13
## [2,] 5 17
## [3,] 9 21Note that the Python result is a special view of the original array data, not a copy. In this case it’s not stored in contiguous memory addresses and is neither really row- nor column-major. This is shown in the array flags:
np.reshape(np.arange(1,25), (4,3,2), "F")[0].flags
## C_CONTIGUOUS : False
## F_CONTIGUOUS : False
## OWNDATA : False
## WRITEABLE : True
## ALIGNED : True
## UPDATEIFCOPY : FalseThe reticulate package lets us easily mix R and Python code and data. Recall that R represents all dense arrays in column-major order but Python/NumPy can represent dense arrays much more generally. That difference warrants attention and can easily lead to confusion!
Remember the following things when working with R and Python arrays, especially n-d arrays with n > 2.
Also worth knowing:
Point number 3. introduces the most potential confusion. Let’s look at some examples to explore these points.
The following example creates a 2x2x2 array in Python using native NumPy row-major ordering and imports it into R. Despite the fact that they print out differently, they are in fact the same.
library(reticulate)
np <- import("numpy", convert=FALSE)
(x <- np$arange(1, 9)$reshape(2L, 2L, 2L))
## [[[ 1. 2.]
## [ 3. 4.]]
##
## [[ 5. 6.]
## [ 7. 8.]]]
(y <- py_to_r(x))
## , , 1
## [,1] [,2]
## [1,] 1 3
## [2,] 5 7
##
## , , 2
## [,1] [,2]
## [1,] 2 4
## [2,] 6 8Wait a minute! They look different! But remember Python’s print order is different. The first “rows” (first index values) are grouped together. Let’s pull out elements with first index of 1 in the R result, with and without dropping the unused dimension to show precisely what we’re indexing here:
y[1,,, drop=FALSE]
## , , 1
## [,1] [,2]
## [1,] 1 3
##
## , , 2
## [,1] [,2]
## [1,] 2 4
y[1,,, drop=TRUE]
## [,1] [,2]
## [1,] 1 2
## [2,] 3 4Note that this is the same as the first block printed in the Python output above! These arrays really are the same in Python and R, respectively. Their apparent differences are merely a result of printing.
Let’s look at this again with another example, this time with an array with different lengths along each dimension to make things even more clear (hopefully). Consider the following 4x3x2 array constructed in Python in row-major order:
np <- import("numpy", convert=FALSE)
(x <- np$reshape(np$arange(1, 25), c(4L, 3L, 2L)))
## [[[ 1. 2.]
## [ 3. 4.]
## [ 5. 6.]]
##
## [[ 7. 8.]
## [ 9. 10.]
## [ 11. 12.]]
##
## [[ 13. 14.]
## [ 15. 16.]
## [ 17. 18.]]
##
## [[ 19. 20.]
## [ 21. 22.]
## [ 23. 24.]]]
(y <- py_to_r(x))
## , , 1
## [,1] [,2] [,3]
## [1,] 1 3 5
## [2,] 7 9 11
## [3,] 13 15 17
## [4,] 19 21 23
##
## , , 2
## [,1] [,2] [,3]
## [1,] 2 4 6
## [2,] 8 10 12
## [3,] 14 16 18
## [4,] 20 22 24Again, they look quite different but the R and Python arrays are really the same. Let’s pick out the sub-array with third index = 0 (Python), equivalently the third index = 1 in R.
np$take(x, 0L, 2L)
## [[ 1. 3. 5.]
## [ 7. 9. 11.]
## [ 13. 15. 17.]
## [ 19. 21. 23.]]
y[, , 1]
## [,1] [,2] [,3]
## [1,] 1 3 5
## [2,] 7 9 11
## [3,] 13 15 17
## [4,] 19 21 23The NumPy take() function is equivalent in this example to the Python notation x[:, :, 0]; that is, entries with third dimension index = 0. See https://numpy.org/doc/stable/reference/generated/numpy.take.html for more information.
The corresponding R notation, y[, , 1], returns the same thing in this example: a 4x3 matrix.
Despite different internal memory ordering, and particularly despite awkward differences in printing arrays, the arrays are the same and are indexed the same way in each language respectively.
The previous examples focused on row-major arrays natively constructed in Python. Let’s see what happens when we start with column-major arrays from R and work with them in Python.
(y <- array(1:24, c(4, 3, 2))) # In R
## , , 1
## [,1] [,2] [,3]
## [1,] 1 5 9
## [2,] 2 6 10
## [3,] 3 7 11
## [4,] 4 8 12
##
## , , 2
## [,1] [,2] [,3]
## [1,] 13 17 21
## [2,] 14 18 22
## [3,] 15 19 23
## [4,] 16 20 24
(x <- np$array(y)) # and now in Python
## [[[ 1 13]
## [ 5 17]
## [ 9 21]]
##
## [[ 2 14]
## [ 6 18]
## [10 22]]
##
## [[ 3 15]
## [ 7 19]
## [11 23]]
##
## [[ 4 16]
## [ 8 20]
## [12 24]]]Note that the Python version takes advantage of NumPy’s extraordinary flexibility and preserves R’s column-major ordering:
x$flags
## C_CONTIGUOUS : False
## F_CONTIGUOUS : True
## OWNDATA : True
## WRITEABLE : True
## ALIGNED : True
## UPDATEIFCOPY : FalseYou can probably tell already from the previous sections that these arrays are the same, and obey the same indexing conventions. The next example selects a subarray such that the third index of each array is 0 (Python) or 1 (R):
y[, , 1]
## [,1] [,2] [,3]
## [1,] 1 5 9
## [2,] 2 6 10
## [3,] 3 7 11
## [4,] 4 8 12
np$take(x, 0L, 2L)
## [[ 1 5 9]
## [ 2 6 10]
## [ 3 7 11]
## [ 4 8 12]]It’s important to remember that the order is preserved from Python when copying an array result back into R:
py_to_r(np$take(x, 0L, 2L))
## [,1] [,2] [,3]
## [1,] 1 5 9
## [2,] 2 6 10
## [3,] 3 7 11
## [4,] 4 8 12The upshot is that arrays created by either R or Python are indexed exactly the same in either language.
Right. That’s just a simple consequence of the default column-major and row-major formats used in R and NumPy respectively. You are always free to use R’s column-major format directly in Python, for example using the “F” flag below (for Fortran):
np$reshape(np$arange(1, 25), c(4L, 3L, 2L), "F")
## [[[ 1. 13.]
## [ 5. 17.]
## [ 9. 21.]]
##
## [[ 2. 14.]
## [ 6. 18.]
## [ 10. 22.]]
##
## [[ 3. 15.]
## [ 7. 19.]
## [ 11. 23.]]
##
## [[ 4. 16.]
## [ 8. 20.]
## [ 12. 24.]]]Note that the result is just like the one we got starting from R above.
Re-arranging R arrays into row-major order requires more work. R is less flexible than Python and we can’t explicitly change R’s memory order representation. When the array is a matrix when we can simply use byrow=TRUE. In the n-d array case, a portion of the problem can be reduced to using byrow=TRUE followed by judicious index permutation with aperm(). Here is one somewhat inefficient example:
y <- aperm(array(matrix(1:24, c(3 * 4, 2), byrow=TRUE),
c(3, 4, 2)), c(2, 1, 3))See the last section below for a different example.
We can verify that the above ugly expression exactly reproduces a NumPy row-major array by subtracting our R array from a native Python one:
np <- import("numpy", convert=FALSE)
o <- import("operator", convert=FALSE)
o$sub(np$arange(1, 25)$reshape(4L, 3L, 2L), np$array(y))
## [[[ 0. 0.]
## [ 0. 0.]
## [ 0. 0.]]
##
## [[ 0. 0.]
## [ 0. 0.]
## [ 0. 0.]]
##
## [[ 0. 0.]
## [ 0. 0.]
## [ 0. 0.]]
##
## [[ 0. 0.]
## [ 0. 0.]
## [ 0. 0.]]]The above NumPy arrays are the same, their element-wise difference is zero.
In R you would typically reshape an array using the dim<-() function. For example:
dim(x) <- c(1000, 28, 28)In R, this operation simply changes the “dim” attribute of the array, effectively re-interpreting the array indices as specified using column-major semantics.
However, the NumPy reshape method uses row-major semantics by default, so if you are mixing R and Python code that reshapes arrays you will find that the reshaping will be inconsistent if you use the R dim<-() function.
To overcome this inconsistency, there is an array_reshape() function which will reshape an R array using row-major semantics (i.e. will fill the new dimensions in row-major rather than col-major order). The example above would be re-written as:
x <- array_reshape(x, c(1000, 28, 28))Here’s a further example to illustrate the difference:
# let's construct a 2x2 array from a vector of 4 elements
x <- 1:4
# rearrange will fill the array row-wise
array_reshape(x, c(2, 2))
# [,1] [,2]
# [1,] 1 2
# [2,] 3 4
# setting the dimensions 'fills' the array col-wise
dim(x) <- c(2, 2)
x
# [,1] [,2]
# [1,] 1 3
# [2,] 2 4It’s worth noting that analogs of R’s apply() function in Python behave differently. The following excellent Mathesaurus reference http://mathesaurus.sourceforge.net/r-numpy.html applies well to matrices and vectors, but is misleading for n-d arrays with n > 2.
In particular, Mathesaurus says that if a is a matrix, then the sum of each column in Python may be computed by a.sum(0), and in R (among other possible ways) by apply(a, 2, sum). Although correct for matrices, this is in general not quite right. A more precise R analog of NumPy’s a.sum(0) is apply(a, seq_along(dim(a))[-1], sum). In other words, a.sum(0) means sum over the first dimension, returning an array of the same dimensions as a but with the first dimension removed.
It’s easy to be confused by this, so let’s see an example using a 4x3x2 array, first in Python:
library(reticulate)
np <- import("numpy", convert=FALSE)
x <- np$arange(1, 25)$reshape(c(4L, 3L, 2L))
x$sum(0) # N. B. a 3x2 matrix!
## [[ 40. 44.]
## [ 48. 52.]
## [ 56. 60.]]
# N. B. A tuple() object is required here (NumPy vectors won't work)
x$sum(tuple(1L, 2L))
## [ 21. 57. 93. 129.]And now the corresponding sums in R:
y <- py_to_r(x)
apply(y, dim(y)[-1], sum)
## [,1] [,2] [,3]
## [1,] 40 48 56
## [2,] 44 52 60
apply(y, 1, sum)
## [1] 21 57 93 129These notes were prepared in response to a tensorflow issue now in the reticulate package https://github.com/rstudio/reticulate/issues/9. The issue directly gets to a common source of confusion with n-d arrays in R and Python and how they are printed and stored. A lightly-edited reproduction of the reference Python code in the issue appears below.
library(tensorflow)
np <- import("numpy", convert=FALSE)
a <- np$arange(1, 9)$reshape(c(2L, 2L, 2L))
b <- np$arange(1, 5)$reshape(c(2L, 2L, 1L))
c <- tf$matmul(tf$constant(a), tf$constant(b))
tf$Session()$run(c)
## , , 1
## [,1] [,2]
## [1,] 5 11
## [2,] 39 53The issue goes on to reproduce the example using R-generated arrays as follows:
A <- list(matrix(1:4, nrow=2, byrow=T), matrix(5:8, nrow=2, byrow=T))
A <- array(unlist(A), dim=c(2,2,2))However, already at this point we see that the R-generated array A is not the same as the above array a by comparing a with np$array(A) below.
However, we can see how it can be easy to make the mistake that they are the same simply because of the way the arrays are printed! The R array looks superficially the same as the printed Python array.
print(a)
## [[[ 1. 2.]
## [ 3. 4.]]
##
## [[ 5. 6.]
## [ 7. 8.]]]
print(np$array(A))
## [[[1 5]
## [2 6]]
##
## [[3 7]
## [4 8]]]
print(A)
## , , 1
## [,1] [,2]
## [1,] 1 2
## [2,] 3 4
##
## , , 2
## [,1] [,2]
## [1,] 5 6
## [2,] 7 8Instead, we need to construct the R array A differently to match the row-major order of Python, discussed in the previous sections. We can use many approaches including:
(A <- np$array(aperm(array(1:8, c(2,2,2)), c(3,2,1))))
## [[[1 2]
## [3 4]]
##
## [[5 6]
## [7 8]]]With similar care ordering the values in the b array we can finish replicating the example in R (with the same result as the reference Python example above).
A <- np$array(aperm(array(1:8, c(2,2,2)), c(3,2,1)))
B <- np$array(aperm(array(1:4, c(2,2,1)), c(2,1,3)))
C <- tf$matmul(tf$constant(A), tf$constant(B))
tf$Session()$run(C)
## , , 1
## [,1] [,2]
## [1,] 5 11
## [2,] 39 53