forum.alglib.net

ALGLIB forum
It is currently Thu Mar 28, 2024 9:06 am

All times are UTC


Forum rules


1. This forum can be used for discussion of both ALGLIB-related and general numerical analysis questions
2. This forum is English-only - postings in other languages will be removed.



Post new topic Reply to topic  [ 10 posts ] 
Author Message
 Post subject: Nonlinear least squares fitting
PostPosted: Sat Apr 02, 2016 12:26 pm 
Offline

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).


Top
 Profile  
 
 Post subject: Re: Nonlinear least squares fitting
PostPosted: Mon Apr 04, 2016 1:09 pm 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 903
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 :)


Top
 Profile  
 
 Post subject: Re: Nonlinear least squares fitting
PostPosted: Mon Apr 04, 2016 6:02 pm 
Offline

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.


Top
 Profile  
 
 Post subject: Re: Nonlinear least squares fitting
PostPosted: Tue Apr 05, 2016 8:38 am 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 903
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?


Top
 Profile  
 
 Post subject: Re: Nonlinear least squares fitting
PostPosted: Tue Apr 05, 2016 9:36 am 
Offline

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. :)


Top
 Profile  
 
 Post subject: Re: Nonlinear least squares fitting
PostPosted: Tue Apr 05, 2016 11:02 am 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 903
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?


Top
 Profile  
 
 Post subject: Re: Nonlinear least squares fitting
PostPosted: Tue Apr 05, 2016 12:19 pm 
Offline

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.


Top
 Profile  
 
 Post subject: Re: Nonlinear least squares fitting
PostPosted: Tue Apr 05, 2016 9:07 pm 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 903
OK, I will try to investigate it tomorrow :) I will post my results here.


Top
 Profile  
 
 Post subject: Re: Nonlinear least squares fitting
PostPosted: Wed Apr 06, 2016 9:11 am 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 903
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.


Top
 Profile  
 
 Post subject: Re: Nonlinear least squares fitting
PostPosted: Wed Apr 06, 2016 7:50 pm 
Offline

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.


Top
 Profile  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 10 posts ] 

All times are UTC


Who is online

Users browsing this forum: No registered users and 60 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to:  
cron
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group