From d95500c56e82dddb9ac5eb2e85abdefcd525f4fd Mon Sep 17 00:00:00 2001 From: daniel Date: Thu, 5 May 2022 13:35:07 +0200 Subject: [PATCH] implemented kpir_approx, wip: kpir_sim --- simulations/kpir_sim.R | 701 ++++++++++++++++++++++++++----- tensorPredictors/R/kpir_approx.R | 128 +++++- tensorPredictors/R/kpir_kron.R | 13 +- tensorPredictors/src/poi.c | 12 - 4 files changed, 713 insertions(+), 141 deletions(-) diff --git a/simulations/kpir_sim.R b/simulations/kpir_sim.R index a05cb7e..c10b83e 100644 --- a/simulations/kpir_sim.R +++ b/simulations/kpir_sim.R @@ -20,9 +20,8 @@ sim <- function(X, Fy, shape, alpha.true, beta.true, max.iter = 500L) { ) cat(sprintf( - "%3d | l = %-12.4f - dist = %-.4e - alpha(%d, %d) = %-.4e - beta(%d, %d) = %-.4e\n", - iter, loss, - dist, + "%s(%3d) | l = %-12.4f - dist = %-.4e - alpha(%d, %d) = %-.4e - beta(%d, %d) = %-.4e\n", + name, iter, loss, dist, nrow(alpha), ncol(alpha), dist.alpha, nrow(beta), ncol(beta), dist.beta )) @@ -30,12 +29,11 @@ sim <- function(X, Fy, shape, alpha.true, beta.true, max.iter = 500L) { } # Initialize logger history targets - hist.base <- hist.new <- hist.momentum <- # hist.kron <- + hist.base <- hist.new <- hist.momentum <- hist.approx <- # hist.kron <- data.frame(iter = seq(0L, max.iter), loss = NA, dist = NA, dist.alpha = NA, dist.beta = NA, norm.alpha = NA, norm.beta = NA ) - hist.kron <- NULL # TODO: fit kron version # Base (old) kpir.base(X, Fy, shape, max.iter = max.iter, logger = logger("base")) @@ -49,14 +47,56 @@ sim <- function(X, Fy, shape, alpha.true, beta.true, max.iter = 500L) { # # Residual Covariance Kronecker product assumpton version # kpir.kron(X, Fy, shape, max.iter = max.iter, logger = logger("kron")) + # Approximated MLE with Nesterov Momentum + kpir.approx(X, Fy, shape, max.iter = max.iter, logger = logger("approx")) + # Add method tags - hist.base$type <- factor("base") - hist.new$type <- factor("new") - hist.momentum$type <- factor("momentum") - # hist.kron$type <- factor("kron") + hist.base$method <- factor("base") + hist.new$method <- factor("new") + hist.momentum$method <- factor("momentum") + # hist.kron$method <- factor("kron") + hist.approx$method <- factor("approx") # Combine results and return - rbind(hist.base, hist.new, hist.momentum, hist.kron) + rbind(hist.base, hist.new, hist.momentum, hist.approx) #, hist.kron +} + +## Plot helper functions +plot.hist <- function(hist, response, ...) { + ggplot(hist, aes(x = iter, color = method, group = interaction(method, repetition))) + + geom_line(aes_(y = as.name(response)), na.rm = TRUE) + + geom_point(data = with(sub <- subset(hist, !is.na(as.symbol(response))), + aggregate(sub, list(method, repetition), tail, 1) + ), aes_(y = as.name(response))) + + labs(...) + + theme(legend.position = "bottom") +} +plot.stats <- function(hist, response, ..., title = "Stats") { + ggplot(hist, aes_(x = quote(iter), y = as.name(response), + color = quote(method), group = quote(method))) + + geom_ribbon(aes(color = NULL, fill = method), alpha = 0.2, + stat = "summary", fun.min = "min", fun.max = "max", na.rm = TRUE) + + geom_ribbon(aes(color = NULL, fill = method), alpha = 0.4, + stat = "summary", na.rm = TRUE, + fun.min = function(y) quantile(y, 0.25), + fun.max = function(y) quantile(y, 0.75)) + + geom_line(stat = "summary", fun = "mean", na.rm = TRUE) + + labs(title = title, ...) + + theme(legend.position = "bottom") +} +plot.mean <- function(hist, response, ..., title = "Mean") { + ggplot(hist, aes_(x = quote(iter), y = as.name(response), + color = quote(method), group = quote(method))) + + geom_line(stat = "summary", fun = "mean", na.rm = TRUE) + + labs(title = title, ...) + + theme(legend.position = "bottom") +} +plot.median <- function(hist, response, ..., title = "Median") { + ggplot(hist, aes_(x = quote(iter), y = as.name(response), + color = quote(method), group = quote(method))) + + geom_line(stat = "summary", fun = "median", na.rm = TRUE) + + labs(title = title, ...) + + theme(legend.position = "bottom") } ## Generate some test data / DEBUG @@ -68,7 +108,10 @@ r <- sample(1:15, 1) # 5 print(c(n, p, q, k, r)) hist <- NULL -for (rep in 1:20) { +reps <- 20 +for (rep in 1:reps) { + cat(sprintf("%4d / %d simulation rep. started\n", rep, reps)) + alpha.true <- alpha <- matrix(rnorm(q * r), q, r) beta.true <- beta <- matrix(rnorm(p * k), p, k) y <- rnorm(n) @@ -87,122 +130,481 @@ for (rep in 1:20) { hist <- rbind(hist, hist.sim) } +# Save simulation results +datetime <- format(Sys.time(), "%Y%m%dT%H%M") +saveRDS(hist, file = sprintf("AR_%s.rds", datetime)) -saveRDS(hist, file = "AR.rds") +# for GGPlot2, as factors for grouping hist$repetition <- factor(hist$repetition) -ggplot(hist, aes(x = iter, color = type, group = interaction(type, repetition))) + - geom_line(aes(y = loss)) + - geom_point(data = with(sub <- subset(hist, !is.na(loss)), - aggregate(sub, list(type, repetition), tail, 1) - ), aes(y = loss)) + - labs( - title = bquote(paste("Optimization Objective: negative log-likelihood ", - l(hat(alpha), hat(beta)))), - subtitle = bquote(paste(Delta[i][j] == 0.25, " * ", 0.5^abs(i - j), ", ", - "20 repetitions, ", n == .(n), ", ", - p == .(p), ", ", q == .(q), ", ", k == .(k), ", ", r == .(r))), - x = "nr. of iterations", - y = expression(l(hat(alpha), hat(beta))), - color = "method" - ) + - theme(legend.position = "bottom") +plot.hist(hist, "loss") +dev.print(png, file = sprintf("sim01_loss_%s.png", datetime), width = 768, height = 768, res = 125) +plot.stats(hist, "loss") +dev.print(png, file = sprintf("sim01_loss_stats_%s.png", datetime), width = 768, height = 768, res = 125) -dev.print(png, file = "sim01_loss.png", width = 768, height = 768, res = 125) +plot.hist(hist, "dist") +dev.print(png, file = sprintf("sim01_dist_%s.png", datetime), width = 768, height = 768, res = 125) +plot.stats(hist, "dist") +dev.print(png, file = sprintf("sim01_dist_stats_%s.png", datetime), width = 768, height = 768, res = 125) + +plot.hist(hist, "dist.alpha") +dev.print(png, file = sprintf("sim01_dist_alpha_%s.png", datetime), width = 768, height = 768, res = 125) +plot.stats(hist, "dist.alpha") +dev.print(png, file = sprintf("sim01_dist_alpha_stats_%s.png", datetime), width = 768, height = 768, res = 125) + +plot.hist(hist, "dist.beta") +dev.print(png, file = sprintf("sim01_dist_beta_%s.png", datetime), width = 768, height = 768, res = 125) +plot.stats(hist, "dist.beta") +dev.print(png, file = sprintf("sim01_dist_beta_stats_%s.png", datetime), width = 768, height = 768, res = 125) + +plot.hist(hist, "norm.alpha") +dev.print(png, file = sprintf("sim01_norm_alpha_%s.png", datetime), width = 768, height = 768, res = 125) +plot.stats(hist, "norm.alpha") +dev.print(png, file = sprintf("sim01_norm_alpha_stats_%s.png", datetime), width = 768, height = 768, res = 125) + +plot.hist(hist, "norm.beta") +dev.print(png, file = sprintf("sim01_norm_beta_%s.png", datetime), width = 768, height = 768, res = 125) +plot.stats(hist, "norm.beta") +dev.print(png, file = sprintf("sim01_norm_beta_stats_%s.png", datetime), width = 768, height = 768, res = 125) -ggplot(hist, aes(x = iter, color = type, group = interaction(type, repetition))) + - geom_line(aes(y = dist)) + - geom_point(data = with(sub <- subset(hist, !is.na(dist)), - aggregate(sub, list(type, repetition), tail, 1) - ), aes(y = dist)) + - labs( - title = bquote(paste("Distance of estimate ", hat(B), " to true ", B == alpha %*% beta)), - subtitle = bquote(paste(Delta[i][j] == 0.25, " * ", 0.5^abs(i - j), ", ", - "20 repetitions, ", n == .(n), ", ", - p == .(p), ", ", q == .(q), ", ", k == .(k), ", ", r == .(r))), - x = "nr. of iterations", - y = expression(abs(B * B^T - hat(B) * hat(B)^T)), - color = "method" - ) + - theme(legend.position = "bottom") - -dev.print(png, file = "sim01_dist.png", width = 768, height = 768, res = 125) -ggplot(hist, aes(x = iter, color = type, group = interaction(type, repetition))) + - geom_line(aes(y = dist.alpha)) + - geom_point(data = with(sub <- subset(hist, !is.na(dist.alpha)), - aggregate(sub, list(type, repetition), tail, 1) - ), aes(y = dist.alpha)) + - labs( - title = bquote(paste("Distance of estimate ", hat(alpha), " to true ", alpha)), - subtitle = bquote(paste(Delta[i][j] == 0.25, " * ", 0.5^abs(i - j), ", ", - "20 repetitions, ", n == .(n), ", ", - p == .(p), ", ", q == .(q), ", ", k == .(k), ", ", r == .(r))), - x = "nr. of iterations", - y = expression(abs(alpha * alpha^T - hat(alpha) * hat(alpha)^T)), - color = "method" - ) + - theme(legend.position = "bottom") -dev.print(png, file = "sim01_dist_alpha.png", width = 768, height = 768, res = 125) +n <- 200 # Sample Size +p <- 11 # sample(1:15, 1) +q <- 3 # sample(1:15, 1) +k <- 7 # sample(1:15, 1) +r <- 5 # sample(1:15, 1) +print(c(n, p, q, k, r)) + +hist <- NULL +reps <- 20 + +Delta.1 <- sqrt(0.5)^abs(outer(seq_len(q), seq_len(q), `-`)) +Delta.2 <- sqrt(0.5)^abs(outer(seq_len(p), seq_len(p), `-`)) +Delta <- kronecker(Delta.1, Delta.2) +for (rep in 1:reps) { + cat(sprintf("%4d / %d simulation rep. started\n", rep, reps)) + + alpha.true <- alpha <- matrix(rnorm(q * r), q, r) + beta.true <- beta <- matrix(rnorm(p * k), p, k) + y <- rnorm(n) + Fy <- do.call(cbind, Map(function(slope, offset) { + sin(slope * y + offset) + }, + head(rep(seq(1, ceiling(0.5 * k * r)), each = 2), k * r), + head(rep(c(0, pi / 2), ceiling(0.5 * k * r)), k * r) + )) + X <- tcrossprod(Fy, kronecker(alpha, beta)) + CVarE:::rmvnorm(n, sigma = Delta) + + hist.sim <- sim(X, Fy, shape = c(p, q, k, r), alpha.true, beta.true) + hist.sim$repetition <- rep + + hist <- rbind(hist, hist.sim) +} + + +# Save simulation results +datetime <- format(Sys.time(), "%Y%m%dT%H%M") +saveRDS(hist, file = sprintf("sim02_%s.rds", datetime)) + +# for GGPlot2, as factors for grouping +hist$repetition <- factor(hist$repetition) + +plot.hist(hist, "loss") +dev.print(png, file = sprintf("sim02_loss_%s.png", datetime), width = 768, height = 768, res = 125) +plot.stats(hist, "loss") +dev.print(png, file = sprintf("sim02_loss_stats_%s.png", datetime), width = 768, height = 768, res = 125) + +plot.hist(hist, "dist") +dev.print(png, file = sprintf("sim02_dist_%s.png", datetime), width = 768, height = 768, res = 125) +plot.stats(hist, "dist") +dev.print(png, file = sprintf("sim02_dist_stats_%s.png", datetime), width = 768, height = 768, res = 125) +plot.mean(hist, "dist") +plot.median(hist, "dist") + +plot.hist(hist, "dist.alpha") +dev.print(png, file = sprintf("sim02_dist_alpha_%s.png", datetime), width = 768, height = 768, res = 125) +plot.stats(hist, "dist.alpha") +dev.print(png, file = sprintf("sim02_dist_alpha_stats_%s.png", datetime), width = 768, height = 768, res = 125) +plot.mean(hist, "dist.alpha") +plot.median(hist, "dist.alpha") + +plot.hist(hist, "dist.beta") +dev.print(png, file = sprintf("sim02_dist_beta_%s.png", datetime), width = 768, height = 768, res = 125) +plot.stats(hist, "dist.beta") +dev.print(png, file = sprintf("sim02_dist_beta_stats_%s.png", datetime), width = 768, height = 768, res = 125) +plot.mean(hist, "dist.beta") +plot.median(hist, "dist.beta") + +plot.hist(hist, "norm.alpha") +dev.print(png, file = sprintf("sim02_norm_alpha_%s.png", datetime), width = 768, height = 768, res = 125) +plot.stats(hist, "norm.alpha") +dev.print(png, file = sprintf("sim02_norm_alpha_stats_%s.png", datetime), width = 768, height = 768, res = 125) + +plot.hist(hist, "norm.beta") +dev.print(png, file = sprintf("sim02_norm_beta_%s.png", datetime), width = 768, height = 768, res = 125) +plot.stats(hist, "norm.beta") +dev.print(png, file = sprintf("sim02_norm_beta_stats_%s.png", datetime), width = 768, height = 768, res = 125) + +plot.hist2 <- function(hist, response, type = "all", ...) { + # Extract final results from history + sub <- na.omit(hist[c("iter", response, "method", "repetition")]) + sub <- aggregate(sub, list(sub$method, sub$repetition), tail, 1) + + # Setup ggplot + p <- ggplot(hist, aes_(x = quote(iter), + y = as.name(response), + color = quote(method), + group = quote(interaction(method, repetition)))) + # Add requested layers + if (type == "all") { + p <- p + geom_line(na.rm = TRUE) + p <- p + geom_point(data = sub) + } else if (type == "mean") { + p <- p + geom_line(alpha = 0.5, na.rm = TRUE, linetype = "dotted") + p <- p + geom_point(data = sub, alpha = 0.5) + p <- p + geom_line(aes(group = method), + stat = "summary", fun = "mean", na.rm = TRUE) + } else if (type == "median") { + p <- p + geom_line(alpha = 0.5, na.rm = TRUE, linetype = "dotted") + p <- p + geom_point(data = sub, alpha = 0.5) + p <- p + geom_line(aes(group = method), + stat = "summary", fun = "median", na.rm = TRUE) + } + # return with theme and annotations + p + labs(...) + theme(legend.position = "bottom") +} + +plot.hist2(hist, "dist.alpha", "all", title = "all") + coord_trans(x = "log1p") +plot.hist2(hist, "dist.alpha", "mean", title = "mean") + coord_trans(x = "log1p") +plot.hist2(hist, "dist.alpha", "median", title = "median") + coord_trans(x = "log1p") + + +################################################################################ +### EEG ### +################################################################################ + +suppressPackageStartupMessages({ + library(pROC) +}) + +# acc: Accuracy. P(Yhat = Y). Estimated as: (TP+TN)/(P+N). +acc <- function(y_true, y_pred) mean(round(y_pred) == y_true) +# err: Error rate. P(Yhat != Y). Estimated as: (FP+FN)/(P+N). +err <- function(y_true, y_pred) mean(round(y_pred) != y_true) +# fpr: False positive rate. P(Yhat = + | Y = -). aliases: Fallout. +fpr <- function(y_true, y_pred) mean((round(y_pred) == 1)[y_true == 0]) +# tpr: True positive rate. P(Yhat = + | Y = +). aliases: Sensitivity, Recall. +tpr <- function(y_true, y_pred) mean((round(y_pred) == 1)[y_true == 1]) +# fnr: False negative rate. P(Yhat = - | Y = +). aliases: Miss. +fnr <- function(y_true, y_pred) mean((round(y_pred) == 0)[y_true == 1]) +# tnr: True negative rate. P(Yhat = - | Y = -). +tnr <- function(y_true, y_pred) mean((round(y_pred) == 0)[y_true == 0]) + +# Load EEG dataset +dataset <- readRDS('eeg_analysis/eeg_data.rds') + +eeg_cross_validation <- function(nrFolds = 10L) { + # Set dimenional parameters. + n <- nrow(dataset) # sample size (nr. of people) + p <- 64L # nr. of predictors (count of sensorce) + t <- 256L # nr. of time points (measurements) + + # Extract dimension names from X. + nNames <- dataset$PersonID + tNames <- as.character(seq(t)) + pNames <- unlist(strsplit(colnames(dataset)[2 + t * seq(p)], '_'))[c(T, F)] + + # Split into X-y. + X <- as.matrix(dataset[, -(1:2)]) + y <- dataset$Case_Control + # Reshape X as 3D tenros of shape (n, t, p) aka. samples, timesteps, predictors. + # (Each of the n rows in X iterate over the time bevore switching sensorce.) + dim(X) <- c(n, t, p) + dimnames(X) <- list(nNames, tNames, pNames) + + # Setup Cross-Validation result + CV <- data.frame( + fold = (seq_len(n) %% nrFolds) + 1L, + y_true = y, + y_pred = NA + ) + + # + +} + +#' @param ppc Number of "p"redictor "p"rincipal "c"omponents. +#' @param tpc Number of "t"ime "p"rincipal "c"omponents. +egg_analysis_reduced <- function(methods, ppc, tpc) { + # Set dimenional parameters. + n <- nrow(dataset) # sample size (nr. of people) + p <- 64L # nr. of predictors (count of sensorce) + t <- 256L # nr. of time points (measurements) + + # Extract dimension names from X. + nNames <- dataset$PersonID + tNames <- as.character(seq(t)) + pNames <- unlist(strsplit(colnames(dataset)[2 + t * seq(p)], '_'))[c(T, F)] + + # Split into X-y. + X <- as.matrix(dataset[, -(1:2)]) + y <- dataset$Case_Control + # Reshape X as 3D tenros of shape (n, t, p) aka. samples, timesteps, predictors. + # (Each of the n rows in X iterate over the time bevore switching sensorce.) + X <- array(X, dim = c(n, t, p), + dimnames = list(nNames, tNames, pNames)) + # Reorder axis to (p, t, n) = (predictors, timesteps, samples). + X <- aperm(X, c(3, 2, 1)) + + # Compute Mean of X. + X_mean <- apply(X, c(1, 2), mean) + X_center <- X - as.vector(X_mean) + + # Compute "left" and "right" cov-matrices. + Sigma_t <- matrix(apply(apply(X_center, 3, crossprod), 1, mean), t, t) + Sigma_p <- matrix(apply(apply(X_center, 3, tcrossprod), 1, mean), p, p) + # Get "left", "right" principal components. + V_p <- svd(Sigma_p, ppc, 0L)$u + V_t <- svd(Sigma_t, tpc, 0L)$u + + # Reduce dimension. + X_reduced <- apply(X_center, 3, function(x) crossprod(V_p, x %*% V_t)) + dim(X_reduced) <- c(ppc, tpc, n) + + # Vectorize to shape of (predictors * timesteps, samples) and transpose to + # (samples, predictors * timesteps). + X_vec <- t(matrix(X_reduced, ppc * tpc, n)) + + loo.cv <- expand.grid(method = names(methods), fold = 1:n) + loo.cv$y_true <- y[loo.cv$fold] + loo.cv$y_pred <- NA + + # Performe LOO cross-validation for each method. + for (i in 1L:n) { + # Print progress. + cat(sprintf("\rCross-Validation (p-PC: %d, t-PC: %d): %4d/%d", + ppc, tpc, i, n)) + # Leave Out the i-th element. + X_train <- X_vec[-i, ] + X_test <- X_vec[i, ] + y_train <- y[-i] + # Center y. + y_train <- scale(y_train, center = TRUE, scale = FALSE) + + # For each method. + for (method.name in names(methods)) { + method <- methods[[method.name]] + # Compute reduction using current method under common API. + sdr <- method(X_train, y_train, ppc, tpc) + B <- kronecker(sdr$alpha, sdr$beta) + # Fit a linear model (which ensures a common sdr direction if possible). + model <- glm(y ~ x, family = binomial(link = "logit"), + data = data.frame(y = y[-i], x = X_train %*% B)) + # Predict out of sample and store in LOO CV data.frame. + y_pred <- predict(model, data.frame(x = X_test %*% B), type = "response") + loo.cv[loo.cv$method == method.name & loo.cv$fold == i, 'y_pred'] <- y_pred + } + } + + for (method.name in names(methods)) { + labels <- loo.cv[loo.cv$method == method.name, 'y_true'] + predictions <- loo.cv[loo.cv$method == method.name, 'y_pred'] + ROC <- roc(unlist(labels), unlist(predictions), quiet = TRUE) + # Combined accuracy, error, ... + cat("\nMethod: ", method.name, "\n", + "acc: ", acc(unlist(labels), unlist(predictions)), "\n", + "err: ", err(unlist(labels), unlist(predictions)), "\n", + "fpr: ", fpr(unlist(labels), unlist(predictions)), "\n", + "tpr: ", tpr(unlist(labels), unlist(predictions)), "\n", + "fnr: ", fnr(unlist(labels), unlist(predictions)), "\n", + "tnr: ", tnr(unlist(labels), unlist(predictions)), "\n", + "auc: ", ROC$auc, "\n", + "auc sd: ", sqrt(var(ROC)), "\n", + sep = '') + } + + loo.cv +} + +methods <- list( + KPIR_LS = function(...) tensor_predictor(..., method = "KPIR_LS"), + KPIR_MLE = function(...) tensor_predictor(..., method = "KPIR_MLE"), + KPFC1 = function(...) tensor_predictor(..., method = "KPFC1"), + KPFC2 = function(...) tensor_predictor(..., method = "KPFC2"), + LSIR = LSIR +) + +# ppc, tpc +# ------------ +params <- list( c( 4, 3) + , c( 15, 15) + , c( 30, 20) +) + +for (param in params) { + c(ppc, tpc) %<-% param + sim <- egg_analysis_reduced(methods, ppc, tpc) + + attr(sim, 'param') <- c(ppc = ppc, tpc = tpc) + + saveRDS(sim, file = sprintf('eeg_analysis_reduced_%d_%d.rds', ppc, tpc)) +} + + + +# plot.hist(hist, "loss", +# title = bquote(paste("Optimization Objective: negative log-likelihood ", +# l(hat(alpha), hat(beta)))), +# subtitle = bquote(paste(Delta[i][j] == 0.25, " * ", 0.5^abs(i - j), ", ", +# "20 repetitions, ", n == .(n), ", ", +# p == .(p), ", ", q == .(q), ", ", k == .(k), ", ", r == .(r))), +# x = "nr. of iterations", +# y = expression(l(hat(alpha), hat(beta))) +# ) +# plot.stats(hist, "loss", +# title = bquote(paste("Optimization Objective: negative log-likelihood ", +# l(hat(alpha), hat(beta)))), +# subtitle = bquote(paste(Delta[i][j] == 0.25, " * ", 0.5^abs(i - j), ", ", +# "20 repetitions, ", n == .(n), ", ", +# p == .(p), ", ", q == .(q), ", ", k == .(k), ", ", r == .(r))), +# x = "nr. of iterations", +# y = expression(l(hat(alpha), hat(beta))) +# ) + + +# dev.print(png, file = sprintf("sim01_loss_stat_%s.png", datetime), +# width = 768, height = 768, res = 125) + + +# ggplot(hist, aes(x = iter, color = method, group = interaction(method, repetition))) + +# geom_line(aes(y = dist)) + +# geom_point(data = with(sub <- subset(hist, !is.na(dist)), +# aggregate(sub, list(method, repetition), tail, 1) +# ), aes(y = dist)) + +# labs( +# title = bquote(paste("Distance of estimate ", hat(B), " to true ", B == alpha %*% beta)), +# subtitle = bquote(paste(Delta[i][j] == 0.25, " * ", 0.5^abs(i - j), ", ", +# "20 repetitions, ", n == .(n), ", ", +# p == .(p), ", ", q == .(q), ", ", k == .(k), ", ", r == .(r))), +# x = "nr. of iterations", +# y = expression(abs(B * B^T - hat(B) * hat(B)^T)), +# color = "method" +# ) + +# theme(legend.position = "bottom") + +# dev.print(png, file = sprintf("sim01_dist_%s.png", datetime), +# width = 768, height = 768, res = 125) + +# ggplot(hist, aes(x = iter, y = dist, color = method, group = method)) + +# geom_ribbon(aes(color = NULL, fill = method), alpha = 0.2, +# stat = "summary", fun.min = "min", fun.max = "max", na.rm = TRUE) + +# geom_ribbon(aes(color = NULL, fill = method), alpha = 0.4, +# stat = "summary", fun.min = function(y) quantile(y, 0.25), +# fun.max = function(y) quantile(y, 0.75), na.rm = TRUE) + +# geom_line(stat = "summary", fun = "mean", na.rm = TRUE) + +# labs( +# title = bquote(paste("Distance of estimate ", hat(B), " to true ", B == alpha %*% beta)), +# subtitle = bquote(paste(Delta[i][j] == 0.25, " * ", 0.5^abs(i - j), ", ", +# "20 repetitions, ", n == .(n), ", ", +# p == .(p), ", ", q == .(q), ", ", k == .(k), ", ", r == .(r))), +# x = "nr. of iterations", +# y = expression(abs(B * B^T - hat(B) * hat(B)^T)), +# color = "method" +# ) + +# theme(legend.position = "bottom") + +# dev.print(png, file = sprintf("sim01_dist_stat_%s.png", datetime), +# width = 768, height = 768, res = 125) + + +# ggplot(hist, aes(x = iter, color = method, group = interaction(method, repetition))) + +# geom_line(aes(y = dist.alpha)) + +# geom_point(data = with(sub <- subset(hist, !is.na(dist.alpha)), +# aggregate(sub, list(method, repetition), tail, 1) +# ), aes(y = dist.alpha)) + +# labs( +# title = bquote(paste("Distance of estimate ", hat(alpha), " to true ", alpha)), +# subtitle = bquote(paste(Delta[i][j] == 0.25, " * ", 0.5^abs(i - j), ", ", +# "20 repetitions, ", n == .(n), ", ", +# p == .(p), ", ", q == .(q), ", ", k == .(k), ", ", r == .(r))), +# x = "nr. of iterations", +# y = expression(abs(alpha * alpha^T - hat(alpha) * hat(alpha)^T)), +# color = "method" +# ) + +# theme(legend.position = "bottom") + +# dev.print(png, file = sprintf("sim01_dist_alpha_%s.png", datetime), +# width = 768, height = 768, res = 125) + + +# ggplot(hist, aes(x = iter, color = method, group = interaction(method, repetition))) + +# geom_line(aes(y = dist.beta)) + +# geom_point(data = with(sub <- subset(hist, !is.na(dist.beta)), +# aggregate(sub, list(method, repetition), tail, 1) +# ), aes(y = dist.beta)) + +# labs( +# title = bquote(paste("Distance of estimate ", hat(beta), " to true ", beta)), +# subtitle = bquote(paste(Delta[i][j] == 0.25, " * ", 0.5^abs(i - j), ", ", +# "20 repetitions, ", n == .(n), ", ", +# p == .(p), ", ", q == .(q), ", ", k == .(k), ", ", r == .(r))), +# x = "nr. of iterations", +# y = expression(abs(beta * beta^T - hat(beta) * hat(beta)^T)), +# color = "method" +# ) + +# theme(legend.position = "bottom") + +# dev.print(png, file = sprintf("sim01_dist_beta_%s.png", datetime), +# width = 768, height = 768, res = 125) + + +# ggplot(hist, aes(x = iter, color = method, group = interaction(method, repetition))) + +# geom_line(aes(y = norm.alpha)) + +# geom_point(data = with(sub <- subset(hist, !is.na(norm.alpha)), +# aggregate(sub, list(method, repetition), tail, 1) +# ), aes(y = norm.alpha)) + +# labs( +# title = expression(paste("Norm of ", hat(alpha))), +# subtitle = bquote(paste(Delta[i][j] == 0.25, " * ", 0.5^abs(i - j), ", ", +# "20 repetitions, ", n == .(n), ", ", +# p == .(p), ", ", q == .(q), ", ", k == .(k), ", ", r == .(r))), +# x = "nr. of iterations", +# y = expression(abs(hat(alpha))[F]), +# color = "method" +# ) + +# theme(legend.position = "bottom") + +# dev.print(png, file = sprintf("sim01_norm_alpha_%s.png", datetime), +# width = 768, height = 768, res = 125) + +# ggplot(hist, aes(x = iter, color = method, group = interaction(method, repetition))) + +# geom_line(aes(y = norm.beta)) + +# geom_point(data = with(sub <- subset(hist, !is.na(norm.beta)), +# aggregate(sub, list(method, repetition), tail, 1) +# ), aes(y = norm.beta)) + +# labs( +# title = expression(paste("Norm of ", hat(beta))), +# subtitle = bquote(paste(Delta[i][j] == 0.25, " * ", 0.5^abs(i - j), ", ", +# "20 repetitions, ", n == .(n), ", ", +# p == .(p), ", ", q == .(q), ", ", k == .(k), ", ", r == .(r))), +# x = "nr. of iterations", +# y = expression(abs(hat(beta))[F]), +# color = "method" +# ) + +# theme(legend.position = "bottom") + +# dev.print(png, file = sprintf("sim01_norm_beta_%s.png", datetime), +# width = 768, height = 768, res = 125) -ggplot(hist, aes(x = iter, color = type, group = interaction(type, repetition))) + - geom_line(aes(y = dist.beta)) + - geom_point(data = with(sub <- subset(hist, !is.na(dist.beta)), - aggregate(sub, list(type, repetition), tail, 1) - ), aes(y = dist.beta)) + - labs( - title = bquote(paste("Distance of estimate ", hat(beta), " to true ", beta)), - subtitle = bquote(paste(Delta[i][j] == 0.25, " * ", 0.5^abs(i - j), ", ", - "20 repetitions, ", n == .(n), ", ", - p == .(p), ", ", q == .(q), ", ", k == .(k), ", ", r == .(r))), - x = "nr. of iterations", - y = expression(abs(beta * beta^T - hat(beta) * hat(beta)^T)), - color = "method" - ) + - theme(legend.position = "bottom") -dev.print(png, file = "sim01_dist_beta.png", width = 768, height = 768, res = 125) -ggplot(hist, aes(x = iter, color = type, group = interaction(type, repetition))) + - geom_line(aes(y = norm.alpha)) + - geom_point(data = with(sub <- subset(hist, !is.na(norm.alpha)), - aggregate(sub, list(type, repetition), tail, 1) - ), aes(y = norm.alpha)) + - labs( - title = expression(paste("Norm of ", hat(alpha))), - subtitle = bquote(paste(Delta[i][j] == 0.25, " * ", 0.5^abs(i - j), ", ", - "20 repetitions, ", n == .(n), ", ", - p == .(p), ", ", q == .(q), ", ", k == .(k), ", ", r == .(r))), - x = "nr. of iterations", - y = expression(abs(hat(alpha))[F]), - color = "method" - ) + - theme(legend.position = "bottom") -dev.print(png, file = "sim01_norm_alpha.png", width = 768, height = 768, res = 125) -ggplot(hist, aes(x = iter, color = type, group = interaction(type, repetition))) + - geom_line(aes(y = norm.beta)) + - geom_point(data = with(sub <- subset(hist, !is.na(norm.beta)), - aggregate(sub, list(type, repetition), tail, 1) - ), aes(y = norm.beta)) + - labs( - title = expression(paste("Norm of ", hat(beta))), - subtitle = bquote(paste(Delta[i][j] == 0.25, " * ", 0.5^abs(i - j), ", ", - "20 repetitions, ", n == .(n), ", ", - p == .(p), ", ", q == .(q), ", ", k == .(k), ", ", r == .(r))), - x = "nr. of iterations", - y = expression(abs(hat(beta))[F]), - color = "method" - ) + - theme(legend.position = "bottom") -dev.print(png, file = "sim01_norm_beta.png", width = 768, height = 768, res = 125) @@ -983,3 +1385,82 @@ dev.print(png, file = "sim01_norm_beta.png", width = 768, height = 768, res = 12 # log(det(kronecker(D.1, D.2))) # p * log(det(D.1)) + q * log(det(D.2)) + + + + d <- ggplot(mtcars, aes(cyl, mpg)) + geom_point() + d + stat_summary(fun.data = "mean_cl_boot", colour = "red", size = 2) + + # Orientation follows the discrete axis + ggplot(mtcars, aes(mpg, factor(cyl))) + + geom_point() + + stat_summary(fun.data = "mean_cl_boot", colour = "red", size = 2) + + # You can supply individual functions to summarise the value at + # each x: + d + stat_summary(fun = "median", colour = "red", size = 2, geom = "point") + d + stat_summary(fun = "mean", colour = "red", size = 2, geom = "point") + d + aes(colour = factor(vs)) + stat_summary(fun = mean, geom="line") + + d + stat_summary(fun = mean, fun.min = min, fun.max = max, colour = "red") + + d <- ggplot(diamonds, aes(cut)) + d + geom_bar() + d + stat_summary(aes(y = price), fun = "mean", geom = "bar") + + # Orientation of stat_summary_bin is ambiguous and must be specified directly + ggplot(diamonds, aes(carat, price)) + + stat_summary_bin(fun = "mean", geom = "bar", orientation = 'y') + + + # Don't use ylim to zoom into a summary plot - this throws the + # data away + p <- ggplot(mtcars, aes(cyl, mpg)) + + stat_summary(fun = "mean", geom = "point") + p + p + ylim(15, 30) + # Instead use coord_cartesian + p + coord_cartesian(ylim = c(15, 30)) + + # A set of useful summary functions is provided from the Hmisc package: + stat_sum_df <- function(fun, geom="crossbar", ...) { + stat_summary(fun.data = fun, colour = "red", geom = geom, width = 0.2, ...) + } + d <- ggplot(mtcars, aes(cyl, mpg)) + geom_point() + # The crossbar geom needs grouping to be specified when used with + # a continuous x axis. + d + stat_sum_df("mean_cl_boot", mapping = aes(group = cyl)) + d + stat_sum_df("mean_sdl", mapping = aes(group = cyl)) + d + stat_sum_df("mean_sdl", fun.args = list(mult = 1), mapping = aes(group = cyl)) + d + stat_sum_df("median_hilow", mapping = aes(group = cyl)) + + # An example with highly skewed distributions: + if (require("ggplot2movies")) { + set.seed(596) + mov <- movies[sample(nrow(movies), 1000), ] + m2 <- + ggplot(mov, aes(x = factor(round(rating)), y = votes)) + + geom_point() + m2 <- + m2 + + stat_summary( + fun.data = "mean_cl_boot", + geom = "crossbar", + colour = "red", width = 0.3 + ) + + xlab("rating") + m2 + # Notice how the overplotting skews off visual perception of the mean + # supplementing the raw data with summary statistics is _very_ important + + # Next, we'll look at votes on a log scale. + + # Transforming the scale means the data are transformed + # first, after which statistics are computed: + m2 + scale_y_log10() + # Transforming the coordinate system occurs after the + # statistic has been computed. This means we're calculating the summary on the raw data + # and stretching the geoms onto the log scale. Compare the widths of the + # standard errors. + m2 + coord_trans(y="log10") + } \ No newline at end of file diff --git a/tensorPredictors/R/kpir_approx.R b/tensorPredictors/R/kpir_approx.R index 28cd564..084f5ad 100644 --- a/tensorPredictors/R/kpir_approx.R +++ b/tensorPredictors/R/kpir_approx.R @@ -4,7 +4,7 @@ #' #' Delta_1 = n^-1 sum_i R_i' R_i, #' Delta_2 = n^-1 sum_i R_i R_i'. -#' +#' #' @export kpir.approx <- function(X, Fy, shape = c(dim(X)[-1], dim(Fy[-1])), max.iter = 500L, max.line.iter = 50L, step.size = 1e-3, @@ -91,14 +91,8 @@ kpir.approx <- function(X, Fy, shape = c(dim(X)[-1], dim(Fy[-1])), S.2 <- n^-1 * tcrossprod(mat(R, 2), mat(R %x_3% Delta.1.inv, 2)) # Evaluate negative log-likelihood (2 pi term dropped) - loss <- -0.5 * n * (p * q * log(s) - p * log(det(Delta.1)) - - q * log(det(Delta.2)) - s * sum(S.1 * Delta.1.inv)) - - # Gradient "generating" tensor - G <- (sum(S.1 * Delta.1.inv) - p * q / s) * R - G <- G + R %x_2% ((diag(q, p, p) - s * (Delta.2.inv %*% S.2)) %*% Delta.2.inv) - G <- G + R %x_3% ((diag(p, q, q) - s * (Delta.1.inv %*% S.1)) %*% Delta.1.inv) - G <- G + s * (R %x_2% Delta.2.inv %x_3% Delta.1.inv) + loss <- -0.5 * (n * (p * q * log(s) - p * log(det(Delta.1)) - + q * log(det(Delta.2))) - s * sum(S.1 * Delta.1.inv)) # Call history callback (logger) before the first iteration if (is.function(logger)) { @@ -125,12 +119,116 @@ kpir.approx <- function(X, Fy, shape = c(dim(X)[-1], dim(Fy[-1])), alpha.moment <- alpha1 + ((a0 - 1) / a1) * (alpha1 - alpha0) beta.moment <- beta1 + ((a0 - 1) / a1) * ( beta1 - beta0) } + + # Extrapolated residuals + R <- X - (Fy %x_3% alpha.moment %x_2% beta.moment) + + # Recompute Covariance Estimates and scaling factor + Delta.1 <- tcrossprod(mat(R, 3)) + Delta.2 <- tcrossprod(mat(R, 2)) + s <- sum(diag(Delta.1)) + + # Inverse Covariances + Delta.1.inv <- solve(Delta.1) + Delta.2.inv <- solve(Delta.2) + + # cross dependent covariance estimates + S.1 <- n^-1 * tcrossprod(mat(R, 3), mat(R %x_2% Delta.2.inv, 3)) + S.2 <- n^-1 * tcrossprod(mat(R, 2), mat(R %x_3% Delta.1.inv, 2)) + + # Gradient "generating" tensor + G <- (sum(S.1 * Delta.1.inv) - p * q / s) * R + G <- G + R %x_2% ((diag(q, p) - s * (Delta.2.inv %*% S.2)) %*% Delta.2.inv) + G <- G + R %x_3% ((diag(p, q) - s * (Delta.1.inv %*% S.1)) %*% Delta.1.inv) + G <- G + s * (R %x_2% Delta.2.inv %x_3% Delta.1.inv) + + # Calculate Gradients + grad.alpha <- tcrossprod(mat(G, 3), mat(Fy %x_2% beta.moment, 3)) + grad.beta <- tcrossprod(mat(G, 2), mat(Fy %x_3% alpha.moment, 2)) + + # Backtracking line search (Armijo type) + # The `inner.prod` is used in the Armijo break condition but does not + # depend on the step size. + inner.prod <- sum(grad.alpha^2) + sum(grad.beta^2) + + # backtracking loop + for (delta in step.size * 0.618034^seq.int(0L, len = max.line.iter)) { + # Update `alpha` and `beta` (note: add(+), the gradients are already + # pointing into the negative slope direction of the loss cause they are + # the gradients of the log-likelihood [NOT the negative log-likelihood]) + alpha.temp <- alpha.moment + delta * grad.alpha + beta.temp <- beta.moment + delta * grad.beta + + # Update Residuals, Covariances, ... + R <- X - (Fy %x_3% alpha.temp %x_2% beta.temp) + Delta.1 <- tcrossprod(mat(R, 3)) + Delta.2 <- tcrossprod(mat(R, 2)) + s <- sum(diag(Delta.1)) + Delta.1.inv <- solve(Delta.1) + Delta.2.inv <- solve(Delta.2) + S.1 <- n^-1 * tcrossprod(mat(R, 3), mat(R %x_2% Delta.2.inv, 3)) + # S.2 not needed + + # Re-evaluate negative log-likelihood + loss.temp <- -0.5 * (n * (p * q * log(s) - p * log(det(Delta.1)) - + q * log(det(Delta.2))) - s * sum(S.1 * Delta.1.inv)) + + # Armijo line search break condition + if (loss.temp <= loss - 0.1 * delta * inner.prod) { + break + } + } + + # Call logger (invoke history callback) + if (is.function(logger)) { + logger(iter, loss.temp, alpha.temp, beta.temp, Delta.1, Delta.2, delta) + } + + # Enforce descent + if (loss.temp < loss) { + alpha0 <- alpha1 + alpha1 <- alpha.temp + beta0 <- beta1 + beta1 <- beta.temp + + # check break conditions + if (mean(abs(alpha1)) + mean(abs(beta1)) < eps) { + break.reason <- "alpha, beta numerically zero" + break # estimates are basically zero -> stop + } + if (inner.prod < eps * (p * q + r * k)) { + break.reason <- "mean squared gradient is smaller than epsilon" + break # mean squared gradient is smaller than epsilon -> stop + } + if (abs(loss.temp - loss) < eps) { + break.reason <- "decrease is too small (slow)" + break # decrease is too small (slow) -> stop + } + + loss <- loss.temp + no.nesterov <- FALSE # always reset + } else if (!no.nesterov) { + no.nesterov <- TRUE # retry without momentum + next + } else { + break.reason <- "failed even without momentum" + break # failed even without momentum -> stop + } + + # update momentum scaling + a0 <- a1 + a1 <- nesterov.scaling(a1, iter) + + # Set next iter starting step.size to line searched step size + # (while allowing it to encrease) + step.size <- 1.618034 * delta + } - # Extrapolated residuals - R <- X - (Fy %x_3% alpha.moment %x_2% beta.moment) - - - - list(loss = loss, alpha = alpha1, beta = beta1, Delta = Delta) + list( + loss = loss, + alpha = alpha1, beta = beta1, + Delta.1 = Delta.1, Delta.2 = Delta.2, tr.Delta = s, + break.reason = break.reason + ) } diff --git a/tensorPredictors/R/kpir_kron.R b/tensorPredictors/R/kpir_kron.R index 7b05985..b8be6db 100644 --- a/tensorPredictors/R/kpir_kron.R +++ b/tensorPredictors/R/kpir_kron.R @@ -126,8 +126,8 @@ kpir.kron <- function(X, Fy, shape = c(dim(X)[-1], dim(Fy[-1])), resid.trans <- resid %x_3% solve(Delta.1) %x_2% solve(Delta.2) # Calculate Gradients - grad.alpha <- tcrossprod(mat(resid.trans, 3), mat(Fy %x_2% beta, 3)) - grad.beta <- tcrossprod(mat(resid.trans, 2), mat(Fy %x_3% alpha, 2)) + grad.alpha <- tcrossprod(mat(resid.trans, 3), mat(Fy %x_2% S.beta, 3)) + grad.beta <- tcrossprod(mat(resid.trans, 2), mat(Fy %x_3% S.alpha, 2)) # Backtracking line search (Armijo type) # The `inner.prod` is used in the Armijo break condition but does not @@ -161,7 +161,7 @@ kpir.kron <- function(X, Fy, shape = c(dim(X)[-1], dim(Fy[-1])), } } - # Call logger (invoce history callback) + # Call logger (invoke history callback) if (is.function(logger)) { logger(iter, loss.temp, alpha.temp, beta.temp, Delta.1, Delta.2, delta) } @@ -207,5 +207,10 @@ kpir.kron <- function(X, Fy, shape = c(dim(X)[-1], dim(Fy[-1])), } - list(loss = loss, alpha = alpha1, beta = beta1, Delta = Delta, break.reason = break.reason) + list( + loss = loss, + alpha = alpha1, beta = beta1, + Delta.1 = Delta.1, Delta.2 = Delta.2, + break.reason = break.reason + ) } diff --git a/tensorPredictors/src/poi.c b/tensorPredictors/src/poi.c index de3dbc1..e175b4c 100644 --- a/tensorPredictors/src/poi.c +++ b/tensorPredictors/src/poi.c @@ -80,15 +80,3 @@ extern SEXP FastPOI_C_sub(SEXP in_B, SEXP in_Delta, UNPROTECT(1); return out_Z; } - -/* List of registered routines (e.g. C entry points) */ -static const R_CallMethodDef CallEntries[] = { - {"FastPOI_C_sub", (DL_FUNC) &FastPOI_C_sub, 5}, - {NULL, NULL, 0} -}; - -/* Restrict C entry points to registered routines. */ -void R_init_tensorPredictors(DllInfo *dll) { - R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); - R_useDynamicSymbols(dll, FALSE); -}