Handling absolute value and other discontinuities in numerical optimization methods that use gradients

Question

Suppose we have difficult peak fitting problems where the the users wish to fit asymmetric peaks to their experimental data by the least squares method. One such function is illustrated below:

Here

• $a_0$ = peak area,
• $a_1$ = peak center,
• $a_2$ = peak width,
• $a_3$ = kurtosis, and finally
• $a_4$ = skew.

(See Reference [1])

As one can see, the function contains absolute value terms and Heaviside theta (step) functions, which can lead to discontinuities in the gradient based optimization method. Since this function is available in a commercial peak fitting software, I contacted the software author as to which algorithm handles these least squares problems. The author informed that Levenberg-Marquardt (LM) is used for solving this least squares with starting from good guesses for the variables.

One might assume this would cause problems for optimization methods relying on gradient information, such as the Levenberg-Marquardt algorithm and he stated that in 30 years LM never gave a problem provided we start with a good guess solution i.e. good guesses for $a_0$ to $a_4$. Now this is purely from a numerical point of view.

In numerical optimization, there are several strategies for handling such cases:

1. Using algorithms that can handle discontinuities in the gradient, such as genetic algorithm. In my experience with MATLAB this is quite slow when the number of data points are > 5000, and it may take hours without a solution.

2. Using finite difference approximations to estimate gradients, even if they are not smooth, which are employed in LM.

Given this background, my question as a non-mathematician is: Is using finite differences for gradient estimation "good enough" for handling the absolute value derivative in the context of the Levenberg-Marquardt algorithm? Given small step size, it might be as good as a smooth approximation such those obtained from logistic functions.

If so, what might be the theoretical justifications for this apparent robustness of LM? Are there particular conditions or characteristics of the peak fitting problem that make finite differences adequate despite the non-smooth nature of the peak functions?

I would appreciate any insights or practical references to relevant literature that could clarify why this approach works effectively in practice, even if it seems problematic from a theoretical standpoint.

Answer 1

I think you are missing at least two other options.

Maybe there are no "true" discontinuities

First thing to note is that the use absolute value and a step function does not necessarily imply discontinuous derivatives that would hinder optimization.

For example the step function in the formula seems to be there to just provide zero density above/below some threshold. Since zero density will never be optimal, you may just treat this is a constraint on the space you optimize over.

The use of the absolute value in your formula is a bit more tricky to understand, but it may easily be that both one sided limits of the derivative at the supposed discontinuity are in fact equal and so there is no problem at all and analytical gradients work great. E.g. $|x^3|$ is a function that is easy to minimize with gradient descent, although absolute value is involved.

Ignore the discontinuities completely

A second option is that there indeed are real discontinuities, but you just ignore them and use analytical gradients anyway.

One way to theoretically justify this approach would be that you split the parameter space into multiple regions which are all smooth and you are trying to find a local optimum within the region that includes your initial values.

I would also assume that most LM implementations are able to detect if they get thwarted by a discontinuity so if that does not happen in the optimization run it is safe-ish to assume you were just moving within a continuous region.

Since it is plausible that good guesses about the parameters can be made directly (e.g. presumably sample mean/median is good guess for $a_1$ and sample skew/kurtosis a good guess for $a_3, a_4$), it is also reasonable to expect that you basically always start in the region that contains the global optimum.