Nonlinear least squares fitting

stl · **Joined:** Fri Apr 01, 2016 9:09 pm **Posts:** 5

Hello,
i am implementing curve fitting in C# and i tried curve fitting by function values, gradient and i am trying Hessian right now. But i am not really sure about Hessian values, I am not so good in math and this is the first time for me implementing some huge math problem.

My function is:
y = s*(1-exp(-a*pow(x,n)))
where s, a, n are unknown variables which i am looking for.

Could someone check my Hessian, please? I hope it is correct, but i need to be sure. Also i have question about tweaking the calculation, by which parameters i can influence precision of the results? Thank you very much!

Code:

 private void button2_Click(object sender, EventArgs e)
        {
            double[,] x = new double[,] { { 0.0000 }, { 0.4900 }, { 1.0000 }, { 1.4700 }, { 1.9600 }, { 2.4500 }, { 2.9400 }, { 3.4300 }, { 3.9200 }, { 4.4100 }, { 4.9000 } };
            double[] y = new double[] { 0.0000, 0.1003, 0.1375, 0.1581, 0.1743, 0.1856, 0.1952, 0.2027, 0.2096, 0.2147, 0.2213, };
            double[] c = new double[] { 0.2, 0.5, 0.7 };
            double epsf = 0;
            double epsx = 0.000001;
            int maxits = 0;
            int info;
            alglib.lsfitstate state;
            alglib.lsfitreport rep;
            double diffstep = 0.0001;
            // Fitting without weights
            //            X       -   array[0..N-1,0..M-1], points (one row = one point)
            //            Y       -   array[0..N-1], function values.
            //            C       -   array[0..K-1], initial approximation to the solution,
            //            N       -   number of points, N>1
            //            M       -   dimension of space
            //            K       -   number of parameters being fitted
            //            DiffStep-   numerical differentiation step; should not be very small or large; large = loss of accuracy small = growth of round-off errors
            alglib.lsfitcreatef(x, y, c, diffstep, out state);
            alglib.lsfitsetcond(state, epsf, epsx, maxits);
            //alglib.lsfitfit(state, function_cx_1_func, null, null);   // using function values only
            //alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, null, null);   // Gradient
            alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, function_cx_1_hess, null, null); // Gradient + Hessian            
            alglib.lsfitresults(state, out info, out c, out rep);

            label2.Text = alglib.ap.format(c, 7);
            label3.Text = info.ToString();
        }

        public static void function_cx_1_func(double[] c, double[] x, ref double func, object obj)
        {
            ////c[0] == s 
            ////c[1] == a
            ////c[2] == n
            ////y = (s (1-exp(-a x^n)))
            //func = (1 - Math.Exp(-c[0] * Math.Pow(x[0], c[1]))); //(1-exp(-a x^n))
            func = c[0] * (1 - Math.Exp(-c[1] * Math.Pow(x[0], c[2]))); // s*(1-exp(-a*pow(x,n)))
        }

        public static void function_cx_1_grad(double[] c, double[] x, ref double func, double[] grad, object obj)
        {
            ////c[0] == s 
            ////c[1] == a
            ////c[2] == n
            ////(d)/(ds) (s (1-exp(-a x^n))) = 1-e^(-a x^n)
            ////(d)/(da) (s (1-exp(-a x^n))) = s x^n e^(-a x^n)
            ////(d)/(dn) (s (1-exp(-a x^n))) = a s x^n log(x) e^(-a x^n)
            grad[0] = 1 - Math.Exp(-c[1] * Math.Pow(x[0],c[2]));//1-e^(-a x^n)
            grad[1] = c[0] * Math.Pow(x[0],c[2]) * Math.Exp(-c[1] * Math.Pow(x[0],c[2])); //s x^n e^(-a x^n)
            grad[2] = c[1] * c[0] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2])); //a s x^n log(x) e^(-a x^n)
        
        }

        public static void function_cx_1_hess(double[] c, double[] x, ref double func, double[] grad, double[,] hess, object obj)
        {
            ////c[0] == s 
            ////c[1] == a
            ////c[2] == n
            ////(d^2)/(ds^2) (s (1-exp(-a x^n))) = 0
            ////(d^2)/(da^2) (s (1-exp(-a x^n))) = s x^(2 n) (-e^(-a x^n))
            ////(d^2)/(dn^2) (s (1-exp(-a x^n))) = -a s x^n log^2(x) e^(-a x^n) (a x^n-1)
            func = c[0] * (1 - Math.Exp(-c[1] * Math.Pow(x[0], c[2])));
            grad[0] = 1 - Math.Exp(-c[1] * Math.Pow(x[0], c[2]));//1-e^(-a x^n)
            grad[1] = c[0] * Math.Pow(x[0], c[2]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2])); //s x^n e^(-a x^n)
            grad[2] = c[1] * c[0] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2])); //a s x^n log(x) e^(-a x^n)
            // [c0,c0]...[c0,c2]
            //(d^2)/(ds^2) (s (1-exp(-a x^n))) = 0
            //(d)/(ds)(d)/(da) (s (1-exp(-a x^n))) = x^n e^(-a x^n)
            //(d)/(ds)(d)/(dn) (s (1-exp(-a x^n))) = a x^n log(x) e^(-a x^n)
            hess[0,0] = 0;
            hess[0, 1] = Math.Pow(x[0], c[2]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2]));
            hess[0, 2] = c[1] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2]));
            //[c1,c0]...[c1,c2]
            //(d)/(da)(d)/(ds) (s (1-exp(-a x^n))) = x^n e^(-a x^n)
            //(d^2)/(da^2) (s (1-exp(-a x^n))) = s x^(2 n) (-e^(-a x^n))
            //(d)/(da)(d)/(dn) (s (1-exp(-a x^n))) = s x^n log(x) (-e^(-a x^n)) (a x^n-1)
            hess[1, 0] = Math.Pow(x[0], c[2]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2]));
            hess[1, 1] = c[0] * Math.Pow(x[0],(2 * c[2])) * (-Math.Exp(-c[1] * Math.Pow(x[0], c[2])));
            hess[1, 2] = c[0] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * (-Math.Exp(-c[1] * Math.Pow(x[0], c[2]))) * (c[1] * Math.Pow(x[0],(c[2] - 1)));
            //[c2,c0]...[c2,c2]
            //(d)/(dn)(d)/(ds) (s (1-exp(-a x^n))) = a x^n log(x) e^(-a x^n)
            //(d)/(dn)(d)/(da) (s (1-exp(-a x^n))) = s x^n log(x) (-e^(-a x^n)) (a x^n-1)
            //(d^2)/(dn^2) (s (1-exp(-a x^n))) = -a s x^n log^2(x) e^(-a x^n) (a x^n-1)
            hess[2, 0] = c[1] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2]));
            hess[2, 1] = c[0] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * (-Math.Exp(-c[1] * Math.Pow(x[0], c[2]))) * (c[1] * Math.Pow(x[0], (c[2] - 1)));
            hess[2, 2] = -c[1] * c[0] * Math.Pow(x[0], c[2]) * Math.Pow(Math.Log10(x[0]),2) * (-Math.Exp(-c[1] * Math.Pow(x[0], c[2]))) * (c[1] * Math.Pow(x[0], (c[2] - 1)));
        }

Thank you for this library, I found it as the easiest to understand and implement unlike other ones. So as a student i am considering small donation to this project (via paypal if it is possible).

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 930

Hello!

First, your code does NOT perform Hessian-based optimization :) you created optimizer object with lsfitcreatef function call, so it works in "function value + numerical differentiation" mode. Even gradient values are not used. You should call some other constructor function (lsfitcreatefg for gradient-based optimization, lsfitcreatefgh for Hessian-based optimization).

Second, it is hard to check your Hessian analytically - it is really huge. However, you can easily perform numerical check. Just calculate Hessian matrix using your code at some arbitrary point (use non-trivial values for coefficients). And then calculate numerical Hessian (use one of numerical differentiation formulae, some of the basic ones), and compare results. If they are dissimilar, then you have an error. If they are similar... well, either everything is correct - or you made "even number of mistakes" - choose a point where incorrect formula gives correct results purely by chance alone.

As for the donation... I really appreciate your intentions! But I think the best donation would be to recommend library to someone else in a search for numerical software :) do it next time you see someone trying to solve some numerical problem, and everyone will be happy :)

stl · **Joined:** Fri Apr 01, 2016 9:09 pm **Posts:** 5

Thank you for help, i forgot about other constructors!
So i modified code and tried to run gradient and hessian fitting and i am getting back same array as initial "approximation to the solution" (c). I found that there was no iterations in algorithm (rep.iterationscount) and as info is "4 gradient norm is no more than EpsG". I checked gradients and they should be correct, i have also searched in examples and manual, but i have no clue how to solve this. What am i missing? Thanks.

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 930

I just noticed that your gradient calculation code does not set func parameter. It should calculate function value + gradient, but it calculates just gradient. It confuses optimizer. Try fixing it. Does it solve your problem?

stl · **Joined:** Fri Apr 01, 2016 9:09 pm **Posts:** 5

Yes, i modified gradient calculation while i was adding new constructors for gradient and hessian. I have tried different values for epsf, epsx, maxits and change CheapFG to true or false. I am still getting same results and same info state (4). Actual code (hessian calculation is not included) is:

Code:

private double[] LeastSquaresFitting()
        {
            double[,] x = new double[,] { { 0.0000 }, { 0.4900 }, { 1.0000 }, { 1.4700 }, { 1.9600 }, { 2.4500 }, { 2.9400 }, { 3.4300 }, { 3.9200 }, { 4.4100 }, { 4.9000 } };
            double[] y = new double[] { 0.0000, 0.1003, 0.1375, 0.1581, 0.1743, 0.1856, 0.1952, 0.2027, 0.2096, 0.2147, 0.2213, };
            double[] c = new double[] { 0.2, 0.5, 0.7 };

            double epsf = 0;
            double epsx = 0.000001; //0.000001;
            int maxits = 0;
            int info;
            alglib.lsfitstate state;
            alglib.lsfitreport rep;
            double diffstep = 0.0001;
            c = new double[] { 0.2, 0.5, 0.7 };
            int Mode = 2; //1- function values, 2- gradient, 3-hessian
            switch (Mode)
            {
                case 1:
                    alglib.lsfitcreatef(x, y, c, diffstep, out state);
                    alglib.lsfitfit(state, function_cx_1_func, null, null);
                    break;
                case 2:
                    alglib.lsfitcreatefg(x, y, c, true, out state);
                    alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, null, null); 
                    break;
                case 3:
                    alglib.lsfitcreatefgh(x, y, c, out state);
                    alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, function_cx_1_hess, null, null);
                    break;
                default:
                    alglib.lsfitcreatef(x, y, c, diffstep, out state);
                    alglib.lsfitfit(state, function_cx_1_func, null, null);
                    break;
            }
            alglib.lsfitsetcond(state, epsf, epsx, maxits);
            alglib.lsfitresults(state, out info, out c, out rep);

            LogEvent("Linearization ended: " + rep.iterationscount.ToString() + " iterations; status " + info.ToString() + ";");

            return c;
        }

        public static void function_cx_1_func(double[] c, double[] x, ref double func, object obj)
        {
            ////c[0] == s   c[1] == a   c[2] == n      y = (s (1-exp(-a x^n)))
            func = c[0] * (1 - Math.Exp(-c[1] * Math.Pow(x[0], c[2]))); // s*(1-exp(-a*pow(x,n)))
        }

        public static void function_cx_1_grad(double[] c, double[] x, ref double func, double[] grad, object obj)
        {
            ////c[0] == s  c[1] == a  c[2] == n
            ////(d)/(ds) (s (1-exp(-a x^n))) = 1-e^(-a x^n)
            ////(d)/(da) (s (1-exp(-a x^n))) = s x^n e^(-a x^n)
            ////(d)/(dn) (s (1-exp(-a x^n))) = a s x^n log(x) e^(-a x^n)
            func = c[0] * (1 - Math.Exp(-c[1] * Math.Pow(x[0], c[2]))); // s*(1-exp(-a*pow(x,n)))
            grad[0] = 1 - Math.Exp(-c[1] * Math.Pow(x[0],c[2])); //1-e^(-a x^n)
            grad[1] = c[0] * Math.Pow(x[0],c[2]) * Math.Exp(-c[1] * Math.Pow(x[0],c[2])); //s x^n e^(-a x^n)
            grad[2] = c[1] * c[0] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2])); //a s x^n log(x) e^(-a x^n)
        
        }

Thank you very much. :)

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 930

One more error in grad[2] - you should call Ln, not Log10, because correct derivative formula uses natural logarithm.

My recommendations - set box constraints for your parameters with lsfitsetbc(). You may get really strange results for n<=0 (optimizer may accidentally step into such a point), so it is better to set n>=0.01 or something like that. It is also good to constrain parameter a, because very large (or very small) values may "catch" optimizer.

Does it help?

stl · **Joined:** Fri Apr 01, 2016 9:09 pm **Posts:** 5

Thank you for pointing logarithm error. I have corrected that and set upper and lower bounds as constraints. I have tried Lower bound from 0.01 to 0.2 and upper from 1.0 to 2.0 but still getting same results and info (4) everytime.

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 930

OK, I will try to investigate it tomorrow :) I will post my results here.

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 930

Hello!

Here are my findings:

1. if you want to experiment with setcond(), you should insert it BEFORE call to lsfitfit(). There is no sense in setting stopping conditions after optimization session is done :)

2. optimization algorithm stalls because one of x[] has exactly zero value. However, derivative for x^n implicitly assumes that x is non-zero, so for zero x it results in attempt to calculate Ln(0). Which, in turn, results in NaN value, which confuses optimizer.

You should either (a) replace zero x by small non-zero value, or (b) fix your derivative formula to handle borderline cases correctly.

3. I do not know why you got no success with finite differentiation version of algorithm. In my experiments it worked correctly because it avoids triggering this Log-related bug. In any case, it is a good option for fast prototyping. I recommend you to use as much as possible.

stl · **Joined:** Fri Apr 01, 2016 9:09 pm **Posts:** 5

Hi Sergey,
thank you very much for your time! When i was adding corect constructors i also moved setcond() before lsfitfit(). I have replaced zero value by 0.0001 and it works. So right now i am able to calculate function parameters by all three linearization modes. Thank you very much for your support.

forum.alglib.net

Forum rules

Nonlinear least squares fitting

Who is online