R/fn.lnfgrad.R

In AdMit: Adaptive Mixture of Student-t Distributions

## Function which builds the objective function and its gradient for the
## optimization of the probabilities, both evaluated at the same 'lambda'
## __input__
## lnK      : [NpxH matrix] of kernel values
## lnD      : [NpxH^2 matrix] of Student-t densities components
## __output__
## [list] with the following components:
## $p       : [function] probabilities corresponding to 'lambda'
## $f       : [function] objective function at 'lambda'
## $g       : [function] gradient of the objective function at 'lambda'
## __20080427__
'fn.lnfgrad' <- function(lnK, lnD)
{
  Np <- nrow(lnK)
  H <- ncol(lnK)

  ## The objective is log(f1) - 2 log(f2) with f1 and f2 weighted moments of
  ## exp(lnK)/mixture, so adding a constant to every element of 'lnK' leaves it
  ## and its gradient unchanged. Working from the raw values, however, makes
  ## exp(2 * lnw) overflow above about +350 and underflow below about -350, and
  ## a log-kernel of a few hundred is ordinary for any posterior with more than
  ## a handful of observations. Centring restores the invariance in practice.
  shift <- if (any(is.finite(lnK))) max(lnK[is.finite(lnK)]) else 0
  lnKvec <- as.double(as.vector(t(lnK-shift)))
  lnDvec <- as.double(as.vector(t(lnD)))

  ## function which transforms the probabilities (positivity and summability),
  ## through a max-shifted softmax so that a large 'lambda' cannot overflow
  'fn.lambdap' <- function(lambda)
  {
    e <- c(lambda,0)
    e <- exp(e-max(e))
    as.vector(e/sum(e))
  }

  ## 'fnlnf_C' returns the objective function and its gradient in one pass.
  ## Both are cached on the value of 'lambda': the optimizer normally asks for
  ## the objective and the gradient at the same point, but nothing guarantees
  ## it, and a gradient left over from another point would silently mislead
  ## the optimizer.
  memo <- new.env(hash = FALSE)
  memo$lambda <- NULL

  'fn.eval' <- function(lambda)
  {
    if (is.null(memo$lambda) || !identical(memo$lambda, lambda))
    {
      r <- .C('fnlnf_C',
              lnp = as.double(log(fn.lambdap(lambda))),
              lnk = lnKvec,
              lnd = lnDvec,
              Np = as.integer(Np),
              H = as.integer(H),
              f = as.double(0),
              grad = vector('double',H),
              PACKAGE = 'AdMit',
              NAOK = TRUE)

      memo$lambda <- lambda
      memo$f <- as.numeric(r$f)
      memo$grad <- r$grad
    }
    invisible(NULL)
  }

  list(
    p = fn.lambdap,

    ## objective function
    f = function(lambda)
    {
      fn.eval(lambda)
      memo$f
    },

    ## gradient of the objective function with respect to 'lambda'
    g = function(lambda)
    {
      fn.eval(lambda) ## the gradient at 'lambda', never one left over
      e <- c(exp(lambda),1)
      s <- sum(e)
      tmp <- -e %*% t(e) / s^2
      diag(tmp) <- (e*s-e^2) / s^2
      gradlambda <- as.matrix(tmp[1:H,1:(H-1)])

      as.vector(t(gradlambda)%*%memo$grad)
    })
}