def compute_fd_hessian(fun, x0, epsilon, anagrad=True):
N = len(x0)<EOL>def compute_gradient(x):<EOL><INDENT>if anagrad:<EOL><INDENT>return fun(x, do_gradient=True)[<NUM_LIT:1>]<EOL><DEDENT>else:<EOL><INDENT>gradient = np.zeros(N, float)<EOL>for i in range(N):<EOL><INDENT>xh = x.copy()<EOL>xh[i] += <NUM_LIT:0.5>*epsilon<EOL>xl = x.copy()<EOL>xl[i] -= <NUM_LIT:0.5>*epsilon<EOL>gradient[i] = (fun(xh)-fun(xl))/epsilon<EOL><DEDENT>return gradient<EOL><DEDENT><DEDENT>hessian = np.zeros((N,N), float)<EOL>for i in range(N):<EOL><INDENT>xh = x0.copy()<EOL>xh[i] += <NUM_LIT:0.5>*epsilon<EOL>xl = x0.copy()<EOL>xl[i] -= <NUM_LIT:0.5>*epsilon<EOL>hessian[i] = (compute_gradient(xh) - compute_gradient(xl))/epsilon<EOL><DEDENT>return <NUM_LIT:0.5>*(hessian + hessian.transpose())<EOL>
Compute the Hessian using the finite difference method

       Arguments:
        | ``fun``  --  the function for which the Hessian should be computed,
                       more info below
        | ``x0``  --  the point at which the Hessian must be computed
        | ``epsilon``  --  a small scalar step size used to compute the finite
                           differences

       Optional argument:
        | ``anagrad``  --  when True, analytical gradients are used
                           [default=True]

       The function ``fun`` takes a mandatory argument ``x`` and an optional
       argument ``do_gradient``:
        | ``x``  --  the arguments of the function to be tested
        | ``do_gradient``  --  When False, only the function value is returned.
                               When True, a 2-tuple with the function value and
                               the gradient are returned [default=False]