CODEWITHSUNDEEP
pythonarraysnumpysqrt
Vahid Azimi
Vahid Azimi

Reputation: 27

Weird behavior when squaring elements in numpy array

I have two numpy arrays of shape (1, 250000):

a = [[  0 254   1 ..., 255   0   1]]
b = [[  1   0 252 ...,   0 255 255]]

I want to create a new numpy array whose elements are the square root of the sum of squares of elements in the arrays a and b, but I am not getting the correct result:

>>> c = np.sqrt(np.square(a)+np.square(b))
>>> print c
[[ 1.          2.          4.12310553 ...,  1.          1.          1.41421354]]

Am I missing something simple here?

Upvotes: 0

Views: 1843

Answers (2)

Warren Weckesser
Warren Weckesser

Reputation: 114831

Presumably your arrays a and b are arrays of unsigned 8 bit integers--you can check by inspecting the attribute a.dtype. When you square them, the data type is preserved, and the 8 bit values overflow, which means the values "wrap around" (i.e. the squared values are modulo 256):

In [7]: a = np.array([[0, 254, 1, 255, 0, 1]], dtype=np.uint8)

In [8]: np.square(a)
Out[8]: array([[0, 4, 1, 1, 0, 1]], dtype=uint8)

In [9]: b = np.array([[1, 0, 252, 0, 255, 255]], dtype=np.uint8)

In [10]: np.square(a) + np.square(b)
Out[10]: array([[ 1,  4, 17,  1,  1,  2]], dtype=uint8)

In [11]: np.sqrt(np.square(a) + np.square(b))
Out[11]: 
array([[ 1.        ,  2.        ,  4.12310553,  1.        ,  1.        ,
         1.41421354]], dtype=float32)

To avoid the problem, you can tell np.square to use a floating point data type:

In [15]: np.sqrt(np.square(a, dtype=np.float64) + np.square(b, dtype=np.float64))
Out[15]: 
array([[   1.        ,  254.        ,  252.00198412,  255.        ,
         255.        ,  255.00196078]])

You could also use the function numpy.hypot, but you might still want to use the dtype argument, otherwise the default data type is np.float16:

In [16]: np.hypot(a, b)
Out[16]: array([[   1.,  254.,  252.,  255.,  255.,  255.]], dtype=float16)

In [17]: np.hypot(a, b, dtype=np.float64)
Out[17]: 
array([[   1.        ,  254.        ,  252.00198412,  255.        ,
         255.        ,  255.00196078]])

You might wonder why the dtype argument that I used in numpy.square and numpy.hypot is not shown in the functions' docstrings. Both of these functions are numpy "ufuncs", and the authors of numpy decided that it was better to show only the main arguments in the docstring. The optional arguments are documented in the reference manual.

Upvotes: 7

kmario23
kmario23

Reputation: 61355

For this simple case, it works perfectly fine:

In [1]: a = np.array([[ 0, 2, 4, 6, 8]])
In [2]: b = np.array([[ 1, 3, 5, 7, 9]])
In [3]: c = np.sqrt(np.square(a) + np.square(b))

In [4]: print(c)
[[  1.     3.60555128   6.40312424   9.21954446  12.04159458]]

You must be doing something wrong.

Upvotes: 0

Related Questions