Optimization in MMA for time resolved spectroscopy

Question

I want to fit some spectroscopic data with a model. I detect the fluorescence of something over different wavelengths. The detection over each wavelength is one data set in time. The lifetimes of fluorescence do not change over different wavelength, but the amplitude of detected signals changes over different wavelengths.

1. We assume that the instrument response function(IRF) is Gaussian and fit the data with some exponential functions. The number of these exponential functions are not known prior to fitting. However, someone might guess on the basis of previous knowledge of the system under study.
2. The detected signal is a convolution of IRF and the model, in this case an exponential decay.
3. There are several detected signals which must be fitted with the same model, but each model for each detected signal have some common parameters and some parameters which are not shared. The common ones are lifetimes of the decay of detected signals and those which are not shared are the amplitudes of each signal.
4. The goal is to subtract each experimental data point from the corresponding point generated by convolution of IRF and model, square it, add all of the squared terms from all the data points, and finally minimize the total term.

$\Delta$ is full width at half maximum and $\mu$ is the position of IRF.

$a_i$ are amplitudes and $\tau_i$ are corresponding lifetimes. Here, I used 3 exponential functions hence 3 lifetimes.

The detected signal is the convolution of IRF and model:

After integration I should fit the resulting function with my experimental data.

At this point I consider one exponential function and one data set. Here are MMA codes:

SetDirectory[NotebookDirectory[]]
data = Import["data.dat"];(* Experimentally measured data*)
taxis = Import["taxis.dat"];(* the time axis*)
numberOFexp = 1;
model[z_] := 
 Sum[ToExpression["a" <> ToString[i] <> ToString[j]]* 
   Exp[-z/ToExpression["\[Tau]" <> ToString[i]]], {i, 1, numberOFexp}]

irf[z_] := (2*Sqrt[2*Log[2]])/(\[CapitalDelta]*Sqrt[2*Pi])*
   Exp[-4*Log[2]*((
      z - ToExpression["\[Mu]" <> ToString[j]])/\[CapitalDelta])^2];

Signal[t_] := 
 Evaluate@ParallelTable[
   Integrate[irf[z]*model[t - z], {z, 0, t}] + 
    ToExpression["b" <> ToString[j]], {j, 1, 1}]

Signal\[TripleDot]points = ParallelMap[Signal, taxis];
diff = Dot[Flatten[Signal\[TripleDot]points] - data[[All, 2]], 
   Flatten[Signal\[TripleDot]points] - data[[All, 2]]];

AbsoluteTiming[
 FindMinimum[
  diff, {{b1, 1.0}, {a11, 300.0}, {\[Tau]1, 10.0}, {\[Mu]1, 
    40.0}, {\[CapitalDelta], 20.0}}, Method -> "LevenbergMarquardt"]]

"LevenbergMarquardt" gives the fastest solution, it give a solution after one second. Other methods take from 3 to 120 seconds. The real problem is much bigger. I must fit at least 104 parameters. I tried to use MMA and used FindMinimum. There are two problems:

• It consumes a huge amount of RAM, 24GB in this case and because I don't have more RAM it crashes. I did not use any particular method for the real problem.
• It takes a lot of time, if it does not crash.

Then, I used NMinimize. But, it takes more than 5 hours to get the result. Sometimes it won't give any result in that time windows.

I fitted the same data, which took 1 seconds in MMA by using FindMinimum and "LevenbergMarquardt", in Matlab and it takes 0.1 seconds. MMA gives a smaller minimum, However, the lifetimes obtained by Matlab is as good as MMA and physically acceptable. So, Matlabs gives me the solution 10 times faster. However, for the real problem I can't use Matlab because I don't know how to implement constrained minimization in Matlab and in general the syntax is not as good as MMA. So, first I want to make my MMA code as fast as that of Matlab and if I failed to do so, learn Matlab and try to solve my problem in Matlab. I should mention that I think fminsearch in Matlab uses 'Nelder-Mead simplex direct search' algorithm.

Here is Matlab code:

the function which must be saved in an .m fi

function myfun = fun(p)
T1 = p(1);
del = p(2);
a11 = p(3);
b1 = p(4);
mu1 = p(5);
t = [];% This is time axis. Get it from below
y1 = [];%This is experimentally measured data. Get it from below

myfun1 = sum((y1 - (b1 + (1/2)*a11*exp ((-16*t*T1 + 16*mu1*T1 + del^2/log (2))/(16*T1^2)).*(erf ((del^2 + 8*mu1*T1*log (2))/(4*del*T1*sqrt (log (2)))) - erf ((del^2 + 8*(-t + mu1)*T1*log (2))/(4*del*T1*sqrt (log (2))))))).^2);
myfun = myfun1;

The other part which must be saved in the same place as the previous function:

format compact
format long
%starting guess
pguess = [10,20,300,1,40];
%options = optimset('PlotFcns',@optimplotfval);
tic
[p,fminres] = fminsearch(@fun,pguess)
toc

data, the signal, for MMA

taxis for MMA which is the same as data[[All,1]] or t for Matlab

y1, the signal, for Matlab which is the same as data[[All,2]]

Is there any way to get a faster solution, possibly as fast as Matlab, in MMA? I don't want to compile the objective function because it takes a lot of time itself.

Answer 1

So I'll show you the utility of singular value decompositions for one-component systems, and show how to make a 1-component fit. Doing this for multiple-component systems is likely to be harder.

First, import the data:

data = N@Import[SystemDialogInput["FileOpen"]][[1 ;; 300]];

Plot the data:

<< Developer`
c = ToPackedArray[{0.3, 1.0, 0.1}];
Image[(data/10)\[TensorProduct]c]

Compute the SVD and show the first component:

{U, S, V} = SingularValueDecomposition[data];
Image[S[[1, 
    1]] U[[;; , 1]]\[TensorProduct]V[[;; , 1]]\[TensorProduct]c/10]

Here's the second component (brightened since it's weak):

Image[3 S[[2, 
    2]] U[[;; , 2]]\[TensorProduct]V[[;; , 2]]\[TensorProduct]c]

You'll notice that the second component looks like noise, which is a strong indication that the data is the fluorescence of a single chemical species. Note that nowhere did I tell the program that the data was a single-component fluorescence trace; that's the magic of the SVD.

We can now isolate the time profile and fit it. Here's what it looks like:

timeProfile = -U[[;; , 1]];
ListLinePlot[timeProfile, PlotRange -> All]

I did a rough model as follows:

fun[\[Tau]_, a_] := 
  a Re@InverseFourier[
      Fourier[Join[ConstantArray[0, 200], 
         Exp[-Range[1400]/\[Tau]]]] Fourier[
        Chop@Exp[-0.00005 Range[-200, 1399]^2]]][[200 ;; 1024 + 199]];

I was having trouble minimizing Norm[fun[\[Tau], a] - timeProfile] using NMinimize (getting error messages, possibly due to bad starting point, but that's another problem and I'm not familiar with minimization, so perhaps ask someone else about that.