R/EnvelopeModels.R

In abnormally-distributed/cvreg: Cross Validation and Robust Estimation Utilities

#' Envelope Regression for Sufficient Dimension Reduction
#'
#' @description This is an adaptation with some improvements on functions in the Renvlp package. This can fit
#' envelope models for dimension reduction to either univariate or multivariate numeric responses with numeric
#' predictors. In the case of a univaraite outcome, a predictor envelope model is fit. For multivariate outcomes,
#' an envelope is fit to both the response matrix and model matrix.
#'
#' @param formula model formula
#' @param data data frame
#' @param rank number of dimensions
#' @param yrank number of dimensions for a multivariate response.
#' @param se should standard errors for the regression basis be calculated? defaults to TRUE.
#' @param init an optional initial covariate dimension reduction subspace. defaults to NULL.
#'
#'
#' @return an sdr object
#' @export
#'
ENV <- function(formula, data, rank = "all", yrank = "all", se = TRUE, init = NULL){

  X = model.matrix(formula, data)[,-1]
  Y = model.frame(formula, data)
  Y = model.response(Y)
  Y <- as.matrix(Y)
  if (rank == "all" || is.null(rank)){
    rank = ncol(X)
  }
  if (ncol(Y)==1){
    ENV.uni(formula, data, rank, se, init)
  }
  else{
    if (yrank == "all" || yrank == ncol(Y)){ENV.mv(formula, data, rank, se, init)}
    else{ENV.mvxy(formula, data, rank, yrank, se, Pinit = init, Ginit = NULL)}
  }
}

ENV.uni <- function(formula, data, rank, se = TRUE, init = NULL) {

  X = model.matrix(formula, data)[,-1]
  Y = model.frame(formula, data)
  Y = model.response(Y)
  Y <- as.matrix(Y)
  X <- as.matrix(X)
  a <- dim(Y)
  n <- a[1]
  r <- a[2]
  p <- ncol(X)
  if (a[1] != nrow(X)) stop("X and Y should have the same number of observations.")
  if (rank > p || rank < 0) stop("rank must be an integer between 0 and p.")
  if(sum(duplicated(cbind(X, Y), MARGIN = 2)) > 0) stop("Some responses also appear in the predictors, or there maybe duplicated columns in X or Y.")

  SigmaY  <- cov(Y) * ((n-1)/n)
  SigmaYX <- cov(Y,X) * ((n-1)/n)
  SigmaX  <- solve(.ridgeinv(cov(X))) * ((n-1)/n)
  TauY <- .ridgeinv(SigmaY)
  eig.SigmaY <- eigen(SigmaY)
  U <- crossprod(SigmaYX, TauY) %*% SigmaYX
  M <- SigmaX - U
  tmp <- envlp_MU(M, U, rank)

  if (!is.null(init)) {
    if (nrow(init) != p || ncol(init) != rank) stop("The initial value should have p rows and rank columns.")
    tmp0 <- qr.Q(qr(init), complete = TRUE)
    tmp$Gammahat <- as.matrix(tmp0[, 1:rank])
    tmp$Gamma0hat <- as.matrix(tmp0[, (rank+1):p])
  }

  Gammahat <- tmp$Gammahat;Gamma0hat <- tmp$Gamma0hat;objfun <- tmp$objfun
  covMatrix <- NULL;std.err <- NULL;ratio <- NULL

  if (rank == p) {
    TauX <- .ridgeinv(SigmaX)
    olscoef <- tcrossprod(TauX, SigmaYX)
    etahat <- olscoef
    Omegahat <- SigmaX
    Omega0hat <- NULL
    muhat <- colMeans(Y) - crossprod(olscoef, colMeans(X))
    betahat <- olscoef
    SigmaXhat <- M + U
    SigmaYcXhat <- SigmaY - SigmaYX %*% olscoef
    log_lik <- - n * (r + p) / 2 * (log(2 * pi) + 1) - n / 2 * (objfun + sum(log(eig.SigmaY$values)))
    if (se == T) {
      covMatrix <- kronecker(SigmaYcXhat, TauX)
      std.err <- matrix(sqrt(diag(covMatrix)), nrow = p)
      ratio <- matrix(1, p, r)
    }
  } else {
    TauX <- .ridgeinv(SigmaX)
    etahat <- crossprod(Gammahat, t(SigmaYX))
    Omegahat <- crossprod(Gammahat, SigmaX) %*% Gammahat
    Omega0hat <- crossprod(Gamma0hat, SigmaX) %*% Gamma0hat
    invOmegahat <- .ridgeinv(Omegahat)
    betahat <- Gammahat %*% invOmegahat %*% etahat
    muhat <- colMeans(Y) - crossprod(betahat, colMeans(X))
    SigmaXhat <- Gammahat %*% tcrossprod(Omegahat, Gammahat) + Gamma0hat %*% tcrossprod(Omega0hat, Gamma0hat)
    olscoef <- tcrossprod(TauX, SigmaYX)
    PGamma <- tcrossprod(Gammahat)
    SigmaYcXhat <- SigmaY - SigmaYX %*% PGamma %*% .ridgeinv(SigmaXhat) %*% tcrossprod(PGamma , SigmaYX)
    log_lik <- - n * (r + p) / 2 * (log(2 * pi) + 1) - n / 2 * (objfun + sum(log(eig.SigmaY$values)))

    if (se == T) {
      covMatrix <- kronecker(SigmaYcXhat, TauX)
      seFm <- matrix(sqrt(diag(covMatrix)), nrow = p)
      TaumaYcXhat <- .ridgeinv(SigmaYcXhat)
      invOmega0hat <- .ridgeinv(Omega0hat)
      temp <- kronecker(etahat %*% tcrossprod(TaumaYcXhat, etahat), Omega0hat) + kronecker(invOmegahat, Omega0hat) + kronecker(Omegahat, invOmega0hat) - 2 * kronecker(diag(rank), diag(p - rank))
      temp2 <- kronecker(t(etahat), Gamma0hat)
      covMatrix <- kronecker(SigmaYcXhat, Gammahat %*% tcrossprod(invOmegahat, Gammahat)) + temp2 %*% tcrossprod(.ridgeinv(temp), temp2)
      std.err <- matrix(sqrt(diag(covMatrix)), nrow = p)
      ratio <- seFm / std.err
    }
  }

  betahat <- as.vector(betahat)
  names(betahat) <- colnames(X)
  predictors <- X %*% Gammahat
  fitted <- as.vector(X %*% betahat)
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  colnames(Gammahat) <- paste0("SP", 1:ncol(predictors))
  rownames(Gammahat) <- colnames(X)
  if (se){
    std.err <- std.err * (1/sqrt(n))
    p.values <- 2*pnorm(-abs((betahat/std.err)))
  } else {
    std.err <- NULL
    p.values <- NULL
  }
  return(structure(list(basis = Gammahat, coef = betahat, fitted = fitted, predictors = predictors, eta = etahat, muhat = muhat, Omega = Omegahat, SigmaYcX = SigmaYcXhat, SigmaX = SigmaXhat, log_lik = log_lik,  n = n, covMatrix = covMatrix, std.err = std.err, p.values = p.values, ratio = ratio,  y.obs = Y, formula = formula, mf = model.frame(formula, data)), class = "sdr", model = "ENV"))
}


ENV.mv <- function(formula, data, rank, se = TRUE, init = NULL) {

  X = model.matrix(formula, data)[,-1];Y = model.frame(formula, data)
  Y = model.response(Y);Y <- as.matrix(Y);X <- as.matrix(X)
  a <- dim(Y);n <- a[1];r <- a[2];p <- ncol(X);R <- rep(1, p)
  if (a[1] != nrow(X)) stop("X and Y should have the same number of observations.")
  if (rank > p || rank < 0) stop("rank must be an integer between 0 and p.")
  if(sum(duplicated(cbind(X, Y), MARGIN = 2)) > 0) stop("Some responses also appear in the predictors, or there maybe duplicated columns in X or Y.")
  SigmaY  <- cov(Y) * ((n-1)/n)
  SigmaYX <- cov(Y,X) * ((n-1)/n)
  SigmaX  <- solve(.ridgeinv(cov(X))) * ((n-1)/n)
  TauY <- .ridgeinv(SigmaY)
  eig.SigmaY <- eigen(SigmaY)
  eig.SigmaX <- eigen(SigmaX)
  U <- crossprod(SigmaYX, TauY) %*% SigmaYX
  M <- SigmaX - U
  q <- length(R)
  tmp <- envlp_MU(M, U, rank)

  if (!is.null(init)) {
    if (nrow(init) != p || ncol(init) != rank) stop("The initial value should have p rows and rank columns.")
    tmp0 <- qr.Q(qr(init), complete = TRUE)
    tmp$Gammahat <- as.matrix(tmp0[, 1:rank])
    tmp$Gamma0hat <- as.matrix(tmp0[, (rank+1):p])
  }

  Gammahat <- tmp$Gammahat
  Gamma0hat <- tmp$Gamma0hat
  objfun <- tmp$objfun
  covMatrix <- NULL
  std.err <- NULL
  ratio <- NULL

  if (rank >= (p - (q - 1) / r)) {
    Gammahat <- diag(p)
    Gamma0hat <- NULL
    TauX <- .ridgeinv(SigmaX)
    olscoef <- tcrossprod(TauX, SigmaYX)
    etahat <- olscoef
    Omegahat <- SigmaX
    Omega0hat <- NULL
    muYhat <- colMeans(Y)
    muXhat <- colMeans(X)
    betahat <- olscoef
    SigmaXhat <- SigmaX
    SigmaYcXhat <- SigmaY - SigmaYX %*% olscoef
    eig.SigmaXhat <- eigen(SigmaXhat)
    eig.SigmaYcXhat <- eigen(SigmaYcXhat)
    objfun <- (sum(log(eig.SigmaXhat$values)) + sum(log(eig.SigmaYcXhat$values)) + p)
    log_lik <- -n * (r + p)/2 * log(2 * pi) - n * r/2 - n/2 * objfun
    if (se == T) {
      covMatrix <- kronecker(SigmaYcXhat, TauX)
      std.err <- matrix(sqrt(diag(covMatrix)), nrow = p)
      ratio <- matrix(1, p, r)
    }
  } else {
    tmp1 <- sxenvlp_MU(X, Y, rank, R)

    if (!is.null(init)) {
      if (nrow(init) != p || ncol(init) != rank) stop("The initial value should have p rows and rank columns.")
      tmp0 <- qr.Q(qr(init), complete = TRUE)
      tmp1$Gammahat <- as.matrix(tmp0[, 1:rank])
      tmp1$Gamma0hat <- as.matrix(tmp0[, (rank+1):p])
      if (length(R) == p) {
        make_objfun <- function(SigmaY, SigmaYX, SigmaX){
          function(d, Gamma, X, Y){
            X <- as.matrix(X);Y <- as.matrix(Y);a <- dim(Y);n <- a[1];r <- a[2];p <- ncol(X)
            TauX <- .ridgeinv(SigmaX);TauY <- .ridgeinv(SigmaY)
            t1 <- crossprod(SigmaYX, TauY);U <- t1 %*% SigmaYX; M <- SigmaX - U
            d1 <- c(1, d);Lambda <- diag(d1);invLambda <- diag(1 / d1)
            m1 <- crossprod(Gamma, Lambda);eigtem1 <- eigen(m1 %*% tcrossprod(TauX, m1))
            m2 <- crossprod(Gamma, invLambda);eigtem2 <- eigen(m2 %*% tcrossprod(M, m2))
            temp1 <- sum(log(eigtem1$values));temp2 <- sum(log(eigtem2$values))
            objfun <- temp1 + temp2
            return(objfun)
          }
        }
        objfun <- make_objfun(SigmaY, SigmaYX, SigmaX)
        d.init <- rep(1, (p - 1))
        k2 <- rep(0, (p - 1))
        tmp.init <- Rsolnp::solnp(pars = d.init, fun = objfun, LB = k2,
                                  .control = list(delta = 1e-6, tol = 0.0005, trace = 0),
                                  Gamma = tmp1$Gammahat, X = X, Y = Y)
        d <- tmp.init$pars
        d1 <- c(1, d)
        tmp1$Lambdahat <- diag(d1)
      } else {
        make_objfun1 <- function(SigmaY, SigmaYX, SigmaX){
          function(d, Gamma, X, Y, R){
            X <- as.matrix(X);Y <- as.matrix(Y)
            a <- dim(Y);n <- a[1];r <- a[2];p <- ncol(X)
            TauX <- .ridgeinv(SigmaX);TauY<-.ridgeinv(SigmaY)
            t1 <- crossprod(SigmaYX, TauY)
            U <- t1 %*% SigmaYX;M <- SigmaX - U;Lambda <- diag(d);invLambda <- diag(1 / d)
            m1 <- crossprod(Gamma, Lambda);eigtem1 <- eigen(m1 %*% tcrossprod(TauX, m1))
            m2 <- crossprod(Gamma, invLambda);eigtem2 <- eigen(m2 %*% tcrossprod(M, m2))
            temp1 <- sum(log(eigtem1$values));temp2 <- sum(log(eigtem2$values));objfun <- temp1 + temp2
            return(objfun)
          }
        }
        objfun1 <- make_objfun1(SigmaY, SigmaYX, SigmaX)
        heq <- function (d, Gamma, X, Y, R){
          iter = sum(R)
          i <- 1
          cont <- NULL
          while (i < iter) {
            C <- matrix(0, (R[i] - 1), sum(R))
            for (j in 2 : (sum(R) - 1)){
              if (R[i] == j) {
                for (k in 1 : (j - 1)){
                  s <- sum(R[0 : (i - 1)]) + 1
                  C[k , s] <- 1
                  C[k , (s + k)] <- -1
                }
              }
            }
            cont <- rbind(cont, C)
            if (i == length(R)) {
              c1 <- c(1, rep(0, (sum(R) - 1)))
              cont <- rbind(c1, cont)
              break
            } else {
              i <- i + 1
            }
          }

          g <- rep(NA, nrow(cont))
          g[1] <- d[1]
          for (m in 2 : nrow(cont)){
            g[m] <- d[which(cont[m, ] == 1)] - d[which(cont[m, ] == -1)]
          }
          g
        }
        d.init <- rep(1, p)
        g1 <- heq(d.init, tmp1$Gammahat, X, Y, R)
        k1 <- c(1, rep(0, length(g1) - 1))
        k2 <- rep(0, p)

        tmp.init <- Rsolnp::solnp(pars = d.init, fun = objfun1, eqfun = heq, eqB = k1, LB = k2,
                                  .control = list(delta = 1e-6, tol = 0.0005, trace = 0),
                                  R = R, Gamma = tmp1$Gammahat, X = X, Y = Y)
        d <- tmp.init$pars
        tmp1$Lambdahat <- diag(d)
      }
    }
    TauX <- .ridgeinv(SigmaX)
    Gammahat <- tmp1$Gammahat
    Gamma0hat <- tmp1$Gamma0hat
    Lambdahat <- tmp1$Lambdahat
    d1 <- diag(Lambdahat)
    invLambdahat <- diag(1 / d1)
    E1 <- crossprod(Gammahat, invLambdahat)
    Omegahat <- E1 %*% tcrossprod(SigmaX, E1)
    E2 <- crossprod(Gamma0hat, invLambdahat)
    Omega0hat <- E2 %*% tcrossprod(SigmaX, E2)
    invOmegahat <- .ridgeinv(Omegahat)
    etahat <- invOmegahat %*% tcrossprod(E1, SigmaYX)
    betahat <- invLambdahat %*% Gammahat %*% etahat
    muYhat <- colMeans(Y)
    muXhat <- colMeans(X)
    E3 <- Lambdahat %*% Gammahat
    sig1 <- E3 %*% tcrossprod(Omegahat, E3)
    E4 <- Lambdahat %*% Gamma0hat
    sig2 <- E4 %*% tcrossprod(Omega0hat, E4)
    SigmaXhat <- sig1 + sig2
    Yc <- Y - tcrossprod(rep(1, nrow(Y)), colMeans(Y))
    Xc <- X - tcrossprod(rep(1, nrow(X)), colMeans(X))
    E5 <- Yc - Xc %*% betahat
    SigmaYcXhat <- crossprod(E5) / n
    TauYcX <- .ridgeinv(mpd(SigmaYcXhat))
    eig.SigmaXhat <- eigen(mpd(SigmaXhat))
    eig.SigmaYcXhat <- eigen(mpd(SigmaYcXhat))
    TauXhat <-.ridgeinv(mpd(SigmaXhat))
    s1 <- TauXhat %*% SigmaX
    eig.invX <- eigen(s1)
    objfun <- (sum(log(eig.SigmaXhat$values)) + sum(log(eig.SigmaYcXhat$values))  + sum(eig.invX$values))
    log_lik <- -n * (r + p)/2 * log(2 * pi) - n * r/2 - n/2 * objfun
    if (se == T) {
      U1 <- SigmaYX %*% TauX
      U2 <- tcrossprod(U1, SigmaYX)
      M2 <- SigmaY - U2
      covMatrix <- kronecker(M2, TauX)
      seFm <- matrix(sqrt(diag(covMatrix)), nrow = p)
      Ep <- .expan(p)
      Eu <- .expan(rank)
      Epu <- .expan((p - rank))
      Cp <- .contr(p)
      j11 <- kronecker(TauYcX, SigmaXhat)
      j222 <- kronecker(TauXhat, TauXhat)
      j22 <- crossprod(Ep, j222) %*% Ep / 2
      J <- cbind(rbind(j11, matrix(0, nrow = nrow(j22), ncol = ncol(j11))), rbind(matrix(0, nrow = nrow(j11), ncol = ncol(j22)), j22))
      diagvec <- function(R){
        E <- NULL
        p <- sum(R)
        for(i in 1 : (length(R) - 1)){
          e <- rep(0, p)
          s <- sum(R[1 : i])
          e[s + 1] <- 1
          t <- kronecker(e, e)
          E <- cbind(E, t)
        }
        E <- matrix(E, nrow = p^2)
      }
      L <- diagvec(R)
      SigmaL <- Gammahat %*% tcrossprod(Omegahat, Gammahat) + Gamma0hat %*% tcrossprod(Omega0hat, Gamma0hat)
      h1 <- crossprod(Gammahat, invLambdahat)
      h2 <- crossprod(etahat, h1)
      h11 <- - kronecker(h2, invLambdahat) %*% L
      h3 <- kronecker(Lambdahat %*% SigmaL, diag(p))
      h21 <- 2 * Cp %*% h3 %*% L
      h12 <- kronecker(diag(r), invLambdahat %*% Gammahat)
      h13 <- kronecker(t(etahat), invLambdahat)
      h4 <- Lambdahat %*% Gammahat %*% Omegahat
      h5 <- kronecker(h4, Lambdahat)
      h6 <- Lambdahat %*% Gammahat
      h7 <- Lambdahat %*% Gamma0hat
      h78 <- h7 %*% tcrossprod(Omega0hat, Gamma0hat)
      h8  <- kronecker(h6, h78)
      h23 <- 2 * Cp %*% (h5 - h8)
      h9 <- Lambdahat %*% Gammahat
      h24 <- Cp %*% kronecker(h9, h9) %*% Eu
      h10 <- Lambdahat %*% Gamma0hat
      h25 <- Cp %*% kronecker(h10, h10) %*% Epu
      H <- rbind(cbind(h11, h12, h13, matrix(0, nrow = nrow(h11), ncol = ncol(h24)), matrix(0, nrow = nrow(h11), ncol = ncol(h25))), cbind(h21, matrix(0, nrow = nrow(h21), ncol = ncol(h12)), h23, h24, h25))
      M1 <- crossprod(H, J) %*% H
      invM1 <- .ridgeinv(M1)
      V <- H %*% tcrossprod(invM1, H)
      covMatrix <- V[1:(p * r), 1:(p * r)]
      std.err <- matrix(sqrt(diag(covMatrix)), nrow = p)
      ratio <- seFm/std.err
    }
  }
  rownames(betahat) <- colnames(X)
  predictors <- X %*% Gammahat
  fitted <- X %*% betahat
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  colnames(Gammahat) <- paste0("SP", 1:ncol(predictors))
  rownames(Gammahat) <- colnames(X)
  if (se){
    std.err <- std.err * (1/sqrt(n))
    p.values <- 2*pnorm(-abs((betahat / std.err)))
  } else {
    std.err <- NULL
    p.values <- NULL
  }
  return(structure(list(basis = Gammahat, coef = betahat, fitted = fitted, predictors = predictors, eta = etahat, muhat = muYhat, Gamma0 = Gamma0hat, Omega = Omegahat, SigmaYcX = SigmaYcXhat, SigmaX = SigmaXhat, log_lik = log_lik,  n = n, covMatrix = covMatrix, std.err = std.err, p.values = p.values, ratio = ratio, formula = formula, mf = model.frame(formula, data), y.obs = Y), class = "sdr", model = "ENV"))
}


ENV.mvxy <- function(formula, data, xrank, yrank, se = TRUE, Pinit = NULL, Ginit = NULL) {

  X = model.matrix(formula, data)[,-1];Y = model.frame(formula, data)
  Y = model.response(Y);Y <- as.matrix(Y);X <- as.matrix(X)
  a <- dim(Y);n <- a[1];r <- a[2];p <- ncol(X)
  if (a[1] != nrow(X))
    stop("X and Y should have the same number of observations.")
  if (yrank > r || yrank < 0)
    stop("yrank must be an integer between 0 and r.")
  if (xrank > p || xrank < 0)
    stop("xrank must be an integer between 0 and p.")
  if (sum(duplicated(cbind(X, Y), MARGIN = 2)) > 0)
    stop("Some responses also appear in the predictors, or there maybe duplicated columns in X or Y.")

  SigmaY  <- cov(Y) * ((n-1)/n)
  SigmaYX <- cov(Y,X) * ((n-1)/n)
  SigmaX  <- solve(.ridgeinv(cov(X))) * ((n-1)/n)
  TauX <- .ridgeinv(SigmaX)
  betaOLS <- SigmaYX %*% TauX
  covMatrix <- NULL
  std.err <- NULL
  if (yrank == 0 || xrank ==0) {
    Gammahat <- NULL
    Gamma0hat <- diag(r)
    Phihat <- NULL
    Phi0hat <- diag(p)
    etahat <- NULL
    Omegahat <- NULL
    Omega0hat <- SigmaY
    Deltahat <- NULL
    Delta0hat <- SigmaX
    muhat <- colMeans(Y)
    betahat <- matrix(0, p, r)
    SigmaYcXhat <- SigmaY
    SigmaXhat <- SigmaX
    tmp1 <- eigen(SigmaX)
    tmp2 <- eigen(SigmaY)
    objfun <- sum(log(tmp1$values)) + sum(log(tmp2$values))
    log_lik <- - n/2 * objfun
  }
  else if (yrank == r && xrank == p) {
    M <- SigmaY - tcrossprod(betaOLS, SigmaYX)
    Gammahat <- diag(r)
    Gamma0hat <- NULL
    Phihat <- diag(p)
    Phi0hat <- NULL
    etahat <- betaOLS
    Omegahat <- M
    Omega0hat <- NULL
    Deltahat <- SigmaX
    Delta0hat <- NULL
    muhat <- colMeans(Y) - betaOLS %*% colMeans(X)
    betahat <- t(betaOLS)
    SigmaYcXhat <- M
    SigmaXhat <- SigmaX
    tmp.e1 <- eigen(Omegahat)
    tmp.e3 <- eigen(Deltahat)
    objfun <- sum(log(tmp.e1$values)) + sum(log(tmp.e3$values))
    log_lik <- -n * (p + r) / 2 - n/2 * objfun
    if (se == T) {
      covMatrix <- kronecker(M, TauX)
      std.err <- matrix(sqrt(diag(covMatrix)), nrow = p)
    }
  }
  else if (xrank == p) {

    tmp.envfun <- function(SigmaX, SigmaY, SigmaYX, TauX, yrank){
      betaOLS <- SigmaYX %*% TauX
      U <- tcrossprod(betaOLS, SigmaYX)
      M <- SigmaY - U
      tmp <- envlp_MU(M, U, yrank)
      Gammahat <- tmp$Gammahat
      Gamma0hat <- tmp$Gamma0hat
      etahat <- crossprod(Gammahat, betaOLS)
      betahat <- Gammahat %*% etahat
      muhat <- colMeans(Y) - betahat %*% colMeans(X)
      Omegahat <- crossprod(Gammahat, M) %*% Gammahat
      Omega0hat <- crossprod(Gamma0hat, sigY) %*% Gamma0hat
      Sigma1 <- Gammahat %*% tcrossprod(Omegahat, Gammahat)
      Sigmahat <- Sigma1 + Gamma0hat %*% tcrossprod(Omega0hat, Gamma0hat)
      return(list(Gamma = Gammahat, Gamma0 = Gamma0hat, Omega = Omegahat, Omega0 = Omega0hat, beta = betahat, Sigma = Sigmahat))
    }

    tmp.env <- tmp.envfun(SigmaX, SigmaY, SigmaYX, TauX, yrank)
    Gammahat <- tmp.env$Gamma
    Gamma0hat <- tmp.env$Gamma0
    Phihat <- diag(p)
    Phi0hat <- NULL
    etahat <- TauX %*% crossprod(SigmaYX, Gammahat)
    Omegahat <- tmp.env$Omega
    Omega0hat <- tmp.env$Omega0
    Deltahat <- SigmaX
    Delta0hat <- NULL
    betahat <- t(tmp.env$beta)
    muhat <- colMeans(Y) - crossprod(betahat, colMeans(X))
    SigmaYcXhat <- tmp.env$Sigma
    SigmaXhat <- SigmaX
    tmp.e1 <- eigen(Omegahat)
    tmp.e2 <- eigen(Omega0hat)
    tmp.e3 <- eigen(Deltahat)
    invome <- .ridgeinv(Omegahat)
    invome0 <- .ridgeinv(Omega0hat)
    invdel <- .ridgeinv(Deltahat)
    t3 <- crossprod(Gammahat, SigmaY) %*% Gammahat %*% invome
    t5 <- crossprod(etahat, SigmaX) %*% etahat %*% invome
    t6 <- t(etahat) %*% crossprod(SigmaYX, Gammahat) %*% invome
    obj1 <- sum(log(tmp.e1$values)) + sum(log(tmp.e2$values))
    obj2 <- sum(log(tmp.e3$values)) + p
    obj3 <- sum(diag(t3)) + (r - yrank)
    obj4 <- sum(diag(t5)) - 2 * sum(diag(t6))
    objfun <- obj1 + obj2 + obj3 + obj4
    log_lik <- - (n/2) * objfun
    if (se == T) {
      M <- SigmaY - tcrossprod(betaOLS, SigmaYX)
      covMatrix <- kronecker(M, TauX)
      seFm <- matrix(sqrt(diag(covMatrix)), nrow = p)
      m1 <- Gammahat %*% tcrossprod(Omegahat, Gammahat)
      temp1 <- kronecker(m1, TauX)
      m2 <- crossprod(etahat, SigmaX) %*% etahat
      m3 <- kronecker(invome0, m2) + kronecker(invome0, Omegahat)
      m4 <- kronecker(Omega0hat, invome) - 2 * kronecker(diag(r - yrank), diag(yrank))
      m5 <- m3 + m4
      invm5 <- .ridgeinv(m5)
      m6 <- kronecker(Gamma0hat, etahat)
      temp2 <- m6 %*% tcrossprod(invm5, m6)
      covMatrix <- temp1 + temp2
      std.err <- matrix(sqrt(diag(covMatrix)), nrow = p)
    }
  }
  else{
    tmp.stenv <- xyenvlp_MU(X, Y, xrank, yrank)
    if (!is.null(Ginit)) {
      if (nrow(Ginit) != r || ncol(Ginit) != yrank) stop("The initial value should have r rows and yrank columns.")
      tmp0 <- qr.Q(qr(Ginit), complete = TRUE)
      tmp.stenv$Gammahat <- as.matrix(tmp0[, 1:yrank])
      tmp.stenv$Gamma0hat <- as.matrix(tmp0[, (yrank+1):r])
    }

    if (!is.null(Pinit)) {
      if (nrow(Pinit) != p || ncol(Pinit) != xrank) stop("The initial value should have p rows and xrank columns.")
      tmp0 <- qr.Q(qr(Pinit), complete = TRUE)
      tmp.stenv$Phihat <- as.matrix(tmp0[, 1:xrank])
      tmp.stenv$Phi0hat <- as.matrix(tmp0[, (xrank+1):p])
    }

    Gammahat <- tmp.stenv$Gammahat
    Gamma0hat <- tmp.stenv$Gamma0hat
    Phihat <- tmp.stenv$Phihat
    Phi0hat <- tmp.stenv$Phi0hat
    Deltahat <- crossprod(Phihat, SigmaX) %*% Phihat
    Delta0hat <- crossprod(Phi0hat, SigmaX) %*% Phi0hat
    invdel <- .ridgeinv(Deltahat)
    et1 <- t(Phihat) %*% crossprod(SigmaYX, Gammahat)
    etahat <- invdel %*% et1
    om1 <- SigmaYX %*% Phihat
    om2 <- om1 %*% tcrossprod(invdel, om1)
    om3 <- SigmaY - om2
    Omegahat <- crossprod(Gammahat, om3) %*% Gammahat
    Omega0hat <- crossprod(Gamma0hat, SigmaY) %*% Gamma0hat
    betahat <- Phihat %*% tcrossprod(etahat, Gammahat)
    muhat <- colMeans(Y) - crossprod(betahat, colMeans(X))
    s1 <- Gammahat %*% tcrossprod(Omegahat, Gammahat)
    SigmaYcXhat <- s1 + Gamma0hat %*% tcrossprod(Omega0hat, Gamma0hat)
    s2 <- Phihat %*% tcrossprod(Deltahat, Phihat)
    SigmaXhat <- s2 + Phi0hat %*% tcrossprod(Delta0hat, Phi0hat)
    tmp.e1 <- eigen(Omegahat)
    tmp.e2 <- eigen(Omega0hat)
    tmp.e3 <- eigen(Deltahat)
    tmp.e4 <- eigen(Delta0hat)
    invome <- .ridgeinv(Omegahat)
    invome0 <- .ridgeinv(Omega0hat)
    invdel0 <- .ridgeinv(Delta0hat)
    t3 <- crossprod(Gammahat, SigmaY) %*% Gammahat %*% invome
    peta <- Phihat %*% etahat
    t5 <- crossprod(peta, SigmaX) %*% peta %*% invome
    t6 <- t(peta) %*% crossprod(SigmaYX, Gammahat) %*% invome
    obj1 <- sum(log(tmp.e1$values)) + sum(log(tmp.e2$values))
    obj2 <- sum(log(tmp.e3$values)) + sum(log(tmp.e4$values)) + p
    obj3 <- sum(diag(t3)) + (r - yrank)
    obj4 <- sum(diag(t5)) - 2 * sum(diag(t6))
    objfun <- obj1 + obj2 + obj3 + obj4
    log_lik <- - (n/2) * objfun
    if (se == T) {
      M <- SigmaY - tcrossprod(betaOLS, SigmaYX)
      covMatrix <- kronecker(M, TauX)
      seFm <- matrix(sqrt(diag(covMatrix)), nrow = p)
      m1 <- Gammahat %*% tcrossprod(Omegahat, Gammahat)
      m2 <- Phihat %*% tcrossprod(invdel, Phihat)
      temp1 <- kronecker(m1, m2)
      d1 <- etahat %*% tcrossprod(invome, etahat)
      m3 <- kronecker(d1, Delta0hat)
      m4 <- kronecker(Deltahat, invdel0) + kronecker(invdel,Delta0hat) - 2 * kronecker(diag(xrank), diag(p - xrank))
      m5 <- m3 + m4
      invm5 <- .ridgeinv(m5)
      m6 <- kronecker(tcrossprod(Gammahat, etahat), Phi0hat)
      temp2 <- m6 %*% tcrossprod(invm5, m6)
      d2 <- crossprod(etahat, Deltahat) %*% etahat
      m8 <- kronecker(invome0, d2) + kronecker(invome0, Omegahat)
      m9 <- kronecker(Omega0hat, invome) - 2 * kronecker(diag(r - yrank), diag(yrank))
      m10 <- m8 + m9
      invm10 <- .ridgeinv(m10)
      m11 <- kronecker(Gamma0hat, Phihat %*% etahat)
      temp3 <- m11 %*% tcrossprod(invm10, m11)
      covMatrix <- temp1 + temp2 + temp3
      std.err <- matrix(sqrt(diag(covMatrix)), nrow = p)
    }
  }

  colnames(betahat) <- colnames(Y)
  rownames(betahat) <- colnames(X)
  predictors <- X %*% Phihat
  y.reduced <- Y %*% Gammahat
  fitted <- X %*% betahat
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  colnames(Phihat) <- paste0("SP", 1:ncol(predictors))
  rownames(Phihat) <- colnames(X)
  if (se){
    std.err <- std.err * (1/sqrt(n))
    p.values <- 2*pnorm(-abs((betahat / std.err)))
  } else {
    std.err <- NULL
    p.values <- NULL
  }

  return(structure(list(coef = betahat, SigmaYcX = SigmaYcXhat, SigmaX = SigmaXhat, Gammahat = Gammahat, Gamma0 = Gamma0hat,
              eta = etahat, Omega = Omegahat, Omega0 = Omega0hat, muhat = muhat, basis = Phihat, Phi0 = Phi0hat,
              Delta = Deltahat, Delta0 = Delta0hat, log_lik = log_lik,  covMatrix = covMatrix, std.err = std.err,
              p.values = p.values, n = n, y.obs = Y, fitted = fitted, y.reduced = y.reduced, predictors = predictors), class = "sdr", model = "ENV"))
}



#' Envelope Estimation of Multivariate Means with Between-Level Heteroscedastic Error
#'
#' @description This model is for fitting a multivariate mean when there is a single categorical predictor. This is
#' similar to MANOVA.
#'
#'
#' @param formula model formula. the response must be multivariate numeric data, and the predictor a single factor.
#' @param data data frame
#' @param rank number of dimensions
#' @param se should standard errors for the regression basis be calculated? defaults to TRUE.
#' @param init an optional initial dimension reduction subspace. defaults to NULL.
#'
#' @return an sdr object
#' @export
#'
#' @references
#' Su, Z. and Cook, R. D. (2013)  Estimation of Multivariate Means with Heteroscedastic Error Using Envelope Models. Statistica Sinica, 23, 213-230.
#'
HENV <- function(formula, data, rank = "all", se = TRUE, init = NULL) {
  fit = TRUE
  levelnames = as.character(levels(model.frame(formula, data)[,2]))
  data = as.data.frame(lapply(data, as.numeric))
  X = model.matrix(formula, data)[,-1]
  if (rank == "all" || is.null(rank)){
    rank = ncol(X)
  }
  Y = model.frame(formula, data)
  Y = model.response(Y)
  Y <- as.matrix(Y)
  groupind <- unique(X)
  XX <- as.factor(X)
  X <- match(XX, levels(XX))
  XX <- as.factor(X)
  Y <- as.matrix(Y)
  a <- dim(Y)
  n <- a[1]
  r <- a[2]
  p <- nlevels(XX)
  u <- r-1
  if(rank < 0 | rank > r){
    stop("rank should be an interger between 0 and r")
  }

   ncumx <- c()
  for (i in 1 : p) {
    ncumx[i] <- length(which(XX == as.numeric(levels(XX)[i])))
  }
  ncum <- cumsum(ncumx)
  ng <- diff(c(0, ncum))
  fracN <- ng / n
  sortx <- sort(X, index.return = T)
  Xs <- sortx$x
  ind <- sortx$ix
  Ys <- Y[ind, ]
  mY <- colMeans(Y)
  sigres = invres <- list(length = p)
  mYg <- matrix(rep(0, r*p), ncol = p)
  for (i in 1 : p) {
    if(i > 1) {
      yy <- Ys[(ncum[i-1] + 1):ncum[i], ]
      sigres[[i]] <- cov(yy) * (ng[i] - 1)/ng[i]
      mYg[ , i] <- colMeans(yy)
      invres[[i]] <- .ridgeinv(sigres[[i]])
    } else {
      yy <- Ys[1:ncum[i], ]
      sigres[[i]]<- cov(yy) * (ng[i] - 1)/ng[i]
      mYg[ , i] <- colMeans(yy)
      invres[[i]] <- .ridgeinv(sigres[[i]])
    }
  }
  SigmaY <- solve(.ridgeinv(cov(Y)))  * (n - 1)/n
  eigtemY <- eigen(SigmaY)$values
  logDetSigmaY <- log(prod(eigtemY[eigtemY > 0]))
  invSigmaY <- .ridgeinv(SigmaY)
  U = M <- list(length = p)
  for (i in 1 : p){
    M[[i]] <- sigres[[i]] * (ng[i] - 1)/ng[i]
    U[[i]] <- SigmaY - M[[i]]
  }
  MU <- SigmaY
  tmp <- henvlp_MU(M, U, MU, rank, n, ng, p)
  if (!is.null(init)) {
    if (nrow(init) != r || ncol(init) != rank) stop("The initial value should have r rows and rank columns.")
    tmp0 <- qr.Q(qr(init), complete = TRUE)
    tmp$Gammahat <- as.matrix(tmp0[, 1:rank])
    tmp$Gamma0hat <- as.matrix(tmp0[, (rank+1):r])
  }
  Gammahat <- tmp$Gammahat
  Gamma0hat <- tmp$Gamma0hat
  covMatrix <- NULL
  std.err <- NULL
  ratio <- NULL
  Yfit <- NULL
  if (rank == r) {
    Sigmahat <- sigres
    muhat <- as.matrix(mY)
    mughat <- mYg
    betahat <- mughat - mY %*% matrix(1, 1, p)
    etahat <- betahat
    Omegahat <- sigres
    Omega0hat <- NULL
    log_lik <- -n * r * (1 + log(2 * pi)) / 2
    for (i in 1 : p) {
      eig <- eigen(sigres[[i]])
      log_lik <- log_lik - ng[i] * sum(log(eig$values)) / 2
    }
    if (se == T) {
      J <- matrix(0, p * r + p * r * (r + 1) / 2, p * r + p * r * (r + 1) / 2)
      for (i in 1 : (p - 1)) {
        for (j in 1 : (p - 1)) {
          aa <- (i - 1) * r + 1
          bb <- i * r
          cc <- (j - 1) * r + 1
          dd <- j * r
          J[aa : bb, cc : dd] <- (fracN[j] * fracN[i] / fracN[p] * diag(r) %*% invres[[p]])
        }
        ee <- (i - 1) * r + 1
        ff <- i * r
        J[ee : ff, ee : ff] <- (fracN[i] * diag(r) %*% invres[[i]] + fracN[i] ^ 2 / fracN[p] * diag(r) %*% invres[[p]])
      }
      for (i in 1 : p) {
        aaa <- r * (p - 1) + 1 + (i - 1) * r * (r + 1) / 2
        bbb <- r * (p - 1) + i * r * (r + 1) / 2
        J[aaa : bbb, aaa : bbb] <- (0.5 * fracN[i] * crossprod(.expan(r), kronecker(invres[[i]], invres[[i]])) %*% .expan(r))
      }
      ggg <- r + 1
      hhh <- p * r + p * r * (r + 1) / 2
      iii <- (p - 1) * r + p * r * (r + 1) / 2
      J[ggg : hhh, ggg : hhh] <- J[1 : iii, 1 : iii]
      J[1 : r, ] <- 0
      J[ggg : hhh, 1 : r] <- 0
      for (i in 1 : p) {
        J[1 : r, 1 : r] <- J[1 : r, 1 : r] + fracN[i] * diag(r) %*% invres[[i]]
      }
      for (i in 1 : (p - 1)) {
        kkk <- i * r + 1
        lll <- (i + 1) * r
        J[1 : r, kkk : lll] <- fracN[i] * (invres[[p]] - invres[[i]])
        J[kkk : lll, 1 : r] <- J[1 : r, kkk : lll]
      }
      temp <- .ridgeinv(J)
      temp1 <- kronecker(matrix(1, 1, (p - 1)), diag(r))
      ccc <- r + 1
      ddd <- r * p
      vargroup <- temp1 %*% tcrossprod(temp[ccc : ddd, ccc : ddd], temp1)
      covMatrix <- matrix(0, r * (p + 1), r * (p + 1))
      covMatrix[1 : ddd, 1 : ddd] <- temp[1 : ddd, 1 : ddd]
      eee <- r * p + 1
      fff <- r * (p + 1)
      covMatrix[eee : fff, eee : fff] <- vargroup
      for (i in 1 : (p - 1)) {
        mmm <- i * r + 1
        nnn <- (i + 1) * r
        covMatrix[eee : fff, mmm : nnn] <- (- temp1 %*% temp[ccc : ddd, mmm : nnn])
      }
      covMatrix[ccc : ddd, eee : fff] <- t(covMatrix[eee : fff, ccc : ddd])
      covMatrix[eee : fff, 1 : r] <- (- temp1 %*% temp[ccc : ddd, 1 : r])
      covMatrix[1 : r, eee : fff] <- t(covMatrix[eee : fff, 1 : r])
      seFm <- matrix(sqrt(diag(covMatrix[ccc : fff, ccc : fff])), r, p)
      std.err <- seFm
      ratio <- matrix(1, r, p)
    }
    if (fit == T) {
      Yfit <- matrix(0, n, r)
      for (i in 1:p) {
        if (i > 1) {
          Yfit[ind[(ncum[i - 1] + 1):ncum[i]], ] <- tcrossprod(matrix(1, ng[i], 1) , mughat[ , i])
        }
        else {
          Yfit[ind[1:ncum[1]], ] <- tcrossprod(matrix(1, ng[1], 1) , mughat[ , 1])
        }
      }
    }
  }
  else {
    muhat <- as.matrix(mY)
    Omega0hat <- crossprod(Gamma0hat, SigmaY) %*% Gamma0hat
    etahat = betahat = etm <- list(length = p)
    for (i in 1:p) {
      etm[[i]] <- mYg[ , i] - mY
      etahat[[i]] <- crossprod(Gammahat, etm[[i]])
      betahat[[i]] <- Gammahat %*% etahat[[i]]
    }
    Sigmahat = Omegahat = mughat = invsig <- list(length = p)
    for (i in 1:p) {
      Omegahat[[i]] <- crossprod(Gammahat, sigres[[i]]) %*% Gammahat
      Sigmahat[[i]] <- Gammahat %*% tcrossprod(Omegahat[[i]], Gammahat) + Gamma0hat %*% tcrossprod(Omega0hat, Gamma0hat)
      mughat[[i]] <- muhat + betahat[[i]]
      invsig[[i]] <- .ridgeinv(Sigmahat[[i]])
    }
    beta <- c()
    for (i in 1 : p) {
      beta <- cbind(beta, betahat[[i]])
    }
    betahat <- beta
    mug <- c()
    for (i in 1 : p) {
      mug <- cbind(mug, mughat[[i]])
    }
    mughat <- mug
    e11 <- crossprod(Gammahat, invSigmaY) %*% Gammahat
    eig2 <- eigen(e11)
    r1 <- sum(log(eig2$values))
    log_lik <- - n * r * (1 + log(2 * pi)) /2 - n * logDetSigmaY / 2 - n * r1 / 2
    for (i in 1:p) {
      e22 <- crossprod(Gammahat, sigres[[i]]) %*% Gammahat
      eig <- eigen(e22)
      log_lik <- log_lik - ng[i] * sum(log(eig$values))/2
    }
    if (se == T) {
      J <- matrix(0, p * r + p * r * (r + 1) / 2, p * r + p * r * (r + 1) / 2)
      for (i in 1 : (p - 1)) {
        for (j in 1 : (p - 1)) {
          aa <- (i - 1) * r + 1
          bb <- i * r
          cc <- (j - 1) * r + 1
          dd <- j * r
          J[aa : bb, cc : dd] <- (fracN[j] * fracN[i] / fracN[p] * diag(r) %*% invsig[[p]])
        }
        ee <- (i - 1) * r + 1
        ff <- i * r
        J[ee : ff, ee : ff] <- (fracN[i] * diag(r) %*% invsig[[i]] + fracN[i] ^ 2 / fracN[p] * diag(r) %*% invsig[[p]])
      }
      for (i in 1 : p) {
        aaa <- r * (p - 1) + 1 + (i - 1) * r * (r + 1) / 2
        bbb <- r * (p - 1) + i * r * (r + 1) / 2
        J[aaa : bbb, aaa : bbb] <- (0.5 * fracN[i] * crossprod(.expan(r), kronecker(invsig[[i]], invsig[[i]])) %*% .expan(r))
      }
      ggg <- r + 1
      hhh <- p * r + p * r * (r + 1) / 2
      iii <- (p - 1) * r + p * r * (r + 1) / 2
      J[ggg : hhh, ggg : hhh] <- J[1 : iii, 1 : iii]
      J[1 : r, ] <- 0
      J[ggg : hhh, 1 : r] <- 0
      for (i in 1 : p) {
        J[1 : r, 1 : r] <- J[1 : r, 1 : r] + fracN[i] * diag(r) %*% invsig[[i]]
      }
      for (i in 1 : (p - 1)) {
        kkk <- i * r + 1
        lll <- (i + 1) * r
        J[1 : r, kkk : lll] <- fracN[i] * (invsig[[p]] - invsig[[i]])
        J[kkk : lll, 1 : r] <- J[1 : r, kkk : lll]
      }

      J1 <- matrix(0, p * r + p * r * (r + 1) / 2, p * r + p * r * (r + 1) / 2)
      for (i in 1 : (p - 1)) {
        for (j in 1 : (p - 1)) {
          aa <- (i - 1) * r + 1
          bb <- i * r
          cc <- (j - 1) * r + 1
          dd <- j * r
          J1[aa : bb, cc : dd] <- (fracN[j] * fracN[i] / fracN[p] * diag(r) %*% invres[[p]])
        }
        ee <- (i - 1) * r + 1
        ff <- i * r
        J1[ee : ff, ee : ff] <- (fracN[i] * diag(r) %*% invres[[i]] + fracN[i] ^ 2 / fracN[p] * diag(r) %*% invres[[p]])
      }
      for (i in 1 : p) {
        aaa <- r * (p - 1) + 1 + (i - 1) * r * (r + 1) / 2
        bbb <- r * (p - 1) + i * r * (r + 1) / 2
        J1[aaa : bbb, aaa : bbb] <- (0.5 * fracN[i] * crossprod(.expan(r), kronecker(invres[[i]], invres[[i]])) %*% .expan(r))
      }
      J1[ggg : hhh, ggg : hhh] <- J[1 : iii, 1 : iii]
      J1[1 : r, ] <- 0
      J1[ggg : hhh, 1 : r] <- 0
      for (i in 1 : p) {
        J1[1 : r, 1 : r] <- J[1 : r, 1 : r] + fracN[i] * diag(r) %*% invres[[i]]
      }
      for (i in 1 : (p - 1)) {
        kkk <- i * r + 1
        lll <- (i + 1) * r
        J1[1 : r, kkk : lll] <- fracN[i] * (invres[[p]] - invres[[i]])
        J1[kkk : lll, 1 : r] <- J1[1 : r, kkk : lll]
      }
      temp1 <- .ridgeinv(J1)
      seFm <- matrix(sqrt(diag(temp1[1 : r * p, 1 : r * p])), r, p)
      aaaa <- p * r + p * r * (r + 1) / 2
      bbbb <- r + rank * (r + p - 1 - rank) + p * rank * (rank + 1) / 2 + (r - rank) * (r - rank + 1) / 2
      H <- matrix(0, aaaa, bbbb)
      for (i in 1 : (p - 1)) {
        cccc <- (i - 1) * r + 1
        dddd <- i * r
        eeee <- (i - 1) * rank + 1
        ffff <- i * rank
        gggg <- (p - 1) * rank + 1
        hhhh <- rank * (r + p - 1 - rank)
        H[cccc : dddd, eeee : ffff] <- Gammahat
        H[cccc : dddd, gggg : hhhh] <- kronecker(t(etahat[[i]]), Gamma0hat)
      }
      for (i in 1 : p) {
        iiii <- r * (p - 1) + (i - 1) * r * (r + 1) / 2 + 1
        jjjj <- r * (p - 1) + i * r * (r + 1) / 2
        kkkk <- (p - 1) * rank + 1
        llll <- rank * (r + p - 1 - rank)
        mmmm <- rank * (r + p - 1 - rank) + (i - 1) * rank * (rank + 1) / 2 + 1
        nnnn <- rank * (r + p - 1 - rank) + i * rank * (rank + 1) / 2
        oooo <- rank * (r + p - 1 - rank) + p * rank * (rank + 1) / 2 + 1
        pppp <- rank * (r + p - 1 - rank) + p * rank * (rank + 1) / 2 + (r - rank) * (r - rank + 1) / 2
        H[iiii : jjjj, kkkk : llll] <- 2 * .contr(r) %*% (kronecker(Gammahat %*% Omegahat[[i]], Gamma0hat) -kronecker(Gammahat, Gamma0hat %*% Omega0hat))
        H[iiii : jjjj, mmmm : nnnn] <- .contr(r) %*% kronecker(Gammahat, Gammahat) %*% .expan(rank)
        H[iiii : jjjj, oooo : pppp] <- .contr(r) %*% kronecker(Gamma0hat, Gamma0hat) %*% .expan(r - rank)
      }
      qqqq <- r + 1
      rrrr <- p * r + p * r * (r + 1) / 2
      ssss <- r + rank * (r + p - 1 - rank) + p * rank * (rank + 1) / 2  + (r - rank) * (r - rank + 1) / 2
      tttt <- (p - 1) * r + p * r * (r + 1) / 2
      uuuu <- rank * (r + p - 1 - rank) + p * rank * (rank + 1) / 2 + (r - rank) * (r - rank + 1) / 2
      wwww <- r * p
      xxxx <- r * p + 1
      yyyy <- r * (p + 1)
      H[qqqq : rrrr, qqqq : ssss] <- H[1 : tttt, 1 : uuuu]
      H[1 : r, qqqq : ssss] <- 0
      H[qqqq : rrrr, 1 : r] <- 0
      H[1 : r, 1 : r] <- diag(r)
      temp3 <- pseudoinverse(.ridgeinv(crossprod(H, J) %*% H),tol=1e-100)
      invtemp3 <- .ridgeinv(temp3)
      temp <- H %*% tcrossprod(invtemp3, H)
      temp2 <- kronecker(matrix(1, 1, (p - 1)), diag(r))
      vargroup <- temp2 %*% tcrossprod(temp[qqqq : wwww, qqqq : wwww], temp2)
      covMatrix <- matrix(0, yyyy, yyyy)
      covMatrix[1 : wwww, 1 : wwww] <- temp[1 : wwww, 1 : wwww]
      covMatrix[xxxx : yyyy, xxxx : yyyy] <- vargroup
      for (i in 1 : (p - 1)) {
        zzzz <- i * r + 1
        zzzzz <- (i + 1) * r
        covMatrix[xxxx : yyyy, zzzz : zzzzz] <- (- temp2 %*% temp[qqqq : wwww, zzzz : zzzzz])
      }
      covMatrix[qqqq : wwww, xxxx : yyyy] <- t(covMatrix[xxxx : yyyy, qqqq : wwww])
      covMatrix[xxxx : yyyy, 1 : r] <- (- temp2 %*% temp[qqqq : wwww, 1 : r])
      covMatrix[1 : r, xxxx : yyyy] <- t(covMatrix[xxxx : yyyy, 1 : r])
      std.err <- matrix(sqrt(diag(covMatrix[qqqq : yyyy, qqqq : yyyy])), r, p)
      ratio <- seFm / std.err
    }
    if (fit == T) {
      Yfit <- matrix(rep(0, n * r), ncol = r)
      for (i in 1:p) {
        if (i > 1) {
          Yfit[ind[(ncum[i - 1] + 1):ncum[i]], ] <- tcrossprod(matrix(1, ng[i], 1), mughat[ , i])
        }
        else {
          Yfit[ind[1:ncum[1]], ] <- tcrossprod(matrix(1, ng[1], 1), mughat[ ,1])
        }
      }
    }
  }
  colnames(Gammahat) <- paste0("SP", 1:ncol(Gammahat))
  predictors <- Y  %*% Gammahat
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  means <- as.data.frame(mughat)
  colnames(means) <- levelnames
  rownames(Gammahat) <- rownames(means)
  colnames(betahat) <- colnames(means)
  rownames(betahat) <- rownames(means)
  if (se){
    std.err <- std.err * (1/sqrt(n))
    p.values <- 2*pnorm(-abs((betahat / std.err)))
  } else {
    std.err <- NULL
    p.values <- NULL
  }
  return(structure(list(basis = Gammahat, coef = betahat, Gamma0 = Gamma0hat,
                        Sigma = Sigmahat, eta = etahat, Omega = Omegahat,
                        Omega0 = Omega0hat, mu = muhat, means = means, log_lik = log_lik,
                        n = n, ng = ng, fitted = Yfit, predictors = predictors, covMatrix = covMatrix,
                        std.err = std.err, p.values = p.values, ratio = ratio, groupInd = groupind), class = "sdr", model = "HENV"))
}




#' Generalized Linear Envelopes for Sufficient Dimension Reduction
#'
#' @description This is an adaptation with several improvements on functions in the Renvlp package. Incorporated
#' here are methods for fitting predictor envelopes for binomial and poisson regression with Firth's biased reduced
#' estimation process, which removes the second-order bias in the maximum likelihood estimate. This also improves
#' the resilience of the estimator when there is class separation, in which case the MLE is degenerate.
#' \cr
#' @param formula model formula
#' @param data data frame
#' @param rank number of dimensions
#' @param se should standard errors for the regression basis be calculated? defaults to TRUE.
#' @param init an optional initial covariate dimension reduction subspace. defaults to NULL.
#'
#' @return an sdr object
#' @export
#'
GLENV <- function(formula, data, rank = "all",  family = c("binomial", "poisson"), se = TRUE, init = NULL){

  family <- match.arg(family)
  X = model.matrix(formula, data)[,-1]
  Y = model.frame(formula, data)
  Y = model.response(Y)
  Y <- as.matrix(Y)
  if (rank == "all" || is.null(rank)){
    rank = ncol(X)
  }
  if (family == "binomial"){
    GLENV.binom(formula, data, rank, se, init)
  }
  else if (family == "poisson"){
    GLENV.poisson(formula, data, rank,  se, init)
  }
}


GLENV.binom <- function(formula, data, rank,se = TRUE,  init = NULL) {

  X = model.matrix(formula, data)[,-1]
  Y = model.frame(formula, data)
  Y = as.matrix(as.numeric(as.factor(model.response(Y)))-1)
  a <- dim(Y)
  n <- a[1]
  r <- a[2]
  p <- ncol(X)
  if (a[1] != nrow(X))
    stop("X and Y should have the same number of observations.")
  if (rank > p || rank < 0)
    stop("rank must be an integer between 0 and p.")
  if (sum(duplicated(cbind(X, Y), MARGIN = 2)) > 0)
    stop("Some responses also appear in the predictors, or there maybe duplicated columns in X or Y.")


  SigmaY <- var(Y)*(n-1)/n
  SigmaYX <- cov(Y,X)*(n-1)/n
  SigmaX <- solve(.ridgeinv(cov(X))*((n-1)/n))
  TauY <- .ridgeinv(SigmaY)
  tmp <- binom.envlp_MU(X, Y, rank)

  if (!is.null(init)) {
    if (nrow(init) != p || ncol(init) != rank) stop("The initial value should have p rows and rank columns.")
    tmp0 <- qr.Q(qr(init), complete = TRUE)
    tmp$Gammahat <- tmp0[, 1:rank]
    tmp$Gamma0hat <- tmp0[, (rank+1):p]
    GX <- X %*% tmp$Gammahat
    fit <- cvreg:::.biasreducedglm(cbind(1,GX), as.vector(Y), family = binomial())
    tmp$muhat <- as.vector(fit$coef[1])
    tmp$etahat <- matrix(fit$coef[2 : (rank + 1)])
    beta <- tmp$Gammahat %*% tmp$etahat
    theta <- matrix(1, n, 1) %*% tmp$muhat + X %*% beta
    c.theta <- - exp(theta) / (1 + exp(theta))^2
    c.theta.mean <- mean(c.theta)
    tmp$weighthat <- c.theta / c.theta.mean
    mu1 <- exp(theta) / (1 + exp(theta))
    tmp$Vhat <- theta + ((Y - mu1) / tmp$weighthat)
    W <- diag(as.vector(tmp$weighthat))
    wx <- W %*% X
    mean.wx <- apply(wx, 2, mean)
    wxx <- X - matrix(1, nrow = n) %*% mean.wx
    Sigmawx <- crossprod(wxx, W) %*% wxx / n
    wv <- W %*% tmp$Vhat
    mean.wv <- mean(wv)
    wvv <- tmp$Vhat - matrix(1, nrow = n) %*% mean.wv
    Sigmawxv <- crossprod(wxx, W) %*% wvv / n
    Tauwx <- .ridgeinv(Sigmawx)
    tmp$avar <- Tauwx / (- c.theta.mean)
    M <- SigmaX
    MU <- SigmaX
    tmp.MU <- eigen(MU)
    invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1/tmp.MU$values, "*") ,tmp.MU$vectors)
    e1 <- eigen(crossprod(tmp$Gammahat, M) %*% tmp$Gammahat)
    e2 <- eigen(crossprod(tmp$Gammahat, invMU) %*% tmp$Gammahat)
    e3 <- eigen(mpd(M))
    e4 <- matrix(1, n, 1) %*% tmp$muhat + X %*% tmp$Gammahat %*% tmp$etahat
    e5 <- crossprod(Y, e4) - colSums(log(1 + exp(e4)))
    tmp$objfun <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))) + e5
  }
  Gammahat <- tmp$Gammahat
  Gamma0hat <- tmp$Gamma0hat
  objfun <- tmp$objfun
  muhat <- tmp$muhat
  etahat <- tmp$etahat
  weighthat <- tmp$weighthat
  Vhat <- tmp$Vhat
  avarhat <- tmp$avar

  covMatrix <- NULL
  std.err <- NULL
  ratio <- NULL
  if (rank == 0) {
    Omegahat <- NULL
    Omega0hat <- SigmaX
    betahat <- matrix(0, p, 1)
    SigmaXhat <- SigmaX
    log_lik <- objfun
    if (se == T)
      ratio <- matrix(1, p, r)
  }
  else if (rank == p) {
    TauX <- .ridgeinv(SigmaX)
    Omegahat <- SigmaX
    Omega0hat <- NULL
    betahat <- etahat
    SigmaXhat <- SigmaX
    log_lik <- objfun
    if (se == T) {
      covMatrix <- avarhat
      std.err <- matrix(sqrt(diag(covMatrix)), nrow = p)
      ratio <- matrix(1, p, r)
    }
  }
  else {
    TauX <- .ridgeinv(SigmaX)
    Omegahat <- crossprod(Gammahat, SigmaX) %*% Gammahat
    Omega0hat <- crossprod(Gamma0hat, SigmaX) %*% Gamma0hat
    betahat <- Gammahat %*% etahat
    SigmaXhat <- Gammahat %*% tcrossprod(Omegahat, Gammahat) + Gamma0hat %*% tcrossprod(Omega0hat, Gamma0hat)
    betaOLS <- tcrossprod(TauX, SigmaYX)
    PGamma <- tcrossprod(Gammahat)
    log_lik <- objfun
    if (se == T) {
      options(warn = -1)
      fit <- cvreg:::.biasreducedglm(cbind(1,X), as.vector(Y), family = binomial())
      mu <- as.vector(fit$coef[1])
      beta <- as.vector(fit$coef[2 : (p + 1)])
      theta <- matrix(1, n, 1) %*% mu + X %*% beta
      c.theta <- - exp(theta) / (1 + exp(theta))^2
      c.theta.mean <- mean(c.theta)
      weight <- c.theta / c.theta.mean
      mu1 <- exp(theta) / (1 + exp(theta))
      V <- theta + ((Y - mu1) / weight)
      W <- diag(as.vector(weight))
      wx <- W %*% X
      mean.wx <- apply(wx, 2, mean)
      wxx <- X - matrix(1, nrow = n) %*% mean.wx
      Sigmawx <- crossprod(wxx, W) %*% wxx / n
      wv <- W %*% V
      mean.wv <- mean(wv)
      wvv <- V - matrix(1, nrow = n) %*% mean.wv
      Sigmawxv <- crossprod(wxx, W) %*% wvv / n
      Tauwx <- .ridgeinv(Sigmawx)
      avar <- Tauwx / (- c.theta.mean)
      covMatrix <- avar
      seFm <- matrix(sqrt(diag(covMatrix)), nrow = p)
      PGamma <- tcrossprod(Gammahat)
      t1 <- avarhat
      T1 <- PGamma %*% t1 %*% PGamma
      invt1 <- .ridgeinv(t1)
      t2 <- kronecker(etahat, t(Gamma0hat))
      invOme <- .ridgeinv(Omegahat)
      invOme0 <- .ridgeinv(Omega0hat)
      t3 <- kronecker(Omegahat, invOme0) + kronecker(invOme, Omega0hat) - 2 * kronecker(diag(rank), diag(p - rank))
      M <- t2 %*% tcrossprod(invt1, t2) + t3
      invM <- .ridgeinv(M)
      T2 <- crossprod(t2, invM) %*% t2
      covMatrix <- T1 + T2
      std.err <- matrix(sqrt(diag(covMatrix)), nrow = p)
      ratio <- seFm/std.err
    }
  }
  rownames(betahat) <- colnames(X)
  predictors <- X %*% Gammahat
  fitted <- X %*% betahat
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  colnames(Gammahat) <- paste0("SP", 1:ncol(predictors))
  rownames(Gammahat) <- colnames(X)
  if (se){
    std.err <- std.err * (1/sqrt(n))
    p.values <- 2*pnorm(-abs((betahat / std.err)))
  } else {
    std.err <- NULL
    p.values <- NULL
  }
  return(structure(list(basis = Gammahat, coef = betahat, fitted = fitted, predictors = predictors, eta = etahat, Omega = Omegahat, SigmaX = SigmaXhat, log_lik = log_lik,  n = n, covMatrix = covMatrix, std.err = std.err, p.values = p.values, ratio = ratio, formula = formula, mf = model.frame(formula, data), y.obs = Y), class = "sdr", model = "ENV"))

}



GLENV.poisson <- function(X, Y, rank, se = TRUE, init = NULL){

  X <- as.matrix(X)
  Y <- as.matrix(Y)
  a <- dim(Y)
  n <- a[1]
  r <- a[2]
  p <- ncol(X)
  if (a[1] != nrow(X))
    stop("X and Y should have the same number of observations.")
  if (rank > p || rank < 0)
    stop("rank must be an integer between 0 and p.")
  if (sum(duplicated(cbind(X, Y), MARGIN = 2)) > 0)
    stop("Some responses also appear in the predictors, or there maybe duplicated columns in X or Y.")

  SigmaY <- var(Y)  * ((n - 1) / n)
  SigmaYX <- cov(Y, X)  * ((n - 1) / n)
  SigmaX <- solve(.ridgeinv(cov(X)))*((n-1)/n)
  TauY <- .ridgeinv(SigmaY)
  tmp <- poisson.envlp_MU(X, Y, rank)

  if (!is.null(init)) {
    if (nrow(init) != p || ncol(init) != rank) stop("The initial value should have p rows and rank columns.")
    tmp0 <- qr.Q(qr(init), complete = TRUE)
    tmp$Gammahat <- tmp0[, 1:rank]
    tmp$Gamma0hat <- tmp0[, (rank+1):p]
    GX <- X %*% tmp$Gammahat
    fit <- cvreg:::.biasreducedglm(cbind(1,GX), as.vector(Y), family = poisson())
    tmp$muhat <- as.vector(fit$coef[1])
    tmp$etahat <- matrix(fit$coef[2 : (rank + 1)])
    beta <- tmp$Gammahat %*% tmp$etahat
    theta <- matrix(1, n, 1) %*% tmp$muhat + X %*% beta
    c.theta <- - exp(theta)
    c.theta.mean <- mean(c.theta)
    tmp$weighthat <- c.theta / c.theta.mean
    tmp$Vhat <- theta + ((Y - c.theta) / tmp$weighthat)
    W <- diag(as.vector(tmp$weighthat))
    wx <- W %*% X
    mean.wx <- apply(wx, 2, mean)
    wxx <- X - matrix(1, nrow = n) %*% mean.wx
    Sigmawx <- crossprod(wxx, W) %*% wxx / n
    wv <- W %*% tmp$Vhat
    mean.wv <- mean(wv)
    wvv <- tmp$Vhat - matrix(1, nrow = n) %*% mean.wv
    Sigmawxv <- crossprod(wxx, W) %*% wvv / n
    Tauwx <- .ridgeinv(Sigmawx)
    tmp$avar <- Tauwx / (- c.theta.mean)
    M <- SigmaX
    MU <- SigmaX
    tmp.MU <- eigen(MU)
    invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1/tmp.MU$values, "*") ,tmp.MU$vectors)
    e1 <- eigen(crossprod(tmp$Gammahat, M) %*% tmp$Gammahat)
    e2 <- eigen(crossprod(tmp$Gammahat, invMU) %*% tmp$Gammahat)
    e3 <- eigen(mpd(M))
    e4 <- matrix(1, n, 1) %*% tmp$muhat + X %*% tmp$Gammahat %*% tmp$etahat
    e5 <- crossprod(Y, e4) - colSums(exp(e4))
    tmp$objfun <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))) + e5
  }

  Gammahat <- tmp$Gammahat
  Gamma0hat <- tmp$Gamma0hat
  objfun <- tmp$objfun
  muhat <- tmp$muhat
  etahat <- tmp$etahat
  weighthat <- tmp$weighthat
  Vhat <- tmp$Vhat
  avarhat <- tmp$avar

  covMatrix <- NULL
  std.err <- NULL
  ratio <- NULL
  if (rank == 0) {
    Omegahat <- NULL
    Omega0hat <- SigmaX
    betahat <- matrix(0, p, 1)
    SigmaXhat <- SigmaX
    log_lik <- objfun
    if (se == T)
      ratio <- matrix(1, p, r)
  }
  else if (rank == p) {
    Omegahat <- SigmaX
    Omega0hat <- NULL
    betahat <- etahat
    SigmaXhat <- SigmaX
    log_lik <- objfun
    if (se == T) {
      covMatrix <- avarhat
      std.err <- matrix(sqrt(diag(covMatrix)), nrow = p)
      ratio <- matrix(1, p, r)
    }
  }
  else {
    TauX<- .ridgeinv(SigmaX)
    Omegahat <- crossprod(Gammahat, SigmaX) %*% Gammahat
    Omega0hat <- crossprod(Gamma0hat, SigmaX) %*% Gamma0hat
    betahat <- Gammahat %*% etahat
    SigmaXhat <- Gammahat %*% tcrossprod(Omegahat, Gammahat) + Gamma0hat %*% tcrossprod(Omega0hat, Gamma0hat)
    PGamma <- tcrossprod(Gammahat)
    log_lik <- objfun
    if (se == T) {
      options(warn = -1)
      fit <- cvreg:::.biasreducedglm(cbind(1,X), as.vector(Y), family = poisson())
      mu <- as.vector(fit$coef[1])
      beta <- as.vector(fit$coef[2 : (p + 1)])
      theta <- matrix(1, n, 1) %*% mu + X %*% beta
      c.theta <- - exp(theta)
      c.theta.mean <- mean(c.theta)
      weight <- c.theta / c.theta.mean
      V <- theta + ((Y - c.theta) / weight)
      W <- diag(as.vector(weight))
      wx <- W %*% X
      mean.wx <- apply(wx, 2, mean)
      wxx <- X - matrix(1, nrow = n) %*% mean.wx
      Sigmawx <- crossprod(wxx, W) %*% wxx / n
      wv <- W %*% V
      mean.wv <- mean(wv)
      wvv <- V - matrix(1, nrow = n) %*% mean.wv
      Sigmawxv <- crossprod(wxx, W) %*% wvv / n
      Tauwx <- .ridgeinv(Sigmawx)
      avar <- Tauwx / (- c.theta.mean)
      covMatrix <- avar
      asyFm <- matrix(sqrt(diag(covMatrix)), nrow = p)
      PGamma <- tcrossprod(Gammahat)
      t1 <- avarhat
      T1 <- PGamma %*% t1 %*% PGamma
      invt1 <- .ridgeinv(t1)
      t2 <- kronecker(etahat, t(Gamma0hat))
      invOme <- .ridgeinv(Omegahat)
      invOme0 <- .ridgeinv(Omega0hat)
      t3 <- kronecker(Omegahat, invOme0) + kronecker(invOme, Omega0hat) - 2 * kronecker(diag(rank), diag(p - rank))
      M <- t2 %*% tcrossprod(invt1, t2) + t3
      invM <- .ridgeinv(M)
      T2 <- crossprod(t2, invM) %*% t2
      covMatrix <- T1 + T2
      std.err <- matrix(sqrt(diag(covMatrix)), nrow = p)
      ratio <- asyFm/std.err
    }
  }
  rownames(betahat) <- colnames(X)
  predictors <- X %*% Gammahat
  fitted <- X %*% betahat
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  colnames(Gammahat) <- paste0("SP", 1:ncol(predictors))
  rownames(Gammahat) <- colnames(X)
  if (se){
    std.err <- std.err * (1/sqrt(n))
    p.values <- 2*pnorm(-abs((betahat / std.err)))
  } else {
    std.err <- NULL
    p.values <- NULL
  }
  return(structure(list(basis = Gammahat, coef = betahat, fitted = fitted, predictors = predictors, eta = etahat, Omega = Omegahat, SigmaX = SigmaXhat, log_lik = log_lik, n = n, covMatrix = covMatrix, std.err = std.err, p.values = p.values, ratio = ratio, formula = formula, mf = model.frame(formula, data), y.obs = Y), class = "sdr", model = "ENV"))
}


## MU functions

envlp_MU <- function(M, U, rank) {

  dimM <- dim(M)
  dimU <- dim(U)
  r <- dimM[1]
  if (dimM[1] != dimM[2] & dimU[1] != dimU[2]) stop("M and U should be square matrices.")
  if (dimM[1] != dimU[1]) stop("M and U should have the same dimension.")
  if (rank > r & rank < 0) stop("rank should be between 0 and r.")
  if (rank == r) {
    Gammahat <- diag(r)
    Gamma0hat <- NULL
    tmp.M <- eigen(mpd(M))
    objfun <- sum(log(tmp.M$values))
  } else if (rank == 1) {
    maxiter = 250
    ftol = 0.0005
    MU <- M + U
    tmp.MU <- eigen(MU)
    invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1 / tmp.MU$values, '*') , tmp.MU$vectors)
    invMU2 <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1 / sqrt(tmp.MU$values), '*') , tmp.MU$vectors)

    midmatrix <- U
    startv <- function(a) crossprod(a, midmatrix) %*% a
    tmp2.MU <- apply(tmp.MU$vectors, 2, startv)
    tmp3.MU <- sort(tmp2.MU, decreasing = TRUE, index.return = TRUE)
    init <- as.matrix(tmp.MU$vectors[, tmp3.MU$ix[1]])

    eig1 <- eigen(crossprod(init,  M) %*% init)
    eig2 <- eigen(crossprod(init, invMU) %*% init)
    obj1 <- sum(log(eig1$values)) + sum(log(eig2$values))

    midmatrix <- invMU2 %*% tcrossprod(U, invMU2)
    tmp2.MU <- apply(tmp.MU$vectors, 2, startv)
    tmp3.MU <- sort(tmp2.MU, decreasing = TRUE, index.return = TRUE)
    init.MU <- as.matrix(tmp.MU$vectors[, tmp3.MU$ix[1]])
    e1 <- eigen(crossprod(init.MU, M) %*% init.MU)
    e2 <- eigen(crossprod(init.MU, invMU) %*% init.MU)
    obj2 <- sum(log(e1$values)) + sum(log(e2$values))
    if (obj2 < obj1) {
      init <- init.MU
      obj1 <- obj2
    }

    tmp.M <- eigen(mpd(M))
    midmatrix <- U
    tmp2.M <- apply(tmp.M$vectors, 2, startv)
    tmp3.M <- sort(tmp2.M, decreasing = TRUE, index.return = TRUE)
    init.M <- as.matrix(tmp.M$vectors[, tmp3.M$ix[1]])
    e1 <- eigen(crossprod(init.M, M) %*% init.M)
    e2 <- eigen(crossprod(init.M, invMU) %*% init.M)
    obj3 <- sum(log(e1$values)) + sum(log(e2$values))
    if (obj3 < obj1) {
      init <- init.M
      obj1 <- obj3
    }

    invM2 <- tcrossprod(sweep(tmp.M$vectors, MARGIN = 2, 1 / sqrt(tmp.M$values), '*') , tmp.M$vectors)
    midmatrix <- invM2 %*% tcrossprod(U, invM2)
    tmp2.M <- apply(tmp.M$vectors, 2, startv)
    tmp3.M <- sort(tmp2.M, decreasing = TRUE, index.return = TRUE)
    init.M <- as.matrix(tmp.M$vectors[, tmp3.M$ix[1]])

    e1 <- eigen(crossprod(init.M, M) %*% init.M)
    e2 <- eigen(crossprod(init.M, invMU) %*% init.M)
    obj4 <- sum(log(e1$values)) + sum(log(e2$values))
    if (obj4 < obj1) {
      init <- init.M
      obj1 <- obj4
    }
    GEidx <- .GaussElim(init)
    Ginit <- init %*% .ridgeinv(init[GEidx[1], ])
    i <- 1
    while (i < maxiter) {
      fobj <- function(x) {
        T1 <- crossprod(x, x)
        T2 <- crossprod(x, M) %*% x
        T3 <- crossprod(x, invMU) %*% x
        -2 * log(T1) + log(T2) + log(T3)
      }

      gobj <- function(x) {
        W1 <- crossprod(x, x)
        T1 <- x / as.vector(W1)
        W2 <- crossprod(x, M) %*% x
        T2 <- M %*% x / as.vector(W2)
        W3 <- crossprod(x, invMU) %*% x
        T3 <- invMU %*% x / as.vector(W3)
        -2 * T1 + T2 + T3
      }

      res <- nlminb(Ginit, fobj, gobj)
      g <- as.matrix(res$par)
      a <- qr(g)
      Gammahat <- qr.Q(a)
      e1 <- eigen(crossprod(Gammahat, M) %*% Gammahat)
      e2 <- eigen(crossprod(Gammahat, invMU) %*% Gammahat)
      obj5 <- sum(log(e1$values)) + sum(log(e2$values))
      if (abs(obj1 - obj5) < ftol * abs(obj1)) {
        break
      } else {
        obj1 <- obj5
        i <- i + 1
      }
    }
    Gamma0hat <- qr.Q(a, complete = TRUE)[, (rank+1):r]
    objfun <- obj5 + sum(log(tmp.MU$values))
    Gammahat <- as.matrix(Gammahat)
    Gamma0hat <- as.matrix(Gamma0hat)

  } else if (rank == r - 1 & rank != 1) {

    maxiter = 250
    ftol = 0.0005
    MU <- M + U
    tmp.MU <- eigen(MU)
    invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1 / tmp.MU$values, '*') , tmp.MU$vectors)
    invMU2 <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1 / sqrt(tmp.MU$values), '*') , tmp.MU$vectors)
    midmatrix <- U
    startv <- function(a) crossprod(a, midmatrix) %*% a
    tmp2.MU <- apply(tmp.MU$vectors, 2, startv)
    tmp3.MU <- sort(tmp2.MU, decreasing = TRUE, index.return = TRUE)
    init <- as.matrix(tmp.MU$vectors[, tmp3.MU$ix[1:rank]])
    eig1 <- eigen(crossprod(init, M) %*% init)
    eig2 <- eigen(crossprod(init, invMU) %*% init)
    obj1 <- sum(log(eig1$values)) + sum(log(eig2$values))
    midmatrix <- invMU2 %*% tcrossprod(U, invMU2)
    tmp2.MU <- apply(tmp.MU$vectors, 2, startv)
    tmp3.MU <- sort(tmp2.MU, decreasing = TRUE, index.return = TRUE)
    init.MU <- as.matrix(tmp.MU$vectors[, tmp3.MU$ix[1:rank]])
    e1 <- eigen(crossprod(init.MU, M) %*% init.MU)
    e2 <- eigen(crossprod(init.MU, invMU) %*% init.MU)
    obj2 <- sum(log(e1$values)) + sum(log(e2$values))
    if (obj2 < obj1) {
      init <- init.MU
      obj1 <- obj2
    }
    tmp.M <- eigen(mpd(M))
    midmatrix <- U
    tmp2.M <- apply(tmp.M$vectors, 2, startv)
    tmp3.M <- sort(tmp2.M, decreasing = TRUE, index.return = TRUE)
    init.M <- as.matrix(tmp.M$vectors[, tmp3.M$ix[1:rank]])
    e1 <- eigen(crossprod(init.M, M) %*% init.M)
    e2 <- eigen(crossprod(init.M, invMU) %*% init.M)
    obj3 <- sum(log(e1$values)) + sum(log(e2$values))
    if (obj3 < obj1) {
      init <- init.M
      obj1 <- obj3
    }

    invM2 <- tcrossprod(sweep(tmp.M$vectors, MARGIN = 2, 1 / sqrt(tmp.M$values), '*') , tmp.M$vectors)
    midmatrix <- invM2 %*% tcrossprod(U, invM2)
    tmp2.M <- apply(tmp.M$vectors, 2, startv)
    tmp3.M <- sort(tmp2.M, decreasing = TRUE, index.return = TRUE)
    init.M <- as.matrix(tmp.M$vectors[, tmp3.M$ix[1:rank]])

    e1 <- eigen(crossprod(init.M, M) %*% init.M)
    e2 <- eigen(crossprod(init.M, invMU) %*% init.M)
    obj4 <- sum(log(e1$values)) + sum(log(e2$values))
    if (obj4 < obj1) {
      init <- init.M
      obj1 <- obj4
    }

    GEidx <- .GaussElim(init)
    Ginit <- init %*% pseudoinverse(init[GEidx[1:rank], ])

    j <- GEidx[r]

    g <- as.matrix(Ginit[j, ])
    t2 <- crossprod(Ginit[-j, ], as.matrix(M[-j, j])) / M[j, j]
    t3 <- crossprod(Ginit[-j, ], as.matrix(invMU[-j, j])) / invMU[j, j]

    GUGt2 <- g + t2
    GUG <- crossprod(Ginit, (M %*% Ginit)) - tcrossprod(GUGt2, GUGt2) * M[j, j]

    GVGt2 <- g + t3
    GVG <- crossprod(Ginit, (invMU %*% Ginit)) - tcrossprod(GVGt2, GVGt2) * invMU[j, j]

    invC1 <- .ridgeinv(GUG);invC2 <- .ridgeinv(GVG)

    fobj <- function(x) {
      tmp2 <- x + t2
      tmp3 <- x + t3
      T2 <- invC1 %*% tmp2
      T3 <- invC2 %*% tmp3
      -2 * log(1 + sum(x^2)) + log(1 + M[j, j] * crossprod(tmp2, T2)) + log(1 + invMU[j, j] * crossprod(tmp3, T3))
    }

    gobj <- function(x) {
      tmp2 <- x + t2
      tmp3 <- x + t3
      T2 <- invC1 %*% tmp2
      T3 <- invC2 %*% tmp3
      -4 	* x %*% .ridgeinv(1 + sum(x^2)) + 2 * T2 / as.numeric(1 / M[j, j] + crossprod(tmp2, T2)) + 2 * T3 / as.numeric(1 / invMU[j, j] + crossprod(tmp3, T3))
    }

    i <- 1
    while (i < maxiter) {
      res <- nlminb(Ginit[j,], fobj, gobj)
      Ginit[j, ] <- res$par
      a <- qr(Ginit)
      Gammahat <- qr.Q(a)
      e1 <- eigen(crossprod(Gammahat, M) %*% Gammahat)
      e2 <- eigen(crossprod(Gammahat, invMU) %*% Gammahat)
      obj5 <- sum(log(e1$values)) + sum(log(e2$values))
      if (abs(obj1 - obj5) < ftol * abs(obj1)) {
        break
      } else {
        obj1 <- obj5
        i <- i + 1
      }
    }

    Gamma0hat <- qr.Q(a, complete = TRUE)[, (rank+1):r, drop = FALSE]
    objfun <- obj5 + sum(log(tmp.MU$values))
    Gammahat <- as.matrix(Gammahat)
    Gamma0hat <- as.matrix(Gamma0hat)

  } else {
    maxiter = 250
    ftol = 0.0005
    MU <- M + U
    tmp.MU <- eigen(MU)
    invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1 / tmp.MU$values, '*') , tmp.MU$vectors)
    invMU2 <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1 / sqrt(tmp.MU$values), '*') , tmp.MU$vectors)
    midmatrix <- U
    startv <- function(a) crossprod(a, midmatrix) %*% a
    tmp2.MU <- apply(tmp.MU$vectors, 2, startv)
    tmp3.MU <- sort(tmp2.MU, decreasing = TRUE, index.return = TRUE)
    init <- as.matrix(tmp.MU$vectors[, tmp3.MU$ix[1:rank]])
    eig1 <- eigen(crossprod(init, M) %*% init)
    eig2 <- eigen(crossprod(init, invMU) %*% init)
    obj1 <- sum(log(eig1$values)) + sum(log(eig2$values))
    midmatrix <- invMU2 %*% tcrossprod(U, invMU2)
    tmp2.MU <- apply(tmp.MU$vectors, 2, startv)
    tmp3.MU <- sort(tmp2.MU, decreasing = TRUE, index.return = TRUE)
    init.MU <- as.matrix(tmp.MU$vectors[, tmp3.MU$ix[1:rank]])
    e1 <- eigen(crossprod(init.MU, M) %*% init.MU)
    e2 <- eigen(crossprod(init.MU, invMU) %*% init.MU)
    obj2 <- sum(log(e1$values)) + sum(log(e2$values))
    if (obj2 < obj1) {
      init <- init.MU
      obj1 <- obj2
    }

    tmp.M <- eigen(mpd(M))
    midmatrix <- U
    tmp2.M <- apply(tmp.M$vectors, 2, startv)
    tmp3.M <- sort(tmp2.M, decreasing = TRUE, index.return = TRUE)
    init.M <- as.matrix(tmp.M$vectors[, tmp3.M$ix[1:rank]])
    e1 <- eigen(crossprod(init.M, M) %*% init.M)
    e2 <- eigen(crossprod(init.M, invMU) %*% init.M)
    obj3 <- sum(log(e1$values)) + sum(log(e2$values))
    if (obj3 < obj1) {
      init <- init.M
      obj1 <- obj3
    }

    invM2 <- tcrossprod(sweep(tmp.M$vectors, MARGIN = 2, 1 / sqrt(tmp.M$values), '*') , tmp.M$vectors)
    midmatrix <- invM2 %*% tcrossprod(U, invM2)
    tmp2.M <- apply(tmp.M$vectors, 2, startv)
    tmp3.M <- sort(tmp2.M, decreasing = TRUE, index.return = TRUE)
    init.M <- as.matrix(tmp.M$vectors[, tmp3.M$ix[1:rank]])
    e1 <- eigen(crossprod(init.M, M) %*% init.M)
    e2 <- eigen(crossprod(init.M, invMU) %*% init.M)
    obj4 <- sum(log(e1$values)) + sum(log(e2$values))
    if (obj4 < obj1) {
      init <- init.M
      obj1 <- obj4
    }
    GEidx <- .GaussElim(init)
    Ginit <- init %*% pseudoinverse(init[GEidx[1:rank], ])
    GUG <- crossprod(Ginit, (M %*% Ginit))
    GVG <- crossprod(Ginit, (invMU %*% Ginit))
    t4 <- crossprod(Ginit[GEidx[(rank+1):r],], Ginit[GEidx[(rank+1):r], ]) + diag(rank)
    i <- 1
    while (i < maxiter) {

      for (j in GEidx[(rank+1):r]) {
        g <- as.matrix(Ginit[j, ])
        t2 <- crossprod(Ginit[-j, ], as.matrix(M[-j, j])) / M[j, j]
        t3 <- crossprod(Ginit[-j, ], as.matrix(invMU[-j, j])) / invMU[j, j]

        GUGt2 <- g + t2
        GUG <- GUG - tcrossprod(GUGt2, GUGt2) * M[j, j]
        GVGt2 <- g + t3
        GVG <- GVG - tcrossprod(GVGt2, GVGt2) * invMU[j, j]
        t4 <- t4 - tcrossprod(g, g)
        invC1 <- .ridgeinv(GUG)
        invC2 <- .ridgeinv(GVG)
        invt4 <- .ridgeinv(t4)

        fobj <- function(x) {
          tmp2 <- x + t2
          tmp3 <- x + t3
          T1 <- invt4 %*% x
          T2 <- invC1 %*% tmp2
          T3 <- invC2 %*% tmp3
          -2 * log(1 + x %*% T1) + log(1 + M[j, j] * crossprod(tmp2, T2)) + log(1 + invMU[j, j] * crossprod(tmp3, T3))
        }

        gobj <- function(x) {
          tmp2 <- x + t2
          tmp3 <- x + t3
          T1 <- invt4 %*% x
          T2 <- invC1 %*% tmp2
          T3 <- invC2 %*% tmp3
          -4 	* T1 / as.numeric(1 + x %*% T1) + 2 * T2 / as.numeric(1 / M[j, j] + crossprod(tmp2, T2)) + 2 * T3 / as.numeric(1 / invMU[j, j] + crossprod(tmp3, T3))
        }
        res <- nlminb(Ginit[j,], fobj, gobj)
        Ginit[j, ] <- res$par
        g <- as.matrix(Ginit[j, ])
        t4 <- t4 + tcrossprod(g, g)
        GUGt2 <- g + t2
        GUG <- GUG + tcrossprod(GUGt2, GUGt2) * M[j, j]
        GVGt2 <- g + t3
        GVG <- GVG + tcrossprod(GVGt2, GVGt2) * invMU[j, j]
      }
      a <- qr(Ginit)
      Gammahat <- qr.Q(a)
      e1 <- eigen(crossprod(Gammahat, M) %*% Gammahat)
      e2 <- eigen(crossprod(Gammahat, invMU) %*% Gammahat)
      obj5 <- sum(log(e1$values)) + sum(log(e2$values))
      if (abs(obj1 - obj5) < ftol * abs(obj1)) {
        break
      } else {
        obj1 <- obj5
        i <- i + 1
      }
    }
    Gamma0hat <- qr.Q(a, complete = TRUE)[, (rank+1):r]
    objfun <- obj5 + sum(log(tmp.MU$values))
    Gammahat <- as.matrix(Gammahat)
    Gamma0hat <- as.matrix(Gamma0hat)
  }
  return(list(Gammahat = Gammahat, Gamma0hat = Gamma0hat, objfun = objfun))
}



sxenvlp_MU <- function(X, Y, rank, R){

  X <- as.matrix(X)
  Y <- as.matrix(Y)
  a <- dim(Y)
  n <- a[1]
  r <- a[2]
  p <- ncol(X)
  if (a[1] != nrow(X))
    stop("X and Y should have the same number of observations.")
  if (rank > p || rank < 0)
    stop("rank must be an integer between 0 and p.")
  if (sum(duplicated(cbind(X, Y), MARGIN = 2)) > 0)
    stop("Some responses also appear in the predictors, or there maybe duplicated columns in X or Y.")

  SigmaY  <- cov(Y)
  SigmaYX <- cov(Y, X)* ((n-1)/n)
  SigmaX  <- solve(.ridgeinv(cov(X))) * ((n-1)/n)
  TauX <- .ridgeinv(SigmaX);TauY <- .ridgeinv(SigmaY);t1 <- crossprod(SigmaYX, TauY)
  U <- t1 %*% SigmaYX;M <- SigmaX - U;tmp <- envlp_MU(M, U, rank);Gamma <- tmp$Gammahat

  if (length(R) == p) {
    make_objfun <- function(SigmaYX, SigmaY, SigmaX){
      function(d, Gamma, X, Y){
        X <- as.matrix(X);Y <- as.matrix(Y);a <- dim(Y);n <- a[1];r <- a[2];p <- ncol(X)
        TauX <- .ridgeinv(SigmaX);TauY <- .ridgeinv(SigmaY)
        t1 <- crossprod(SigmaYX, TauY);U <- t1 %*% SigmaYX;M <- SigmaX - U
        d1 <- c(1, d);Lambda <- diag(d1);invLambda <- diag(1 / d1)
        m1 <- crossprod(Gamma, Lambda);eigtem1 <- eigen(m1 %*% tcrossprod(TauX, m1))
        m2 <- crossprod(Gamma, invLambda);eigtem2 <- eigen(m2 %*% tcrossprod(M, m2))
        temp1 <- sum(log(eigtem1$values));temp2 <- sum(log(eigtem2$values));objfun <- temp1 + temp2
        return(objfun)
      }
    }
    objfun <- make_objfun(SigmaYX, SigmaY, SigmaX);d.init <- rep(1, (p - 1));k2 <- rep(0, (p - 1))
    tmp.init <- Rsolnp::solnp(pars = d.init, fun = objfun, LB = k2,
                              .control = list(delta = 1e-6, tol = 0.0005, trace = 0),
                              Gamma = Gamma, X = X, Y = Y)
    d <- tmp.init$pars;obj1 <- objfun(d, Gamma, X, Y);maxiter = 250;ftol = 0.0005;i <- 1
    while (i < maxiter) {
      d1 <- c(1, d)
      L1 <- diag(d1)
      invL1 <- diag(1 / d1)
      M1 <- invL1 %*% M %*% invL1
      U1 <- invL1 %*% U %*% invL1
      tmp1 <- envlp_MU(M1, U1, rank)
      Gamma <- tmp1$Gammahat
      temp1 <- Rsolnp::solnp(pars = d, fun = objfun, LB = k2,
                             .control = list(delta = 1e-6, tol = 0.0005, trace = 0),
                             Gamma = Gamma, X = X, Y = Y)
      d <- temp1$pars
      obj2 <- objfun(d, Gamma, X, Y)
      if (abs(obj1 - obj2) < ftol) {
        break
      }
      else {
        obj1 <- obj2
        i <- i + 1
      }
    }
    Gamma0 <- tmp1$Gamma0hat
    d1 <- c(1, d)
    Lambda <- diag(d1)
  } else {

    make_objfun1 <- function(SigmaYX, SigmaY, SigmaX){
      function(d, Gamma, X, Y, R){
        X <- as.matrix(X);Y <- as.matrix(Y);a <- dim(Y);n <- a[1];r <- a[2];p <- ncol(X)
        TauX <- .ridgeinv(SigmaX);TauY <- .ridgeinv(SigmaY);t1 <- crossprod(SigmaYX, TauY)
        U <- t1 %*% SigmaYX;M <- SigmaX - U;Lambda <- diag(d);invLambda <- diag(1 / d)
        m1 <- crossprod(Gamma, Lambda);eigtem1 <- eigen(m1 %*% tcrossprod(TauX, m1))
        m2 <- crossprod(Gamma, invLambda);eigtem2 <- eigen(m2 %*% tcrossprod(M, m2))
        temp1 <- sum(log(eigtem1$values));temp2 <- sum(log(eigtem2$values));objfun <- temp1 + temp2
        return(objfun)
      }
    }

    objfun1 <- make_objfun1(SigmaYX, SigmaY, SigmaX)

    heq <- function (d, Gamma, X, Y, R){
      iter = sum(R)
      i <- 1
      cont <- NULL
      while (i < iter) {
        C <- matrix(0, (R[i] - 1), sum(R))
        for (j in 2 : (sum(R) - 1)){
          if (R[i] == j) {
            for (k in 1 : (j - 1)){
              s <- sum(R[0 : (i - 1)]) + 1
              C[k , s] <- 1
              C[k , (s + k)] <- -1
            }
          }
        }
        cont <- rbind(cont, C)
        if (i == length(R)) {
          c1 <- c(1, rep(0, (sum(R) - 1)))
          cont <- rbind(c1, cont)
          break
        } else {
          i <- i + 1
        }
      }

      g <- rep(NA, nrow(cont))
      g[1] <- d[1]
      for (m in 2 : nrow(cont)){
        g[m] <- d[which(cont[m, ] == 1)] - d[which(cont[m, ] == -1)]
      }
      g
    }
    d.init <- rep(1, p);g1 <- heq(d.init, Gamma, X, Y, R);k1 <- c(1, rep(0, length(g1) - 1));k2 <- rep(0, p)

    tmp.init <- Rsolnp::solnp(pars = d.init, fun = objfun1, eqfun = heq, eqB = k1,
                              LB = k2, .control = list(delta = 1e-6, tol = 0.0005, trace = 0),
                              R = R, Gamma = Gamma, X = X, Y = Y)

    d <- tmp.init$pars;obj1 <- objfun1(d, Gamma, X, Y, R);maxiter = 250;ftol = 0.0005;i <- 1

    while (i < maxiter) {
      L1 <- diag(d)
      invL1 <- diag(1 / d)
      M1 <- invL1 %*% M %*% invL1
      U1 <- invL1 %*% U %*% invL1
      tmp1 <- envlp_MU(M1, U1, rank)
      Gamma <- tmp1$Gammahat
      temp1 <- Rsolnp::solnp(pars = d, fun = objfun1, eqfun = heq, eqB = k1,
                             LB = k2, .control = list(delta = 1e-6, tol = 0.0005, trace = 0),
                             R = R, Gamma = Gamma, X = X, Y = Y)
      d <- temp1$pars
      obj2 <- objfun1(d, Gamma, X, Y, R)
      if (abs(obj1 - obj2) < ftol) {
        break
      }
      else {
        obj1 <- obj2
        i <- i + 1
      }
    }
    Gamma0 <- tmp1$Gamma0hat
    Lambda <- diag(d)
  }

  return(list(Gammahat = Gamma, Gamma0hat = Gamma0, Lambdahat = Lambda, objfun = obj2))
}


xyenvlp_MU <- function(X, Y, xrank, yrank){

  X <- as.matrix(X);Y <- as.matrix(Y)
  a <- dim(Y);n <- a[1];r <- a[2];p <- ncol(X)
  if (a[1] != nrow(X))
    stop("X and Y should have the same number of observations.")
  if (yrank > r || yrank < 0)
    stop("yrank must be an integer between 0 and r.")
  if (xrank > p || xrank < 0)
    stop("xrank must be an integer between 0 and p.")
  if (sum(duplicated(cbind(X, Y), MARGIN = 2)) > 0)
    stop("There appears to be duplicated columns in X and Y.")

  SigmaY <- solve(.ridgeinv(cov(Y))) * ((n-1)/n)
  SigmaYX <- cov(Y, X)* ((n-1)/n)
  SigmaX <- solve(.ridgeinv(cov(X))) * ((n-1)/n)
  TauX <- .ridgeinv(SigmaX);TauY <- .ridgeinv(SigmaY)
  betaOLS <- SigmaYX %*% TauX
  U1 <- crossprod(SigmaYX, TauY) %*% SigmaYX
  M1 <- SigmaX - U1
  tmp1 <- envlp_MU(M1, U1, xrank)
  Phi <- tmp1$Gammahat
  Phi0 <- tmp1$Gamma0hat
  E1 <- crossprod(Phi, SigmaX) %*% Phi
  invE1 <- .ridgeinv(E1)
  E2 <- SigmaYX %*% Phi
  U2 <- E2 %*% tcrossprod(invE1, E2)
  M2 <- SigmaY - U2
  tmp2 <- envlp_MU(M2, U2, yrank)
  Ga <- tmp2$Gammahat
  Ga0 <- tmp2$Gamma0hat
  m1 <- crossprod(Phi, SigmaX) %*% Phi
  m2 <- crossprod(Phi0, SigmaX) %*% Phi0
  m3 <- crossprod(Ga0, SigmaY) %*% Ga0
  m4 <- crossprod(Ga, M2) %*% Ga
  e1 <- eigen(m1);e2 <- eigen(m2);e3 <- eigen(m3);e4 <- eigen(m4)
  obj1 <- (sum(log(e1$values)) + sum(log(e2$values))  + sum(log(e3$values)) + sum(log(e4$values)))

  maxiter = 250
  ftol = 0.0005
  i <- 1
  while(i < maxiter){
    E3 <- crossprod(Ga, SigmaY) %*% Ga
    invE3 <- .ridgeinv(E3)
    E4 <- crossprod(SigmaYX, Ga)
    U3 <- E4 %*% tcrossprod(invE3, E4)
    M3 <- SigmaX - U3
    tmp3 <- envlp_MU(M3, U3, xrank)
    Phi <- tmp3$Gammahat
    Phi0 <- tmp3$Gamma0hat
    E5 <- crossprod(Phi, SigmaX) %*% Phi
    invE5 <- .ridgeinv(E5)
    E6 <- SigmaYX %*% Phi
    U4 <- E6 %*% tcrossprod(invE5, E6)
    M4 <- SigmaY - U4
    tmp4 <- envlp_MU(M4, U4, yrank)
    Ga <- tmp4$Gammahat
    Ga0 <- tmp4$Gamma0hat
    m5 <- crossprod(Phi, SigmaX) %*% Phi
    m6 <- crossprod(Phi0, SigmaX) %*% Phi0
    m7 <- crossprod(Ga0, SigmaY) %*% Ga0
    m8 <- crossprod(Ga, M4) %*% Ga
    e5 <- eigen(m5);e6 <- eigen(m6);e7 <- eigen(m7);e8 <- eigen(m8)
    obj2 <- (sum(log(e5$values)) + sum(log(e6$values)) + sum(log(e7$values)) + sum(log(e8$values)))
    if (abs(obj1 - obj2) < ftol) {
      break
    }
    else {
      obj1 <- obj2
      i <- i + 1
    }
  }
  return(list(Gammahat = Ga, Gamma0hat = Ga0,Phihat = Phi, Phi0hat = Phi0, objfun = obj2))
}




binom.envlp_MU <- function(X, Y, rank){
  X <- as.matrix(X);Y <- as.matrix(Y)
  a <- dim(Y);n <- a[1];r <- a[2];p <- ncol(X)
  if (a[1] != nrow(X))
    stop("X and Y should have the same number of observations.")
  if (rank > p || rank < 0)
    stop("rank must be an integer between 0 and r.")
  if (sum(duplicated(cbind(X, Y), MARGIN = 2)) > 0)
    stop("Some responses also appear in the predictors, or there maybe duplicated columns in X or Y.")

  options(warn = -1)
  fit <- cvreg:::.biasreducedglm(cbind(1,X), as.vector(Y), family = binomial())
  mu <- as.vector(fit$coef[1])
  beta <- as.vector(fit$coef[2 : (p + 1)])
  theta <- matrix(1, n, 1) %*% mu + X %*% beta
  c.theta <- - exp(theta) / (1 + exp(theta))^2
  c.theta.mean <- mean(c.theta)
  weight <- c.theta / c.theta.mean
  mu1 <- exp(theta) / (1 + exp(theta))
  V <- theta + ((Y - mu1) / weight)
  W <- diag(as.vector(weight))
  wx <- W %*% X
  mean.wx <- apply(wx, 2, mean)
  wxx <- X - matrix(1, nrow = n) %*% mean.wx
  Sigmawx <- crossprod(wxx, W) %*% wxx / n
  wv <- W %*% V
  mean.wv <- mean(wv)
  wvv <- V - matrix(1, nrow = n) %*% mean.wv
  Sigmawxv <- crossprod(wxx, W) %*% wvv / n
  Tauwx <- .ridgeinv(Sigmawx)

  M.init <- Tauwx / (- c.theta.mean)
  U.init <- tcrossprod(beta)
  tmp1 <- envlp_MU(M.init, U.init, rank)
  gamma.init <- tmp1$Gammahat
  gamma0.init <- tmp1$Gamma0hat

  SigmaX <- stats::cov(X) * ((n - 1) / n)
  TauX <- .ridgeinv(SigmaX)
  SigmaX <- solve(SigmaX)

  if (rank == 0) {
    Gammahat <- NULL
    Gamma0hat <- diag(r)
    tmp.M <- eigen(SigmaX)
    mu <- log(mean(Y) / (1 - mean(Y)))
    muh <- matrix(1, n, 1) %*% mu
    c1 <- log(1 + exp(muh))
    Cn1 <- crossprod(Y, muh) - colSums(c1)
    eta <- NULL
    var <- NULL
    weight <- rep(1, n)
    V <- Y
    objfun <- Cn1 - n/2 * sum(log(tmp.M$values))
  } else if (rank == p) {
    Gammahat <- diag(r)
    Gamma0hat <- NULL
    tmp.M <- eigen(SigmaX)
    eta <- as.matrix(beta)
    var <- Tauwx / (- c.theta.mean)
    c1 <- log(1 + exp(theta))
    Cn1 <- crossprod(Y, theta) - colSums(c1)
    objfun <- Cn1 - n/2 * sum(log(tmp.M$values))
  } else if (rank == p - 1) {
    maxiter = 250
    ftol = 0.0005

    M <- SigmaX
    MU <- SigmaX
    tmp.MU <- eigen(MU)
    invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1/tmp.MU$values, "*"), tmp.MU$vectors)
    c1 <- log(1 + exp(theta))
    temp1 <- crossprod(Y, theta) - colSums(c1)
    e1 <- eigen(crossprod(gamma.init, M) %*% gamma.init)
    e2 <- eigen(crossprod(gamma.init, invMU) %*% gamma.init)
    e3 <- eigen(SigmaX)
    temp2 <- sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))
    obj1 <- temp1 - (n / 2) * temp2
    GiX <- X %*% gamma.init
    fit <- cvreg:::.biasreducedglm(cbind(1,GiX), as.vector(Y), family = binomial())
    eta <- as.matrix(fit$coef[2 : (rank + 1)])
    init <- gamma.init
    GEidx <- .GaussElim(init)
    Ginit <- init %*% pseudoinverse(init[GEidx[1:rank], ])
    j <- GEidx[p]
    g <- as.matrix(Ginit[j, ])
    g1 <- Ginit[-j, ]
    t2 <- crossprod(g1, as.matrix(M[-j, j]))/M[j, j]
    t3 <- crossprod(g1, as.matrix(invMU[-j, j]))/invMU[j, j]
    GUGt2 <- g + t2
    GUG <- crossprod(Ginit, (M %*% Ginit)) - tcrossprod(GUGt2, GUGt2) * M[j, j]
    GVGt2 <- g + t3
    GVG <- crossprod(Ginit, (invMU %*% Ginit)) - tcrossprod(GVGt2, GVGt2) * invMU[j, j]
    invC1 <- .ridgeinv(GUG)
    invC2 <- .ridgeinv(GVG)

    X1 <- as.matrix(X[ , j])
    X2 <- as.matrix(X[ , -j])
    t5 <- tcrossprod(X1,g) %*% eta
    t6 <- matrix(1, n, 1) %*% mu + X2 %*% g1 %*% eta
    t7 <- t6 + t5
    t8 <- crossprod(Y, t7)
    et6 <- log(1 + exp(t7))
    t9 <- colSums(et6)

    fobj <- function(x) {
      tmp2 <- x + t2
      tmp3 <- x + t3
      tmp4 <- tcrossprod(X1,x) %*% eta + t6
      T2 <- invC1 %*% tmp2
      T3 <- invC2 %*% tmp3
      T4 <- log(1 + exp(tmp4))
      n /2 * (- 2 * log(1 + sum(x^2)) + log(1 + M[j, j] * crossprod(tmp2,T2)) + log(1 + invMU[j, j] * crossprod(tmp3,T3))) - t(Y) %*% tmp4 + colSums(T4)
    }

    gobj <- function(x) {
      tmp2 <- x + t2
      tmp3 <- x + t3
      tmp4 <- tcrossprod(X1,x) %*% eta + t6
      Cn <- Y - (exp(tmp4) / (1 + exp(tmp4)))
      T2 <- invC1 %*% tmp2
      T3 <- invC2 %*% tmp3
      A1 <- crossprod(g1, Sigmawx[-j, -j]) %*% g1 + as.matrix(x) %*% Sigmawx[j, -j] %*% g1 + tcrossprod(crossprod(g1, Sigmawx[-j, j]), x) + tcrossprod(x) * Sigmawx[j , j]
      invC3 <- .ridgeinv(A1)
      T5 <- tcrossprod(crossprod(X, Cn),eta)
      tmp5 <- tcrossprod(X1,x) + X2 %*% g1
      T6 <- tcrossprod(Sigmawxv, Cn) %*% tmp5 %*% invC3
      tmp6 <- as.matrix(Sigmawx[ , -j]) %*% g1 + tcrossprod(Sigmawx[ , j], x)
      T7 <- tmp6 %*% invC3 %*% tcrossprod(crossprod(tmp5, Cn), eta)
      T8 <- tmp6 %*% eta %*% crossprod(Cn, tmp5) %*% invC3
      r1 <- rep(0, p)
      r1[j] <- 1
      T9 <- crossprod(r1, (T5 + T6 - T7 - T8))
      - t(T9) + n * (- 4 * x %*% solve(1 + sum(x^2)) + 2 * T2/as.numeric(1/M[j, j] + crossprod(tmp2, T2)) + 2 * T3/as.numeric(1/invMU[j, j] + crossprod(tmp3, T3)))
    }

    i <- 1
    while (i < maxiter) {
      res <- nlminb(Ginit[j, ], fobj, gobj)
      Ginit[j, ] <- res$par
      GiX <- X %*% Ginit
      fit <- cvreg:::.biasreducedglm(cbind(1,GX), as.vector(Y), family = binomial())
      mu <- as.vector(fit$coef[1])
      eta <- as.matrix(fit$coef[2 : (rank + 1)])
      beta <- Ginit %*% eta
      theta <- matrix(1, n, 1) %*% mu + X %*% beta
      c.theta <- - exp(theta) / (1 + exp(theta))^2
      c.theta.mean <- mean(c.theta)
      weight <- c.theta / c.theta.mean
      mu1 <- exp(theta) / (1 + exp(theta))
      V <- theta + ((Y - mu1) / weight)
      W <- diag(as.vector(weight))
      wx <- W %*% X
      mean.wx <- apply(wx, 2, mean)
      wxx <- X - matrix(1, nrow = n) %*% mean.wx
      Sigmawx <- crossprod(wxx, W) %*% wxx / n
      wv <- W %*% V
      mean.wv <- mean(wv)
      wvv <- V - matrix(1, nrow = n) %*% mean.wv
      Sigmawxv <- crossprod(wxx, W) %*% wvv / n
      Tauwx <- .ridgeinv(Sigmawx)
      e1 <- eigen(crossprod(Ginit, M) %*% Ginit)
      e2 <- eigen(crossprod(Ginit, invMU) %*% Ginit)
      e3 <- eigen(crossprod(Ginit))
      e4 <- matrix(1, n, 1) %*% mu + X %*% beta
      e5 <- crossprod(Y, e4) - colSums(log(1 + exp(e4)))
      obj2 <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))) + e5
      if (abs(obj1 - obj2) < ftol) {
        break
      }
      else {
        obj1 <- obj2
        i <- i + 1
      }
    }
    a <- qr(Ginit)
    Gammahat <- qr.Q(a)
    GX <- X %*% Gammahat
    fit <- cvreg:::.biasreducedglm(cbind(1,GX), as.vector(Y), family = binomial())
    mu <- as.vector(fit$coef[1])
    eta <- matrix(fit$coef[2 : (rank + 1)])
    beta <- Gammahat %*% eta
    theta <- matrix(1, n, 1) %*% mu + X %*% beta
    c.theta <- - exp(theta) / (1 + exp(theta))^2
    c.theta.mean <- mean(c.theta)
    weight <- c.theta / c.theta.mean
    mu1 <- exp(theta) / (1 + exp(theta))
    V <- theta + ((Y - mu1) / weight)
    W <- diag(as.vector(weight))
    wx <- W %*% X
    mean.wx <- apply(wx, 2, mean)
    wxx <- X - matrix(1, nrow = n) %*% mean.wx
    Sigmawx <- crossprod(wxx, W) %*% wxx / n
    wv <- W %*% V
    mean.wv <- mean(wv)
    wvv <- V - matrix(1, nrow = n) %*% mean.wv
    Sigmawxv <- crossprod(wxx, W) %*% wvv / n
    Tauwx <- .ridgeinv(Sigmawx)
    var <- Tauwx / (- c.theta.mean)
    e1 <- eigen(crossprod(Gammahat, M) %*% Gammahat)
    e2 <- eigen(crossprod(Gammahat, invMU) %*% Gammahat)
    e3 <- eigen(mpd(M))
    e4 <- matrix(1, n, 1) %*% mu + X %*% Gammahat %*% eta
    e5 <- crossprod(Y, e4) - colSums(log(1 + exp(e4)))
    Gamma0hat <- qr.Q(a, complete = TRUE)[, (rank + 1):p, drop = FALSE]
    objfun <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))) + e5
    Gammahat <- as.matrix(Gammahat)
    Gamma0hat <- as.matrix(Gamma0hat)
  } else {
    maxiter = 250
    ftol = 0.0005
    M <- SigmaX
    MU <- SigmaX
    tmp.MU <- eigen(MU)
    invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1/tmp.MU$values, "*"), tmp.MU$vectors)
    c1 <- log(1 + exp(theta))
    temp1 <- crossprod(Y, theta) - colSums(c1)
    e1 <- eigen(crossprod(gamma.init, M) %*% gamma.init)
    e2 <- eigen(crossprod(gamma.init, invMU) %*% gamma.init)
    e3 <- eigen(SigmaX)
    temp2 <- sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))
    obj1 <- temp1 - (n / 2) * temp2
    options(warn = -1)
    GiX <- X %*% gamma.init
    fit <- cvreg:::.biasreducedglm(cbind(1,GiX), as.vector(Y), family = binomial())
    eta <- as.matrix(fit$coef[2 : (rank + 1)])
    init <- gamma.init
    GEidx <- .GaussElim(init)
    Ginit <- init %*% pseudoinverse(init[GEidx[1:rank], ])
    GUG <- crossprod(Ginit, (M %*% Ginit))
    GVG <- crossprod(Ginit, (invMU %*% Ginit))
    t4 <- crossprod(Ginit[GEidx[(rank + 1):p], ], Ginit[GEidx[(rank + 1):p], ]) + diag(rank)
    i <- 1
    while (i < maxiter) {
      for (j in GEidx[(rank + 1):p]) {
        g <- as.matrix(Ginit[j, ])
        g1 <- as.matrix(Ginit[-j, ])
        t2 <- crossprod(g1, as.matrix(M[-j, j]))/M[j, j]
        t3 <- crossprod(g1, as.matrix(invMU[-j, j]))/invMU[j, j]
        GUGt2 <- g + t2
        GUG <- GUG - tcrossprod(GUGt2, GUGt2) * M[j, j]
        GVGt2 <- g + t3
        GVG <- GVG - tcrossprod(GVGt2, GVGt2) * invMU[j, j]
        t4 <- t4 - tcrossprod(g, g)
        invC1 <- .ridgeinv(GUG)
        invC2 <- .ridgeinv(GVG)
        invt4 <- .ridgeinv(t4)
        X1 <- X[ , j]
        X2 <- X[ , -j]
        t5 <- tcrossprod(X1,g) %*% eta
        t6 <- matrix(1, n, 1) %*% mu + X2 %*% g1 %*% eta
        t7 <- t5 + t6
        t8 <- crossprod(Y, t7)
        et6 <- log(1 + exp(t7))
        t9 <- colSums(et6)

        fobj <- function(x) {
          tmp2 <- x + t2
          tmp3 <- x + t3
          tmp4 <- tcrossprod(X1, x) %*% eta + t6
          T1 <- invt4 %*% x
          T2 <- invC1 %*% tmp2
          T3 <- invC2 %*% tmp3
          T4 <- log(1 + exp(tmp4))
          n /2 * (-2 * log(1 + x %*% T1) + log(1 + M[j, j] * crossprod(tmp2, T2)) + log(1 + invMU[j, j] *crossprod(tmp3, T3))) - t(Y) %*% tmp4 + colSums(T4)
        }

        gobj <- function(x) {
          tmp2 <- x + t2
          tmp3 <- x + t3
          tmp4 <- tcrossprod(X1, x) %*% eta + t6
          Cn <- Y - (exp(tmp4) / (1 + exp(tmp4)))
          T1 <- invt4 %*% x
          T2 <- invC1 %*% tmp2
          T3 <- invC2 %*% tmp3
          A1 <- crossprod(g1, Sigmawx[-j, -j]) %*% g1 + as.matrix(x) %*% Sigmawx[j, -j] %*% g1 + tcrossprod(crossprod(g1, Sigmawx[-j, j]) , x) + tcrossprod(x) * Sigmawx[j , j]
          invC3 <- .ridgeinv(A1)
          T5 <-  tcrossprod(crossprod(X, Cn), eta)
          tmp5 <- tcrossprod(X1, x) + X2 %*% g1
          T6 <- tcrossprod(Sigmawxv, Cn) %*% tmp5 %*% invC3
          tmp6 <- Sigmawx[ , -j] %*% g1 + tcrossprod(Sigmawx[ , j], x)
          T7 <- tmp6 %*% invC3 %*% tcrossprod(crossprod(tmp5, Cn), eta)
          T8 <- tmp6 %*% eta %*% crossprod(Cn, tmp5) %*% invC3
          r1 <- rep(0, p)
          r1[j] <- 1
          T9 <- crossprod(r1, (T5 + T6 - T7 - T8))
          - t(T9) + n * (- 4 * T1/as.numeric(1 + x %*% T1) + 2 * T2/as.numeric(1/M[j, j] + crossprod(tmp2, T2)) + 2 * T3/as.numeric(1/invMU[j, j] + crossprod(tmp3, T3)))
        }
        res <- nlminb(Ginit[j, ], fobj, gobj)
        Ginit[j, ] <- res$par
        g <- as.matrix(Ginit[j, ])
        t4 <- t4 + tcrossprod(g, g)
        GUGt2 <- g + t2
        GUG <- GUG + tcrossprod(GUGt2, GUGt2) * M[j, j]
        GVGt2 <- g + t3
        GVG <- GVG + tcrossprod(GVGt2, GVGt2) * invMU[j, j]
      }

      GiX <- X %*% Ginit
      fit <- cvreg:::.biasreducedglm(cbind(1,GiX), as.vector(Y), family = binomial())
      mu <- as.vector(fit$coef[1])
      eta <- as.matrix(fit$coef[2 : (rank + 1)])
      beta <- Ginit %*% eta
      theta <- matrix(1, n, 1) %*% mu + X %*% beta
      c.theta <- - exp(theta) / (1 + exp(theta))^2
      c.theta.mean <- mean(c.theta)
      weight <- c.theta / c.theta.mean
      mu1 <- exp(theta) / (1 + exp(theta))
      V <- theta + ((Y - mu1) / weight)
      W <- diag(as.vector(weight))
      wx <- W %*% X
      mean.wx <- apply(wx, 2, mean)
      wxx <- X - matrix(1, nrow = n) %*% mean.wx
      Sigmawx <- crossprod(wxx, W) %*% wxx / n
      wv <- W %*% V
      mean.wv <- mean(wv)
      wvv <- V - matrix(1, nrow = n) %*% mean.wv
      Sigmawxv <- crossprod(wxx, W) %*% wvv / n
      Tauwx <- .ridgeinv(Sigmawx)
      e1 <- eigen(crossprod(Ginit, M) %*% Ginit)
      e2 <- eigen(crossprod(Ginit, invMU) %*% Ginit)
      e3 <- eigen(crossprod(Ginit))
      e4 <- matrix(1, n, 1) %*% mu + X %*% beta
      e5 <- crossprod(Y, e4) - colSums(log(1 + exp(e4)))
      obj2 <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))) + e5
      if (abs(obj1 - obj2) < ftol) {
        break
      }
      else {
        obj1 <- obj2
        i <- i + 1
      }
    }
    a <- qr(Ginit)
    Gammahat <- qr.Q(a)
    GX <- X %*% Gammahat
    fit <- cvreg:::.biasreducedglm(cbind(1,GX), as.vector(Y), family = binomial())
    mu <- as.vector(fit$coef[1])
    eta <- matrix(fit$coef[2 : (rank + 1)])
    beta <- Gammahat %*% eta
    theta <- matrix(1, n, 1) %*% mu + X %*% beta
    c.theta <- - exp(theta) / (1 + exp(theta))^2
    c.theta.mean <- mean(c.theta)
    weight <- c.theta / c.theta.mean
    mu1 <- exp(theta) / (1 + exp(theta))
    V <- theta + ((Y - mu1) / weight)
    W <- diag(as.vector(weight))
    wx <- W %*% X
    mean.wx <- apply(wx, 2, mean)
    wxx <- X - matrix(1, nrow = n) %*% mean.wx
    Sigmawx <- crossprod(wxx, W) %*% wxx / n
    wv <- W %*% V
    mean.wv <- mean(wv)
    wvv <- V - matrix(1, nrow = n) %*% mean.wv
    Sigmawxv <- crossprod(wxx, W) %*% wvv / n
    Tauwx <- .ridgeinv(Sigmawx)
    var <- Tauwx / (- c.theta.mean)
    e1 <- eigen(crossprod(Gammahat, M) %*% Gammahat)
    e2 <- eigen(crossprod(Gammahat, invMU) %*% Gammahat)
    e3 <- eigen(mpd(M))
    e4 <- matrix(1, n, 1) %*% mu + X %*% Gammahat %*% eta
    e5 <- crossprod(Y, e4) - colSums(log(1 + exp(e4)))
    Gamma0hat <- qr.Q(a, complete = TRUE)[, (rank + 1):p]
    objfun <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))) + e5
    Gammahat <- as.matrix(Gammahat)
    Gamma0hat <- as.matrix(Gamma0hat)

  }
  return(list(Gammahat = Gammahat, Gamma0hat = Gamma0hat, muhat = mu,
              etahat = eta, weighthat = weight, Vhat = V,
              avar = var, objfun = objfun))
}


poisson.envlp_MU <- function(X, Y, rank){

  X <- as.matrix(X)
  Y <- as.matrix(Y)
  a <- dim(Y)
  n <- a[1]
  r <- a[2]
  p <- ncol(X)
  if (a[1] != nrow(X))
    stop("X and Y should have the same number of observations.")
  if (rank > p || rank < 0)
    stop("rank must be an integer between 0 and r.")
  if (sum(duplicated(cbind(X, Y), MARGIN = 2)) > 0)
    stop("Some responses also appear in the predictors, or there maybe duplicated columns in X or Y.")

  options(warn = -1)
  fit <- cvreg:::.biasreducedglm(cbind(1,X), as.vector(Y), family = poisson())
  mu <- as.vector(fit$coef[1])
  beta <- as.vector(fit$coef[2 : (p + 1)])
  theta <- matrix(1, n, 1) %*% mu + X %*% beta
  c.theta <- -exp(theta)
  c.theta.mean <- mean(c.theta)
  weight <- c.theta / c.theta.mean
  V <- theta + ((Y - c.theta) / weight)
  W <- diag(as.vector(weight))
  wx <- W %*% X
  mean.wx <- apply(wx, 2, mean)
  wxx <- X - matrix(1, nrow = n) %*% mean.wx
  Sigmawx <- crossprod(wxx, W) %*% wxx / n
  wv <- W %*% V
  mean.wv <- mean(wv)
  wvv <- V - matrix(1, nrow = n) %*% mean.wv
  Sigmawxv <- crossprod(wxx, W) %*% wvv / n
  Tauwx <- .ridgeinv(Sigmawx)

  M.init <- Tauwx / (- c.theta.mean)
  U.init <- tcrossprod(beta)
  tmp1 <- envlp_MU(M.init, U.init, rank)
  gamma.init <- tmp1$Gammahat
  gamma0.init <- tmp1$Gamma0hat
  SigmaX <- solve(.ridgeinv(stats::cov(X))) * ((n-1)/n)
  TauX <- .ridgeinv(SigmaX)

  if (rank == 0) {
    Gammahat <- NULL
    Gamma0hat <- diag(r)
    tmp.M <- eigen(SigmaX)
    mu <- log(mean(Y))
    muh <- matrix(1, n, 1) %*% mu
    Cn1 <- crossprod(Y, muh) - colSums(exp(muh))
    eta <- NULL
    var <- NULL
    weight <- rep(1, n)
    V <- Y
    objfun <- Cn1 - n/2 * sum(log(tmp.M$values))
  } else if (rank == p) {
    Gammahat <- diag(r)
    Gamma0hat <- NULL
    tmp.M <- eigen(SigmaX)
    eta <- beta
    var <- Tauwx / (- c.theta.mean)
    Cn1 <- crossprod(Y, theta) - colSums(exp(theta))
    objfun <- Cn1 - n/2 * sum(log(tmp.M$values))
  } else if (rank == p - 1) {
    maxiter = 250
    ftol = 0.0005

    M <- SigmaX
    MU <- SigmaX
    tmp.MU <- eigen(MU)
    invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1/tmp.MU$values, "*") ,tmp.MU$vectors)

    Cn1 <- crossprod(Y, theta) - colSums(exp(theta))
    e1 <- eigen(crossprod(gamma.init, M) %*% gamma.init)
    e2 <- eigen(crossprod(gamma.init, invMU) %*% gamma.init)
    e3 <- eigen(SigmaX)
    temp1 <- sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))
    obj1 <- Cn1 - (n / 2) * temp1

    GiX <- X %*% gamma.init
    fit <- cvreg:::.biasreducedglm(cbind(1,GiX), as.vector(Y), family = poisson())
    eta <- as.matrix(fit$coef[2 : (rank + 1)])

    init <- gamma.init
    GEidx <- .GaussElim(init)
    Ginit <- init %*% pseudoinverse(init[GEidx[1:rank], ])
    j <- GEidx[p]
    g <- as.matrix(Ginit[j, ])
    g1 <- Ginit[-j, ]
    t2 <- crossprod(g1, as.matrix(M[-j, j]))/M[j, j]
    t3 <- crossprod(g1, as.matrix(invMU[-j, j]))/invMU[j, j]
    GUGt2 <- g + t2
    GUG <- crossprod(Ginit, (M %*% Ginit)) - tcrossprod(GUGt2, GUGt2) * M[j, j]
    GVGt2 <- g + t3
    GVG <- crossprod(Ginit, (invMU %*% Ginit)) - tcrossprod(GVGt2, GVGt2) * invMU[j, j]
    invC1 <- .ridgeinv(GUG)
    invC2 <- .ridgeinv(GVG)

    X1 <- X[ , j]
    X2 <- X[ , -j]
    t5 <- tcrossprod(X1,g) %*% eta
    t6 <- matrix(1, n, 1) %*% mu + X2 %*% g1 %*% eta
    t7 <- t6 + t5
    t8 <- crossprod(Y, t7)
    et6 <- exp(t7)
    t9 <- colSums(et6)

    fobj <- function(x) {
      tmp2 <- x + t2
      tmp3 <- x + t3
      tmp4 <- tcrossprod(X1, x) %*% eta + t6
      T2 <- invC1 %*% tmp2
      T3 <- invC2 %*% tmp3
      T4 <- exp(tmp4)
      n /2 * (- 2 * log(1 + sum(x^2)) + log(1 + M[j, j] * crossprod(tmp2,T2)) + log(1 + invMU[j, j] * crossprod(tmp3,T3))) - t(Y) %*% tmp4 + colSums(T4)
    }

    gobj <- function(x) {
      tmp2 <- x + t2
      tmp3 <- x + t3
      tmp4 <- tcrossprod(X1, x) %*% eta + t6
      Cn <- Y - exp(tmp4)
      T2 <- invC1 %*% tmp2
      T3 <- invC2 %*% tmp3
      A1 <- crossprod(g1, Sigmawx[-j, -j]) %*% g1 + as.matrix(x) %*% Sigmawx[j, -j] %*% g1 + tcrossprod(crossprod(g1, Sigmawx[-j, j]),x) + tcrossprod(x) * Sigmawx[j , j]
      invC3 <- .ridgeinv(A1)
      T5 <-  tcrossprod(crossprod(X, Cn), eta)
      tmp5 <- tcrossprod(X1, x) + X2 %*% g1
      T6 <- tcrossprod(Sigmawxv, Cn) %*% tmp5 %*% invC3
      tmp6 <- as.matrix(Sigmawx[ , -j]) %*% g1 + tcrossprod(Sigmawx[ , j], x)
      T7 <- tmp6 %*% invC3 %*% tcrossprod(crossprod(tmp5, Cn), eta)
      T8 <- tmp6 %*% eta %*% crossprod(Cn, tmp5) %*% invC3
      r1 <- rep(0, p)
      r1[j] <- 1
      T9 <- crossprod(r1, (T5 + T6 - T7 - T8))
      - t(T9) + n * (- 4 * x %*% .ridgeinv(1 + sum(x^2)) + 2 * T2/as.numeric(1/M[j, j] + crossprod(tmp2, T2)) + 2 * T3/as.numeric(1/invMU[j, j] + crossprod(tmp3, T3)))
    }

    i <- 1
    while (i < maxiter) {
      res <- nlminb(Ginit[j, ], fobj, gobj)
      Ginit[j, ] <- res$par

      GX <- X %*% Ginit
      fit <- cvreg:::.biasreducedglm(cbind(1,GX), as.vector(Y), family = poisson())
      mu <- as.vector(fit$coef[1])
      eta <- as.matrix(fit$coef[2 : (rank + 1)])
      beta <- Ginit %*% eta
      theta <- matrix(1, n, 1) %*% mu + X %*% beta
      c.theta <- - exp(theta)
      c.theta.mean <- mean(c.theta)
      weight <- c.theta / c.theta.mean
      V <- theta + ((Y - c.theta) / weight)
      W <- diag(as.vector(weight))
      wx <- W %*% X
      mean.wx <- apply(wx, 2, mean)
      wxx <- X - matrix(1, nrow = n) %*% mean.wx
      Sigmawx <- crossprod(wxx, W) %*% wxx / n
      wv <- W %*% V
      mean.wv <- mean(wv)
      wvv <- V - matrix(1, nrow = n) %*% mean.wv
      Sigmawxv <- crossprod(wxx, W) %*% wvv / n
      Tauwx <- .ridgeinv(Sigmawx)
      e1 <- eigen(crossprod(Ginit, M) %*% Ginit)
      e2 <- eigen(crossprod(Ginit, invMU) %*% Ginit)
      e3 <- eigen(crossprod(Ginit))
      e4 <- matrix(1, n, 1) %*% mu + X %*% beta
      e5 <- crossprod(Y, e4) - colSums(exp(e4))
      obj2 <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))) + e5
      if (abs(obj1 - obj2) < ftol) {
        break
      }
      else {
        obj1 <- obj2
        i <- i + 1
      }
    }
    a <- qr(Ginit)
    Gammahat <- qr.Q(a)
    GX <- X %*% Gammahat
    fit <- cvreg:::.biasreducedglm(cbind(1,X), as.vector(Y), family = poisson())
    mu <- as.vector(fit$coef[1])
    eta <- matrix(fit$coef[2 : (rank + 1)])
    beta <- Gammahat %*% eta
    theta <- matrix(1, n, 1) %*% mu + X %*% beta
    theta <- matrix(1, n, 1) %*% mu + X %*% beta
    c.theta <- - exp(theta)
    c.theta.mean <- mean(c.theta)
    weight <- c.theta / c.theta.mean
    V <- theta + ((Y - c.theta) / weight)
    W <- diag(as.vector(weight))
    wx <- W %*% X
    mean.wx <- apply(wx, 2, mean)
    wxx <- X - matrix(1, nrow = n) %*% mean.wx
    Sigmawx <- crossprod(wxx, W) %*% wxx / n
    wv <- W %*% V
    mean.wv <- mean(wv)
    wvv <- V - matrix(1, nrow = n) %*% mean.wv
    Sigmawxv <- crossprod(wxx, W) %*% wvv / n
    Tauwx <- .ridgeinv(Sigmawx)
    var <- Tauwx / (- c.theta.mean)
    e1 <- eigen(crossprod(Gammahat, M) %*% Gammahat)
    e2 <- eigen(crossprod(Gammahat, invMU) %*% Gammahat)
    e3 <- eigen(mpd(M))
    e4 <- matrix(1, n, 1) %*% mu + X %*% Gammahat %*% eta
    e5 <- crossprod(Y, e4) - colSums(exp(e4))
    Gamma0hat <- qr.Q(a, complete = TRUE)[, (rank + 1):r, drop = FALSE]
    objfun <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))) + e5
    Gammahat <- as.matrix(Gammahat)
    Gamma0hat <- as.matrix(Gamma0hat)

  } else {
    maxiter = 250
    ftol = 0.0005

    M <- SigmaX
    MU <- SigmaX
    tmp.MU <- eigen(MU)
    invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1/tmp.MU$values, "*") ,tmp.MU$vectors)

    Cn1 <- crossprod(Y, theta) - colSums(exp(theta))
    e1 <- eigen(crossprod(gamma.init, M) %*% gamma.init)
    e2 <- eigen(crossprod(gamma.init, invMU) %*% gamma.init)
    e3 <- eigen(SigmaX)
    temp1 <- sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))
    obj1 <- Cn1 - (n / 2) * temp1

    GiX <- X %*% gamma.init
    fit <- cvreg:::.biasreducedglm(cbind(1,GiX), as.vector(Y), family = poisson())
    eta <- as.matrix(fit$coef[2 : (rank + 1)])

    init <- gamma.init
    GEidx <- .GaussElim(init)
    Ginit <- init %*% pseudoinverse(init[GEidx[1:rank], ])
    GUG <- crossprod(Ginit, (M %*% Ginit))
    GVG <- crossprod(Ginit, (invMU %*% Ginit))
    t4 <- crossprod(Ginit[GEidx[(rank + 1):p], ], Ginit[GEidx[(rank + 1):p], ]) + diag(rank)

    i <- 1
    while (i < maxiter) {
      for (j in GEidx[(rank + 1):p]) {
        g <- as.matrix(Ginit[j, ])
        g1 <- as.matrix(Ginit[-j, ])
        t2 <- crossprod(g1, as.matrix(M[-j, j]))/M[j, j]
        t3 <- crossprod(g1, as.matrix(invMU[-j, j]))/invMU[j, j]
        GUGt2 <- g + t2
        GUG <- GUG - tcrossprod(GUGt2, GUGt2) * M[j, j]
        GVGt2 <- g + t3
        GVG <- GVG - tcrossprod(GVGt2, GVGt2) * invMU[j, j]
        t4 <- t4 - tcrossprod(g)
        invC1 <-.ridgeinv(GUG)
        invC2 <-.ridgeinv(GVG)
        invt4 <-.ridgeinv(t4)

        X1 <- X[ , j]
        X2 <- X[ , -j]
        t5 <- tcrossprod(X1,g) %*% eta
        t6 <- matrix(1, n, 1) %*% mu + X2 %*% g1 %*% eta
        t7 <- t5 + t6
        t8 <- crossprod(Y, t7)
        et7 <- exp(t7)
        t9 <- colSums(et7)

        fobj <- function(x) {
          tmp2 <- x + t2
          tmp3 <- x + t3
          tmp4 <- tcrossprod(X1, x) %*% eta + t6
          T1 <- invt4 %*% x
          T2 <- invC1 %*% tmp2
          T3 <- invC2 %*% tmp3
          T4 <- exp(tmp4)
          n /2 * (-2 * log(1 + x %*% T1) + log(1 + M[j, j] * crossprod(tmp2, T2)) + log(1 + invMU[j, j] *crossprod(tmp3, T3))) - t(Y) %*% tmp4 + colSums(T4)
        }

        gobj <- function(x) {
          tmp2 <- x + t2
          tmp3 <- x + t3
          tmp4 <- tcrossprod(X1, x) %*% eta + t6
          Cn <- Y - exp(tmp4)
          T1 <- invt4 %*% x
          T2 <- invC1 %*% tmp2
          T3 <- invC2 %*% tmp3
          A1 <- crossprod(g1, Sigmawx[-j, -j]) %*% g1 + as.matrix(x) %*% Sigmawx[j, -j] %*% g1 + tcrossprod(crossprod(g1, Sigmawx[-j, j]), x) + tcrossprod(x) * Sigmawx[j , j]
          invC3 <-.ridgeinv(A1)
          T5 <-  tcrossprod(crossprod(X, Cn), eta)
          tmp5 <- tcrossprod(X1, x) + X2 %*% g1
          T6 <- tcrossprod(Sigmawxv, Cn) %*% tmp5 %*% invC3
          tmp6 <- as.matrix(Sigmawx[ , -j]) %*% g1 + tcrossprod(Sigmawx[ , j], x)
          T7 <- tmp6 %*% invC3 %*% tcrossprod(crossprod(tmp5, Cn), eta)
          T8 <- tmp6 %*% eta %*% crossprod(Cn, tmp5) %*% invC3
          r1 <- rep(0, p)
          r1[j] <- 1
          T9 <- crossprod(r1, (T5 + T6 - T7 - T8))
          - t(T9) + n * (- 4 * T1/as.numeric(1 + x %*% T1) + 2 * T2/as.numeric(1/M[j, j] + crossprod(tmp2, T2)) + 2 * T3/as.numeric(1/invMU[j, j] + crossprod(tmp3, T3)))
        }
        res <- nlminb(Ginit[j, ], fobj, gobj)
        Ginit[j, ] <- res$par
        g <- as.matrix(Ginit[j, ])
        t4 <- t4 + tcrossprod(g, g)
        GUGt2 <- g + t2
        GUG <- GUG + tcrossprod(GUGt2, GUGt2) * M[j, j]
        GVGt2 <- g + t3
        GVG <- GVG + tcrossprod(GVGt2, GVGt2) * invMU[j, j]
      }

      GX <- X %*% Ginit
      fit <-  cvreg:::.biasreducedglm(cbind(1,GX), as.vector(Y), family = poisson())
      mu <- as.vector(fit$coef[1])
      eta <- as.matrix(fit$coef[2 : (rank + 1)])
      beta <- Ginit %*% eta
      theta <- matrix(1, n, 1) %*% mu + X %*% beta
      c.theta <- - exp(theta)
      c.theta.mean <- mean(c.theta)
      weight <- c.theta / c.theta.mean
      V <- theta + ((Y - c.theta) / weight)
      W <- diag(as.vector(weight))
      wx <- W %*% X
      mean.wx <- apply(wx, 2, mean)
      wxx <- X - matrix(1, nrow = n) %*% mean.wx
      Sigmawx <- crossprod(wxx, W) %*% wxx / n
      wv <- W %*% V
      mean.wv <- mean(wv)
      wvv <- V - matrix(1, nrow = n) %*% mean.wv
      Sigmawxv <- crossprod(wxx, W) %*% wvv / n
      Tauwx <-.ridgeinv(Sigmawx)

      e1 <- eigen(crossprod(Ginit, M) %*% Ginit)
      e2 <- eigen(crossprod(Ginit, invMU) %*% Ginit)
      e3 <- eigen(crossprod(Ginit))
      e4 <- matrix(1, n, 1) %*% mu + X %*% beta
      e5 <- crossprod(Y, e4) - colSums(exp(e4))
      obj2 <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))) + e5
      if (abs(obj1 - obj2) < ftol) {
        break
      }
      else {
        obj1 <- obj2
        i <- i + 1
      }
    }
  }
  a <- qr(Ginit)
  Gammahat <- qr.Q(a)
  GX <- X %*% Gammahat
  fit <-  cvreg:::.biasreducedglm(cbind(1,GX), as.vector(Y), family = poisson())
  mu <- as.vector(fit$coef[1])
  eta <- matrix(fit$coef[2 : (rank + 1)])
  beta <- Gammahat %*% eta
  theta <- matrix(1, n, 1) %*% mu + X %*% beta
  c.theta <- - exp(theta)
  c.theta.mean <- mean(c.theta)
  weight <- c.theta / c.theta.mean
  V <- theta + ((Y - c.theta) / weight)
  W <- diag(as.vector(weight))
  wx <- W %*% X
  mean.wx <- apply(wx, 2, mean)
  wxx <- X - matrix(1, nrow = n) %*% mean.wx
  Sigmawx <- crossprod(wxx, W) %*% wxx / n
  wv <- W %*% V
  mean.wv <- mean(wv)
  wvv <- V - matrix(1, nrow = n) %*% mean.wv
  Sigmawxv <- crossprod(wxx, W) %*% wvv / n
  Tauwx <- .ridgeinv(Sigmawx)
  var <- Tauwx / (- c.theta.mean)
  e1 <- eigen(crossprod(Gammahat, M) %*% Gammahat)
  e2 <- eigen(crossprod(Gammahat, invMU) %*% Gammahat)
  e3 <- eigen(mpd(M))
  e4 <- matrix(1, n, 1) %*% mu + X %*% Gammahat %*% eta
  e5 <- crossprod(Y, e4) - colSums(exp(e4))
  Gamma0hat <- qr.Q(a, complete = TRUE)[, (rank + 1):p]
  objfun <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))) + e5
  Gammahat <- as.matrix(Gammahat)
  Gamma0hat <- as.matrix(Gamma0hat)

  return(list(Gammahat = Gammahat, Gamma0hat = Gamma0hat, muhat = mu,
              etahat = eta, weighthat = weight, Vhat = V,
              avar = var, objfun = objfun))
}


henvlp_MU<- function(M, U, MU, rank, n, ng, L){
  p <- L
  dimM <- dim(M[[1]])
  dimU <- dim(U[[1]])
  r <- dimM[1]
  auxinit <- function(M, U, MU, rank, ng, n, p){
    tmp.MU <- eigen(MU)
    invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1 / tmp.MU$values, '*'), tmp.MU$vectors)
    invMU2 <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1 / sqrt(tmp.MU$values), '*') , tmp.MU$vectors)
    midmatrix <- invMU2
    startv4 <- function(a) crossprod(a,midmatrix) %*% a
    tmp2.MU <- apply(tmp.MU$vectors, 2, startv4)
    tmp3.MU <- sort(tmp2.MU, decreasing = TRUE, index.return = TRUE)
    init <- as.matrix(tmp.MU$vectors[, tmp3.MU$ix[1:rank]])
    obj1a <- c()
    for(i in 1:p){
      eig1 <- eigen(crossprod(init, M[[i]]) %*% init)
      eig2 <- eigen(crossprod(init, invMU) %*% init)
      obj1a[i] <- sum(log(eig1$values))*ng[i]/n + sum(log(eig2$values))
    }
    obj1 <- sum(obj1a)
    return(list(init = init, obj1 = obj1, invMU = invMU))
  }
  auxf1 <- function(M1, U1, rank, n, ng, p, init, x, r){
    M <- M1
    U <- U1
    MU <- M + U
    tmp.MU <- eigen(MU)
    invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1 / tmp.MU$values, '*') , tmp.MU$vectors)
    GEidx <- .GaussElim(init)
    Ginit <- init %*% pseudoinverse(init[GEidx[1:rank], ])
    j <- GEidx[r]
    g <- as.matrix(Ginit[j, ])
    t2 <- crossprod(Ginit[-j, ], as.matrix(M[-j, j])) / M[j, j]
    t3 <- crossprod(Ginit[-j, ], as.matrix(invMU[-j, j])) / invMU[j, j]
    GUGt2 <- g + t2
    GUG <- crossprod(Ginit, (M %*% Ginit)) - tcrossprod(GUGt2, GUGt2) * M[j, j]
    GVGt2 <- g + t3
    GVG <- crossprod(Ginit, (invMU %*% Ginit)) - tcrossprod(GVGt2, GVGt2) * invMU[j, j]
    invc1 <- .ridgeinv(GUG)
    invc2 <- .ridgeinv(GVG)
    tmp2 <- x + t2
    tmp3 <- x + t3
    T2 <- invc1 %*% tmp2
    T3 <- invc2 %*% tmp3
    out <- ng * log(1 + M[j, j] * crossprod(tmp2, T2))/n + log(1 + invMU[j, j] * crossprod(tmp3, T3))/p
    return(out)
  }

  auxf2 <- function(M1, U1, t2, t3, invc1, invc2, ng, n, p, x, j){
    M <- M1
    U <- U1
    MU <- M + U
    tmp.MU <- eigen(MU)
    invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1 / tmp.MU$values, '*') , tmp.MU$vectors)
    tmp2 <- x + t2
    tmp3 <- x + t3
    T2 <- invc1 %*% tmp2
    T3 <- invc2 %*% tmp3
    out <- ng * log(1 + M[j, j] * crossprod(tmp2, T2))/n + log(1 + invMU[j, j] * crossprod(tmp3, T3))/p
    return(out)
  }

  auxg1 <- function(M1, U1, rank, n, ng, p, init, x, r){
    M <- M1
    U <- U1
    MU <- M + U
    tmp.MU <- eigen(MU)
    invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1 / tmp.MU$values, '*') , tmp.MU$vectors)
    GEidx <- .GaussElim(init)
    Ginit <- init %*% pseudoinverse(init[GEidx[1 : rank], ])
    j <- GEidx[r]
    g <- as.matrix(Ginit[j, ])
    t2 <- crossprod(Ginit[-j, ], as.matrix(M[-j, j])) / M[j, j]
    t3 <- crossprod(Ginit[-j, ], as.matrix(invMU[-j, j])) / invMU[j, j]
    GUGt2 <- g + t2
    GUG <- crossprod(Ginit, (M %*% Ginit)) - tcrossprod(GUGt2, GUGt2) * M[j, j]
    GVGt2 <- g + t3
    GVG <- crossprod(Ginit, (invMU %*% Ginit)) - tcrossprod(GVGt2, GVGt2) * invMU[j, j]
    invC1 <- .ridgeinv(GUG)
    invC2 <- .ridgeinv(GVG)
    tmp2 <- x + t2
    tmp3 <- x + t3
    T2 <- invC1 %*% tmp2
    T3 <- invC2 %*% tmp3
    out <-  2 * ng * T2 / (n * as.numeric(1 / M[j, j] + crossprod(tmp2, T2))) + 2 * T3 /(p * as.numeric(1 / invMU[j, j] + crossprod(tmp3, T3)))
    return(out)
  }

  auxg2 <- function(M1, U1, t2, t3, invc1, invc2, ng, n, p, x, j){
    M <- M1
    U <- U1
    MU <- M + U
    tmp.MU <- eigen(MU)
    invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2,  1 / tmp.MU$values, '*') , tmp.MU$vectors)
    tmp2 <- x + t2
    tmp3 <- x + t3
    T2 <- invc1 %*% tmp2
    T3 <- invc2 %*% tmp3
    out <-  2 * T2 * ng / (n * as.numeric(1 / M[j, j] + crossprod(tmp2, T2))) + 2 * T3 / (p * as.numeric(1 / invMU[j, j] + crossprod(tmp3, T3)))
    return(out)

  }
  if(rank == 0){
    Gammahat <- NULL
    Gamma0hat <- diag(r)
  }else if (rank == r){
    Gammahat <- diag(r)
    Gamma0hat <- NULL
  }else if (rank == r - 1){
    maxiter = 250
    ftol = 0.001
    initout <- auxinit(M, U, MU, rank, ng, n, p)
    init <- initout$init
    obj1 <- initout$obj1
    invMU <- initout$invMU
    GEidx <- .GaussElim(init)
    Ginit <- init %*% pseudoinverse(init[GEidx[1 : rank], ])
    j <- GEidx[r]
    fobj <- function(x) {
      res <- -2 * log(1 + sum(x^2))
      for(i in 1:p){
        res <- res + auxf1(M[[i]], U[[i]], rank, n, ng[i], p, init, x, r)
      }
      return(res)
    }
    gobj <- function(x) {
      res <- -4 * x %*% solve(1 + sum(x^2))
      for(i in 1:p){
        res <- res + auxg1(M[[i]], U[[i]], rank, n, ng[i], p, init, x, r)
      }
      return(res)
    }
    i <- 1
    while (i < maxiter) {
      if (rank == 1){
        res <- nlminb(Ginit[j,], fn = fobj, gr = gobj)
      } else {
        res <- nlminb(Ginit[j,], fn = fobj, gr = gobj)
      }
      Ginit[j, ] <- res$par
      a <- qr(Ginit)
      Gammahat <- qr.Q(a)
      obj5a <- c()
      for(i in 1:p){
        e1 <- eigen(crossprod(Gammahat, (M[[i]] %*% Gammahat)))
        e2 <- eigen(crossprod(Gammahat, (invMU %*% Gammahat)))
        obj5a[i] <- sum(log(e1$values))*ng[i]/n + sum(log(e2$values))/p
      }
      obj5 <- sum(obj5a)
      if (abs(obj1 - obj5) < ftol * abs(obj1)) {
        break
      } else {
        obj1 <- obj5
        i <- i + 1
      }
    }
    Gamma0hat <- qr.Q(a, complete = TRUE)[, (rank + 1) : r, drop = FALSE]
    Gammahat <- as.matrix(Gammahat)
    Gamma0hat <- as.matrix(Gamma0hat)
  }else{
    maxiter <- 250
    ftol <- 0.001
    initout <- auxinit(M, U, MU, rank, ng, n, p)
    init <- initout$init
    obj1 <- initout$obj1
    GEidx <- .GaussElim(init)
    Ginit <- init %*% pseudoinverse(init[GEidx[1 : rank], ])
    GUG = GVG <- list(length = p)
    for (k in 1 : p){
      MU <- M[[k]] + U[[k]]
      tmp.MU <- eigen(MU)
      invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1 / tmp.MU$values, '*') , tmp.MU$vectors)
      GUG[[k]] <- crossprod(Ginit, (M[[k]] %*% Ginit))
      GVG[[k]] <- crossprod(Ginit, (invMU %*% Ginit))
    }
    t4 <- crossprod(Ginit[GEidx[(rank + 1):r],], Ginit[GEidx[(rank + 1):r], ]) + diag(rank)
    i <- 1
    while (i < maxiter) {
      for (j in GEidx[(rank+1):r]) {
        g <- as.matrix(Ginit[j, ])
        t4 <- t4 - tcrossprod(g, g)
        invt4 <- .ridgeinv(t4)
        t2 = t3 = GUGt2 = GVGt2 = invc1 = invc2 <- list(length=p)
        for(k in 1:p){
          Maux <- M[[k]]
          MU <- M[[k]] + U[[k]]
          tmp.MU <- eigen(MU)
          invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1 / tmp.MU$values, '*'),tmp.MU$vectors)
          t2[[k]] <- crossprod(Ginit[-j, ], as.matrix(Maux[-j, j])) / Maux[j, j]
          t3[[k]] <- crossprod(Ginit[-j, ], as.matrix(invMU[-j, j])) / invMU[j, j]
          GUGt2[[k]] <- g + t2[[k]]
          GUG[[k]] <- GUG[[k]] - tcrossprod(GUGt2[[k]], GUGt2[[k]]) * Maux[j, j]
          GVGt2[[k]] <- g + t3[[k]]
          GVG[[k]] <- GVG[[k]] - tcrossprod(GVGt2[[k]], GVGt2[[k]]) * invMU[j, j]
          invc1[[k]] <- .ridgeinv(GUG[[k]])
          invc2[[k]] <- .ridgeinv(GVG[[k]])
        }
        fobj <- function(x) {
          res <- -2 * log(1 + x %*% invt4 %*% x)
          for(k in 1:p){
            res <- res + auxf2(M[[k]], U[[k]], t2[[k]], t3[[k]], invc1[[k]], invc2[[k]], ng[k], n, p, x, j)
          }
          return(res)
        }
        gobj <- function(x) {
          res <- -4	* invt4 %*% x / as.numeric(1 + x %*% invt4 %*% x)
          for(k in 1:p){
            res <- res + auxg2(M[[k]], U[[k]], t2[[k]], t3[[k]], invc1[[k]], invc2[[k]], ng[k], n, p, x, j)
          }
          return(res)
        }
        if (rank == 1){
          res <- nlminb(Ginit[j,], fn = fobj, gr = gobj)
        } else{
          res <- nlminb(Ginit[j,], fn = fobj, gr = gobj)
        }
        Ginit[j, ] <- res$par
        g <- as.matrix(Ginit[j, ])
        t4 <- t4 + tcrossprod(g, g)
        for(k in 1:p){
          Maux <- M[[k]]
          MU <- M[[k]] + U[[k]]
          tmp.MU <- eigen(MU)
          invMU <- tcrossprod(sweep(tmp.MU$vectors, MARGIN = 2, 1 / tmp.MU$values, '*'), tmp.MU$vectors)
          GUGt2[[k]] <- g + t2[[k]]
          GUG[[k]] <- GUG[[k]] + tcrossprod(GUGt2[[k]], GUGt2[[k]]) * Maux[j, j]

          GVGt2[[k]] <- g + t3[[k]]
          GVG[[k]] <- GVG[[k]] + tcrossprod(GVGt2[[k]], GVGt2[[k]]) * invMU[j, j]
        }
      }
      a <- qr(Ginit)
      Gammahat <- qr.Q(a)
      obj5a <- c()
      for(i in 1 : p){
        e1 <- eigen(crossprod(Gammahat, (M[[i]] %*% Gammahat)))
        e2 <- eigen(crossprod(Gammahat, (invMU %*% Gammahat)))
        obj5a[i] <- sum(log(e1$values))*ng[i]/n + sum(log(e2$values))/p
      }
      obj5 <- sum(obj5a)
      if (abs(obj1 - obj5) < ftol * abs(obj1)) {
        break
      } else {
        obj1 <- obj5
        i <- i + 1
      }
    }
    Gamma0hat <- qr.Q(a, complete = TRUE)[, (rank + 1) : r]
    Gammahat <- as.matrix(Gammahat)
    Gamma0hat <- as.matrix(Gamma0hat)
  }
  return(list(Gammahat = Gammahat, Gamma0hat = Gamma0hat))
}


.GaussElim <- function(A) {
  a <- dim(A);n <- a[1];p <- a[2];idx <- rep(0, p);res.idx <- 1:n;i <- 1
  while (i <= p) {
    tmp <- max(abs(A[res.idx, i]))
    Stmp <- setdiff(which(abs(A[, i]) == tmp), idx)
    idx[i] <- Stmp[1]
    res.idx <- setdiff(res.idx, idx[i])
    for (j in 1:(n-i)) {
      A[res.idx[j], ] <- A[res.idx[j], ] - A[res.idx[j], i] / A[idx[i], i] * A[idx[i], ]
    }
    i <- i + 1
  }
  c(idx, res.idx)
}

.ridgeinv <- function(m){
  cvreg:::ridgepow(m, -1)
}

.expan <- function(d) {
  E <- matrix(0, d^2,d*(d+1)/2)
  for (i in 1:d) {
    for (j in 1:d) {
      if (j >= i) {
        E[d*(i-1)+j,(2*d-i)/2*(i-1)+j] <- 1
      } else {
        E[d*(i-1)+ j,(2*d-j)/2*(j-1)+i] <- 1
      }
    }
  }
  return(E)
}
.contr <- function(d) {
  C <- matrix(0, d * (d + 1) / 2, d ^ 2)
  for (i in 1:d) {
    for (j in 1:d) {
      if (j == i) {
        C[(2*d-i)/2*(i-1)+j,d*(i-1)+j] <- 1
      } else if (j > i) {
        C[(2*d-i)/2*(i-1)+j,d*(i-1)+j] <- 0.50
      } else {
        C[(2*d-j)/2*(j-1)+i,d*(i-1)+j] <- 0.50
      }
    }
  }
  return(C)
}