# Ceres Library: basic use, taking handwritten Gauss Newton method as an example

Posted by webwiese on Thu, 17 Feb 2022 14:13:54 +0100

# Ceres Library

## brief introduction

Ceres library is an open source C + + nonlinear optimization library developed by Google, which is widely used to solve least squares problems.

## install

git clone https://ceres-solver.googlesource.com/ceres-solver


Second, you need to install dependencies:

sudo apt install liblapack-dev libsuitesparse-dev libcxsparse3 libgflags-dev libgoogle-glog-dev libgtest-dev


Then, enter the source directory for compilation and installation:

cd ceres-solver
mkdir build && cd build
cmake ..
make -j12
sudo make install


Default installation file locations: \ usr\local\include\ceres and \ usr\local\lib

## Engineering configuration

Import the library file, cmakelists. In the project root directory Txt fill in as follows:

# Ceres-solver
find_package(Ceres REQUIRED)
include_directories(${CERES_INCLUDE_DIRS})  Cmakelists under src Txt file is filled as follows: # link ceres-solver target_link_libraries(Study${CERES_LIBRARIES})


include header file main H fill in the following:

//  Ceres-solver
#include <ceres/ceres.h>


## Problem description

In the Cere library, the least squares problem is modeled as follows (kernel least squares with boundary):
min ⁡ x 1 2 ∑ i ρ i ( ∥ f i ( x i 1 , x i 2 , ⋯   , x i k ) ∥ 2 ) s . t . l j ≤ x j ≤ u j \min_x\frac{1}{2}\sum_i\rho_i\Bigl(\begin{Vmatrix}f_i(x_{i_1},x_{i_2},\cdots,x_{i_k})\end{Vmatrix}^2\Bigr)\\ s.t.\quad l_j\leq x_j \leq u_j xmin​21​i∑​ρi​(∥∥​fi​(xi1​​,xi2​​,⋯,xik​​)​∥∥​2)s.t.lj​≤xj​≤uj​
Among the questions expressed above, x 1 , ⋯   , x n x_1,\cdots,x_n x1,..., xn are called optimization variables, also known as Parameter blocks; f i ( ⋅ ) f_i(\cdot) fi (⋅) is called Cost function; ρ i ( ⋅ ) \rho_i(\cdot) ρ i (⋅) is called Loss function, also known as kernel function; ρ i ( ∥ f i ( x i 1 , x i 2 , ⋯   , x i k ) ∥ 2 ) \rho_i\Bigl(\begin{Vmatrix}f_i(x_{i_1},x_{i_2},\cdots,x_{i_k})\end{Vmatrix}^2\Bigr) ρ i (‖ ‖ fi (xi1, xi2,..., xik) ‖ ‖ 2) is called Residual block; parameter l j l_j lj # and u j u_j uj is called the lower and upper bounds of the parameter.

When the loss function is taken as a constant value: ρ i ( x ) = x \rho_i(x)=x ρ i (x)=x, at this time, the upper and lower limits of the parameters are established as l j = − ∞ , u j = ∞ l_j=-\infin,u_j=\infin lj = − ∞, uj = ∞, that is, it is constructed as a least squares problem with unbounded constraints. At this time, the cost function is equivalent to the loss function, that is, the residual block.

# Ceres library API

Using Ceres library to solve the least squares problem is mainly composed of three parts: the construction of cost function, the construction of least squares problem and the solution of least squares problem.

In order to facilitate the learning of Ceres library, the handwritten Gauss Newton method is still used here as an example. The problem is described as follows: https://blog.csdn.net/weixin_45929038/article/details/122973583

Similarly, extract points as observation data sets:

int main()
{
/*--------  Initial parameter configuration--------*/
//  Actual curve parameters
double ar = 1.0, br = 1.0, cr = 1.0;
//  Initial value of estimated curve parameters
double ae = 2.0, be = 1.0, ce = 5.0;
//  Number of sampling observation data points
int N = 100;
//  Noise standard deviation and its reciprocal
double w_sigma = 1.0;
//  Random number generator
cv::RNG rng;

/*--------  Observation data generation--------*/
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*w_sigma));
}
//  Parameter block
double abc = {ae, be, ce};

return 0;
}


## Constructing cost function

The cost function is constructed by using the imitation function of C + +. Its essence is the structure struct or class, and it overloads () to make it have the same call mode as the function, so it is called imitation function. Struct is used for definition here, which mainly includes constructor and overloaded operator ():

For example, the cost function for solving curve parameters is constructed as follows:

/*--------  Build CostFunction--------*/
struct CURVE_FITTING_COST{
//  Constructor
CURVE_FITTING_COST(double x, double y): _x(x), _y(y){}

template<typename T>
bool operator()(const T *const abc, T *residual) const
{
//  y - exp(ax^2 + bx + c)
residual = T(_y) - ceres::exp(abc * T(_x) * T(_x) + abc * T(_x) + abc);
return true;
}

//  Member variables that store values
const double _x, _y;
};


### Constructor

The constructor is used to pass the value of known parameters (sensor observation data) into the cost function. When there are no known parameters in the optimized problem, there is no need to write the constructor.

In this example, 100 points on the curve are used as input.

The operator () is a template method whose return value is of type bool. The accepted parameters are parameter block and residual block.

#### input

Parameters should be passed into the array of variables at one time using pointers. For example, for the above problems, the optimized parameters include a , b , c a,b,c a. b and c are three parameters, so the transmission parameter block is * const abc, which corresponds to an array of three parameters.

The above problem is an unconstrained least squares problem, so the residual block is equivalent to the cost function y i − e x p ( a x i 2 + b x i + c ) y_i-exp(ax_i^2+bx_i+c) yi​−exp(axi2​+bxi​+c).

#### Template type T

Since different open source algorithm libraries have their own defined data types (Vector in Eigen corresponds to Point in Opencv and Matrix corresponds to Matrix), the input and output of operators are unified as template type T.

#### const keyword

In the process of solving the least squares, we don't want to modify the contents of overloaded operator () and parameter block due to abnormal operation, so we use const to modify it.

## Constructing least squares problem

The least squares problem is mainly constructed by Ceres::Problem class:

//  instantiation
ceres::Problem problem;


After instantiating this class, the least squares problem is constructed by calling the corresponding member function.

### Residual block

Use the following function to pass the residuals:

problem.AddResidualBlock(CostFunction* cost_function,
LossFunction* loss_function,
const vector<double *> parameter_blocks);

• cost_function: cost function
• loss_function: loss function. If the constructed problem is an unconstrained least squares problem, the input parameter is nullptr
• parameter_blocks: parameter blocks, that is, the parameters of the final solution

The cost function is described as follows:

### Derivative model of cost function

Ceres provides three template classes for derivation with different accuracy and efficiency:

• AutoDiffCostFunction: automatic derivation class. The derivation method is determined by Ceres library. It is the longest used class.
• NumericDiffCostFunction: numerical derivative class. The derivative numerical solution form is written manually by the user. It is usually used when automatic derivation cannot be called directly. Its accuracy and efficiency are lower than that of automatic derivation
• Analytical derivatives: analytical derivatives. It is used when derivatives have closed analytical form. It needs to manage residuals and Jacobian matrix by itself.

Ceres officially recommends automatic derivation:

ceres::AutoDiffCostFunction<CostFunctor, int residualDim, int paramDim>(CostFunctor* functor);

• CostFunctor: the cost function built in the first step
• residualDim: residual block latitude, that is, the output latitude of the model
• paramDim: parameter block latitude, that is, the input latitude of the model
• Function: pointer, which is the construction of cost function

For example, the problems in handwritten Gauss Newton method can be constructed as follows:

ceres::AutoDiffCostFunction<CURVE_FITTING_COST, 1, 3>(new CURVE_FITTING_COST(x_data[i], y_data[i]))


Among them, the residual latitude is 1, that is, the error in the problem e i e_i ei​； The parameter latitude is 3, that is, the latitude of the curve parameter a , b , c a,b,c a,b,c； The passed in parameters are used to build the cost function of each time.

For this problem, the configuration of the construction residual term is as follows:

for(int i = 0; i < 100; i++){
//  Residual term configuration
problem.AddResidualBlock(new ceres::AutoDiffCostFunction<CURVE_FITTING_COST, 1, 3>(new CURVE_FITTING_COST(x_data[i], y_data[i])),
nullptr,
abc);
}


Here, the loss function (unconstrained) is not used, so it is passed into nullptr; The parameter block is a three-dimensional matrix abc

### other

See the reference blog: Ceres detailed explanation (I) Problem class

## Solving the least squares problem

Use ceres::Solve to solve, and its function prototype is as follows:

void Solve(const Solver::Options& options, Problem* problem, Solver::Summary* summary)

• Options: configuration of solver and configuration options of solution
• Problem: the problem to be solved, that is, the least squares problem we build
• summary: the optimization information of the solution, which is used to store the optimization information in the solution process

The configuration of the solver is described as follows:

### Solver::Options

• minimizer_type: iterative solution method. Options are as follows:
• TRUST_REGION: trust region method, default value
• LINEAR_SEARCH: linear search method
• Parameters usually remain at their default values
• trust_ region_ strategy_ The trust domain policy is as follows:
• LEVENBERG_MARQUARDT, lewenberg Marquardt method, default
• DOGLEG: dog leg method, commonly known as DOGLEG method
• linear_solver_type: method of solving linear equations
• DENSE_QR: QR decomposition method, default value, used for small-scale least squares solution
• DENSE_NORMAL_CHOLESKY and SPARSE_NORMAL_CHOLESKY: CHolesky decomposition, which is used for solving large-scale nonlinear least squares with sparsity
• CGNR: conjugate gradient method for solving sparse equations
• DENSE_ SCUR and spark_ SCHUR: SCHUR decomposition, used to solve BA problems
• ITERATIVE_SCHUR: conjugate gradient SCHUR, used to solve BA problem
• minimizer_progress_to_stdout: whether to output the optimization information to the terminal
• bool type, false by default. If set to true, the output is different according to the iteration method:
• Trust region method
• cost: the value of the current objective function
• cost_change: objective function change caused by current parameter change
• |step|: parameter variation
• tr_ratio: the ratio of the actual change of the objective function to the change of the objective function in the trust region
• ls_iter: number of iterations of the linear solver
• iter_time: current iteration time
• total_time: total optimization time
• Linear search method
• f: The value of the current objective function
• d: Objective function change caused by current parameter change
• g: Module length of current gradient
• h: Parameter variation
• s: Linear search optimal step size
• it: current iteration time
• tt: total optimization time

In the handwritten Gauss Newton method problem, the following configuration is carried out:

//  Options
ceres::Solver::Options options;
//  cholesky decomposition
options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY;
//  Output optimization information
options.minimizer_progress_to_stdout = true;


### Solver::Summary

Solver::Summary is the information of solver and various variables. The common member functions are as follows:

• Brief report (): output a simple summary of a single line;
• FullReport(): output a complete summary of multiple lines.

## Ceres Library: handwritten Gauss Newton method

### source code

The code is as follows:

#include "main.h"
#include <chrono>

/*--------  Build CostFunction--------*/
struct CURVE_FITTING_COST{
//  Constructor
CURVE_FITTING_COST(double x, double y): _x(x), _y(y){}

template<typename T>
bool operator()(const T *const abc, T *residual) const
{
//  y - exp(ax^2 + bx + c)
residual = T(_y) - ceres::exp(abc * T(_x) * T(_x) + abc * T(_x) + abc);
return true;
}

//  Member variables that store values
const double _x, _y;
};

int main()
{
/*--------  Initial parameter configuration--------*/
//  Actual curve parameters
double ar = 1.0, br = 1.0, cr = 1.0;
//  Initial value of estimated curve parameters
double ae = 2.0, be = 1.0, ce = 5.0;
//  Number of sampling observation data points
int N = 100;
//  Noise standard deviation and its reciprocal
double w_sigma = 1.0;
//  Random number generator
cv::RNG rng;

/*--------  Observation data generation--------*/
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 * w_sigma));
}
//  Parameter block
double abc = {ae, be, ce};

/*--------  Constructing least squares problem--------*/
//  instantiation
ceres::Problem problem;
for(int i = 0; i < N; i++){
//  Residual term configuration
problem.AddResidualBlock(new ceres::AutoDiffCostFunction<CURVE_FITTING_COST, 1, 3>(new CURVE_FITTING_COST(x_data[i], y_data[i])),
nullptr,
abc);
}

/*--------  Solving the least squares problem--------*/
//  Options
ceres::Solver::Options options;
//  cholesky decomposition
options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY;
//  Output optimization information
options.minimizer_progress_to_stdout = true;

//  summary
ceres::Solver::Summary summary;
//  Solve start timing t1
//  solve
ceres::Solve(options, &problem, &summary);
//  Solution end timing t2
//  Total solution time
chrono::duration<double> time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
cout << "solve time cost = " << time_used.count() << " seconds. " << endl;

/*--------  Result output--------*/
cout << summary.BriefReport() << endl;
cout << "estimated a,b,c = ";
for(auto a:abc) cout << a << " ";
cout << endl;

return 0;
}


### output

The output contents are as follows:

iter      cost      cost_change  |gradient|   |step|    tr_ratio  tr_radius  ls_iter  iter_time  total_time
0  4.631766e+07    0.00e+00    9.35e+07   0.00e+00   0.00e+00  1.00e+04        0    1.30e-03    1.76e-03
1  6.190334e+06    4.01e+07    1.27e+07   9.68e-01   8.66e-01  1.65e+04        1    1.29e-03    3.80e-03
2  8.093942e+05    5.38e+06    1.73e+06   9.43e-01   8.69e-01  2.76e+04        1    1.03e-03    4.89e-03
3  9.975023e+04    7.10e+05    2.38e+05   8.72e-01   8.77e-01  4.83e+04        1    1.19e-03    6.12e-03
4  1.050054e+04    8.92e+04    3.31e+04   7.14e-01   8.95e-01  9.53e+04        1    1.34e-03    7.51e-03
5  7.599604e+02    9.74e+03    4.54e+03   5.14e-01   9.32e-01  2.69e+05        1    1.13e-03    8.70e-03
6  6.636318e+01    6.94e+02    4.62e+02   5.17e-01   9.78e-01  8.07e+05        1    1.20e-03    9.95e-03
7  5.102078e+01    1.53e+01    1.40e+01   2.57e-01   9.98e-01  2.42e+06        1    1.13e-03    1.11e-02
8  5.100209e+01    1.87e-02    1.87e-02   1.44e-02   9.99e-01  7.26e+06        1    9.99e-04    1.22e-02
solve time cost = 0.0123565 seconds.
Ceres Solver Report: Iterations: 9, Initial cost: 4.631766e+07, Final cost: 5.100209e+01, Termination: CONVERGENCE
estimated a,b,c = 0.877573 1.21245 0.931249


Topics: C++