DFP算法

拟牛顿法中的DFP(Davidon-Fletcher-Powell)算法为：

DFP算法：初始化，取正定矩阵 $H _ 0$ ，初始点 $\boldsymbol x _ 0$ ， $k=0$ 。

令 $\boldsymbol d _ k=-H _ k\nabla f(\boldsymbol x _ k)$ .
优化 $\min _ {\alpha\geqslant0}f(\boldsymbol x _ {k}+\alpha \boldsymbol d _ k)$ 得到 $\alpha _ k,\boldsymbol x _ {k+1}$ 及 $\boldsymbol s _ k=\alpha _ k\boldsymbol d _ k$ .
令 $\boldsymbol y _ k=\nabla f(\boldsymbol x _ {k+1})-\nabla f(\boldsymbol x _ k)$ ，以及
$H _ {k+1}=H _ k-\frac{H _ k\boldsymbol y _ k\boldsymbol y _ k^\mathsf TH _ k}{\boldsymbol y _ k^\mathsf TH _ k\boldsymbol y _ k}+\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}.$ 更新 $k$ ，若未满足终止条件，回到步骤1。

Contents

1 DFP算法
2 个人推导
3 代数方法

令 $B _ k=H _ k^{-1}$ ，可以证明所有 $B _ k,H _ k$ 都是正定的（可见https://gaomj.cn/optimization-2/）。步骤三中的迭代公式中， $B _ {k+1}:=H _ {k+1}^{-1}$ 据称是如下优化问题的解：
$\begin{aligned} \min _ B{}&\quad \|B-B _ k\| _ W\\ \text{s.t. }&\quad B^\mathsf T=B,\\ &\quad B\boldsymbol s _ k=\boldsymbol y _ k. \end{aligned}$ 其中范数 $\|\cdot\| _ W$ 为加权Frobenius范数，定义为 $\|A\| _ W:=\|W^\frac12AW^\frac12\| _ F$ ，这里 $\|\cdot\| _ F$ 则为熟知的Frobenius范数，定义为 $\|A\|^2 _ F=\mathsf{tr}(A^\mathsf TA)=\sum _ {i,j}a _ {ij}^2$ ；这里我们选择权重矩阵 $W$ 使得 $W\boldsymbol y _ k=\boldsymbol s _ k$ 且 $W\succ 0$ ，这是可以做到的，因为可以证明总有 $\boldsymbol y _ k^\mathsf T\boldsymbol s _ k>0$ ，即夹角是锐角，我们考虑一个旋转变换（可以参见https://gaomj.cn/linearalgebra6/，属于正交变换的基本结果）的矩阵 $U$ 使得旋转后 $U\boldsymbol y,U\boldsymbol s$ 的分量都是正数，这样就有一个正对角线的对角矩阵 $\Lambda$ 使得 $\Lambda U\boldsymbol y=U\boldsymbol s$ ，即 $U^\mathsf T\Lambda U\boldsymbol y=\boldsymbol s$ ，就得到了正定矩阵 $U^\mathsf T\Lambda U=:W$ 。

下面给出该公式的个人推导。在最后还将给出查阅得到的其他人给出的推导。

个人推导

思路是采取微积分的方法（Lagrange乘数法）证明，但是也要涉及大量的矩阵代数的运算，十分地不容易。（“容错”也不高，稍微剑走偏锋就很容易证不了）

第一步

先做一个代换， $B\leftarrow W^\frac12BW^\frac12$ ， $B _ k\leftarrow W^\frac12B _ kW^\frac12$ ，代换是可逆的。现在问题变成了最小化 $\|B-B _ k\| _ F$ 使得 $B$ 是对称矩阵且 $W^{-\frac12}BW^{-\frac12}\boldsymbol s _ k=\boldsymbol y _ k$ ，即 $BW^{-\frac12}\boldsymbol s _ k=W^\frac12\boldsymbol y _ k$ ，然后再做一个代换： $\boldsymbol s _ k\leftarrow W^{-\frac12}\boldsymbol s _ k$ ， $y _ k\leftarrow W^\frac12\boldsymbol y _ k$ ，这个代换也是可逆的。注意到 $W$ 的选择已经使得 $W\boldsymbol y _ k=\boldsymbol s _ k$ ，故 $W^\frac12\boldsymbol y _ k=W^{-\frac12}\boldsymbol s _ k$ ，所以代换后的 $\boldsymbol s _ k,\boldsymbol y _ k$ 是相同的，现在问题变成了
$\begin{aligned} \min _ B{}&\quad \|B-B _ k\| _ F\\ \text{s.t. }&\quad B^\mathsf T=B,\\ &\quad B\boldsymbol s _ k=\boldsymbol s _ k. \end{aligned}$ 这个问题更为简单，我们把范数变为了普通的Frobenius范数且 $\boldsymbol y _ k$ 变为了 $\boldsymbol s _ k$ 。

考虑一个“更大的”优化问题：
$\begin{aligned} \min _ B{}&\quad \|B-B _ k\| _ F\\ \text{s.t. }&\quad \boldsymbol s _ k^\mathsf TB=\boldsymbol s _ k^\mathsf T,\\ &\quad B\boldsymbol s _ k=\boldsymbol s _ k. \end{aligned}$ 显然前面的问题是这个问题的子问题，如果这个问题的最优解 $B$ 是对称矩阵，那么我们就知道前面的子问题的最优解也是这个矩阵 $B$ 。

注：如果考虑另一个“更大的”优化问题： $\min _ B\|B-B _ k\| _ F$ 使得 $B\boldsymbol s _ k=\boldsymbol s _ k$ ，即去掉 $B$ 对称这一限制，按同样的思路如果解出一个对称的 $B$ 那么它就是原问题的解，可惜的是，最终解出来的 $B$ 不是对称矩阵。

第二步

把这个矩阵优化问题变为向量的问题。

将矩阵拉直为向量后要用到矩阵的Kronecker积的性质，定义和部分性质可参见另一篇矩阵梯度下的文章。下面要用到的有 $(A\otimes B)^\mathsf T=A^\mathsf T\otimes B^\mathsf T$ ， $AC\otimes BD=(A\otimes B)(C\otimes D)$ ， $(AXB)^\mathsf V=(B^\mathsf T\otimes A)\, X^\mathsf V$ ， $k\otimes A=A\otimes k=kA$ 。矩阵的Frobenius范数恰为矩阵拉直后的向量的范数。所以为方便用微积分的方法解决问题，将上面的问题改写成
$\begin{aligned} \min _ B{}&\quad\frac12\|B^\mathsf V-B _ k^\mathsf V\|^2,\\ \text{s.t. }&\quad (I _ n\otimes\boldsymbol s _ k^\mathsf T)\, B^\mathsf V=\boldsymbol s _ k,\\ &\quad (\boldsymbol s _ k^\mathsf T\otimes I _ n)\, B^\mathsf V=\boldsymbol s _ k. \end{aligned}$ 用Lagrange乘数法解决。

第三步

引入Lagrange函数
$\frac12\|B^\mathsf V-B _ k^\mathsf V\|^2-\boldsymbol \lambda^\mathsf T[(\boldsymbol s _ k^\mathsf T\otimes I _ n)\, B^\mathsf V-\boldsymbol s _ k]-\boldsymbol \mu^\mathsf T[(I _ n\otimes\boldsymbol s _ k^\mathsf T)\, B^\mathsf V-\boldsymbol s _ k].$
对 $B^\mathsf V$ 求梯度并置零得到
$B^\mathsf V-B _ k^\mathsf V=(\boldsymbol s _ k\otimes I _ n)\, \boldsymbol \lambda+(I _ n\otimes\boldsymbol s _ k)\, \boldsymbol \mu.$ 注意到约束条件，左乘 $(\boldsymbol s _ k^\mathsf T\otimes I _ n)$ 或 $(I _ n\otimes\boldsymbol s _ k^\mathsf T)$ 分别得到
$\begin{aligned} \boldsymbol s _ k-(\boldsymbol s _ k^\mathsf T\otimes I _ n)\, B _ k^\mathsf V&=(\boldsymbol s _ k^\mathsf T\boldsymbol s _ k\otimes I _ n)\, \boldsymbol \lambda+(\boldsymbol s _ k^\mathsf T\otimes \boldsymbol s _ k)\boldsymbol \mu,\\ \boldsymbol s _ k-(I _ n\otimes\boldsymbol s _ k^\mathsf T)\, B _ k^\mathsf V&=(\boldsymbol s _ k\otimes \boldsymbol s _ k^\mathsf T)\, \boldsymbol \lambda+(I _ n\otimes \boldsymbol s _ k^\mathsf T\boldsymbol s _ k)\boldsymbol \mu. \end{aligned}$ 这里注意到Kronecker积的另一个性质：
$\boldsymbol s _ k\otimes \boldsymbol s _ k^\mathsf T=\boldsymbol s _ k^\mathsf T\otimes \boldsymbol s _ k=\boldsymbol s _ k\boldsymbol s _ k^\mathsf T.$ 原式化简为
$\begin{aligned} \boldsymbol s _ k-(\boldsymbol s _ k^\mathsf T\otimes I _ n)\, B _ k^\mathsf V&=\boldsymbol s _ k^\mathsf T\boldsymbol s _ k\boldsymbol \lambda+\boldsymbol s _ k\boldsymbol s _ k^\mathsf T\boldsymbol \mu,\\ \boldsymbol s _ k-(I _ n\otimes\boldsymbol s _ k^\mathsf T)\, B _ k^\mathsf V&=\boldsymbol s _ k\boldsymbol s _ k^\mathsf T\boldsymbol \lambda+\boldsymbol s _ k^\mathsf T\boldsymbol s _ k\boldsymbol \mu. \end{aligned}$

第四步

考虑用上式消去 $B^\mathsf V-B _ k^\mathsf V=(\boldsymbol s _ k\otimes I _ n)\, \boldsymbol \lambda+(I _ n\otimes\boldsymbol s _ k)\, \boldsymbol \mu$ 中的 $\boldsymbol \lambda,\boldsymbol \mu$ 。首先求 $\boldsymbol \lambda$ ，
$\begin{gathered} \boldsymbol \lambda=\frac{1}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}[\boldsymbol s _ k-(\boldsymbol s _ k^\mathsf T\otimes I _ n)\, B _ k^\mathsf V-\boldsymbol s _ k\boldsymbol s _ k^\mathsf T\boldsymbol \mu],\\ (\boldsymbol s _ k\otimes I _ n)\, \boldsymbol \lambda=\frac{1}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}[(1-\boldsymbol s _ k^\mathsf T\boldsymbol \mu)(\boldsymbol s _ k\otimes\boldsymbol s _ k)-(\boldsymbol s _ k\boldsymbol s _ k^\mathsf T\otimes I _ n)\, B _ k^\mathsf V]. \end{gathered}$ $\boldsymbol \mu$ 同样可得
$(I _ n\otimes\boldsymbol s _ k)\, \boldsymbol \mu=\frac{1}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}[(1-\boldsymbol s _ k^\mathsf T\boldsymbol \lambda)(\boldsymbol s _ k\otimes\boldsymbol s _ k)-(I _ n\otimes\boldsymbol s _ k\boldsymbol s _ k^\mathsf T)\, B _ k^\mathsf V].$ 这样两式相加就得到
$B^\mathsf V-B _ k^\mathsf V=\frac{1}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}[(2-\boldsymbol s _ k^\mathsf T(\boldsymbol\lambda+\boldsymbol \mu))(\boldsymbol s _ k\otimes\boldsymbol s _ k)-(\boldsymbol s _ k\boldsymbol s _ k^\mathsf T\otimes I _ n)\, B _ k^\mathsf V-(I _ n\otimes\boldsymbol s _ k\boldsymbol s _ k^\mathsf T)\, B _ k^\mathsf V].$
只需要再计算 $2-\boldsymbol s _ k^\mathsf T(\boldsymbol\lambda+\boldsymbol \mu)$ ，只需要用 $\boldsymbol s _ k-(\boldsymbol s _ k^\mathsf T\otimes I _ n)\, B _ k^\mathsf V=\boldsymbol s _ k^\mathsf T\boldsymbol s _ k\boldsymbol \lambda+\boldsymbol s _ k\boldsymbol s _ k^\mathsf T\boldsymbol \mu$ 左乘 $\boldsymbol s _ k^\mathsf T$ 计算，于是
$1-\frac{1}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}(\boldsymbol s _ k^\mathsf T\otimes\boldsymbol s _ k^\mathsf T)\, B _ k^\mathsf V=\boldsymbol s _ k^\mathsf T\boldsymbol \lambda+\boldsymbol s _ k^\mathsf T\boldsymbol \mu.$ 这样 $2-\boldsymbol s _ k^\mathsf T(\boldsymbol\lambda+\boldsymbol \mu)=1+\frac{1}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}(\boldsymbol s _ k^\mathsf T\otimes\boldsymbol s _ k^\mathsf T)\, B _ k^\mathsf V$ ，再注意到 $\boldsymbol s _ k^\mathsf TB _ k\boldsymbol s _ k=(\boldsymbol s _ k^\mathsf TB _ k\boldsymbol s _ k)^\mathsf V=(\boldsymbol s _ k^\mathsf T\otimes\boldsymbol s _ k^\mathsf T)\, B _ k^\mathsf V$ ，

代入 $B^\mathsf V$ 的式子
$B^\mathsf V=B _ k^\mathsf V+\frac{1}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}\Big[\Big(1+\frac{\boldsymbol s _ k^\mathsf TB _ k\boldsymbol s _ k}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}\Big)(\boldsymbol s _ k\otimes\boldsymbol s _ k)-(\boldsymbol s _ k\boldsymbol s _ k^\mathsf T\otimes I _ n)\, B _ k^\mathsf V-(I _ n\otimes\boldsymbol s _ k\boldsymbol s _ k^\mathsf T)\, B _ k^\mathsf V\Big].$

第五步(B)

首先，将向量变回矩阵，这样做的依据是两个相同维度的矩阵的拉直向量相同那么这两个矩阵相等。
$B _ {k+1}=B _ k+\frac{1}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}\Big[\Big(1+\frac{\boldsymbol s _ k^\mathsf TB _ k\boldsymbol s _ k}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}\Big)\boldsymbol s _ k\boldsymbol s _ k^\mathsf T-B _ k\boldsymbol s _ k\boldsymbol s _ k^\mathsf T-\boldsymbol s _ k\boldsymbol s _ k^\mathsf TB _ k\Big].$
将含 $B _ k$ 的项写在一起，
$\begin{aligned} B _ {k+1}&=\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}+B _ k-B _ k\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}-\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}B _ k+\frac{\boldsymbol s _ k^\mathsf TB _ k\boldsymbol s _ k}{(\boldsymbol s _ k^\mathsf T\boldsymbol s _ k)^2}\boldsymbol s _ k\boldsymbol s _ k^\mathsf T\\ &=\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}+B _ k-B _ k\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}-\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}B _ k+\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}B _ k\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}\\ &=\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}+\Big(I _ n-\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}\Big)B _ k\Big(I _ n-\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}\Big). \end{aligned}$
其次，将第一步的变换换回去， $B\rightarrow W^\frac12BW^\frac12$ ， $B _ k\rightarrow W^\frac12B _ kW^\frac12$ ， $\boldsymbol s _ k\rightarrow W^{-\frac12}\boldsymbol s _ k=W^\frac12\boldsymbol y _ k$ ，
$\begin{aligned} W^\frac12B _ {k+1}W^\frac12&=\frac{W^\frac12\boldsymbol y _ k\boldsymbol y _ k^\mathsf TW^\frac12}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}\\ &\quad{}+\Big(W^{\frac12}W^{-\frac12}-\frac{W^\frac12\boldsymbol y _ k\boldsymbol s _ k^\mathsf TW^{-\frac12}}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}\Big)W^\frac12B _ kW^\frac12\Big(W^{-\frac12}W^\frac12-\frac{W^{-\frac12}\boldsymbol s _ k\boldsymbol y _ k^\mathsf TW^\frac12}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}\Big). \end{aligned}$ 求得最初的优化问题结果为（由于是对称矩阵，这就是最终结果）
$B _ {k+1}=\frac{\boldsymbol y _ k\boldsymbol y _ k^\mathsf T}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}+\Big(I _ n-\frac{\boldsymbol y _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}\Big)B _ k\Big(I _ n-\frac{\boldsymbol s _ k\boldsymbol y _ k^\mathsf T}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}\Big).$ 此即DFP算法对 $B$ 的更新公式。

第五步(H)

我们也可以用 $B _ k=H _ k^{-1},B _ {k+1}=H _ {k+1}^{-1}$ 得到关于 $H$ 的更新公式。类似Sherman–Morrison公式的推导（可见https://gaomj.cn/sherman-morrison/），这个公式形式为：
$(A+\boldsymbol u\boldsymbol v^\mathsf T)^{-1}=A^{-1}-\frac{A^{-1}\boldsymbol u\boldsymbol v^\mathsf TA^{-1}}{1+\boldsymbol v^\mathsf TA^{-1}\boldsymbol u}.$
我们再次拆开 $B _ {k+1}$ 的括号（或者跳过第五步( $B$ )，直接利用第五步( $B$ )开头的式子，就不用前面合并后这里又拆开了）
$\begin{aligned} B _ {k+1}&=B _ k+\frac{\boldsymbol y _ k\boldsymbol y _ k^\mathsf T}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}-\frac{\boldsymbol y _ k\boldsymbol s _ k^\mathsf TB _ k}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}-\frac{B _ k\boldsymbol s _ k\boldsymbol y _ k^\mathsf T}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}+\frac{\boldsymbol y _ k\boldsymbol s _ k^\mathsf TB _ k\boldsymbol s _ k\boldsymbol y _ k^\mathsf T}{(\boldsymbol y _ k^\mathsf T\boldsymbol s _ k)^2}\\ &=B _ k+\frac1{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}\Big(\boldsymbol y _ k\boldsymbol y _ k^\mathsf T-\boldsymbol y _ k\boldsymbol s _ k^\mathsf TB _ k-B _ k\boldsymbol s _ k\boldsymbol y _ k^\mathsf T+\frac{\boldsymbol y _ k\boldsymbol s _ k^\mathsf TB _ k\boldsymbol s _ k\boldsymbol y _ k^\mathsf T}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}\Big)\\ &=B _ k+\frac1{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}\Big[\Big(1+\frac{\boldsymbol s _ k^\mathsf TB _ k\boldsymbol s _ k}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}\Big)\boldsymbol y _ k\boldsymbol y _ k^\mathsf T-\boldsymbol y _ k\boldsymbol s _ k^\mathsf TB _ k-B _ k\boldsymbol s _ k\boldsymbol y _ k^\mathsf T\Big]\\ &=B _ k+\frac1{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}\begin{bmatrix}\boldsymbol y _ k& B _ k\boldsymbol s _ k\end{bmatrix} \begin{bmatrix} 1+\frac{\boldsymbol s _ k^\mathsf TB _ k\boldsymbol s _ k}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}&-1\\ -1&0 \end{bmatrix} \begin{bmatrix} \boldsymbol y _ k^\mathsf T\\ \boldsymbol s _ k^\mathsf TB _ k \end{bmatrix} . \end{aligned}$ 要对右边求逆，记
$U=\begin{bmatrix} \boldsymbol y _ k&B _ k\boldsymbol s _ k \end{bmatrix},\quad V=\frac1{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}\begin{bmatrix} 1+\frac{\boldsymbol s _ k^\mathsf TB _ k\boldsymbol s _ k}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}&-1\\ -1&0 \end{bmatrix}.$ 我们就要求 $B _ k+UVU^\mathsf T$ 的逆，考虑矩阵方程
$\begin{bmatrix} B _ k& U\\ U^\mathsf T& -V^{-1} \end{bmatrix} \begin{bmatrix} X\\ Y \end{bmatrix} =\begin{bmatrix} I _ n\\ 0 \end{bmatrix}\implies \begin{cases} B _ kX+UY=I _ n,\\ U^\mathsf TX-V^{-1}Y=0. \end{cases}$ 消去 $Y$ 得
$B _ kX+UVU^\mathsf TX=I _ n\implies X=(B _ k+UVU^\mathsf T)^{-1}.$ 就是说 $X$ 即为我们要求的。上面方程组的第一条式子又表明 $X=B _ k^{-1}(I _ n-UY)$ ，代入方程组的第二条式子， $U^\mathsf TX=V^{-1}Y$ 得 $U^\mathsf TB _ k^{-1}(I _ n-UY)=V^{-1}Y$ ，由此又可以求出 $Y$ ，
$U^\mathsf TB _ k^{-1}=(V^{-1}+U^\mathsf TB _ k^{-1}U)Y\implies Y=(V^{-1}+U^\mathsf TB _ k^{-1}U)^{-1}U^\mathsf TB _ k^{-1}.$ 这个 $Y$ 可稍作计算，
$\begin{aligned} V^{-1}+U^\mathsf TB _ k^{-1}U&=\boldsymbol y _ k^\mathsf T\boldsymbol s _ k\begin{bmatrix} 0&-1\\ -1&-1-\frac{\boldsymbol s _ k^\mathsf TB _ k\boldsymbol s _ k}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k} \end{bmatrix}+\begin{bmatrix} \boldsymbol y _ k^\mathsf TB _ k^{-1}\boldsymbol y _ k& \boldsymbol y _ k^\mathsf T\boldsymbol s _ k\\ \boldsymbol s _ k^\mathsf T\boldsymbol y _ k&\boldsymbol s _ k^\mathsf TB _ k^{-1}\boldsymbol s _ k \end{bmatrix}\\ &=\begin{bmatrix} \boldsymbol y _ k^\mathsf TB _ k^{-1}\boldsymbol y _ k& 0\\ 0&-\boldsymbol y _ k^\mathsf T\boldsymbol s _ k \end{bmatrix}, \end{aligned}$ 从而
$Y=\begin{bmatrix} 1\mathbin/\boldsymbol y _ k^\mathsf TB _ k^{-1}\boldsymbol y _ k& 0\\ 0&-1\mathbin/\boldsymbol y _ k^\mathsf T\boldsymbol s _ k \end{bmatrix} \begin{bmatrix} \boldsymbol y _ k^\mathsf TB _ k^{-1}\\ \boldsymbol s _ k^\mathsf T \end{bmatrix}.$ 将此 $Y$ 又代回原方程组的第一条 $B _ kX+UY=I _ n$ ，就得到
$\begin{aligned} (B _ k+UVU^\mathsf T)^{-1}=X&=B _ k^{-1}-B _ k^{-1}UY\\ &=B _ k^{-1}-\begin{bmatrix} B _ k^{-1}\boldsymbol y _ k&\boldsymbol s _ k \end{bmatrix} \begin{bmatrix} 1\mathbin/\boldsymbol y _ k^\mathsf TB _ k^{-1}\boldsymbol y _ k& 0\\ 0&-1\mathbin/\boldsymbol y _ k^\mathsf T\boldsymbol s _ k \end{bmatrix} \begin{bmatrix} \boldsymbol y _ k^\mathsf TB _ k^{-1}\\ \boldsymbol s _ k^\mathsf T \end{bmatrix}\\ B _ {k+1}^{-1}&=B _ k^{-1}-\frac{B _ k^{-1}\boldsymbol y _ k\boldsymbol y _ k^\mathsf TB _ k^{-1}}{\boldsymbol y _ k^\mathsf TB _ k^{-1}\boldsymbol y _ k}+\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}. \end{aligned}$ 这正是
$H _ {k+1}=H _ k-\frac{H _ k\boldsymbol y _ k\boldsymbol y _ k^\mathsf TH _ k}{\boldsymbol y _ k^\mathsf TH _ k\boldsymbol y _ k}+\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol y _ k^\mathsf T\boldsymbol s _ k}.$

代数方法

最后给出查阅得到的其他人给出的方法。这个方法用几乎纯代数的方式解前面提到的矩阵优化问题，范围对应上文第二到第四步。
$\begin{aligned} \min _ B{}&\quad \|B-B _ k\| _ F\\ \text{s.t. }&\quad B^\mathsf T=B,\\ &\quad B\boldsymbol s _ k=\boldsymbol s _ k. \end{aligned}$ 注意到第二条限制给出了 $B$ 的一个特征向量和特征值 $1$ 。设单位向量 $\boldsymbol u=\frac{\boldsymbol s _ k}{\|\boldsymbol s _ k\|}$ ，设正交矩阵 $U$ 的第一列为 $\boldsymbol u$ ， $U=\begin{pmatrix}\boldsymbol u& U _ \perp\end{pmatrix}$ ，则 $\boldsymbol u^\mathsf TB\boldsymbol u=1$ ， $U _ \perp^\mathsf TB\boldsymbol u=\boldsymbol 0$ 。注意到 $\|U^\mathsf TAU\| _ F^2=\mathsf{tr}(U^\mathsf TA^\mathsf TAU)=\mathsf{tr}(AA^\mathsf T)=\|A\| _ F^2$ ，
$\begin{aligned} \|B-B _ k\| _ F^2&=\|U^\mathsf TBU-U^\mathsf TB _ kU\| _ F^2\\ &=\left\|\begin{bmatrix} 1&\boldsymbol 0^\mathsf T\\ \boldsymbol 0& U _ \perp^\mathsf TBU _ \perp \end{bmatrix}- \begin{bmatrix} \boldsymbol u^\mathsf TB _ k\boldsymbol u&\boldsymbol u^\mathsf TB _ kU _ \perp\\ U _ \perp^\mathsf TB _ k\boldsymbol u&U _ \perp^\mathsf TB _ kU _ \perp \end{bmatrix}\right\| _ F^2\\ &=\left\| \begin{bmatrix} \boldsymbol u^\mathsf TB _ k\boldsymbol u-1&\boldsymbol u^\mathsf TB _ kU _ \perp\\ U _ \perp^\mathsf TB _ k\boldsymbol u&U _ \perp^\mathsf TB _ kU _ \perp -U _ \perp^\mathsf TBU _ \perp\end{bmatrix}\right\| _ F^2\\ &=(\boldsymbol u^\mathsf TB _ k\boldsymbol u-1)^2+\|\boldsymbol u^\mathsf TB _ kU _ \perp\| _ F^2+\|U _ \perp^\mathsf TB _ k\boldsymbol u\| _ F^2+\|U _ \perp^\mathsf TB _ kU _ \perp -U _ \perp^\mathsf TBU _ \perp\| _ F^2. \end{aligned}$ 于是 $B$ 使得 $\|B-B _ k\| _ F$ 最小化就是使得 $\|U _ \perp^\mathsf TB _ kU _ \perp -U _ \perp^\mathsf TBU _ \perp\| _ F$ 最小化，故 $B _ {k+1}$ 应满足
$U _ \perp^\mathsf TB _ {k+1}U _ \perp=U _ \perp^\mathsf TB _ kU _ \perp.$ 故
$\begin{aligned} B _ {k+1}&=U\begin{bmatrix} 1&\boldsymbol 0^\mathsf T\\ \boldsymbol 0& U _ \perp^\mathsf TB _ {k+1}U _ \perp \end{bmatrix}U^\mathsf T=U\begin{bmatrix} 1&\boldsymbol 0^\mathsf T\\ \boldsymbol 0& U _ \perp^\mathsf TB _ {k}U _ \perp \end{bmatrix}U^\mathsf T\\ &=\boldsymbol u\boldsymbol u^\mathsf T+U _ \perp U _ \perp^\mathsf TB _ kU _ \perp U _ \perp^\mathsf T\\ &=\boldsymbol u\boldsymbol u^\mathsf T+(I _ n-\boldsymbol u\boldsymbol u^\mathsf T)B _ k(I _ n-\boldsymbol u\boldsymbol u^\mathsf T)\\ &=\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}+\Big(I _ n-\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}\Big)B _ k\Big(I _ n-\frac{\boldsymbol s _ k\boldsymbol s _ k^\mathsf T}{\boldsymbol s _ k^\mathsf T\boldsymbol s _ k}\Big). \end{aligned}$ 再次得到了这个结果。

这个方法更为简单，前面的方法用Lagrange乘数，需要用复杂的代数运算以求解出 $\boldsymbol \lambda,\boldsymbol \mu$ ，稍有偏差便求解不得，本人推导是不少波折。这个一比真是简单到令人惊叹。

第五步本人对逆的求解没有直接用到Sherman-Morrison公式，而是直接构造矩阵方程解出来。按照各资料的说法，由于要求逆的矩阵是一个矩阵的二秩矫正(rank two correction)，所以可以用两次Sherman-Morrison公式解出来。感觉操作上可能不会很简单，就不再尝试这种做法。

DFP(Davidon-Fletcher-Powell)算法的个人推导

DFP算法

个人推导

第一步

第二步

第三步

第四步

第五步(B)

第五步(H)

代数方法

评论

发表回复取消回复

DFP(Davidon-Fletcher-Powell)算法的个人推导

DFP算法

个人推导

第一步

第二步

第三步

第四步

第五步(B)

第五步(H)

代数方法

评论

发表回复 取消回复

发表回复取消回复