正定二次目标函数

$$min f(x) = \frac 1 2 x^T G x + r^T x + \sigma$$

非线性最小二乘

$$\begin{align*} f(x) & = \frac 1 2 \sum_{i=1}^m r^2_i(x) \\ & = \frac 1 2 r(x) r^T(x) \end{align*}$$

• $r_i(x)$：通常为非线性函数
• $r(x) = (r_1(x), \cdots, r_n(x))^T$
• $x \in R^n, m \geq n$

则

$$\begin{align*} \nabla f(x) & = \sum_{i=1}^m \nabla r_i(x) r_i(x) \\ & = J(x)^T r(x) \\ \nabla^2 f(x) & = \sum_{i=1}^m \nabla r_i(x) r_i(x) + \sum_{i=1}^m r_i \nabla^2 r_i(x) \\ & = J(x)^T J(x) + \sum_{i=1}^m r_i(x) \nabla^2 r_i(x) \end{align*}$$

• $$J(x) = \begin{bmatrix} \frac {\partial r_1} {\partial x_1} & \frac {\partial r_1} {\partial x_2} & \cdots & \frac {\partial r_1} {\partial x_n} \\ \frac {\partial r_2} {\partial x_1} & \frac {\partial r_2} {\partial x_2} & \cdots & \frac {\partial r_2} {\partial x_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac {\partial r_m} {\partial x_1} & \frac {\partial r_m} {\partial x_2} & \cdots & \frac {\partial r_m} {\partial x_n} \end{bmatrix} = \begin{bmatrix} \nabla r_1(x)^T \\ \nabla r_2(x)^T \\ \vdots \\ \nabla r_m(x)^T \\ \end{bmatrix}$$ 为$r(x)$的Jacobi矩阵

Gauss-Newton法

Newton法中为简化计算，略去其Hesse矩阵中 $\sum_{i=1}^m r_i(x) \nabla^2 r_i(x)$项，即直接求解 方程组

$$J(x^{(k)})^T J(x^{(k)}) d = -J(x^{(k)})^T r(x^{(k)})$$

算法

同Newton法，仅求解Newton方程改为求解以上方程组

特点

• 实际问题中

  • 局部解$x^{ }$对应的目标函数值$f(x^{ })$接近0 时，$|r(x^{(k)})|$较小
  • 曲线$r_i(x)$接近直线， $\nabla^2 r_i(x) \approx 0$

  采用Gauss-Newton法效果较好，否则效果一般

• 矩阵$J(x^{(k)})^T J(x^{(k)})$是半正定矩阵，当Jacobi矩阵 列满秩时为正定矩阵，此时虽然$d^{(k)}$是下降方向，但仍需 类似修正牛顿法增加一维搜索策略保证目标函数值不上升

Levenberg-Marquardt方法

但$J(x^{(k)})$中各列线性相关、接近线性相关，则求解 Newton-Gauss方法中的方程组会出现困难，可以改为求解

$$(J(x^{(k)})^T J(x^{(k)}) + vI) d = -J(x^{(k)})^T r(x^{(k)})$$

• $v$：迭代过程中需要调整的参数，LM方法的关键即如何调整

定理1

• 若$d(v)$是以上方程组的解，则$|d(v)|^2$是$v$的连续下降 函数，且$v \rightarrow +\infty, |d(v)| \rightarrow 0$

• $J(x^{(k)})^T J(x^{(k)})$是对称半正定矩阵，则存在正交阵

  $$(P^{(k)})^T J(x^{(k)})^T J(x^{(k)}) P^{(k)} = \Lambda^{(k)}$$

• 则可以解出$|d(v)|^2$

• 增大$v$可以限制$|d^{(k)}|$，所以LM方法也被称为阻尼最小 二乘法

定理2

• 若$d(v)$是以上方程的解，则$d(v)$是$f(x)$在$x^{(k)}$处的 下降方向，且$v \rightarrow + \infty$时，$d(v)$的方向与 $-J(x^{(k)})^T r(x^{(k)})$方向一致

• 下降方向：$\nabla f(x^{(k)}) d(v) < 0$即可
• 方向一致：夹角余弦

• $v$充分大时，LM方法产生的搜索方向$d^{(k)}$和负梯度方向 一致

参数调整方法

使用梯度、近似Hesse矩阵定义二次函数

$$q(d) = f(x^{(k)}) + (J(x^{(k)})^T r(x^{(k)}))^T d + \frac 1 2 d^T (J(x^{(k)})^T J(x^{(k)})) d$$

• 其增量为

  $$\begin{align*} \Delta q^{(k)} & = q(d^{(k)}) - q(0) \\ & = (J(x^{(k)})^T r(x^{(k)}))^T d^{(k)} + \frac 1 2 (d^{(k)})^T (J(x^{(k)})^T J(x^{(k)})) d^{(k)} \end{align*}$$

• 目标函数增量

  $$\begin{align*} \Delta f^{(k)} & = f(x^{(k)} + d^{(k)}) - f(x^{(k)}) \\ & = f(x^{(k+1)}) - f(x^{(k)}) \end{align*}$$

• 定义$\gamma_k = \frac {\Delta f^{(k)}} {\Delta q^{(k)}}$

• $\gamma_k$接近1说明$\Delta f^{(k)}$接近$\Delta q^{(k)}$

  • 即$f(x^{(k)} + d^{(k+1)})$接近$q(d^{(k)})$
  • 即$f(x)$在$x^{(k)}$附近接近二次函数
  • 即使用Gauss-Newton方法求解最小二乘问题效果较好
  • 即LM方法求解时$v$参数应该较小

• $\gamma_k$接近0说明$\Delta f^{(k)}$与$\Delta q^{(k)}$ 近似程度不好

  • $d^{(k)}$不应取得过大，应减少$d^{(k)}$得模长
  • 应该增加参数$v$进行限制
  • 迭代方向趋近于负梯度方向

• $\gamma_k$适中时，认为参数$v$选取合适，不做调整

  • 临界值通常为0.25、0.75

算法

1. 初始点$x^{(1)}$、初始参数$v$（小值）、精度要求$\epsilon$ ，置k=k+1

2. 若$|J(x^{(k)})^T r(x^{(k)})| < \epsilon$，则停止计算， 得到问题解$x^{(k)}$，否则求解线性方程组

   $$(J(x^{(k)})^T J(x^{(k)}) + v_kI) d = -J(x^{(k)})^T r(x^{(k)})$$

   得到$d^{(k)}$

3. 置$x^{(k+1)} = x^{(k)} + d^{(k)}$，计算$\gamma_k$

4. 若

   • $\gamma < 0.25$，置$v_{k+1} = 4 v_k$
   • $\gamma > 0.75$，置$v_{k+1} = v_k / 2$
   • 否则置$v_{k+1} = v_k$

5. 置k=k+1，转2

Newton法

思想

• 若$x^{ * }$是无约束问题局部解，则有

  $$\nabla f(x^{ * }) = 0$$

  可求解此问题，得到无约束问题最优解

• 原始问题是非线性，考虑求解其线性逼近，在初始点$x^{(1)}$ 处泰勒展开

  $$\nabla f(x) \approx \nabla f(x^{(1)}) + \nabla^2 f(x^{(1)})(x - x^{(1)})$$

  解得

  $$x^{(2)} = x^{(1)} - (\nabla^2 f(x^{(1)}))^{-1} \nabla f(x^{(1)})$$

  作为$x^{ * }$的第二次近似

• 不断迭代，得到如下序列

  $$x^{(k+1)} = x^{(k)} + d^{(k)}$$

  • $d^{(k)}$：Newton方向，即以下方程解 $$\nabla^2 f(x^{(k)}) d = -\nabla f(x^{(k)})$$

算法

• 初始点 $x^{(1)}$、精度要求 $\epsilon$，置 $k=1$

• 考虑 $|\nabla f(x^{(k)})| \leq \epsilon$

  • 若满足，则停止计算，得到最优解 $x^{(k)}$

  • 否则求解如下方程，得到 $d^{(k)}$

    $$\nabla^2 f(x^{(k)}) d = -\nabla f(x^{(k)})$$

• 如下设置，并转2

  $$x^{(k+1)} = x^{(k)} + d^{(k)}, k = k+1$$

特点

• 优点

  • 产生点列 ${x^{k}}$ 若收敛，则具有二阶收敛速率
  • 具有二次终止性，事实上对正定二次函数，一步即可收敛

• 缺点

  • 可能会在某步迭代时目标函数值上升
  • 当初始点 $x^{(1)}$ 距离最优解 $x^{ * }$ 时，产生的点列 可能不收敛，或者收敛到鞍点
  • 需要计算 Hesse 矩阵
    • 计算量大
    • Hesse 矩阵可能不可逆，算法终止
    • Hesse 矩阵不正定，Newton 方向可能不是下降方向

阻尼/修正 Newton 法

• 克服 Newton 法目标函数值上升的缺点
• 一定程度上克服点列可能不收敛缺点

算法

• 初始点 $x^{(1)}$、精度要求 $\epsilon$，置 $k=1$

• 考虑 $|\nabla f(x^{(k)})| \leq \epsilon$

  • 若满足，停止计算，得到最优解 $x^{(k)}$
  • 否则求解如下方程得到 $d^{(k)}$ $$\nabla^2 f(x^{(k)}) d = -\nabla f(x^{(k)})$$

• 一维搜索，求解一维问题

  $$\arg\min_{\alpha} \phi(\alpha) = f(x^{(k)} + \alpha d^{(k)})$$

  得到 $\alpha_k$，如下设置并转2

  $$x^{(k+1)} = x^{(k)} + \alpha_k d^{(k)}, k = k+1$$

其他改进

• Newton 法、修正 Newton 法的改进方向
  • 结合最速下降方向修正迭代方向
  • Hesse 矩阵不正定情形下的替代

结合最速下降方向

• 将 Newton 方向和最速下降方向结合

• 设 $\theta_k$ 是 $$ 之间夹角，显然希望 $\theta < \frac \pi 2$

• 则置限制条件 $\eta$，取迭代方向

  $$d^{(k)} = \left \{ \begin{array}{l} d^{(k)}, & cos\theta_k \geq \eta \\ -\nabla f(x^{(k)}), & 其他 \end{array} \right.$$

Negative Curvature

• 当 Hesse 矩阵非正定时，选择负曲率下降方向 $d^{(k)}$（一定存在）

• Hesse 矩阵非正定时，一定存在负特征值、相应特征向量 $u$

  • 取负曲率下降方向作为迭代方向

    $$d^{(k)} = -sign(u^T \nabla f(x^{(k)})) u$$

  • $x^{(k)}$ 处负曲率方向 $d^{(k)}$ 满足

    $${d^{(k)}}^T \nabla^2 f(x^{(k)}) d^{(k)} < 0$$

修正 Hesse 矩阵

• 取迭代方向 $d^{(k)}$ 为以下方程的解

  $$(\nabla^2 f(x^{(k)}) + v_k I) d = -\nabla f(x^{k})$$

• $v_k$：大于 $\nabla^2 f(x^{(k)})$ 最大负特征值绝对值

综述

拟Newton法/变度量法：不需要求解Hesse矩阵，使用一阶导构造 二阶信息的近似矩阵

• 使用迭代过程中信息，创建近似矩阵$B^{(k)}$代替Hesse矩阵

• 用以下方程组替代Newton方程，其解$d^{(k)}$作为搜索方向

  $$B^{(k)} d = - \triangledown f(x^{(k)})$$

思想

• 考虑$\triangledown f(x)$在$x^{(k+1)}$处泰勒展开

  $$\triangledown f(x) \approx \triangledown f(x^{(k+1)}) + \triangledown^2 f(x^{(k+1)})(x - x^{(k+1)})$$

• 取$x = x^{(k)}$，有

  $$\begin{align*} \triangledown f(x^{(k+1)}) - \triangledown f(x^{(k)}) & \approx \triangledown^2 f(x^{(x+1)}) (x^{(k+1) } - x^{(k)}) \\ \triangledown^2 f(x^{k+1}) s^{(k)} & \approx y^{(k)} \end{align*}$$

  • $s^{(k)} = x^{(k+1)} - x^{(k)}$
  • $y^{(k)} = \triangledown f(x^{(k+1)}) - \triangledown f(x^{(k)})$

• 要求$B^{(k)}$近似$\triangledown^2 f(x^{(k)})$，带入并将 $\approx$改为$=$，得到拟Newton方程

  $$B^{(k+1)} s^{(k)} = y^{(k)}$$

  并假设$B^{(k)}$对称

• 拟Newton方程不能唯一确定$B^{(k+1)}$，需要附加条件，自然 的想法就是$B^{(k+1)}$可由$B^{(k)}$修正得到，即

  $$B^{(k+1)} = B^{(k)} + \Delta B^{(k)}$$

  且修正项$\Delta B^{(k)}$具有“简单形式”

Hesse矩阵修正

对称秩1修正

认为简单指矩阵秩小：即认为$\Delta B^{(k)}$秩为最小值1

• 设$\Delta B^{(k)} = u v^T$，带入有

  $$\begin{align*} y^{(k)} & = B^{(k+1)} s^{(k)} \\ & = B^{(k)} s^{(k)} + (v^T s^{(k)}) u \\ y^{(k)} - B^{(k)} s^{(k)} & = (v^T s^{(k)}) u \end{align*}$$
  • 这里有的书会设$\Delta B^{(k)} = \alpha u v^T$， 其实对向量没有必要
  • $v^T s^{(k)}$是数，所以$u$必然与共线，同理也没有必要 考虑系数，直接取相等即可
  • 而且系数不会影响最终结果

• 可取$u = y^{(k)} - B^{(k)} s{(k)}$，取$v$满足 $v^T s^{(k)} = 1$

• 由$B^{(k)}$的对称性、并希望$B^{(k+1)}$保持对称，需要 $u, v$共线，则有

  $$\begin{align*} v & = \lambda u = \lambda (y^{(k)} - B^{(k)} s^{(k)}) \\ 1 & = \lambda (y^{(k)} - B^{(k)} s^{(k)})^T s^{(k)} \end{align*}$$

• 得到$B^{(k)}$的对称秩1修正公式

  $$B^{(k+1)} = B^{(k)} + \frac {(y^{(k) - B^{(k)} s^{(k)}}) (y^{(k)} - B^{(k)} s^{(k)})^T} {(y^{(k)} - B^{(k)} s^{(k)})^T s^{(k)}}$$

算法

1. 初始点$x^{(1)}$、初始矩阵$B^{(1)} = I$、精度要求 $\epsilon$、置$k=1$

2. 若$|\triangledown f(x^{(k)})| \leq \epsilon$，停止计算 ，得到解$x^{(k)}$，否则求解以下方程得到$d^{(k)}$

   $$B^{(k)} d = -\triangle f(x^{(k)})$$

3. 一维搜索，求解

   $$\arg\min_{\alpha} \phi(\alpha)=f(x^{(k)} + \alpha d^{(k)})$$

   得到$\alpha_k$，置$x^{(k+1)}=x^{(k)} + \alpha_k d^{(k)}$

4. 修正$B^{(k)}$

   $$\begin{align*} s^{(k)} & = x^{(k+1)} - x^{(k)} \\ y^{(k)} & = \triangledown f(x^{(k+1)}) - \triangledown f(x^{(k)}) \\ B^{(k+1)} & = B^{(k)} + \frac {(y^{(k) - B^{(k)} s^{(k)}}) (y^{(k)} - B^{(k)} s^{(k)})^T} {(y^{(k)} - B^{(k)} s^{(k)})^T s^{(k)}} \end{align*}$$

5. 置$k = k+1$，转2

特点

• 缺点

  • 要求$(y^{(k)} - B^{(k)} s^{(k)})^T s^{(k)} \neq 0$， 否则无法继续计算

  • 不能保证正定性传递，只有 $(y^{(k)} - B^{(k)} s^{(k)})^T s^{(k)} > 0$才能保证 $B^{(k+1)}$也正定

  • 即使$(y^{(k)} - B^{(k)} s^{(k)})^T s^{(k)} > 0$， 也可能很小，容易产生较大的舍入误差

对称秩2修正

• 为克服秩1修正公式缺点，考虑$\Delta B^{(k)}$秩为2，设

  $$\Delta B^{(k)} = u^{(1)} (v^{(1)})^T + u^{(2)} (v^{(2)})^T$$

• 带入拟Newton方程有

  $$B^{(k)} s^{(k)} + ((v^{(1)})^T s^{(k)}) u^{(1)} + ((v^{(2)})^T s^{(k)}) u^{(2)} = y^{(k)}$$

• 类似的取

  $$\left \{ \begin{array}{l} u^{(1)} = y^{(k)} \\ (v^{(1)})^T s^{(k)} = 1 \end{array} \right.$$ $$\left \{ \begin{array}{l} u^{(2)} = -B^{(k)} s^{(k)} \\ (v^{(2)})^T s^{(k)} = 1 \end{array} \right.$$

• 同秩1公式保持对称性推导，得到对称秩2修正公式/BFGS公式

  $$B^{(k+1)} = B^{(k)} - \frac {B^{(k)} s^{(k)} (s^{(k)})^T B^{(k)}} {(s^{(k)})^T B^{(k)} s^{(k)}} + \frac {y^{(k)} (y^{(k)})^T} {(y^{(k)})^T s^{(k)}}$$

BFGS算法

类似同秩1修正算法，仅第4步使用对称秩2修正公式

Hesse逆修正

对称秩2修正

• 考虑直接构造近似于$(\triangledown^2 f(x^{(k)}))^{-1}$的 矩阵$H^{(k)}$

• 这样无需求解线性方程组，直接计算

  $$d^{(k)} = -H^{(k)} \triangledown f(x^{(k)})$$

• 相应拟Newton方程为

  $$H^{(k+1)} y^{(k)} = s^{(k)}$$

• 可得$H^{(k)}$的对称秩1修正公式

  $$H^{(k+1)} = H^{(k)} + \frac {(s^{(k)} - H^{(k)} y^{(k)}) (s^{(k)} - H^{(k)} y^{(k)})T} {(s^{(k)} - H^{(k)} y^{(k)})^T y^{(k)}}$$

• 可得$H^{(k)}$的对称秩2修正公式/DFP公式

  $$H^{(k+1)} = H^{(k)} - \frac {H^{(k)} y^{(k)} (y^{(k)})^T H^{(k)}} {(y^{(k)})^T H^{(k)} y^{(k)}} + \frac {s^{(k)} (s^{(k)})^T} {(s^{(k)})^T y^{(k)}}$$

DFP算法

类似BFGS算法，只是

• 使用$H^{(k)}$计算更新方向
• 使用$H^{(k)}$的对称秩2修正公式修正

• 对正定二次函数，BFGS算法和DFP算法效果相同
• 对一般可微（非正定二次函数），一般认为BFGS算法在收敛性质 、数值计算方面均由于DFP算法

Hesse逆的BFGS算法

• 考虑

  $$\begin{align*} B^{(k+1)} & = B^{(k)} + u^{(1)} (v^{(1)})^T + u^{(2)} (v^{(2)})^T \\ H^{(k+1)} & = (B^{(k+1)})^{-1} \\ & = (B^{(k)} + u^{(1)} (v^{(1)})^T + u^{(2)} (v^{(2)})^T)^{-1} \\ \end{align*}$$

• 两次利用Sherman-Morrison公式，可得

  $$H^{(k+1)} = (I - \frac {s^{(k)} (y^{(k)})^T} {(y^{(k)})^T s^{(k)}}) H^{(k)} (I - \frac {s^{(k)} (y^{(k)})^T} {(y^{(k)})^T s^{(k)}})^T + \frac {s^{(k)} (s^{(k)})^T} {(y^{(k)})^T s^{(k)}}$$

todo

• 还可以进一步展开

  $$H^{(k+1)} = H^{(k)} + (\frac 1 {(s^{(k)})^T y^{(k)}} + \frac {(y^{(k)})^T H^{(k)} y^{(k)}} {((s^{(k)})^T y^{(k)})^2}) s^{(k)} (s^{(k)})^T - \frac 1 {(s^{(k)})^T y^{(k)}} (H^{(k)} y^{(k)} (s^{(k)})^T + s^{(k)} (y^{(k)})^T H^{(k)})$$

变度量法的基本性质

算法的下降性

定理1

• 设$B^{(k)}$（$H^{(k)}$）是正定对称矩阵，且有 $(s^{(k)})^T y^{(k)} > 0$，则由BFGS（DFS）公式构造的 $B^{(k+1)}$（$H^{(k+1)}$）是正定对称的

• 考虑$B^{(k)}$对称正定，有 $B^{(k)} = (B^{(k)})^{1/2} (B^{(k)})^{1/2}$

• 带入利用柯西不等式即可证

• 中间插入正定矩阵的向量内积不等式也称为广义柯西不等式

定理2

• 若$d^{(k)}$v是下降方向，且一维搜索是精确的，设 $B^{(k)}$（$H^{(k)}$）是正定对称矩阵，则有BFGS（DFP） 公式构造的$B^{(k+1)}$（$H^{(k+1)}$）是正定对称的

• 精确一维搜索$(d^{(k)})^T \triangledown f(x^{(k+1)}) = 0$
• 则有$(s^{(k)})^T y^{(k)} > 0$

定理3

• 若用BFGS算法（DFP算法）求解无约束问题，设初始矩阵 $B^{(1)}$（$H^{(1)}$）是正定对称矩阵，且一维搜索是精确的 ，若$\triangledown f(x^{(k)}) \neq 0$，则产生搜索方向 $d^{(k)}$是下降方向

• 结合上2个结论，数学归纳法即可

总结

• 若每步迭代一维搜索精确，或满足$(s^{(k)})^T y^{(k)} > 0$

  • 停止在某一稳定点
  • 或产生严格递减的序列${f(x^{(k)})}$

• 若目标函数满足一定条件我，可以证明变度量法产生的点列 ${x^{(k)}}$收敛到极小点，且收敛速率超线性

搜索方向共轭性

• 用变度量法BFGS（DFP）算法求解正定二次函数

$$
min f(x) = \frac 1 2 x^T G x + r^T x + \sigma
$$

若一维搜索是精确的，假设已经进行了m次迭代，则

• 搜索方向$d^{(1)}, \cdots, d^{(m)}$是m个非零的G共轭方向

• 对于$j = 1, 2, \cdots, m$，有

$$
B^{(m+1)} s^{(j)} = y^{(j)}
(H^{(m+1)} y^{(j)} = s^{(j)})
$$

且$m = n$时有吧

$$
B^{(n+1)} = G(H^{(n+1)} = G^{-1})
$$

变度量法二次终止

• 若一维搜索是精确的，则变度量法（BFGS、DFP）具有二次终止

• 若$\triangle f(x^{(k)}) = 0, k \leq n$，则得到最优解 $x^{(k)}$

• 否则得到的搜索方向是共轭的，由扩展空间子定理， $x^{(n+1)}$是最优解