Lab 4 - 非线性最小二乘

本次实验我们实现 Gauss Newton 法求解简单的非线性最小二乘问题。
本次实验内容为课程作业，计算成绩。你需要将源代码、可执行程序（注明运行环境）和实验文档上传至学在浙大，标题为Lab4-学号-姓名，压缩包和文件夹名称都为Lab4-学号-姓名。
本次作业提交截止时间：2021年04月14日 23:59:59，逾期将要扣分。

知识回顾

Gauss Newton 法来源于 Newton Rapson 法，后者利用二次型逼近目标函数。为了进行二次逼近，NR 法需要计算目标函数的 Hessian，但这并不是个容易完成的任务。GN 方法利用最小二乘问题本身的特点，通过计算 Jacobian 来近似 Hessian 。

NR 法采用二阶 Taylor 展开近似目标函数：

$F(x+\Delta x) \approx F(x)+J_F\Delta x+\frac12\Delta x^TH_F\Delta x$

最小二乘问题具有 $F(x) = \|R(x)\|_2^2$ 的形式，对 $R$ 进行一阶 Taylor 展开，可以得到：

$\begin{aligned} F(x+\Delta x)&=\|R(x+\Delta x)\|_2^2\\ &\approx \|R(x)+J_R\Delta x\|_2^2\\ &=\|R(x)\|_2^2+2R^TJ_R\Delta x+\Delta x^TJ_R^TJ_R\Delta x\\ &= F(x)+J_F\Delta x+\Delta x^TJ_R^TJ_R\Delta x \end{aligned}$

（请务必自己推导一次，理解全部过程, Idea: Line search）

将两种展开方式进行对比，可知在最小二乘问题里 $H_F\approx 2J_R^TJ_R$ 。

继续依照 NR 法的思路，每一步迭代中，我们求解如下的标准方程：

$J_R^TJ_R\Delta x+J_R^TR=0$

这一标准方程对应于求 $J_R\Delta x=-R$ 的（线性）最小二乘解，可以采用共轭梯度法求解。本次实验中，可以直接采用 OpenCV 提供的矩阵算法完成。

得到 $\Delta x$ 后，以它作为下降方向，我们进行线性搜索求出合适的步长$\alpha$，更新 $x$ ：

$x\leftarrow x+\alpha\Delta x$

这样就得到了 GN 法，GN 法的算法框架如下：

• $x\leftarrow x_0$
• $n\leftarrow 0$
• $\text{while }n < n_{\max} :$
  • $\Delta x \leftarrow\text{Solution of }J_R\Delta x=-R :$
    • $\text{Conjugate Gradient or Other}$
  • $\text{if }\|R\|_\infty\leq\varepsilon_r \vee \|\Delta x\|_\infty\leq\varepsilon_g\text{ return }x$
  • $\alpha\leftarrow\arg\min_\alpha \{x+\alpha\Delta x\}$
  • $x\leftarrow x+\alpha\Delta x$
  • $n\leftarrow n+1$

实验要求

本次实验你需要你实现 Gauss Newton 求解最小二乘问题。

接口要求

你需要实现名为 Solver#### 的类，其中 #### 代表你学号的后四位。
这个类需要继承自接口类 GaussNewtonSolver ，我们提供了头文件，请下载 hw3_gn.h 并包含在你的优化器类中。
你需要严格按照接口要求实现，否则会扣分。

GaussNewtonSolver 的定义如下：

class GaussNewtonSolver {
public:
    virtual double solve(
               ResidualFunction *f, // 目标函数
               double *X,           // 输入作为初值，输出作为结果
               GaussNewtonParams param = GaussNewtonParams(), // 优化参数
               GaussNewtonReport *report = nullptr // 优化结果报告
               ) = 0;
};

• 当你的类的 solve 函数被调用时，它将采用 Gauss Newton 法最小化函数 f 。
• 输入时数组 X 内包含初始点，并且在优化后将被修改为最优点，维度需要与目标函数定义的维度一致。
• param 是 GN 算法运行的各种相关参数。
• 如果 report 不是空指针，优化完成后应当把相关的报告记录到 report 。
• 函数返回值是优化得到的目标函数的最优值。

目标函数通过继承接口类 ResidualFunction 进行定义：

class ResidualFunction {
public:
    virtual int nR() const = 0;
    virtual int nX() const = 0;
    virtual void eval(double *R, double *J, double *X) = 0;
};

nR 函数需要返回余项向量的维度。类似的，nX 函数需要返回变量 X 的维度。
eval 运行时读入 X 将计算得到的余项和Jacobian写入 R 和 J 。

• 计算得到的余项向量需要写入到 R ，R 需要预先分配好，长度为 nR。
• 计算得到的 Jacobian 需要写入到 J ，J 需要预先分配好，大小为 nR*nX。
• 输入的 X 是一个长度为 nX 的数组，包含所有的变量。

注意：你的优化器需要负责在优化开始前分配好 R 和 J 需要的空间，在结束后销毁。X 作为输入，由用户（调用 solve 的程序）分配并填充好初始值。

优化的参数通过如下的结构体定义，各个参数的含义见注释。

struct GaussNewtonParams{
    GaussNewtonParams() :
       exact_line_search(false),
       gradient_tolerance(1e-5),
       residual_tolerance(1e-5),
       max_iter(1000),
       verbose(false)
    {}
    bool exact_line_search; // 使用精确线性搜索还是近似线性搜索
    double gradient_tolerance; // 梯度阈值，当前梯度小于这个阈值时停止迭代
    double residual_tolerance; // 余项阈值，当前余项小于这个阈值时停止迭代
    int max_iter; // 最大迭代步数
    bool verbose; // 是否打印每步迭代的信息
};

优化结果的详细信息储存的结构体定义如下，含义见注释。

struct GaussNewtonReport {
    enum StopType {
        STOP_GRAD_TOL,       // 梯度达到阈值
        STOP_RESIDUAL_TOL,   // 余项达到阈值
        STOP_NO_CONVERGE,    // 不收敛
        STOP_NUMERIC_FAILURE // 其它数值错误
    };
    StopType stop_type; // 优化终止的原因
    double n_iter;      // 迭代次数
};

测试问题

作为测试，我们求解从三维点云拟合椭球的问题，我们的椭球模型是

$\frac{x^2}{A^2}+\frac{y^2}{B^2}+\frac{z^2}{C^2} = 1$

下载测试数据文件 ellipse753.txt ，其中包含了 753 个含有噪音的点，你需要建立目标函数，然后利用你的优化器最小化目标函数，从这个点云恢复出参数 AA BB CC 。

注意：这个测试问题同样可以用线性最小二乘求解，作业中必须使用非线性最小二乘的形式，但同时你可以用线性最小二乘求解来进行验证。

参考资料：K.Madsen et.al., Methods for Non-Linear Least Square Problems. 2nd Ed., 2014. PDF