NumPy/Scipy中非线性函数的数值梯度

4
我正在尝试在numpy中实现数值梯度计算,以用作cyipopt中梯度的回调函数。 我对numpy梯度函数的理解是,它应返回基于有限差分逼近法在某一点计算的梯度。
我不明白如何使用此模块实现非线性函数的梯度。所提供的示例问题似乎是一个线性函数。
>>> 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 ])

我的代码片段如下:

import numpy as np

# Hock & Schittkowski test problem #40
x = np.mgrid[0.75:0.85:0.01, 0.75:0.8:0.01, 0.75:0.8:0.01, 0.75:0.8:0.01]
# target is evaluation at x = [0.8, 0.8, 0.8, 0.8]
f = -x[0] * x[1] * x[2] * x[3]
g = np.gradient(f)

print g

另一个缺点是我必须在多个点上评估x(并且它会在多个点返回梯度)。在numpy/scipy中是否有更好的选项,可以在单个点上数值地评估梯度,以便我可以将其实现为回调函数?
2个回答

6

首先,一些警告:

  • 数值优化很难做好
  • ipopt 是非常复杂的软件
    • 将 ipopt 与数值微分结合起来听起来像是在自讨苦吃,但当然这取决于你的问题
    • ipopt 几乎总是基于自动微分工具而不是数值微分

还有一些:

  • 由于这是一个复杂的任务,而且 python + ipopt 的状态不如其他一些语言(例如julia + JuMP),所以需要花费一些努力

还有一些备选方案:

  • 使用pyomo,它包装了ipopt并具有自动微分功能
  • 使用casadi,它也包装了ipopt并具有自动微分功能
  • 使用autograd自动计算numpy代码子集上的梯度
    • 然后使用cyipopt添加这些
  • 使用求解器SLSQP或COBYLA的scipy.minimize,它可以为您完成所有工作(SLSQP可以使用等式和不等式约束;COBYLA仅使用不等式约束,通过x >= y+x <= y模拟等式约束也可以起作用)

使用您的工具来处理任务

您的完整示例问题在非线性规划代码的测试示例中进行了定义:

enter image description here

这里是基于数值微分的代码,解决了您的测试问题,包括官方设定(函数、梯度、起始点、边界等)。

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

输出:

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.8, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:       12
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:        4
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        3
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

Objective value at iteration #0 is - -0.4096
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0 -4.0960000e-01 2.88e-01 2.53e-02   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
Objective value at iteration #1 is - -0.255391
   1 -2.5539060e-01 1.28e-02 2.98e-01 -11.0 2.51e-01    -  1.00e+00 1.00e+00h  1
Objective value at iteration #2 is - -0.249299
   2 -2.4929898e-01 8.29e-05 3.73e-01 -11.0 7.77e-03    -  1.00e+00 1.00e+00h  1
Objective value at iteration #3 is - -0.25077
   3 -2.5076955e-01 1.32e-03 3.28e-01 -11.0 2.46e-02    -  1.00e+00 1.00e+00h  1
Objective value at iteration #4 is - -0.250025
   4 -2.5002535e-01 4.06e-05 1.93e-02 -11.0 4.65e-03    -  1.00e+00 1.00e+00h  1
Objective value at iteration #5 is - -0.25
   5 -2.5000038e-01 6.57e-07 1.70e-04 -11.0 5.46e-04    -  1.00e+00 1.00e+00h  1
Objective value at iteration #6 is - -0.25
   6 -2.5000001e-01 2.18e-08 2.20e-06 -11.0 9.69e-05    -  1.00e+00 1.00e+00h  1
Objective value at iteration #7 is - -0.25
   7 -2.5000000e-01 3.73e-12 4.42e-10 -11.0 1.27e-06    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 7

                                   (scaled)                 (unscaled)
Objective...............:  -2.5000000000225586e-01   -2.5000000000225586e-01
Dual infeasibility......:   4.4218750883118219e-10    4.4218750883118219e-10
Constraint violation....:   3.7250202922223252e-12    3.7250202922223252e-12
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   4.4218750883118219e-10    4.4218750883118219e-10


Number of objective function evaluations             = 8
Number of objective gradient evaluations             = 8
Number of equality constraint evaluations            = 8
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 8
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 0
Total CPU secs in IPOPT (w/o function evaluations)   =      0.016
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
[ 0.79370053  0.70710678  0.52973155  0.84089641]
{'x': array([ 0.79370053,  0.70710678,  0.52973155,  0.84089641]), 'g': array([  3.72502029e-12,  -3.93685085e-13,   5.86974913e-13]), 'obj_val': -0.25000000000225586, 'mult_g': array([ 0.49999999, -0.47193715,  0.35355339]), 'mult_x_L': array([ 0.,  0.,  0.,  0.]), 'mult_x_U': array([ 0.,  0.,  0.,  0.]), 'status': 0, 'status_msg': b'Algorithm terminated successfully at a locally optimal point, satisfying the convergence tolerances (can be specified by options).'}

代码备注

  • 我们使用scipy的approx_fprime,这基本上是为了所有基于梯度的优化器在scipy.optimize中添加的。
  • 如来源所述;我没有考虑到ipopt需要海森矩阵,而我们使用了ipopt的hessian-approximation
  • 我忽略了ipopt对约束雅可比稀疏结构的需求
    • 默认假设:使用默认的海森矩阵结构为下三角矩阵,我不会对此作出任何保证(性能差 vs. 破坏一切)

需要进行分析。当存在唯一解时(这应该是情况;我提到了一个可能的问题和一个参数来调整!),不应该发生这种情况!人们在定义这些时经常犯错误。我不确定你是否犯了错误,但我只是提一下。你的退出状态是否正常? - sascha
我的输出与你的完全相同。我唯一收到的警告与缺少库有关:无法加载'libmetis.so' - 使用备用的AMD代码 - stagermane
这个警告只涉及性能。嗯...鉴于这些信息,我们没有太多可以分析的。 - sascha
1
好的,进一步检查后,我发现在我的显式雅可比计算中有一个打字错误。现在,使用两种方法得到了相同的结果,根据此网站的预期Objective=-0.25000000008228956 - stagermane
1
很高兴听到这个消息!只有一件事情需要注意:当不使用LBFGS时,要小心你的显式雅可比方法!你不需要纯雅可比矩阵,而是需要一个带有拉格朗日项的雅可比矩阵(根据IPOPT文档)。现在,在许多情况下,这个最后一项为零,但这取决于你的问题(我认为numdifftools不关心这个)! - sascha
显示剩余2条评论

0

我觉得你对数学函数的定义以及其数值实现有一些误解。

你应该将你的函数定义为:

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的定义隐藏,因此不允许在不同点进行多个梯度评估。此外,使用您的方法使您依赖于您离散化步骤。


在我看来,当具有函数时,导数发生;但是,如果只有一组变量(尚未适应任何函数),我们可以计算该数组中所有变量的偏导数,就像是“从这一点”开始,它们本质上定义了雅可比矩阵(梯度)...顺便说一句,用于拟合函数的系数...因此,为了计算任何目标或任何约束的梯度,可以使用numerical_grad = optimize.approx_fprime(inputs0_, func, epsilon)...而所有这些梯度一起形成了方程系统(包括目标和约束)的雅可比矩阵 - JeeyCi
方程系统:jac = lambda x0, *args, **kwargs: approx_fprime(x0, fun, eps, *args, **kwargs)来自示例)... 并使用np.allclosenp.testing.assert_almost_equal()相对于解析计算进行测试。 - JeeyCi

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接