R/ingUtils.r

In SCGLR: Supervised Component Generalized Linear Regression

# @title Projected Iterated Normed Gradient
# @author F. mortier, C. Trottier, G. Cornu and X. Bry
# @param Z matrix of the working variables
# @param X matrix of the normalized covariates
# @param AX matrix of the suplementary covariates
# @param W matrix of weights
# @param u vector of loadings
# @param method structural relevance criterion. Object of class "method.SCGLR"  built by  \code{\link{methodSR}} for Structural Relevance.
# @return u_new the updated loading vectors

ping <- function(Z, X, AX, W, F, u, method) {
  u_cur <- u
  h_cur <- hFunct(Z=Z, X=X, AX=AX, W=W, u=u_cur, method=method)$h
  u_new <- update_u(Z=Z, X=X, AX=AX, W=W, F=F, u=u_cur, method=method)
  h_new <- hFunct(Z=Z, X=X, AX=AX, W=W, u=u_new, method=method)$h
  ing_iter <- 0
  while((abs((h_new-h_cur)/h_cur)>method$epsilon) && (ing_iter<=method$maxiter)) {
    u_cur <- u_new
    h_cur <- h_new
    u_new <- update_u(Z=Z, X=X, AX=AX, W=W, F=F, u=u_cur, method=method)
    h_new <- hFunct(Z=Z, X=X, AX=AX, W=W, u=u_new, method=method)$h
    ing_iter <- ing_iter+1
  }
  return(u_new)
}

update_u <- function(Z, X, AX, W, F, u, method) {
  out_h <- hFunct(Z=Z, X=X, AX=AX, W=W, u=u, method=method)
  
  if(is.null(F)) {
    m <- out_h$gradh/sqrt(sum(out_h$gradh^2))
  } else {
    C <- (crossprod(X,F))#/nrow(X))
    proj_C_ortho <- out_h$gradh - C %*% solve(crossprod(C), crossprod(C, out_h$gradh))
    m <- c(proj_C_ortho / sqrt(sum(proj_C_ortho^2)))
  }
  
  h_m <- hFunct(Z=Z, X=X, AX=AX, W=W, u=m, method=method)$h
  h_u <- out_h$h
  k <- 1
  
  while((h_m<h_u) && (k<method$bailout)) {
    m <- c(u)+m
    m <- m/sqrt(sum(m^2))
    h_m <- hFunct(Z=Z, X=X, AX=AX, W=W, u=m, method=method)$h
    k <- k+1
  }
  
  if(k>method$bailout) print("ARGH") # TODO do something else :-)
  u_new <- m
  return(u_new)
}

hFunct<- function(Z, X, AX, W, u, method)
{
  f <- c(X %*% u)
  psi <- 0
  gradpsi <- rep(0, length(u))
  
  if(!is.null(AX)) {
    
    for(k in 1:ncol(Z)) {
      
      AXtWkAX <- crossprod(AX,W[,k]*AX)
      projWkfAX <- c(AX %*% solve(AXtWkAX, crossprod(AX,W[,k]*f)))
      projWkforthoAX <- f - projWkfAX
      Zk <- wtScale(Z[,k],W[,k])#Z[,k] - sum(W[,k]*Z[,k])
      WZk <- W[,k]*Zk
      projWkZAX <- AX %*% solve(AXtWkAX, crossprod(AX,WZk))
      
      # compute psi
      scalsqpfZ <- sum(c(projWkforthoAX)*WZk)^2
      scalsqpfpf <- sum(c(projWkforthoAX)^2*W[,k])
      term1psi <- sum(scalsqpfZ/(scalsqpfpf))
      term2psi <- sum(WZk*projWkZAX)
      psi <- psi+term1psi+term2psi
      
      # compute grad(psi)
      PiorthoPrimeWkZ <- WZk -  W[,k]*AX %*% solve(AXtWkAX,crossprod(AX,WZk))
      XprimeprojorthoWZ <- crossprod(X,PiorthoPrimeWkZ)
      term1 <- c(XprimeprojorthoWZ %*% crossprod(XprimeprojorthoWZ,u))/(scalsqpfpf)
      
      WprojWkOrthof <- W[,k]*projWkforthoAX
      term2 <-  scalsqpfZ*c(crossprod(X,WprojWkOrthof))/(scalsqpfpf^2)
      gradpsi <- gradpsi +(term1-term2)
    }
    
    gradpsi <- 2*gradpsi
    
  } else {
    
    for(k in 1:ncol(Z)){
      
      Zk <- wtScale(Z[,k],W[,k])#Z[,k] - sum(W[,k]*Z[,k])
      WZk <- W[,k]*Zk
      scalsqpfZ <- sum(c(f)*WZk)^2
      scalsqpfpf <- sum(c(f)^2*W[,k])
      
      # compute psi
      psi <- psi+sum(scalsqpfZ/(scalsqpfpf))
      
      # compute grad(psi)
      XprimeWZ <- crossprod(X,WZk) #X'W_k Z_k
      term1 <- c(XprimeWZ %*% crossprod(XprimeWZ,u))/(scalsqpfpf)
      term2 <-  scalsqpfZ*c(crossprod(X,W[,k]*f))/(scalsqpfpf^2)
      gradpsi <- gradpsi +(term1-term2)
    }
    
    gradpsi <- 2*gradpsi
  }
  
  n <- nrow(X)
  
  # calcul phi Component Variance: cv
  if(method$phi=="cv") {
    
    phi <- c(crossprod(f))/n
    
    # calcul grad phi
    gradphi <- c(2*crossprod(X,f/n))
    
  } else {
    
    ### autre calcul de phi avec l>=1 : vpi: Variable Powered Inertia
    scalsqfX <- colSums(f*X/n)
    XtWX <- crossprod(X)/n
    phi <- (sum((scalsqfX^2)^method$l))^(1/method$l)
    
    # calcul de grad phi
    gradphi <- 2*phi^(1-method$l)*rowSums(XtWX%*%diag(scalsqfX)^(2*method$l-1))
  }
  
  # calcul de h (s in R+)
  #h = log(psi)+method$s*log(phi)
  #gradh=gradpsi/psi+method$s*gradphi/phi
  
  # calcul de h (s in [0..1])
  h <- (1-method$s)*log(psi)+method$s*log(phi)
  gradh <- (1-method$s)*gradpsi/psi+method$s*gradphi/phi
  
  return(list(h=h, gradh=gradh,psi=psi,gradpsi=gradpsi,phi=phi,gradphi=gradphi))
}