[R] Optimal use of optim?

[Leo Mada]

Dear R-Users,

I am interested in the optimal strategy for optim:

Q: How to formulate the optimization problem?
Q1: Are there benefits for abs(f(x)) vs (f(x))^2?
Q2: Are there any limitations for using abs(...)?

Regarding point 1: my feeling is that the gradients should be more robust with the abs(...) method; but I did not test it.

Regarding point 2: Obviously, abs(f(x)) is not differentiable. This is a limitation from a mathematical point of view, but the gradient should be still computable using a finite difference method. Are there any other limitations?

The question is based on the following example:
solving gamma(x) = y for some given y.

Function multiroot in package rootSolve fails when y < 1. The version using optim fares much better. However, I had to tweak somewhat the parameters in order to get a higher precision. This made me curious about the optimal strategy; unfortunately, I do not have much experience with optimizations.

The example is also available on GitHub:
https://github.com/discoleo/R/blob/master/Math/Integrals.Gamma.Inv.R

Sincerely,

Leonard

### Example:
gamma.invN = function(x, x0, lim, ...,
            ndeps = 1E-5, rtol = 1E-9, gtol = 0.0001, digits = 10) {
      v = x;
      if(length(lim) == 1) lim = c(lim - 1 + gtol, lim - gtol);
      cntr = c(list(...), ndeps = ndeps, pgtol = rtol);
      # abs(): should provide greater precision than ()^2 when computing the gradients;
      # param ndeps: needed to improve accuracy;
      r = optim(x0, \(x) { abs(gamma(x) - v); }, lower = lim[1], upper =lim[2],
            method = "L-BFGS-B", control = cntr);
      if( ! is.null(digits)) print(r$par, digits);
      return(r$par);
}

### Example
Euler   = 0.57721566490153286060651209008240243079;
x = gamma.invN(Euler, -3.5, lim = -3)
gamma(x) - Euler
x = gamma.invN(Euler, -3.9, lim = -3)
gamma(x) - Euler


	[[alternative HTML version deleted]]