非线性优化

习题一:矩阵微分

在优化中经常会遇到矩阵微分的问题。例如,当自变量为向量$x$,求标量函数$u(x)$对$x$的导数时,即为矩阵微分。通常线性代数教材不会深入探讨此事,这往往是矩阵论的内容。回答下列问题: 设变量为$x\in{R^N}$,那么:

  1. 矩阵$A\in{R^{N×N}}$,那么$d(Ax)/dx$是什么?
  2. 矩阵$A\in{R^{N×N}}$,那么$d(x^TAx)/dx$是什么?
  3. 证明:

    $$ x^TA^Tx = tr(Axx^T ). $$

    分析

  4. 涉及到矩阵和向量的求导,总共有五种类型,注意相关的顺序以及于雅各比矩阵的区别,再强调一次,雅各比矩阵的列数为自变量的个数,这一点要记清楚。 A.向量对标量; B.标量对向量; C.向量对向量; D.矩阵对标量; E.标量对矩阵; 他们分别对应的形式如下所示:
    $$ 向量对标量: \frac{\partial{Y}}{\partial{x}} = \begin{bmatrix} \frac{\partial{y_1}}{\partial{x}} \ \frac{\partial{y_2}}{\partial{x}} \ \vdots \ \frac{\partial{y_m}}{\partial{x}} \end{bmatrix} $$

$$ 标量对向量: \frac{\partial{y}}{\partial{X}} = \begin{bmatrix} \frac{\partial{y}}{\partial{x_1}} & \frac{\partial{y}}{\partial{x_2}} & \cdots & \frac{\partial{y}}{\partial{x_n}} \end{bmatrix} $$

$$ 向量对向量: \frac{\partial{Y}}{\partial{X}} = \begin{bmatrix} \frac{\partial{y_1}}{\partial{x_1}} & \frac{\partial{y_2}}{\partial{x_1}} & \cdots & \frac{\partial{y_m}}{\partial{x_1}} \\ \frac{\partial{y_1}}{\partial{x_2}} & \frac{\partial{y_2}}{\partial{x_2}} & \cdots & \frac{\partial{y_m}}{\partial{x_2}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial{y_1}}{\partial{x_n}} & \frac{\partial{y_2}}{\partial{x_n}} & \cdots & \frac{\partial{y_m}}{\partial{x_n}} \\ \end{bmatrix} $$

$$ 矩阵对标量: \frac{\partial{Y}}{\partial{x}} = \begin{bmatrix} \frac{\partial{y_{11}}}{\partial{x}} & \frac{\partial{y_{21}}}{\partial{x}} & \cdots & \frac{\partial{y_{m1}}}{\partial{x}} \\\\ \frac{\partial{y_{12}}}{\partial{x}} & \frac{\partial{y_{22}}}{\partial{x}} & \cdots & \frac{\partial{y_{m2}}}{\partial{x}} \\\\ \vdots & \vdots & \ddots & \vdots \\\\ \frac{\partial{y_{1n}}}{\partial{x}} & \frac{\partial{y_{2n}}}{\partial{x}} & \cdots & \frac{\partial{y_{mn}}}{\partial{x}} \\\\ \end{bmatrix} $$

$$ 标量对矩阵: \frac{\partial{y}}{\partial{X}} = \begin{bmatrix} \frac{\partial{y}}{\partial{x_{11}}} & \frac{\partial{y}}{\partial{x_{12}}} & \cdots & \frac{\partial{y}}{\partial{x_{1n}}} \\ \frac{\partial{y}}{\partial{x_{21}}} & \frac{\partial{y}}{\partial{x{22}}} & \cdots & \frac{\partial{y}}{\partial{x_{2n}}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial{y}}{\partial{x_{m1}}} & \frac{\partial{y}}{\partial{x_{m2}}} & \cdots & \frac{\partial{y}}{\partial{x_{mn}}} \\ \end{bmatrix} $$

我们常用的矩阵求导公式有:

$$ Y = A * X –> \frac{\partial{Y}}{\partial{X}} = A^T $$ $$ Y = X * A –> \frac{\partial{Y}}{\partial{X}} = A
$$ $$ Y = A^T * X * B –> \frac{\partial{Y}}{\partial{X}} = A * B^T
$$ $$ Y = A^T * X^T * B –> \frac{\partial{Y}}{\partial{X}} = B * A^T
$$ $$ \frac{\partial{X^TX}}{\partial{X}} = X $$

答案:

1、设$Y=AX$,则$Y\in{R^{N*1}}$,所以$\frac{\partial{Y}}{\partial{X}}$属于向量对向量的形式,其中Y向量我们可以计算得出:

$$ Y = \begin{bmatrix} a_{11}*x_1+a_{12}*x_2+ \cdots + a_{1n}*x_n\\ a_{21}*x_1+a_{22}*x_2+ \cdots + a_{2n}*x_n \\ \vdots \\ a_{n1}*x_1+a_{n2}*x_2+ \cdots + a_{nn}*x_n \end{bmatrix} $$

则我们由上述向量对向量的式子可以得到:

$$ \frac{\partial{Y}}{\partial{X}} = \begin{bmatrix} \frac{\partial{y_1}}{\partial{x_1}} & \frac{\partial{y_2}}{\partial{x_1}} & \cdots & \frac{\partial{y_m}}{\partial{x_1}} \\ \frac{\partial{y_1}}{\partial{x_2}} & \frac{\partial{y_2}}{\partial{x_2}} & \cdots & \frac{\partial{y_m}}{\partial{x_2}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial{y_1}}{\partial{x_n}} & \frac{\partial{y_2}}{\partial{x_n}} & \cdots & \frac{\partial{y_m}}{\partial{x_n}} \\ \end{bmatrix} $$

带入YX变量可得:

$$ \frac{\partial{Y}}{\partial{X}} = \begin{pmatrix} a_ {11} & a_ {21} & \cdots & a_ {n1} \\ a_ {12} & a_ {22} & \cdots & a_ {n2}\\ \vdots & \vdots & \ddots & \vdots \\ a_ {1n} & a_ {2n} & \cdots & a_ {nn} \\ \end{pmatrix} = A^T $$

2、由第一问我们可以知道,$Y=AX$,我们设 $ Z = x^TY$,则我们可以知道$Z\in{R^{1*1}}$,则

$$ Z = \sum_ {i=1}^n\sum_ {j=1}^n a_ {ij}*x_i*x_j $$

则我们由上述标量对向量的式子可以得到:

$$ \frac{\partial{z}}{\partial{X}} = \begin{bmatrix} \frac{\partial{z}}{\partial{x_1}} & \frac{\partial{z}}{\partial{x_2}} & \cdots & \frac{\partial{z}}{\partial{x_n}} \end{bmatrix} $$

我们计算出Z对与 $x_l$的导数为:

$$ \frac{\partial{z}}{\partial{x_l}} = \sum_{i=1}^nx_i*a_{il}+\sum_{j=1}^nx_j*a_{lj} = (A^T+A)x $$

3、根据上述的第二问,我们可以将$x^TA^Tx$展开,同时可以将$Axx^T$展开,最终我们能够证明$x^TA^Tx = tr(Axx^T)$,这里不再展开叙述。

习题二:高斯牛顿法的曲线拟合实验

我们在课上演示了用Ceresg2o进行曲线拟合的实验,可以看到优化框架给我们带来了诸多便利。本题中你需要自己实现一遍高斯牛顿的迭代过程,求解曲线的参数。我们将原题复述如下。设有曲线满足以下方程:

$$ y = exp(ax^2+bx+c)+w $$

其中a,b,c为曲线参数,w为噪声。现有N个数据点(x,y),希望通过此N个点来拟合a,b,c。实验中取$N = 100$。
那么,定义误差为 $ e_i = y_i-exp(ax_i^2+bx_i+c)$,于是$(a,b,c)$的最优解可通过解以下最小二乘获得:

$$ min_{a b c} = \frac{1}{2}\sum_{i=1}^N {||y_i - exp(ax_i^2+bx_i+c)||^2} $$

现在请你书写 Gauss-Newton 的程序以解决此问题。

分析

  1. 注意雅克比矩阵和矩阵微分的形式不同,雅克比矩阵列数为自变量的个数,例如本题中,雅克比矩阵的大小应为 $1*3$,雅克比矩阵和梯度矩阵的关系$A=J^T$。
  2. 程序计算的时候,100个H矩阵相加,这一点可以考虑为100数据都这样的大矩阵,分模块相乘就是各个小模块相乘再累加的过程。
  3. 一阶梯度法就是快速下降法,但是快速下降法过于贪心,容易出现锯齿路线,反而增加了迭代的次数。
  4. 二阶梯度法就是牛顿法,缺点是需要计算H矩阵,在大规模的问题中这是很难实现的,我们一般避免计算H矩阵。
  5. 一阶和二阶的方法,都是将$f(x)^2$在$x$附近泰勒展开,但是高斯牛顿的方法是将$f(x)$进行一阶泰勒展开。
  6. 对于列文伯格-马夸尔特的方法,一般认为比高斯牛顿更为健壮,但是它收敛的速度可能会慢于高斯牛顿法,但却在SLAM中大量应用。

答案

程序最终结果如图所示:

完整程序链接:https://github.com/XLMaverick/Visual-Localization-Percessing/tree/master/vslam_fourteen_lectures/gauessnewton
代码如下:

    #include <iostream>
    #include <opencv2/opencv.hpp>
    #include <Eigen/Core>
    #include <Eigen/Dense>

    using namespace std;
    using namespace Eigen;

    int main(int argc, char **argv) 
    {
        double ar = 1.0, br = 2.0, cr = 1.0;         // 真实参数值
        double ae = 2.0, be = -1.0, ce = 5.0;        // 估计参数值
        int N = 100;                                 // 数据点
        double w_sigma = 1.0;                        // 噪声Sigma值
        cv::RNG rng;                                 // OpenCV随机数产生器

        vector<double> x_data, y_data;      // 数据
        for (int i = 0; i < N; i++) 
        {
            double x = i / 100.0;
            x_data.push_back(x);
            y_data.push_back(exp(ar * x * x + br * x + cr) + rng.gaussian(w_sigma));
        }

        // 开始Gauss-Newton迭代
        int iterations = 100;    // 迭代次数
        double cost = 0, lastCost = 0;  // 本次迭代的cost和上一次迭代的cost

        for (int iter = 0; iter < iterations; iter++) 
        {
            Matrix3d H = Matrix3d::Zero();             // Hessian = J^T J in Gauss-Newton
            Vector3d b = Vector3d::Zero();             // bias
            cost = 0;

            for (int i = 0; i < N; i++) 
            {
                double xi = x_data[i], yi = y_data[i];  // 第i个数据点
                double error = 0;   // 第i个数据点的计算误差
                error =yi -exp(ae*xi*xi+be*xi+ce); // 填写计算error的表达式
                Vector3d J; // 雅可比矩阵
                J[0] = -exp(ae*xi*xi+be*xi+ce)*xi*xi;  // de/da
                J[1] = -exp(ae*xi*xi+be*xi+ce)*xi;  // de/db
                J[2] = -exp(ae*xi*xi+be*xi+ce);  // de/dc

                H += J * J.transpose(); // GN近似的H
                b += -error * J;

                cost += error * error;
            }

            // 求解线性方程 Hx=b,建议用ldlt
            Vector3d dx;
            dx = H.ldlt().solve(b);

            if (isnan(dx[0])) 
            {
                cout << "result is nan!" << endl;
                break;
            }

            if (iter > 0 && cost > lastCost) {
                // 误差增长了,说明近似的不够好
                cout << "cost: " << cost << ", last cost: " << lastCost << endl;
                break;
            }

            // 更新abc估计值
            ae += dx[0];
            be += dx[1];
            ce += dx[2];

            lastCost = cost;

            cout << "total cost: " << cost << endl;
        }

        cout << "estimated abc = " << ae << ", " << be << ", " << ce << endl;
        return 0;
    }