GNOME Bugzilla – Bug 690045
Implement least square fitting
Last modified: 2018-05-22 13:55:46 UTC
When applying least squares technics, one tries to minimize a least squares cost function and gnumeric solver is able to do so. However, usually, one needs not only the optimal values of the parameters, but also for the confidence intervals, one needs the covariance matrix, which is the inverse of the Hessian matrix of the cost function. This Hessian matrix is internally computed by gnumeric's solver when performing Newton iterations and a simple option to output it somewhere in the workbook would be very valuable. I wrote a little patch (below) to output the Hessian on the command line at every Newton iteration. It is dirty but useful. root@LGGE-PC144:/usr/local/gnumeric# diff plugins/nlsolve/gnm-nlsolve.c.orig plugins/nlsolve/gnm-nlsolve.c 318a319,321 > /* Frédéric Parrenin - for Hessian outputing */ > int ifp,jfp; > /* Frédéric Parrenin - end of for Hessian outputing */ 321a325,332 > /* Frédéric Parrenin - Hessian outputing */ > printf("\n"); > for (ifp=0;ifp<n;ifp++) > { > for (jfp=0;jfp<n;jfp++) {printf("%g ",H[ifp][jfp]);} > printf("\n"); > } > /* Frédéric Parrenin - end of Hessian outputing */
https://docs.google.com/open?id=0BzX8dPORePBsNWdsSC02ZGF3UW8 In this file, I minimized 'Inversion'D2 by changing 'Inversion'D5:D9 and I put the Hessian in the 'P^-1' sheet. What is strange here is that the bottom-right term of the Hessian, which is a diagonal term, is negative, while it should be positive because we are close to a minimum (the function to minimize should be convex near the minimum). Does anybody has a clue why this happens? Any help would be greatly appreciated.
This does not look right: printf("%g ",H[ifp][jfp]); H is not necessarily a double, so using the %g format is a bad idea. GNM_FORMAT_g would be much more appropriate.
The Hessian is computed numerically. It is not surprising that you would see strange numbers now and then.
Is it not a problem that the Hessian is not positive definite?
> Is it not a problem that the Hessian is not positive definite? Well, obviously we would rather have the numbers match the theory, but it's a fact of life that sometimes we don't get to have that. There are many causes: * Using a secant instead of a proper derivative * Improper objective function (non-continuous, for example) * Rounding error * Underflow * Outright bugs, perhaps The algorithms used should be fairly robust in the face of numbers that don't make sense.
This type of algorithm often (usually?) only use an approximate Hessian since it is not feasible to calculate the true Hessian. That's why looking at the matrix used can be misleading.
I would again advocate for implementing a least-squares fitting method, more specific than the solver. Optimizing a least-squares cost function only requires to estimate the gradient of the model to fit, which should be far more precise than numerically calculating the Hessian of the cost function. Moreover, least-squares problems are the most often used optimization problems.
least square fitting (ie. fitting model parameters to the data to minimizre the sum of the squares of the deviation) has absolutely nothing to do with "solving". These are two completely distinct ideas. In fact noting stops you of performing least square fits.
Least square fitting certainly has to do with solving. In least square fitting, we minimize a cost function which has a particular form. Moreover, the least square theory allows to evaluate the confidence interval on the inferred model parameters (and also on the model output), which the current gnumeric solver cannot do. So what stops me at performing least square is a convenient and efficient way to do it in Gnumeric, including the reconstruction of confidence intervals. I will try to write my own specific code but I believe a generic solution inside gnumeric would be a great added value to the software.
I would like to change the title of this bug to "Implement least square fitting" but I don't see how to do it (and maybe I don't have the permissions to do so). I can ask some money to pay somebody to implement least square fitting. (good chance that it would be accepted but not 100% sure) Would some of gnumeric's dev be interested in this kind of contract? If so, how much would you ask for that? If no gnumeric dev is interested, could you suggest somebody else?
Without some reference to a specific algorithm you would like to see implemented, a generic "Implement least square fitting" seems to be pretty much an open ended suggestion.
Dear Andreas, Please see this reference: https://docs.google.com/open?id=0BzX8dPORePBsTFhsVnA4UkxWSEk The least square problem is described section 3 starting p. 57. The classical quasi-Newton algorithm is described p. 69. Please tell me how much you would ask to implement this feature. Best regards, Frédéric Parrenin
Hi Frédéric, I am not in a position at this time to commit to any specific implementation task. Andreas
Frédéric, I had a brief look at your reference and it seems similar to the materials I know on this issue: They are not appropriate to be passed to a programmer to be implemented. Andreas
Dear Andreas, It is indeed a classical material. Why do you think it is not an appropriate reference to be passed to a programmer? Of course just the general ideas are described. I can find more specific references to the implementation of the quasi-Newton algorithm if needed. For example, there are references given by the m1qn3 devs: https://who.rocq.inria.fr/Jean-Charles.Gilbert/modulopt/optimization-routines/m1qn3/m1qn3.html Best regards, Frédéric Parrenin
We definitely seem to be speaking different languages: m1qn3 uses L-BFGS, L-BFGS is a solver, ie. it finds the minimum of sum function. This is different from least-square approximation. Of course one can formulate a least square approximation to a given model as a minimization problem and then solve it using a solver, but these are still 2 distinct steps. I don't see anything really wrong with the current solver. Andreas PS: Why do I think it is not an appropriate reference to be passed to a programmer? Programmers are typically not mathematicians. The text does not even try to indicate how to translate the mathematics into an algorithm. But of course you may know programmers who could handle that...
Dear Andreas, Concerning comment #2, what would be the syntax? Frédéric Parrenin
Concerning comment #16, it is my bad, I read too quickly. M1qn3 is indeed a solver which needs to evaluate the function and its gradient.
So to be more precise of what I would like to see in this quasi-Newton least-square fitting algorithm: - the algorithm would define the model as some input parameters/cells to be modified (same than the current solver) + some output parameters/cells (the current solver only take one cell to be minimized or maximized) - the algorithm would check if the input parameters are not the results of calculations - the algorithm would define the \mu parameter defined in eq. (3.51) of Tarantola-2005 (constant in this case). - the algorithm would take as input the covariance matrices on input parameters (C_M) and output parameters (C_M) and check that their dimensions are equal to the numbers of input and output parameters respectively. For simplicity, the algorithm could offer an option where C_M (resp. C_D) are diagonal and therefore take as input only the vector of the diagonal terms. - the algorithm would also take as input the places in the workbook where to ouput the posterior covariance matrices {\tilde C}_M and {\tilde C}_D. - the algorithm would take as input the number of iterations to be performed. - the aglorithm would then apply the iterative walk defined in eq. (3.51) of Tarantola-2005, calculating the gradient of the model at each iteration and displaying the result of the cost function defined, for example, in eq. (3.35) of Tarantola-2005. The operator could interrupt the process at every time (like in the current solver). - when the max number of iterations is reached or when the operator has interrupted the process, the algorithm would populate the input cells with the values reached in the last iterations and would output the posterior covariance matrices where it has been defined. Please tell me if there is still something unclear in this RFE.
Well, "the algorithm would define the \mu parameter defined in eq. (3.51) of Tarantola-2005 (constant in this case)." I don't know what it means that an algorithm 'defines' a parameter. Moreover eq 3.51 does not define the parameter but uses it. More importantly: "the aglorithm would then apply the iterative walk defined in eq. (3.51) of Tarantola-2005, calculating the gradient of the model at each iteration". How should it "calculate" the gradient. At best such an algorithm could try to estimate the gradient by calculating difference quotient. Similarly {\tilde C}_M and {\tilde C}_D can't be calculated in this setup but as best estimated. Tarantola assumes that you have a model that is given in a form that it can be analytically analysed. In your case you have a spreadsheet formulas which do not lend themselves to that treatment. good luck of finding somebody interested in trying to implement something like this.
-- GitLab Migration Automatic Message -- This bug has been migrated to GNOME's GitLab instance and has been closed from further activity. You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.gnome.org/GNOME/gnumeric/issues/206.