Source code for statsmodels.tools.numdiff

"""numerical differentiation function, gradient, Jacobian, and Hessian

Author : josef-pkt
License : BSD

Notes
-----
These are simple forward differentiation, so that we have them available
without dependencies.

* Jacobian should be faster than numdifftools because it doesn't use loop over
  observations.
* numerical precision will vary and depend on the choice of stepsizes
"""

# TODO:
# * some cleanup
# * check numerical accuracy (and bugs) with numdifftools and analytical
#   derivatives
#   - linear least squares case: (hess - 2*X'X) is 1e-8 or so
#   - gradient and Hessian agree with numdifftools when evaluated away from
#     minimum
#   - forward gradient, Jacobian evaluated at minimum is inaccurate, centered
#     (+/- epsilon) is ok
# * dot product of Jacobian is different from Hessian, either wrong example or
#   a bug (unlikely), or a real difference
#
#
# What are the conditions that Jacobian dotproduct and Hessian are the same?
#
# See also:
#
# BHHH: Greene p481 17.4.6,  MLE Jacobian = d loglike / d beta , where loglike
# is vector for each observation
#    see also example 17.4 when J'J is very different from Hessian
#    also does it hold only at the minimum, what's relationship to covariance
#    of Jacobian matrix
# http://projects.scipy.org/scipy/ticket/1157
# http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm
#    objective: sum((y-f(beta,x)**2),   Jacobian = d f/d beta
#    and not d objective/d beta as in MLE Greene
#    similar: http://crsouza.blogspot.com/2009/11/neural-network-learning-by-levenberg_18.html#hessian
#
# in example: if J = d x*beta / d beta then J'J == X'X
#    similar to http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm
from __future__ import print_function
from statsmodels.compat.python import range
import numpy as np

# NOTE: we only do double precision internally so far
EPS = np.MachAr().eps

_hessian_docs = """
    Calculate Hessian with finite difference derivative approximation

    Parameters
    ----------
    x : array_like
       value at which function derivative is evaluated
    f : function
       function of one array f(x, `*args`, `**kwargs`)
    epsilon : float or array-like, optional
       Stepsize used, if None, then stepsize is automatically chosen
       according to EPS**(1/%(scale)s)*x.
    args : tuple
        Arguments for function `f`.
    kwargs : dict
        Keyword arguments for function `f`.
    %(extra_params)s

    Returns
    -------
    hess : ndarray
       array of partial second derivatives, Hessian
    %(extra_returns)s

    Notes
    -----
    Equation (%(equation_number)s) in Ridout. Computes the Hessian as::

      %(equation)s

    where e[j] is a vector with element j == 1 and the rest are zero and
    d[i] is epsilon[i].

    References
    ----------:

    Ridout, M.S. (2009) Statistical applications of the complex-step method
        of numerical differentiation. The American Statistician, 63, 66-74
"""


def _get_epsilon(x, s, epsilon, n):
    if epsilon is None:
        h = EPS**(1. / s) * np.maximum(np.abs(x), 0.1)
    else:
        if np.isscalar(epsilon):
            h = np.empty(n)
            h.fill(epsilon)
        else:  # pragma : no cover
            h = np.asarray(epsilon)
            if h.shape != x.shape:
                raise ValueError("If h is not a scalar it must have the same"
                                 " shape as x.")
    return h


[docs]def approx_fprime(x, f, epsilon=None, args=(), kwargs={}, centered=False): ''' Gradient of function, or Jacobian if function f returns 1d array Parameters ---------- x : array parameters at which the derivative is evaluated f : function `f(*((x,)+args), **kwargs)` returning either one value or 1d array epsilon : float, optional Stepsize, if None, optimal stepsize is used. This is EPS**(1/2)*x for `centered` == False and EPS**(1/3)*x for `centered` == True. args : tuple Tuple of additional arguments for function `f`. kwargs : dict Dictionary of additional keyword arguments for function `f`. centered : bool Whether central difference should be returned. If not, does forward differencing. Returns ------- grad : array gradient or Jacobian Notes ----- If f returns a 1d array, it returns a Jacobian. If a 2d array is returned by f (e.g., with a value for each observation), it returns a 3d array with the Jacobian of each observation with shape xk x nobs x xk. I.e., the Jacobian of the first observation would be [:, 0, :] ''' n = len(x) # TODO: add scaled stepsize f0 = f(*((x,)+args), **kwargs) dim = np.atleast_1d(f0).shape # it could be a scalar grad = np.zeros((n,) + dim, np.promote_types(float, x.dtype)) ei = np.zeros((n,), float) if not centered: epsilon = _get_epsilon(x, 2, epsilon, n) for k in range(n): ei[k] = epsilon[k] grad[k, :] = (f(*((x+ei,) + args), **kwargs) - f0)/epsilon[k] ei[k] = 0.0 else: epsilon = _get_epsilon(x, 3, epsilon, n) / 2. for k in range(len(x)): ei[k] = epsilon[k] grad[k, :] = (f(*((x+ei,)+args), **kwargs) - f(*((x-ei,)+args), **kwargs))/(2 * epsilon[k]) ei[k] = 0.0 return grad.squeeze().T
[docs]def approx_fprime_cs(x, f, epsilon=None, args=(), kwargs={}): ''' Calculate gradient or Jacobian with complex step derivative approximation Parameters ---------- x : array parameters at which the derivative is evaluated f : function `f(*((x,)+args), **kwargs)` returning either one value or 1d array epsilon : float, optional Stepsize, if None, optimal stepsize is used. Optimal step-size is EPS*x. See note. args : tuple Tuple of additional arguments for function `f`. kwargs : dict Dictionary of additional keyword arguments for function `f`. Returns ------- partials : ndarray array of partial derivatives, Gradient or Jacobian Notes ----- The complex-step derivative has truncation error O(epsilon**2), so truncation error can be eliminated by choosing epsilon to be very small. The complex-step derivative avoids the problem of round-off error with small epsilon because there is no subtraction. ''' # From Guilherme P. de Freitas, numpy mailing list # May 04 2010 thread "Improvement of performance" # http://mail.scipy.org/pipermail/numpy-discussion/2010-May/050250.html n = len(x) epsilon = _get_epsilon(x, 1, epsilon, n) increments = np.identity(n) * 1j * epsilon # TODO: see if this can be vectorized, but usually dim is small partials = [f(x+ih, *args, **kwargs).imag / epsilon[i] for i, ih in enumerate(increments)] return np.array(partials).T
[docs]def approx_hess_cs(x, f, epsilon=None, args=(), kwargs={}): '''Calculate Hessian with complex-step derivative approximation Parameters ---------- x : array_like value at which function derivative is evaluated f : function function of one array f(x) epsilon : float stepsize, if None, then stepsize is automatically chosen Returns ------- hess : ndarray array of partial second derivatives, Hessian Notes ----- based on equation 10 in M. S. RIDOUT: Statistical Applications of the Complex-step Method of Numerical Differentiation, University of Kent, Canterbury, Kent, U.K. The stepsize is the same for the complex and the finite difference part. ''' # TODO: might want to consider lowering the step for pure derivatives n = len(x) h = _get_epsilon(x, 3, epsilon, n) ee = np.diag(h) hess = np.outer(h, h) n = len(x) for i in range(n): for j in range(i, n): hess[i, j] = (f(*((x + 1j*ee[i, :] + ee[j, :],) + args), **kwargs) - f(*((x + 1j*ee[i, :] - ee[j, :],)+args), **kwargs)).imag/2./hess[i, j] hess[j, i] = hess[i, j] return hess
approx_hess_cs.__doc__ = (("Calculate Hessian with complex-step derivative " "approximation\n") + "\n".join(_hessian_docs.split("\n")[1:]) % dict(scale="3", extra_params="", extra_returns="", equation_number="10", equation=("1/(2*d_j*d_k) * " "imag(f(x + i*d[j]*e[j] + " "d[k]*e[k]) -\n" " " "f(x + i*d[j]*e[j] - d[k]*e[k]))\n")) )
[docs]def approx_hess1(x, f, epsilon=None, args=(), kwargs={}, return_grad=False): n = len(x) h = _get_epsilon(x, 3, epsilon, n) ee = np.diag(h) f0 = f(*((x,)+args), **kwargs) # Compute forward step g = np.zeros(n) for i in range(n): g[i] = f(*((x+ee[i, :],)+args), **kwargs) hess = np.outer(h, h) # this is now epsilon**2 # Compute "double" forward step for i in range(n): for j in range(i, n): hess[i, j] = (f(*((x + ee[i, :] + ee[j, :],) + args), **kwargs) - g[i] - g[j] + f0)/hess[i, j] hess[j, i] = hess[i, j] if return_grad: grad = (g - f0)/h return hess, grad else: return hess
approx_hess1.__doc__ = _hessian_docs % dict(scale="3", extra_params="""return_grad : bool Whether or not to also return the gradient """, extra_returns="""grad : nparray Gradient if return_grad == True """, equation_number="7", equation="""1/(d_j*d_k) * ((f(x + d[j]*e[j] + d[k]*e[k]) - f(x + d[j]*e[j]))) """)
[docs]def approx_hess2(x, f, epsilon=None, args=(), kwargs={}, return_grad=False): # n = len(x) # NOTE: ridout suggesting using eps**(1/4)*theta h = _get_epsilon(x, 3, epsilon, n) ee = np.diag(h) f0 = f(*((x,)+args), **kwargs) # Compute forward step g = np.zeros(n) gg = np.zeros(n) for i in range(n): g[i] = f(*((x+ee[i, :],)+args), **kwargs) gg[i] = f(*((x-ee[i, :],)+args), **kwargs) hess = np.outer(h, h) # this is now epsilon**2 # Compute "double" forward step for i in range(n): for j in range(i, n): hess[i, j] = (f(*((x + ee[i, :] + ee[j, :],) + args), **kwargs) - g[i] - g[j] + f0 + f(*((x - ee[i, :] - ee[j, :],) + args), **kwargs) - gg[i] - gg[j] + f0)/(2 * hess[i, j]) hess[j, i] = hess[i, j] if return_grad: grad = (g - f0)/h return hess, grad else: return hess
approx_hess2.__doc__ = _hessian_docs % dict(scale="3", extra_params="""return_grad : bool Whether or not to also return the gradient """, extra_returns="""grad : nparray Gradient if return_grad == True """, equation_number="8", equation = """1/(2*d_j*d_k) * ((f(x + d[j]*e[j] + d[k]*e[k]) - f(x + d[j]*e[j])) - (f(x + d[k]*e[k]) - f(x)) + (f(x - d[j]*e[j] - d[k]*e[k]) - f(x + d[j]*e[j])) - (f(x - d[k]*e[k]) - f(x))) """)
[docs]def approx_hess3(x, f, epsilon=None, args=(), kwargs={}): n = len(x) h = _get_epsilon(x, 4, epsilon, n) ee = np.diag(h) hess = np.outer(h,h) for i in range(n): for j in range(i, n): hess[i, j] = (f(*((x + ee[i, :] + ee[j, :],) + args), **kwargs) - f(*((x + ee[i, :] - ee[j, :],) + args), **kwargs) - (f(*((x - ee[i, :] + ee[j, :],) + args), **kwargs) - f(*((x - ee[i, :] - ee[j, :],) + args), **kwargs)) )/(4.*hess[i, j]) hess[j, i] = hess[i, j] return hess
approx_hess3.__doc__ = _hessian_docs % dict(scale="4", extra_params="", extra_returns="", equation_number="9", equation = """1/(4*d_j*d_k) * ((f(x + d[j]*e[j] + d[k]*e[k]) - f(x + d[j]*e[j] - d[k]*e[k])) - (f(x - d[j]*e[j] + d[k]*e[k]) - f(x - d[j]*e[j] - d[k]*e[k]))""") approx_hess = approx_hess3 approx_hess.__doc__ += "\n This is an alias for approx_hess3" if __name__ == '__main__': #pragma : no cover import statsmodels.api as sm from scipy.optimize.optimize import approx_fhess_p import numpy as np data = sm.datasets.spector.load() data.exog = sm.add_constant(data.exog, prepend=False) mod = sm.Probit(data.endog, data.exog) res = mod.fit(method="newton") test_params = [1,0.25,1.4,-7] llf = mod.loglike score = mod.score hess = mod.hessian # below is Josef's scratch work def approx_hess_cs_old(x, func, args=(), h=1.0e-20, epsilon=1e-6): def grad(x): return approx_fprime_cs(x, func, args=args, h=1.0e-20) #Hessian from gradient: return (approx_fprime(x, grad, epsilon) + approx_fprime(x, grad, -epsilon))/2. def fun(beta, x): return np.dot(x, beta).sum(0) def fun1(beta, y, x): #print(beta.shape, x.shape) xb = np.dot(x, beta) return (y-xb)**2 #(xb-xb.mean(0))**2 def fun2(beta, y, x): #print(beta.shape, x.shape) return fun1(beta, y, x).sum(0) nobs = 200 x = np.arange(nobs*3).reshape(nobs,-1) x = np.random.randn(nobs,3) xk = np.array([1,2,3]) xk = np.array([1.,1.,1.]) #xk = np.zeros(3) beta = xk y = np.dot(x, beta) + 0.1*np.random.randn(nobs) xk = np.dot(np.linalg.pinv(x),y) epsilon = 1e-6 args = (y,x) from scipy import optimize xfmin = optimize.fmin(fun2, (0,0,0), args) print(approx_fprime((1,2,3),fun,epsilon,x)) jac = approx_fprime(xk,fun1,epsilon,args) jacmin = approx_fprime(xk,fun1,-epsilon,args) #print(jac) print(jac.sum(0)) print('\nnp.dot(jac.T, jac)') print(np.dot(jac.T, jac)) print('\n2*np.dot(x.T, x)') print(2*np.dot(x.T, x)) jac2 = (jac+jacmin)/2. print(np.dot(jac2.T, jac2)) #he = approx_hess(xk,fun2,epsilon,*args) print(approx_hess_old(xk,fun2,1e-3,args)) he = approx_hess_old(xk,fun2,None,args) print('hessfd') print(he) print('epsilon =', None) print(he[0] - 2*np.dot(x.T, x)) for eps in [1e-3,1e-4,1e-5,1e-6]: print('eps =', eps) print(approx_hess_old(xk,fun2,eps,args)[0] - 2*np.dot(x.T, x)) hcs2 = approx_hess_cs(xk,fun2,args=args) print('hcs2') print(hcs2 - 2*np.dot(x.T, x)) hfd3 = approx_hess(xk,fun2,args=args) print('hfd3') print(hfd3 - 2*np.dot(x.T, x)) import numdifftools as nd hnd = nd.Hessian(lambda a: fun2(a, y, x)) hessnd = hnd(xk) print('numdiff') print(hessnd - 2*np.dot(x.T, x)) #assert_almost_equal(hessnd, he[0]) gnd = nd.Gradient(lambda a: fun2(a, y, x)) gradnd = gnd(xk)