牛顿法与拟牛顿法

牛顿法&拟牛顿法

设无约束优化问题：
$$\min f(x),x\in R^n$$

1 牛顿法

基本思想，通过泰勒二阶展开，通过对泰勒展开求导，并令其等于0，从而求得极小值。

将$f(x)$在$x_k$处进行泰勒展开：
$$f(x)\approx f(x_k)+\mathrm{\Delta }f(x_k)(x-x_k)+\frac{1}{2}{\left(x-x_k\right)}^T{\mathrm{\Delta }}^{2}f(x_k)(x-x_k)+0(||x-x_k|{|}^2)$$

对泰勒式及进行求导：
$$\nabla f(x)\approx \nabla f(x_k)+{\nabla }^{2}f(x_k)(x-x_k)$$

令上次为0，便可求得极小值位置：
$$\nabla f(x_k)=-{\nabla }^{2}f(x_k)(x-x_k)$$

得：
$$x=x_k-{\nabla }^{2}f(x_k{)}^{-1}\nabla f(x_k)$$
我们便可得下降到最小值的方向：$-{\nabla }^{2}f(x_k{)}^{-1}\nabla f(x_k)$

牛顿法流程：
设定初始点：$x_0$，迭代步数$k=0$
开始迭代：
1）计算梯度，Hessian矩阵，从而得到牛顿迭代方向：
$$\nabla f(x_k{)}^2,f(x_k),d_k=-{\nabla }^{2}f(x_k{)}^{-1}\nabla f(x_k)$$

2）判断当前迭代步长是否足够小或达到迭代次数，是则停止，否，继续下一步

3）更新：
$$\begin{array}{l}x=x_k+d_k\\ k++\end{array}$$
转步骤1）

Tips:
此外，牛顿法，是局部收敛的。因为是基于求解极值的方法，因而会局限于局部极小值，因而要求初始点选取在全局最优解附近。
牛顿法有一次收敛性。但由于${\nabla }^{2}f(x_k)$不能保证一定正定，会出现秩为0的情况，此时，牛顿法会失效。因而基于修正hessian矩阵的算法出现，如LM算法。

2 拟牛顿法

拟牛顿法是一种基于牛顿法的变形，提供一种思想，包含多种解决方案。
基本思想：牛顿法虽收敛速度快，但要计算Hessian矩阵及其逆矩阵，且Hessian矩阵不一定正定，为了减少计算量，拟牛顿法被提出，其基本思想是使用目标函数的一阶信息近似牛顿法的Hessian矩阵。

拟牛顿法的条件：使用k时刻的一阶二阶近似出$k+1$时刻的二阶Hessian矩阵。
我们将目标函数$f(x)$在$x_{k+1}$处进行泰勒展开：
$$f(x)\approx f(x_{k+1})+\mathrm{\Delta }f(x_{k+1})(x-x_{k+1})+\frac{1}{2}{\left(x-x_{k+1}\right)}^T{\mathrm{\Delta }}^{2}f(x_{k+1})(x-x_{k+1})+0(||x-x_{k+!}|{|}^2)$$

对泰勒式及进行求导：
$$\nabla f(x)\approx \nabla f(x_{k+1})+{\nabla }^{2}f(x_{k+1})(x-x_{k+1})$$

令$x=x_k,g_k=\nabla f(x_{k+1})$，则得：
$$g_k-g_{k+1}\approx {\nabla }^{2}f(x_{k+1})(x_k-x_{k+1})$$$$x_k-x_{k+1}\approx {\nabla }^{2}f(x_{k+1})(g_k-g_{k+1})$$

令$q_k=g_k-g_{k+1},p_k=x_k-x_{k+1},H_{k+1}={\left({\nabla }^{2}f(x_{k+1})\right)}^{-1}$，则得拟牛顿法的条件：
$$p_k\approx H_{k+1}{q}_k$$

Tips：
此处H表示的是Hessian矩阵的逆。
关于这个公式的理解，可以将x想为时间t，即，两帧(两次迭代之间$k->k+1$)的速度变化，等于$k+1$时刻的加速度*两帧之间的时差。
此外，如果我们对目标函数在$x_k$处展开，再将$x_{k+1}$代入，可以得到 。这里，引入思考，泰勒展开式只保证在展开点的一定邻域内有效，当$x_k$与$x_{k+1}$相差较大时，以上拟牛顿法的条件是否也有一定误差。

牛顿法是使$H_{k+1}=H_k+\Delta H_k$来拟合 $H_{k+1}$。

2.1 对称秩1校正

这个思想是，保证 的半正定及对称性，我们令：
$$\Delta H_k=v_k{v}_k^T，v_k\ne 0$$

由于$v_k{v}_k^T$使得$\Delta H_k$的秩为1，故而此方法名为秩1校正。
代入拟牛顿法条件
$$p_k\approx H_{k+1}{q}_k$$$$\begin{array}{l}{p}_k\approx H_{k+1}{q}_k=\left(H_k+\Delta H_k\right)q_k\\ p_k=\left(H_k+v_k{v}_k^T\right)q_k=H_k{q}_k+v_k{v}_k^T{q}_k\end{array}$$

从而得：$v_k{v}_k^T{q}_k=p_k-H_k{q}_k$
由于$q_k$为向量，上式可以看成$Ax=b$，$A$为对称对阵，求$A$的问题。
等式两侧都乘以${\left(p_k-H_k{q}_k\right)}^T$，得：
$$v_k{v}_k^T{q}_k{\left(p_k-H_k{q}_k\right)}^T=\left(p_k-H_k{q}_k\right){\left(p_k-H_k{q}_k\right)}^T$$
进而得：
$$v_k{v}_k^T=\frac{\left(p_k-H_k{q}_k\right){\left(p_k-H_k{q}_k\right)}^T}{q_k{\left(p_k-H_k{q}_k\right)}^T}$$

因而我们可以发现，秩1校正成立的条件 $q_k{\left(p_k-H_k{q}_k\right)}^T>0$

2.2 DFP

基本思想，我们可以从秩1校正中发现$\Delta H_k$的是由$p_k$和$H_k{q}_k$组成的，因而我么可以设计以下方法来完成校正：a和b为系数。
$$\mathrm{\Delta }{H}_k=a\left(p_k{p}_k^T\right)+b(H_k{q}_k)(H_k{q}_k{)}^T=a\left(P_k{P}_k^T\right)+b{H}_k{q}_k{q}_k^T{H}_k^T$$

代入拟牛顿法条件
$$p_k\approx H_{k+1}{q}_k$$$$\begin{array}{l}{p}_k\approx H_{k+1}{q}_k=\left(H_k+\Delta H_k\right)q_k\\ p_k=\left(H_k+a\left(p_k{p}_k^T\right)+b{H}_k{q}_k{q}_k^T{H}_k^T\right)q_k=H_k{q}_k+a{p}_k{p}_k^T{q}_k+b{H}_k{q}_k{q}_k^T{H}_k^T{q}_k\end{array}$$

目标转为求系数a和b。根据乘法的结合率，我们可以知道：$p_k^T{q}_k$和$q_k^T{H}_k^T{q}_k$为常数，因此，可得：
$$\begin{array}{l}\left(H_k+b{q}_k^T{H}_k^T{q}_k{H}_k\right)q_k+\left(a{p}_k^T{q}_k-1\right)p_k=0\\ \Rightarrow \left(1+b{q}_k^T{H}_k^T{q}_k\right)H_k{q}_k+\left(a{p}_k^T{q}_k-1\right)p_k=0\end{array}$$
设$p_k$和$H_k{q}_k$线性无关，则要使上式$=0$成立，则需：
$$\{\begin{array}{c}1+b{q}_k^T{H}_k^T{q}_k=0\\ a{p}_k^T{q}_k-1=0\end{array}\Rightarrow \{\begin{array}{c}a=\frac{1}{p_k^T{q}_k}\\ b=-\frac{1}{q_k^T{H}_k{q}_k}\end{array}$$
从而可得：
$$\mathrm{\Delta }{H}_k=\frac{p_k{p}_k^T}{p_k^T{q}_k}-\frac{H_k{q}_k\cdot q_k^T{H}_k}{q_k^T{H}_k{q}_k}$$
$H$为对称矩阵$H^T=H$

DFP流程：
设定初始点：$x_0$，迭代步数$k=0$， H为单位矩阵
开始迭代：
1）计算梯度，信息矩阵，从而得到牛顿迭代方向：
$$\nabla f(x_k),{\nabla }^{2}f(x_k),d_k=-H_k\cdot \nabla f(x_k)$$
2）判断当前迭代步长是否足够小或达到迭代次数，是则停止，否，继续下一步

3）更新：
$$\begin{array}{l}{x}_{k+1}=x_k+d_k;g_{k+1}=\mathrm{\Delta }f(x_{k+1}),p_k=x_{k+1}-x_k,q_k=g_{k+1}-g_k;\\ H_{k+1}=H_k+\frac{p_k{p}_k^T}{p_k^T{q}_k}-\frac{H_k{q}_k\cdot q_k^T{H}_k}{q_k^T{H}_k{q}_k}\\ k++;\end{array}$$
转步骤1）
示例代码：

#include <iostream>
#include <cmath>

using namespace std;

// 目标函数
double objectiveFunction(double x)
{
    return pow(x, 2) - 4 * x + 4;
}

// 目标函数的一阶导数
double derivativeFunction(double x)
{
    return 2 * x - 4;
}

// 拟牛顿法（BFGS算法）求解函数
double bfgsMethod(double initialGuess, double tolerance)
{
    double x = initialGuess;
    double dx = 0.0;
    double H = 1.0; // 初始Hessian逆矩阵
    float s = 0.1;
    do {
        double g = derivativeFunction(x);
        double step = -H * g;
        x = x + s*step;

        double dg = derivativeFunction(x) - g;
        
        // 更新Hessian逆矩阵
        double deltaX = step;
        double deltaY = dg;
        H += (deltaX * deltaX) / (deltaX * deltaY) - (H * deltaY * deltaY * H) / (deltaY * H * deltaY);

        dx = abs(step);
    } while (dx > tolerance);

    return x;
}

int main()
{
    double initialGuess = 1.0; // 初始猜测值
    double tolerance = 0.0001; // 精度

    double result = bfgsMethod(initialGuess, tolerance);

    cout << "Solution: " << result << endl;

    return 0;
}

2.3 BFGS

基本思想：由BFGS拟合Hessian阵来完成迭代。
$$x_k-x_{k+1}\approx {\nabla }^{2}f(x_{k+1}{)}^{-1}(g_k-g_{k+1}).\mapsto g_k-g_k\approx {\nabla }^{2}f(x_{k+1})(x_k-x_{k+1})$$
即$p_k\approx H_{k+1}{q}_k\mapsto q_k\approx B_{k+1}{p}_k$
求$\Delta B_k$的思路，即可以把$\mathrm{\Delta }{H}_k=\frac{p_k{p}_k^T}{p_k^T{q}_k}-\frac{H_k{q}_k\cdot q_k^T{H}_k}{q_k^T{H}_k{q}_k}$中的$p_k$换成$q_k$，将$H_k{q}_k$替换为$B_k{p}_k$，最后得BFGS校正公式：
$$\mathrm{\Delta }{B}_k=\frac{q_k{q}_k^T}{q_k^T{p}_k}-\frac{B_k{p}_k\cdot p_k^T{B}_k}{p_k^T{B}_k{p}_k}$$
该公式是较为有效率的一个校正公式。

此外，我们对$\Delta B_k$(DFP)的校正公式进行求逆，可以得到另一个的校正公式：
$$\mathrm{\Delta }{B}_k=(I-\frac{q_k{p}_k^T}{q_k^T{p}_k})B_k(I-\frac{p_k{q}_k^T}{q_k^T{p}_k})+\frac{q_k{q}_k^T}{q_k^T{p}_k}$$
示例代码

#include <iostream>
#include <cmath>

using namespace std;

// 目标函数
double objectiveFunction(double x)
{
    return pow(x, 2) - 4 * x + 4;
}

// 目标函数的一阶导数
double derivativeFunction(double x)
{
    return 2 * x - 4;
}

// 拟牛顿法（BFGS算法）求解函数
double bfgsMethod(double initialGuess, double tolerance)
{
    double x = initialGuess;
    double dx = 0.0;
     double B = 1.0; // 初始Hessian逆矩阵
     float s = 0.1;
     do {
         double g = derivativeFunction(x);
         double step = -g / B;
         x = x + s*step;

         double dg = derivativeFunction(x) - g;

         // 更新Hessian逆矩阵
         double deltaX = step;
         double deltaY = dg;
#if 1
        B += (1-(deltaY * deltaX) / (deltaY * deltaX)) * B *
             (1-(deltaX * deltaY) / (deltaY * deltaX))
           + (deltaY * deltaY) / (deltaY * deltaX);
#else
        B += (deltaY * deltaY) / (deltaY * deltaX)
           - (B * deltaX * deltaX * B) / (deltaX * B * deltaX);
#endif

        dx = abs(step);
    } while (dx > tolerance);

    return x;
}

int main()
{
    double initialGuess = 1.0; // 初始猜测值
    double tolerance = 0.0001; // 精度

    double result = bfgsMethod(initialGuess, tolerance);

    cout << "Solution: " << result << endl;

    return 0;
}

2维加步长确定：

#include <iostream>
#include <cmath>
#include <Eigen/Dense>
#include <Eigen/Core>

using namespace Eigen;
using namespace std;

// Objective function to minimize (example: quadratic function)
double objectiveFunction(const VectorXd& x) {
    return x[0]*x[0] + 2*x[1]*x[1];
}

// Gradient of the objective function
VectorXd gradient(const VectorXd& x) {
    VectorXd grad(2);
    grad[0] = 2*x[0];
    grad[1] = 4*x[1];
    return grad;
}

// BFGS optimization algorithm
void bfgsOptimization(VectorXd& x0, double epsilon, int maxIterations) {
    int n = x0.size();
    MatrixXd B = MatrixXd::Identity(n, n);  // Initial approximation of the inverse Hessian matrix
    VectorXd x = x0;

    for (int iter = 0; iter < maxIterations; ++iter) {
        VectorXd grad = gradient(x);
        double alpha = 1;  // Line search step size

        // Line search to determine the step size alpha
        //步长太大，缩一半
#if 1
        //if (f(x_k+1) > f(x_k) - 0.5*f'(x_k)*d_k)
        if (objectiveFunction(x - alpha * B * grad) > objectiveFunction(x) - 0.1 * alpha * grad.transpose() * B * grad) {
            alpha *= 0.1;
        }
        VectorXd s = -alpha * B.inverse() * grad;

#else

        if (objectiveFunction(x - alpha * grad) > objectiveFunction(x) - 0.1 * alpha * grad.transpose() * grad) {
            alpha *= 0.1;
        }
        VectorXd s = -alpha * grad;
#endif

        VectorXd xNext = x + s;
        VectorXd gradNext = gradient(xNext);
        VectorXd y = gradNext - grad;

        // BFGS update formula
        B = B + (y * y.transpose()) / (y.transpose() * s) - (B * s * s.transpose() * B) / (s.transpose() * B * s);

        x = xNext;

        cout << "iter: " << iter+1 << " x: " << x.transpose() << " cost: " << objectiveFunction(x) << endl;
        // Check for convergence
        if (grad.norm() < epsilon) {
            std::cout << "Converged after " << iter + 1 << " iterations." << std::endl;
            break;
        }
    }

    std::cout << "Optimal solution: " << x.transpose() << std::endl;
    std::cout << "Optimal value: " << objectiveFunction(x) << std::endl;
}

int main() {
    VectorXd x0(2);
    x0 << 1.0, 1.0;  // Initial guess

    double epsilon = 1e-6;  // Tolerance for convergence
    int maxIterations = 100;  // Maximum number of iterations

    bfgsOptimization(x0, epsilon, maxIterations);

    return 0;

拟牛顿法与线搜索

由于拟牛顿法使用的是近似的hessian矩阵，因此步长的选择非常重要，通常与线搜索方法结合。

#include <iostream>
#include <cmath>
#include <Eigen/Dense>

using namespace Eigen;
using namespace std;

// 定义目标函数 f(x)
double f(Vector2d x) {
    return x[0] * x[0] + 2 * x[0] * x[1] + x[1] * x[1];
}

// 计算目标函数在 x 处的梯度
Vector2d grad_f(Vector2d x) {
    Vector2d res;
    res[0] = 2 * x[0] + 2 * x[1];
    res[1] = 2 * x[0] + 2 * x[1];
    return res;
}

// BFGS 算法
Vector2d BFGS(Vector2d x0, double tol, int maxiter) {
    Vector2d x = x0;
    Matrix2d B = Matrix2d::Identity();
    int k = 0;
    while (k < maxiter) {
        Vector2d g = grad_f(x);
        if (g.norm() < tol) {
            break;
        }
        Vector2d p = -B * g;
        double alpha = 1.0;
        double rho = 0.9;
        double c = 0.0001;
        double phi = f(x);
        Vector2d grad_phi = grad_f(x);
        while (f(x + alpha * p) > phi + c * alpha * p.dot(grad_phi)) {
            alpha *= rho;
        }
        Vector2d s = alpha * p;
        Vector2d y = grad_f(x + s) - g;
        double rho_inv = 1 / y.dot(s);
        Matrix2d term1 = Matrix2d::Identity() - rho_inv * s * y.transpose();
        Matrix2d term2 = rho_inv * s * s.transpose();
        B = term1 * B * term1.transpose() + term2;
        x += s;
        k++;
    }
    return x;
}

int main() {
    Vector2d x0(1, 2); // 初始参数值
    double tol = 1e-6; // 迭代停止容差
    int maxiter = 100; // 最大迭代次数
    Vector2d res = BFGS(x0, tol, maxiter);
    cout << "Optimal parameters: (" << res[0] << ", " << res[1] << ")" << endl;
    cout << "Optimal function value: " << f(res) << endl;
    return 0;
}

参考链接
https://www.bilibili.com/video/BV1e64y1Y7Sr/