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 2
The 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 2
Consider 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 2
Efficient 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 2
And 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 24
import 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 ® 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 21
np.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 21
Note 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 : False
The 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 8
Wait 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 4
Note 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 24
Again, 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 23
The NumPy take() function is equivalent in this example to the Python
notation x[:, :, 0]; that is, entries with third dimension index = 0. (See
https://docs.scipy.org/doc/numpy/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 : False
You 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 ®:
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 12
The 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 4
It'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 129
These 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 53
The 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 8
Instead, 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