首页 > 其他分享 >cereas学习(2) 鲍威尔函数 多项式

cereas学习(2) 鲍威尔函数 多项式

时间:2023-01-03 00:55:07浏览次数:68  
标签:03 const 鲍威尔 cereas 多项式 x2 x3 x1 x4

 

 

struct F4 {
  template <typename T>
  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;
  }
};

  

 

  同样,我们可以定义类F1F2F3评估 f1(x1,x2),f2(x3,x4)和f3(x2,x3) 分别。

struct F1 {
  template <typename T>
  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 <typename T>
  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 <typename T>
  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 <typename T>
  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;
  }
};

  

 

使用这些,问题可以构造如下:

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<F1, 1, 1, 1>(new F1), nullptr, &x1, &x2);
problem.AddResidualBlock(
  new AutoDiffCostFunction<F2, 1, 1, 1>(new F2), nullptr, &x3, &x4);
problem.AddResidualBlock(
  new AutoDiffCostFunction<F3, 1, 1, 1>(new F3), nullptr, &x2, &x3);
problem.AddResidualBlock(
  new AutoDiffCostFunction<F4, 1, 1, 1>(new F4), nullptr, &x1, &x4);

  请注意,每个ResidualBlock仅取决于相应残差对象所依赖的两个参数,而不是所有四个参数。编译和运行examples/powell.cc 给我们:

 

Initial x1 = 3, x2 = -1, x3 = 0, x4 = 1
iter      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    4.95e-04    2.30e-03
   1  5.036190e+00    1.02e+02    2.00e+01   2.16e+00   9.53e-01  3.00e+04       1    4.39e-05    2.40e-03
   2  3.148168e-01    4.72e+00    2.50e+00   6.23e-01   9.37e-01  9.00e+04       1    9.06e-06    2.43e-03
   3  1.967760e-02    2.95e-01    3.13e-01   3.08e-01   9.37e-01  2.70e+05       1    8.11e-06    2.45e-03
   4  1.229900e-03    1.84e-02    3.91e-02   1.54e-01   9.37e-01  8.10e+05       1    6.91e-06    2.48e-03
   5  7.687123e-05    1.15e-03    4.89e-03   7.69e-02   9.37e-01  2.43e+06       1    7.87e-06    2.50e-03
   6  4.804625e-06    7.21e-05    6.11e-04   3.85e-02   9.37e-01  7.29e+06       1    5.96e-06    2.52e-03
   7  3.003028e-07    4.50e-06    7.64e-05   1.92e-02   9.37e-01  2.19e+07       1    5.96e-06    2.55e-03
   8  1.877006e-08    2.82e-07    9.54e-06   9.62e-03   9.37e-01  6.56e+07       1    5.96e-06    2.57e-03
   9  1.173223e-09    1.76e-08    1.19e-06   4.81e-03   9.37e-01  1.97e+08       1    7.87e-06    2.60e-03
  10  7.333425e-11    1.10e-09    1.49e-07   2.40e-03   9.37e-01  5.90e+08       1    6.20e-06    2.63e-03
  11  4.584044e-12    6.88e-11    1.86e-08   1.20e-03   9.37e-01  1.77e+09       1    6.91e-06    2.65e-03
  12  2.865573e-13    4.30e-12    2.33e-09   6.02e-04   9.37e-01  5.31e+09       1    5.96e-06    2.67e-03
  13  1.791438e-14    2.69e-13    2.91e-10   3.01e-04   9.37e-01  1.59e+10       1    7.15e-06    2.69e-03

Ceres Solver v1.12.0 Solve Report
----------------------------------
                                     Original                  Reduced
Parameter blocks                            4                        4
Parameters                                  4                        4
Residual blocks                             4                        4
Residual                                    4                        4

Minimizer                        TRUST_REGION

Dense linear algebra library            EIGEN
Trust region strategy     LEVENBERG_MARQUARDT

                                        Given                     Used
Linear solver                        DENSE_QR                 DENSE_QR
Threads                                     1                        1
Linear solver threads                       1                        1

Cost:
Initial                          1.075000e+02
Final                            1.791438e-14
Change                           1.075000e+02

Minimizer iterations                       14
Successful steps                           14
Unsuccessful steps                          0

Time (in seconds):
Preprocessor                            0.002

  Residual evaluation                   0.000
  Jacobian evaluation                   0.000
  Linear solver                         0.000
Minimizer                               0.001

Postprocessor                           0.000
Total                                   0.005

Termination:                      CONVERGENCE (Gradient tolerance reached. Gradient max norm: 3.642190e-11 <= 1.000000e-10)

Final x1 = 0.000292189, x2 = -2.92189e-05, x3 = 4.79511e-05, x4 = 4.79511e-05

  

 

 

//
// An example program that minimizes Powell's singular function.
//
//   F = 1/2 (f1^2 + f2^2 + f3^2 + f4^2)
//
//   f1 = x1 + 10*x2;
//   f2 = sqrt(5) * (x3 - x4)
//   f3 = (x2 - 2*x3)^2
//   f4 = sqrt(10) * (x1 - x4)^2
//
// The starting values are x1 = 3, x2 = -1, x3 = 0, x4 = 1.
// The minimum is 0 at (x1, x2, x3, x4) = 0.
//
// From: Testing Unconstrained Optimization Software by Jorge J. More, Burton S.
// Garbow and Kenneth E. Hillstrom in ACM Transactions on Mathematical Software,
// Vol 7(1), March 1981.
#include <vector>
#include "ceres/ceres.h"
#include "gflags/gflags.h"
#include "glog/logging.h"
using ceres::AutoDiffCostFunction;
using ceres::CostFunction;
using ceres::Problem;
using ceres::Solve;
using ceres::Solver;
struct F1 {
  template <typename T>
  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 <typename T>
  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 <typename T>
  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 <typename T>
  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<F1, 1, 1, 1>(new F1), nullptr, &x1, &x2);
  problem.AddResidualBlock(
      new AutoDiffCostFunction<F2, 1, 1, 1>(new F2), nullptr, &x3, &x4);
  problem.AddResidualBlock(
      new AutoDiffCostFunction<F3, 1, 1, 1>(new F3), nullptr, &x2, &x3);
  problem.AddResidualBlock(
      new AutoDiffCostFunction<F4, 1, 1, 1>(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;
}

  

标签:03,const,鲍威尔,cereas,多项式,x2,x3,x1,x4
From: https://www.cnblogs.com/gooutlook/p/17020928.html

相关文章

  • 多项式乘法学习笔记
    多项式乘法给定两个多项式\(A(x)=\sum_{i=0}^{n-1}{a_ix^i}\)\(B(x)=\sum^{m-1}_{i=0}{b_ix^i}\)求\(C(x)=A(x)\timesB(x)=\sum_{i=0}^{n+m-1}{c_ix^i}\)......
  • HDU 6801 Game on a Circle 题解 (推式子,多项式)
    题目链接首先注意到我们对这个环的扫描是一轮一轮进行的,每轮都会从左到右对每个没被删除的元素以p的概率删除。如果我们能对每个\(t(t\in[0,\infin],t是整数)和i\)求出c......
  • 多项式笔记
    多项式,清空,封装,边界,是讨厌的。NTT懒得讲原理,喵。 llA[maxn],B[maxn],pos[maxn],S[maxn]; constintmod=998244353,g=3,ginv=332748118; llksm(llx,inty) { ll......
  • [NOIP2009 普及组] 多项式输出
    [NOIP2009普及组]多项式输出题目描述一元$n$次多项式可用如下的表达式表示:$$f(x)=a_nxn+a_{n-1}x{n-1}+\cdots+a_1x+a_0,a_n\ne0$$其中,$a_ix^i$称为$i$次项,$a......
  • 省选07. 多项式
    P3338[ZJOI2014]力\[\begin{aligned}E_i&=\sum_{j=1}^{i-1}\frac{q_j}{(i-j)^2}-\sum_{j=i+1}^n\frac{q_j}{(i-j)^2}\end{aligned}\]设\(f(x)=q_x\),\(g(x)=x^2\),\(h(......
  • luoguP5383 普通多项式转下降幂多项式 题解
    学习了bztMinamoto大佬的做法,希望这篇题解可以使得那个方法更加易于理解。既然下降幂多项式转普通多项式可以采取分治\(\operatorname{NTT}\),那么可以猜测逆过来也可以......
  • 【顺序结构】计算多项式的值
    前言在上一次发布【顺序结构】甲流疫情死亡率后,今天我们继续来感受一下【顺序结构】计算多项式的值的做法。正文题目描述:对于多项式f(x)=ax3+bx2+cx+d和给定的x,a,b,......
  • 下标模意义下的多项式乘法及其应用
    已知两个数组\(a,b\),求一个数组\(c\),满足\(c_i=\sum_{j+k\equivi(MOD\K)}a_jb_k(MOD\M)\)。这里我们把这个东西称为下标模意义下的多项式乘法。那么这个东西怎么做呢?先......
  • 尝试把多项式板子分段喂给 ChatGPT 辨认
    #include<bits/stdc++.h>#defineendl'\n'#definerep(i,a,b)for(inti=(a);i<=(b);++i)#defineRep(i,a,b)for(inti=(a);i>=(b);--i)usingnamespacestd;const......
  • 【视频】什么是非线性模型与R语言多项式回归、局部平滑样条、 广义相加GAM分析工资数
    最近我们被客户要求撰写关于非线性模型的研究报告,包括一些图形和统计输出。在这文中,我将介绍非线性回归的基础知识。非线性回归是一种对因变量和一组自变量之间的非线性关系......