fitting models that are mostly linear

65 views
Skip to first unread message

Nicholas Zobrist

unread,
Jan 2, 2020, 4:14:00 PM1/2/20
to lmfit-py
Hi,

I have a suggestion for a new feature, and would be happy to work on a PR if it seems useful to people:

Models that are linear in their parameters are easy to fit since there is an analytic least squares solution. However, as soon as a model has at least one non-linear parameter, a generic solver is typically required, like those used in lmfit. Having to use a generic solver for these problems can be unfortunate since it can introduce convergence issues and dramatically increase the number of required function evaluations.

To avoid this problem, it is possible to fit only the non-linear parameters with a generic solver while solving directly for the linear parameters at each function evaluation. This way the generic solver has a much lower dimensional space to search for the solution.

My proposed addition is to include a fit method that implements this algorithm.

---

I learned about this algorithm in a data analysis course a long time ago and use it quite often in my own work. However, I've never seen any out-of-the-box implementations. I think lmfit would be a good place to add such a method (in my opinion) because it would benefit from the flexibility of the parameters and model objects.

Let me know if there is any interest.

-Nicholas

Matt Newville

unread,
Jan 3, 2020, 10:23:04 AM1/3/20
to lmfit-py
Hi Nicholas, 

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/c0a6bb35-1d5c-4681-ab6b-7dbf4799f78e%40googlegroups.com.


Sure: what do you have in mind?  Like, how "non-linear" could a parameter be? 

--Matt Newville 

Nicholas Zobrist

unread,
Jan 3, 2020, 1:31:08 PM1/3/20
to lmfit-py
Hi Matt,

I mean linear parameters in the sense that when all of the other parameters are fixed, the linear parameters can be solved for with scipy's linear_lsq(). For example, in y = a + b * x + c * log(x) + d * exp(e * x) the parameters a, b, c, and d are linear, but e is not. I typically use this kind of algorithm when fitting multiple peak-like models to data. In that case, all of the peak amplitudes are linear.

To give you a better idea of what the implementation might be like, I've attached an example using the double exponential problem from the documentation. Both the new and standard fitting methods converge with the starting point in the documentation, but the new method gets to the solution with many less evaluations. I've also included some commented out starting points for which the new method converges but the standard one does not.

-N


On Friday, January 3, 2020 at 7:23:04 AM UTC-8, Matt Newville wrote:
Hi Nicholas, 

On Thu, Jan 2, 2020 at 1:14 PM Nicholas Zobrist <nzob...@ucsb.edu> wrote:
Hi,

I have a suggestion for a new feature, and would be happy to work on a PR if it seems useful to people:

Models that are linear in their parameters are easy to fit since there is an analytic least squares solution. However, as soon as a model has at least one non-linear parameter, a generic solver is typically required, like those used in lmfit. Having to use a generic solver for these problems can be unfortunate since it can introduce convergence issues and dramatically increase the number of required function evaluations.

To avoid this problem, it is possible to fit only the non-linear parameters with a generic solver while solving directly for the linear parameters at each function evaluation. This way the generic solver has a much lower dimensional space to search for the solution.

My proposed addition is to include a fit method that implements this algorithm.

---

I learned about this algorithm in a data analysis course a long time ago and use it quite often in my own work. However, I've never seen any out-of-the-box implementations. I think lmfit would be a good place to add such a method (in my opinion) because it would benefit from the flexibility of the parameters and model objects.

Let me know if there is any interest.

-Nicholas

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfi...@googlegroups.com.
partial_linear_fitting.py

Kenneth Spears

unread,
Jan 4, 2020, 5:30:48 PM1/4/20
to lmfi...@googlegroups.com
Hello Nicholas and Matt,
    There is quite a literature on algorithms for solving problems with modest numbers of non-linear parameters (<20) and huge numbers of linear parameters.  The best algorithm that I know about is called the Variable Projection Algorithm.  I have written code in python 5 years ago and used a subset of the LMFIT code involving the parameter class and parameter bounds.  The essential point is that most fitting problems with LMFIT models are not sufficiently complex nor have the kind of over-determined data sets that can use the matrix methods that really speed up calculations. This algorithm is not like earlier alternative schemes of bouncing between non-linear and linear fitting steps, and uses a clever method for iterations.   I attach below an excerpt from one of my manuscripts that describes the context for time resolved spectroscopy applications (x-ray data fitting is another) with references.  These data sets may have 300x400 points with 6-10 non-linear parameters and 300-400 linear parameters.  I have an old test program that distills the essence of the algorithm in a toy example model, but it has not been updated to work with LMFIT 0.9 series.  Always on my to-do list, which gets longer every day, as the lab keeps using the version 8.3 LMFIT with Python 3 until I get to updating.
Ken Spears
Currently at University of Michigan
Emeritus Professor Chemistry, Northwestern University

From Supplement to Manuscript:

Data Analysis and Fitting

Time resolved spectroscopy data consist of a large number of wavelengths and time points that need to be modeled mathematically in a form that can allow subsequent kinetic modeling of species versus time.  The most useful  mathematical function is a sum of exponentials in time, each with a unique rate constant and each weighted by amplitudes at all wavelengths.  This type of model is often called global analysis and was discussed in an earlier section where the model outcome is named a decay associated difference spectra or DADS.  The global analysis takes advantage of many data points in both wavelength and time to reduce the effects of noise in the data set.  This occurs because there are many hundreds of wavelengths and time points to provide an overdetermined data set that is fit with a relatively few non-linear parameters in exponential mathematical functions that are weighted by amplitudes at each wavelength (linear parameters).  The fitting of such a data set is efficient, as it can be formulated as a matrix math problem with a variable projection algorithm developed by Golub and Pereyra6,7 and improved by Kaufman.8  The algorithm effectively uses a least squares method that is a mix of a very fast linear algorithm for the linear parameter fitting and much slower algorithm for the small number of non-linear parameters such as rate constants in exponentials.  The variable projection algorithm was used in an early kinetic study by Nagle9 and many other studies are discussed in a review by vanStokkum.10,11  VanStokkum and his co-workers have developed computer code called TIMP in the R computer language12 and Snellenburg in this group has provided free access to a graphical front end in JAVA called Glotaran (glotaran.org) whose acronym is derived from Global and Target Analysis (GloTarAn),13 where  experimental system modeling, such as kinetic modeling of species, is called target analysis by vanStokkum and co-workers. 

In our work we have developed our own global fitting code in the Python programming language to allow more complex fitting situations such as pump and repump experiments.  We have compared basic fitting results with Glotaran and our Python code and find equivalent parameter outputs. The code follows the mathematical matrix model suggested by Mullen and vanStokkum.12  In our case the variable projection algorithm is created from matrix and non-linear least squares fitting modules in SciPy and NumPy.  We also use the parameter class definition in a non-linear fitting package called LMFIT, which enables output parameter statistical analysis and confidence interval analysis.  LMFIT uses SciPy algorithms for non-linear fitting and adds a method of bounding parameters in the fitting process, which can be useful in complex cases. 

The general approach to an unfamiliar raw data set is to manually make a chirp corrected data set and then run a Singular Value Decomposition (SVD) to get initial guesses for the number of exponentials and their rate constants.  The data fitting with our variable projection code usually uses the slower simplex method of Nelder-Mead to ensure not being in a local minimum, and then uses the faster Levenberg-Marquardt algorithm to obtain error estimates.  When necessary there is an option to fit the raw data to obtain an optimized chirp correction simultaneous with parameter fitting.  The code models the excitation pulse as a Gaussian and the response functions are convoluted with the excitation pulse shape.  It is possible to add non-linear response functions or Raman functions at very short time scales, although this is rarely required for our data.  In addition, the code can model the pump and delayed repump experiments described in this manuscript.  There are a number of options for convenient control of complex fitting cases, where different weighting schemes or parameter bounding can allow better convergence and/or error estimates. The fitting of several hundred wavelengths and 5-6 exponentials can take under a minute on a typical laptop machine with the Nelder-Mead method, and the time is much shorter when working with repeat data sets and the Levenberg-Marquardt algorithm.  The Python environment is the Anaconda implementation with the Spyder editor.

(6)        Golub, G. H.; Pereyra, V., Differentiation of pseudo-inverses and nonlinear least-squares problems whose variables separate. Siam J. Numer. Anal. 1973, 10, 413-432.

(7)        Gene, G.; Victor, P., Separable nonlinear least squares: the variable projection method and its applications. Inverse Probl. 2003, 19, R1-R26.

(8)        Kaufman, L., A variable projection method for solving separable nonlinear least squares problems. BIT Num. Math. 1975, 15, 49-57.

(9)        Nagle, J. F.; Parodi, L. A.; Lozier, R. H., Procedure for testing kinetic-models of the photocycle of bacteriorhodopsin. Biophys. J. 1982, 38, 161-174.

(10)      van Stokkum, I. H. M.; Larsen, D. S.; van Grondelle, R., Global and target analysis of time-resolved spectra (vol 1657, pg 82, 2004). Biochim. Biophys. Acta, Bioenergetics 2004, 1658, 262-262.

(11)      van Stokkum, I. H. M.; Larsen, D. S.; van Grondelle, R., Global and target analysis of time-resolved spectra. Biochim. Biophys. Acta, Bioenergetics 2004, 1657, 82-104.

(12)      Mullen, K. M.; van Stokkum, I. H. M., TIMP: An R package for modeling multi-way spectroscopic measurements. J. Stat. Soft. 2007, 18, http://dx.doi.org/10.18637/jss.v18018.i18603.

(13)      Snellenburg, J. J.; Laptenok, S. P.; Seger, R.; Mullen, K. M.; van Stokkum, I. H. M., Glotaran: A Java-Based Graphical User Interface for the R Package TIMP. J. Stat. Soft. 2012, 49, http://dx.doi.org/10.18637/jss.v18049.i18603.


To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CA%2B7ESbpcBg3ZUa%3DQ5R%3Diz-0e--zB7jZ3kZNZ12TKqNBYMDnovw%40mail.gmail.com.

Nicholas Zobrist

unread,
Jan 4, 2020, 7:24:15 PM1/4/20
to lmfi...@googlegroups.com
Hi Ken,

Interesting, it looks like I need to do some reading before/if I try to implement anything. I’m not surprised that there is a smarter way to go about this, and it would be useful, I think, to have an algorithm with a strong basis in the literature. 

I agree that there is probably some limit where the complexity of a problem warrants looking for an alternative to lmfit, but I’m not familiar with where that boundary lies. I still think it would be beneficial if there was a way to take advantage of model linearity in some way in the package. 

As a side note, it would also be convenient to be able to solve a fully linear problem in lmfit, under the guiding principle of having a unified interface to all fitting problems. 

-N

Kenneth Spears

unread,
Jan 4, 2020, 8:03:53 PM1/4/20
to lmfi...@googlegroups.com
I believe he commercial packages in Matlab have some fitting option using the Variable Projection Algorithm, but I do not know that language and find it too expensive unless tied to a university computer system.  Therefore I really do not know how it functions and if it could be made similar in LMFIT.  My gut feeling is that when the total parameters get to be above 20-50 with nonlinear in the 10-15 range then conventional non-linear fitting becomes cumbersome and separating the linear might be useful.

To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CA%2B_JPuR%3D8hM%3DeWQuHcNTX%2B%3DF2vwrVKa8WXtw%2BhQsfAddCaMEiw%40mail.gmail.com.

Matt Newville

unread,
Jan 6, 2020, 12:59:29 AM1/6/20
to lmfit-py
Hi Nicholas, Kenneth,

Thanks very much -- this all looks interesting.  I'm not sure I fully understand Nicholas' example, but it and Variable Projection Algorithm will get some more study over the next month or so.  

FWIW, I definitely have problems where a series of lineshapes (say, Gaussian or sine waves) that each have an Amplitude parameter, which would be "linear" here.   Many times, that amplitude / line-strength is really an important part of the result.  I hesitate to say how "most people" use lmfit, but I suspect this is very common.   Also, to agree with Kenneth,  it does appear that when you get to more than 20 to 50 variables, many fitting problems need -- and can really benefit from -- extra help.  So, a method to speed up fits that have some linear variables would be very helpful.    

For some problems, I spend some effort on automating good guesses.   Perhaps that is not too different from these methods, but if they can be generalized better, that would be great.  

Just as an aside:  currently, lmfit does not have any compiled code.  We have relied on scipy to do all that for us.   I'm not necessarily opposed to supporting compiled code (I do it for other projects), but that would probably need its own discussion.



Till Stensitzki

unread,
Feb 10, 2020, 5:19:59 PM2/10/20
to lmfit-py
Variable projection is "just"  a name for doing linear regression after building the basis columns from the non-linear variables.
Basically this makes the linear parameters almost free in my problem, where the same basis columns are applied to different
measurements.

Kenneth did't I sent you some code six years ago? Some acknowledgement would have been nice.

greetings
Till

-N

Hi Nicholas, 

To unsubscribe from this group and stop receiving emails from it, send an email to lmfi...@googlegroups.com.


Sure: what do you have in mind?  Like, how "non-linear" could a parameter be? 

--Matt Newville 

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfi...@googlegroups.com.

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfi...@googlegroups.com.

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfi...@googlegroups.com.

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfi...@googlegroups.com.


--Matt 

Kenneth Spears

unread,
Feb 10, 2020, 7:40:03 PM2/10/20
to lmfi...@googlegroups.com
Hi Till,  
Yes, and I have used it with direct  acknowledgements in a simple test code based a lot on your model.  I did not send that code out in this response.  The published references are what came to mind for this response as the algorithm sources seemed most of interest. Sorry for unintentional omission.  I keep intending to update the simple example code to work with current lmfit and you are well documented there.  

 Ken

On Feb 10, 2020, at 5:19 PM, Till Stensitzki <mail...@gmx.de> wrote:


To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/1fed0c8d-71b3-4529-a4fd-678ee7e75071%40googlegroups.com.

Nicholas Zobrist

unread,
Feb 11, 2020, 1:08:06 PM2/11/20
to lmfi...@googlegroups.com
Hi Ken and Till, 

Are either of you willing to share your code to use as a starting point for integrating a variable projection algorithm into lmfit? It seems like this would be a useful addition for the community, and I have the time to work on a PR.

Best,
Nicholas

Kenneth Spears

unread,
Feb 11, 2020, 8:45:57 PM2/11/20
to lmfi...@googlegroups.com
I am traveling and perhaps between trips in a couple of weeks I can update code to work with new lmfit and send it.  Code has comments, and as I recall my implementation of the algorithm had some differences with Till, that I also documented.  The structure of method and matrix formulation was described in detail in the references I sent, and Dutch group has a lot of detail.  When back I could send those pdf files.

 Ken

On Feb 11, 2020, at 1:08 PM, Nicholas Zobrist <nzob...@physics.ucsb.edu> wrote:

Hi Ken and Till, 

Matt Newville

unread,
Feb 13, 2020, 12:04:26 AM2/13/20
to lmfit-py
Hi Ken, Nicholas, Till, 

Thanks for looking into this -- it definitely seems worth pursuing to me.  

From reading (but not trying on a real-world case) the example code posted a while back, I have a sneaking suspicion that it would not actually be much faster than just doing a simple fit, at least for problems with a handful of parameters. I guess if you get to 20 or more parameters, it might really give some savings.   Does anyone have any benchmarks on this?

Thanks,

--Matt


Nicholas Zobrist

unread,
Feb 13, 2020, 1:30:42 PM2/13/20
to lmfi...@googlegroups.com
Hi Matt,

I definitely agree that benchmarking is a good idea. It would be nice to point a user to example problems where this method is definitively better. Unfortunately, I have not yet fully understood all of the differences between the code I posted and the varpro routine Ken and Till have discussed, so I won’t comment on that.

As for my example code, it is a bit sloppy, using more function evaluations than strictly necessary. With a bit of care it is possible to re-write the algorithm so that it does not introduce any extra evaluations per fit iteration. In the double exponential example case, the number of function evaluations decreases from 609 to 106. 

I have seen similar savings in my own work where the function evaluations are the most time consuming step---even when the number of parameters is <10. However, I don’t claim that this is the case for every problem.

-Nicholas

Reply all
Reply to author
Forward
0 new messages