本文只考虑最简单的一维问题,包括求根和求最小值点两个子问题:

  • 非线性方程求根 \(f(x)=0\)
  • 非线性方程求最小值点 \(x_0 = \text{argmin} f(x)\)
    • 精确一维搜索
    • 非精确一维搜索

由于线性问题是平凡的,因此默认\(f(x)\)是非线性的。

非线性方程求根

对分法

假设\(f(x)\)在有限区间\([a,b]\)满足:\(f(a)f(b)<0\) 异号,那么\(f(x)\)\((a,b)\)至少有一个零点 \(x^*\)

算法:

  1. 起始的计算区间 \([x_l,x_r]=[a,b]\),此时满足两端异号 \(f(x_l)f(x_r)<0\)
  2. 选取中点 \(x_m = \frac12(x_l+x_r)\) ,计算\(f(x_m)\)
  3. 如果中点 \(|f(x_m)| \le \varepsilon\),或者区间长度\(|x_r-x_l|\le \varepsilon\) 已经足够小,则认为中点已经足够精确,已经找到零点,求解完成。
  4. 否则左右两个区间 \([x_l,x_m]\)\([x_m,x_r]\) 恰有一个满足两端异号,选取它作为新的计算区间\([x_l',x_r']\),回到第二步。

收敛性:进行 \(k\) 次对分之后,计算区间长度为 \(2^{-k}(b-a)\),因此至多进行 \(N\) 次对分。 \[ \frac{b-a}{2^N} \le \varepsilon \Rightarrow N = \left[\log_2(\frac{b-a}{\varepsilon})\right] \]

误差:第 k 次对分之后中点 \(x_m^{(k)}\) 满足 \[ \Vert x_m^{(k)} - x^*\Vert \le \frac{b-a}{2^{k+1}} \]

Remark:算法的前提条件保证了至少有一个零点,但如果在区间中有多个零点,则可能收敛到其中任一零点。

牛顿法

假设 \(f(x)\)\(x\) 附近有一个零点 \(x^*\) ,牛顿法的流程是:从 \(x^{(0)}=x\) 出发,构造迭代点列 \(\{x^{(k)}\}\),希望极限趋于零点 \[ \lim_{k \to \infty} x^{(k)} = x^* \] 具体的构造方法为:在\((x^{(k)},f(x^{(k)}))\)处作 \(f\) 的切线,切线与 x轴交点记作 \(x^{(k+1)}\)\[ \begin{aligned} 0 - f(x^{(k)}) &= f'(x^{(k)})(x^{(k+1)}-x^{(k)})\\ x^{(k+1)} &= x^{(k)} - \frac{f(x^{(k)})}{f'(x^{(k)})} \end{aligned} \] 迭代的终止条件有两个: \[ \begin{aligned} |f(x^{(k)})| \le \varepsilon_1 \\ \Vert x^{(k)} - x^{(k-1)})\Vert \le \varepsilon_2 \end{aligned} \]

Remark:

  • 牛顿法对迭代的初值很敏感,要求初值与真正的零点足够靠近(靠近程度也与 \(f\) 自身性质有关);否则可能跳转并收敛到另外的零点,或者直接不收敛。
  • 在初值合适选取时,并且 \(f \in C^2\)\(x^*\) 有单根时,收敛速度可以达到二阶。(如果在 \(x^*\) 是重根,可能需要对格式进行修改)

\[ \Vert x^{(k+1)}-x^*\Vert \le C \Vert x^{(k)}-x^*\Vert^2 \]

  • 对于单调且凸的函数,才足够保证牛顿法的全局收敛性,例如 \(f(x)=x^2\) 对应的是经典的平方根算法。

收敛速度的证明: 记 \(e^{(k)} = x^{(k)} - x^*\),则有 \[ \begin{aligned} e^{(k+1)} &= x^{(k+1)} - x^* = x^{(k)} - \frac{f(x^{(k)})}{f'(x^{(k)})} - x^* \\ &= e^{(k)} - \frac{f(x^{(k)})}{f'(x^{(k)})} = \frac{e^{(k)} f'(x^{(k)}) - f(x^{(k)})}{f'(x^{(k)})} \end{aligned} \]\(f(x^{*}) = f(x^{(k)}-e^{(k)})\)\(x^{(k)}\) 处进行泰勒展开 \[ 0 = f(x^*) = f(x^{(k)}) - f'(x^{(k)}) e^{(k)} + \frac12 f''(\xi) (e^{(k)})^2 \] 因此有 \[ e^{(k+1)} = \frac12 \frac{f''(\xi)}{f'(x^{(k)})} (e^{(k)})^2 \approx C (e^{(k)})^2 \]

例:经典的平方根算法:计算 \(f(x)=x^2-R\) 的零点 \[ x^{(k+1)} = x^{(k)} - \frac{(x^{(k)})^2-R}{2x^{(k)}} = \frac12 \left( x^{(k)} + \frac{R}{x^{(k)}}\right) \]

Remark: 关于非线性方程求根的牛顿法,可以直接推广到非线性方程组上,求 \(X\in R^n\),满足 \(F(X) = 0 \in R^n\),这里为了简化问题,取约束个数和自变量维数相同,并且假设问题的解存在。最简单是线性情形 \(F(X)=AX\)\(A \in R^{n\times n}\) 是常系数方阵。

迭代公式为 \[ X^{(k+1)} = X^{(k)} - J(X^{(k)})^{-1}F(X^{(k)}) \] 其中 \(J(X)\) 为 Jacobi 矩阵。 \[ J(X)_{ij} = \frac{\partial F_i}{\partial x_j}(X) \] 这里要求 \(F(X)\) 足够光滑,并且假设 Jacobi 矩阵始终可逆。

割线法

牛顿法需要使用导函数信息,有时对于非线性函数 \(f(x)\) 直接求导很困难,可以使用最近两个点的割线斜率进行近似。

\[ \begin{aligned} x^{(k+1)} &= x^{(k)} - \frac{f(x^{(k)})(x^{(k)}-x^{(k-1)})}{f(x^{(k)})-f(x^{(k-1)})} \end{aligned} \] 收敛的速度会略有降低(阶数\(\frac{1+\sqrt{5}}2\approx 1.62\)),此时为了启动格式,我们需要提供两个初始值\(x_0,x_1\)

不动点法

对于方程 \(f(x)=0\) 转化为如下形式 \[ x = F(x) \]\(f\) 的零点 \(x\)\(F\) 的不动点。然后可以据此设计迭代格式 \[ x^{(k+1)} = F(x^{(k)}) \]

基于压缩映射原理,我们希望在区间\([a,b]\)上的 \(F\) 是压缩映射,要求\(F:[a,b]\to [a,b]\) ,并且存在 \(0<L<1\) 使得 \[ \Vert F(x) - F(y) \Vert \le L \Vert x-y \Vert \]

牛顿法实质是一种不动点法: \[ F_N(x) = x - \frac{f(x)}{f'(x)} \] 割线法实质也是一种不动点法: \[ F_G(x) = x - \frac{f(x)^2}{f(x+f(x))-f(x)} \]

不动点法提供了非常自由的构造方式:针对具体问题,设计不同的 \(F\) 并检验压缩映射条件,就可以给出很多不动点迭代格式。

精确一维搜索

精确一维搜索通常要求 \(f\) 在区间 \([a,b]\) 上为单峰的函数,已知存在极小值,否则会计算失败。

黄金分割法

这个方法只需要计算 \(f(x)\) 的值,不需要导数信息。 主要的步骤是不断压缩计算区间,假设初始计算区间 \([x_l,x_r]\),区间长度 \(x_r-x_l\),比例系数 \(\rho \in (0,\frac12)\) 待定。

取内部左右两侧的两个点 \(a_l,a_r\)\[ \begin{aligned} a_l = x_l + \rho (x_r-x_l)\\ a_r = x_r - \rho (x_r-x_l)\\ \end{aligned} \]

  • 如果 \(f(a_l)<f(a_r)\),那么极小点位于 \(a_r\) 左侧,选取新的计算区间为 \([x_l',x_r']=[x_l,a_r]\)
  • 如果 \(f(a_l)\ge f(a_r)\),那么极小点位于 \(a_l\) 右侧,选取新的计算区间为 \([x_l',x_r']=[a_l,x_r]\)

新的计算区间长度为 \((1-\rho)(x_r-x_l)\)

我们希望可以充分利用计算的点值,以第一种情况为例,希望 \(a_l\) 可以在新计算区间 \([x_l',x_r']\) 也被利用,作为右侧点,需要满足 \[ \frac{1-2\rho}{1-\rho} = \rho \] 解得 \(\rho = \frac{3-\sqrt5}2 \approx 0.382\)

(黄金分割法的思路是固定 \(\rho\), 另一种思路是动态调整 \(\rho\),取一个与斐波那契数列相关的系数 \(\rho_k\),得到斐波那契数列法)

每次计算区间的压缩比例为 \(1-\rho \approx 0.618\) 即黄金分割比。 除了第一次压缩区间之外,每一次压缩区间只需要选择一个内部点进行求值。

算法:

  1. 起始的计算区间 \([x_l,x_r]=[a,b]\)
  2. 选取左侧点 \(a_l = x_l + \rho(x_r-x_l)\) ,右侧点 \(a_r = x_r - \rho(x_r-x_l)\),计算 \(f(a_l)\)\(f(a_r)\)
  3. 压缩区间:
    1. 如果 \(f(a_l)<f(a_r)\),那么极小点位于 \(a_r\) 左侧,选取新的计算区间为 \([x_l',x_r']=[x_l,a_r]\)\(a_l\) 可以作为下一步的 \(a_r'\)
    2. 如果 \(f(a_l)\ge f(a_r)\),那么极小点位于 \(a_l\) 右侧,选取新的计算区间为 \([x_l',x_r']=[a_l,x_r]\)\(a_r\) 可以作为下一步的 \(a_l'\)
  4. 如果计算区间已经足够小,则求解结束,得到最小值点,否则回到第二步。

对分法

对分法主要依赖 \(f'(x)\),和求根的对分法基本类似。

算法:

  1. 起始的计算区间 \([x_l,x_r]=[a,b]\)
  2. 选取中点 \(x_m = \frac12(x_l+x_r)\) ,计算 \(f'(x_m)\)
  3. 如果中点 \(|f'(x_m)| \le \varepsilon\),或者区间长度 \(|x_r-x_l|\le \varepsilon\) 已经足够小,则认为已经找到极值点,求解完成。
  4. 否则压缩区间:若 \(f'(x_m) > 0\),取左侧为新的计算区间 \([x_l',x_r'] = [x_l,x_m]\) ;若 \(f'(x_m) < 0\),取右侧为新的计算区间 \([x_l',x_r'] = [x_m,x_r]\)

牛顿法

牛顿法需要依赖 \(f'(x),f''(x)\),并且假设 \(f''>0\) 在计算区域内恒成立,否则难以保证收敛性。

本质是通过二阶泰勒展开将 \(f\)\(x^{(k)}\) 局部近似为二次函数 \(q(h)\) ,然后更新到二次函数 \(q(h)\) 的极小值点,得到 \(x^{(k+1)}\)

\[ q(h)= f(x^{(k)}) + f'(x^{(k)}) h + \frac{f''(x^{(k)})}2 h^2 \approx f(x^{(k)}+h) \]

极小值点 \[ h^* = -\frac{f'(x^{(k)})}{f''(x^{(k)})} \] 因此迭代格式为 \[ x^{(k+1)} = x^{(k)} + h^* = x^{(k)}-\frac{f'(x^{(k)})}{f''(x^{(k)})} \]

Remark: 关于多维的非线性方程最小值问题:自变量 \(X\in R^n\) ,目标函数 \(f: R^n \to R\) 足够光滑,假设目标函数存在极小值,那么牛顿迭代法的公式为 \[ X^{(k+1)} = X^{(k)} - H(X^{(k)})^{-1} \nabla f(X^{(k)}) \] 其中 \(\nabla f\) 是梯度向量,\(H = \nabla^2 f\) 是 Hesse 方阵。 \[ \begin{aligned} \nabla f(X)_i &= \frac{\partial f}{\partial X_i}(X)\\ H(X)_{ij} &= \nabla f(X)_{ij} = \frac{\partial^2 f}{\partial x_i \partial x_j}(X) \end{aligned} \] 要求 Hesse 方阵始终可逆,计算才可以继续,事实上我们还希望 Hesse 方阵始终对称正定,此时二阶近似才存在最小值点,否则可能不收敛。为了确保计算可以继续,从牛顿法发展出了拟牛顿法:用一族对称正定矩阵 \(H_k\) 来近似实际的 Hesse 方阵 \(H(X^{(k)})\) 参与计算。

非精确一维搜索

非精确一维搜索不要求直接求出最小值点,而是给定一个下降准则,求出一个点使得函数有足够的下降即可。

常用的下降准则有

Armijo准则:取定 \(\rho \in (0,1)\) ,要求 \(t\) 使得 \[ f(t) \le f(0) + \rho f'(0) t \] Goldstein准则:取定 \(\rho \in (0,\frac12)\) ,要求 \(t\) 使得 \[ \begin{aligned} f(t) &\le f(0) + \rho f'(0)t\\ f(t) &\ge f(0) + (1-\rho)f'(0) t \end{aligned} \] 注意:真正的最小值点也可能不满足Goldstein准则!

Wolfe准则:取定 \(\rho \in (0,\frac12)\) 以及 \(\sigma \in (\rho,1)\) ,要求 \(t\) 使得 \[ \begin{aligned} f(t) &\le f(0) + \rho f'(0)t \\ f'(t) &\ge \sigma f'(0) \end{aligned} \] 第二个条件通常会被加强为双边条件,称为强Wolfe准则 \[ |f'(t)| \le \sigma |f'(0)| = -\sigma f'(0) \]

基于这些准则可以设计很多具体的非精确一维搜索算法,主要的过程都是不断使用二分法、插值法或者步长乘以比例系数的方法得到新的试探点,判定准则是否成立,成立即可结束搜索,不成立则缩减搜索区间,继续尝试。

以下各有一个例子。

满足 Armijo 准则的回退法

基于Armijo准则设计的常用算法是回退法: 初始选取一个步长 \(t_0 = a\),选取一个指数衰减因子 \(\gamma \in (0,1)\),检查是否满足准则。 如果不满足准则,选取更小的步长 \(t_{k+1} = \gamma t_{k} = \gamma^{k} a\) 继续尝试,直到最终使用足够小的步长\(t\),可以满足准则。

\[ t = \min_{k \in N}\{\gamma^k a : f(\gamma^k a) \le f(0) + \rho f'(0) \gamma^k a \} \]

满足 Goldstein 准则的对分法

取定初始的搜索区间 \([0,a]\),以及两个参数 \(\rho \in (0,\frac12),\delta >0\),计算\(f(0),f'(0)\) 用于准则的判定。

  1. \(t=0\)出发,当前搜索区间是 \([x_l,x_r]=[0,a]\),选取中点作为试探点 \(x =\frac12 (x_l+x_r)\),计算 \(f(x)\)
  2. 判定:验证第一条准则: \(f(x) \le f(0) + \rho f'(0)x\)
    1. 若满足第一条准则,验证第二条准则: \(f(x) \ge f(0) + (1-\rho)f'(0) x\)
      1. 若满足第二条准则,则搜索完成,输出 \(t=x\)
      2. 若不满足第二条准则,搜索区间改为右侧子区间 \([x_l',x_r'] = [x,x_r]\) ,下一次的试探点选取:如果右端点 \(x_r' < a\) 取中点 \(x' = \frac12(x_l'+x_r')\),否则取 \(x' = (1+\delta)x\) 。继续判定
    2. 若不满足第一条准则,搜索区间改为左侧子区间 \([x_l',x_r'] = [x_l,x]\),下一次的试探点仍然取作中点 \(x' = \frac12(x_l'+x_r')\) 。继续判定

满足 Wolfe 准则的算法

取定初始的搜索区间 \([0,a]\),以及两个参数 \(\rho \in (0,\frac12),\sigma \in (\rho,1)\),计算\(f(0),f'(0)\)用于后续判定。

  • \(t=0\)出发,当前搜索区间是 \([x_l,x_r]=[0,a]\),选取适当的试探点 \(x \in (x_l,x_r)\),计算 \(f(x)\)
  • 判定:验证第一条准则 \(f(x) \le f(0) + \rho f'(0)x\)
    1. 若满足第一条准则,验证第二条准则 \(f'(x) \ge \sigma f'(0)\)
      1. 若满足第二条准则,则搜索完成,输出 \(t=x\)
      2. 若不满足第二条准则,搜索区间改为右侧子区间 \([x_l',x_r'] = [x,x_r]\) ,下一次的试探点 \(x'\) 取为一个两点二次插值多项式的极小值点,继续判定
    2. 若不满足第一条准则,搜索区间改为左侧子区间 \([x_l',x_r'] = [x_l,x]\),下一次的试探点 \(x'\) 取为一个两点二次插值多项式的极小值点,继续判定

更具体的细节:

当调整到左侧子区间\([x_l',x_r'] = [x_l,x]\),取两点二次插值多项式 \(p(s),s\in [x_l,x]\), 基于如下信息(我们还没有计算 \(f'(x)\) 的值) \[ \begin{aligned} p(x_l) &= f(x_l)\\ p'(x_l) &= f'(x_l)\\ p(x) &= f(x) \end{aligned} \] 得到 \(p(s)\) 的表达式 \[ \begin{aligned} p(s) &= p(x_l)\frac{(s-x)(s+x-2x_l)}{-(x_l-x)^2} + p'(x_l)\frac{(s-x_l)(s-x)}{x_l-x} + p(x)\frac{(s-x_l)^2}{(x-x_l)^2}\\ \end{aligned} \] 取新的试探点为极小值点 \[ x' = x_l +\frac12 \frac{(x_l-x)^2 p'(x_l)}{(p(x_l)-p(x))-(x_l-x)p'(x_l)} \]

当调整到右侧子区间\([x_l',x_r'] = [x,x_r]\),取两点二次插值多项式 \(p(s),s\in [0,1]\), 基于如下信息(我们已经计算了 \(f'(x)\) 的值) \[ \begin{aligned} p(x) &= f(x)\\ p'(x) &= f'(x)\\ p'(x_l) &= f'(x_l) \end{aligned} \] 得到 \(p(s)\) 的表达式 \[ \begin{aligned} p(s) &= p'(x_l)\frac{(s-x)^2}{2(x_l-x)} + p(x) + p'(x)\frac{(s-x)(s-2x_l+x)}{2(x-x_l)}\\ \end{aligned} \] 取新的试探点为极小值点 \[ x' = x - \frac{(x_l-x)p'(x)}{p'(x_l)-p'(x)} \]