网站建设建站经验,网站底部 图标,wordpress引入php文件,营销型网站的目标是文章目录1、Curve Fitting1.1、残差定义1.2、 Problem问题构造1.3、完整代码1.4、运行结果2、Robust Curve Fitting1、Curve Fitting
到目前为止#xff0c;我们看到的示例都是没有数据的简单优化问题。最小二乘和非线性最小二乘分析的原始目的是对数据进行曲线拟合。
以一个…
文章目录1、Curve Fitting1.1、残差定义1.2、 Problem问题构造1.3、完整代码1.4、运行结果2、Robust Curve Fitting1、Curve Fitting
到目前为止我们看到的示例都是没有数据的简单优化问题。最小二乘和非线性最小二乘分析的原始目的是对数据进行曲线拟合。
以一个简单的曲线拟合的问题为例。采样点是根据曲线 ye0.3x0.1y e^{0.3x 0.1}ye0.3x0.1 生成并且添加标准差 σ0.2σ0.2σ0.2 的高斯噪声。我们用下列带未知参数的方程来拟合这些采样点yemxc.y e^{mx c}.yemxc.
1.1、残差定义
首先定义一个模板对象来计算残差。每一个观察值(采样点)都有一个残差
struct ExponentialResidual {ExponentialResidual(double x, double y): x_(x), y_(y) {}template typename Tbool operator()(const T* const m, const T* const c, T* residual) const {residual[0] y_ - exp(m[0] * x_ c[0]);return true;}private:// Observations for a sample.const double x_;const double y_;
};1.2、 Problem问题构造
假设观测数据是一个名为data的2n大小的的数组为每一个观察值创建一个CostFunction的问题(problem)构造是一个简单的事。
double m 0.0;
double c 0.0;Problem problem;
for (int i 0; i kNumObservations; i) {CostFunction* cost_function new AutoDiffCostFunctionExponentialResidual, 1, 1, 1(new ExponentialResidual(data[2 * i], data[2 * i 1]));problem.AddResidualBlock(cost_function, nullptr, m, c);
}
*/与Hello World的f(x)10−x对比
struct CostFunctor {template typename Tbool operator()(const T* const x, T* residual) const {residual[0] T(10.0) - x[0];return true;}
};
CostFunction* cost_function new AutoDiffCostFunctionCostFunctor, 1, 1(new CostFunctor);
problem.AddResidualBlock(cost_function, NULL, x);对比结果:
1.在Hello World中CostFunctor中是没有显式构造函数的也就同样没有了初始值。所以在构造对象时可以直接New CostFunctor。而在本节的例子中构造对象时还要加上初始值即 new ExponentialResidual(data[2 * i], data[2 * i 1]));2.在AutoDiffCostFunction的模板中本例中一共有三个1而在Hello World中只有两个1即residual和x的维度。注意先是残差后是输入参数而且一一对应。
1.3、完整代码
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2015 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS AS IS
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwalgoogle.com (Sameer Agarwal)#include ceres/ceres.h
#include glog/logging.husing ceres::AutoDiffCostFunction;
using ceres::CostFunction;
using ceres::Problem;
using ceres::Solver;
using ceres::Solve;// Data generated using the following octave code.
// randn(seed, 23497);
// m 0.3;
// c 0.1;
// x[0:0.075:5];
// y exp(m * x c);
// noise randn(size(x)) * 0.2;
// y_observed y noise;
// data [x, y_observed];const int kNumObservations 67;
const double data[] {0.000000e00, 1.133898e00,7.500000e-02, 1.334902e00,1.500000e-01, 1.213546e00,2.250000e-01, 1.252016e00,3.000000e-01, 1.392265e00,3.750000e-01, 1.314458e00,4.500000e-01, 1.472541e00,5.250000e-01, 1.536218e00,6.000000e-01, 1.355679e00,6.750000e-01, 1.463566e00,7.500000e-01, 1.490201e00,8.250000e-01, 1.658699e00,9.000000e-01, 1.067574e00,9.750000e-01, 1.464629e00,1.050000e00, 1.402653e00,1.125000e00, 1.713141e00,1.200000e00, 1.527021e00,1.275000e00, 1.702632e00,1.350000e00, 1.423899e00,1.425000e00, 1.543078e00,1.500000e00, 1.664015e00,1.575000e00, 1.732484e00,1.650000e00, 1.543296e00,1.725000e00, 1.959523e00,1.800000e00, 1.685132e00,1.875000e00, 1.951791e00,1.950000e00, 2.095346e00,2.025000e00, 2.361460e00,2.100000e00, 2.169119e00,2.175000e00, 2.061745e00,2.250000e00, 2.178641e00,2.325000e00, 2.104346e00,2.400000e00, 2.584470e00,2.475000e00, 1.914158e00,2.550000e00, 2.368375e00,2.625000e00, 2.686125e00,2.700000e00, 2.712395e00,2.775000e00, 2.499511e00,2.850000e00, 2.558897e00,2.925000e00, 2.309154e00,3.000000e00, 2.869503e00,3.075000e00, 3.116645e00,3.150000e00, 3.094907e00,3.225000e00, 2.471759e00,3.300000e00, 3.017131e00,3.375000e00, 3.232381e00,3.450000e00, 2.944596e00,3.525000e00, 3.385343e00,3.600000e00, 3.199826e00,3.675000e00, 3.423039e00,3.750000e00, 3.621552e00,3.825000e00, 3.559255e00,3.900000e00, 3.530713e00,3.975000e00, 3.561766e00,4.050000e00, 3.544574e00,4.125000e00, 3.867945e00,4.200000e00, 4.049776e00,4.275000e00, 3.885601e00,4.350000e00, 4.110505e00,4.425000e00, 4.345320e00,4.500000e00, 4.161241e00,4.575000e00, 4.363407e00,4.650000e00, 4.161576e00,4.725000e00, 4.619728e00,4.800000e00, 4.737410e00,4.875000e00, 4.727863e00,4.950000e00, 4.669206e00,
};struct ExponentialResidual {ExponentialResidual(double x, double y): x_(x), y_(y) {}template typename T bool operator()(const T* const m,const T* const c,T* residual) const {residual[0] y_ - exp(m[0] * x_ c[0]);return true;}private:const double x_;const double y_;
};int main(int argc, char** argv) {google::InitGoogleLogging(argv[0]);double m 0.0;double c 0.0;Problem problem;for (int i 0; i kNumObservations; i) {problem.AddResidualBlock(new AutoDiffCostFunctionExponentialResidual, 1, 1, 1(new ExponentialResidual(data[2 * i], data[2 * i 1])),NULL,m, c);}Solver::Options options;options.max_num_iterations 25;options.linear_solver_type ceres::DENSE_QR;options.minimizer_progress_to_stdout true;Solver::Summary summary;Solve(options, problem, summary);std::cout summary.BriefReport() \n;std::cout Initial m: 0.0 c: 0.0 \n;std::cout Final m: m c: c \n;return 0;
}1.4、运行结果
iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter iter_time total_time0 1.211734e02 0.00e00 3.61e02 0.00e00 0.00e00 1.00e04 0 1.84e-03 2.30e-031 2.334822e03 -2.21e03 0.00e00 7.52e-01 -1.87e01 5.00e03 1 7.51e-04 3.41e-032 2.331438e03 -2.21e03 0.00e00 7.51e-01 -1.86e01 1.25e03 1 3.35e-04 3.82e-033 2.311313e03 -2.19e03 0.00e00 7.48e-01 -1.85e01 1.56e02 1 3.31e-04 4.22e-034 2.137268e03 -2.02e03 0.00e00 7.22e-01 -1.70e01 9.77e00 1 3.32e-04 4.62e-035 8.553131e02 -7.34e02 0.00e00 5.78e-01 -6.32e00 3.05e-01 1 3.30e-04 5.02e-036 3.306595e01 8.81e01 4.10e02 3.18e-01 1.37e00 9.16e-01 1 1.95e-03 7.04e-037 6.426770e00 2.66e01 1.81e02 1.29e-01 1.10e00 2.75e00 1 2.03e-03 9.14e-038 3.344546e00 3.08e00 5.51e01 3.05e-02 1.03e00 8.24e00 1 4.12e-03 1.34e-029 1.987485e00 1.36e00 2.33e01 8.87e-02 9.94e-01 2.47e01 1 2.04e-03 1.55e-0210 1.211585e00 7.76e-01 8.22e00 1.05e-01 9.89e-01 7.42e01 1 1.96e-03 1.76e-0211 1.063265e00 1.48e-01 1.44e00 6.06e-02 9.97e-01 2.22e02 1 1.97e-03 1.96e-0212 1.056795e00 6.47e-03 1.18e-01 1.47e-02 1.00e00 6.67e02 1 1.97e-03 2.17e-0213 1.056751e00 4.39e-05 3.79e-03 1.28e-03 1.00e00 2.00e03 1 1.96e-03 2.40e-02
Ceres Solver Report: Iterations: 14, Initial cost: 1.211734e02, Final cost: 1.056751e00, Termination: CONVERGENCE
Initial m: 0 c: 0
Final m: 0.291861 c: 0.131439参数的初始值为m0c0初始代价函数为121.173。最后的解是m0.291861c0.131439代价是1.05675。 这个结果和期望解m0.3c0.1有一些细微的差别但是是合理的(因为添加了高斯噪声)。当从噪声数据重建曲线时我们预计会看到这样的偏差。 实际上如果您要对m0.3,c0.1的目标函数进行评估那么当目标函数值为1.082425时拟合效果会更差。下图说明了适合度。
2、Robust Curve Fitting
现在假设我们给出的数据有一些异常值离群点/外点也就是说我们有一些点不服从噪声模型继续使用上面的代码来拟合这些数据我们将得到拟合曲线是偏离实际预期。 处理异常值的标准方法是使用LossFunction。损失函数降低了残差高的残差块(通常是与异常值对应的残差块)的影响。为了将损失函数与残差块联系起来我们改变
// problem.AddResidualBlock(cost_function, nullptr , m, c);problem.AddResidualBlock(cost_function, new CauchyLoss(0.5) , m, c);这里使用SizedCostFunction类需要填写计算jacobians。完整代码如下
#include ceres/ceres.h
#include glog/logging.h// Data generated using the following octave code.
// randn(seed, 23497);
// m 0.3;
// c 0.1;
// x[0:0.075:5];
// y exp(m * x c);
// noise randn(size(x)) * 0.2;
// outlier_noise rand(size(x)) 0.05;
// y_observed y noise outlier_noise;
// data [x, y_observed];const int kNumObservations 67;
const double data[] {
0.000000e00, 1.133898e00,
7.500000e-02, 1.334902e00,
1.500000e-01, 1.213546e00,
2.250000e-01, 1.252016e00,
3.000000e-01, 1.392265e00,
3.750000e-01, 1.314458e00,
4.500000e-01, 1.472541e00,
5.250000e-01, 1.536218e00,
6.000000e-01, 1.355679e00,
6.750000e-01, 1.463566e00,
7.500000e-01, 1.490201e00,
8.250000e-01, 1.658699e00,
9.000000e-01, 1.067574e00,
9.750000e-01, 1.464629e00,
1.050000e00, 1.402653e00,
1.125000e00, 1.713141e00,
1.200000e00, 1.527021e00,
1.275000e00, 1.702632e00,
1.350000e00, 1.423899e00,
1.425000e00, 5.543078e00, // Outlier point
1.500000e00, 5.664015e00, // Outlier point
1.575000e00, 1.732484e00,
1.650000e00, 1.543296e00,
1.725000e00, 1.959523e00,
1.800000e00, 1.685132e00,
1.875000e00, 1.951791e00,
1.950000e00, 2.095346e00,
2.025000e00, 2.361460e00,
2.100000e00, 2.169119e00,
2.175000e00, 2.061745e00,
2.250000e00, 2.178641e00,
2.325000e00, 2.104346e00,
2.400000e00, 2.584470e00,
2.475000e00, 1.914158e00,
2.550000e00, 2.368375e00,
2.625000e00, 2.686125e00,
2.700000e00, 2.712395e00,
2.775000e00, 2.499511e00,
2.850000e00, 2.558897e00,
2.925000e00, 2.309154e00,
3.000000e00, 2.869503e00,
3.075000e00, 3.116645e00,
3.150000e00, 3.094907e00,
3.225000e00, 2.471759e00,
3.300000e00, 3.017131e00,
3.375000e00, 3.232381e00,
3.450000e00, 2.944596e00,
3.525000e00, 3.385343e00,
3.600000e00, 3.199826e00,
3.675000e00, 3.423039e00,
3.750000e00, 3.621552e00,
3.825000e00, 3.559255e00,
3.900000e00, 3.530713e00,
3.975000e00, 3.561766e00,
4.050000e00, 3.544574e00,
4.125000e00, 3.867945e00,
4.200000e00, 4.049776e00,
4.275000e00, 3.885601e00,
4.350000e00, 4.110505e00,
4.425000e00, 4.345320e00,
4.500000e00, 4.161241e00,
4.575000e00, 4.363407e00,
4.650000e00, 4.161576e00,
4.725000e00, 4.619728e00,
4.800000e00, 4.737410e00,
4.875000e00, 4.727863e00,
4.950000e00, 4.669206e00
};using ceres::AutoDiffCostFunction;
using ceres::CostFunction;
using ceres::CauchyLoss;
using ceres::Problem;
using ceres::Solve;
using ceres::Solver;class ExponentialResidual : public ceres::SizedCostFunction1, 2 {
public:ExponentialResidual(double x, double y): x_(x), y_(y) { }virtual bool Evaluate(const double* const* mc,double* residual,double** jacobians) const override{double tmp_y exp(mc[0][0] * x_ mc[0][1]);residual[0] tmp_y - y_; // r exp(mxc) - y if(jacobians jacobians[0]) {jacobians[0][0] x_ * tmp_y; // dr / dmjacobians[0][1] tmp_y; // dr / dc}return true;}private:const double x_, y_;
};int main(int argc, char** argv) {google::InitGoogleLogging(argv[0]);double mc[] {0,0};double m mc[0];double c mc[1];Problem problem;for (int i 0; i kNumObservations; i) {problem.AddResidualBlock(new ExponentialResidual(data[2 * i], data[2 * i 1]),new CauchyLoss(0.5), // 使用损失LossFunction不再是nullptrmc);}Solver::Options options;options.linear_solver_type ceres::DENSE_QR;options.minimizer_progress_to_stdout true;Solver::Summary summary;Solve(options, problem, summary);std::cout summary.BriefReport() \n;std::cout Initial m: 0.0 c: 0.0 \n;std::cout Final m: m c: c \n;return 0;
}CauchyLoss 是Ceres Solver中的一个损失函数参数0.5指定了损失函数的尺度。重新运行程序得到的优化的结果相比不使用损失函数的结果更优。
iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter iter_time total_time0 1.815138e01 0.00e00 2.04e01 0.00e00 0.00e00 1.00e04 0 6.93e-04 1.14e-031 2.259471e01 -4.44e00 0.00e00 5.48e-01 -7.74e-01 5.00e03 1 7.61e-04 2.26e-032 2.258929e01 -4.44e00 0.00e00 5.48e-01 -7.73e-01 1.25e03 1 4.05e-04 2.74e-033 2.255683e01 -4.41e00 0.00e00 5.48e-01 -7.68e-01 1.56e02 1 4.51e-04 5.16e-034 2.225747e01 -4.11e00 0.00e00 5.41e-01 -7.16e-01 9.77e00 1 3.37e-04 5.58e-035 1.784270e01 3.09e-01 8.54e01 4.72e-01 5.44e-02 5.72e00 1 8.33e-04 6.48e-036 7.557353e00 1.03e01 1.13e02 1.06e-01 2.07e00 1.72e01 1 8.29e-04 7.38e-037 2.674796e00 4.88e00 7.69e01 5.41e-02 1.78e00 5.15e01 1 8.22e-04 8.28e-038 1.946177e00 7.29e-01 1.61e01 5.00e-02 1.23e00 1.54e02 1 8.22e-04 9.17e-039 1.904587e00 4.16e-02 2.20e00 2.64e-02 1.11e00 4.63e02 1 8.23e-04 1.01e-0210 1.902929e00 1.66e-03 2.28e-01 6.94e-03 1.12e00 1.39e03 1 8.93e-04 1.10e-0211 1.902884e00 4.51e-05 1.49e-02 1.27e-03 1.14e00 4.17e03 1 9.34e-04 1.21e-02
Ceres Solver Report: Iterations: 12, Initial cost: 1.815138e01, Final cost: 1.902884e00, Termination: CONVERGENCE
Initial m: 0 c: 0
Final m: 0.287605 c: 0.151213 文章转载自： http://www.morning.bqts.cn.gov.cn.bqts.cn http://www.morning.mtrz.cn.gov.cn.mtrz.cn http://www.morning.sjzsjsm.com.gov.cn.sjzsjsm.com http://www.morning.lgnz.cn.gov.cn.lgnz.cn http://www.morning.a3e2r.com.gov.cn.a3e2r.com http://www.morning.spfq.cn.gov.cn.spfq.cn http://www.morning.mnslh.cn.gov.cn.mnslh.cn http://www.morning.rnygs.cn.gov.cn.rnygs.cn http://www.morning.xcdph.cn.gov.cn.xcdph.cn http://www.morning.jqwpw.cn.gov.cn.jqwpw.cn http://www.morning.zkqsc.cn.gov.cn.zkqsc.cn http://www.morning.skqfx.cn.gov.cn.skqfx.cn http://www.morning.ksjnl.cn.gov.cn.ksjnl.cn http://www.morning.zpxwg.cn.gov.cn.zpxwg.cn http://www.morning.hwhnx.cn.gov.cn.hwhnx.cn http://www.morning.bangaw.cn.gov.cn.bangaw.cn http://www.morning.bsgfl.cn.gov.cn.bsgfl.cn http://www.morning.qnzld.cn.gov.cn.qnzld.cn http://www.morning.sskkf.cn.gov.cn.sskkf.cn http://www.morning.ccdyc.cn.gov.cn.ccdyc.cn http://www.morning.fdxhk.cn.gov.cn.fdxhk.cn http://www.morning.ldgqh.cn.gov.cn.ldgqh.cn http://www.morning.ummpdl.cn.gov.cn.ummpdl.cn http://www.morning.qjxkx.cn.gov.cn.qjxkx.cn http://www.morning.jgnjl.cn.gov.cn.jgnjl.cn http://www.morning.sbjhm.cn.gov.cn.sbjhm.cn http://www.morning.dmtld.cn.gov.cn.dmtld.cn http://www.morning.fqljq.cn.gov.cn.fqljq.cn http://www.morning.lqljj.cn.gov.cn.lqljj.cn http://www.morning.fplqh.cn.gov.cn.fplqh.cn http://www.morning.hyfrd.cn.gov.cn.hyfrd.cn http://www.morning.cnvlog.cn.gov.cn.cnvlog.cn http://www.morning.qmsbr.cn.gov.cn.qmsbr.cn http://www.morning.mhcft.cn.gov.cn.mhcft.cn http://www.morning.gydth.cn.gov.cn.gydth.cn http://www.morning.xtrzh.cn.gov.cn.xtrzh.cn http://www.morning.xdjwh.cn.gov.cn.xdjwh.cn http://www.morning.lkkkf.cn.gov.cn.lkkkf.cn http://www.morning.sfdky.cn.gov.cn.sfdky.cn http://www.morning.dwxqf.cn.gov.cn.dwxqf.cn http://www.morning.fydsr.cn.gov.cn.fydsr.cn http://www.morning.nwljj.cn.gov.cn.nwljj.cn http://www.morning.ntwfr.cn.gov.cn.ntwfr.cn http://www.morning.ngcw.cn.gov.cn.ngcw.cn http://www.morning.jmtrq.cn.gov.cn.jmtrq.cn http://www.morning.kpzbf.cn.gov.cn.kpzbf.cn http://www.morning.qjmnl.cn.gov.cn.qjmnl.cn http://www.morning.rjmd.cn.gov.cn.rjmd.cn http://www.morning.nrbcx.cn.gov.cn.nrbcx.cn http://www.morning.yldgw.cn.gov.cn.yldgw.cn http://www.morning.fmgwx.cn.gov.cn.fmgwx.cn http://www.morning.mzhhr.cn.gov.cn.mzhhr.cn http://www.morning.dpbgw.cn.gov.cn.dpbgw.cn http://www.morning.zpnfc.cn.gov.cn.zpnfc.cn http://www.morning.lynmt.cn.gov.cn.lynmt.cn http://www.morning.brkc.cn.gov.cn.brkc.cn http://www.morning.nckzt.cn.gov.cn.nckzt.cn http://www.morning.qysnd.cn.gov.cn.qysnd.cn http://www.morning.lxfdh.cn.gov.cn.lxfdh.cn http://www.morning.ydmml.cn.gov.cn.ydmml.cn http://www.morning.sgfgz.cn.gov.cn.sgfgz.cn http://www.morning.pgcmz.cn.gov.cn.pgcmz.cn http://www.morning.myzfz.com.gov.cn.myzfz.com http://www.morning.bpmfq.cn.gov.cn.bpmfq.cn http://www.morning.btmwd.cn.gov.cn.btmwd.cn http://www.morning.wdrxh.cn.gov.cn.wdrxh.cn http://www.morning.pyswr.cn.gov.cn.pyswr.cn http://www.morning.tthmg.cn.gov.cn.tthmg.cn http://www.morning.gwxwl.cn.gov.cn.gwxwl.cn http://www.morning.srrzb.cn.gov.cn.srrzb.cn http://www.morning.sjwzl.cn.gov.cn.sjwzl.cn http://www.morning.pdghl.cn.gov.cn.pdghl.cn http://www.morning.qlkzl.cn.gov.cn.qlkzl.cn http://www.morning.pjyrl.cn.gov.cn.pjyrl.cn http://www.morning.cfqyx.cn.gov.cn.cfqyx.cn http://www.morning.pggkr.cn.gov.cn.pggkr.cn http://www.morning.rpth.cn.gov.cn.rpth.cn http://www.morning.tgfjm.cn.gov.cn.tgfjm.cn http://www.morning.yhywr.cn.gov.cn.yhywr.cn http://www.morning.fslxc.cn.gov.cn.fslxc.cn
