求根

求解 $\frac{1}{2}(10-x^2)$ 的根：

#include <iostream>
#include <ceres/ceres.h>

/* 1.构建代价函数 */
class CostFunc{
public:
    CostFunc()=default;
    template<typename T>
    bool operator()(const T* x, T* residual) const{
        *residual = 0.5 * pow(T(10.0) - *x, 2);
        return true;
    }
};

int main(int argc, char** argv){
    /* 2.构建待求解的优化问题 */
    ceres::Problem problem;
    double x = 9.0; // 待优化参数，初始为9
    // 使用自动求导, <代价函数, 残差维度, 待优化参数的维度>
    ceres::CostFunction* cost_function 
        = new ceres::AutoDiffCostFunction<CostFunc, 1, 1>(new CostFunc);
    problem.AddResidualBlock(   cost_function,  // 代价函数
                             	nullptr,        // 核函数
                             	&x              // 待优化参数
                            );
	/* 3.配置求解器参数 */
    ceres::Solver::Options options;
    options.linear_solver_type = ceres::DENSE_QR; // 增量方程的解法
    options.minimizer_progress_to_stdout = true;
    ceres::Solver::Summary summary;
    ceres::Solve(options, &problem, &summary);
    /* 求解 */
    std::cout << summary.BriefReport() << std::endl;
    std::cout << "result: " << x << std::endl;
}

结果为 $x = 9.99951$

曲线拟合

以 $y = 3 ? s i n (x) + 2 ? c o s (x) + 1$ 为真实值生成观测数据，如下图所示，原始数据见附录。
在这里插入图片描述

最小二乘

使用最小二乘法求解，令观测模型为：
$y = a ? s i n (x) + b ? c o s (x) + c$
令观测向量 $h_i=\begin{bmatrix}sin(x_i)&cos(x_i)&1\end{bmatrix}^T$ ，待优化参数 $p=\begin{bmatrix}a&b&c\end{bmatrix}^T$ ，则 $y_i=h_i^Tp$ 。将观测向量 $h_i$ 按行排列得到观测矩阵：
$H=\begin{bmatrix}h_1^T\\h_2^T\\\vdots\\h_n^T\end{bmatrix}$
观测方程如下：
$Y=\begin{bmatrix}y_1\\y_2\\\vdots\\y_n\end{bmatrix}=Hp$
两端同时乘以观测矩阵的转置 $H^T$ ，得到：
$H^T=H^THp$
由于矩阵乘以矩阵的转置其结果一定可逆，所以解得：
$p=(H^TH)^{-1}H^TY$
这就是待求解参数的最小二乘解。Matlab求解程序如下：

x = 0 : 0.2 : 10;
x = x';
y_true = 3 * sin(x) + 2 * cos(x) + 1; % 真实值

% 生成噪声
w = wgn(length(x), 1, 10*log(0.4));
% 产生观测值
z = y_true + w;
% 生成观测矩阵
H = [sin(x), cos(x), x.^0];
% 计算最小二乘估计
p = (H' * H)^(-1) * H' * z;
% 计算均方误差
y = H * p;
mse = (y - y_true)' * (y - y_true) / length(y);
fprintf("mse: %f\n", mse);
% 可视化
plot(x, y_true)
hold on 
plot(x, z)
hold on 
plot(x, y)
legend("真实值", "观测值", "拟合值")

结果如下，与真实值基本一致。

Ceres求解

使用ceres求解，程序如下：

#include <iostream>
#include <ceres/ceres.h>

#include <fstream>
#include <istream>

/* 1.构建代价函数 */
class CostFunc{
public:
    CostFunc()=default;
    CostFunc(double x, double y) : x_(x), y_(y){ }
    template<typename T>
    bool operator()(const T* abc, T* residual) const{
        // residual = y - f(x)
        // f(x) = a·sin(x) + b·cos(x) + c
        *residual = T(y_) - (abc[0] * sin(T(x_)) + abc[1] * cos(T( x_)) + abc[2]);
        return true;
    }
private:
    double x_, y_;
};

int main(int argc, char** argv){
    ceres::Problem problem;
    double abc[] = {0, 0, 0}; // 待优化参数，全部初始化为0
    /* 2.构建待求解的优化问题 */
    std::fstream data_fs("lse_data.txt"); // 从文件读取观测数据
    double x, y; // 观测数据
    // 添加代价项
    for(int i = 0; i < 51; i++){
        data_fs >> x >> y;
        // 使用自动求导, <代价函数, 残差维度, 待优化参数的维度>
        ceres::CostFunction* cost_function 
            = new ceres::AutoDiffCostFunction<CostFunc, 1, 3>(new CostFunc(x, y));
        problem.AddResidualBlock(   cost_function,  // 代价函数
                                    nullptr,        // 核函数
                                    abc             // 待优化参数
                                );
    }
    /* 3.配置求解器参数 */
    ceres::Solver::Options options; 
    options.linear_solver_type = ceres::DENSE_QR; // 增量方程的解法
    options.minimizer_progress_to_stdout = true;
    ceres::Solver::Summary summary; // 求解结果摘要
    /* 求解 */
    ceres::Solve(options, &problem, &summary);
    std::cout << summary.BriefReport() << std::endl;
    std::cout << "result: " << abc[0] << " " << abc[1] << " " << abc[2] << std::endl;
}

结果如下：

iter      cost      cost_change  |gradient|   |step|    tr_ratio  tr_radius  ls_iter  iter_time  total_time
   0  2.147215e+02    0.00e+00    8.29e+01   0.00e+00   0.00e+00  1.00e+04        0    2.19e-05    4.89e-05
   1  1.659848e-06    2.15e+02    7.21e-03   3.74e+00   1.00e+00  3.00e+04        1    2.10e-05    9.11e-05
   2  1.511573e-15    1.66e-06    2.15e-07   3.39e-04   1.00e+00  9.00e+04        1    1.60e-05    1.11e-04
Ceres Solver Report: Iterations: 3, Initial cost: 2.147215e+02, Final cost: 1.511573e-15, Termination: CONVERGENCE
result: 3 2 1

附录

曲线拟合观测数据

0                  3
0.200000000000000  3.55614114806767
0.400000000000000  4.01037701493172
0.600000000000000  4.34459865000446
0.800000000000000  4.54548169139290
1                  4.60501756615997
1.20000000000000   4.52083276685503
1.40000000000000   4.29628347576586
1.60000000000000   3.94032176452194
1.80000000000000   3.46713870324841
2                  2.89559860738276
2.20000000000000   2.24848697694808
2.40000000000000   1.55160211057096
2.60000000000000   0.832726608726498
2.80000000000000   0.120519769130397
3                  -0.556624969021289
3.20000000000000   -1.17171198187225
3.40000000000000   -1.70021969123942
3.60000000000000   -2.12107816255285
3.80000000000000   -2.41750909665699
4                  -2.57769472765101
4.20000000000000   -2.59524895992216
4.40000000000000   -2.46947196162539
4.60000000000000   -2.20537806477050
4.80000000000000   -1.81349585962863
5                  -1.30944845306296
5.20000000000000   -0.713330624559708
5.40000000000000   -0.0489077107826954
5.60000000000000   0.657331843403535
5.80000000000000   1.37723249564137
6                  2.08209407870395
6.20000000000000   2.74381598559394
6.40000000000000   3.33601745206787
6.60000000000000   3.83508927445719
6.80000000000000   4.22113503411548
7                  4.47876430484298
7.20000000000000   4.59770622061197
7.40000000000000   4.57321894258366
7.60000000000000   4.40627870125897
7.80000000000000   4.10354087724912
8                  3.67707467225292
8.20000000000000   3.14388194807165
8.40000000000000   2.52521941603147
8.60000000000000   1.84575119898232
8.80000000000000   1.13256555055197
9                  0.414094931955916
9.20000000000000   -0.281017500507584
9.40000000000000   -0.925059807710340
9.60000000000000   -1.49235605525719
9.80000000000000   -1.96028993196529
10                 -2.31020639082101