Inconsistency between tf.matmul in Python and tf$matmul in R about reticulate HOT 35 CLOSED

I don't know numpy and R matrix functions well enough to say why, but I can tell you that tf$constant in R is literally call Python's tf.constant (there is no intermediate code which transforms the input and/or would cause the functions to behave differently.

The numpy constants you are creating are:

[1 2 3 4 5 6 7 8]
[1 2 3 4]

The R matrices you are creating are (when converted to NumPy):

[[[1 5]
  [3 7]]

 [[2 6]
  [4 8]]]

[[[1]
  [3]]

So I think the R matrices are a red-herrring (you don't need to be doing that much munging around with matrices!). If we simplify your code to this it yields the expected result:

library(tensorflow)
A_tf = tf$constant(c(1:8), dtype="float64", shape=c(2,2,2))
B_tf = tf$constant(c(1:4), dtype="float64", shape=c(2,2,1))
sess = tf$Session()
sess$run(A_tf)
sess$run(B_tf)
sess$run(tf$matmul(A_tf,B_tf))

from reticulate.

From @bwlewis on February 14, 2017 14:2

No, it's not quite the same B_tf above, and consequentially the result
is a matrix instead of an array of vectors as desired.

R code:

library(tensorflow)
import("numpy") -> np
tf$Session()$run(tf$constant(np$arange(1, 5, dtype=np$int32), shape=c(2,2,1)))

#, , 1
#     [,1] [,2]
#[1,]    1    2
#[2,]    3    4

Python equivalent:

tf.Session().run(tf.constant(np.arange(1, 5, dtype=np.int32),shape=[2, 2, 1]) )

#array([[[1],
#        [2]],
#
#       [[3],
#        [4]]], dtype=int32)

There is indeed something funny going on here, possibly related to avoidance of native Python lists; trying to figure it out...

Also note, as @G-Lynn notes above:

tf$Session()$run(tf$constant(np$arange(1, 5, dtype=np$int32), shape=c(2,1,2)))
#, , 1
#
#     [,1]
#[1,]    1
#[2,]    3
#
#, , 2
#
#     [,1]
#[1,]    2
#[2,]    4

(sorry to keep editing this), https://github.com/bwlewis/python is also wrong! Curiously, the error is different--it's a transposition error instead of a dimensional error. Closer, but still wrong. Wow.

library(python)
 import("tensorflow") -> tf
 import("numpy") -> np
tf$Session()$run(tf$constant(np$arange(1, 5, dtype=np$int32), shape=c(2,2,1)))
#[[[1]
#  [2]]
#
# [[3]
#  [4]]]

from reticulate.

From @bwlewis on February 14, 2017 14:29

Type of shape is not defined in the doc: https://www.tensorflow.org/api_docs/python/constant_op/constant_value_tensors

bummer. But it seems to work with either a Python list or a numpy array:

>>> tf.Session().run(tf.constant(np.arange(1,5), shape=np.array([2,2,1])))
array([[[1],
        [2]],

       [[3],
        [4]]])
>>> tf.Session().run(tf.constant(np.arange(1,5), shape=[2,2,1]))
array([[[1],
        [2]],

       [[3],
        [4]]])

from reticulate.

The R shape function will yield a Python list: https://rstudio.github.io/tensorflow/using_tensorflow_api.html#tensor_shapes

Does that help?

from reticulate.

From @bwlewis on February 14, 2017 15:20

It doesn't, curiously enough:

tf$Session()$run(tf$constant(np$arange(1, 5, dtype=np$int32), shape=shape(2,2,1)))
#, , 1
#
#     [,1] [,2]
#[1,]    1    2
#[2,]    3    4

from reticulate.

Even explicitly passing an R array (which will become a numpy array on the other side) gives the same result:

tf$Session()$run(tf$constant(np$arange(1, 5, dtype=np$int32), shape=array(c(2,2,1))))

from reticulate.

This code yields the expected output:

py <- import("__builtin__")
py$print(tf$Session()$run(tf$constant(np$arange(1, 5, dtype=np$int32), shape=c(2,2,1))))

This means that the problem is NOT on the input side, it's on the marshaling of a 3-d array back into R. That code is here:

https://github.com/rstudio/reticulate/blob/master/src/python.cpp#L516

I do this to "cast" the NumPy array into a Fortran compatible array:

array = (PyArrayObject*)PyArray_CastToType(array, descr, NPY_ARRAY_FARRAY);

Then I just access the array data in memory order. Clearly this isn't working. I had assumed this was working properly based on this test (https://github.com/rstudio/reticulate/blob/master/tests/testthat/test-python-numpy.R#L18-L22) however looking closely at the test it's invalid because the result compared to also goes through the faulty marshaling layer.

from reticulate.

From @bwlewis on February 14, 2017 16:12

Aha. Note also a row/column memory ordering problem in the tensor:

import tensorflow as tf
import numpy as np
tf.Session().run(tf.constant(np.arange(1, 9, dtype=np.int32),shape=[2, 2, 2]))
array([[[1, 2],
        [3, 4]],

       [[5, 6],
        [7, 8]]], dtype=int32)

versus

library(tensorflow)
import("numpy") -> np
tf$constant(np$arange(1, 9, dtype=np$int32),shape=c(2,2,2))
tf$Session()$run(tf$constant(np$arange(1, 9, dtype=np$int32),shape=c(2,2,2)))
, , 1
     [,1] [,2]
[1,]    1    3
[2,]    5    7
, , 2
     [,1] [,2]
[1,]    2    4
[2,]    6    8

again, as @jjallaire shows, this is on the R marshaling side because we see the right order in Python:

py <- import("__builtin__")
py$print(tf$Session()$run(tf$constant(np$arange(1, 9, dtype=np$int32),shape=c(2,2,2))))
[[[1 2]
  [3 4]]

 [[5 6]
  [7 8]]]

from reticulate.

I just exported an r_to_py function which enabled me to write this test code:

> library(reticulate)

> a <- array(c(1:12), dim = c(2,3,2))
> a
, , 1

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

, , 2

     [,1] [,2] [,3]
[1,]    7    9   11
[2,]    8   10   12

 
> r_to_py(a)
[[[ 1  7]
  [ 3  9]
  [ 5 11]]

 [[ 2  8]
  [ 4 10]
  [ 6 12]]]

@bwlewis I'm not thinking clearly in multiple Fortran vs. C dimensions right now! Do these results surprise you? (it seems to me like the arrays are different however it could be that I'm missing a subtlety of the distinct dimension ordering).

from reticulate.

From @bwlewis on February 14, 2017 16:36

@jjallaire It's a memory row/column order problem. I need to carefully look at the numpy source. I suspect this:

• since default numpy is row-ordered, that means it probably links to a "c" blas instead of blas.
• results from that cblas are probably ordered by row, maybe not respecting the NP_ARRAY_FARRAY_RO flag?
• even given NP_ARRAY_FARRAY_RO flag, if numpy is lazy then most buffers are copied and re-ordered anyway. Even if numpy is clever that is not avoidable in many cases. The effort to make zero-copy representations between R and Python is wasted in this case. Needs further investigation...

Anyway, will carefully look at numpy source tonight.

from reticulate.

From @G-Lynn on February 14, 2017 17:0

@jjallaire @bwlewis Thank you for looking into this!

from reticulate.

Okay, @bwlewis I am eternally grateful for your help on this!

from reticulate.

From @eddelbuettel on February 14, 2017 17:13

Department of foggy memory here, but IIRC in package RcppCNPy we always end up transposing (for R) what we get (via the cnpy library) as NumPy data.

from reticulate.

Would one possibility be that tensorflow is taking the NumPy array and accessing it's buffer directly with the assumption of row-major ordering? Under the hood tensorflow is using Eigen and I'm guessing they don't make copies of the NumPy arrays.

from reticulate.

I think there are problems going into Python as well. Here's a simple example where on the R side we have 4 2x3 arrays but on the Python side we end up with 2 3x4 arrays:

> library(reticulate)

> a <- array(c(1:24), dim = c(2,3,4))

> a
, , 1

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

, , 2

     [,1] [,2] [,3]
[1,]    7    9   11
[2,]    8   10   12

, , 3

     [,1] [,2] [,3]
[1,]   13   15   17
[2,]   14   16   18

, , 4

     [,1] [,2] [,3]
[1,]   19   21   23
[2,]   20   22   24


> py_a <- r_to_py(a)

> py_a
[[[ 1  7 13 19]
  [ 3  9 15 21]
  [ 5 11 17 23]]

 [[ 2  8 14 20]
  [ 4 10 16 22]
  [ 6 12 18 24]]]

from reticulate.

From @bwlewis on February 14, 2017 18:22

in these examples "transpose" needs to apply to a 3d array.
On Feb 14, 2017 12:13, "Dirk Eddelbuettel" [email protected] wrote:

Department of foggy memory here, but IIRC in package RcppCNPy we always
end up transposing (for R) what we get (via the cnpy library) as NumPy data.

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
rstudio/tensorflow#88 (comment),
or mute the thread
https://github.com/notifications/unsubscribe-auth/AAIsnisK1rXOuJLi5fELQ5HgrxhKqT1Sks5rceDGgaJpZM4L_3Ni
.

from reticulate.

From @bwlewis on February 14, 2017 18:23

possible! would perhaps then be a python/tensorflow bug.
On Feb 14, 2017 13:11, "JJ Allaire" [email protected] wrote:

Would one possibility be that tensorflow is taking the NumPy array and
accessing it's buffer directly with the assumption of row-major ordering?
Under the hood tensorflow is using Eigen and I'm guessing they don't make
copies of the NumPy arrays.

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
rstudio/tensorflow#88 (comment),
or mute the thread
https://github.com/notifications/unsubscribe-auth/AAIsnrrnBc78Zxul-bvQfTrca5NEZ6iXks5rce5YgaJpZM4L_3Ni
.

from reticulate.

But then again it seems like my simple repro above that has nothing to do with tensorflow indicates we definitely have problems of our own.

from reticulate.

Here's another example which I believe illustrates that the order of the third dimension differs between R and NumPy:

> library(reticulate)
> a <- array(c(1:24), dim = c(2,3,4))
> a
, , 1

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

, , 2

     [,1] [,2] [,3]
[1,]    7    9   11
[2,]    8   10   12

, , 3

     [,1] [,2] [,3]
[1,]   13   15   17
[2,]   14   16   18

, , 4

     [,1] [,2] [,3]
[1,]   19   21   23
[2,]   20   22   24

> py_a <- r_to_py(aperm(a, c(3,1,2)))
> py_a
[[[ 1  3  5]
  [ 2  4  6]]

 [[ 7  9 11]
  [ 8 10 12]]

 [[13 15 17]
  [14 16 18]]

 [[19 21 23]
  [20 22 24]]]
> 

The printed arrays are the same because ee used aperm to move the 3rd dimension to the first slot.

This is sufficiently unfamiliar territory to me that I could be mistaken. Please advise :-)

from reticulate.

This stack overflow thread indicates that NumPy prints 3-d arrays in a way that is somewhat unintuitive to those of us familiar with R and Matlab:

http://stackoverflow.com/questions/22981845/3-dimensional-array-in-numpy

Could this be the source of all the problems/confusion? (i.e. the arrays are correct just printed incomprehensibly).

from reticulate.

This example proves that when 3-d arrays are marshaled back and forth they preserve their equality:

library(reticulate)

a <- array(c(1:24), dim = c(2,3,4))

py_a <- r_to_py(a)

r_a <- py_to_r(py_a)

identical(a, r_a)

from reticulate.

Another example which illustrates that NumPy receives the 3d array correctly:

> library(reticulate)

> a <- array(c(1:24), dim = c(2,3,4))

> py_a <- r_to_py(a)

> py_a$shape
[[1]]
[1] 2

[[2]]
[1] 3

[[3]]
[1] 4


> py_a$flatten(order = 'F')
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

> r_a <- py_to_r(py_a)

> identical(a, r_a)
[1] TRUE

from reticulate.

Here's another yet example showing that we get different NumPy arrays depending on whether we create them in R or NumPy

> library(reticulate)

> np <- import("numpy")

> # create a numpy 3d array from an R array
> a <- array(c(1:24), dim = c(2,3,4))

> r_to_py(a)
[[[ 1  7 13 19]
  [ 3  9 15 21]
  [ 5 11 17 23]]

 [[ 2  8 14 20]
  [ 4 10 16 22]
  [ 6 12 18 24]]]

> # create a numpy 3d array via Python
> py <- py_run_string("import numpy; a = numpy.arange(1,25).reshape(2,3,4)")

> py_get_attr(py, "a")    
[[[ 1  2  3  4]
  [ 5  6  7  8]
  [ 9 10 11 12]]

 [[13 14 15 16]
  [17 18 19 20]
  [21 22 23 24]]]
> 

I'm so confused! I feel like rather than 2 hard problems in computer science (naming things and caching invalidation) we now need to add a 3rd (in-memory array representations)!.

from reticulate.

The reason they are different is that the higher order dimensions go first in Python (e.g. it should be reshape(4, 2, 3) to match the dim = c(2,3,4) used in R). So there may not be a problem here.

from reticulate.

The NumPy C API is also reporting dims with higher order dimensions first! This means that we are getting the dim attribute wrong in R (the reason that we round trip successfully is that we get it wrong going in and going out in the same way). I'll keep experimenting.

from reticulate.

EDIT on above CBLAS comment: I forgot that cblas is just a lightweight wrapper over blas that inserts transposes as needed, see for instance http://www.netlib.org/blas/blast-forum/cblas.tgz -- so my comments on copying are wrong above, no copies are needed.

from reticulate.

@jjallaire Consider your example above:

library(reticulate)
np <- import("numpy")
a <- array(c(1:24), dim = c(2,3,4))
p <- r_to_py(a)
p$flags
#  C_CONTIGUOUS : False
#  F_CONTIGUOUS : True
#  OWNDATA : False
#  WRITEABLE : False
#  ALIGNED : True
#  UPDATEIFCOPY : False

On the other hand:

py_run_string("import numpy; a = numpy.arange(1,25).reshape(2,3,4); print  a.flags")
#  C_CONTIGUOUS : True
#  F_CONTIGUOUS : False
#  OWNDATA : False
#  WRITEABLE : True
#  ALIGNED : True
#  UPDATEIFCOPY : False

Thus:

1. we need to check numpy array attributes on return to R to see if they are column or row ordered and adjust accordingly on the return.
2. This may cause some inefficiencies in binary operators with mixed order arguments...in that one case, one of those arrays needs to be copied internally by Python. I don't think this is a super-serious problem.

This leads to the question, why isn't this line of C++ code translating the array for us? https://github.com/rstudio/reticulate/blob/master/src/python.cpp#L549:

array = (PyArrayObject*)PyArray_CastToType(array, descr, NPY_ARRAY_FARRAY);

Getting back to the original problem, there may be cases in tensorflow where row order is assumed. But one thing at a time.

from reticulate.

Maybe a Python n-d array bug in PyArray_CastToType? It works fine in 2-d. Here using the old bwlewis/python package to help see what's going on:

x <- matrix(1:25.0, 5)
x %*% x
#     [,1] [,2] [,3] [,4] [,5]
#[1,]  215  490  765 1040 1315
#[2,]  230  530  830 1130 1430
#[3,]  245  570  895 1220 1545
#[4,]  260  610  960 1310 1660
#[5,]  275  650 1025 1400 1775

library(python)
import("numpy") -> np
p <- np$array(x)
p$dot(p)$flags
#  C_CONTIGUOUS : True          ## <<<---- NOTE
#  F_CONTIGUOUS : False
#  OWNDATA : True
#  WRITEABLE : True
#  ALIGNED : True
#  UPDATEIFCOPY : False

# But, it comes back to R in this example correctly!
R(p$dot(p))
#     [,1] [,2] [,3] [,4] [,5]
#[1,]  215  490  765 1040 1315
#[2,]  230  530  830 1130 1430
#[3,]  245  570  895 1220 1545
#[4,]  260  610  960 1310 1660
#[5,]  275  650 1025 1400 1775

I think the problems are limited to the n-d case, n > 2.

from reticulate.

Also, just a nit on the use of flatten above, from the Python API doc:

order : {‘C’, ‘F’, ‘A’, ‘K’}, optional
‘C’ means to flatten in row-major (C-style) order.
‘F’ means to flatten in column-major (Fortran- style) order.
‘A’ means to flatten in column-major order if a is Fortran contiguous in memory, row-major order otherwise.
‘K’ means to flatten a in the order the elements occur in memory.
The default is ‘C’.

Really, shouldn't the default be K or, arguably, A? Seems nuts to me otherwise.

from reticulate.

And just to verify that your example above is not merely a print/display problem, these things are in fact different within Python. So there is a problem in the assigment of the Python objects from R?

library(reticulate)
np <- import("numpy")
a <- array(c(1:24), dim = c(2,3,4))
pf <- r_to_py(a)
pc <- py_run_string("import numpy; a = numpy.arange(1,25).reshape(2,3,4)")
pc <- py_get_attr(pc, "a")

o <- import("operator")
o$sub(pf, pc)
#, , 1
#
#     [,1] [,2] [,3]
#[1,]    0   -2   -4
#[2,]  -11  -13  -15

# ...  (should all be zero)

from reticulate.

Maybe it's not us. Forget about R for a minute:

library(reticulate)
np <- import("numpy")
pf <- py_run_string("import numpy; a = numpy.reshape(numpy.arange(1,25), (2,3,4), 'F')")
pf <- py_get_attr(pf, "a")
pc <- py_run_string("import numpy; b = numpy.arange(1,25).reshape(2,3,4)")
pc <- py_get_attr(pc, "b")

o <- import("operator")
o$sub(pf, pc)
#, , 1
#
#     [,1] [,2] [,3]
#[1,]    0   -2   -4
#[2,]  -11  -13  -15   # ...  (not all zero)

EDIT: I'm wrong, and have been mis-interpreting the order flag above. It's not simply saying how to fill data structures by rows first ("C_CONTIGUOUS") or columns first ("F_CONTIGUOUS"), but also a signal to all the operators on the element order for arithmetic. In the above example, the data layout are identical in memory -- see this with pc$flatten('K') and pf$flatten('K'). Despite this, the subtract operator sees the flags and adjusts element order to match.

from reticulate.

OK, getting back to the original problem, here it is in a nutshell:

library(reticulate)
np <- import("numpy")
# 2-d
a <- py_run_string("import numpy; a = numpy.arange(1, 5).reshape(2, 2)")
a <- py_get_attr(a, "a")

# 3-d not consistent
b <- py_run_string("import numpy; b = numpy.arange(1, 9).reshape(2, 2, 2)")
b <- py_get_attr(b, "b")

a
#[[1 2]
# [3 4]]
py_to_r(a)
#     [,1] [,2]
#[1,]    1    2
#[2,]    3    4

b
#[[[1 2]
#  [3 4]]
#
# [[5 6]
#  [7 8]]]
py_to_r(b)
#, , 1
#     [,1] [,2]
#[1,]    1    3
#[2,]    5    7
#
#, , 2
#     [,1] [,2]
#[1,]    2    4
#[2,]    6    8

There is an inconsistency in how C_CONTIGUOUS arrays are brought into R from Python. Things work as expected in the 2-d case, but not so in the 3-d case. This difference is unexpected.

from reticulate.

from reticulate.

from reticulate.

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.
The R array you construct in the example is not the same as the
one constructed in the reference Python code, but it's really easy
to confuse them!

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

An extended version of this discussion will appear in a package vignette.

from reticulate.