R/FitRDeltaQSym.R

In Correlplot: A Collection of Functions for Graphing Correlation Matrices

FitRDeltaQSym <- function(R,W=NULL,nd=2,eps=1e-6,delta.init=0,
                   q.init=rep(0,ncol(R)),itmax=1000,verbose=FALSE) {
  losshistory <- NULL
  rmsehistory <- NULL
  
  p <- ncol(R)
  J <- matrix(1,p,p)
  I <- diag(p)
  bone <- rep(1,p)
  if(is.null(W)) {
    W <- J - I
  } 
  W <- W/sum(W)
  
  sd <- eigen(R)
  
  if(nd > 1) {
    V  <- sd$vectors[,1:nd]
    Dl <- diag(sd$values[1:nd])  
  } else {
    V  <- as.matrix(sd$vectors[,1],ncol=1)
    Dl <- matrix(sd$values[1],1,1)
  }
  
  Rh <- V%*%Dl%*%t(V)
  E     <- R - Rh
  
  oldloss  <- tr(E%*%(W*E))
  i <- 1

  repeat {
    Radj <- R - delta.init*J - 0.5*bone%o%q.init - 0.5*q.init%o%bone
  
    out.wals <- wAddPCA(Radj, W, add = "nul", p = nd, verboseout = FALSE, 
                      epsout = 1e-8, bnd = "opt", epsin = 1e-8,
                      itmaxout = 5000, itmaxin = 20000)
  
    out.sd <- eigen(out.wals$y)
  
    if(any(round(out.sd$values,8) < 0)) {
      cat("warning: adjusted correlation matrix has negative eigenvalues.\n")
      print(out.sd$values)
    }
    
    V  <- out.sd$vectors[,1:nd]
    
    if(nd > 1) {
      V  <- out.sd$vectors[,1:nd]
      Dl <- diag(out.sd$values[1:nd])  
      Ds <- sqrt(Dl)
      G  <- V[,1:nd]%*%Ds[1:nd,1:nd]
    } else {
      V  <- out.sd$vectors[,1]
      Dl <- out.sd$values[1]
      Ds <- sqrt(Dl)
      G  <- V*Ds
      G  <- matrix(G,ncol=1)
    }
    
    GGt <- G%*%t(G)

    delta <- (tr(W%*%(R-GGt)) - p*sum(W*(R-GGt)))/(1-p)
    q     <- rowSums(p*W*(R-GGt)) - delta*bone
    E     <- R - delta*J - bone%o%q - GGt
    loss  <- tr(E%*%(W*E))
    delta.init <- delta
    q.init     <- q
 
    losshistory <- c(losshistory,loss)
    Rhat <- delta*J + bone%o%q + GGt
    rmsehistory <- c(rmsehistory,rmse(R,Rhat))
  
    decreaseloss <- oldloss - loss
    
    if(verbose) {
      cat("Iteration ", formatC(i, digits = 0, width = 5, format = "d"), 
          " loss", formatC(loss, digits = 6, width = 10, format = "f"),"\n")
    }
    
    if ((decreaseloss < eps) || (i == itmax)) {
      if (i == itmax) {
        cat("Maximum number of iterations (", itmax, ") reached.\n")
        cat("The algorithm may not have converged.\n")
      }
      break
    }
  
    oldloss <- loss
    i <- i + 1
  }
  
  Rhat <- delta*J + bone%o%q + GGt

  fit.rmse <- rmse(R,Rhat)
  
  imin <- which.min(losshistory)
  
  if(imin < length(losshistory)) {
    cat("Smallest loss does not occur at the last iteration but at iteration ",imin,".\n")
  }
  
  losshistory <- cbind(1:length(losshistory),losshistory)
  
  if(verbose) {
    cat("RMSE of the fit:\n")
    cat(formatC(fit.rmse, digits = 4, width = 10, format = "f"),"\n")
    cat("Final column adjustments:\n")
    final <- delta + q
    cat(formatC(final, digits = 4, width = 10, format = "f"),"\n")
  }
 
  return(list(delta=delta,q=q,G=G,fit.rmse=fit.rmse,
              losshistory=losshistory,
              rmsehistory=rmsehistory,
              Rhat=Rhat,
              eps=eps,nd=nd)) 
}