用numdifftools计算hessian矩阵

Question

我想学习如何使用numdifftools制作函数的Hessian矩阵

例如，对于平方误差梯度函数：

import numpy as np

def predict(x_i, beta):
    """assumes that the first element of each x_i is 1"""
    return np.dot(x_i, beta)

def error(x_i, y_i, beta):
    return y_i - predict(x_i, beta)

def squared_error(x_i, y_i, beta):
    return error(x_i, y_i, beta) ** 2

def squared_error_gradient(x_i, y_i, beta):
    """The gradient (with respect to beta)
    correspond to the ith squared error term"""
    return [-2 * x_ij * error(x_i, y_i, beta)
           for x_ij in x_i]

我想从这个函数中得到Hessian矩阵。实际上，代价函数的最佳学习率严格小于2λ，其中λ是Hessian函数的最大特征值。所以我试着：

import numdifftools as nd

Hfun = nd.Hessian(squared_error_gradient)
h = Hfun([1, 49, 4, 2500], # X
         [2000], # y
         [random.random() for x in [1, 49, 4, 2500]]) # beta

但它的回报是：

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-22-576c9ee22a18> in <module>
      1 import numdifftools as nd
      2 Hfun = nd.Hessian(squared_error_gradient)
----> 3 h = Hfun([1, 49, 4, 2500],[2000], [random.random() for x in [1, 49, 4, 2500]])
      4 h

/opt/conda/lib/python3.7/site-packages/numdifftools/core.py in __call__(self, x, *args, **kwds)
    847 
    848     def __call__(self, x, *args, **kwds):
--> 849         return super(Hessdiag, self).__call__(np.atleast_1d(x), *args, **kwds)
    850 
    851 

/opt/conda/lib/python3.7/site-packages/numdifftools/core.py in __call__(self, x, *args, **kwds)
    374     def __call__(self, x, *args, **kwds):
    375         xi = np.asarray(x)
--> 376         results = self._derivative(xi, args, kwds)
    377         derivative, info = self._extrapolate(*results)
    378         if self.full_output:

/opt/conda/lib/python3.7/site-packages/numdifftools/core.py in _derivative_nonzero_order(self, xi, args, kwds)
    258         steps, step_ratio = self._get_steps(xi)
    259         fxi = self._eval_first(f, xi)
--> 260         results = [diff(f, fxi, xi, h) for h in steps]
    261 
    262         self.set_richardson_rule(step_ratio, self.richardson_terms)

/opt/conda/lib/python3.7/site-packages/numdifftools/core.py in <listcomp>(.0)
    258         steps, step_ratio = self._get_steps(xi)
    259         fxi = self._eval_first(f, xi)
--> 260         results = [diff(f, fxi, xi, h) for h in steps]
    261 
    262         self.set_richardson_rule(step_ratio, self.richardson_terms)

/opt/conda/lib/python3.7/site-packages/numdifftools/core.py in _central_even(f, fx, x, h)
    888         n = len(x)
    889         ee = np.diag(h)
--> 890         dtype = np.result_type(fx, float)  # make sure it is at least float64
    891         hess = np.empty((n, n), dtype=dtype)
    892         np.outer(h, h, out=hess)

<__array_function__ internals> in result_type(*args, **kwargs)

TypeError: Field elements must be 2- or 3-tuples, got 'array([-1411.17618406])'