Search code examples
pythonarraysnumpyrounding-error

Avoid incorrect rounding when part of a Numpy Array is replaced


In one part of my code I have used the following code with numpy arrays:

import numpy as np

A=np.array([[1,-8],[2,-1],[2,14]])

k=2
x = np.array([[0],[15]])       

v=np.array([[1/np.sqrt(2)],[1/np.sqrt(2)]])

temp = x - 2*v @ v.T @ x

print("temp = ",temp)
I=np.array([[0,0],[0,0],[0,0]])
print("I = ",I)

I[k-1:3,k-1:2]= temp
print("new I =",I)

This is the output that I see:

temp =  
[[-1.50000000e+01]
 [ 3.55271368e-15]]
I =  
[[0 0]
 [0 0]
 [0 0]]
new I = 
[[  0   0]
 [  0 -14]
 [  0   0]]

However -1.50000000e+01 = -15 so I am not sure why this element is being replaced with -14.


Solution

  • If I modify your code to create float dtype arrays, especially I:

    In [88]: A=np.array([[1,-8],[2,-1],[2,14]], float)    ###
        ...:  
        ...: k=2 
        ...: x = np.array([[0],[15]])        
        ...:  
        ...: v=np.array([[0.70710678],[0.70710678]]) 
        ...:  
        ...: temp = x - 2*v @ v.T @ x 
        ...:  
        ...: print("temp = ",temp) 
        ...: I=np.array([[0,0],[0,0],[0,0]], float)      ####
        ...: print("I = ",I) 
        ...:  
        ...: I[k-1:3,k-1:2]= temp 
        ...: print("new I =",I)                                                                    
    temp =  [[-1.49999999e+01]
     [ 5.03409456e-08]]
    I =  [[0. 0.]
     [0. 0.]
     [0. 0.]]
    new I = [[ 0.00000000e+00  0.00000000e+00]
     [ 0.00000000e+00 -1.49999999e+01]
     [ 0.00000000e+00  5.03409456e-08]]
    

    round gives me more control than truncation:

    In [90]: np.round(I)                                                                           
    Out[90]: 
    array([[  0.,   0.],
           [  0., -15.],
           [  0.,   0.]])
    In [91]: I.astype(int)                                                                         
    Out[91]: 
    array([[  0,   0],
           [  0, -14],
           [  0,   0]])