R/Function_evaluation.R

In MMAD: An R Package of Minorization-Maximization Algorithm via the Assembly--Decomposition Technology

#' Function Evaluation
#'
#' This function allows the user to compute the value, gradients and the Hessian matrix of the target function \eqn{f(\boldsymbol{\theta})} at a given point \eqn{\boldsymbol{\theta}^*}.
#'
#' @param Function_obj An R list depicting the target function.
#' @param input The value \eqn{\boldsymbol{\theta}^*} at which the target function is evaluated.
#'
#' @return
#' An R list recording the value, gradients and the Hessian matrix of the target function \eqn{f(\boldsymbol{\theta})} at a given point \eqn{\boldsymbol{\theta}^*}.
#'
#' @export
Function_evaluation<-function(Function_obj,input)
{
  if(length(input)!=Function_obj$dimension)
  {
    print("Dimension mismatch!")
  }
  else
  {
    p<-length(input)
    value_output<-rep(0,length(Function_obj$components))
    gradient_output<-matrix(0,length(Function_obj$components),p)
    Hessian_output<-matrix(0,length(Function_obj$components)*p,p)
    for(i in 1:length(Function_obj$components))
    {
      ## Put value in the function argument
      update_parameter<-rep(NA,length(Function_obj$components[[i]]$parameter))

      zero_which<-which(Function_obj$components[[i]]$parameter==0)
      neg<-which(Function_obj$components[[i]]$parameter<0)
      neg_value<--Function_obj$components[[i]]$parameter[neg]
      pos<-which(Function_obj$components[[i]]$parameter>0)
      pos_value<-Function_obj$components[[i]]$parameter[pos]

      update_parameter[neg]<-input[neg_value]
      update_parameter[pos]<-value_output[pos_value]
      update_parameter[zero_which]<-1

      ## Compute the value
      value_output[i]<-get(Function_obj$components[[i]]$functions)(coefficient = Function_obj$components[[i]]$coefficient,parameter = update_parameter)
      ## Compute the gradient
      update_gradient<-matrix(0,length(Function_obj$components[[i]]$parameter),p)
      update_gradient[neg,]<-(matrix(rep(1:p,each=length(neg)),ncol=p)==neg_value)
      update_gradient[pos,]<-gradient_output[pos_value,]
      gradient_mid<-get(paste0("Gradient_",Function_obj$components[[i]]$functions))(coefficient = Function_obj$components[[i]]$coefficient,parameter = update_parameter)
      gradient_output[i,]<-colSums(update_gradient*gradient_mid)
      ## Compute the Hessian
      Hessian_output[(i-1)*p+1:p,]<-t(update_gradient)%*%(get(paste0("Hessian_",Function_obj$components[[i]]$functions))(coefficient = Function_obj$components[[i]]$coefficient,parameter = update_parameter))%*%update_gradient
      if(length(pos_value)>0)
      {
        for(j in 1:length(pos_value))
        {
          Hessian_output[(i-1)*p+1:p,]<-gradient_mid[pos[j]]*Hessian_output[(pos_value[j]-1)*p+1:p,]+Hessian_output[(i-1)*p+1:p,]
        }
      }
    }
    return(list(Value=value_output[i],Gradient=gradient_output[i,],Hessian=Hessian_output[(i-1)*p+1:p,]))
  }
}