Calculation of gradients

I write codes for calculating gradients of real valued function.
But the result is not desired one. For example, gradient of f(x,y) = x^2 + y^2 at (1,1) should be (2,2) but in my codes, it is (1,1). Where do I misunderstand?

Here, gradient is a vector whose i-th component is given by

{ (f(x1,x2,…,xi+h,…) – f(x1,x2,…,xi,…) ) /h

+(f(x1,x2,…,xi-h,…) – f(x1,x2,…,xi,…) ) /(-h) }/2.

from numpy import *

def grad(f,x):
    h=1e-4
    print("h =",h)
    grad = zeros_like(x)
    print( " grad = zeros_like(x):"  ,grad )
    for i in range(x.size):
        print("i=", i)
        grad[i] = (f(x[i] + h) - f(x[i] - h)) / (2*h)
        print(" grad[",i,"] =", grad[i] )
    return grad






def f(x):
    return sum(x**2)

print("should be (2,2) but it is ",  grad(f,array([1,1])) )

>Solution :

Use

    grad = zeros_like(x, dtype=float)

All the x values you use are integer type. The division operation will convert the result to float, but since grad is set to be like x, it will have an integer type, and the floating point results will be truncated to integer before assigned to grad[i]. Hence you get a truncated floating point value as integer, which produces the result of what you see.

You may expect the result of the gradient calculation to be 2 when x[i] == 1. But the floating point division is imprecise (in decimal), thus it may be something like 1.999998, which truncates to 1 before it is assigned. That’s why your result is 1, not 2.

The worst part is that if you have set h to something like h=0.015625, which is 2**-6 and can be represented exactly in binary, the result would have been (2, 2), since the division results in precisely 2.0 in binary, and thus the truncation to integer values would result in 2. This can lead to (hidden) bugs which will only appear weeks or months later, if at all (but may lead to incorrect results in published results, with all kind of possible consequences. Just as an example of a worst-case scenario.)