无约束优化——最速下降算法及Matlab实现

本文主要带大家回顾了一下梯度概念，并讲解了最速下降法的有趣点，最后教如何利用matlab 构建代码模块完成最速下降法的代码编写。本文与其他文章的对比的创新点就是提供了一个可视化的动态寻优图，还有就是本文的代码是经过多次测试是非常适用二元n维的函数。...

文章共57,400字 · 阅读需要大约192分钟

一键AI生成摘要，助你高效阅读

问答

cugautozp

15833人浏览 · 2022-05-22 20:25:00

cugautozp · 2022-05-22 20:25:00 发布

0 引言

最近太忙了，更新的有点慢了。本篇就更新一个通常用的较多的算法——最速下降法，其实数值分析中会讲到这个方法。此篇先回顾梯度概念，再讲述最速下降法，最后通过例子理解算法以及实现代码编写。

希望大家看完这篇文章之后能有所收获！！！！
下一篇文章：Matlab 自编雅可比矩阵 (jacobi) 函数与官方的Jacobian matrix（雅可比矩阵）函数对比及创新

这次需要提前必备一些 matlab 函数的使用

matlab 版本：2020a ，2022a (2022.6.5本人Matlab版本更新到了2022a)
matlab2022b版本中suptitle不能使用直接把suptitle改为sgtitle

句柄函数（function_handle）的使用¹——f=@(变量名) 函数表达式
matlabFunction 的使用²
范数函数 norm
绘制三维：mesh，等高线：contour，网格生成：meshgrid
获取每一帧图像：getframe，更新：drawnow，暂停延迟：pause （这些可以作为了解）

1 数学知识

了解这个最速下降法，需要一些数学知识的，所以先来回顾一些数学知识。

必看：如果数学知识没有掌握好的同学，我建议还是老老实实的看完这第一节内所有内容，否则后面有些关于公式推导会看的云里雾里飘。

1.1 梯度（gradient) 概念

这个概念参考中³ 有详细的解释说明，这里讲述一下其中梯度的概念

梯度：是一个向量，表示某一函数在该点处的方向导数沿着该方向取得最大值。（1）有大小：变化率最大（为该梯度的模）；（2）有方向：函数在该点处沿着该方向（此梯度的方向）变化最快⁴。

定义： 设 $z = f (x, y)$ 在点 $P_0(x_0,y_0)$ 处存在偏导数 $f'_x(x_0,y_0)$ 和 $f'_y(x_0,y_0)$ , 则称向量 ${f'_x(x_0,y_0) ,f'_y(x_0,y_0)\}$ 为 $f (x, y)$ 在 $P_0(x_0,y_0)$ 的梯度，记作 $\left.\nabla f\right|_{P_{0}},\left.\nabla z\right|_{P_{0}},\left.\operatorname{gradf}\right|_{P_{0}} \text { 或 }\left.\operatorname{gradz}\right|_{P_{0}}$
$\therefore \left.\nabla f\right|_{P_{0}}=\left.\operatorname{gradf}\right|_{P_{0}}=\{f'_x(x_0,y_0) ,f'_y(x_0,y_0)\}$

其中 $\nabla$ (Nabla)算子： $\nabla = \frac{\partial}{\partial x}i + \frac{\partial}{\partial y}j$

1.2 梯度函数

定义： 若 $f (x, y)$ 在 $D$ 内处处存在偏导数，则称 $\nabla f= \{f'_x(x,y) ,f'_y(x,y)\}$ 为 $f (x, y)$ 在D内的梯度函数

（1）梯度的大小： $|\nabla f|=\sqrt{\left[f_{x}^{\prime}(x, y)\right]^{2}+\left[f_{y}^{\prime}(x, y)\right]^{2}}$
（2）梯度的方向：

设 $v=\{v_1,v_2\} (|v|=1)$ 是任一给定方向，则对 $\nabla f$ 与 $v$ 的夹角 $\theta$ 有： $\begin{aligned} \left.\frac{\partial f}{\partial v}\right|_{P_{0}} &=f_{x}^{\prime}\left(x_{0}, y_{0}\right) v_{1}+f_{y}^{\prime}\left(x_{0}, y_{0}\right) v_{2} \\ &=\left\{f_{x}^{\prime}\left(x_{0}, y_{0}\right), f_{y}^{\prime}\left(x_{0}, y_{0}\right)\right\} \bullet\left\{v_{1}, v_{2}\right\} \\ &=\left.\nabla f\right|_{P_{0}} \bullet v=|\nabla f|_{P_{0}}|\cdot |v| \cos \theta \end{aligned}$
由于 $∣ v ∣ = 1$ 那么上式为： $\frac{\partial f}{\partial v}|_{P_{0}} =|\nabla f|_{P_{0}}|\cos \theta$

补充知识(高中内容)：已知两个向量 $\overrightarrow{v_1} =(a,b),\overrightarrow{v_2} =(c,d)$ ，那么两个向量内积： $v_1v_2=ac+bd=|\overrightarrow{v_1}||\overrightarrow{v_2}|\cos\theta$ ，其中： $\cos\theta=\frac{ac+bd}{|\overrightarrow{v_1}||\overrightarrow{v_2}|}$

$\therefore$ 当该点 $P_0$ 的方向 $v$ 与梯度方向 $\nabla f|_{P_{0}}$ 夹角为 $0$ 时取得最大，这就证明了沿着梯度方向函数值变化的最快，这公式也证明了⁵ ，梯度的方向是函数在指定点上升最快的方向，那么梯度的反方向自然是下降最快的方向了⁶：

如果点 $P_0$ 继续沿着梯度方向走，函数值增大
当 $\theta = 180°$ 即沿着梯度反方向走，函数值减小
垂直于梯度方向，函数值不变，即为：等高线

以上三个总结之后会用到，接下来开启最速下降法的算法讲解。

1.3 Hessian 矩阵

梯度的作用就是给了表示它变化的向量场。所以我觉得梯度跟导数更像一家人。继续来看二次近似： $Q(x)=f\left(x_{0}\right)+f^{\prime}\left(x_{0}\right)\left(x-x_{0}\right)+\frac{1}{2} f^{\prime \prime}\left(x_{0}\right)\left(x-x_{0}\right)^{2}$ 如果我们想把它推广到多元函数，那么我们会想应该也是类似的状况，除了前面的这一部分 $L(x,y)=a+b(x_0,y_0)+c(x_0,y_0)$ ，我们会想要再加上二次的部分： $d(x-x_0)^2+e(x-x_0)(y-y_0)+f(y-y_0)^2$
$\begin{array}{c} Q(x, y)=f\left(x_{0}, y_{0}\right)+f_{x}\left(x_{0}, y_{0}\right)\left(x-x_{0}\right)+f_{y}\left(x_{0}, y_{0}\right)\left(y-y_{0}\right) \\ +\frac{1}{2} f_{x x}\left(x_{0}, y_{0}\right)\left(x-x_{0}\right)^{2}+f_{x y}\left(x-x_{0}\right)\left(y-y_{0}\right)+\frac{1}{2} f_{y y}\left(x_{0}, y_{0}\right)\left(y-y_{0}\right)^{2} \end{array}$
通过类比或者计算，我们都可以得到以上的二次近似的结果。然后我们可以继续引出二次近似的向量形式： $Q(\mathbf{x})=f\left(\mathbf{x}_{0}\right)+\nabla f\left(\mathbf{x}_{\mathbf{0}}\right) \cdot\left(\mathbf{x}-\mathbf{x}_{\mathbf{0}}\right)+\frac{1}{2}\left(\mathbf{x}-\mathbf{x}_{0}\right)^{T} H\left(\mathbf{x}_{0}\right)\left(\mathbf{x}-\mathbf{x}_{0}\right)$
其中 $H$ 是 Hessian 矩阵： $\nabla (\nabla f)=H=\left[\begin{array}{cc} \frac{\partial f}{\partial x^{2}} & \frac{\partial f}{\partial x \partial y} \\ \frac{\partial f}{\partial y \partial x} & \frac{\partial f}{\partial y^{2}} \end{array}\right]$
当然我们也可以推广到更高的维度： $\mathbf{H}=\left[\begin{array}{cccc} \frac{\partial^{2} f}{\partial x_{1}^{2}} & \frac{\partial^{2} f}{\partial x_{1} \partial x_{2}} & \cdots & \frac{\partial^{2} f}{\partial x_{1} \partial x_{n}} \\ \frac{\partial^{2} f}{\partial x_{2} \partial x_{1}} & \frac{\partial^{2} f}{\partial x_{2}^{2}} & \cdots & \frac{\partial^{2} f}{\partial x_{2} \partial x_{n}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^{2} f}{\partial x_{n} \partial x_{1}} & \frac{\partial^{2} f}{\partial x_{n} \partial x_{2}} & \cdots & \frac{\partial^{2} f}{\partial x_{n}^{2}} \end{array}\right]$ Hessian 矩阵必定是一个对称矩阵

2 最速下降法

通过上面梯度知识的回顾，再来看这个算法就不会觉的很茫然。最速梯度下降法：解决的问题是无约束优化问题，而所谓的无约束优化问题就是对目标函数的求解，没有任何的约束限制的优化问题，比如求下方最小值： $min f (x)$

其中的函数 $f:R^n \to R$ .

求解这类的问题可以分为两大类：最优条件法和迭代法

2.1 概念⁷

最速下降法（Steepest descent）是梯度下降法的一种更具体实现形式，其理念为在每次迭代中选择合适的步长 $\alpha_k$ ，使得目标函数值能够得到最大程度的减少。

每一次迭代，沿梯度的反方向，我们总可以找到一个 $x^{(k+1)} = x^k - \alpha_k \cdot \nabla f(x^{(k)})$ ,使得在这个方向 $f(x^{(k+1)})$ 取最小值，即 $\alpha_{(k)}=\operatorname{argmin}f\left(x^{k}-\alpha \cdot \nabla f\left(x^{(k)}\right)\right)$

可以发下一个很有趣的事情就是：最速下降法每次更新的轨迹都和上一次垂直（证明过程在2.3中会讲解）

2.2 算法步骤⁷

下面先给出最速下降法的计算步骤：

第1步：选取初始点 $x^0$ ，给定终止误差 $\varepsilon >0$ ，令 $k = 0$ ；

第2步：计算 $\nabla f(x^k)$ ，若 $||\nabla f(x^k)|| \le \varepsilon$ ，停止迭代，输出 $x^k$ 。否则进行第三步；

第3步：取 $p^k = -\nabla f(x^k)$ ；

第4步：进行一维搜索，求 $\lambda_k$ ，使得 $f\left(x^{k}+\lambda_{k} p^{k}\right)=\min _{\lambda \geq 0} f\left(x^{k}+\lambda p^{k}\right)$ 令 $x^{k+1} = x^k + \lambda_k p^k，k=k+1，\text{转第 2 步}$

2.3 原理详解⁸

设 $x=x^{(k)}+\lambda p^{(k)}$ 在 $x^{(k)}$ 点处作一阶泰勒展开 $f(x)=f(x^{(k)}+\lambda p^{(k)}) = f(x^{(k)})+\lambda \nabla f(x^{(k)})^{T} p^{(k)}+O(||\lambda p^{(k)}||)$ 其中 $O(||\lambda p^{(k)}||)=O(||\lambda||)$ 是比 $\lambda$ 高阶的无穷小量，而 $p^{(k)}||=1$ .

因此 $f(x^{(k+1)})-f(x^{(k)}) \approx \lambda \nabla f(x^{(k)})^{T} p^{(k)}$
$\nabla f(x^{(k)})^{T} p^{(k)}=||\nabla f(x^{(k)})||\cdot||p^{(k)}|| \cos \theta=||\nabla f(x^{(k)})||\cdot\cos \theta$

搜索方向 $p^{(k)}$ ：在第1小节中讲述了取 $\theta = 180°$ 梯度负方向原因。 $\therefore$ $p^{(k)}=-\nabla f(x^{(k)})$

求梯度向量模的值： $||\nabla f(x^{(k)})^{T}||$ .

若 $||\nabla f(x^{(k)})^{T}|| < \varepsilon$ ，停止计算，输出 $x^{(k)}$ 作为极小点近似值。
若 $||\nabla f(x^{(k)})^{T}|| > \varepsilon$ ，转下一步

很多讲解都没有对这一个地方进行详细的说明，其实一开始本人也是不太理解为什么只需要证明了梯度向量的模长 $||\nabla f(x^{(k)})^{T}|| < \varepsilon$ 就可以了。

证明：
其实可以结合公式( $p^{(k)}||=1$ )： $f(x^{(k+1)})-f(x^{(k)}) \approx \lambda \nabla f(x^{(k)})^{T} p^{(k)}=||\nabla f(x^{(k)})||\cdot\cos \theta$ ，其中的 $f(x^{(k+1)})$ 表示的是下一个寻优值， $f(x^{(k)})$ 表示当前的迭代值：
左边式子的含义：就是当前值与未来值之差越来越小，那就说明变化的速度越来越慢了，速度慢下来了，那不就意味着离这个最优的极小点越来越近了。
右边式子含义：左边两个式子差值最终等于一个梯度的大小 $||\nabla f(x^{(k)})||$ × $\cos \theta$ ，大家都知道 $\cos \theta \in [0,1]$ ， $\therefore$ 当 $||\nabla f(x^{(k)})^{T}|| < \varepsilon$ 时是不是说明找的这个值是不是可以近似看作一个极小值，且 $x^{(k)}$ 为极小点。

最优步长 $\lambda_k$ ：

设 $f (x)$ 具有二阶连续偏导数，将它在 $x^{(k)}$ 点作二阶泰勒展开： $\begin{aligned} f\left(x^{(k)}-\lambda \nabla f\left(x^{(k)}\right)\right) &=f\left(x^{(k)}\right)+\nabla f\left(x^{(k)}\right)^{T}\left(-\lambda \nabla f\left(x^{(k)}\right)\right) \\ &+\frac{1}{2}\left(-\lambda \nabla f\left(x^{(k)}\right)\right)^{T} H\left(x^{(k)}\right)\left(-\lambda \nabla f\left(x^{(k)}\right)\right) \\ &+O\left(\| \lambda \nabla f\left(x^{(k)} \|^{2}\right)\right. \end{aligned}$
将上式主要部分记作 $H(\lambda)$ $H(\lambda)= f\left(x^{(k)}\right)+\nabla f\left(x^{(k)}\right)^{T}\left(-\lambda \nabla f\left(x^{(k)}\right)\right)+\frac{1}{2}\left(-\lambda \nabla f\left(x^{(k)}\right)\right)^{T} H\left(x^{(k)}\right)\left(-\lambda \nabla f\left(x^{(k)}\right)\right)$
函数 $H(\lambda)$ 的唯一的极值点（令 $H'(\lambda)=0$ ，即可求出 $\lambda$ ）为： $\lambda_{k}=\frac{\nabla f\left(x^{(k)}\right)^{T} \nabla f\left(x^{(k)}\right)}{\nabla f\left(x^{(k)}\right)^{T} H\left(x^{(k)}\right) \nabla f\left(x^{(k)}\right)}$
其中 $H(x^{(k)})$ 在第1节说过是 Hessian 矩阵

这里证明前面说的垂直关系，证明如下：

将点 $x^{(k+1)}=x^{(k)}+\lambda p^{k}$ 代入函数 $f (x)$ 中，再令 $f^{'} (x) = 0$ 可以得到： $f'(x^{(k+1)})=f'(x^{(k)}+\lambda p^{k})=\nabla f(x^{(k)}+\lambda p^{k})^{T}p^{(k)}=0$
将 $p^{(k)}=-\nabla f(x^{(k)})$ 代入上式： $-\nabla f(x^{(k)}+\lambda p^{k})^{T}\nabla f(x^{(k)})=-\nabla f(x^{(k+1)})^{T}\nabla f(x^{(k)})=0$
这不就证明出来了，前一个点与后一个点的梯度是垂直关系

下面讲下这个步长 $\lambda$ 取大取小的影响，以一元二次函数为例，如图 1 所示（这个图本人用matlab画出来的，代码会放大最后）：

$\text{图 1}$
可以从这个图中看到最优步长为 $\lambda_{op}$ ，当步长 $\lambda \in(0,\lambda_{op}) \cup (\lambda_{op},2\lambda_{op})$ 迭代次数很多收敛速度很慢。很有意思的是当 $\lambda =2\lambda _{op}$ 时不迭代了直接在等高线上来回迭代。当 $\lambda \in(2\lambda_{op},\infty)$ 时反向迭代了，远离了极小值点。

2.4 缺点⁷

某点的负梯度方向，通常只是在该点附近才具有这种最速下降的性质。

最速下降法利用目标函数一阶梯度进行下降求解，易产生锯齿现象（如下图），在开始几步，目标函数下降较快；但在接近极小点时，收敛速度长久不理想了。特别适当目标函数的等值线为比较扁平的椭圆时，收敛就更慢了。
在这里插入图片描述
因此，在实用中常用最速下降法和其他方法联合应用，在前期使用最速下降法，而在接近极小值点时，可改用收敛较快的其他方法。

3 示例说明

这个例子就用这篇文章的示例⁹ 虽然用了他的例子本人代码和他代码完全不一样。如果大家感兴趣也可以去看这个人的代码。好了废话不多说上题目： $min f(x) = x_1^{2}+2x_2^2-2x_1x_2-2x_2$ 其中 $x=(x_1,x_2)^T$ ， $x^{0}=(0,0)^T$ 。

(1) 求梯度函数 $\nabla f(x)=\begin{pmatrix} 2x_1-2x_2\\ 4x_2-2x_1-2\end{pmatrix}$ Hessian 矩阵 $H(x)=\nabla(\nabla f(x))=\begin{pmatrix} 2&-2\\ -2&4\end{pmatrix}$

(2) 将 $x^{(0)}=(0,0)^T$ 代入上式 (1) 中，得到 $\nabla f(x^{(0)})=\begin{pmatrix} 0\\ -2\end{pmatrix}$ ，即 $p^{(0)}=-\nabla f(x^{(0)})=\begin{pmatrix} 0\\ 2\end{pmatrix}$

(3) $H(x^{(0)})=\begin{pmatrix} 2&-2\\ -2&4\end{pmatrix}$ 和 $\nabla f(x^{(0)})$ 代入 $\lambda_{k}=\frac{\nabla f\left(x^{(k)}\right)^{T} \nabla f\left(x^{(k)}\right)}{\nabla f\left(x^{(k)}\right)^{T} H\left(x^{(k)}\right) \nabla f\left(x^{(k)}\right)}$
得到 $\lambda=\frac{\begin{pmatrix} 0& -2\end{pmatrix} \begin{pmatrix} 0\\ -2\end{pmatrix}}{\begin{pmatrix} 0& -2\end{pmatrix}\begin{pmatrix} 2&-2\\ -2&4\end{pmatrix} \begin{pmatrix} 0\\ -2\end{pmatrix}}=\frac{1}{4}$
(4) $x^{(1)}=x^{(0)}+\lambda p^{(0)}=\begin{pmatrix} 0\\ \frac{1}{2}\end{pmatrix}$

同理转回（2）可以一直迭代回去一直到满足条件为止，得到最优解为 $x^*=(1,1)^T，y^*=-1$

4 代码实现

根据上面的方法步骤来看，其实编写代码的时候只需要几个模块就行，（1）求 $\nabla f(x)$ 的函数，（2）求 Hessian 矩阵 $H (x)$ 的模块；（3）将前面两个函数放在一起构造一个求最速下法的函数即可；（4）为了可视化，本人又增加了一个生成GIF的动态最速下降法的画图函数。

声明：

构造了：df=nabla_f(fun,x)，求 $\nabla f(x)$ 的函数
构造了：H=Hesse(df,x,x0,n)，求 Hessian 矩阵 $H (x)$ 的模块
构造了：[x0,i] = GD(f,x,x0,epsilon)，为最速下降法函数
构造了：GDplot(f,x0,i,x1,x2,GIFname)，最速下降画图函数
本人代码现在适用于 二元n次的函数 进行梯度下降寻极小值点。（(2022.5.30)现已经解决了限制，可以实现n元n次函数，GDplot()画图函数只能绘制三维图）

注意：接下来的几个模快代码里有会用到在开头就提醒大家得必备一些函数的应用。

先放结果：
在这里插入图片描述
生成的GIF图：

4.1 $\nabla f(x)$ 的函数模块代码实现

这个函数的输入有两个：（1）fun：表示函数表达式，matlab中一般用句柄函数形式表达。句柄函数写法在文章的开头就写过，或者不了解的同学可以去看参考文献的第一个。（2）x：表示变量字符串，这里一般是 $x_1,x_2$ ，下面也给出了例子。
输出：df：是一个句柄函数

代码如下：

function df=nabla_f(fun,x)
%  ∇f 求梯度
    df=[];
    x=str2sym(x);
    for i=1:length(x)
        df1 = diff(fun,x(i));
        df=[df;df1];
    end
    df= matlabFunction(df);
end

例如输入以下代码运行后可以得到一个df的句柄函数：

f=@(x1,x2) x1.^2+2*x2.^2-2*x1.*x2-2*x2;
x ='[x1,x2]';
df=nabla_f(f,x)

output：
在这里插入图片描述

4.2 Hessian 矩阵 $H (x)$ 的模块代码实现

这个函数的输入有4个：（1）df：表示为 $\nabla f$ ，即为函数nabla_f 的结果。（2）x：表示变量字符串，这里一般是 $x_1,x_2$ 。（3）x0：表示初始的 $x^{(0)}$ 。（4）n：表示当前迭代次数。
输出：df：是一个句柄函数

代码如下：

function H = Hesse(df,x,x0,n)
%  Hesse矩阵 H(x)=∇(∇f) 
%  df 为 ∇f，即为函数nabla_f的结果 ∇f
    H=[];
    x=str2sym(x);
    for i=1:length(x)
        df1 = diff(df,x(i));
        H=[H,df1];
    end
%     H = matlabFunction(H);
    s=char(H);
    if find(s=='x')
        H = matlabFunction(H);
        H = H(x0{n}(1),x0{n}(2));
    else
        H = double(H);
    end
end

4.3 最速下降法函数模块代码实现

输入：（1）f：为上面定义的句柄函数例如：f= @(x1,x2) 2*x1.^2+2*x2.^2+2*x1.*x2+x1-x2。（2） x：为表示变量字符串，这里一般是 $x_1,x_2$ .（3）x0：表示起始点 $x^{(0)}$ 。（4）epsilon 为 $ε$ 精度

代码如下：

function [x0,i] = GD(f,x,x0,epsilon)
% 梯度下降法 也可以称：最速下降法
% input：f 为上面定义的句柄函数 例如：f= @(x1,x2) 2*x1.^2+2*x2.^2+2*x1.*x2+x1-x2;
%        x 为表示变量字符串，这里一般是 x1,x2
%        x0 表示起始点
%        epsilon 为ε精度
% output：x0 为GD函数出来的从起始点x0寻优到极小值点的所有点集合
%        i 为GD函数寻优过程得到点的个数即：i==length(x0)
%        可以打印出最小值 min{f}，以及极小值值点 x
%

    i=1;
    df = nabla_f(f,x);
    H = Hesse(df,x,x0,i);
    dfx=df(x0{1}(1),x0{1}(2));
   
    er=norm(dfx);
    while er > epsilon                
            p = -dfx;
            lambda = dfx'*dfx/(dfx'*H*dfx);
            i=i+1;
            x0{i} = x0{i-1}+lambda*p;
            dfx = df(x0{i}(1),x0{i}(2));
            H = Hesse(df,x,x0,i);
            er = norm(dfx);     
    end
    fmin = f(x0{i}(1),x0{i}(2));
    disp('极小值点：');
    disp(['x1 = ',num2str(x0{i}(1))]);
    disp(['x2 = ',num2str(x0{i}(2))]);
    fprintf('\nf最小值：\n min{f}=%f\n',fmin);
end

*4.4 最速下降法画图函数模块代码实现

这个函数代码感兴趣的同学可以看看，个人只把代码放在这里就不加解释。主要是写了这么长太累偷个懒不过分吧！

创新点：1. 可以根据 f的表达式自动生成图的标题，意思就是不需要自己手动打图标题，输入任意的 $f$ 可以实现任意 $f$ 的图标题显示，这都得益于写了一个 f2t() 函数，该函数实现了将句柄函数格式 (function_handle) 转化为 $T e x$ 格式的字符串。

更新2022-11-6已经将suptitle改为了sgtitle，并上传了f2s函数

代码如下：

function GDplot(f,x0,i,x1,x2,GIFname)
% input：f 为上面定义的句柄函数 例如：f= @(x1,x2) 2*x1.^2+2*x2.^2+2*x1.*x2+x1-x2;
%        x0 为GD函数出来的从起始点x0寻优到极小值点的所有点集合
%        i 为GD函数寻优过程得到点的个数即：i==length(x0)
%        x1 为f函数中第一个变量的取值范围
%        x2 为f函数中第二个变量的取值范围
% 
% output：生成一个gif

        [x1,x2]=meshgrid(x1,x2);
        z=f(x1,x2);
        figure('color','w')
        sgtitle(['\it f=',f2s(f)])
        subplot(211) 
        mesh(x1,x2,z)
        axis off
        view([-35,45])
        hold on
        subplot(212)
        contour(x1,x2,z,20)
        zlim([0,0.5])
        set(gca,'ZTick',[],'zcolor','w')
        axis off
        view([-35,45])
        hold on

        pic_num = 1;

        for j =1:i-1    
            a=[x0{j}(1),x0{j}(2),f(x0{j}(1),x0{j}(2))];
            b=[x0{j+1}(1),x0{j+1}(2),f(x0{j+1}(1),x0{j+1}(2))];
            c=[a',b'];
            a1=[x0{j}(1),x0{j}(2)];
            b1=[x0{j+1}(1),x0{j+1}(2)];
            c1=[a1',b1'];
    
            subplot(211)
            plot3(x0{j}(1),x0{j}(2),f(x0{j}(1),x0{j}(2)),'r.','MarkerSize',10)
            subplot(212)
            plot(x0{j}(1),x0{j}(2),'r.','MarkerSize',10)
            drawnow
            F(j)=getframe(gcf);
            pause(0.5)

            subplot(211)
            plot3(c(1,:),c(2,:),c(3,:),'r--')    
            subplot(212)
            plot(c1(1,:),c1(2,:),'r--')
            drawnow 
            F(2*j)=getframe(gcf);
            pause(0.5)
            
            % 绘制并保存gif
            I=frame2im(F(j));
            [I,map]=rgb2ind(I,256);
            I1=frame2im(F(2*j));
            [I1,map1]=rgb2ind(I1,256);
            if pic_num == 1
                imwrite(I,map, GIFname ,'gif', 'Loopcount',inf,'DelayTime',0.5);
            else
                imwrite(I,map, GIFname ,'gif','WriteMode','append','DelayTime',0.5);
                imwrite(I1,map1, GIFname ,'gif','WriteMode','append','DelayTime',0.5);
            end
            pic_num = pic_num + 1;
        end
        subplot(211)
        plot3(x0{i}(1),x0{i}(2),f(x0{i}(1),x0{i}(2)),'r.','MarkerSize',9)
        subplot(212)
        plot(x0{i}(1),x0{i}(2),'r.','MarkerSize',9)
        F(2*i-1)=getframe(gcf);
        I=frame2im(F(2*i-1));
        [I,map]=rgb2ind(I,256);
        imwrite(I,map, GIFname ,'gif','WriteMode','append','DelayTime',0.5);
end

function s = f2s(fun)
% 句柄函数的转换为字符串 
% 主要是用来画图title用的，可自动将句柄函数转为字符串函数
% 在绘制画图时可以自动生成函数表达式 title（并且以Latex形式显示出来），避免手动敲击函数公式
    s = func2str(fun);
    s = char(s);
    c = strfind(s,')');
    s(1:c(1))=[];
    c1 = strfind(s,'.');
    s(c1)=[];
    c2 = strfind(s,'*');
    s(c2)=[];
    c3 = strfind(s,'x');
    for  i = 1:length(c3)
       s = insertAfter(s,c3(i)+i-1,'_'); 
    end   
end

4.5 主函数

代码如下：

clc
clear all
close all
% f= @(x1,x2) 2*x1.^2+2*x2.^2+2*x1.*x2+x1-x2;
f= @(x1,x2) x1.^2+2*x2.^2-2*x1.*x2-2*x2;      % 句柄函数表达式
% f= @(x1,x2) (x1-1).^2+(x2-1).^2;
% f=@(x1,x2) x1.^4+3*x1.^2*x2-x2.^4;
x0{1}=[0;0];
x ='[x1,x2]';      % 函数变量字符串
epsilon=0.001;     % ε为误差精度
[x0,i] = GD(f,x,x0,epsilon);
x1=0:0.01:2;
x2=x1;
GIFname = 'f1.gif';
GDplot(f,x0,i,x1,x2,GIFname)

4.6 上面图1的代码

图1代码直接下载文件链接：点击此处下载代码文件

4.7 完整代码文件

如果是小白的同学可以下载这个完整代码（包含本文章的代码+图1的代码）：
代码已经放入个人GitHub中：https://github.com/cug-auto-zp/CSDN/tree/main/gradient
如果进不了GitHub的读者，这里还提供了资源下载链接：资源链接

5 总结

写这篇文章还是要感谢这些博主写的文章个人感觉看了他们的文章收获非常多。给大家代码一方面是方便大家学习，另一方面是我也是通过大佬们的文章和公开的代码发现matlab有很个人不知道的函数，边看边学代码，有很多人问我是这么学的matlab，其实也不用刻意去学习这个东西，我都是在边看边积累的过程不断地扩充了自己的matlab的知识。所以个人坚持将这种开源的代码的思想一直进行下去。