tests/online_stochastic_pca.R

In dfleis/morpca: Manifold Optimization for Robust PCA

#library(expm)
library(matrixcalc) # for is.positive.semi.definite() to check PSD while testing
library(Rcpp)
library(microbenchmark)
sourceCpp("test/logmat.cpp")
sourceCpp("test/expmat.cpp")

#===== functions =====#
project_re <- function(M, d, k) {
  eig <- eigen(M)
  V <- eig$vectors
  s <- eig$values
  
  S_new <- sapply(0:(d - 1), function(k_prime) {
    s[1:k_prime] <- 1/(d - k_prime) 
    s[(k_prime + 1):d] <- 
      (1 - sum(s[1:k_prime])) * s[(k_prime + 1):d]/sum(s[(k_prime + 1):d])
    s
  })
  
  S_all_idx <- which(apply(S_new, 2, function(x) sum(x <= 1/(d - k))) == d)
  S_idx <- max(S_all_idx)
  s_new <- S_new[,S_idx]
  
  V %*% tcrossprod(diag(s_new), V)
}

#===== parameters =====#
set.seed(124)
n <- 100
d <- 5
k <- 3

maxiter <- n

#===== generate data =====#
X <- matrix(rnorm(n * d), nrow = n)
I <- diag(d)

#===== initialize =====#
# generate PSD matrix M0 such that the largest eigenvalue 
# is <= 1/(d - k) and tr(M0) = 1
INIT_TMP <- matrix(rnorm(d^2), nrow = d, ncol = d)
SVD <- svd(INIT_TMP)
U0 <- SVD$u

#maxsig <- 1/(d - k)
#nb_maxsig <- 1/maxsig
#SIG0 <- diag(c(rep(maxsig, nb_maxsig), rep(0, d - nb_maxsig)))

SIG0 <- diag(abs(rnorm(d)))
M0 <- round(U0 %*% SIG0 %*% t(U0), 6)

M_list <- vector(mode = 'list', length = maxiter + 1)
M_list[[1]] <- M0

#===== start calculations =====#
pt <- proc.time()
for (i in 1:maxiter) {
  eta <- sqrt(1/i)
  x <- X[i,]
  
  xxt <- tcrossprod(x)
  #M_tmp <- expmat_cpp(logmat_cpp(M_list[[i]]) - eta * xxt)
  M_tmp <- expmat_cpp(expm::logm(M_list[[i]], method = "Eigen") - eta * xxt)
  
  M_list[[i + 1]] <- project_re(M_tmp, d, k)
}
proc.time() - pt

# calculate objective in (5)
obj <- sapply(M_list, function(M) {
  obj_all <- apply(X, 1, function(x) {
    sum(diag(M %*% tcrossprod(x)))
  })
  mean(obj_all)
})


M <- M_list[[maxiter + 1]]
Q <- qr.Q(qr(M))
R <- qr.R(qr(M))
chol(I - M)
rankMatrix(M)


#===== figures =====#
plot(obj[2:length(obj)], type = 'o', pch = 21, cex = 0.75, bg = 'white')

sapply(M_list, rankMatrix)



eig <- eigen(M0)
V <- eig$vectors
s <- rep(1, length(eig$values))
s[1] <- .Machine$double.eps

X <- V %*% diag(s) %*% t(V)
V %*% diag(exp(s)) %*% t(V)
expm(X)

round(V %*% diag(log(s)) %*% t(V) - logm(X), 6)
V %*% diag(log(s)) %*% t(V)
logm(X)