Levenberg-Marquadt Algorithm (GSL implementation)

  • Thread starter pisuke
  • Start date
  • Tags
    Algorithm
In summary: I mean for example you can do an exponential fitting for: y = a*exp(b*t) + c, and you will get a good result, but if you do an exponential fitting for y = a*exp(b*t) + c*exp(d*t) you will get a bad result.I consider a good approximation for the function y = a*exp(-b*t) + c*exp(-d*t) a solution for which:1. The sum of the squares of the deviations between the curve and the data is a minimum2. The parameters are physically sound, that means that the decay times are positive, that the maximum value is coherent
  • #1
pisuke
3
0
Hello there,

I have to fit a one term exponential curve: y = a*exp(-b*t), and I have to implement the fitting procedure in a C program. Some days ago I discovered the GSL (GNU Scientific Library), and I told myself better later than never ;-)

When reading the GSL documentation and the implementation examples at http://www.gnu.org/software/gsl/man...-for-Nonlinear-Least_002dSquares-Fitting.html I gladly found that one of the examples is an exponential curve of the type y = a* exp(-b*t) + c, it means that with just a slight modification of the code I will get what I am looking for.

However, in the code there's something that I don't understand, the fitting procedure adds a term called sigma with a value of 0.1 in every calculation of the function, here is the explanation on the documentation:
The following example program fits a weighted exponential model with background to experimental data, Y = A \exp(-\lambda t) + b. The first part of the program sets up the functions expb_f and expb_df to calculate the model and its Jacobian. The appropriate fitting function is given by,

f_i = ((A \exp(-\lambda t_i) + b) - y_i)/\sigma_i

Here is the code doing the stuff:
Code:
...
sigma[i] = 0.1;
...

for (i = 0; i < n; i++)
         {
           /* Model Yi = A * exp(-lambda * i) + b */
           double t = i;
           double Yi = A * exp (-lambda * t) + b;
           gsl_vector_set (f, i, (Yi - y[i])/sigma[i]);
         }
...

I don't know the reason and what is more important I don't know if I have to use such term in my own program.

I am not a mathematician, I come from the biologic field and I need to fit data which comes from light absorption measurements.
Hopefully can someone here give me a hand on that issue.

Thanks,
Pisuke
 
Mathematics news on Phys.org
  • #2
While waiting for somebody more experienced for a proper answer, here's a snack:

You probably don't need the L-M algorithm at all; your equation can be fitted using logarithmic regression.
That is, given y = a* exp(-b*t), applying a logarithm to both sides, you get: ln y = ln a - b * t, which is a simple linear equation, which you can fit by the least squares method, for example.
Just use, instead of your (t,y) data pairs, the pairs (x', y') = (t, ln y), fit the linear equation y' = p + q * x', and then use b = -q and a = exp(p) in your original equation.

The "small" difference of a similar equation with an added constant +c in the L-M case example you cite, is due precisely because, after adding a constant, log regression is no longer applicable.
 
Last edited:
  • #3
Hello Dodo,

Thanks for your answer, but applying the natural logarithm to both sides of the equation does not solve my problem, what I have really is a double exponential decay of the form:

y = a*exp(-b*t) + c*exp(-d*t)

I approximate the double exponential to one exponential and I accept the deviation I get, because what I need is the factor, which multiplies the t in a single exponential decay.

y = a*exp(-b*t)

I have already tried your ln approach and I have compared the results with a direct LM fitting of the "deviated" one term exponential curve (I have used MATLAB fitting tool), the LM result is 10% better than the regression line result, that's why I need to fit directly the exponential decay.

Actually what would be ideal, would be a fitting procedure using the double exponential, and then somehow converge the parameters b and d in one, but that as long as I know is not possible.

Pisuke
 
  • #4
As stated in the text that accompanies the example

The main part of the program sets up a Levenberg-Marquardt solver and some simulated random data. The data uses the known parameters (1.0,5.0,0.1) combined with gaussian noise (standard deviation = 0.1) over a range of 40 timesteps

So I think the extra sigma term in the example is the noise with a standard deviation of 0.1.

Hope it helps!
 
  • #5
hy,
i've read your question about the two exponential fitting of the function:

y = a*exp(-b*t) + c*exp(-d*t)


i've the same problem but i must implement the solving algorithm in excel

can you help me finding code example or other material?

you've write:

"I approximate the double exponential to one exponential and I accept the deviation I get, because what I need is the factor, which multiplies the t in a single exponential decay.

y = a*exp(-b*t)

I have already tried your ln approach and I have compared the results with a direct LM fitting of the "deviated" one term exponential curve (I have used MATLAB fitting tool), the LM result is 10% better than the regression line result, that's why I need to fit directly the exponential decay"

is it really a good approximation?

my email adress is :
Lucabellesi@yahoo.it

best regards
 
  • #6
Hi!

lucbel said:
hy,
i've read your question about the two exponential fitting of the function:

y = a*exp(-b*t) + c*exp(-d*t)


i've the same problem but i must implement the solving algorithm in excel

can you help me finding code example or other material?


In your case I would implement the algorithm using visual basic for applications. you can find the source code of a C implementation on the GSL (GNU Scientific Library): http://linux.duke.edu/~mstenner/free-docs/gsl-ref-1.1/gsl-ref_36.html [Broken]

Unfortunately I have not experience in visual basic for applications, so I can't give you hints on translating the code to basic.


lucbel said:
is it really a good approximation?


Buff... that depends on the function parameters, in my case it could be better but it could also be worse. However you should realize that excel, even offering a exponential fitting does not do an actual non linear fitting procedure. Excel does a linearization applying the natural logarithm to your data and afterwards a linear fitting procedure. And the results are radically different if your function is a complex exponential and not a single one.

I recommend you to play a little bit with Matlab's curve fitting tool (it has a fantastic help). Matlab also offers an Excel link tool, which allow you to link both programs, you can generate a function in matlab, which is a high level language very very flexible and then call that function from Excel. The MATLAB function can get as input the values which are in some cells of an excel sheet, and write the output in another cell of the same excel sheet. The shortcoming is that you need both programs, the advantage is that is easier than using visual basic, specially if you don't have previous experience.
Again the MATLAB help for the excel link tool is outstanding.

Regards,
Pisuke
 
Last edited by a moderator:
  • #7
The 2nd attachment (Nonlinear Least Squares Example) in the below link might help.

https://www.physicsforums.com/showthread.php?t=97391 [Broken]

There are two ways this can be solved with Excel. The easiest is to set up the equations within the spreadsheet and solve using built-in Excel matrix functions. The second is to write VBA routine, which would still use the built-in Excel matrix functions.
 
Last edited by a moderator:
  • #8
Attached is an example of how one might solve such a problem in Excel. The data I just made up. In this example, the matrix calculations are done with Excel built-in spreadsheet functions (i.e. MMULT, Transpose, Minverse) and the iterative portion is handled by some simple VBA code that does a copy/paste operation on the solution vector and adjusts the L-M parameter (Lambda) based on how the sum of sqaures is behaving. Of course, one can get much fancier by using matrix decomposition techniques such as QR instead of calculating a matrix inverse and using a line search to determine updated parameters instead of a fixed multiplication factor (mu), but additional VBA code would be needed. This particular example is tricky because it is ill conditioned and actually is unstable.
 

Attachments

  • NONLIN.gif
    NONLIN.gif
    55.2 KB · Views: 671

1. What is the Levenberg-Marquadt algorithm?

The Levenberg-Marquadt algorithm is an optimization method used to find the minimum of a function by iteratively adjusting the parameters of the function. It is commonly used in nonlinear least squares fitting problems.

2. What is the GSL implementation of the Levenberg-Marquadt algorithm?

The GSL (GNU Scientific Library) implementation of the Levenberg-Marquadt algorithm is a software library that provides a robust and efficient implementation of the algorithm. It is written in the C programming language and can be used for a variety of optimization problems.

3. How does the Levenberg-Marquadt algorithm work?

The algorithm works by minimizing a cost function, which is a measure of the difference between the expected and actual values of a set of data. It does this by iteratively adjusting the parameters of the function being optimized until the cost function is minimized.

4. What are the advantages of using the Levenberg-Marquadt algorithm?

The Levenberg-Marquadt algorithm is often preferred over other optimization methods because it is robust, efficient, and can handle a wide range of problems. It also has good convergence properties, meaning it can find a solution even when the initial guesses for the parameters are far from the true values.

5. Are there any limitations to the Levenberg-Marquadt algorithm?

Although the Levenberg-Marquadt algorithm is a powerful and widely used optimization method, it does have some limitations. It may not always find the global minimum of a function and can be sensitive to initial parameter values. It also requires the user to have a good understanding of the problem and choose appropriate starting values for the parameters.

Similar threads

  • General Math
Replies
2
Views
666
  • MATLAB, Maple, Mathematica, LaTeX
Replies
9
Views
4K
  • Programming and Computer Science
Replies
1
Views
2K
  • Set Theory, Logic, Probability, Statistics
Replies
1
Views
1K
Replies
10
Views
5K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
2K
  • Programming and Computer Science
Replies
2
Views
840
  • Mechanical Engineering
Replies
2
Views
926
  • Programming and Computer Science
Replies
1
Views
1K
  • General Math
Replies
1
Views
1K
Back
Top