[R] Optimize function in R: unable to find maximum of logistic function

[William Dunlap]

> One can do better yet by using the logspace_add(logx, logy) function
> available in the C API to R, but there is currently not an R-language interface
> to that (or logspace_subtract(logx,logy)).  They calculate log(exp(logx) +- exp(logy))
> with minimal roundoff error.
> 
> You could probably also use the built-in plogis() function, with its log=TRUE and
> perhaps lower.tail=FALSE arguments.

Here are functions, "f2" and "f3", which compute log(f()), via the above algorithms,
along with your original "f" and my first suggestion "f1".  f2 and f3 give almost
identical results over a wide range of inputs (pos. and neg.).  f1 agrees with them
for positive inputs but falls apart for larger negative inputs.  Any of f1, f2, or f3 will
do your optimization, although I don't know why you want to optimize a monotonic
function.

f <-
function (k, T_s = 20) {
    2 - 2/(1 + exp(-2*T_s*k))
}
f1 <-
function(k, T_s = 20) {
    # log(f(k))
    log(2) - 2*T_s*k - log1p(exp(-2*T_s*k))
}
f2 <-
function(k, T_s = 20) {
    # log(f(k))
    log(2) - logspace_add(2*T_s*k, 0)
}
logspace_add <-
function(logx, logy)
{
    # log(exp(logx) + exp(logy))
    mn <- pmin(logx, logy)
    mx <- pmax(logx, logy)
    mx + log1p(exp(mn - mx))
}
f3 <-
function(k, T_s = 20) {
    # log(f(k))
    log(2) + plogis(k, lower.tail=FALSE, log=TRUE, scale=1/(T_s*2))
}

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com


> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf
> Of William Dunlap
> Sent: Monday, October 28, 2013 3:03 PM
> To: Rolf Turner; Shantanu MULLICK
> Cc: r-help at r-project.org
> Subject: Re: [R] Optimize function in R: unable to find maximum of logistic function
> 
> Note that f(x) returns exactly zero for all x above 1, hence its estimated
> derivative is 0 everywhere in that region.   You are asking optimize to find
> the non-zero part of a function which is 0 in more than 80% of its domain.
> 
> You can put a trace on f() to see where optimize() looks:
>    > trace(f, quote(cat("k=", deparse(k), "\n")), print=FALSE)
>   [1] "f"
>   > optimize(f, c(0, 5), tol = 1e-14, maximum= TRUE)
>   k= 1.90983005625053
>   k= 3.09016994374947
>   k= 3.81966011250105
>   k= 4.27050983124842
>   k= 4.54915028125263
>   k= 4.72135954999579
>   k= 4.82779073125683
>   k= 4.89356881873896
>   ... eventually homing in on 5 ...
> It looks in a few places, but f() is 0 in all of them so it figures that
> is the mininum and the maximum value it ever takes.  I suppose
> you could ask that optimize look in more places for a place with
> a non-zero derivative but it would be very expensive to look in all
> (c. 2^50) places.
> 
> You probably should have it maximize the log of your objective
> function, which gives you more range to play with, and you will
> have to do some numerical analysis to minimize underflow and
> roundoff error.  E.g., one could use
>     f1 <- function (k)  {
>         # log(f(k))
>         T_s <- 20
>         log(2) - 2 * T_s * k - log1p(exp(-2 * T_s * k))
>    }
> where log1p(x) calculates log(1+x) more accurately than a computer's
> 'log' and '+' can.  (Remember that 1+10^-17 exactly equals 1 here.)
> 
> One can do better yet by using the logspace_add(logx, logy) function
> available in the C API to R, but there is currently not an R-language interface
> to that (or logspace_subtract(logx,logy)).  They calculate log(exp(logx) +- exp(logy))
> with minimal roundoff error.
> 
> You could probably also use the built-in plogis() function, with its log=TRUE and
> perhaps lower.tail=FALSE arguments.
> 
> Bill Dunlap
> Spotfire, TIBCO Software
> wdunlap tibco.com
> 
> 
> > -----Original Message-----
> > From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf
> > Of Rolf Turner
> > Sent: Monday, October 28, 2013 2:01 PM
> > To: Shantanu MULLICK
> > Cc: r-help at r-project.org
> > Subject: Re: [R] Optimize function in R: unable to find maximum of logistic function
> >
> >
> > This could be described as a bug, perhaps.  Or it could be described as
> > an indication that numerical optimization is inevitably tricky. Notice that
> > if you narrow down your search interval from [0,5] to [0,0.5] you get the
> > right answer:
> >
> >  > optimize(f, c(0, 0.5), maximum= TRUE,tol=1e-10)
> > $maximum
> > [1] 4.192436e-11
> >
> > $objective
> > [1] 1
> >
> > I guess there's a problem with finding a gradient that is effectively
> > (numerically) zero when "k" is equal to 5.
> >
> >
> >      cheers,
> >
> >      Rolf Turner
> >
> >
> > On 10/29/13 06:00, Shantanu MULLICK wrote:
> > > Hello Everyone,
> > >
> > > I want to perform a 1-D optimization by using the optimize() function. I
> > > want to find the maximum value of a "logistic" function. The optimize()
> > > function gives the wrong result.
> > >
> > > My code:
> > > f= function (k) {
> > > T_s = 20
> > > result = (2- 2/(1+ exp(-2*T_s*k)))
> > > return(result)
> > > }
> > > optimize(f, c(0, 5), tol = 0.0000000000001, maximum= TRUE)
> > >
> > > The maximum value for the function happens at k=0, and the maximum value is
> > > 1. Yet  the optimise function, says that the maximum value happens at k=
> > > 4.9995, and the maximum value is 0.
> > >
> >
> > ______________________________________________
> > R-help at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.