无约束优化方法

一、最速下降法(梯度下降法)

基本思想:

  1. 初始点 x0Rnx^{0} \in \mathbb{R}^{n}0ϵ10 \leq \epsilon \ll 1k=0k=0
  2. 计算 f(xk)\nabla f\left(\boldsymbol{x}^{k}\right) ,如果 f(xk)2ϵ\left\|\nabla f\left(\boldsymbol{x}^{k}\right)\right\|_{2} \leq \epsilon ,停止迭代:否则,令 dk=f(xk)\boldsymbol{d}^{k}=-\nabla f\left(\boldsymbol{x}^{k}\right)
  3. 沿 dkd^{k} 方向确定步长 αk\alpha_{k}
  4. xk+1=xk+αkdk\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}+\alpha_{k} \boldsymbol{d}^{k}k=k+1k=k+1,转 Step 2。

精确一维搜索最速下降法:
若最速下降法 Step 3 中确定的步长 αk\alpha_{k} 满足:αk=argminα>0f(xk+αdk)\alpha_{k}=\arg \min _{\alpha>0} f\left(\boldsymbol{x}^{k}+\alpha \boldsymbol{d}^{k}\right),即 f(xk+αkdk)=minα>0f(xk+αdk)f\left(\boldsymbol{x}^{k}+\alpha_{k} \boldsymbol{d}^{k}\right)=\min _{\alpha>0} f\left(\boldsymbol{x}^{k}+\alpha \boldsymbol{d}^{k}\right) 。也即,若 αk\alpha_{k} 经精确一维搜索得到,则相应的最速下降算法称为精确一维(线)搜索最速下降法,且称 αk\alpha_{k} 为精确步长因子。

最速下降法特点:
dk\boldsymbol{d}^{k}dk+1\boldsymbol{d}^{k+1} 正交,即相邻两次迭代中搜索方向是正交的。最速下降法向极小点逼近是曲折前进的,这种现象称为锯齿现象。该现象会使得迭代点逼近极小点的过程是"之"字形。而在实际优化过程中,对于任意初始点,最速下降法都可以很快到达极小点附近,但是越靠近极小点步长越小,导致最速下降法的收敛速度很慢。

证明:由于 αk=argminα>0f(xk+αdk)\alpha_{k}=\arg \min _{\alpha>0} f\left(\boldsymbol{x}^{k}+\alpha \boldsymbol{d}^{k}\right) ,则 df(xk+αdk)dαα=αk=0\left.\frac{d f\left(\boldsymbol{x}^{k}+\alpha \boldsymbol{d}^{k}\right)}{d \alpha}\right|_{\alpha=\alpha_{k}}=0 。即 (dk)Tf(xk+αkdk)=(dk)Tf(xk+1)=0\left(\boldsymbol{d}^{k}\right)^{\mathrm{T}} \nabla f\left(\boldsymbol{x}^{k}+\alpha_{k} \boldsymbol{d}^{k}\right)=\left(\boldsymbol{d}^{k}\right)^{\mathrm{T}} \nabla f\left(\boldsymbol{x}^{k+1}\right)=0
dk+1=f(xk+1)\boldsymbol{d}^{k+1}=-\nabla f\left(\boldsymbol{x}^{k+1}\right) ,故 (dk)Tdk+1=0\left(\boldsymbol{d}^{k}\right)^{\mathrm{T}} \boldsymbol{d}^{k+1}=0 ,即正交。

最速下降法向极小点逼近是曲折前进的,这种现象称为锯齿现象。该现象会使得迭代点逼近极小点的过程是"之"字形。而在实际优化过程中,对于任意初始点,最速下降法都可以很快到达极小点附近,但是越靠近极小点步长越小,导致最速下降法的收敛速度很慢。

收敛性定理:
f(x)f(\boldsymbol{x}) 连续可微,且水平集有界L={xRnf(x)f(x0)}L=\left\{\boldsymbol{x} \in \mathbb{R}^{n} \mid f(\boldsymbol{x}) \leq f\left(\boldsymbol{x}^{0}\right)\right\},则对于某个 kk ,由精确一维搜索最速下降法产生的序列满足:f(xk)2=0\left\|\nabla f\left(\boldsymbol{x}^{k}\right)\right\|_{2}=0,或limk+f(xk)2=0\lim _{k \rightarrow+\infty}\left\|\nabla f\left(\boldsymbol{x}^{k}\right)\right\|_{2}=0

优点:

  • 计算复杂度低,存储量小;
  • 具有全局收敛性,对初始点要求不严格。

缺点:
-收敛速度慢,尤其在极小点附近,该方法产生的点列逼近函数的极小点越来越慢。

导致缺点的原因:

  • dk=f(xk)\boldsymbol{d}^{k}=-\nabla f\left(\boldsymbol{x}^{k}\right) 仅反映 f(x)f(\boldsymbol{x})xk\boldsymbol{x}^{k} 处的局部性质;
  • dk\boldsymbol{d}^{k}dk+1\boldsymbol{d}^{k+1} 正交,即相邻两次迭代的搜索方向正交,致使优化"远快近慢"。

为了解决最速下降法中相邻两步搜索方向正交的所导致的收敛速度慢的问题,可使用以下改进方法:

  1. 选择不同初始点:有利于避免迭代过程中出现锯齿现象。例如对于 f(x)=x12+25x22f(\boldsymbol{x})=x_{1}^{2}+25 x_{2}^{2} ,分别取 x0=(2,2)T\boldsymbol{x}^{0}=(2,2)^{\mathrm{T}}x0=(100,0)T\boldsymbol{x}^{0}=(100,0)^{\mathrm{T}} 。虽然后一初始点较前一初始点离 xx^{*} 更远,但迭代中不出现锯齿现象。
  2. 采用不精确的一维搜索:可使相邻两个迭代点处的梯度不正交,从而改变收敛性。对于最速下降法,有时为了减少计算工作量,采用固定步长,称为固定步长最速下降法。
  3. 采用加速梯度法:负梯度方向和 pk=xkxk2\boldsymbol{p}^{k}=\boldsymbol{x}^{k}-\boldsymbol{x}^{k-2} 相结合。即,由于最速下降法在极小点附近呈"锯齿"状,因此第 kk 步的搜索方向取为 pk\boldsymbol{p}^{k} ,而第 k+1k+1 和第 k+2k+2 步继续使用负梯度方向。

二、共轭梯度法

1.问题背景和基本思想:

最速下降法计算复杂度低,但锯齿现象致使其收敛速度慢。因此,为避免锯齿现象出现,可利用当前迭代点 xk\boldsymbol{x}^{k} 处的负梯度方向 f(xk)-\nabla f\left(\boldsymbol{x}^{k}\right) 与上一步的搜索方向 dk1\boldsymbol{d}^{k-1} ,通过寻找 f(xk)-\nabla f\left(\boldsymbol{x}^{k}\right)dk1\boldsymbol{d}^{k-1} 两者的适当线性组合,从而构造出搜索方向 dk\boldsymbol{d}^{k} ,而非每次迭代的搜索方向都为负梯度方向。

由 Taylor 公式知,函数在某点附近的性质与二次函数十分接近。因此,先针对正定二次函数建立简单、高效、收敛速度快的求解算法,然后再推广到一般函数上去。

2.共轭及其性质

d1,d2,,dn\boldsymbol{d}_{1}, \boldsymbol{d}_{2}, \cdots, \boldsymbol{d}_{n}Rn\mathbb{R}^{n} 中任意一组非零向量, QRn×n\boldsymbol{Q} \in \mathbb{R}^{n \times n} 是 n 阶对称正定矩阵。若

diTQdj=0(ij)\boldsymbol{d}_{i}^{\mathrm{T}} \boldsymbol{Q} \boldsymbol{d}_{j}=0(i \neq j)

则称 di\boldsymbol{d}_{i}dj\boldsymbol{d}_{j}Q\boldsymbol{Q} 共轭的,向量组 d1,d2,,dn\boldsymbol{d}_{1}, \boldsymbol{d}_{2}, \cdots, \boldsymbol{d}_{n}Q\boldsymbol{Q} 共轭向量组。
Q=I\boldsymbol{Q}=\boldsymbol{I}diTIdj=0(ij)\boldsymbol{d}_{i}^{\mathrm{T}} \boldsymbol{I} \boldsymbol{d}_{j}=0(i \neq j) ,则 d1,d2,,dn\boldsymbol{d}_{1}, \boldsymbol{d}_{2}, \cdots, \boldsymbol{d}_{n} 是正交的,因此共轭是正交的推广,共轭向量组是正交向量组的推广。

Q\boldsymbol{Q}nn 阶对称正定矩阵,非零向量组 d1,d2,,dn\boldsymbol{d}_{1}, \boldsymbol{d}_{2}, \cdots, \boldsymbol{d}_{n} 关于 Q\boldsymbol{Q} 共轭,则该向量组必线性无关。

Q\boldsymbol{Q}nn 阶对称正定矩阵,非零向量组 d1,d2,,dn\boldsymbol{d}_{1}, \boldsymbol{d}_{2}, \cdots, \boldsymbol{d}_{n} 关于 Q\boldsymbol{Q} 共轭,则它们构成 Rn\mathbb{R}^{n} 的一组基。

Q\boldsymbol{Q}nn 阶对称正定矩阵,非零向量组 d1,d2,,dn\boldsymbol{d}_{1}, \boldsymbol{d}_{2}, \cdots, \boldsymbol{d}_{n} 关于 Q\boldsymbol{Q} 共轭。若向量 v\boldsymbol{v}d1\boldsymbol{d}_{1}, d2\boldsymbol{d}_{2}, \cdots, dn\boldsymbol{d}_{n} 均正交,则 v=0\boldsymbol{v}=\mathbf{0}

3.为什么用共轭方向?

对于下列正定二次函数求极小问题:

\min f(\boldsymbol{x})=\frac{1}{2} \boldsymbol{x}^{\mathrm{T}} \boldsymbol{Q} \boldsymbol{x}+\boldsymbol{b}^{\mathrm{T}} \boldsymbol{x}+c$$以二元二次函数为例 ($n=2$) ,取任意初始点 $\boldsymbol{x}^{1}$ ,沿某个下降方向 $\boldsymbol{d}^{1}$ 作精确—维搜索得$x^{2}=x^{1}+\alpha_{1} d^{1}$。 由精确—维搜索可知: $\nabla f\left(\boldsymbol{x}^{2}\right)^{\mathrm{T}} \boldsymbol{d}^{1}=0$ 。若要提升迭代算法收敛速度,则不能取 $\boldsymbol{d}^{2}=-\nabla f\left(\boldsymbol{x}^{2}\right)$ 。否则引发锯齿现象,致使收敛速度变慢。 如果希望下次迭代就可以求得极小点 $\boldsymbol{x}^{*}$ ,则搜索方向 $\boldsymbol{d}^{2}$ 该如何选取? 如果能够选定这样的搜索方向,那么对于二元二次函数只需依次沿搜索方向 $d^{1}, d^{2}$ 各进行次精确一维搜索就可以求到极小点 $x^{*}$ ,即:$\boldsymbol{x}^{*}=\boldsymbol{x}^{2}+\alpha_{2} \boldsymbol{d}^{2}$。 注意到 $x^{*}$ 是 $f(x)$ 的极小点,故 $x^{*}$ 是 $f(x)$ 的驻点:$\nabla f\left(\boldsymbol{x}^{*}\right)=\boldsymbol{Q} \boldsymbol{x}^{*}+\boldsymbol{b}=\mathbf{0}$。 展开整理:$$\begin{array}{c l}\boldsymbol{0}&=\nabla f\left(\boldsymbol{x}^{*}\right)\\\boldsymbol{0}&=\boldsymbol{Q}\left(\boldsymbol{x}^{2}+\alpha_{2} \boldsymbol{d}^{2}\right)+\boldsymbol{b}\\\boldsymbol{0}&=\boldsymbol{Q} \boldsymbol{x}^{2}+\boldsymbol{b}+\alpha_{2} \boldsymbol{Q} \boldsymbol{d}^{2}\\\boldsymbol{0}&=\nabla f\left(\boldsymbol{x}^{2}\right)+\alpha_{2} \boldsymbol{Q} \boldsymbol{d}^{2}\end{array}$$等式两边同时左乘 $\left(d^{1}\right)^{\mathrm{T}}$ 可得:$$\begin{array}{c l}\boldsymbol{0}&=\left(\boldsymbol{d}^{1}\right)^{\mathrm{T}}\nabla f\left(\boldsymbol{x}^{2}\right)+\alpha_{2}\left(\boldsymbol{d}^{1}\right)^{\mathrm{T}} \boldsymbol{Q} \boldsymbol{d}^{2}\\\boldsymbol{0}&=0+\alpha_{2}\left(\boldsymbol{d}^{1}\right)^{\mathrm{T}} \boldsymbol{Q} \boldsymbol{d}^{2}\\\boldsymbol{0}&=\left(\boldsymbol{d}^{1}\right)^{\mathrm{T}} \boldsymbol{Q} \boldsymbol{d}^{2}\end{array}$$这说明,搜索方向 $d^{1}$ 和 $d^{2}$ 关于$\boldsymbol{Q}$ 共轭,其中 $d^{1}$ 是某个搜索方向。因此,用迭代法求解正定二次函数极小的问题:$\min f(\boldsymbol{x})=\frac{1}{2} \boldsymbol{x}^{\mathrm{T}} \boldsymbol{Q x}+\boldsymbol{b}^{\mathrm{T}} \boldsymbol{x}+c$ 可转化为构造或寻找 $n$ 个关于 $\boldsymbol{Q}$ 共轭的方向的问题。 基本定理:设 $\boldsymbol{Q}$ 为 $n$ 阶对称正定矩阵,向量组 $\boldsymbol{d}_{1}, \boldsymbol{d}_{2}, \cdots, \boldsymbol{d}_{n}$ 关于 $\boldsymbol{Q}$ 共轭,对正定二次函数$\min f(\boldsymbol{x})=\frac{1}{2} \boldsymbol{x}^{\mathrm{T}} \boldsymbol{Q x}+\boldsymbol{b}^{\mathrm{T}} \boldsymbol{x}+c$,给定任意初始点 $\boldsymbol{x}^{1}$ ,依次沿 $\boldsymbol{d}_{i}(i=1, \cdots, n)$ 进行精确线搜索,则至多经过 $n$ 次迭代即可求得 $f(\boldsymbol{x})$ 的极小值点。 **二次收敛性(二次终止性)**: 给定任意初始点,若某算法经有限次迭代即可达到正定二次函数的极小点,称该算法具有二次收敛性。 ##### 4.求解正定二次函数的共轭梯度法 (1)基于**待定系数法**构造共轭方向 基本思想:结合当前迭代点处的负梯度方向以及上一个迭代点处的下降方向来构造共轭方向。同时,迭代格式仍为 $x^{k+1}=x^{k}+\alpha_{k} d^{k}$ ,步长 $\alpha_{k}$ 由精确线搜索获得,且第一步迭代时取 $\boldsymbol{d}^{1}=-\nabla f\left(\boldsymbol{x}^{1}\right)$ 。 根据待定系数法,令新的搜索方向为 $\boldsymbol{d}^{k+1}=-\nabla f\left(\boldsymbol{x}^{k+1}\right)+\lambda_{k} \boldsymbol{d}^{k}$,其中 $\lambda_{k}$ 末知且为待定组合系数。要求 $\boldsymbol{d}^{k}$ 与 $\boldsymbol{d}^{k+1}$ 是 $\boldsymbol{Q}$ 共轭的,则左乘 $\left(\boldsymbol{d}^{k}\right)^{\mathrm{T}} \boldsymbol{Q}$ ,并令 $\left(\boldsymbol{d}^{k}\right)^{\mathrm{T}} \boldsymbol{Q} \boldsymbol{d}^{k+1}=0$ ,解之得: $$\lambda_{k}=\frac{\left(\boldsymbol{d}^{k}\right)^{\mathrm{T}} \boldsymbol{Q} \nabla f\left(\boldsymbol{x}^{k+1}\right)}{\left(\boldsymbol{d}^{k}\right)^{\mathrm{T}} \boldsymbol{Q} \boldsymbol{d}^{k}}

对于上述方式产生的下降方向 dk\boldsymbol{d}^{k} ,可以证明:(di)TQdk+1=0,i=1,2,,k\left(\boldsymbol{d}^{i}\right)^{\mathrm{T}} \boldsymbol{Q} \boldsymbol{d}^{k+1}=0, i=1,2, \cdots, k

也即,按上述方法得到的 nn 个方向 d1,d2,,dn\boldsymbol{d}^{1}, \boldsymbol{d}^{2}, \cdots, \boldsymbol{d}^{n} 关于 Q\boldsymbol{Q} 是两两共轭的。
此外,上述共轭梯度法迭代过程中,序列 {xk}\left\{x^{k}\right\} 所对应的梯度向量组 f(xk),k=1,2,,n\nabla f\left(\boldsymbol{x}^{k}\right), k=1,2, \cdots, n 是正交向量组。

(2)适用于求解正定二次函数的共轭梯度法算法

  1. 初始化 x1Rnx^{1} \in \mathbb{R}^{n}d1=f(x1)\boldsymbol{d}^{1}=-\nabla f\left(\boldsymbol{x}^{1}\right)0ϵ<10 \leq \epsilon<1k=1k=1
  2. 如果 f(xk)2ϵ\left\|\nabla f\left(\boldsymbol{x}^{k}\right)\right\|_{2} \leq \epsilon ,停止迭代;否则沿 dkd^{k} 精确线搜索求步长:$$\alpha_{k}=\arg \min _{\alpha>0} f\left(\boldsymbol{x}^{k}+\alpha \boldsymbol{d}{k}\right)=-\frac{\left(\boldsymbol{d}{k}\right)^{\mathrm{T}} \nabla f\left(\boldsymbol{x}{k}\right)}{\left(\boldsymbol{d}{k}\right)^{\mathrm{T}} \boldsymbol{Q} \boldsymbol{d}^{k}}

    令 $\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}+\alpha_{k} \boldsymbol{d}^{k}$ ,并计算: $\boldsymbol{d}^{k+1}=-\nabla f\left(\boldsymbol{x}^{k+1}\right)+\lambda_{k} \boldsymbol{d}^{k}$ ,其中 $\lambda_{k}=\frac{\left(\boldsymbol{d}^{k}\right)^{\mathrm{T}} \boldsymbol{Q} \nabla f\left(\boldsymbol{x}^{k+1}\right)}{\left(\boldsymbol{d}^{k}\right)^{\mathrm{T}} \boldsymbol{Q} \boldsymbol{d}^{k}}$ 。

注意:αk\alpha_{k}是步长,λk\lambda_{k}是组合系数,二者不同。

针对正定二次函数迭代过程中最优步长 αk\alpha_{k} 的推导过程:
不失一般性,设 f(x)=12xTQx+bTx+cf(\boldsymbol{x})=\frac{1}{2} \boldsymbol{x}^{\mathrm{T}} \boldsymbol{Q} \boldsymbol{x}+\boldsymbol{b}^{\mathrm{T}} \boldsymbol{x}+c ,其中 Q\boldsymbol{Q} 为对称正定矩阵,且 f(x)=Qx+b\nabla f(\boldsymbol{x})=\boldsymbol{Q} \boldsymbol{x}+\boldsymbol{b} 。给定任意迭代点 xk\boldsymbol{x}^{k} ,则最佳步长 αk\alpha_{k} 满足αk=argminα>0f(xk+αdk)\alpha_{k}=\arg \min _{\alpha>0} f\left(\boldsymbol{x}^{k}+\alpha \boldsymbol{d}^{k}\right),即

0=f(xk+αdk)=f(xk+αdk)Tdk=(Q(xk+αdk)+b)Tdk=(f(xk)+αQdk)Tdk\begin{aligned} 0 & =f^{\prime}\left(\boldsymbol{x}^{k}+\alpha \boldsymbol{d}^{k}\right)=\nabla f\left(\boldsymbol{x}^{k}+\alpha \boldsymbol{d}^{k}\right)^{\mathrm{T}} \boldsymbol{d}^{k} \\ & =\left(\boldsymbol{Q}\left(\boldsymbol{x}^{k}+\alpha \boldsymbol{d}^{k}\right)+\boldsymbol{b}\right)^{\mathrm{T}} \boldsymbol{d}^{k}=\left(\nabla f\left(\boldsymbol{x}^{k}\right)+\alpha \boldsymbol{Q} \boldsymbol{d}^{k}\right)^{\mathrm{T}} \boldsymbol{d}^{k} \end{aligned}

解得:$$\alpha_{k}=-\frac{\left(\boldsymbol{d}{k}\right){T} \nabla f\left(\boldsymbol{x}{k}\right)}{\left(\boldsymbol{d}{k}\right)^{T} \boldsymbol{Q} \boldsymbol{d}^{k}}$$

(3)收敛性
综上所述,若中途不停机的话,这样得到的 nn 个方向 d1,d2,,dn\boldsymbol{d}^{1}, \boldsymbol{d}^{2}, \cdots, \boldsymbol{d}^{n} 关于 Q\boldsymbol{Q} 是两两共轭的。再结合基本定理,那么 xn+1x^{n+1} 一定是所求的最优解。

5.组合系数的选取

上面的算法中,组合系数的选取为如下 Hestenes-Stiefel 公式(HS):

λk=(dk)TQf(xk+1)(dk)TQdk\lambda_{k}=\frac{\left(\boldsymbol{d}^{k}\right)^{\mathrm{T}} \boldsymbol{Q} \nabla f\left(\boldsymbol{x}^{k+1}\right)}{\left(\boldsymbol{d}^{k}\right)^{\mathrm{T}} \boldsymbol{Q} \boldsymbol{d}^{k}}

注意到该公式中含有 Hessian 矩阵 Q ,计算复杂度高且所需存储空间较大,故不能直接应用于一般非线性函数求极值的问题。

组合系数的其他形式一正定二次函数:

  1. Fletcher-Reeve 公式(FR)

    λk=f(xk+1)22f(xk)22\lambda_{k}=\frac{\left\|\nabla f\left(\boldsymbol{x}^{k+1}\right)\right\|_{2}^{2}}{\left\|\nabla f\left(\boldsymbol{x}^{k}\right)\right\|_{2}^{2}}

  2. Dixon-Myers 公式(DM)

    λk=f(xk+1)22(dk)Tf(xk)\lambda_{k}=-\frac{\left\|\nabla f\left(\boldsymbol{x}^{k+1}\right)\right\|_{2}^{2}}{\left(\boldsymbol{d}^{k}\right)^{\mathrm{T}} \nabla f\left(\boldsymbol{x}^{k}\right)}

  3. Polak-Ribiere-Polyak 公式(PRP)

    λk=f(xk+1)T(f(xk+1)f(xk))f(xk)22\lambda_{k}=\frac{\nabla f\left(\boldsymbol{x}^{k+1}\right)^{\mathrm{T}}\left(\nabla f\left(\boldsymbol{x}^{k+1}\right)-\nabla f\left(\boldsymbol{x}^{k}\right)\right)}{\left\|\nabla f\left(\boldsymbol{x}^{k}\right)\right\|_{2}^{2}}

    对于正定二次函数,HS、FR、DM、PRP 等四个公式是等价的;对应的四种共轭梯度法也是等价的。对于非二次函数,产生的搜索方向不再相同,常利用 FR、DM、PRP 等三个公式,通常不用 HS 公式 (其中含有 Hessian 矩阵)。经验上 PRP 效果较好。
6.求解一般函数极值的 FR 共轭梯度法
  1. 给定 x1Rn\boldsymbol{x}^{1} \in \mathbb{R}^{n}d1=f(x1)\boldsymbol{d}^{1}=-\nabla f\left(\boldsymbol{x}^{1}\right)0ϵ<10 \leq \epsilon<1k=1k=1
  2. 如果 f(xk)2ϵ\left\|\nabla f\left(\boldsymbol{x}^{k}\right)\right\|_{2} \leq \epsilon ,停止迭代,取 x=xk\boldsymbol{x}^{*}=\boldsymbol{x}^{k} ;否则,由精确线搜索求步长:αk=argminα>0f(xk+αdk)\alpha_{k}=\arg \min _{\alpha>0} f\left(\boldsymbol{x}^{k}+\alpha \boldsymbol{d}^{k}\right),则 xk+1=xk+αkdk\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}+\alpha_{k} \boldsymbol{d}^{k}dk+1=f(xk+1)+λkdk\boldsymbol{d}^{k+1}=-\nabla f\left(\boldsymbol{x}^{k+1}\right)+\lambda_{k} \boldsymbol{d}^{k} ,其中 λk=f(xk+1)22f(xk)22\lambda_{k}=\frac{\left\|\nabla f\left(\boldsymbol{x}^{k+1}\right)\right\|_{2}^{2}}{\left\|\nabla f\left(\boldsymbol{x}^{k}\right)\right\|_{2}^{2}}
  3. k=k+1k=k+1 ,返回 Step 2。

收敛性定理:
设凸函数 f(x)f(\boldsymbol{x}) 存在一阶连续偏导数,且水平集有界L={xRnf(x)f(x1)}L=\left\{\boldsymbol{x} \in \mathbb{R}^{n} \mid f(\boldsymbol{x}) \leq f\left(\boldsymbol{x}^{1}\right)\right\},则由共轭梯度法得到的点列 {xk}\left\{\boldsymbol{x}^{k}\right\} 有如下性质:

  1. f(xk)f\left(\boldsymbol{x}^{k}\right) 严格单调下降,且 limk+f(xk)\lim _{k \rightarrow+\infty} f\left(\boldsymbol{x}^{k}\right) 存在;
  2. {xk}\left\{\boldsymbol{x}^{k}\right\} 的任意驻点 x\boldsymbol{x}^{*} 都是 f(x)f(\boldsymbol{x}) 的极小点。

优点:

  • 克服了最速下降法收敛速度慢的缺点,共轭梯度法的收敛速率不坏于最速下降法;
  • 建立在二次模型之上,具有二次终止性,即求解 n 元正定二次函数(严格凸二次函数)极值,最多 n 步迭代便可求出最优解;
  • 不用求 Hessian 矩阵或者逆矩阵,算法所需存储空间小,是求解大规模问题的重要方法;
  • 对凸函数全局收敛。

缺点:

  • 误差可能会使 n 步迭代得不到正定二次函数的极小点;
  • n 维空间中共轭方向最多有 n 个, n 步后构造的搜索方向不再是共轭的,会降低收敛速度;
  • 对于非二次函数,由于目标函数的 Hessian 矩阵不再是常数矩阵,因而产生的方向不再是共轭方向。
7.周期性的共轭梯度法

对于 FR 算法和 PRP 算法,如果初始方向不取负梯度方向或不使用精确线搜索求最优步长,即使对于二次函数也很难产生 nn 个共轭方向。

重开始技术:
用这两个方法迭代 nn 步后,将 xn+1x^{n+1} 作为初始点 x1x^{1} 重新开始迭代。具体而言,当迭代 nn 步后,将 xn+1x^{n+1} 的负梯度方向设定为搜索方向,然后用共轭梯度法继续求解以产生 nn 个新点列。不断重复上述迭代过程,直至收敛或满足给定的终止条件。

周期性 FR 共轭梯度法算法:

  1. 初始化 x0Rnx^{0} \in \mathbb{R}^{n}0ϵ10 \leq \epsilon \ll 1k=0k=0
  2. 计算 f(xk)\nabla f\left(\boldsymbol{x}^{k}\right) ,如果 f(xk)2ϵ\left\|\nabla f\left(\boldsymbol{x}^{k}\right)\right\|_{2} \leq \epsilon ,停止迭代;否则转 Step 3;
  3. kknn 的倍数,则令 dk=f(xk)\boldsymbol{d}^{k}=-\nabla f\left(\boldsymbol{x}^{k}\right) (重开始);否则令 dk=f(xk)+f(xk)22f(xk1)22dk1\boldsymbol{d}^{k}=-\nabla f\left(\boldsymbol{x}^{k}\right)+\frac{\left\|\nabla f\left(\boldsymbol{x}^{k}\right)\right\|_{2}^{2}}{\left\|\nabla f\left(\boldsymbol{x}^{k-1}\right)\right\|_{2}^{2}} \boldsymbol{d}^{k-1}
  4. 由精确线搜索求 αk\alpha_{k} ,令 xk+1=xk+αkdk\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}+\alpha_{k} \boldsymbol{d}^{k}k=k+1k=k+1 ,返回 Step 2。

三、牛顿法(牛顿法和阻尼牛顿法)

1.牛顿法

基本思想:
f(x)f(\boldsymbol{x}) 是二阶连续可微的,由略去高阶项的 Taylor 公式f(x)f(xk)+f(xk)T(xxk)+12(xxk)T2f(xk)(xxk)f(\boldsymbol{x}) \approx f\left(\boldsymbol{x}^{k}\right)+\nabla f\left(\boldsymbol{x}^{k}\right)^{\mathrm{T}}\left(\boldsymbol{x}-\boldsymbol{x}^{k}\right)+\frac{1}{2}\left(\boldsymbol{x}-\boldsymbol{x}^{k}\right)^{\mathrm{T}} \nabla^{2} f\left(\boldsymbol{x}^{k}\right)\left(\boldsymbol{x}-\boldsymbol{x}^{k}\right)
两边对 xx 求导,且令其等于 0 :f(xk)+2f(xk)(xxk)=0\nabla f\left(\boldsymbol{x}^{k}\right)+\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\left(\boldsymbol{x}-\boldsymbol{x}^{k}\right)=\mathbf{0},若 2f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right) 可逆,则有x=xk[2f(xk)]1f(xk)\boldsymbol{x}=\boldsymbol{x}^{k}-\left[\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\right]^{-1} \nabla f\left(\boldsymbol{x}^{k}\right)。若它不是 f(x)f(x) 的极值点,则令它作为新的迭代点 xk+1x^{k+1} ,即xk+1=xk[2f(xk)]1f(xk)\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}-\left[\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\right]^{-1} \nabla f\left(\boldsymbol{x}^{k}\right)

二次终止性
对于 Hessian 矩阵正定的二次函数,给定任意初始点,牛顿法只需迭代一次就可以得到其极小点。

证明
不失一般性,设 f(x)=12xTQx+bTx+cf(\boldsymbol{x})=\frac{1}{2} \boldsymbol{x}^{\mathrm{T}} \boldsymbol{Q x}+\boldsymbol{b}^{\mathrm{T}} \boldsymbol{x}+c ,其中 Q\boldsymbol{Q} 为对称正定矩阵,则 f(x)f(x) 的极小值点为 x=Q1bx^{*}=-\boldsymbol{Q}^{-1} \boldsymbol{b}
给定初始点 x0x^{0} ,则:f(x0)=Qx0+b\nabla f\left(\boldsymbol{x}^{0}\right)=\boldsymbol{Q} \boldsymbol{x}^{0}+\boldsymbol{b}2f(x0)=Q\nabla^{2} f\left(\boldsymbol{x}^{0}\right)=\boldsymbol{Q} ,故

x1=x0[2f(x0)]1f(x0)=x0Q1(Qx0+b)=Q1b=x\begin{aligned} \boldsymbol{x}^{1} & =\boldsymbol{x}^{0}-\left[\nabla^{2} f\left(\boldsymbol{x}^{0}\right)\right]^{-1} \nabla f\left(\boldsymbol{x}^{0}\right) \\ & =\boldsymbol{x}^{0}-\boldsymbol{Q}^{-1}\left(\boldsymbol{Q} \boldsymbol{x}^{0}+\boldsymbol{b}\right) \\ & =-\boldsymbol{Q}^{-1} \boldsymbol{b}\\ &=\boldsymbol{x}^{*} \end{aligned}

牛顿法算法:

  1. 初始化 x0Rnx^{0} \in \mathbb{R}^{n}0ϵ10 \leq \epsilon \ll 1k=0k=0
  2. 计算 f(xk)\nabla f\left(\boldsymbol{x}^{k}\right) ,如果 f(xk)2ϵ\left\|\nabla f\left(\boldsymbol{x}^{k}\right)\right\|_{2} \leq \epsilon ,停止迭代;否则计算 2f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right) ,并令

    dk=[2f(xk)]1f(xk)\boldsymbol{d}^{k}=-\left[\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\right]^{-1} \nabla f\left(\boldsymbol{x}^{k}\right)

  3. xk+1=xk+dk,k=k+1\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}+\boldsymbol{d}^{k}, k=k+1 ,返回 Step 2。

如何计算可逆矩阵 A 的逆矩阵?
(AI)(A \mid I) 进行初等行变换。当 (AI)(\boldsymbol{A} \mid \boldsymbol{I})A\boldsymbol{A} 变换为 I\boldsymbol{I} 时,则原 I\boldsymbol{I} 变换为 A1\boldsymbol{A}^{-1} ,即 (IA1)\left(I \mid A^{-1}\right)

优点:

  • 对于正定二次函数,牛顿法从任意初始点迭代一次即可得到极小点;
  • 若最优解 x\boldsymbol{x}^{*}2f(x)\nabla^{2} f\left(\boldsymbol{x}^{*}\right) 正定,且初始点选取合适,算法很快收敛。

缺点:

  • 要求函数至少二阶连续可微;
    -收敛性对初始点的依赖性很强,即不具备全局收敛性,有可能收敛到鞍点或极大值点;
  • 每次迭代需要计算 Hessian 矩阵 2f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right) ,计算复杂度高;
  • 每次迭代需要求解方程组 2f(xk)dk=f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right) \boldsymbol{d}^{k}=-\nabla f\left(\boldsymbol{x}^{k}\right) ,当方程组为奇异(矩阵对应的行列式值为零,无法求逆矩阵)或病态(指系数矩阵或常数项的微小扰动会引起解剧烈变化的线性方程组)时,无法求得 dk\boldsymbol{d}^{k}dk\boldsymbol{d}^{k} 不是下降方向。
2.阻尼牛顿法

基本思想:
针对牛顿法的第一个缺点,在求新迭代点时,以 dk\boldsymbol{d}^{k} 作为搜索方向进行一维搜索求步长 αk\alpha_{k} ,而不是将步长固定为 αk=1\alpha_{k}=1αk=argminα>0f(xk+αdk)\alpha_{k}=\arg \min _{\alpha>0} f\left(\boldsymbol{x}^{k}+\alpha \boldsymbol{d}^{k}\right),即利用精确一维搜索为牛顿法搜索得最优步长 αk\alpha_{k} ,然后令 xk+1=xk+αkdk=xkαk[2f(xk)]1f(xk)\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}+\alpha_{k} \boldsymbol{d}^{k}=\boldsymbol{x}^{k}-\alpha_{k}\left[\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\right]^{-1} \nabla f\left(\boldsymbol{x}^{k}\right)

二次终止性:
对于 Hessian 矩阵正定的二次函数,给定任意初始点,阻尼牛顿法只需迭代一次就可以得到其极小点。

阻尼牛顿法算法

  1. 初始化 x0Rnx^{0} \in \mathbb{R}^{n}0ϵ10 \leq \epsilon \ll 1k=0k=0
  2. 计算 f(xk)\nabla f\left(\boldsymbol{x}^{k}\right) ,若 f(xk)2ϵ\left\|\nabla f\left(\boldsymbol{x}^{k}\right)\right\|_{2} \leq \epsilon ,停止迭代;否则计算 2f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right) ,并令dk=[2f(xk)]1f(xk)\boldsymbol{d}^{k}=-\left[\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\right]^{-1} \nabla f\left(\boldsymbol{x}^{k}\right)
  3. 沿 dk\boldsymbol{d}^{k} 进行精确线搜索,获得最优步长 αk\alpha_{k}
  4. xk+1=xk+αkdk\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}+\alpha_{k} \boldsymbol{d}^{k}k=k+1k=k+1 ,转 Step 2。

收敛性定理:
f(x)f(\boldsymbol{x}) 二阶连续可微, 2f(x)\nabla^{2} f(\boldsymbol{x}) 正定。记 {xk}\left\{\boldsymbol{x}^{k}\right\} 是由阻尼牛顿法所得迭代点列,若水平集 L={xRnf(x)f(x0)}L=\left\{\boldsymbol{x} \in \mathbb{R}^{n} \mid f(\boldsymbol{x}) \leq f\left(\boldsymbol{x}^{0}\right)\right\} 有界,则:

  1. {f(xk)}\left\{f\left(\boldsymbol{x}^{k}\right)\right\} 为严格单调下降数列;
  2. {xk}\left\{\boldsymbol{x}^{k}\right\} 必有聚点,且任意聚点 x\boldsymbol{x}^{*} 满足 f(x)=0\nabla f\left(\boldsymbol{x}^{*}\right)=\mathbf{0}

针对每次迭代都要计算 Hessian 矩阵、计算量大的改进方法:

  • 为减小工作量,可令每 mm (正整数)次迭代使用同一个 Hessian 矩阵。但收敛速度随 mm 的增大而下降。具体而言,随 m+m \rightarrow+\infty ,算法收敛速度逐渐降低为线性收敛;
  • 采用拟牛顿法,随着迭代的进程逐渐近似 Hessian 矩阵 2f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right)

针对 2f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right) 非正定和奇异的改进方法:

  • dk\boldsymbol{d}^{k} 为函数上升方向时,可沿其负方向搜索。但可能出现在 ±dk\pm \boldsymbol{d}^{k} 方向上函数值都不变化的情况;
  • μ>0\mu>0 足够小,且使得 2f(xk)+μI\nabla^{2} f\left(\boldsymbol{x}^{k}\right)+\mu \boldsymbol{I} 正定,而后用 2f(xk)+μI\nabla^{2} f\left(\boldsymbol{x}^{k}\right)+\mu \boldsymbol{I} 取代 2f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right)进行迭代,且该策略为全局二阶收敛;
  • 2f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right) 正定,则 dk=[2f(xk)]1f(xk)\boldsymbol{d}^{k}=-\left[\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\right]^{-1} \nabla f\left(\boldsymbol{x}^{k}\right) ;否则 dk=f(xk)\boldsymbol{d}^{k}=-\nabla f\left(\boldsymbol{x}^{k}\right) 。同时,采用非精确一维搜索搜索步长 αk\alpha_{k} 。在一定条件下,该方法全局收敛,但当迭代过程中 2f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right) 非正定情况较多时,收敛速度降为接近线性;
  • 采用拟牛顿法,随着迭代不断近似 Hessian 矩阵的逆 [2f(xk)]1\left[\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\right]^{-1}

四、拟牛顿法(变尺度法)

1.拟牛顿法基本思想

牛顿方向:dk=[2f(xk)]1f(xk)\boldsymbol{d}^{k}=-\left[\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\right]^{-1} \nabla f\left(\boldsymbol{x}^{k}\right),可以看到,牛顿法在每次迭代时都需要计算 Hessian 矩阵 2f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right) 及其逆矩阵,计算量和存储量都较大。1959 年,Davidon 提出用迭代过程中得到的梯度信息来近似 Hessian 矩阵或其逆矩阵,由此产生了一类非常成功的算法——拟牛顿法。

目前最常用的两个拟牛顿法:

  • Davidon-Fletcher-Powell(DFP)算法
  • Broyden-Fletcher-Goldfarb-Shanno(BFGS)算法
2.DFP 算法

Hk\boldsymbol{H}_{k}[2f(xk)]1\left[\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\right]^{-1} 的近似矩阵,则新算法中每次的搜索方向为dk=Hkf(xk)\boldsymbol{d}^{k}=-\boldsymbol{H}_{k} \nabla f\left(\boldsymbol{x}^{k}\right)。由此产生的方法称为变尺度法, Hk\boldsymbol{H}_{k} 称为尺度矩阵。
要使新算法比牛顿法计算和存储量小且同时整体收敛性好,关键在于:如何构造近似矩阵阵列 {Hk}\left\{\boldsymbol{H}_{k}\right\}

要求: {Hk}\left\{\boldsymbol{H}_{k}\right\} 既能逐步逼近 [2f(xk)]1\left[\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\right]^{-1} ,又无需计算二阶导数。同时保证算法的二次收敛性和稳定性。

{Hk}\left\{\boldsymbol{H}_{k}\right\} 需满足以下条件:
C1Hk\boldsymbol{H}_{k} 是对称正定矩阵;
C2Hk\boldsymbol{H}_{k}Hk1\boldsymbol{H}_{k-1} 简单修正而得:Hk=Hk1+Ek\boldsymbol{H}_{k}=\boldsymbol{H}_{k-1}+\boldsymbol{E}_{k}
C3Hk\boldsymbol{H}_{k} 满足拟牛顿方程(与牛顿法类似)Hkγk=δk\boldsymbol{H}_{k} \boldsymbol{\gamma}_{k}=\boldsymbol{\delta}_{k},其中 δk=xkxk1\boldsymbol{\delta}_{k}=\boldsymbol{x}^{k}-\boldsymbol{x}^{k-1}γk=f(xk)f(xk1)\boldsymbol{\gamma}_{k}=\nabla f\left(\boldsymbol{x}^{k}\right)-\nabla f\left(\boldsymbol{x}^{k-1}\right)

以下给出 {Hk}\left\{\boldsymbol{H}_{k}\right\} 需满足 C3 的原因:

f(x)f(\boldsymbol{x}) 是二阶连续可微的,由 Taylor 公式

f(x)f(xk)+f(xk)T(xxk)+12(xxk)T2f(xk)(xxk),f(x)f(xk)+2f(xk)(xxk).\begin{array}{l} f(\boldsymbol{x}) \approx f\left(\boldsymbol{x}^{k}\right)+\nabla f\left(\boldsymbol{x}^{k}\right)^{\mathrm{T}}\left(\boldsymbol{x}-\boldsymbol{x}^{k}\right)+\frac{1}{2}\left(\boldsymbol{x}-\boldsymbol{x}^{k}\right)^{\mathrm{T}} \nabla^{2} f\left(\boldsymbol{x}^{k}\right)\left(\boldsymbol{x}-\boldsymbol{x}^{k}\right), \\ \nabla f(\boldsymbol{x}) \approx \nabla f\left(\boldsymbol{x}^{k}\right)+\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\left(\boldsymbol{x}-\boldsymbol{x}^{k}\right) . \end{array}

x=xk1\boldsymbol{x}=\boldsymbol{x}^{k-1} ,则:2f(xk)(xkxk1)f(xk)f(xk1)\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\left(\boldsymbol{x}^{k}-\boldsymbol{x}^{k-1}\right) \approx \nabla f\left(\boldsymbol{x}^{k}\right)-\nabla f\left(\boldsymbol{x}^{k-1}\right)
δk=xkxk1\boldsymbol{\delta}_{k}=\boldsymbol{x}^{k}-\boldsymbol{x}^{k-1}, γk=f(xk)f(xk1)\boldsymbol{\gamma}_{k}=\nabla f\left(\boldsymbol{x}^{k}\right)-\nabla f\left(\boldsymbol{x}^{k-1}\right) ,因此 2f(xk)δkγk\nabla^{2} f\left(\boldsymbol{x}^{k}\right) \boldsymbol{\delta}_{k} \approx \gamma_{k} (对二次函数为等号)。
2f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right) 可逆,则有 [2f(xk)]1γkδk\left[\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\right]^{-1} \gamma_{k} \approx \boldsymbol{\delta}_{k}
Hk\boldsymbol{H}_{k} 代替 [2f(xk)]1\left[\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\right]^{-1} 并强制上式取等号,得到拟牛顿方程。

对称秩 1 校正:

设校正矩阵的形式为:Hk=Hk1+βkukukT\boldsymbol{H}_{k}=\boldsymbol{H}_{k-1}+\beta_{k} \boldsymbol{u}_{k} \boldsymbol{u}_{k}^{\mathrm{T}},其中 βk0\beta_{k} \neq 0uk=(uk1,uk2,,ukn)T0\boldsymbol{u}_{k}=\left(u_{k 1}, u_{k 2}, \cdots, u_{k n}\right)^{\mathrm{T}} \neq \mathbf{0}
代入拟牛顿方程 Hkγk=δk\boldsymbol{H}_{k} \boldsymbol{\gamma}_{k}=\boldsymbol{\delta}_{k} 得:
Hk1γk+βkukukTγk=δk\boldsymbol{H}_{k-1} \boldsymbol{\gamma}_{k}+\beta_{k} \boldsymbol{u}_{k} \boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k}=\boldsymbol{\delta}_{k},即: βkuk(ukTγk)=δkHk1γk\beta_{k} \boldsymbol{u}_{k}\left(\boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k}\right)=\boldsymbol{\delta}_{k}-\boldsymbol{H}_{k-1} \boldsymbol{\gamma}_{k}
注意到 βk\beta_{k}ukTγk\boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k} 都是标量,所以向量 uk\boldsymbol{u}_{k}δkHk1γk\boldsymbol{\delta}_{k}-\boldsymbol{H}_{k-1} \boldsymbol{\gamma}_{k} 成比例。因此不失一般性,不妨取 βk(ukTγk)=1\beta_{k}\left(\boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k}\right)=1 ,则有 uk=δkHk1γk\boldsymbol{u}_{k}=\boldsymbol{\delta}_{k}-\boldsymbol{H}_{k-1} \boldsymbol{\gamma}_{k}

此时βk=1ukTγk=1γkTuk=1γkT(δkHk1γk)\beta_{k}=\frac{1}{\boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k}}=\frac{1}{\boldsymbol{\gamma}_{k}^{\mathrm{T}} \boldsymbol{u}_{k}}=\frac{1}{\boldsymbol{\gamma}_{k}^{\mathrm{T}}\left(\boldsymbol{\delta}_{k}-\boldsymbol{H}_{k-1} \boldsymbol{\gamma}_{k}\right)}

代入到Hk=Hk1+βkukukT\boldsymbol{H}_{k}=\boldsymbol{H}_{k-1}+\beta_{k} \boldsymbol{u}_{k} \boldsymbol{u}_{k}^{\mathrm{T}}式,得秩1校正公式

Hk=Hk1+(δkHk1γk)(δkHk1γk)TγkT(δkHk1γk)\boldsymbol{H}_{k}=\boldsymbol{H}_{k-1}+\frac{\left(\boldsymbol{\delta}_{k}-\boldsymbol{H}_{k-1} \boldsymbol{\gamma}_{k}\right)\left(\boldsymbol{\delta}_{k}-\boldsymbol{H}_{k-1} \boldsymbol{\gamma}_{k}\right)^{\mathrm{T}}}{\boldsymbol{\gamma}_{k}^{\mathrm{T}}\left(\boldsymbol{\delta}_{k}-\boldsymbol{H}_{k-1} \boldsymbol{\gamma}_{k}\right)}

基于秩 1 校正的拟牛顿算法:

  1. 初始化 x0Rnx^{0} \in \mathbb{R}^{n}H0=I\boldsymbol{H}_{0}=\boldsymbol{I}0ϵ10 \leq \epsilon \ll 1k=0k=0 ,取 d0=f(x0)\boldsymbol{d}^{0}=-\nabla f\left(\boldsymbol{x}^{0}\right)
  2. f(xk)2ϵ\left\|\nabla f\left(\boldsymbol{x}^{k}\right)\right\|_{2} \leq \epsilon ,停止运算,取 x=xk\boldsymbol{x}^{*}=\boldsymbol{x}^{k} ;否则,利用精确线搜索求步长 αk\alpha_{k} ,即 αk=argminf(xk+αdk)\alpha_{k}=\arg \min f\left(\boldsymbol{x}^{k}+\alpha \boldsymbol{d}^{k}\right) ,并令 xk+1=xk+αkdk\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}+\alpha_{k} \boldsymbol{d}^{k}
  3. 计算

    δk+1=xk+1xkγk+1=f(xk+1)f(xk)Hk+1=Hk+(δk+1Hkγk+1)(δk+1Hkγk+1)Tγk+1T(δk+1Hkγk+1)dk+1=Hk+1f(xk+1)\begin{array}{l} \boldsymbol{\delta}_{k+1}=\boldsymbol{x}^{k+1}-\boldsymbol{x}^{k} \\ \boldsymbol{\gamma}_{k+1}=\nabla f\left(\boldsymbol{x}^{k+1}\right)-\nabla f\left(\boldsymbol{x}^{k}\right) \\ \boldsymbol{H}_{k+1}=\boldsymbol{H}_{k}+\frac{\left(\boldsymbol{\delta}_{k+1}-\boldsymbol{H}_{k} \boldsymbol{\gamma}_{k+1}\right)\left(\boldsymbol{\delta}_{k+1}-\boldsymbol{H}_{k} \gamma_{k+1}\right)^{\mathrm{T}}}{\boldsymbol{\gamma}_{k+1}^{\mathrm{T}}\left(\boldsymbol{\delta}_{k+1}-\boldsymbol{H}_{k} \gamma_{k+1}\right)} \\ \boldsymbol{d}^{k+1}=-\boldsymbol{H}_{k+1} \nabla f\left(\boldsymbol{x}^{k+1}\right) \end{array}

  4. 令 k=k+1 ,返回 Step 2。

优点:

  • 算法简洁,具有二次终止性。
    缺点:
  • Hk1\boldsymbol{H}_{k-1} 正定时,由秩 1 校正公式确定的 Hk\boldsymbol{H}_{k} 不一定正定,因此不能保证搜索方向 dk=Hkf(xk)\boldsymbol{d}^{k}=-\boldsymbol{H}_{k} \nabla f\left(\boldsymbol{x}^{k}\right) 是下降方向;
  • 秩 1 校正公式中分母可能为零或者近似于零,导致算法不稳定。

对称秩 2 校正:

设校正矩阵的形式为:Hk=Hk1+λkukukT+βkvkvkT\boldsymbol{H}_{k}=\boldsymbol{H}_{k-1}+\lambda_{k} \boldsymbol{u}_{k} \boldsymbol{u}_{k}^{\mathrm{T}}+\beta_{k} \boldsymbol{v}_{k} \boldsymbol{v}_{k}^{\mathrm{T}},且 βk0\beta_{k} \neq 0λk0\lambda_{k} \neq 0uk=(uk1,uk2,,ukn)T0\boldsymbol{u}_{k}=\left(u_{k 1}, u_{k 2}, \cdots, u_{k n}\right)^{\mathrm{T}} \neq \mathbf{0}vk=(vk1,vk2,,vkn)T0\boldsymbol{v}_{k}=\left(v_{k 1}, v_{k 2}, \cdots, v_{k n}\right)^{\mathrm{T}} \neq \mathbf{0}

代入拟牛顿方程 Hkγk=δk\boldsymbol{H}_{k} \boldsymbol{\gamma}_{k}=\boldsymbol{\delta}_{k} 得到Hk1γk+λkukukTγk+βkvkvkTγk=δk\boldsymbol{H}_{k-1} \boldsymbol{\gamma}_{k}+\lambda_{k} \boldsymbol{u}_{k} \boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k}+\beta_{k} \boldsymbol{v}_{k} \boldsymbol{v}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k}=\boldsymbol{\delta}_{k},即 λkukukTγk+βkvkvkTγk=δkHk1γk\lambda_{k} \boldsymbol{u}_{k} \boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k}+\beta_{k} \boldsymbol{v}_{k} \boldsymbol{v}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k}=\boldsymbol{\delta}_{k}-\boldsymbol{H}_{k-1} \boldsymbol{\gamma}_{k}

λkuk(ukTγk)=δk\lambda_{k} \boldsymbol{u}_{k}\left(\boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k}\right)=\boldsymbol{\delta}_{k}βkvk(vkTγk)=Hk1γk\beta_{k} \boldsymbol{v}_{k}\left(\boldsymbol{v}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k}\right)=-\boldsymbol{H}_{k-1} \boldsymbol{\gamma}_{k} ,则上式成立。
注意到 λk\lambda_{k}ukTγk\boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k}βk\beta_{k}vkTγk\boldsymbol{v}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k} 都是标量,那么若令λk(ukTγk)=1\lambda_{k}\left(\boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k}\right)=1βk(vkTγk)=1\beta_{k}\left(\boldsymbol{v}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k}\right)=-1

则有 uk=δk\boldsymbol{u}_{k}=\boldsymbol{\delta}_{k}, vk=Hk1γk\boldsymbol{v}_{k}=\boldsymbol{H}_{k-1} \boldsymbol{\gamma}_{k}
将这些代入Hk=Hk1+λkukukT+βkvkvkT\boldsymbol{H}_{k}=\boldsymbol{H}_{k-1}+\lambda_{k} \boldsymbol{u}_{k} \boldsymbol{u}_{k}^{\mathrm{T}}+\beta_{k} \boldsymbol{v}_{k} \boldsymbol{v}_{k}^{\mathrm{T}},得DFP公式(秩 2 校正公式):

Hk=Hk1+δkδkTδkTγkHk1γkγkTHk1γkTHk1γk\boldsymbol{H}_{k}=\boldsymbol{H}_{k-1}+\frac{\boldsymbol{\delta}_{k} \boldsymbol{\delta}_{k}^{\mathrm{T}}}{\boldsymbol{\delta}_{k}^{\mathrm{T}} \boldsymbol{\gamma}_{k}}-\frac{\boldsymbol{H}_{k-1} \boldsymbol{\gamma}_{k} \boldsymbol{\gamma}_{k}^{\mathrm{T}} \boldsymbol{H}_{k-1}}{\boldsymbol{\gamma}_{k}^{\mathrm{T}} \boldsymbol{H}_{k-1} \boldsymbol{\gamma}_{k}}

基于秩 2 校正的拟牛顿算法(DFP):

  1. 初始化 x0Rn\boldsymbol{x}^{0} \in \mathbb{R}^{n}H0=I\boldsymbol{H}_{0}=\boldsymbol{I}0ϵ10 \leq \epsilon \ll 1k=0k=0 ,取 d0=f(x0)\boldsymbol{d}^{0}=-\nabla f\left(\boldsymbol{x}^{0}\right)
  2. f(xk)2ϵ\left\|\nabla f\left(\boldsymbol{x}^{k}\right)\right\|_{2} \leq \epsilon ,停止迭代,取 x=xk\boldsymbol{x}^{*}=\boldsymbol{x}^{k} ;否则,利用精确线搜索求步长 αk\alpha_{k} ,即 αk=argminf(xk+αdk)\alpha_{k}=\arg \min f\left(\boldsymbol{x}^{k}+\alpha \boldsymbol{d}^{k}\right) ,并令 xk+1=xk+αkdk\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}+\alpha_{k} \boldsymbol{d}^{k}
  3. 计算

    δk+1=xk+1xk,γk+1=f(xk+1)f(xk),Hk+1=Hk+δk+1δk+1Tδk+1Tγk+1Hkγk+1γk+1THkγk+1THkγk+1,dk+1=Hk+1f(xk+1)\begin{array}{l} \boldsymbol{\delta}_{k+1}=\boldsymbol{x}^{k+1}-\boldsymbol{x}^{k}, \\ \boldsymbol{\gamma}_{k+1}=\nabla f\left(\boldsymbol{x}^{k+1}\right)-\nabla f\left(\boldsymbol{x}^{k}\right), \\ \boldsymbol{H}_{k+1}=\boldsymbol{H}_{k}+\frac{\boldsymbol{\delta}_{k+1} \boldsymbol{\delta}_{k+1}^{\mathrm{T}}}{\boldsymbol{\delta}_{k+1}^{\mathrm{T}} \boldsymbol{\gamma}_{k+1}}-\frac{\boldsymbol{H}_{k} \boldsymbol{\gamma}_{k+1} \boldsymbol{\gamma}_{k+1}^{\mathrm{T}} \boldsymbol{H}_{k}}{\boldsymbol{\gamma}_{k+1}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\gamma}_{k+1}}, \\ \boldsymbol{d}^{k+1}=-\boldsymbol{H}_{k+1} \nabla f\left(\boldsymbol{x}^{k+1}\right) 。 \end{array}

  4. k=k+1k=k+1 ,返回 Step 2。

收敛定理
f(x)f(\boldsymbol{x}) 一阶连续可微, x0Rn\boldsymbol{x}^{0} \in \mathbb{R}^{n} ,且水平集 L={xRnf(x)f(x0)}L=\left\{\boldsymbol{x} \in \mathbb{R}^{n} \mid f(\boldsymbol{x}) \leq f\left(\boldsymbol{x}^{0}\right)\right\} 有界,则由 DFP 算法产生的点列 {xk}\left\{\boldsymbol{x}^{k}\right\} 满足:

  • {xk}\left\{\boldsymbol{x}^{k}\right\} 是有穷点列,则 xk\boldsymbol{x}^{k}f(x)f(\boldsymbol{x}) 的驻点。
  • {xk}\left\{\boldsymbol{x}^{k}\right\} 是无穷点列,它必有极限且任一极限是 f(x)f(\boldsymbol{x}) 的驻点。

优点:

  • 若函数 f(x)f(\boldsymbol{x}) 是二次严格凸函数,则 DFP 算法具有二次终止性;
  • 若函数 f(x)f(\boldsymbol{x}) 是至少一阶连续可微的严格凸函数,则 DFP 算法是全局收敛的;
  • Hk\boldsymbol{H}_{k} 正定, f(xk)0\nabla f\left(\boldsymbol{x}^{k}\right) \neq \mathbf{0} ,则 Hk+1\boldsymbol{H}_{k+1} 正定;
  • DFP 算法克服了秩 1 校正公式中分母为 0 或近似为 0 的缺点;
  • 对非二次函数,DFP 算法的效果也很好,它比最速下降法和共轭梯度法要有效的多,收敛速度也是超线性的。

缺点:

  • 需要的存储量较大,大约需要 O(n2)O\left(n^{2}\right) 个存储单元,因此对于大规模问题,与共轭梯度法相比,计算量、存储量较大;
  • 实际运算中,由于舍入误差的存在以及一维搜索的不精确,DFP算法的稳定性和计算效率都会受到很大的影响,数值计算的稳定性不如后面将介绍的 BFGS 算法;
  • 若求步长时不采用精确一维搜索,则计算效率不如 BFGS 算法。
3.BFGS 算法

在对牛顿方向修正时,若不考虑构造 [2f(xk)]1\left[\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\right]^{-1} 的近似矩阵,而是考虑 Hessian 矩阵 2f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right) 的近似矩阵 Bk\boldsymbol{B}_{k} ,同样可构造出另一类拟牛顿算法,其中 BFGS 算法是其中之一。
Bk\boldsymbol{B}_{k}2f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right) 的近似矩阵,则新算法中每次搜索方向满足Bkdk+f(xk)=0\boldsymbol{B}_{k} \boldsymbol{d}^{k}+\nabla f\left(\boldsymbol{x}^{k}\right)=\mathbf{0}

思考:要使新算法比牛顿法计算和存储量小,且整体收敛性好,关键在于:

如何构造近似矩阵阵列 {Bk}\left\{\boldsymbol{B}_{k}\right\}
要求: {Bk}\left\{\boldsymbol{B}_{k}\right\} 的选取既能逐步逼近 2f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right) ,又无需计算二阶导数。同时保证算法的二次收敛性和稳定性。

需具备以下条件:
C1Bk\boldsymbol{B}_{k} 是对称正定矩阵;
C2Bk\boldsymbol{B}_{k}Bk1\boldsymbol{B}_{k-1} 简单修正而得:Bk=Bk1+Ck\boldsymbol{B}_{k}=\boldsymbol{B}_{k-1}+\boldsymbol{C}_{k}
C3Bk\boldsymbol{B}_{k} 满足拟牛顿方程:Bkδk=γk\boldsymbol{B}_{k} \boldsymbol{\delta}_{k}=\gamma_{k}

其中 δk=xkxk1\boldsymbol{\delta}_{k}=\boldsymbol{x}^{k}-\boldsymbol{x}^{k-1}γk=f(xk)f(xk1)\boldsymbol{\gamma}_{k}=\nabla f\left(\boldsymbol{x}^{k}\right)-\nabla f\left(\boldsymbol{x}^{k-1}\right)

以下给出需满足 C3 的原因:

f(x)f(\boldsymbol{x}) 是二阶连续可微的,由 Taylor 公式

f(x)f(xk)+f(xk)T(xxk)+12(xxk)T2f(xk)(xxk)f(x)f(xk)+2f(xk)(xxk)\begin{array}{l} f(\boldsymbol{x}) \approx f\left(\boldsymbol{x}^{k}\right)+\nabla f\left(\boldsymbol{x}^{k}\right)^{\mathrm{T}}\left(\boldsymbol{x}-\boldsymbol{x}^{k}\right)+\frac{1}{2}\left(\boldsymbol{x}-\boldsymbol{x}^{k}\right)^{\mathrm{T}} \nabla^{2} f\left(\boldsymbol{x}^{k}\right)\left(\boldsymbol{x}-\boldsymbol{x}^{k}\right) \\ \nabla f(\boldsymbol{x}) \approx \nabla f\left(\boldsymbol{x}^{k}\right)+\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\left(\boldsymbol{x}-\boldsymbol{x}^{k}\right) \end{array}

x=xk1\boldsymbol{x}=\boldsymbol{x}^{k-1} ,则:2f(xk)(xkxk1)f(xk)f(xk1)\nabla^{2} f\left(\boldsymbol{x}^{k}\right)\left(\boldsymbol{x}^{k}-\boldsymbol{x}^{k-1}\right) \approx \nabla f\left(\boldsymbol{x}^{k}\right)-\nabla f\left(\boldsymbol{x}^{k-1}\right)
δk=xkxk1\boldsymbol{\delta}_{k}=\boldsymbol{x}^{k}-\boldsymbol{x}^{k-1}γk=f(xk)f(xk1)\boldsymbol{\gamma}_{k}=\nabla f\left(\boldsymbol{x}^{k}\right)-\nabla f\left(\boldsymbol{x}^{k-1}\right),因此 2f(xk)δkγk\nabla^{2} f\left(\boldsymbol{x}^{k}\right) \boldsymbol{\delta}_{k} \approx \gamma_{k} (对二次函数为等式)。

Bk\boldsymbol{B}_{k} 代替上式中的 2f(xk)\nabla^{2} f\left(\boldsymbol{x}^{k}\right) 并强制上式取等号,得到拟牛顿方程。

设校正矩阵的形式为:Bk=Bk1+λkukukT+βkvkvkT\boldsymbol{B}_{k}=\boldsymbol{B}_{k-1}+\lambda_{k} \boldsymbol{u}_{k} \boldsymbol{u}_{k}^{\mathrm{T}}+\beta_{k} \boldsymbol{v}_{k} \boldsymbol{v}_{k}^{\mathrm{T}}

λk0\lambda_{k} \neq 0βk0\beta_{k} \neq 0uk=(uk1,uk2,,ukn)T0\boldsymbol{u}_{k}=\left(u_{k 1}, u_{k 2}, \cdots, u_{k n}\right)^{\mathrm{T}} \neq \mathbf{0}vk=(vk1,vk2,,vkn)T0\boldsymbol{v}_{k}=\left(v_{k 1}, v_{k 2}, \cdots, v_{k n}\right)^{\mathrm{T}} \neq \mathbf{0}
代入拟牛顿方程 Bkδk=γk\boldsymbol{B}_{k} \boldsymbol{\delta}_{k}=\boldsymbol{\gamma}_{k}
得到 Bk1δk+λkukukTδk+βkvkvkTδk=γk\boldsymbol{B}_{k-1} \boldsymbol{\delta}_{k}+\lambda_{k} \boldsymbol{u}_{k} \boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{\delta}_{k}+\beta_{k} \boldsymbol{v}_{k} \boldsymbol{v}_{k}^{\mathrm{T}} \boldsymbol{\delta}_{k}=\boldsymbol{\gamma}_{k}
λkukukTδk+βkvkvkTδk=γkBk1δk\lambda_{k} \boldsymbol{u}_{k} \boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{\delta}_{k}+\beta_{k} \boldsymbol{v}_{k} \boldsymbol{v}_{k}^{\mathrm{T}} \boldsymbol{\delta}_{k}=\boldsymbol{\gamma}_{k}-\boldsymbol{B}_{k-1} \boldsymbol{\delta}_{k}

λkukukTδk+βkvkvkTδk=γkBk1δk\lambda_{k} \boldsymbol{u}_{k} \boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{\delta}_{k}+\beta_{k} \boldsymbol{v}_{k} \boldsymbol{v}_{k}^{\mathrm{T}} \boldsymbol{\delta}_{k}=\boldsymbol{\gamma}_{k}-\boldsymbol{B}_{k-1} \boldsymbol{\delta}_{k}
可见若 λkuk(ukTδk)=γk\lambda_{k} \boldsymbol{u}_{k}\left(\boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{\delta}_{k}\right)=\boldsymbol{\gamma}_{k}βkvk(vkTδk)=Bk1δk\beta_{k} \boldsymbol{v}_{k}\left(\boldsymbol{v}_{k}^{\mathrm{T}} \boldsymbol{\delta}_{k}\right)=-\boldsymbol{B}_{k-1} \boldsymbol{\delta}_{k} 上式成立。
注意到 λk\lambda_{k}ukTδk\boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{\delta}_{k}βk\beta_{k}vkTδk\boldsymbol{v}_{k}^{\mathrm{T}} \boldsymbol{\delta}_{k} 都是标量,那么若令 λk(ukTδk)=1\lambda_{k}\left(\boldsymbol{u}_{k}{ }^{\mathrm{T}} \boldsymbol{\delta}_{k}\right)=1βk(vkTδk)=1\beta_{k}\left(\boldsymbol{v}_{k}{ }^{\mathrm{T}} \boldsymbol{\delta}_{k}\right)=-1 ,则有 uk=γk\boldsymbol{u}_{k}=\boldsymbol{\gamma}_{k}vk=Bk1δk\boldsymbol{v}_{k}=\boldsymbol{B}_{k-1} \boldsymbol{\delta}_{k}

将这些代入Bk=Bk1+λkukukT+βkvkvkT\boldsymbol{B}_{k}=\boldsymbol{B}_{k-1}+\lambda_{k} \boldsymbol{u}_{k} \boldsymbol{u}_{k}^{\mathrm{T}}+\beta_{k} \boldsymbol{v}_{k} \boldsymbol{v}_{k}^{\mathrm{T}},得BFGS 校正公式:

Bk=Bk1+γkγkTγkTδkBk1δkδkTBk1δkTBk1δk\boldsymbol{B}_{k}=\boldsymbol{B}_{k-1}+\frac{\boldsymbol{\gamma}_{k} \boldsymbol{\gamma}_{k}^{\mathrm{T}}}{\boldsymbol{\gamma}_{k}^{\mathrm{T}} \boldsymbol{\delta}_{k}}-\frac{\boldsymbol{B}_{k-1} \boldsymbol{\delta}_{k} \boldsymbol{\delta}_{k}^{\mathrm{T}} \boldsymbol{B}_{k-1}}{\boldsymbol{\delta}_{k}^{\mathrm{T}} \boldsymbol{B}_{k-1} \boldsymbol{\delta}_{k}}

BFGS 算法:

  1. 初始化 x0Rn\boldsymbol{x}^{0} \in \mathbb{R}^{n}B0=I\boldsymbol{B}_{0}=\boldsymbol{I}0ϵ10 \leq \epsilon \ll 1k=0k=0,取 d0=f(x0)\boldsymbol{d}^{0}=-\nabla f\left(\boldsymbol{x}^{0}\right)
  2. f(xk)2ϵ\left\|\nabla f\left(\boldsymbol{x}^{k}\right)\right\|_{2} \leq \epsilon ,停止运算,取 x=xk\boldsymbol{x}^{*}=\boldsymbol{x}^{k} ;否则,利用精确线搜索求步长 αk\alpha_{k} ,即 αk=argminf(xk+αdk)\alpha_{k}=\arg \min f\left(\boldsymbol{x}^{k}+\alpha \boldsymbol{d}^{k}\right) ,并令 xk+1=xk+αkdk\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}+\alpha_{k} \boldsymbol{d}^{k}
  3. 计算

    δk+1=xk+1xk,γk+1=f(xk+1)f(xk),Bk+1=Bk+γk+1γk+1T1k+1TBkδk+1δk+1TBkδk+1TBkδk+1\begin{array}{l} \boldsymbol{\delta}_{k+1}=\boldsymbol{x}^{k+1}-\boldsymbol{x}^{k}, \\ \boldsymbol{\gamma}_{k+1}=\nabla f\left(\boldsymbol{x}^{\boldsymbol{k}+1}\right)-\nabla f\left(\boldsymbol{x}^{\boldsymbol{k}}\right), \\ \boldsymbol{B}_{k+1}=\boldsymbol{B}_{k}+\frac{\gamma_{k+1}}{\gamma_{k+1}^{\mathrm{T}} \boldsymbol{1}_{k+1}^{\mathrm{T}}}-\frac{\boldsymbol{B}_{k} \boldsymbol{\delta}_{k+1} \boldsymbol{\delta}_{k+1}^{\mathrm{T}} \boldsymbol{B}_{k}}{\boldsymbol{\delta}_{k+1}^{\mathrm{T}} \boldsymbol{B}_{k} \boldsymbol{\delta}_{k+1}} 。 \end{array}

    解方程组 Bk+1dk+1+f(xk+1)=0\boldsymbol{B}_{k+1} \boldsymbol{d}^{k+1}+\nabla f\left(\boldsymbol{x}^{k+1}\right)=\mathbf{0}dk+1\boldsymbol{d}^{k+1}
  4. k=k+1k=k+1 ,返回 Step 2。

BFGS 算法与 DFP 算法的比较:
BFGS 算法与 DFP 算法都属于拟牛顿算法,二者性质有一定类似性。 BFGS 算法的数值稳定性比 DFP 算法好,而且求步长若不用精确一维搜索,仍然是超线性收敛的。

五、最小二乘法

1.最小二乘问题

在某些最优化问题中,比如曲线拟合问题,目标函数由若干个函数的平方和构成。一般可以形式化为:F(x)=i=1mfi2(x)\mathcal{F}(\boldsymbol{x})=\sum_{i=1}^{m} f_{i}^{2}(\boldsymbol{x})
其中 xRn\boldsymbol{x} \in \mathbb{R}^{n} 。一般假设 mnm \geq n 。我们称如下无约束优化问题

minxRnF(x)\min _{\boldsymbol{x} \in \mathbb{R}^{n}} \mathcal{F}(\boldsymbol{x})

为最小二乘问题。特别地,当 fi(x)f_{i}(\boldsymbol{x}) 全为 x\boldsymbol{x} 的线性函数时,称其为线性最小二乘问题。当 fi(x)f_{i}(\boldsymbol{x}) 不全为 x\boldsymbol{x} 的线性函数,即至少一个 fi(x)f_{i}(\boldsymbol{x})x\boldsymbol{x} 的非线性函数时,称其为非线性最小二乘问题。

由于目标函数 F(x)\mathcal{F}(\boldsymbol{x}) 具有若干个函数平方和这种特殊形式,即目标函数具有光滑性质(Smoothness),因此给问题的求解带来便利。对于这类问题,除了能够运用前面介绍的迭代求解方法外,如最速下降法、共轭梯度法、牛顿法、拟牛顿法等,还可以利用解析法求解。

2.线性最小二乘问题

在带参数的回归模型中,最简单的模型是线性回归模型。设 (wi,bi)\left(\boldsymbol{w}_{i}, b_{i}\right) 为观测到的自变量和响应变量,且不同数据点相互独立,则对每个数据点有如下式子成立:

bi=wi,1x1+wi,2x2++wi,n1xn1+xn+εi,i=1,,mb_{i}=w_{i, 1} x_{1}+w_{i, 2} x_{2}+\cdots+w_{i, n-1} x_{n-1}+x_{n}+\varepsilon_{i}, \quad i=1, \cdots, m

其中 xix_{i} 是需要确定的参数, εi\varepsilon_{i} 是某种噪声且不同数据点之间满足独立同分布性质。令 ai=(wiT,1)T\boldsymbol{a}_{i}=\left(\boldsymbol{w}_{i}^{\mathrm{T}}, 1\right)^{\mathrm{T}}x=(x1,x2,,xn)TRn\boldsymbol{x}=\left(x_{1}, x_{2}, \cdots, x_{n}\right)^{\mathrm{T}} \in \mathbb{R}^{n} ,则线性回归模型可以简写为 bi=aiTx+εib_{i}=\boldsymbol{a}_{i}^{\mathrm{T}} \boldsymbol{x}+\varepsilon_{i}

将训练集中的输入特征写成一个矩阵 A ,将标签 b_{i} 和噪声 \varepsilon_{i} 写成向量形式,即

A=(a1TamT)Rm×n,b=(b1bm)Rm,ε=(ε1εm)Rm\boldsymbol{A}=\left(\begin{array}{c} \boldsymbol{a}_{1}^{\mathrm{T}} \\ \vdots \\ \boldsymbol{a}_{m}^{\mathrm{T}} \end{array}\right) \in \mathbb{R}^{m \times n}, \boldsymbol{b}=\left(\begin{array}{c} b_{1} \\ \vdots \\ b_{m} \end{array}\right) \in \mathbb{R}^{m}, \boldsymbol{\varepsilon}=\left(\begin{array}{c} \varepsilon_{1} \\ \vdots \\ \varepsilon_{m} \end{array}\right) \in \mathbb{R}^{m}

则线性回归模型可整理为如下矩阵向量形式:b=Ax+ε0\boldsymbol{b}=\boldsymbol{A} \boldsymbol{x}+\varepsilon_{0}

假设 εi\varepsilon_{i} 是高斯白噪声,即 εiN(0,σ2)\varepsilon_{i} \sim \mathcal{N}\left(0, \sigma^{2}\right) 。那么我们有 p(biai;x)=12πσ2exp((biaiTx)22σ2)p\left(b_{i} \mid \boldsymbol{a}_{i} ; \boldsymbol{x}\right)=\frac{1}{\sqrt{2 \pi} \sigma^{2}} \exp \left(-\frac{\left(b_{i}-\boldsymbol{a}_{i}^{\mathrm{T}} \boldsymbol{x}\right)^{2}}{2 \sigma^{2}}\right)

则对数似然函数为(x)=ln[i=1mp(biai;x)]=m2ln(2π)mlnσ12σ2i=1m(biaiTx)2\ell(\boldsymbol{x})=\ln \left[\prod_{i=1}^{m} p\left(b_{i} \mid \boldsymbol{a}_{i} ; \boldsymbol{x}\right)\right]=-\frac{m}{2} \ln (2 \pi)-m \ln \sigma-\frac{1}{2 \sigma^{2}} \sum_{i=1}^{m}\left(b_{i}-\boldsymbol{a}_{i}^{\mathrm{T}} \boldsymbol{x}\right)^{2}

最大似然估计则是极大化对数似然函数,去除掉常数项之后我们得到了如下最小二乘问题

maxxRn(x)minxRnAxb22\max _{\boldsymbol{x} \in \mathbb{R}^{n}} \ell(\boldsymbol{x}) \Leftrightarrow \min _{\boldsymbol{x} \in \mathbb{R}^{n}}\|\boldsymbol{A} \boldsymbol{x}-\boldsymbol{b}\|_{2}^{2}

由上述,最小二乘问题中的目标函数 F(x)\mathcal{F}(x) 可以表示为:

F(x)=Axb22=(bAx)T(bAx)=bTb(Ax)TbbTAx+xTATAx=bTb((Ax)Tb)TbTAx+xTATAx=bTb2bTAx+xTATAx.\begin{aligned} \mathcal{F}(\boldsymbol{x}) & =\|\boldsymbol{A} \boldsymbol{x}-\boldsymbol{b}\|_{2}^{2}=(\boldsymbol{b}-\boldsymbol{A} \boldsymbol{x})^{\mathrm{T}}(\boldsymbol{b}-\boldsymbol{A} \boldsymbol{x}) \\ & =\boldsymbol{b}^{\mathrm{T}} \boldsymbol{b}-(\boldsymbol{A} \boldsymbol{x})^{\mathrm{T}} \boldsymbol{b}-\boldsymbol{b}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{x}+\boldsymbol{x}^{\mathrm{T}} \boldsymbol{A}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{x} \\ & =\boldsymbol{b}^{\mathrm{T}} \boldsymbol{b}-\left((\boldsymbol{A} \boldsymbol{x})^{\mathrm{T}} \boldsymbol{b}\right)^{\mathrm{T}}-\boldsymbol{b}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{x}+\boldsymbol{x}^{\mathrm{T}} \boldsymbol{A}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{x} \\ & =\boldsymbol{b}^{\mathrm{T}} \boldsymbol{b}-2 \boldsymbol{b}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{x}+\boldsymbol{x}^{\mathrm{T}} \boldsymbol{A}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{x} . \end{aligned}

注意:

  • r=Axb\boldsymbol{r}=\boldsymbol{A} \boldsymbol{x}-\boldsymbol{b} ,则 r\boldsymbol{r} 称为残差(Residual);
  • εi\varepsilon_{i} 不是高斯白噪声时,需要借助似然函数来构造其他噪声所对应的目标函数。例如,在某些噪声下构造出的模型实际上为最小一乘问题:minxRnAxb1\min _{\boldsymbol{x} \in \mathbb{R}^{n}}\|\boldsymbol{A} \boldsymbol{x}-\boldsymbol{b}\|_{1}

易知, F(x)=2ATAx2ATb\nabla \mathcal{F}(\boldsymbol{x})=2 \boldsymbol{A}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{x}-2 \boldsymbol{A}^{\mathrm{T}} \boldsymbol{b} 以及 2F(x)=2ATA\nabla^{2} \mathcal{F}(\boldsymbol{x})=2 \boldsymbol{A}^{\mathrm{T}} \boldsymbol{A} ,因而本章所学习的无约束优化方法可用来求解上述线性最小二乘问题。
另一方面,通过寻找最小二乘问题的驻点也可以获得其解析解。令F(x)=2ATAx2ATb=0\nabla \mathcal{F}(\boldsymbol{x})=2 \boldsymbol{A}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{x}-2 \boldsymbol{A}^{\mathrm{T}} \boldsymbol{b}=\mathbf{0}

也即,最小二乘问题的驻点满足法方程(normal equation)$$\boldsymbol{A}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{x}=\boldsymbol{A}^{\mathrm{T}} \boldsymbol{b}$$
AA 列满秩, ATAA^{\mathrm{T}} A 为对称正定矩阵,由此得到 F(x)\mathcal{F}(x) 的驻点为最小二乘解:$$\boldsymbol{x}{*}=\left(\boldsymbol{A}{\mathrm{T}} \boldsymbol{A}\right)^{-1} \boldsymbol{A}^{\mathrm{T}} \boldsymbol{b}$$由于 F(x)\mathcal{F}(x) 是凸函数, xx^{*} 必是全局最优解。

3.正则化最小二乘法

Tikhonov 正则化
为了平衡模型的拟合性质和解的光滑性,Tikhonov 正则化或岭回归 (ridge regression)添加 2\ell_{2} 范数平方为正则项。假设 εi\varepsilon_{i} 是高斯白噪声,则带 2\ell_{2} 范数平方正则项的线性回归模型为如下问题:minxRnAxb22+λx22\min _{\boldsymbol{x} \in \mathbb{R}^{n}}\|\boldsymbol{A} \boldsymbol{x}-\boldsymbol{b}\|_{2}^{2}+\lambda\|\boldsymbol{x}\|_{2}^{2}
由于正则项的存在,该问题的目标函数是强凸函数,解的性质得到改善。同时,该模型是对 xx 较大的分量产生惩罚,故另一种形式为:minxRnAxb22 s.t. x2μ\min _{\boldsymbol{x} \in \mathbb{R}^{n}}\|\boldsymbol{A} \boldsymbol{x}-\boldsymbol{b}\|_{2}^{2} \text { s.t. }\|\boldsymbol{x}\|_{2} \leq \mu

由于上述两个模型最优性条件类似,当参数 λ\lambdaμ\mu 满足一定关系时,它们的解是相同的。

LASSO 正则化
如果希望得到的解 xx 是稀疏的,则可考虑添加 1\ell_{1} 范数平方为正则项,也即转化为如下 LASSO 问题:minxRnAxb22+λx1\min _{\boldsymbol{x} \in \mathbb{R}^{n}}\|\boldsymbol{A} \boldsymbol{x}-\boldsymbol{b}\|_{2}^{2}+\lambda\|\boldsymbol{x}\|_{1 \bullet}
LASSO 问题通过惩罚参数 λ\lambda 来控制解的稀疏性。如果 xx 是稀疏的,那么预测值 bib_{i} 只和 ai\boldsymbol{a}_{i} 的部分元素相关。同时,也可以考虑:minxRnAxb22 s.t. x1τ\min _{\boldsymbol{x} \in \mathbb{R}^{n}}\|\boldsymbol{A} \boldsymbol{x}-\boldsymbol{b}\|_{2}^{2} \text { s.t. }\|\boldsymbol{x}\|_{1} \leq \tau。或者考虑到噪声 εi\varepsilon_{i} 的存在,转化为如下形式minxRnx1 s.t. Axb22δ\min _{\boldsymbol{x} \in \mathbb{R}^{n}}\|\boldsymbol{x}\|_{1} \text { s.t. }\|\boldsymbol{A} \boldsymbol{x}-\boldsymbol{b}\|_{2}^{2} \leq \delta

Author

秦宇春

Posted on

2025-11-15

Updated on

2025-11-15

Licensed under