LM拟合 C++
未完成
#include <iostream>
#include <vector>
#include <array>
#include <ctime>
#include <random>
using namespace std;
void Calc_J_fx(vector<array <double, 2> >& data,
double& k,
vector<double>& Js,
vector<double>& fx)
{
for (int i = 0; i < data.size(); ++i)
{
double x = data[i][0];
double y = data[i][1];
Js[i] = 2 * x;
fx[i] = k * x * x - y;
}
}
void Calc_H_g(vector<double>& Js,
vector<double>& fx,
double& H,
double& g)
{
H = 0;
g = 0;
for (int i = 0; i < data.size(); ++i)
{
H += Js[i] * Js[i];
g += -Js[i] * fx[i]
}
}
void Calc_Cost(vector<double>& fx,
double& cost)
{
cost = 0
for (int i = 0; i < fx.size(); ++i)
{
cost += fx[i] * fx[i]
}
}
void Calc_F(vector<array <double, 2> >& data,
double& k,
double& F)
{
for (int i = 0; i < data.size(); ++i)
{
double x = data[i][0];
double y = data[i][1];
double f = k * x * x - y;
F += 0.5 * f;
}
}
void Calc_L(double& mu,
double& g,
double& dk,
double& dL)
{
dL = 0.5 * dk * (mu * dk - g);
}
int main()
{
int N = 30;
vector< array<double, 2> > data;
std::default_random_engine e;
std::uniform_real_distribution<double> u(0,5);
std::default_random_engine e_;
std::uniform_real_distribution<double> u_(-0.5,0.5);
e.seed(time(0));
for (int i = 0; i < N; ++i)
{
double x = u(e);
double y = k * x * x + double(u_(e_));
data[i][0] = x;
data[i][1] = y;
}
k = 0;
double
int max_iter_num = 200;
double epsilon2 = 1e-6;
double epsilon1 = 1e-6;
vector<double> Js;
vector<double> fx
double H = 0;
double k = 2;
double g = 0;
double cost = 0;
Cala_J_fx(data, k, Js, fx);
Cala_H_g(Js, fx, H, g);
double mu = H;
double vu = 2.0;
data.reserve(N);
Js.reserve(N);
fx.reserve(N);
for (int it = 0; it < max_iter_num; ++it)
{
Cala_J_fx(data, k, Js, fx);
Cala_H_g(Js, fx, H, g);
// 梯度较小
if (sqrt(g * g) < epsilon1)
{
break
}
double dk = 1./H * g;
// 增量少
if (sqrt(dk * dk) <= epsilon2 * sqrt(k * k) + epsilon2)
{
break;
}
else
{
double new_k = k + dk;
double F0, F1, rho;
Calc_F(data, k, F0);
Calc_F(data, new_k, F1);
Calc_L(mu, g, dk, dL);
rho = (F0 - F1)/dL;
if (rho > 0)
{
k = new_k;
mu = max(1.0 / 3.0, 1 - pow(2 * roi - 1, 3));
vu = 2.0;
}
else
{
mu = mu * vu;
vu *= 2;
}
}
}
return 0;
}