Gauss-Newton Method (c + + implementation)

#include <iostream> #include <cmath> #include <Eigen/Eigen> using namespace Eigen; #define ITER_STEP (1e-5) #define ITER_CNT () typedef void (*FUNC_PTR) (const VECTORXD &input, const VECTORXD &

Amp;params, Vectorxd &output); void People_func (const vectorxd &input, const vectorxd &params, Vectorxd &output) {Double A = params (0, 0
    );

    Double B = params (1, 0);
    for (int i = 0; i < input.rows (); i++) {output (i, 0) = A * EXP (input (i, 0) * B); } void Get_jacobian (func_ptr func, const VECTORXD &input, const VECTORXD &params, Matri
    Xxd &output) {int m = input.rows ();

    int n = params.rows ();
    Vectorxd out0 (M, 1);
    Vectorxd OUT1 (M, 1); 
    Vectorxd param0 (n, 1); 

    Vectorxd param1 (n, 1);

    Output.resize (M, n);
        for (int j = 0; J < N; j +) {param0 = params;
        param1 = params;
        Param0 (j, 0)-= iter_step; Param1 (j, 0) += Iter_step;
        Func (input, Param0, out0);
        Func (input, param1, OUT1);
    Output.block (0, J, m, 1) = (out1-out0)/(2 * iter_step); } void Gauss_newton (func_ptr func, const VECTORXD &inputs, const VECTORXD &output, Vect
    Orxd &params) {int m = inputs.rows ();

    int n = params.rows ();
    Jacobian matrixxd Jmat (m, n);
    Vectorxd R (M, 1);

    Vectorxd tmp (M, 1);
    Double Pre_mse = 0.0;

    Double MSE;
        for (int i = 0; i < iter_cnt i++) {MSE = 0.0;
        Func (inputs, params, TMP);
        R = output-tmp;
        Get_jacobian (func, inputs, params, jmat);
        MSE = R.transpose () * r;
        MSE/= m;
        if (Fabs (MSE-PRE_MSE) < 1e-8) {break;
        } PRE_MSE = MSE;
        Vectorxd Delta = (jmat.transpose () * jmat). Inverse () * Jmat.transpose () * r;
        printf ("i =%d, MSE%lf.\n", I, MSE);
    params = Delta; } std::cout <<


"Params:" << params.transpose () << Std::endl;}
    int People_fit () {//A * EXP (B * x) VECTORXD Inputs (8, 1);
    Inputs << 1, 2, 3, 4, 5, 6, 7, 8;
    Vectorxd output (8, 1);

    Output << 8.3, 11.0, 14.7, 19.7, 26.7, 35.2, 44.4, 55.9;
    Vectorxd params (2, 1);
    Params << 8, 0.7;
    Gauss_newton (people_func, inputs, output, params);
return 0; } void Tri_func (const vectorxd &input, const vectorxd &params, Vectorxd &output) {Double A = params (0, 0
    );
    Double B = params (1, 0);
    Double C = params (2, 0);

    Double D = params (3, 0);
    for (int i = 0; i < input.rows (); i++) {output (i, 0) = A*sin (B*input (i, 0)) + c*cos (d*input (i, 0));
    } int Tri_func_fit () {//A * sin (Bx) + C * cos (Dx) Double A = 5;
    Double B = 1;
    Double C = 10;

    Double D = 2;
    const int smp_cnt = 100;
    Vectorxd input (smp_cnt, 1);
    Vectorxd output (smp_cnt, 1);

    VECTORXD params (4, 1); Params &LT;&LT

    1, 1, 9, 1;
        for (int i = 0; i < smp_cnt i++) {input (i, 0) = i;
    Output (i, 0) = A*sin (B*input (i, 0)) +c*cos (D*input (i, 0)) + (rand ()% 255)/255.0;
    } Gauss_newton (Tri_func, input, output, params);
return 0;
    int main () {people_fit ();
Tri_func_fit ();

 }