# Optimization in MMA for time resolved spectroscopy

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.

• Random comic sans code – shrx Oct 29 '14 at 21:46
• One question: you say "The common ones are lifetimes of the decay of detected signals and those which are not shared are the amplitudes", which implies that $\tau_1=\tau_2=\tau_3:=\tau_\text{eff}$. However, that seems sort of nonsensical from a fitting perspective, since your model can then be simplified to $m(t)=(a_1+a_2+a_3)\exp(-t/\tau_\text{eff})=a_\text{eff}\exp(-t/\tau_\text{eff})$, which means that your model has two redundant degrees of freedom. Am I missing something? – DumpsterDoofus Oct 29 '14 at 21:49
• @DumpsterDoofus. I did not mean that tau1 = tau2. I mean if I calculate the convolution for each data set, the tau1 is the same in all those calculated convolution, and likewise for tau2 and tau3. But, a1 is different for each data set and hence when I calculate convolution for each data set a1 is different for every calculation of convolution. From physical point of view, if I detect signal over different wavelengths, the lifetimes do not change but the amplitudes in different signals detected at different wavelengths can be different. – MOON Oct 29 '14 at 21:54
• In other words, because the 2D spectrum has the form $$\psi(\lambda,t)=\sum_{j=1}^{N_\text{comp}}c_j(t)\epsilon_j(\lambda),$$ the 2D dataset will admit a rank-1 decomposition $$D=\sum_{j=1}^N\sigma_j\mathbf{u}_j\otimes\mathbf{v}_j$$ which is precisely the singular value decomposition of the dataset. By selecting the $N_\text{comp}$ largest principle components, you can directly determine the $c_j(t)$ and $\epsilon_j(\lambda)$ and bypass the fitting process altogether. Then you can just deconvolve the IRF from the $c_j$, and do $N_\text{comp}$ 1-dimensional simple-exponential fits. – DumpsterDoofus Oct 29 '14 at 23:25
• Ugh, never mind about what I said about directly decomposing the data via SVD, I was wrong; SVD has orthogonal bases, and the decaying exponentials are not orthogonal to each other, so SVD gives something different. The article does say that SVD can give estimates for the number of components, though. – DumpsterDoofus Oct 31 '14 at 16:49

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.

• Thank you. I must spend time to learn your code. To me it seems that ListLinePlot[timeProfile] is the spectra of the sample not its time profile. However, I am nit sure currently. – MOON Nov 2 '14 at 0:14
• Time profile information is in U matrix and spectra information is in V` matrix. When you use timeProfile = Table[-100*U[[All, i]], {i, 1, 3}]; ListPlot[timeProfile, PlotRange -> All] You see that there is just one decay curve and others are noise so it suggests for one lifetime to describe the data. It is simply beautiful! – MOON Nov 2 '14 at 0:42
• @yashar: Oh sorry, is the vertical axis the time axis and the horizontal axis the frequency axis? I must have them backwards, apologies. I fixed my answer to reflect this. – DumpsterDoofus Nov 2 '14 at 0:46