It would be nice if the C functions Rf_logspace_sum, Rf_logspace_add, and Rf_logspace_sub were available as R functions. (I wish the '_sub' were '_subtract' because 'sub' means too many things in R.)
I think Rf_logspace_sum in R could be a little better. E.g., using the C code #include <R.h> #include <Rinternals.h> #include <Rmath.h> SEXP Call_logspace_sum(SEXP x) { if (TYPEOF(x) != REALSXP) { Rf_error("'x' must be a numeric vector"); } return ScalarReal(Rf_logspace_sum(REAL(x), length(x))); } and the R functions logspace_sum <- function (x) .Call("Call_logspace_sum", as.numeric(x)) and test <- function (x) { x <- as.numeric(x) rbind(Rmpfr = as.numeric(log(sum(exp(Rmpfr::mpfr(x, precBits=5000))))), Rf_logspace_sum = logspace_sum(x), subtract_xmax = log(sum(exp(x - max(x)))) + max(x), naive = log(sum(exp(x)))) } R-3.3.2 on Linux gives, after options(digits=17) > test(c(0, -50)) [,1] Rmpfr 1.9287498479639178e-22 Rf_logspace_sum 1.9287498479639178e-22 subtract_xmax 0.0000000000000000e+00 naive 0.0000000000000000e+00 which is nice, but also the not so nice > test(c(0, -50, -50)) [,1] Rmpfr 3.8574996959278356e-22 Rf_logspace_sum 0.0000000000000000e+00 subtract_xmax 0.0000000000000000e+00 naive 0.0000000000000000e+00 With TERR the second test has Rmpfr==Rf_logspace_sum for that example. Bill Dunlap TIBCO Software wdunlap tibco.com On Mon, Nov 7, 2016 at 3:08 AM, Martin Maechler <maech...@stat.math.ethz.ch> wrote: > >>>>> William Dunlap via R-help <r-help@r-project.org> > >>>>> on Sun, 6 Nov 2016 20:53:17 -0800 writes: > > > Perhaps the C function Rf_logspace_sum(double *x, int n) would help > in > > computing log(b). It computes log(sum(exp(x_i))) for i in 1..n, > avoiding > > unnecessary under- and overflow. > > Indeed! > > I had thought more than twice to also export it to the R level > notably as we have been using two R level versions in a package > I maintain ('copula'). They are vectorized there in a way that > seemed particularly useful to our (Marius Hofert and my) use cases. > > More on this -- making these available in R, how exactly? -- > probably should move to the R-devel list. > > Thank you Bill for bringing it up! > Martin > > > Bill Dunlap > > TIBCO Software > > wdunlap tibco.com > > > On Sun, Nov 6, 2016 at 5:25 PM, Rolf Turner <r.tur...@auckland.ac.nz> > wrote: > > >> On 07/11/16 13:07, William Dunlap wrote: > >> > >>> Have you tried reparameterizing, using logb (=log(b)) instead of b? > >>> > >> > >> Uh, no. I don't think that that makes any sense in my context. > >> > >> The "b" values are probabilities and must satisfy a "sum-to-1" > >> constraint. To accommodate this constraint I re-parametrise via a > >> "logistic" style parametrisation --- basically > >> > >> b_i = exp(z_i)/[sum_j exp(z_j)], j = 1, ... n > >> > >> with the parameters that the optimiser works with being z_1, ..., > z_{n-1} > >> (and with z_n == 0 for identifiability). The objective function is > of the > >> form sum_i(a_i * log(b_i)), so I transform back > >> from the z_i to the b_i in order calculate the value of the > objective > >> function. But when the z_i get moderately large-negative, the b_i > become > >> numerically 0 and then log(b_i) becomes -Inf. And the optimiser > falls over. > >> > >> cheers, > >> > >> Rolf > >> > >> > >>> Bill Dunlap > >>> TIBCO Software > >>> wdunlap tibco.com <http://tibco.com> > >>> > >>> On Sun, Nov 6, 2016 at 1:17 PM, Rolf Turner < > r.tur...@auckland.ac.nz > >>> <mailto:r.tur...@auckland.ac.nz>> wrote: > >>> > >>> > >>> I am trying to deal with a maximisation problem in which it is > >>> possible for the objective function to (quite legitimately) return > >>> the value -Inf, which causes the numerical optimisers that I have > >>> tried to fall over. > >>> > >>> The -Inf values arise from expressions of the form "a * log(b)", > >>> with b = 0. Under the *starting* values of the parameters, a must > >>> equal equal 0 whenever b = 0, so we can legitimately say that a * > >>> log(b) = 0 in these circumstances. However as the maximisation > >>> algorithm searches over parameters it is possible for b to take the > >>> value 0 for values of > >>> a that are strictly positive. (The values of "a" do not change > during > >>> this search, although they *do* change between "successive > searches".) > >>> > >>> Clearly if one is *maximising* the objective then -Inf is not a > value > >>> of > >>> particular interest, and we should be able to "move away". But the > >>> optimising function just stops. > >>> > >>> It is also clear that "moving away" is not a simple task; you can't > >>> estimate a gradient or Hessian at a point where the function value > >>> is -Inf. > >>> > >>> Can anyone suggest a way out of this dilemma, perhaps an optimiser > >>> that is equipped to cope with -Inf values in some sneaky way? > >>> > >>> Various ad hoc kludges spring to mind, but they all seem to be > >>> fraught with peril. > >>> > >>> I have tried changing the value returned by the objective function > >>> from > >>> "v" to exp(v) --- which maps -Inf to 0, which is nice and finite. > >>> However this seemed to flatten out the objective surface too much, > >>> and the search stalled at the 0 value, which is the antithesis of > >>> optimal. > >>> > >>> The problem arises in a context of applying the EM algorithm where > >>> the M-step cannot be carried out explicitly, whence numerical > >>> optimisation. > >>> I can give more detail if anyone thinks that it could be relevant. > >>> > >>> I would appreciate advice from younger and wiser heads! :-) > >>> > >>> cheers, > >>> > >>> Rolf Turner > >>> > >>> -- > >>> Technical Editor ANZJS > >>> Department of Statistics > >>> University of Auckland > >>> Phone: +64-9-373-7599 ext. 88276 <tel:%2B64-9-373-7599%20ext.%2 > 088276> > >>> > >>> ______________________________________________ > >>> R-help@r-project.org <mailto:R-help@r-project.org> mailing list -- > >>> To UNSUBSCRIBE and more, see > >>> https://stat.ethz.ch/mailman/listinfo/r-help > >>> <https://stat.ethz.ch/mailman/listinfo/r-help> > >>> PLEASE do read the posting guide > >>> http://www.R-project.org/posting-guide.html > >>> <http://www.R-project.org/posting-guide.html> > >>> and provide commented, minimal, self-contained, reproducible code. > >>> > >>> > >>> > >> > >> -- > >> Technical Editor ANZJS > >> Department of Statistics > >> University of Auckland > >> Phone: +64-9-373-7599 ext. 88276 > >> > > > [[alternative HTML version deleted]] > > > ______________________________________________ > > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see > > 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. > [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.