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;
}