numpy/scipy中非线性函数的数值梯度

Question

我试图在numpy中实现一个数值梯度计算，作为cyipot中梯度的回调函数。我对numpy梯度函数的理解是，它应该返回基于finite different approximation的点计算的梯度。在

我不明白如何用这个模块来实现非线性函数的梯度。给出的样本问题似乎是一个线性函数。在

>>> f = np.array([1, 2, 4, 7, 11, 16], dtype=np.float)
>>> np.gradient(f)
array([ 1. ,  1.5,  2.5,  3.5,  4.5,  5. ])
>>> np.gradient(f, 2)
array([ 0.5 ,  0.75,  1.25,  1.75,  2.25,  2.5 ])

我的代码片段如下：

^{pr2}$

另一个缺点是，我必须在几个点上计算x（它返回几个点的梯度） 在numpy/scipy中，是否有更好的选项可以在一个单个点处对梯度进行数值计算，这样我就可以将其作为回调函数来实现了？在

Answer 1

我想你对什么是数学函数和它的数值实现有一些误解。在

您应该将您的功能定义为：

def func(x1, x2, x3, x4):
    return -x1*x2*x3*x4

现在您需要在特定的点计算函数，可以使用您提供的np.mgrid来实现。在

如果要计算渐变，请使用copy.misc.derivative（https://docs.scipy.org/doc/scipy/reference/generated/scipy.misc.derivative.html）（注意dx的默认参数通常不好，请将其更改为1e-5。在数值计算中，线性梯度和非线性梯度没有区别，只是对于非线性函数，梯度并非处处相同。在

使用np.gradient实际上是从数组中的点计算梯度，函数的定义被f的定义所隐藏，因此不允许在不同的点进行多个梯度计算。同样，使用你的方法会使你依赖于你的离散化步骤。在

Answer 2

首先，一些警告：

• 数值优化很难做对
• ipopt是非常复杂的软件
  • 将ipopt与数值微分相结合听起来像是在自找麻烦，但这当然取决于你的问题
  • ipopt几乎总是基于automatic-differentiation tools，而不是numerical-differentiation！在

还有更多：

• 由于这是一项复杂的任务，而且python+ipopt的状态没有其他语言中的好（例如julia+JuMP），所以需要做一些工作

还有一些替代品：

• 使用pyomo，它包装了ipopt并具有自动区分功能
• 使用casadi，它也包装ipopt并具有自动区分功能
• 使用autograd自动计算numpy代码子集上的渐变
  • 然后使用cyipot添加这些
• scipy.minimize with solvers SLSQP or COBYLA，可以为您做任何事情（SLSQP可以使用等式和不等式约束；只有COBYLA的不等式约束，其中通过x >= y+x <= y可以来模拟等式约束）

使用工具完成任务

完整的示例问题在Test Examples for Nonlinear Programming Codes中定义：

这里有一些代码，基于数值微分，解决你的测试问题，包括官方设置（函数，梯度，起点，边界…）

import numpy as np
import scipy.sparse as sps
import ipopt
from scipy.optimize import approx_fprime


class Problem40(object):
    """ # Hock & Schittkowski test problem #40
            Basic structure  follows:
            - cyipopt example from https://pythonhosted.org/ipopt/tutorial.html#defining-the-problem
            - which follows ipopt's docs from: https://www.coin-or.org/Ipopt/documentation/node22.html
            Changes:
            - numerical-diff using scipy for function & constraints
            - removal of hessian-calculation
              - we will use limited-memory approximation
                - ipopt docs: https://www.coin-or.org/Ipopt/documentation/node31.html
              - (because i'm too lazy to reason about the math; lagrange and co.)
    """
    def __init__(self):
        self.num_diff_eps = 1e-8  # maybe tuning needed!

    def objective(self, x):
        # callback for objective
        return -np.prod(x)  # -x1 x2 x3 x4

    def constraint_0(self, x):
        return np.array([x[0]**3 + x[1]**2 -1])

    def constraint_1(self, x):
        return np.array([x[0]**2 * x[3] - x[2]])

    def constraint_2(self, x):
        return np.array([x[3]**2 - x[1]])

    def constraints(self, x):
        # callback for constraints
        return np.concatenate([self.constraint_0(x),
                               self.constraint_1(x),
                               self.constraint_2(x)])

    def gradient(self, x):
        # callback for gradient
        return approx_fprime(x, self.objective, self.num_diff_eps)

    def jacobian(self, x):
        # callback for jacobian
        return np.concatenate([
            approx_fprime(x, self.constraint_0, self.num_diff_eps),
            approx_fprime(x, self.constraint_1, self.num_diff_eps),
            approx_fprime(x, self.constraint_2, self.num_diff_eps)])

    def hessian(self, x, lagrange, obj_factor):
        return False  # we will use quasi-newton approaches to use hessian-info

    # progress callback
    def intermediate(
            self,
            alg_mod,
            iter_count,
            obj_value,
            inf_pr,
            inf_du,
            mu,
            d_norm,
            regularization_size,
            alpha_du,
            alpha_pr,
            ls_trials
            ):

        print("Objective value at iteration #%d is - %g" % (iter_count, obj_value))

# Remaining problem definition; still following official source:
# http://www.ai7.uni-bayreuth.de/test_problem_coll.pdf

# start-point -> infeasible
x0 = [0.8, 0.8, 0.8, 0.8]

# variable-bounds -> empty => np.inf-approach deviates from cyipopt docs!
lb = [-np.inf, -np.inf, -np.inf, -np.inf]
ub = [np.inf, np.inf, np.inf, np.inf]

# constraint bounds -> c == 0 needed -> both bounds = 0
cl = [0, 0, 0]
cu = [0, 0, 0]

nlp = ipopt.problem(
            n=len(x0),
            m=len(cl),
            problem_obj=Problem40(),
            lb=lb,
            ub=ub,
            cl=cl,
            cu=cu
            )

# IMPORTANT: need to use limited-memory / lbfgs here as we didn't give a valid hessian-callback
nlp.addOption(b'hessian_approximation', b'limited-memory')
x, info = nlp.solve(x0)
print(x)
print(info)

# CORRECT RESULT & SUCCESSFUL STATE

输出：

^{pr2}$

关于代码

• 我们使用scipy's approx_fprime，它基本上是为scipy.optimize中所有基于梯度的优化器添加的
• 正如源代码中所述；我没有考虑ipopt对hessian的需求，我们使用ipopts hessian-approximation
  • 基本思想在wiki: LBFGS中描述
• 我忽略了ipopts对约束Jacobian的稀疏结构的需要
  • 一个默认假设：默认的hessian结构是一个较低的三角形矩阵，使用了，我不会对这里可能发生的情况给出任何保证（性能差vs.破坏一切）