R/linear_DSPP.R

Defines functions dspp_compute_Sw dspp_distnorm dspp_compute_wi do.dspp

Documented in do.dspp

#' Discriminative Sparsity Preserving Projection
#'
#' Discriminative Sparsity Preserving Projection (DSPP) is a supervised dimension reduction method
#' that employs sparse representation model to adaptively build both intrinsic adjacency graph and
#' penalty graph. It follows an integration of global within-class structure into manifold learning
#' under exploiting discriminative nature provided from label information.
#'
#'
#'
#' @param X an \eqn{(n\times p)} matrix or data frame whose rows are observations.
#' @param label a length-\eqn{n} vector of data class labels.
#' @param ndim an integer-valued target dimension.
#' @param preprocess  an additional option for preprocessing the data.
#' Default is "center". See also \code{\link{aux.preprocess}} for more details.
#' @param lambda regularization parameter for constructing sparsely weighted network.
#' @param rho a parameter for balancing the local and global contribution.
#'
#'
#' @return a named list containing
#' \describe{
#' \item{Y}{an \eqn{(n\times ndim)} matrix whose rows are embedded observations.}
#' \item{trfinfo}{a list containing information for out-of-sample prediction.}
#' \item{projection}{a \eqn{(p\times ndim)} whose columns are basis for projection.}
#' }
#'
#' @examples
#' \dontrun{
#' ## use iris data
#' data(iris)
#' set.seed(100)
#' subid = sample(1:150, 50)
#' X     = as.matrix(iris[subid,1:4])
#' label = as.factor(iris[subid,5])
#'
#' ## try different rho values
#' out1 <- do.dspp(X, label, ndim=2, rho=0.01)
#' out2 <- do.dspp(X, label, ndim=2, rho=0.1)
#' out3 <- do.dspp(X, label, ndim=2, rho=1)
#'
#' ## visualize
#' opar <- par(no.readonly=TRUE)
#' par(mfrow=c(1,3))
#' plot(out1$Y, main="rho=0.01", col=label, pch=19)
#' plot(out2$Y, main="rho=0.1",  col=label, pch=19)
#' plot(out3$Y, main="rho=1",    col=label, pch=19)
#' par(opar)
#' }
#'
#' @references
#' \insertRef{gao_discriminative_2015}{Rdimtools}
#'
#' @author Kisung You
#' @rdname linear_DSPP
#' @concept linear_methods
#' @export
do.dspp <- function(X, label, ndim=2, preprocess=c("center","scale","cscale","decorrelate","whiten"),
                    lambda=1.0, rho=1.0){
  #------------------------------------------------------------------------
  ## PREPROCESSING
  #   1. data matrix
  aux.typecheck(X)
  n = nrow(X)
  p = ncol(X)
  #   2. label : check and return a de-factored vector
  #   For this example, there should be no degenerate class of size 1.
  label  = check_label(label, n)
  ulabel = unique(label)
  for (i in 1:length(ulabel)){
    if (sum(label==ulabel[i])==1){
      stop("* do.dspp : no degerate class of size 1 is allowed.")
    }
  }
  if (any(is.na(label))||(any(is.infinite(label)))){stop("* Supervised Learning : any element of 'label' as NA or Inf will simply be considered as a class, not missing entries.")  }
  #   3. ndim
  ndim = as.integer(ndim)
  if (!check_ndim(ndim,p)){stop("* do.dspp : 'ndim' is a positive integer in [1,#(covariates)).")}
  #   4. preprocess
  if (missing(preprocess)){
    algpreprocess = "center"
  } else {
    algpreprocess = match.arg(preprocess)
  }
  #   5. lambda : positive real number
  lambda = as.double(lambda)
  if (!check_NumMM(lambda, 0, 1e+10, compact=FALSE)){stop("* do.dspp : 'lambda' is a positive real number.")}
  #   6. rho    : nonnegative real number
  rho = as.double(rho)
  if (!check_NumMM(rho, 0, 1e+10, compact=TRUE)){stop(" do.dspp : 'rho' is a nonnegative real number.")}

  #------------------------------------------------------------------------
  ## COMPUTATION 1. Preliminary Computation
  tmplist = aux.preprocess.hidden(X,type=algpreprocess,algtype="linear")
  trfinfo = tmplist$info
  pX      = tmplist$pX
  #------------------------------------------------------------------------
  ## COMPUTATION 2. Three Important Matrices
  #   1. W : adjacency, intrinsic
  W = array(0,c(n,n))
  for (i in 1:n){
    #   1-1. target vector
    xi = matrix(pX[i,],ncol=1)
    #   1-2. target matrix
    tgtidx = setdiff(which(label==label[i]),i)
    Xi = t(pX[tgtidx,])
    #   1-3. compute wi and assign it as column vector
    W[tgtidx,i] = dspp_compute_wi(xi,Xi,lambda)
  }
  #   2. B : adjacency, penalty
  #   2-1. compute radius
  radius = rep(0,n)
  for (i in 1:n){
    tgtrow  = as.vector(W[i,])
    idxmax  = which(tgtrow==max(tgtrow))
    if (length(idxmax)>1){
      idxmax = as.integer(idxmax[1])
    }
    maxdist = as.vector(pX[i,])-as.vector(pX[idxmax,])
    radius[i] = sqrt(sum(maxdist*maxdist))
  }
  #   2-2. let's iterate
  B = array(0,c(n,n))
  for (i in 1:n){
    veci = as.vector(pX[i,])
    for (j in 1:n){
      vecj = as.vector(pX[j,])
      if (i!=j){
        if (dspp_distnorm(veci,vecj)<=radius[i]){
          if (label[i]!=label[j]){
            B[i,j] = 1
          }
        }
      }
    }
  }
  #   3. compute class-wise mean
  classmean = array(0,c(length(ulabel), p))
  for (i in 1:length(ulabel)){
    partidx = which(label==ulabel[i])
    classmean[i,] = colMeans(pX[partidx,])
  }
  #   4. Sw : does not have the name on it.
  Sw = array(0,c(p,p))
  for (i in 1:n){
    #   4-1. select one vector
    cvec = pX[i,]
    clabel = label[i]
    #   4-2. select target
    tgtmat = pX[-i,]
    tgtlabel = label[-i]
    #   4-3. compute inner summation
    Smat = dspp_compute_Sw(cvec,clabel,tgtmat,tgtlabel,ulabel,classmean,as.double(radius[i]))
    Sw   = Sw+Smat
  }
  #------------------------------------------------------------------------
  ## COMPUTATION 3. set geigen setting
  #   0. symmetrize
  W = (W+t(W))/2
  B = (B+t(B))/2
  Sw= (Sw+t(Sw))/2
  #   1. Lw, Lb
  Lw = diag(rowSums(W))-W
  Lb = diag(rowSums(B))-B
  #   2. LHS and RHS
  LHS = ((t(pX)%*%Lw%*%pX) + (rho*Sw))
  RHS = (t(pX)%*%Lb%*%pX)
  #   3. use with SMALLEST ones
  projection = aux.geigen(LHS,RHS,ndim,maximal=FALSE)

  #------------------------------------------------------------------------
  ## RETURN OUTPUT

  result = list()
  result$Y = pX%*%projection
  result$trfinfo = trfinfo
  result$projection = projection
  return(result)
}



#  ------------------------------------------------------------------------
#' @keywords internal
#' @noRd
dspp_compute_wi <- function(xi, Xi, lambda){
  n   = ncol(Xi)
  si  = CVXR::Variable(n)
  obj = (CVXR::p_norm(xi-(Xi%*%si),p=2) + lambda*CVXR::p_norm(si,p=1))
  constr = list(si>=0)
  prob   = CVXR::Problem(Minimize(obj), constr)
  result = solve(prob)
  return(as.vector(result$getValue(si)))
}
#' @keywords internal
#' @noRd
dspp_distnorm <- function(vec1, vec2){
  vecdiff = as.vector(vec1)-as.vector(vec2)
  output  = as.double(sqrt(sum(vecdiff*vecdiff)))
  return(output)
}
#' @keywords internal
#' @noRd
dspp_compute_Sw <- function(cvec, clabel, tgtmat, tgtlabel, ulabel, classmean, radius){
  p = length(cvec)
  Smat = array(0,c(p,p))
  ntgt = nrow(tgtmat)
  for (i in 1:ntgt){
    # 1. within radius one
    vecdiff = as.vector(cvec)-as.vector(tgtmat[i,])
    if (sqrt(sum(vecdiff*vecdiff))<radius){
      # 2. different class only
      if (tgtlabel[i]!=clabel){
        tgtvec  = (tgtmat[i,])
        tgtidx  = which(ulabel==clabel)
        meanvec = (classmean[tgtidx,])
        mvdiff  = as.vector(tgtvec)-as.vector(meanvec)
        # 3. update Smat
        Smat = Smat + outer(mvdiff,mvdiff)
      }
    }
  }
  return(Smat)
}
kisungyou/Rdimtools documentation built on Jan. 2, 2023, 9:55 a.m.