# how to compute the gradient of a function at an extremal point

I am writing a computer program that searches for the minimum of a multivariate function $f: \mathbb{R}^n \to \mathbb{R}$. This function is in fact the sum of many functions:

$$f(x) = \sum_{i=1}^m f_i(x)$$

The computation starts from some guess $x_0$ of a minimum point, and at each iteration $j$ it improves the current guess $x_j$ into a better one $x_{j+1}$ (using some variant of Newton's method or steepest descent). As $x_j$ approaches the minimum, the gradient becomes closer to 0. The problem is, that to compute the gradient, I need to sum up many gradients, like this

$$\nabla f = \sum_{i=1}^m \nabla f_i$$

and the $\nabla f_i$'s are nonzero. The fact that their sum is close to zero means that if we compute this sum in floating point arithmetics, we get a big accumulated error and the computation is very inaccurate. When I try to approximate the gradient by evaluating $f$ at different points near $x_j$, the approximate gradient comes out indeed quite different from the computed gradient. (In comparison, at points where $\nabla f$ is far from zero, the computed and approximate gradients are very close).

Does anyone know of a good method to compute the gradient near extremal points?

-
Did you try using conjugate gradient methods? –  yohBS Sep 23 '12 at 12:38
i just looked it up in wikipedia and this method seems to use the gradient. but in my case the gradient is inaccurate, so how can this method help? –  roel Sep 23 '12 at 13:17
If the the minimum you're looking for is quadratic, that is, if the hessian matrix at the minimum is non-singular, I think it should work. Look up Numerical Recipies Chapter 10.8 (I have the 3rd edition), it is very well explained there. –  yohBS Sep 23 '12 at 20:56
There are techniques for accurately computing the sum of many floating-point numbers without too much cancellation error; see Higham's "Accurate and Efficient Floating Point Summation". (I presume you're using double precision already, of course -- if not, forget everything and do that first.) –  Rahul Jun 7 '13 at 0:12

There are basically three mathods to provide the gradients to an optimization algorithm:

1. You provide a function gradF(x1,x2...) in which you implement the analytic form of the gradient function.

2. You use finite-differences to numerically approximate the gradient.

3. You can use symbolic (automatic) differentiation (i.e. to generate the gradient function implementation given the implementation of your objective function).

Also, you have to understand the behaviour of your function and you have to adjust the termination conditions such as function tolerance (i.e. $f(x_i)-f(x_{i-1})$, where $i$ is the step number) and the value tolerance (i.e. $x_i-x_{i-1}$).

-
I'm using method #1 - gradF(...) is implemented analytically. To test gradF for correctness, I use method #2 to get the "observed" gradient, and compare the results. As for the error tolerance, the fact that gradF becomes increasingly inaccurate in the neighborhood of the minimum gives me very bad error tolerance. –  roel Sep 24 '12 at 8:04
It could be possible that the resulting function from adding your $n$ functions will have a flat region in which the optimization algorithm can't do any progress. What are the tolerances that you are using? –  Haider Sep 24 '12 at 8:36
I stop when the gradient norm becomes less than 1e-5. The approximated |gradient| at the stop point is ~ 0.5. The function value at that point is roughly -400000. I estimate that the minimum can be improved by order of 0.1. –  roel Sep 24 '12 at 9:53
use value tolarance $x_i - x_{i-1}$ where $i$ is the step number, as a termination criteria –  Haider Sep 24 '12 at 10:09
If you have a computer program that already computes the values of $f_i$, and you need to compute derivatives, I would strongly avoid finite difference methods, if possible.