R/acmtf_fun.R

In CMTFtoolbox: Create (Advanced) Coupled Matrix and Tensor Factorization Models

#' Calculate function value of ACMTF
#'
#' @param x Vectorized parameters of the CMTF model.
#' @param Z Z object as generated by [setupCMTFdata()].
#' @param alpha Alpha value of the loss function as specified by Acar et al., 2014
#' @param beta Beta value of the loss function as specified by Acar et al., 2014
#' @param epsilon Epsilon value of the loss function as specified by Acar et al., 2014
#' @param manual Manual calculation of each loss term (default FALSE)
#'
#' @return Scalar of the loss function value (when manual=FALSE), otherwise a list containing all loss terms
#' @export
#'
#' @examples
#' A = array(rnorm(108*2), c(108, 2))
#' B = array(rnorm(100*2), c(100, 2))
#' C = array(rnorm(10*2), c(10, 2))
#' D = array(rnorm(100*2), c(100,2))
#' E = array(rnorm(10*2), c(10,2))
#'
#' df1 = reinflateTensor(A, B, C)
#' df2 = reinflateTensor(A, D, E)
#' datasets = list(df1, df2)
#' modes = list(c(1,2,3), c(1,4,5))
#' Z = setupCMTFdata(datasets, modes, normalize=FALSE)
#'
#' init = initializeACMTF(Z, 2, output="vect")
#' f = acmtf_fun(init, Z)
acmtf_fun = function(x, Z, alpha=1, beta=rep(1e-3, length(Z$object)), epsilon=1e-8, manual=FALSE){

  numDatasets = length(Z$object)
  numModes = max(unlist(Z$modes))
  Fac = vect_to_fac(x, Z)
  numComponents = ncol(Fac[[1]])
  reinflatedBlocks = reinflateFac(Fac, Z, returnAsTensor=TRUE)

  # Penalty for fit on X
  f_per_term = rep(NA, numDatasets)
  for(p in 1:numDatasets){
    modes = Z$modes[[p]]
    reinflatedBlock = reinflatedBlocks[[p]]
    residuals = Z$object[[p]] - reinflatedBlock
    residuals = Z$missing[[p]] * residuals

    Fnorm = rTensor::fnorm(residuals) # verified to work for matrices too
    f_per_term[p] = 0.5 * Fnorm^2
  }

  # Penalty to make the solution norm 1
  f_norm = matrix(NA, nrow=numModes, ncol=numComponents)
  for(i in 1:numModes){
    for(j in 1:numComponents){
      f_norm[i,j] = 0.5 * alpha * (norm(as.matrix(Fac[[i]][,j]), "2")-1)^2
    }
  }

  # Penalty on the lambdas
  f_lambda = matrix(NA, nrow=numComponents, ncol=numDatasets)
  for(i in 1:numComponents){
    for(p in 1:numDatasets){
      lambda_r = Fac[[numModes+1]][p,i]
      f_lambda[i,p] =  0.5 * beta[p] * (sqrt(lambda_r^2 + epsilon))
    }
  }

  if(manual == FALSE){
    f = sum(f_per_term) + sum(f_norm) + sum(f_lambda)
    return(f)
  } else{
    return(list(f_per_term, f_norm, f_lambda))
  }
}