How to speed up Voigt profile fits?

Question

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[[1]]/tg[[1]]

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.

Answer 1

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[2]],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.

Answer 2

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[26]= "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.