R/misc.R

In koekvall/mpredcc: Fully Likelihood-based Principal Components Regression

grad_wrap <- function(psi, opts){
  p <- ncol(opts$X)
  eg <- eigen(matrix(psi, p, p), symmetric = TRUE)
  return(grad_Psi(U = eg$vectors, d = eg$values, alpha = opts$alpha,
                  ssx = opts$ssx, Y = opts$Y, X = opts$X,
                  X_tilde = opts$X_tilde, Y_tilde = opts$Y_tilde))
}

prox_proj <- function(psi, t, opts)
{
  p <- ncol(opts$X)
  k <- opts$k
  Psi <- matrix(psi, p, p)
  ed <- eigen(Psi, symmetric = TRUE)
  ed$values[ed$values < 0] <- 0
  if(k < p){
    ed$values[(k + 1):p] <- 0
  }
  true_dim <- sum(ed$values > 0)
  if(true_dim == 0){
    return(matrix(0, p, p))
  } else{
    return(as.vector(ed$vectors[, 1:true_dim, drop = FALSE] %*% 
                       sweep(t(ed$vectors[, 1:true_dim, drop = FALSE]), 1,
                             ed$values[1:true_dim], FUN = "*")))  
  }
}

proj_evals <- function(d, k)
{
  d[d < 0] <- 0
  if(k < length(d)) d[(k + 1):length(d)] <- 0
  return(d)
}

project <- function(Psi, k){
  ed <- eigen(Psi, symmetric = TRUE)
  ed$values <- proj_evals(ed$values, k = k)
  return(ed$vectors %*% sweep(t(ed$vectors), 1, ed$values, FUN = "*"))
}

vec_norm <- function(x, q){
  return(sum(x^q)^(1 / q))
}

get_X_tilde<- function(X, H, phi)
{
  n <- nrow(X)
  X <- t(matrix(as.vector(t(X)) - H %*% phi, ncol = n))
  return(X)
}