Ceres简单应用-求解(Powell's Function)鲍威尔函数最小值

【资料图】

Ceres 求解 Powell’s function 的最小化

\(\quad\)现在考虑一个稍微复杂一点的例子—鲍威尔函数的最小化。

\(\quad{}\) \(x=[x_1,x_2,x_3,x_4]\) 并且

\[\begin{array}{l}f_{1}(x)=x_{1}+10 x_{2} \\f_{2}(x)=\sqrt{5}\left(x_{3}-x_{4}\right) \\f_{3}(x)=\left(x_{2}-2 x_{3}\right)^{2} \\f_{4}(x)=\sqrt{10}\left(x_{1}-x_{4}\right)^{2} \\ \\ F(x)=\left[f_{1}(x), f_{2}(x), f_{3}(x), f_{4}(x)\right]\end{array}\]

\(\quad{}\) \(F(x)\) 是一个拥有四个参数的函数，并且拥有四个残差块，我们希望找到一个 \(x\) 使得下面的式子得到最小化：

\[\frac{1}{2}||F(x)||^2\]

\(\quad{}\)同样，第一步是定义对目标函子中的项进行评估的函子。 这是评估的代码，翻译过来不太准确，或者说这其实是在定义残差块。对 \(f_4(x_1,x_4)\) 的定义如下：

struct F4 {  template   bool operator()(const T* const x1, const T* const x4, T* residual) const {    residual[0] = sqrt(10.0) * (x1[0] - x4[0]) * (x1[0] - x4[0]);    return true;  }};

\(\quad{}\) 同样，我们也可以类似的定义 \(F_1\), \(F_2\) 和 \(F_3\) 去评估计算 \(f_1(x_1,x_2)\)，\(f_2(x_3,x_4)\)，\(f_3(x_2,x_3)\)。这样的话，问题将会被定义为如下的格式：

// 赋予初值double x1 =  3.0; double x2 = -1.0; double x3 =  0.0; double x4 = 1.0;Problem problem;// Add residual terms to the problem using the autodiff// wrapper to get the derivatives automatically.problem.AddResidualBlock(  new AutoDiffCostFunction(new F1), nullptr, &x1, &x2);problem.AddResidualBlock(  new AutoDiffCostFunction(new F2), nullptr, &x3, &x4);problem.AddResidualBlock(  new AutoDiffCostFunction(new F3), nullptr, &x2, &x3);problem.AddResidualBlock(  new AutoDiffCostFunction(new F4), nullptr, &x1, &x4);

\(\quad{}\) 请注意，每个 ResidualBlock仅依赖于对应残差对象所依赖的两个参数，而不依赖于所有四个参数。

所有的程序如下所示：

#include #include "ceres/ceres.h"#include "gflags/gflags.h"#include "glog/logging.h"#include "ceres/internal/port.h"using ceres::AutoDiffCostFunction;using ceres::CostFunction;using ceres::Problem;using ceres::Solve;using ceres::Solver;struct F1{  template   bool operator()(const T *const x1, const T *const x2, T *residual) const  {    // f1 = x1 + 10 * x2;    residual[0] = x1[0] + 10.0 * x2[0];    return true;  }};struct F2{  template   bool operator()(const T *const x3, const T *const x4, T *residual) const  {    // f2 = sqrt(5) (x3 - x4)    residual[0] = sqrt(5.0) * (x3[0] - x4[0]);    return true;  }};struct F3{  template   bool operator()(const T *const x2, const T *const x3, T *residual) const  {    // f3 = (x2 - 2 x3)^2    residual[0] = (x2[0] - 2.0 * x3[0]) * (x2[0] - 2.0 * x3[0]);    return true;  }};struct F4{  template   bool operator()(const T *const x1, const T *const x4, T *residual) const  {    // f4 = sqrt(10) (x1 - x4)^2    residual[0] = sqrt(10.0) * (x1[0] - x4[0]) * (x1[0] - x4[0]);    return true;  }};DEFINE_string(minimizer,              "trust_region",              "Minimizer type to use, choices are: line_search & trust_region");int main(int argc, char **argv){  GFLAGS_NAMESPACE::ParseCommandLineFlags(&argc, &argv, true);  google::InitGoogleLogging(argv[0]);  double x1 = 3.0;  double x2 = -1.0;  double x3 = 0.0;  double x4 = 1.0;  Problem problem;  // Add residual terms to the problem using the autodiff  // wrapper to get the derivatives automatically. The parameters, x1 through  // x4, are modified in place.  problem.AddResidualBlock(      new AutoDiffCostFunction(new F1), nullptr, &x1, &x2);  problem.AddResidualBlock(      new AutoDiffCostFunction(new F2), nullptr, &x3, &x4);  problem.AddResidualBlock(      new AutoDiffCostFunction(new F3), nullptr, &x2, &x3);  problem.AddResidualBlock(      new AutoDiffCostFunction(new F4), nullptr, &x1, &x4);  Solver::Options options;  // LOG_IF(FATAL,  //        !ceres::StringToMinimizerType(CERES_GET_FLAG(FLAGS_minimizer),  //                                      &options.minimizer_type))  //     << "Invalid minimizer: " << CERES_GET_FLAG(FLAGS_minimizer)  //     << ", valid options are: trust_region and line_search.";  options.max_num_iterations = 100;  options.linear_solver_type = ceres::DENSE_QR;  options.minimizer_progress_to_stdout = true;  // clang-format off  std::cout << "Initial x1 = " << x1            << ", x2 = " << x2            << ", x3 = " << x3            << ", x4 = " << x4            << "\n";  // clang-format on  // Run the solver!  Solver::Summary summary;  Solve(options, &problem, &summary);  std::cout << summary.FullReport() << "\n";  // clang-format off  std::cout << "Final x1 = " << x1            << ", x2 = " << x2            << ", x3 = " << x3            << ", x4 = " << x4            << "\n";  // clang-format on  return 0;}

\(\quad{}\) 编译并运行可得到如下输出：

Initial x1 = 3, x2 = -1, x3 = 0, x4 = 1iter      cost      cost_change  |gradient|   |step|    tr_ratio  tr_radius  ls_iter  iter_time  total_time   0  1.075000e+02    0.00e+00    1.55e+02   0.00e+00   0.00e+00  1.00e+04        0    5.89e-05    5.71e-04   1  5.036190e+00    1.02e+02    2.00e+01   2.16e+00   9.53e-01  3.00e+04        1    1.80e-04    7.80e-04   2  3.148168e-01    4.72e+00    2.50e+00   6.23e-01   9.37e-01  9.00e+04        1    3.90e-05    8.32e-04   3  1.967760e-02    2.95e-01    3.13e-01   3.08e-01   9.37e-01  2.70e+05        1    3.65e-05    8.79e-04   4  1.229900e-03    1.84e-02    3.91e-02   1.54e-01   9.37e-01  8.10e+05        1    3.57e-05    9.23e-04   5  7.687123e-05    1.15e-03    4.89e-03   7.69e-02   9.37e-01  2.43e+06        1    3.54e-05    9.67e-04   6  4.804625e-06    7.21e-05    6.11e-04   3.85e-02   9.37e-01  7.29e+06        1    3.52e-05    1.01e-03   7  3.003028e-07    4.50e-06    7.64e-05   1.92e-02   9.37e-01  2.19e+07        1    3.53e-05    1.05e-03   8  1.877006e-08    2.82e-07    9.54e-06   9.62e-03   9.37e-01  6.56e+07        1    3.53e-05    1.10e-03   9  1.173223e-09    1.76e-08    1.19e-06   4.81e-03   9.37e-01  1.97e+08        1    3.58e-05    1.15e-03  10  7.333425e-11    1.10e-09    1.49e-07   2.40e-03   9.37e-01  5.90e+08        1    3.54e-05    1.19e-03  11  4.584044e-12    6.88e-11    1.86e-08   1.20e-03   9.37e-01  1.77e+09        1    3.51e-05    1.24e-03  12  2.865573e-13    4.30e-12    2.33e-09   6.02e-04   9.37e-01  5.31e+09        1    3.55e-05    1.28e-03  13  1.791438e-14    2.69e-13    2.91e-10   3.01e-04   9.37e-01  1.59e+10        1    3.55e-05    1.32e-03  14  1.120029e-15    1.68e-14    3.64e-11   1.51e-04   9.37e-01  4.78e+10        1    3.54e-05    1.37e-03Solver Summary (v 1.13.0-eigen-(3.3.4)-lapack-suitesparse-(5.1.2)-cxsparse-(3.1.9)-openmp)                                     Original                  ReducedParameter blocks                            4                        4Parameters                                  4                        4Residual blocks                             4                        4Residual                                    4                        4Minimizer                        TRUST_REGIONDense linear algebra library            EIGENTrust region strategy     LEVENBERG_MARQUARDT                                        Given                     UsedLinear solver                        DENSE_QR                 DENSE_QRThreads                                     1                        1Linear solver threads                       1                        1Linear solver ordering              AUTOMATIC                        4Cost:Initial                          1.075000e+02Final                            1.120029e-15Change                           1.075000e+02Minimizer iterations                       15Successful steps                           15Unsuccessful steps                          0Time (in seconds):Preprocessor                           0.0005  Residual evaluation                  0.0000  Jacobian evaluation                  0.0005  Linear solver                        0.0001Minimizer                              0.0009Postprocessor                          0.0000Total                                  0.0014Termination:                      CONVERGENCE (Gradient tolerance reached. Gradient max norm: 3.642190e-11 <= 1.000000e-10)Final x1 = 0.000146222, x2 = -1.46222e-05, x3 = 2.40957e-05, x4 = 2.40957e-05

\(\quad{}\) 其实我们很容易看出函数的最优解在 \(x_1=0, x_2=0,x_3=0,x_4=0\) ，而最终的计算结果虽然并不是完全为\(0\) 但是也基本接近。