Fitting a Multi-Exponential Decay With Constraints

Question

I am having trouble finding a way to fit my data in Mathematica to a multi-exponential decay, with the critical condition that all contributions to the decay sum to 1.

To start, I can fit my data to single exponential decay just fine using:

model = Exp[-t/k];
fit = FindFit[cordata,model,k,t];

This gets me a decent fit, and I can also fit my data to a multi-exponential decay using:

modelx = m1 * Exp[-t/k1] + m2 * Exp[-t/k2];
fit = FindFit[cordata,modelx,{{m1,.5},{k1,x},{m2,.5},{k2,x}},t];

where x is the value of k returned by the single exponential fitting. This multi-exponential fit is even better (so long as I provide reasonable starting values for m1, m2, k1 and k2. However, if I attempt to constrain the fitting so that all contributions sum to 1:

modelx = m1 * Exp[-t/k1] + m2 * Exp[-t/k2];
fit = FindFit[cordata,{modelx,m1+m2==1},{{m1,.5},{k1,x},{m2,.5},{k2,x}},t];

The fitting no longer works. m1, m2, k1 and k2 are just returned as the starting values and thus the fit is no better than for a single exponential decay.

Just for the record, my data basically is an exponential decay that starts at 1 and decays to zero. Fitting it to a multi-exponential decay where all components sum to 1 should be totally reasonable (I think). I need to lock the sum of the components because it is very important to the interpretation and application of the fit. If anyone has any advice on how I can do this, I would greatly appreciate it. Sorry if the answer is obvious. I am new to Mathematica, but I have tried a variety of options to get around this with no success.

Also, if it is somehow possible to take a normal multi-exponential decay and adjust it so that all contributions sum to 1, I would be open to that as well. Sorry if I missed this, my math is not the best.

Edit: Although I showed a second order exponential decay to be concise, I actually want to apply this method to a sixth order decay:

modelx = m1 * Exp[-t/k1] + m2 * Exp[-t/k2] + m3 * Exp[-t/k3] + m4 * Exp[-t/k4] + m5 * Exp[-t/k5] + m6 * Exp[-t/k6]

Where all factors sum to 1. So I cannot define one factor using the other.

Answer 1

There is much to exponential fitting that is obscured by non-linear fits to a number of discrete components. Unless you have thousands of densely sampled points (the vector t in your notation) you will never truly resolve 6 components. If 6 spikes appear, it is merely a convenient fiction. Beware. A true model containing a wide block of continuous signal has a least-squares fit composed of a few delta-functions. Nothing like the input continuous model. Remember, exponentials are hopelessly non-orthogonal functions.

It is vitally important that you abandon a non-linear approach to finding multiple components. It is far better to use the Mathematica implementation of non-negative least-squares by Lawson and Hanson. See my answer and the comments to this question. Provide a list of k values (in your notation) and let a function calculate amplitudes at these many locations. A non-linear simultaneous calculation of both k values and amplitudes will never converge in practice.

First a function for log sampling of the relaxation time axis, T2 in magnetic resonance, otherwise k in your notation.

T2sample[n_Integer, min_, max_, ilog_Integer] :=
   If[ilog == 0,
      min + Range[0., n - 1]*(max - min)/(n - 1.), 
      min * 10.^(Range[0, n - 1]*Log[10, max/min]/(n - 1.))]

The following t2l1con1 uses linear programming to return a list of discrete components at the input list of relaxation times T2, or k values in your notation. There is a constraint to unit total amplitude.

• t: list of measurement times
• dat: list of data amplitudes at the measurement times
• stn: list of standard deviations of the data dat
• t2: list of relaxation times at which amplitudes are to be calculated

The code.

t2l1con1[t_List, dat_List, stn_List, t2_List, opts___] :=
   Block[{nt = Length[t], nt2 = Length[t2], c, m, rhs, soln},
      (* objective function *)
      c = Flatten[{Table[0., {nt2}], Table[1., {2*nt}]}];
      (* data constraints *)
      m = Table[0., {nt + 1}, {nt2 + 2*nt}];
      m[[Range[nt],Range[nt2]]] = DiagonalMatrix[1/stn].E^Outer[Times,-t,1/t2];
      Do[m[[i, {nt2 + 2*i - 1, nt2 + 2*i}]] = {1., -1.}, {i, 1, nt}];
      (* weighting requiring unit sum of amplitudes *)
      m[[nt + 1, Range[nt2]]] = 1.;
      (* right hand side *)
      rhs = Join[Transpose[{dat/stn, Table[0., {nt}]}], {{1., 0.}}];
      (* find solution *)
      soln = Chop[LinearProgramming[c, m, rhs, opts]];
      Transpose[{t2, soln[[Range[nt2]]]}]
 ]

A bi-exponential model with no noise. Note the split peak due to the model relaxation time (0.100) lying between two input T2 values. Simply add the two amplitudes. Split peaks may be merged by re-running with an additional T2 value inserted into the input list. Running a linear method many times will still be faster than a non-linear attempt.

Block[{t, dat, stn, t2},
   t = Range[0.01, 0.5, 0.005];
   dat = 0.3*E^(-t/0.05) + 0.7*E^(-t/0.100);
   stn = ConstantArray[0.01, Length[dat]];
   t2 = T2sample[250, 0.010, 1.00, 1];
   ListLogLinearPlot[
      t2l1con1[t, dat, stn, t2],
      Filling -> Axis, Frame -> True, 
      FrameLabel -> {"T2 Relaxation Time", "Amplitude"}]
]

The same bi-exponential model with about 5% noise. Note the variation from the true, input model. Such is life with exponential fitting.

Block[{t, dat, stn, t2},
   t = Range[0.01, 0.5, 0.005];
   dat = 0.3*E^(-t/0.05) + 0.7*E^(-t/0.100) +
               0.05*RandomReal[{-1, 1}, Length[t]];
   stn = ConstantArray[0.01, Length[dat]];
   t2 = T2sample[250, 0.010, 1.00, 1];
   ListLogLinearPlot[
      t2l1con1[t, dat, stn, t2],
      Filling -> Axis, Frame -> True, 
      FrameLabel -> {"T2 Relaxation Time", "Amplitude"}]
]