最优化之最速下降法.

最优化问题——梯度下降法1、⽆约束最优化问题求解此问题的⽅法⽅法分为两⼤类：最优条件法和迭代法。

2、最优条件法我们常常就是通过这个必要条件去求取可能的极⼩值点，再验证这些点是否真的是极⼩值点。

当上式⽅程可以求解的时候，⽆约束最优化问题基本就解决了。

实际中，这个⽅程往往难以求解。

这就引出了第⼆⼤类⽅法：迭代法。

最优条件法：最⼩⼆乘估计3、迭代法（1）梯度下降法（gradient descent），⼜称最速下降法（steepest descent）梯度下降法是求解⽆约束最优化问题的⼀种最常⽤的⽅法。

梯度下降法是迭代算法，每⼀步需要求解⽬标函数的梯度向量。

必备条件：函数f(x)必须可微，也就是说函数f(x)的梯度必须存在优点：实现简单缺点：最速下降法是⼀阶收敛的，往往需要多次迭代才能接近问题最优解。

算法A.1（梯度下降法）输⼊：⽬标函数f(x)，梯度函数g(x)=▽f(x)，计算精度ε；输出：f(x)的极⼩点x*总结：选取适当的初值x(0),不断迭代，更新x的值，进⾏⽬标函数的极⼩化，直到收敛。

由于负梯度⽅向是使函数值下降最快的⽅向，在迭代的每⼀步，以负梯度⽅向更新x的值，从⽽达到减少函数值的⽬的。

λk叫步长或者学习率；梯度⽅向g k=g(x(k)）是x=x(k)时⽬标函数f(x)的⼀阶微分值。

学习率/步长λ的确定：当f(x)的形式确定，我们可以通过求解这个⼀元⽅程来获得迭代步长λ。

当此⽅程形式复杂，解析解不存在，我们就需要使⽤“⼀维搜索”来求解λ了。

⼀维搜索是⼀些数值⽅法，有0.618法、Fibonacci法、抛物线法等等，这⾥不详细解释了。

在实际使⽤中，为了简便，也可以使⽤⼀个预定义的常数⽽不⽤⼀维搜索来确定步长λ。

这时步长的选择往往根据经验或者通过试算来确定。

步长过⼩则收敛慢，步长过⼤可能震荡⽽不收敛。

如下图：当⽬标函数是凸函数时，梯度下降法的解是全局最优解。

但是，⼀般情况下，往往不是凸函数，所以其解不保证是全局最优解。

最优化Armijo算法确定步长的最速下降法资料最速下降法是最优化算法中最简单、最基础的一种方法，但其收敛速度较慢且容易陷入局部最优解。

因此，在最速下降法的基础上，可以通过引入步长的方法来提高算法的收敛速度。

而Armijo算法就是一种常见的用于确定步长的方法。

最速下降法基础假设我们要最小化目标函数f(x)，那么最速下降法的思路就是从一个初始点x0开始，不断朝着负梯度方向进行迭代，直到找到最优解x∗，即：$x_{k+1} = x_k - \\alpha_k \ abla f(x_k)$其中，ablaf(x k)是f(x)在x k处的梯度，$\\alpha_k$ 是步长（也称为学习率），表示每次迭代的步长大小。

但这里还有一个问题：如何确定每次迭代的步长呢？Armijo算法Armijo算法是一种基于梯度下降法的步长确定方法。

它的思路是，每次迭代的步长不应该过大，否则容易导致超出收敛区域。

同时，步长也不应该过小，否则收敛速度会变得非常缓慢。

因此，步长的大小应该恰到好处，即在一定范围内找到一个最优的步长大小。

具体地，Armijo算法通过二分搜索的方法，在可行步长范围内找到一个最优的步长 $\\alpha_k$。

具体过程如下：1.首先初始化 $\\alpha_0$，并设定一些参数，如尝试步长大小t、可行步长下界 $\\tau$ 和函数下降的最小比例 $\\gamma$。

2.计算目标函数f(x k−t ablaf(x k))，以及根据一定准则确定下一个$\\alpha$。

3.如果 $f(x_k - \\alpha_k\ abla f(x_k))$ 函数值比f(x k)减小了一些比例$\\gamma$，则认为当前 $\\alpha_k$ 是可行的步长。

4.如果当前 $\\alpha_k$ 不是可行的步长，则将其折半，即 $\\alpha_k\\leftarrow \\alpha_k/2$，直到找到一个可行的步长为止。

最优化算法最速下降法、⽜顿法、拟⽜顿法Python实现---------------------------------------2020.9.23更新---------------------------------把 BFGS(x)改写了⼀下，变简洁了def BFGS(x): #拟⽜顿法epsilon, h, maxiter = 10**-5, 10**-5, 10**4Bk = np.eye(x.size)for iter1 in range(maxiter):grad = num_grad(x, h)if np.linalg.norm(grad) < epsilon:return xdk = -np.dot((np.linalg.inv(Bk)), grad)ak = linesearch(x, dk)x = x + dk*akyk = num_grad(x, h) -gradsk = ak*dkif np.dot(yk, sk) > 0:Bs = np.dot(Bk,sk)ys = np.dot(yk,sk)sBs = np.dot(np.dot(sk,Bk),sk)Bk = Bk - 1.0*Bs.reshape((n,1))*Bs/sBs + 1.0*yk.reshape((n,1))*yk/ysreturn x---------------------------------------2020.9.23更新---------------------------------只⽤到了numpy这⼀个库，只要安装有这个库应该都可以直接运⾏import numpy as npdef f(x): #⽬标函数x1 = x[0]x2 = x[1]y = 100*((x2 - x1**2)**2) + (x1-1)**2return ydef num_grad(x, h): #求梯度df = np.zeros(x.size)for i in range(x.size):x1, x2 = x.copy(), x.copy() #这⾥需要⽤到复制，⽽不能⽤赋值号（=），原因是Python⾥⾯=号只是取别名，不是复制（c/c++⾥⾯是）x1[i] = x[i] - hx2[i] = x[i] + hy1, y2 = f(x1), f(x2)df[i] = (y2-y1)/(2*h)return dfdef num_hess(x, h): #求hess矩阵hess = np.zeros((x.size, x.size))for i in range(x.size):x1 = x.copy()x1[i] = x[i] - hdf1 = num_grad(x1, h)x2 = x.copy()x2[i] = x[i] + hdf2 = num_grad(x2, h)d2f = (df2 - df1) / (2 * h)hess[i] = d2freturn hessdef linesearch(x, dk): #求步长ak = 1for i in range(20):newf, oldf = f(x + ak * dk), f(x)if newf < oldf:return akelse:ak = ak / 4 #迭代更新步长，步长可随意变换，保证newf⽐oldf⼩就可以了（如改为: ak=ak/2 也是可以的）return akdef steepest(x): #最速下降法epsilon, h, maxiter = 10**-5, 10**-5, 10**4for iter1 in range(maxiter):grad = num_grad(x, h)if np.linalg.norm(grad) < epsilon:return xdk = -gradak = linesearch(x, dk)x = x + ak * dkreturn xdef newTonFuction(x): #⽜顿法epsilon, h1, h2, maxiter = 10**-5, 10**-5, 10**-5, 10**4for iter1 in range(maxiter):grad = num_grad(x, h1)if np.linalg.norm(grad) < epsilon:return xhess = num_hess(x, h2)dk = -np.dot((np.linalg.inv(hess)), grad)x = x + dkreturn xdef BFGS(x): #拟⽜顿法epsilon, h, maxiter = 10**-5, 10**-5, 10**4Bk = np.eye(x.size)for iter1 in range(maxiter):grad = num_grad(x, h)if np.linalg.norm(grad) < epsilon:return xdk = -np.dot((np.linalg.inv(Bk)), grad)ak = linesearch(x, dk)x = x + dk*akyk = num_grad(x, h) -gradsk = ak*dkif np.dot(yk.reshape(1, grad.shape[0]), sk) > 0:'''第⼀种分步计算实现t0 = np.dot(Bk, sk)t1 = np.dot(t0.reshape(sk.shape[0], 1), sk.reshape(1, sk.shape[0]))temp0 = np.dot(t1, Bk)temp1 = np.dot(np.dot(sk.reshape(1, sk.shape[0]), Bk), sk)tmp0 = np.dot(yk.reshape(yk.shape[0], 1), yk.reshape(1, yk.shape[0]))tmp1 = np.dot(yk.reshape(1, yk.shape[0]), sk)Bk = Bk - temp0 / temp1 + tmp0 / tmp1'''#第⼆种直接写公式实现Bk = Bk - np.dot(np.dot(np.dot(Bk, sk).reshape(sk.shape[0], 1), sk.reshape(1, sk.shape[0])), Bk)/np.dot(np.dot(sk.reshape(1, sk.shape[0]), Bk), sk) + np.dot(yk.reshape(yk.shape[0], 1), yk.reshape(1, yk.shape[0])) / np.dot(yk.reshape(1, yreturn x#x0 = np.array([0.999960983973235, 0.999921911551354]) #初始解x0 = np.array([0.7, 0.9]) #初始解x = steepest(x0) #调⽤最速下降法print("最速下降法最后的解向量：",x)print("最速下降法最后的解：",f(x))print('')x = newTonFuction(x0) #调⽤⽜顿法print("⽜顿法最后的解向量：", x)print("⽜顿法最后的解：", f(x))print('')x = BFGS(x0) #调⽤拟⽜顿法print("拟⽜顿法最后的解向量：", x)print("拟⽜顿法最后的解：", f(x))print('')结果如下拟⽜顿法感觉弄⿇烦了，暂时也没想法改，先就这样吧。

随着人工智能、模糊控制、模式识别、人工网络等新技术的应用和发展。

可以让它们与广义预测控制相结合,建立高精度、多模态的预测模型。

使广义预测控制在异常情况下可以稳定运行,推进广义预测控制的进一步发展。

2.2.1最速下降法最速下降法是无约束最优化中是比较有效的方法,它是以d}=一可(x})作为下降方向的算法。

其迭代格式为xx+i=xx一。

*Of (xk)上式中,一般通过精确线搜索准则求得步长因子。

*,当然也不排除可以利用非精确线搜索准则来求得步长因子。

*。

不管最速下降法采取何种线搜索准则,它均具有全局收敛性,但是这也不能直接就认为最速下降算法就是一个良好的优化算法。

在实际试验中,有很多优化问题利用最速下降法并不是下降的特快,反而下将的十分缓慢。

这是因为出现了锯齿现象:就是在计算过程中,最速下降法开始几步还是挺快的,但是当目标函数f (x)的等高线接近于一个球的时候,就出现了类似锯齿现象,前进十分缓慢,降低了算法的效能。

2.2.12.2.2牛顿法牛顿法也是无约束最优化问题中的一种经典算法,它是利用目标函数.f (x)的二次泰勒展开式,并将二次泰勒展开式进行极小化。

其迭代格式为x}+}=xA十d}(2-5)其中步长因子。

、=l} d、为02f (x} )d + Of (xA ) = 0的解。

当目标函数f(x)是正定二次函数的时候,牛顿法可以一步达到最优解;当目标函数f (x)是非二次函数的时候,牛顿法经过有限次迭代之后就不能确保求得目标函数f (x)的最优解。

我们知道目标函数f (x)在极小点附近是很接近于二次函数的,所以,假如初始点非常靠近无约束最优化问题((1-1)的最优解x的时候,并且}Z.f (x.)正定的时候,那么牛顿法就会有很快的收敛速度,而由此算法产生的点列也具有了超线性收敛速度,同时还在一定条件下具有二次收敛性;假如初始点与无约束最优化问题(1-1)的最优解x’相距比较远的时候,这时的}Z.}(x})就不一定是正定的了,也就存在了一个问题,那就是此时的牛顿方向就不一定是下降方向,有可能是上升方向,此时由此算法产生的点列可能也就不收敛于无约束最优化问题((1-1)的最优解了。