# How to speed up Voigt profile fits?

Voigt profile fits in Mathematica seem terribly slow. For an example data set with 81 points a corresponding fitting procedures for Voigt fits is >1000 times slower than for Gaussian or Lorentzian profile. How can Voigt profile fits be speed up?

Here is what I do. First we define the Voigt and Gaussian profile. As Mathematica often complains about δ and σ being <0, I use Abs instead of constraining the model as this proved faster. In addition I use the compile command to make the function even faster.

voigtprofile = Compile[{{δ, _Real, 0}, {σ, _Real, 0}, {A, _Real,0}, {ν0,_Real, 0}, {ν, _Real, 0}}, A PDF[VoigtDistribution[Abs@δ, Abs@σ], ν - ν0]];
gaussian[σ_, A_, ν0_, ν_] := Return[A PDF[NormalDistribution[ν0, σ], ν]];


Then I create noisy example data sets with a bit of noise:

noisyDataV = {#,
voigtprofile[0.15, 0.1, 1, 0, #] + RandomReal[{-0.1, 0.1}]} & /@ Range[-2, 2, 0.05];
noisyDataG = {#, gaussian[0.1, 1, 0, #] + RandomReal[{-0.1, 0.1}]} & /@
Range[-2, 2, 0.05];


And we use NonlinearModelFit to fit the data with excellent start parameters:

tv = AbsoluteTiming[vfit = NonlinearModelFit[noisyDataV,
voigtprofile[δ, σ,
A, ν0, ν], {{δ, 0.15}, {σ, 0.1}, {A,
1}, {ν0, 0}}, ν];]
tg = AbsoluteTiming[gfit = NonlinearModelFit[noisyDataG,
gaussian[σ,
A, ν0, ν], {{σ, 0.1}, {A, 1}, {ν0,
0}}, ν];]


And if we compare the required time for fit:

tv[]/tg[]


I get values between 1000 and 6000, which is terrible. In addition, selecting a fit Method e.g. NMinimize or other does at best yield the same result.

As this minimal example is just a very simple example, times scale up to unbearable long times for more realistic scenarios with real data. I'm glad for any hint on how to speed this simple example up.

Thanks to J. M. for the suggestion. Instead of using the computationally intensive VoigtDistribution, we can use the Pseudo-Voigt numeric approximation from this post which exhibits very low errors compared to the analytical Voigt (see here or references in here) and speeds up the fitting procedure by approximately a factor of 100 or more.

pseudoVoigtProfile = With[{n = 24, τ = 12}, With[{d = N[Range[n] π/τ],b = N[Exp[-(Range[n] π/τ)^2]], s = N[PadRight[{}, n, {-1, 1}]], sq = N[Sqrt],sp = N[Sqrt[2 π]]},
Compile[{{δ, _Real}, {σ, _Real}, {x, _Real}},
Module[{z = (x + I δ)/(σ sq), e},
e = Exp[I τ z];
Re[(I (1 -e)/(τ z) + (2 I z/τ)b.((e s -1)/((d + z) (d - z))))]/(σ sp)],
RuntimeAttributes -> {Listable}]]];


and fit it to the noisy test data

tpv = AbsoluteTiming[pvfit = NonlinearModelFit[noisyDataV,
pseudoVoigtProfile[δ, σ, ν - ν0], {{δ, 0.15}, {σ, 0.1}, {A, 1}, {ν0, 0}}, ν];]


which leads to a ratio tpv/tg ~ 10 on my computer.

Note: Although the Pseudo-Voigt profile is a numerical approximation it will do the job for most daily life applications where the noise on the e.g. optical spectrum is at least an order of magnitude larger than the numeric error of the approximation.

• As Voigt is Gaussian convolved by Lorentzian another approach regularly used in spectrometry is to work in Fourier space: see sites.math.washington.edu/~morrow/papers/GVmathThesis2010.pdf for instance (sorry I do not have all the details in mind for a complete answer here) – Picaud Vincent Mar 15 at 15:11
• The compiled implementation you got from my answer is a genuine implementation of the Voigt profile, nothing "pseudo" about it. – J. M. is away Mar 15 at 15:17
• Thanks J.M. to point that imprecision out. It's a numerical approximation of the Voigt profile and not a Pseudo-Voigt profile. – Adrian Beckert Mar 18 at 12:40

You can use CompilePrint to check what Compile actually did:

Needs["CompiledFunctionTools"];
voigtprofile1 =
Compile[{{δ, _Real, 0}, {σ, _Real, 0}, {A, _Real,
0}, {ν0, _Real, 0}, {ν, _Real, 0}},
A PDF[VoigtDistribution[Abs@δ, Abs@σ], ν - ν0]
];
CompilePrint[voigtprofile1]

Out= "5 arguments
7 Real registers
Underflow checking off
Overflow checking off
Integer overflow checking on
RuntimeAttributes -> {}
R0 = A1
R1 = A2
R2 = A3
R3 = A4
R4 = A5
Result = R6
1   R5 = MainEvaluate[ Function[{δ, σ, A, ν0, ν},
PDF[VoigtDistribution[Abs[δ], Abs[σ]], ν - ν0]][ R0, R1, R2, R3, R4]]
2   R6 = R2 * R5
3   Return"


As you can see, all the function really does is call MainEvaluate and then multiply two numbers together. Basically, you haven't compiled anything.

Instead, evaluate the PDF inside of the Compile to get a symbolic expression that is compilable:

voigtprofile2 = Compile[
{{δ, _Real, 0}, {σ, _Real, 0}, {A, _Real, 0}, {ν0, _Real, 0}, {ν, _Real, 0}},
Evaluate[
A PDF[VoigtDistribution[Abs@δ, Abs@σ], ν - ν0]
],
CompilationOptions -> {"ExpressionOptimization" -> True},
RuntimeAttributes -> {Listable}
];


This should be significantly faster.

• The problem with your voigtprofile2 is that Erfc[] only works for real arguments within a compiled function, so it gets wrapped in MainEvaluate. That was one reason why I wrote the routine in the thread I linked to in the comments. – J. M. is away Mar 15 at 16:32
• Ok, that's handy to know. That's a specific problem of VoigtDistribution I didn't anticipate. I still wanted to point out the importance of using CompilePrint`, though. – Sjoerd Smit Mar 15 at 16:39