2.7. Mathematical optimization: finding minima of functions

作者GaëlVaroquaux

数学优化处理查找函数的数字最小值(或最大值或零)的问题。在这种情况下,该函数称为成本函数目标函数能量

在这里,我们有兴趣使用scipy.optimize进行黑盒优化:我们不依赖于我们正在优化的函数的数学表达式。请注意,此表达式通常可用于更高效,非黑盒,优化。

先决条件

  • Numpy,Scipy
  • matplotlib

也可以看看

参考

数学优化非常...数学。如果你想要表演,真的很想读书:

2.7.1. 知道您的问题

不是所有的优化问题都是一样的了解您的问题使您能够选择正确的工具。

问题的维度

优化问题的规模几乎由问题的维度设置,即执行搜索的标量变量的数量。

2.7.1.1. 凸/非优化

凸函数

  • f高于其所有切线。
  • 对于两个点A,B,f(C)位于段[f(A),f(B])之下],如果A
非凸函数

优化凸函数很容易。优化非凸函数可能非常困难。

注意

可以证明,对于凸函数,局部最小值也是全局最小值。然后,在某种意义上,最小值是唯一的。

2.7.1.2. 平滑和不平滑的问题

平滑函数

渐变是随处定义的,是一个连续函数

非平滑函数

优化平滑函数更容易(在黑盒优化的上下文中为true,否则线性规划分段线性函数)。

2.7.1.3. 嘈杂与精确成本函数

嘈杂(蓝色)和无噪声(绿色)功能

噪音梯度

许多优化方法依赖于目标函数的梯度。如果未给出梯度函数,则它们被数值计算,这导致误差。在这种情况下,即使目标函数没有噪声,基于梯度的优化也可能是有噪声的优化。

2.7.1.4. 约束

约束下的优化

这里:

2.7.2. 对不同优化程序的回顾

2.7.2.1. 入门:1D优化

使用scipy.optimize.brent()可最小化1D函数。它将包围策略与抛物线近似相结合。

布林特对二次函数的方法:它在3次迭代中收敛,因为二次近似是精确的。
Brent的非凸函数方法:注意,优化器避免了局部最小值的事实是运气问题。
>>> from scipy import optimize
>>> def f(x):
... return -np.exp(-(x - .7)**2)
>>> x_min = optimize.brent(f) # It actually converges in 9 iterations!
>>> x_min
0.699999999...
>>> x_min - .7
-2.1605...e-10

注意

Brent的方法也可以用于优化使用scipy.optimize.fminbound()限制到间隔

注意

在scipy 0.11中,scipy.optimize.minimize_scalar()给出了1D标量最小化的通用接口

2.7.2.2. 基于渐变的方法

2.7.2.2.1. Some intuitions about gradient descent

这里我们关注直觉,而不是代码。代码随后。

梯度下降基本上包括沿梯度方向,即最速下降的方向采取小步长。

Fixed step gradient descent
良好条件二次函数。

病态条件二次函数。

梯度方法对病态问题的核心问题是梯度往往不指向最小值的方向。

我们可以看到,各向异性(病态)函数难以优化。

回家消息:条件数字和预处理

如果你知道你的变量的自然缩放,预定标它们,使它们的行为相似。这与预条件相关。

此外,显然可采取更大的步骤是有利的。这是使用线搜索在梯度下降代码中完成的。

Adaptive step gradient descent
良好的二次函数。
一个病态条件二次函数。
病态非二次函数。
病态非常非二次函数。

函数看起来像一个二次函数(椭圆等曲线),它越容易优化。

2.7.2.2.2. Conjugate gradient descent

上述梯度下降算法是不用于真实问题的玩具。

从上述实验可以看出,简单梯度下降算法的问题之一是它倾向于跨越谷而振荡,每次沿着梯度的方向,使其穿过谷。共轭梯度通过添加摩擦项来解决这个问题:每个步骤取决于梯度的最后两个值,并且急转弯减少。

Conjugate gradient descent
病态非二次函数。
病态非常非二次函数。

基于共轭梯度的方法在scipy中命名为'cg'使函数最小化的简单共轭梯度法是scipy.optimize.fmin_cg()

>>> def f(x):   # The rosenbrock function
... return .5*(1 - x[0])**2 + (x[1] - x[0]**2)**2
>>> optimize.fmin_cg(f, [2, 2])
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 13
Function evaluations: 120
Gradient evaluations: 30
array([ 0.99998968, 0.99997855])

这些方法需要函数的梯度。他们可以计算它,但会执行更好,如果你可以传递他们的渐变:

>>> def fprime(x):
... return np.array((-2*.5*(1 - x[0]) - 4*x[0]*(x[1] - x[0]**2), 2*(x[1] - x[0]**2)))
>>> optimize.fmin_cg(f, [2, 2], fprime=fprime)
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 13
Function evaluations: 30
Gradient evaluations: 30
array([ 0.99999199, 0.99998336])

注意,该函数仅被评估30次,与没有梯度的120相比。

2.7.2.3. 牛顿和准牛顿法

2.7.2.3.1. Newton methods: using the Hessian (2nd differential)

牛顿方法使用局部二次近似来计算跳转方向。为此,它们依赖于函数的二阶导数:gradientHessian

病态条件二次函数:

注意,由于二次近似是精确的,牛顿法是快速的

病态非二次函数:

这里我们优化高斯,它总是低于其二次近似。结果,牛顿法过调并导致振荡。

病态非常非二次函数:

在scipy中,用于优化的牛顿方法在scipy.optimize.fmin_ncg()中实现(cg这里指的是内部操作,即Hessian的反演由共轭梯度执行的事实)。scipy.optimize.fmin_tnc() can be use for constraint problems, although it is less versatile:

>>> def f(x):   # The rosenbrock function
... return .5*(1 - x[0])**2 + (x[1] - x[0]**2)**2
>>> def fprime(x):
... return np.array((-2*.5*(1 - x[0]) - 4*x[0]*(x[1] - x[0]**2), 2*(x[1] - x[0]**2)))
>>> optimize.fmin_ncg(f, [2, 2], fprime=fprime)
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 9
Function evaluations: 11
Gradient evaluations: 51
Hessian evaluations: 0
array([ 1., 1.])

注意,与共轭梯度(上面)相比,牛顿法需要更少的函数评价,但更多的梯度评价,因为它使用它来近似Hessian。让我们计算Hessian并将其传递给算法:

>>> def hessian(x): # Computed with sympy
... return np.array(((1 - 4*x[1] + 12*x[0]**2, -4*x[0]), (-4*x[0], 2)))
>>> optimize.fmin_ncg(f, [2, 2], fprime=fprime, fhess=hessian)
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 9
Function evaluations: 11
Gradient evaluations: 19
Hessian evaluations: 9
array([ 1., 1.])

注意

在非常高维度上,Hessian的反演可能是昂贵且不稳定的(大规模> 250)。

注意

牛顿优化器不应该与牛顿的根查找方法混淆,基于相同的原则,scipy.optimize.newton()

2.7.2.3.2. Quasi-Newton methods: approximating the Hessian on the fly

BFGS:BFGS(Broyden-Fletcher-Goldfarb-Shanno算法)在每个步骤精炼Hessian的近似。

病态条件二次函数:

在一个精确的二次函数上,BFGS不像牛顿的方法那么快,但仍然非常快。

病态非二次函数:

这里BFGS比牛顿更好,因为它的曲率的经验估计好于由Hessian给出的。

病态非常非二次函数:
>>> def f(x):   # The rosenbrock function
... return .5*(1 - x[0])**2 + (x[1] - x[0]**2)**2
>>> def fprime(x):
... return np.array((-2*.5*(1 - x[0]) - 4*x[0]*(x[1] - x[0]**2), 2*(x[1] - x[0]**2)))
>>> optimize.fmin_bfgs(f, [2, 2], fprime=fprime)
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 16
Function evaluations: 24
Gradient evaluations: 24
array([ 1.00000017, 1.00000026])

L-BFGS:有限内存BFGS位于BFGS和共轭梯度之间:在非常高的维数(> 250)中,Hessian矩阵太贵而不能计算和反转。L-BFGS保持低秩版本。此外,scipy版本scipy.optimize.fmin_l_bfgs_b()包括框边界:

>>> def f(x):   # The rosenbrock function
... return .5*(1 - x[0])**2 + (x[1] - x[0]**2)**2
>>> def fprime(x):
... return np.array((-2*.5*(1 - x[0]) - 4*x[0]*(x[1] - x[0]**2), 2*(x[1] - x[0]**2)))
>>> optimize.fmin_l_bfgs_b(f, [2, 2], fprime=fprime)
(array([ 1.00000005, 1.00000009]), 1.4417677473011859e-15, {...})

注意

如果您没有指定梯度到L-BFGS求解器,您需要添加approx_grad = 1

2.7.2.4. 无渐变的方法

2.7.2.4.1. A shooting method: the Powell algorithm

几乎是一个渐变的方法

病态条件二次函数:

鲍威尔的方法对于低维度的局部不适应不太敏感

病态非常非二次函数:

2.7.2.4.2. Simplex method: the Nelder-Mead

Nelder-Mead算法是对高维空间的二分法的泛化。该算法通过细化单纯形(将间隔和三角形广义化为高维空间)以将最小值括起来来工作。

强点:它对噪声是鲁棒的,因为它不依赖于计算梯度。因此,它可以工作在不是局部平滑的功能,例如实验数据点,只要它们显示大规模钟形行为。然而,它比基于梯度的平滑,无噪声函数的方法慢。

病态非二次函数:
病态非常非二次函数:

在scipy中,scipy.optimize.fmin()实现Nelder-Mead方法:

>>> def f(x):   # The rosenbrock function
... return .5*(1 - x[0])**2 + (x[1] - x[0]**2)**2
>>> optimize.fmin(f, [2, 2])
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 46
Function evaluations: 91
array([ 0.99998568, 0.99996682])

2.7.2.5. 全局优化程序

如果您的问题不允许唯一的本地最小值(除非函数是凸的,否则很难测试),并且您没有先前的信息来初始化靠近解决方案的优化,您可能需要一个全局优化器。

2.7.3. 使用scipy 的实用指南

2.7.3.1. 选择方法

../../_images/sphx_glr_plot_compare_optimizers_001.png
没有梯度的知识:
 
有梯度知识:
 
  • BFGS(scipy.optimize.fmin_bfgs())或L-BFGS(scipy.optimize.fmin_l_bfgs_b())。
  • BFGS的计算开销大于L-BFGS,其本身大于共轭梯度的。另一方面,BFGS通常比CG需要更少的功能评估。因此,共轭梯度法在优化计算廉价函数方面优于BFGS。
与黑森州:
 
如果您有噪声测量:
 

2.7.3.2. 让您的优化工具更快速

  • 选择正确的方法(见上文),做分析计算梯度和Hessian,如果可以。
  • 尽可能使用preconditionning
  • 明智地选择您的初始化点。例如,如果您正在运行许多类似的优化,请使用另一个的结果热启动一个。
  • 如果不需要精度,请放宽公差

2.7.3.3. 计算梯度

计算梯度,甚至更多的Hessians,是非常乏味,但值得的努力。使用Sympy的符号计算可能会派上用场。

警告

非常优化的不收敛好的公共来源是梯度计算中的人为误差。您可以使用scipy.optimize.check_grad()来检查渐变是否正确。它返回给定梯度之间的不同的范数,以及数字计算的梯度:

>>> optimize.check_grad(f, fprime, [2, 2])
2.384185791015625e-07

另请参见scipy.optimize.approx_fprime()以查找错误。

2.7.3.4. 合成练习

../../_images/sphx_glr_plot_exercise_ill_conditioned_001.png

练习:一个简单的(?)二次函数

优化以下函数,使用K [0]作为起点:

np.random.seed(0)
K = np.random.normal(size=(100, 100))
def f(x):
return np.sum((np.dot(K, x - 1))**2) + np.sum(x**2)**2

时间你的方法。找到最快的方法。为什么BFGS不能很好地工作?

练习:局部平坦的最小值

考虑函数exp(-1 /(1 * x ** 2 + y ** 2)此函数允许(0,0)中的最小值。从在(1,1)处的初始化开始,尝试在该最小点的1e-8内。

flat_min_0 flat_min_1

2.7.4. 特殊情况:非线性最小二乘法

2.7.4.1. 最小化向量函数的范数

最小二乘问题最小化向量函数的范数,具有可以在scipy.optimize.leastsq()中实现的Levenberg-Marquardt算法中使用的特定结构。

让我们尝试最小化以下向量函数的范数:

>>> def f(x):
... return np.arctan(x) - np.arctan(np.linspace(0, 1, len(x)))
>>> x0 = np.zeros(10)
>>> optimize.leastsq(f, x0)
(array([ 0. , 0.11111111, 0.22222222, 0.33333333, 0.44444444,
0.55555556, 0.66666667, 0.77777778, 0.88888889, 1. ]), 2)

这需要67个函数求值(用'full_output = 1'检查)。如果我们自己计算规范,并使用一个好的通用优化器(BFGS)

>>> def g(x):
... return np.sum(f(x)**2)
>>> optimize.fmin_bfgs(g, x0)
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 11
Function evaluations: 144
Gradient evaluations: 12
array([ -7.4...-09, 1.1...e-01, 2.2...e-01,
3.3...e-01, 4.4...e-01, 5.5...e-01,
6.6...e-01, 7.7...e-01, 8.8...e-01,
1.0...e+00])

BFGS需要更多的函数调用,并给出不太精确的结果。

注意

leastsq仅当输出向量的维数较大时比BFGS有意义,并且大于要优化的参数的数量。

警告

如果函数是线性的,这是一个线性代数问题,应该用scipy.linalg.lstsq()来解决。

2.7.4.2. 曲线拟合

../../_images/sphx_glr_plot_curve_fit_001.png

最小二乘问题在将非线性拟合到数据时经常发生。虽然可以自己构造我们的优化问题,但是scipy为此提供了一个帮助函数:scipy.optimize.curve_fit()

>>> def f(t, omega, phi):
... return np.cos(omega * t + phi)
>>> x = np.linspace(0, 3, 50)
>>> y = f(x, 1.5, 1) + .1*np.random.normal(size=50)
>>> optimize.curve_fit(f, x, y)
(array([ 1.51854577, 0.92665541]), array([[ 0.00037994, -0.00056796],
[-0.00056796, 0.00123978]]))

练习

对omega = 3也一样。有什么困难?

2.7.5. 使用约束的优化

2.7.5.1. 框边界

框边界对应于限制优化的每个单独参数。注意,一些最初不写为框边界的问题可以通过变量的改变来重写。

../../_images/sphx_glr_plot_constraints_002.png

2.7.5.2. 一般约束

指定为函数的等式和不等式约束:f(x)= 0g(x) 0

  • scipy.optimize.fmin_slsqp()顺序最小二乘法规划:等式和不等式约束:

    ../../_images/sphx_glr_plot_non_bounds_constraints_001.png
    >>> def f(x):
    
    ... return np.sqrt((x[0] - 3)**2 + (x[1] - 2)**2)
    >>> def constraint(x):
    ... return np.atleast_1d(1.5 - np.sum(np.abs(x)))
    >>> optimize.fmin_slsqp(f, np.array([0, 0]), ieqcons=[constraint, ])
    Optimization terminated successfully. (Exit mode 0)
    Current function value: 2.47487373504
    Iterations: 5
    Function evaluations: 20
    Gradient evaluations: 5
    array([ 1.25004696, 0.24995304])
  • scipy.optimize.fmin_cobyla()线性近似的约束优化:仅限不等式约束:

    >>> optimize.fmin_cobyla(f, np.array([0, 0]), cons=constraint)
    
    array([ 1.25009622, 0.24990378])

警告

上述问题被称为统计中的Lasso问题,并且它存在非常有效的求解器(例如在scikit-learn中)。一般来说,当特定的存在时,不要使用通用解算器。

拉格朗日乘子

如果你准备做一点数学,许多约束优化问题可以使用称为拉格朗日乘子的数学技巧转换为非约束优化问题。