方程的极值问题

詹令
lealzhan@126.com
2017.8.14
last modified 2017.10.30

Unconstrained

多维无约束优化问题:

min ⁡ x ∈ R n f ( x ) \min_{x \in R^n} f(x) xRnminf(x)

最优化方法主要是基于迭代方法。 迭代方法要求迭代过程中目标函数是下降的:

f ( x k + 1 ) < f ( x k ) f(x_k+1)<f(x_k) f(xk+1)<f(xk)

迭代方法的基本要素

  • 搜索方向
  • 步长
  • 迭代格式
  • 终止条件

迭代方法的基本流程

两大类迭代方法

  • 线搜索
    先计算搜索方向,再计算步长
    最速下降,共轭梯度,牛顿法,伪牛顿法,
    又分为精确线搜索和非精确线搜索。

    • 精确线搜索一般不用
    • 非精确线搜索又分为好多种。
    • Arjimo
    • Goldstein
  • 信赖域
    先计算步长,再计算搜索方向

1. Line Search 线搜索

1.1 基于解析梯度

1.1.1 Gradient Descent/Stepest Descent:
1.1.1.1 特性
  • 多维无约束基于梯度的优化算法中最简单,可靠的方法。
  • 收敛速度较慢,一次收敛
  • 越接近越接近目标,步长越小。步长选择对结果影响很大。BB步长?
  • 可用于求解非线性方程组?
1.1.1.2 算法

第k步迭代时的搜索方向

min ⁡ d k f ( x k + d k ) \min_{d_k} f(x_k + d_k) dkminf(xk+dk)

f ( x k + d k ) f(x_k + d_k) f(xk+dk) 一阶泰勒展开:

f ( x k + d k ) = f ( x k ) + ∇ f ( x k ) T d k f(x_k + d_k) = f(x_k) + \nabla f(x_k)^T d_k f(xk+dk)=f(xk)+f(xk)Tdk

等价于求

min ⁡ d k ( f ( x k ) + ∇ f ( x k ) T d k ) \min_{d_k} \left(f(x_k) + \nabla f(x_k)^T d_k \right) dkmin(f(xk)+f(xk)Tdk)

==

d k d_k dk使其在 ∇ f ( x k ) ) \nabla f(x_k)) f(xk))上的投影最小. 显然,当 d k = − ∇ f ( x k ) \color{red}{d_k = -\nabla f(x_k)} dk=f(xk) 时, 投影最小(-1)。

x的搜索方向

d k = − ∇ f ( x k ) \color{red}{d_k= -\nabla f(x_k) } dk=f(xk)

步长 α k \alpha_k αk 由一维搜索来确定:

f ( x k + α k d k ) = min ⁡ α > 0 f ( x k + α d k ) f(x_k + \alpha_k d_k) = \min_{\alpha > 0} f(x_k + \alpha d_k) f(xk+αkdk)=α>0minf(xk+αdk)

迭代格式

x k + 1 = x k + α k d k x_{k+1} = x_k + \alpha_k d_k xk+1=xk+αkdk

终止条件

Δ f ( x ) < ε \Delta f(x)<\varepsilon Δf(x)<ε

或者

∣ ∣ d k ∣ ∣ 2 < ε ||d_k||_2 < \varepsilon dk2<ε

?以上的梯度如果是基于f的解析梯度求得,则是最原始的梯度下降。如果是基于f(x)的样本点求得,则分为以下三种情况:批量梯度下降法, 随机梯度下降,小批量梯度下降法。?

1.1.1.3 困难点及解决方案
1.1.1.4 变种
1.1.1.4.1 批量梯度下降法

(Batch Gradient Descent, BGD)

1.1.1.4.2 随机梯度下降

(Stochastic gradient descent, SGD)
http://www.cnblogs.com/maybe2030/p/5089753.html
相对于MGD更常用些。
最早是在Tom Mitchel的机器学习[]的讲人工神经网络时看到的.

1.1.1.4.3 小批量梯度下降法

(Mini-batch Gradient Descent,MBGD)

1.1.1.5 并行的问题
1.1.2 Linear Conjugate Gradient
1.1.2.1 特性
1.1.2.2 算法
1.1.2.3 困难点
1.1.2.4 解决方案

目标函数是凸/正定二次函数的情形,求以下函数 f ( x ) f(x) f(x)的最小值:

f ( x ) =   g T x + 1 2 x T H x f(x) = \ g^Tx + \frac{1}{2}x^THx f(x)= gTx+21xTHx

其中 g ∈ R n , H ∈ R n × n 对 称 正 定 g \in R^n, H \in R^{n \times n}对称正定 gRn,HRn×n
等价于求 H x + g = 0 Hx+g=0 Hx+g=0.

let g k = ∇ f ( x k ) = H x k + g g_k = \nabla f(x_k) = Hx_k+g gk=f(xk)=Hxk+g

步长

α k = − g k T d k d k T H d k \alpha_k = -\frac{g_k^Td_k}{d_k^THd_k} αk=dkTHdkgkTdk

迭代格式

x k + 1 = x k + α k d k x_{k+1} = x_k + \alpha_k d_k xk+1=xk+αkdk

搜索方向

d k + 1 = − g k + 1 + β k d k d_{k+1} = - g_{k+1} + \beta_k d_{k} dk+1=gk+1+βkdk

其中 β k = g k + 1 T H d k d k T H d k \beta_k = \frac{g_{k+1}^T Hd_k}{d_k^THd_k} βk=dkTHdkgk+1THdk
或者 β k = g k + 1 T H d k d k T H d k \beta_k = \frac{g_{k+1}^T Hd_k}{d_k^THd_k} βk=dkTHdkgk+1THdk

preconditioned conjugate gradient?

1.1.3 None-linear Conjugate Gradient methods

目标函数是的一般函数的情形, 求以下函数 f ( x ) f(x) f(x)的最小值:

f ( x ) f(x) f(x)

β k \beta_k βk的选取和Linear CG不同

storage of order N

1.1.3.1 Fletcher-Reeves algorithm

the original none-linear cg

1.1.3.2 Polak-Ribiere algorithm

need 1st derivative
new gamma_i
(scipy.optimize.minimize “CG”)

搜索方向

1.1.4 Newton Raphson Method
特性
  • 要求目标函数二次可微
  • 要求Hessian矩阵正定
  • 迭代求求解时 二次收敛
算法

第k步迭代时的搜索方向

min ⁡ d k f ( x k + d k ) \min_{d_k} f(x_k + d_k) dkminf(xk+dk)

f ( x k + d k ) f(x_k + d_k) f(xk+dk) 一阶泰勒展开:

f ( x k + d k ) = f ( x k ) + ∇ f ( x k ) T d k f(x_k + d_k) = f(x_k) + \nabla f(x_k)^T d_k f(xk+dk)=f(xk)+f(xk)Tdk

等价于求

min ⁡ d k ( f ( x k ) + ∇ f ( x k ) T d k ) \min_{d_k} \left(f(x_k) + \nabla f(x_k)^T d_k \right) dkmin(f(xk)+f(xk)Tdk)

==

d k d_k dk使其在 ∇ f ( x k ) ) \nabla f(x_k)) f(xk))上的投影最小. 显然,当 d k = − ∇ f ( x k ) \color{red}{d_k = -\nabla f(x_k)} dk=f(xk) 时, 投影最小(-1)。

f ( x k + d k ) f(x_k + d_k) f(xk+dk) 二阶泰勒展开:

f ( x k + d k ) = f ( x k ) + ∇ f ( x k ) T d k + 1 2 d k T ∇ 2 f ( x k ) d k f(x_k + d_k) = f(x_k) + \nabla f(x_k)^T d_k + \frac{1}{2}d_k^T\nabla^2f(x_k)d_k f(xk+dk)=f(xk)+f(xk)Tdk+21dkT2f(xk)dk

等价于求

min ⁡ d k ( f ( x k ) + ∇ f ( x k ) T d k + 1 2 d k T ∇ 2 f ( x k ) d k ) \min_{d_k} \left( f(x_k) + \nabla f(x_k)^T d_k + \frac{1}{2}d_k^T\nabla^2f(x_k)d_k \right) dkmin(f(xk)+f(xk)Tdk+21dkT2f(xk)dk)

上式中的目标函数对 d k d_k dk的梯度是零:

∇ f ( x k ) T + ∇ 2 f ( x k ) d k = 0 \nabla f(x_k)^T + \nabla^2f(x_k)d_k = 0 f(xk)T+2f(xk)dk=0

d k = − [ ∇ 2 f ( x k ) ] − 1 ∇ f ( x k ) T \color{red}{d_k = -[\nabla^2f(x_k)]^{-1}\nabla f(x_k)^T } dk=[2f(xk)]1f(xk)T


求目标数 f ( x ) f(x) f(x)的最小值。基于迭代法,等价于求第k步迭代时目标函数在哪个方向上

min ⁡ d ∈ R n f ( x k + d ) \min_{d \in R^n} {f(x_k+d)} dRnminf(xk+d)

假定 f ( x ) f(x) f(x) 二次可微, 且 Hessian矩阵 ∇ 2 f ( x ) \nabla^2 f(x) 2f(x) 正定。
f ( x ) f(x) f(x)在当前迭代点 x k x_k xk上泰勒展开,

f ( x k + d ) ≈   f ( x k ) + d T ∇ f ( x k ) + 1 2 d T ∇ 2 f ( x k ) d f(x_k+d) \approx \ f(x_k) + d^T \nabla f(x_k) + \frac{1}{2}d^T \nabla^2 f(x_k)d f(xk+d) f(xk)+dTf(xk)+21dT2f(xk)d

⟹ \Longrightarrow

min ⁡ d ∈ R n ( f ( x k + d ) ) \min_{d \in R^n} {(f(x_k+d))} dRnmin(f(xk+d)) 等价于泰勒展开的右边项 f ( x k ) + d T ∇ f ( x k ) + 1 2 d T ∇ 2 f ( x k ) d f(x_k) + d^T \nabla f(x_k) + \frac{1}{2}d^T \nabla^2 f(x_k)d f(xk)+dTf(xk)+21dT2f(xk)d d d d的梯度为零

∇ f ( x k ) + ∇ 2 f ( x k ) d = 0 \nabla f(x_k) + \nabla^2 f(x_k)d = 0 f(xk)+2f(xk)d=0

因为 ∇ 2 f ( x ) \nabla^2 f(x) 2f(x) 正定[], 可得搜索方向

d k = − [ ∇ 2 f ( x k ) ] − 1 ∇ f ( x k ) \color{red}{d_k=-[\nabla^2 f(x_k)]^{-1}\nabla f(x_k) } dk=[2f(xk)]1f(xk)

步长 α k \alpha_k αk 由一维/线性搜索来确定:
待定
f ( x k + α k d k ) = min ⁡ α > 0 f ( x k + α d k ) f(x_k + \alpha_k d_k) = \min_{\alpha > 0} f(x_k + \alpha d_k) f(xk+αkdk)=α>0minf(xk+αdk)

迭代格式

x k + 1 = x k + α k d k x_{k+1} = x_k + \alpha_k d_k xk+1=xk+αkdk

终止条件
待定

Δ f ( x ) < ε \Delta f(x)<\varepsilon Δf(x)<ε

或者

∣ ∣ d k ∣ ∣ 2 < ε ||d_k||_2 < \varepsilon dk2<ε

困难点
  • Hessian矩阵难以求出。什么情况下难以求出?鞍点?
  • Hessian矩阵可以求出,但是Hessian矩阵的逆可能数值不稳定.
解决方案
  • 截尾牛顿法
  • 修正牛顿法
  • 拟牛顿法
并行的问题
1.1.5 Quasi-Newton/Pseudo-Newton/拟牛顿法
  • 超线性收敛 superlinear convergence

∇ 2 f ( x ) \nabla^2 f(x) 2f(x) 其实不一定正定(马鞍面?),所以Newton Raphson Method的 d k d_k dk不一定是最优方向。但是可以构造出一个正定的 ∇ 2 f ( x ) \nabla^2 f(x) 2f(x), 便有了Quasi-Newton?

f ( x k + d ) ≈   f ( x k ) + d T g k + 1 2 d T B k d f(x_k+d) \approx \ f(x_k)+d^T g_k + \frac{1}{2}d^T B_k d f(xk+d) f(xk)+dTgk+21dTBkd

其中 g k = ∇ f ( x k ) , B k R n × n g_k= \nabla f(x_k), B_k R^{n \times n} gk=f(xk),BkRn×n是对称矩阵。
可以发现 B k B_k Bk被用于近似Newton Raphson Method中的目标函数的海森矩阵 ∇ 2 f ( x k ) \nabla^2 f(x_k) 2f(xk).

mostly better than Conjugate gradient methods, storage of order N2

1.1.5.1 Davidon-Fletcher-Powell (DFP) algorithm
1.1.5.2 Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm
  • even good for non-smooth optimization
    • need 1st derivative(gradient) (analytically or numerically)
    • 2nd derivative(hessian matrix) is approximated and stored; =>? storage of order N2
    • better than DFP;
    • superlinear convergence(< quadritic)
    • (scipy.optimize.minimize “BFGS”)
    • 变种: L-BFGS (Low-Storage BFGS) (Limited-Memory BFGS)
      • [LINK] [Numerical Optimization 2nd 9.1]
      • 用于大规模问题
      • 只存储 少量长度为n的向量 来近似 n*n的Hessian矩阵 =>? storage of order N
      • linear convergency
      • (scipy.optimize.minimize “L-BFGS-B”)

H = ∇ 2 f ( x k ) ← B H = \nabla^2 f(x_k) \leftarrow B H=2f(xk)B

其中

B 0 = I B_0 = I B0=I

B k + 1 = B k + Δ k x Δ k T x Δ k T x Δ k θ + B k Δ k x ( B k Δ k x ) T Δ k T x B k Δ k x B_{k+1} = B_k + \frac{\Delta_kx \Delta_k^Tx}{\Delta_k^Tx \Delta_k\theta} + \frac{B_k \Delta_kx(B_k\Delta_kx)^T}{\Delta_k^Tx B_k \Delta_kx} Bk+1=Bk+ΔkTxΔkθΔkxΔkTx+ΔkTxBkΔkxBkΔkx(BkΔkx)T

搜索方向

d k = − B k − 1 ∇ f ( x k ) \color{red}{d_k= - B_k^{-1}\nabla f(x_k)} dk=Bk1f(xk)

拟牛顿法和牛顿法一样,都会碰到一个工程上的挑战:当x的维度很高时,H就会变得非常大。比方说x是一百万维,H就有一万亿个元素,一般计算机的内存可对付不了这么大的矩阵。于是,还要把H做个近似,将其分解为两个条状矩阵和一个对角阵,这样规模就大大降低了。如果在BFGS方法上应用这一技巧,就是著名的L-BFGS方法,这里的“L”,指的是“Limited-Memory”。
http://www.sohu.com/a/154817861_647918

梯度下降法和牛顿法/拟牛顿法相比,两者都是迭代求解,不过梯度下降法是梯度求解,而牛顿法/拟牛顿法是用二阶的海森矩阵的逆矩阵或伪逆矩阵求解。相对而言,使用牛顿法/拟牛顿法收敛更快。但是每次迭代的时间比梯度下降法长。 http://www.cnblogs.com/pinard/p/5970503.html

1.2 基于目标函数的有限差分式计算数值梯度

  • Quasi-Newton or variable metric
  • Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm
    (scipy.optimize.minimize “BFGS”, if deritive func is not given, use numerical deritive)

1.3 无梯度/直接法

无约束优化问题

min ⁡ x ∈ R n f ( x ) \min_{x \in R^n} f(x) xRnminf(x)

等价于求 ∇ f ( x ) = 0 ? \nabla f(x)=0 ? f(x)=0?.

1.3.1 Nelder-Mead/Downhill Simplex Method
1.3.1.1 特性
  • slow, robust, storage of order N2
  • (scipy.optimize.minimize “Nelder-Mead”)
  • 二维空间 三个点, 而爬N维空间里的山峰,就得N+1个人配合了
1.3.2 Powell’s Method/Direction Set Method
1.3.2.1 特性
  • quadratic convergence(?), robust, storage of order N2
  • (scipy.optimize.minimize “Powell”)

储存复杂度 N2和N,N为空间维度。(不是很理解) N代表优化变量的数目?

2. Trust Region 置信域法

说到这儿,细心的读者可能会发现一个问题,不论是梯度法、牛顿法、BFGS还是L-BFGS,都可以总结为“先方向,后步长”的方法。但是在找方向的过程中,又是编二阶导,又是近似的,这个方向找的还靠谱么?我只能说,大致上靠谱,不过确实也打了不小的折扣。那么如果不做这些妥协和近似行不行呢,办法也是存在的,不过这就要要变成“先步长,后方向”的思路。
先方向后步长方法的典型代表,是置信域法(Trust-region method),与前面的拟牛顿法相比,它不用近似一个正定的H——这里的妙处在于,只要用金箍棒划了一个圈儿,不管圈儿里是抛物面还是马鞍面,最高点是可以求出来的,而且有个接近闭式的解!
不管是先步长还是先方向,步长到底选多大合适呢?说实话,这基本上是个玄学问题(如果不是神学问题的话):步子迈小了,得多怎才能爬上去;步子迈大了,象毗湿奴那样一下子迈过三界的话,什么梯度、Hession矩阵就全乱了,优化过程也变得不可控。但是,规律还是有的,那就是步长总是要按一定的规律逐渐变小的,好比你罚点球助跑的时候,前面大步流星,后面碎步紧倒,才能准确找到球。 http://www.sohu.com/a/154817861_647918

Constrained:

有约束的优化问题,可以通过拉格朗日乘数法转化成无约束的方程组优化问题

Test Code

Test Code 0

# python 3.5
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize
from scipy.optimize import minimize, rosen, rosen_der
import numpy as np
#x0 = [1.2, 1.5, 0.39, 1.33]
#res = minimize(rosen, x0, method='Nelder-Mead', tol=1e-4)
#print res.x


def func(x0):
	return (x0[0]-1)*(x0[0]-1) + x0[1]*x0[1] + x0[2]*x0[2] + np.abs(x0[3]*x0[0]) + np.abs(x0[2]*x0[0])

def test():
	x0 = [6.9, 1.1, 1.1, 0.1]

	res = minimize(rosen, x0, method='Nelder-Mead', tol=1e-8)
	print ("Nelder-Mead")
	print ("result", res.x)
	print ("iteration numb", res.nit)
	print ("----------")

	res = minimize(rosen, x0, method='Powell', tol=1e-8)
	print ("Powell")
	print ("result", res.x)
	print ("iteration numb", res.nit)
	print ("----------")

	res = minimize(rosen, x0, method='CG', tol=1e-8)
	print ("CG")
	print ("result", res.x)
	print ("iteration numb", res.nit)
	print ("----------")

	res = minimize(rosen, x0, method='BFGS', tol=1e-8)
	print ("BFGS")
	print ("result", res.x)
	print ("iteration numb", res.nit)
	print ("----------")

	#need jacobian
	#res = minimize(rosen, x0, method='Newton-CG', tol=1e-8)
	#print ("Newton-CG", "result", res.x, "iteration numb", res.nit

	res = minimize(rosen, x0, method='L-BFGS-B', tol=1e-8)
	print ("L-BFGS-B")
	print ("result", res.x)
	print ("iteration numb", res.nit)
	print ("----------")

	res = minimize(rosen, x0, method='TNC', tol=1e-8)
	print ("TNC")
	print ("result", res.x)
	print ("iteration numb", res.nit)
	print ("----------")

	res = minimize(rosen, x0, method='COBYLA', tol=1e-8)
	print ("COBYLA")
	print ("result", res.x)
	print ("iteration numb")#, res.nit\
	print ("----------")

	res = minimize(rosen, x0, method='SLSQP', tol=1e-8)
	print ("SLSQP")
	print ("result", res.x)
	print ("iteration numb", res.nit)
	print ("----------")

	#need jacobian
	#res = minimize(rosen, x0, method='dogleg', tol=1e-8)
	#print ("dogleg", "result", res.x, "iteration numb"

	#need jacobian
	#res = minimize(rosen, x0, method='trust-ncg', tol=1e-8)
	#print ("trust-ncg", "result", res.x, "iteration numb"

if __name__=='__main__':
	test()

Note

Gradient, Jacobian, Hessian

正定矩阵 可逆?

多元函数极值和最小二乘问题的关系


Reference

[] Nocedal, Jorge, and S. J. Wright. “Numerical Optimization.” Springer 9.4(1999):1556-1556.
[] python调用 [[LINK](https://docs.scipy.org/doc/scipy/reference/generated/ scipy.optimize.minimize.html#scipy.optimize.minimize)]
[] Powell’s Method[LINK]
[] NLopt [LINK]
[] nonelinear conjugate gradient methods [wiki]
[] ZOOpt 无梯度优化工具包 [github]
[] stockoverflow : CG vs BFGS [LINK]
[] Comparison of derivative-free optimization algorithms [LINK]

[] 梯度下降法和共轭梯度法有何异同?[LINK]
[] http://www.sohu.com/a/154817861_647918

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值