R/Function_minorization.R

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

#' Function Minorization
#'
#' This function allows the user to minorize 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 minorized.
#'
#' @return
#' An R list recording the surrogate function \eqn{S(\boldsymbol{\theta}|\boldsymbol{\theta}^*)=S_C+\sum_{i=1}^pS_i(\theta_{i}|\boldsymbol{\theta}^*)} where
#' \itemize{
#'   \item the first \eqn{p} objects (named Surrogate_1, Surrogate_2, ..., Surrogate_p) record \eqn{S_1(\theta_{1}|\boldsymbol{\theta}^*),\ldots,S_p(\theta_{p}|\boldsymbol{\theta}^*)};
#'   \item the final object (named Constant) records the constant \eqn{S_C}.
#' }
#'
#' @export
Function_minorization<-function(Function_obj,input)
{
  if(length(input)!=Function_obj$dimension)
  {
    print("Dimension mismatch!")
  }
  else
  {
    p<-length(input)
    mid_value<-rep(0,length(Function_obj$components))
    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]<-mid_value[pos_value]
      update_parameter[zero_which]<-1

      ## Compute the value
      mid_value[i]<-get(Function_obj$components[[i]]$functions)(coefficient = Function_obj$components[[i]]$coefficient,parameter = update_parameter)
    }

    Surrogate_function_mid<-list()
    Surrogate_function_final<-list()
    Surrogate_function_final_coefficient<-list()
    for(i in 1:Function_obj$dimension)
    {
      Surrogate_function_final[[i]]<-list()
      Surrogate_function_final[[i]]$dimension<-1
      Surrogate_function_final[[i]]$components<-list()
      Surrogate_function_final_coefficient[[i]]<-list(coefficient=NULL,parameter=NULL,functions="Linear_combination")
    }
    Surrogate_constant<-0

    coefficients<-1
    Surrogate_function_mid[[1]]<-Function_obj$components[[length(Function_obj$components)]]

    while (length(Surrogate_function_mid)>0)
    {
      Function_component<-Surrogate_function_mid[[1]]
      if((length(Function_component$parameter)==1)&(Function_component$parameter[1]<0))
      {
        ii<--Function_component$parameter

        add_position<-length(Surrogate_function_final[[ii]]$components)+1
        add_function<-list(coefficient=Function_component$coefficient,parameter=-1,functions=Function_component$functions)
        Surrogate_function_final[[ii]]$components[[add_position]]<-add_function
        Surrogate_function_final_coefficient[[ii]]$coefficient<-c(Surrogate_function_final_coefficient[[ii]]$coefficient,coefficients[1])
        Surrogate_function_final_coefficient[[ii]]$parameter<-c(Surrogate_function_final_coefficient[[ii]]$parameter,add_position)
      }
      else
      {
        convexity_indicator<-convexity(Function_component$functions,coefficients[1])
        if((Function_component$functions=="Power")&(abs(Function_component$coefficient[length(Function_component$coefficient)]-1/2)<1/2))
        {
          convexity_indicator<-3-convexity_indicator
        }

        if(convexity_indicator==0)
        {
          for(j in 1:length(Function_component$parameter))
          {
            if(Function_component$parameter[j]<0)
            {
              ii<--Function_component$parameter[j]
              Surrogate_function_final_coefficient[[ii]]$coefficient<-c(Surrogate_function_final_coefficient[[ii]]$coefficient,coefficients[1]*Function_component$coefficient[j])
              Surrogate_function_final_coefficient[[ii]]$parameter<-c(Surrogate_function_final_coefficient[[ii]]$parameter,-1)
            }
            if(Function_component$parameter[j]==0)
            {
              Surrogate_constant<-Surrogate_constant+coefficients[1]*Function_component$coefficient[j]
            }
            if(Function_component$parameter[j]>0)
            {
              Surrogate_function_mid[[length(Surrogate_function_mid)+1]]<-Function_obj$components[[Function_component$parameter[j]]]
              coefficients<-c(coefficients,coefficients[1]*Function_component$coefficient[j])
            }
          }
        }
        if(convexity_indicator==1)
        {
          update_parameter<-rep(1,length(Function_component$parameter))
          update_parameter[Function_component$parameter<0]<-input[-Function_component$parameter[Function_component$parameter<0]]
          update_parameter[Function_component$parameter>0]<-mid_value[Function_component$parameter[Function_component$parameter>0]]
          GG<-get(paste0("Gradient_",Function_component$functions))(coefficient = Function_component$coefficient,parameter = update_parameter)
          VV<-get(Function_component$functions)(coefficient = Function_component$coefficient,parameter = update_parameter)
          Surrogate_constant<-Surrogate_constant+coefficients[1]*(VV-sum(GG*update_parameter))
          power_indicator<-(Function_component$functions=="Power")
          for(j in 1:length(Function_component$parameter))
          {
            if(Function_component$parameter[j]==0)
            {
              Surrogate_constant<-Surrogate_constant+coefficients[1]*GG[j]*update_parameter[j]
            }
            if(Function_component$parameter[j]<0)
            {
              ii<--Function_component$parameter[j]

              Surrogate_function_final_coefficient[[ii]]$coefficient<-c(Surrogate_function_final_coefficient[[ii]]$coefficient,coefficients[1]*GG[j])
              Surrogate_function_final_coefficient[[ii]]$parameter<-c(Surrogate_function_final_coefficient[[ii]]$parameter,-1)
            }
            if(Function_component$parameter[j]>0)
            {
              Surrogate_function_mid[[length(Surrogate_function_mid)+1]]<-Function_obj$components[[Function_component$parameter[j]]]
              coefficients<-c(coefficients,coefficients[1]*GG[j])
            }
          }
        }
        if(convexity_indicator==2)
        {
          update_parameter<-rep(1,length(Function_component$parameter))
          update_parameter[Function_component$parameter<0]<-input[-Function_component$parameter[Function_component$parameter<0]]
          update_parameter[Function_component$parameter>0]<-mid_value[Function_component$parameter[Function_component$parameter>0]]
          linear_component<-update_parameter*Function_component$coefficient[1:length(update_parameter)]
          linear_sum<-sum(linear_component)
          power_indicator<-(Function_component$functions=="Power")
          for(j in 1:length(Function_component$parameter))
          {
            if(Function_component$parameter[j]==0)
            {
              Surrogate_constant<-Surrogate_constant+coefficients[1]*linear_component[j]/linear_sum*get(Function_component$functions)(coefficient = linear_sum,parameter = 1)
            }
            if(Function_component$parameter[j]<0)
            {
              ii<--Function_component$parameter[j]

              add_position<-length(Surrogate_function_final[[ii]]$components)+1
              add_function<-list(coefficient=Function_component$coefficient[j]*linear_sum/linear_component[j],parameter=-1,functions=Function_component$functions)
              if(power_indicator)
              {
                add_function$coefficient<-c(add_function$coefficient,Function_component$coefficient[length(Function_component$coefficient)])
              }

              Surrogate_function_final[[ii]]$components[[add_position]]<-add_function

              Surrogate_function_final_coefficient[[ii]]$coefficient<-c(Surrogate_function_final_coefficient[[ii]]$coefficient,coefficients[1]*linear_component[j]/linear_sum)
              Surrogate_function_final_coefficient[[ii]]$parameter<-c(Surrogate_function_final_coefficient[[ii]]$parameter,add_position)
            }
            ###### Check point
            if(Function_component$parameter[j]>0)
            {
              outer_function<-Function_component$functions
              internal_function<-Function_obj$components[[Function_component$parameter[j]]]

              coefficient_general<-coefficients[1]*linear_component[j]/linear_sum
              coefficient_out<-Function_component$coefficient[j]*linear_sum/linear_component[j]
              coefficient_in<-internal_function$coefficient

              parameter_in<-internal_function$parameter

              if(outer_function=="Logarithmic")
              {
                if(internal_function$functions=="Power")
                {
                  add_function<-list(coefficient=coefficient_in[-length(coefficient_in)],parameter=internal_function$parameter,functions="Logarithmic")
                  Surrogate_function_mid[[length(Surrogate_function_mid)+1]]<-add_function
                  coefficients<-c(coefficients,coefficient_general*coefficient_in[length(coefficient_in)])
                  Surrogate_constant<-Surrogate_constant+coefficient_general*log(coefficient_out)
                }
              }

              if(outer_function=="Logarithmic")
              {
                if(internal_function$functions=="Exponential")
                {
                  add_function<-list(coefficient=coefficient_in,parameter=internal_function$parameter,functions="Linear_combination")
                  Surrogate_function_mid[[length(Surrogate_function_mid)+1]]<-add_function
                  coefficients<-c(coefficients,coefficient_general)
                  Surrogate_constant<-Surrogate_constant+coefficient_general*log(coefficient_out)
                }
              }

              if(outer_function=="Exponential")
              {
                if(internal_function$functions=="Logarithmic")
                {
                  if(coefficient_out!=1)
                  {
                    add_function<-list(coefficient=c(coefficient_in,coefficient_out),parameter=internal_function$parameter,functions="Power")
                    Surrogate_function_mid[[length(Surrogate_function_mid)+1]]<-add_function
                    coefficients<-c(coefficients,coefficient_general)
                  }
                  if(coefficient_out==1)
                  {
                    add_function<-list(coefficient=c(coefficient_in),parameter=internal_function$parameter,functions="Linear_combination")
                    Surrogate_function_mid[[length(Surrogate_function_mid)+1]]<-add_function
                    coefficients<-c(coefficients,coefficient_general)
                  }
                }
              }
            }
          }
        }
        if(convexity_indicator==3)
        {
          update_parameter<-rep(1,length(Function_component$parameter))
          update_parameter[Function_component$parameter<0]<-input[-Function_component$parameter[Function_component$parameter<0]]
          update_parameter[Function_component$parameter>0]<-mid_value[Function_component$parameter[Function_component$parameter>0]]
          linear_component<-update_parameter*Function_component$coefficient[1:length(update_parameter)]
          linear_sum<-sum(linear_component)

          VV<-get(Function_component$functions)(coefficient = Function_component$coefficient,parameter = update_parameter)
          Surrogate_constant<-Surrogate_constant+coefficients[1]*(VV+linear_sum/(1-linear_sum)*log(linear_sum))

          Surrogate_function_mid[[length(Surrogate_function_mid)+1]]<-list(coefficient=Function_component$coefficient,parameter=Function_component$parameter,functions="Logarithmic")
          coefficients<-c(coefficients,-coefficients[1]*linear_sum/(1-linear_sum))
        }
      }
      Surrogate_function_mid<-Surrogate_function_mid[-1]
      coefficients<-coefficients[-1]
      # print(length(Surrogate_function_mid))
    }

    for(i in 1:Function_obj$dimension)
    {
      add_position<-length(Surrogate_function_final[[i]]$components)+1
      Surrogate_function_final[[i]]$components[[add_position]]<-Surrogate_function_final_coefficient[[i]]
    }


    Surrogate_function_final$Constant<-Surrogate_constant
    names(Surrogate_function_final)[1:Function_obj$dimension]<-paste0("Surrogate_",1:Function_obj$dimension)
    return(Surrogate_function_final)
  }
}