From 06be69a2b7478910560a49dee33890cb39ed4c37 Mon Sep 17 00:00:00 2001 From: Julien Chiquet Date: Wed, 22 Sep 2021 17:01:39 +0200 Subject: [PATCH] Now use n * p variational parameters for the spherical covariance case, fixing #79 --- R/PLNLDAfit-class.R | 2 +- R/PLNfit-class.R | 5 +- R/PLNmixturefamily-class.R | 6 +-- R/PLNmixturefit-class.R | 5 +- inst/case_studies/oaks_tree.R | 6 ++- src/optimize.cpp | 96 ++--------------------------------- src/optimize_ve.cpp | 24 ++++----- 7 files changed, 25 insertions(+), 119 deletions(-) diff --git a/R/PLNLDAfit-class.R b/R/PLNLDAfit-class.R index dfb61946..08bb7283 100644 --- a/R/PLNLDAfit-class.R +++ b/R/PLNLDAfit-class.R @@ -92,7 +92,7 @@ PLNLDAfit <- R6Class( covariates <- cbind(covariates, model.matrix( ~ grouping + 0)) super$postTreatment(responses, covariates, offsets) rownames(private$B) <- colnames(private$B) <- colnames(responses) - if (private$covariance != "spherical") colnames(private$S2) <- 1:self$q + colnames(private$S2) <- 1:self$q self$setVisualization() }, diff --git a/R/PLNfit-class.R b/R/PLNfit-class.R index 9ac28df5..3d6a8e2e 100644 --- a/R/PLNfit-class.R +++ b/R/PLNfit-class.R @@ -97,7 +97,7 @@ PLNfit <- R6Class( private$Theta <- do.call(rbind, lapply(LMs, coefficients)) residuals <- do.call(cbind, lapply(LMs, residuals)) private$M <- residuals - private$S2 <- matrix(0.1, n, ifelse(control$covariance == "spherical", 1, p)) + private$S2 <- matrix(0.1,n,p) if (control$covariance == "spherical") { private$Sigma <- diag(sum(residuals^2)/(n*p), p, p) } else if (control$covariance == "diagonal") { @@ -173,8 +173,7 @@ PLNfit <- R6Class( ## Initialize the variational parameters with the appropriate new dimension of the data optim_out <- VEstep_optimizer( - list(M = matrix(0, n, p), - S = matrix(sqrt(0.1), n, ifelse(self$vcov_model == "spherical", 1, p))), + list(M = matrix(0, n, p), S = matrix(sqrt(0.1), n, p)), responses, covariates, offsets, diff --git a/R/PLNmixturefamily-class.R b/R/PLNmixturefamily-class.R index fddb70db..0a3c4836 100644 --- a/R/PLNmixturefamily-class.R +++ b/R/PLNmixturefamily-class.R @@ -126,11 +126,7 @@ PLNmixturefamily <- myPLN <- PLNfit$new(responses, covariates, offsets, rep(1, nrow(responses)), formula, xlevels, control) myPLN$optimize(responses, covariates, offsets, rep(1, nrow(responses)), control) - if(control$covariance == 'spherical') - Sbar <- c(myPLN$var_par$S2) * myPLN$p - else - Sbar <- rowSums(myPLN$var_par$S2) - + Sbar <- rowSums(myPLN$var_par$S2) D <- sqrt(as.matrix(dist(myPLN$var_par$M)^2) + outer(Sbar,rep(1,myPLN$n)) + outer(rep(1, myPLN$n), Sbar)) if (is.numeric(control$init_cl)) { diff --git a/R/PLNmixturefit-class.R b/R/PLNmixturefit-class.R index 49bd02b3..c190c085 100644 --- a/R/PLNmixturefit-class.R +++ b/R/PLNmixturefit-class.R @@ -46,7 +46,6 @@ PLNmixturefit <- M <- private$comp %>% map("var_par") %>% map("M") S2 <- private$comp %>% map("var_par") %>% map("S2") - if(private$covariance == "spherical") S2 <- map(S2, ~outer(as.numeric(.x), rep(1, self$p) )) mu <- private$comp %>% map(coef) %>% map(~outer(rep(1, self$n), as.numeric(.x))) @@ -252,7 +251,7 @@ PLNmixturefit <- #' @return a [`ggplot`] graphic plot_clustering_data = function(main = "Expected counts reorder by clustering", plot = TRUE, log_scale = TRUE) { M <- private$mix_up('var_par$M') - S2 <- switch(private$covariance, "spherical" = private$mix_up('var_par$S2') %*% rbind(rep(1, ncol(M))), private$mix_up('var_par$S2')) + S2 <- private$mix_up('var_par$S2') mu <- self$posteriorProb %*% t(self$group_means) A <- exp(mu + M + .5 * S2) p <- plot_matrix(A, 'samples', 'variables', self$memberships, log_scale) @@ -346,7 +345,7 @@ PLNmixturefit <- #' @field entropy_latent Entropy of the variational distribution of the latent vector (Gaussian) entropy_latent = function() { .5 * (sum(map_dbl(private$comp, function(component) { - sum( diag(component$weights) %*% log(component$var_par$S2) * ifelse(component$vcov_model == "spherical", self$p, 1) ) + sum( diag(component$weights) %*% log(component$var_par$S2) ) })) + self$n * self$p * log(2*pi*exp(1))) }, #' @field entropy Full entropy of the variational distribution (latent vector + clustering) diff --git a/inst/case_studies/oaks_tree.R b/inst/case_studies/oaks_tree.R index 4be679db..7d220f8a 100644 --- a/inst/case_studies/oaks_tree.R +++ b/inst/case_studies/oaks_tree.R @@ -32,6 +32,9 @@ plot(myLDA_tree_diagonal) otu.family <- factor(rep(c("fungi", "E. aphiltoides", "bacteria"), c(47, 1, 66))) plot(myLDA_tree, "variable", var_cols = otu.family) ## TODO: add color for arrows to check +myLDA_tree_spherical <- PLNLDA(Abundance ~ 1 + offset(log(Offset)), grouping = tree, data = oaks, control = list(covariance = "spherical")) +plot(myLDA_tree_spherical) + ## One dimensional check of plot myLDA_orientation <- PLNLDA(Abundance ~ 1 + offset(log(Offset)), grouping = orientation, data = oaks) plot(myLDA_orientation) @@ -91,8 +94,7 @@ data.frame( ggplot(aes(x = nb_components, y = value, colour = score)) + geom_line() + theme_bw() + labs(y = "clustering similarity", x = "number of components") ## Mixture model to recover tree structure - with covariates -system.time(my_mixtures <- PLNmixture(Abundance ~ 0 + tree + distTOground + offset(log(Offset)), data = oaks, - control_main = list(covariance = "spherical"))) +system.time(my_mixtures <- PLNmixture(Abundance ~ 0 + tree + distTOground + offset(log(Offset)), data = oaks)) plot(my_mixtures, criteria = c("loglik", "ICL", "BIC"), reverse = TRUE) diff --git a/src/optimize.cpp b/src/optimize.cpp index 3838bc93..e35d37ea 100644 --- a/src/optimize.cpp +++ b/src/optimize.cpp @@ -123,97 +123,7 @@ Rcpp::List cpp_optimize_spherical( // Conversion from R, prepare optimization const auto init_Theta = Rcpp::as(init_parameters["Theta"]); // (p,d) const auto init_M = Rcpp::as(init_parameters["M"]); // (n,p) - const auto init_S = Rcpp::as(init_parameters["S"]); // (n) - - const auto metadata = tuple_metadata(init_Theta, init_M, init_S); - enum { THETA_ID, M_ID, S_ID }; // Names for metadata indexes - - auto parameters = std::vector(metadata.packed_size); - metadata.map(parameters.data()) = init_Theta; - metadata.map(parameters.data()) = init_M; - metadata.map(parameters.data()) = init_S; - - auto optimizer = new_nlopt_optimizer(configuration, parameters.size()); - if(configuration.containsElementNamed("xtol_abs")) { - SEXP value = configuration["xtol_abs"]; - if(Rcpp::is(value)) { - set_uniform_xtol_abs(optimizer.get(), Rcpp::as(value)); - } else { - auto per_param_list = Rcpp::as(value); - auto packed = std::vector(metadata.packed_size); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["Theta"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["M"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["S"]); - set_per_value_xtol_abs(optimizer.get(), packed); - } - } - - const double w_bar = accu(w); - - // Optimize - auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &w_bar](const double * params, double * grad) -> double { - const arma::mat Theta = metadata.map(params); - const arma::mat M = metadata.map(params); - const arma::mat S = metadata.map(params); - - arma::mat S2 = S % S; - const arma::uword p = Y.n_cols; - arma::mat Z = O + X * Theta.t() + M; - arma::mat A = exp(Z + 0.5 * S2 * arma::ones(p).t()); - double sigma2 = arma::as_scalar(dot(w, sum(pow(M, 2), 1) + double(p) * S2) / (double(p) * w_bar) ); - double objective = accu(w.t() * (A - Y % Z)) - 0.5 * double(p) * dot(w, log(S2/sigma2)) ; - - metadata.map(grad) = (A - Y).t() * (X.each_col() % w); - metadata.map(grad) = diagmat(w) * (M / sigma2 + A - Y); - metadata.map(grad) = w % (double(p) * S / sigma2 + S % sum(A, 1) - double(p) * pow(S, -1)); - - return objective; - }; - OptimizerResult result = minimize_objective_on_parameters(optimizer.get(), objective_and_grad, parameters); - - // Variational parameters - arma::mat M = metadata.copy(parameters.data()); - arma::mat S = metadata.copy(parameters.data()); // vec(n) -> mat(n, 1) - arma::mat S2 = S % S; - // Regression parameters - arma::mat Theta = metadata.copy(parameters.data()); - // Variance parameters - const arma::uword p = Y.n_cols; - const double sigma2 = arma::as_scalar(dot(w, sum(pow(M, 2), 1) + double(p) * S2)) / (double(p) * w_bar); - arma::sp_mat Sigma(p,p); Sigma.diag() = arma::ones(p) * sigma2; - arma::sp_mat Omega(p,p); Omega.diag() = arma::ones(p) * pow(sigma2, -1); - // Element-wise log-likelihood - arma::mat Z = O + X * Theta.t() + M; - arma::mat A = exp(Z + 0.5 * S2 * arma::ones(p).t()); - arma::mat loglik = sum(Y % Z - A - 0.5 * pow(M, 2) / sigma2, 1) - 0.5 * double(p) * S2 / sigma2 + - 0.5 * double(p) * log(S2 / sigma2) + ki(Y); - - return Rcpp::List::create( - Rcpp::Named("status", static_cast(result.status)), - Rcpp::Named("iterations", result.nb_iterations), - Rcpp::Named("Theta", Theta), - Rcpp::Named("M", M), - Rcpp::Named("S", S), - Rcpp::Named("Z", Z), - Rcpp::Named("A", A), - Rcpp::Named("Sigma", Sigma), - Rcpp::Named("Omega", Omega), - Rcpp::Named("loglik", loglik)); -} - -// [[Rcpp::export]] -Rcpp::List cpp_optimize_spherical_2( - const Rcpp::List & init_parameters, // List(Theta, M, S) - const arma::mat & Y, // responses (n,p) - const arma::mat & X, // covariates (n,d) - const arma::mat & O, // offsets (n,p) - const arma::vec & w, // weights (n) - const Rcpp::List & configuration // List of config values -) { - // Conversion from R, prepare optimization - const auto init_Theta = Rcpp::as(init_parameters["Theta"]); // (p,d) - const auto init_M = Rcpp::as(init_parameters["M"]); // (n,p) - const auto init_S = Rcpp::as(init_parameters["S"]); // (n,p) + const auto init_S = Rcpp::as(init_parameters["S"]); // (n,p) const auto metadata = tuple_metadata(init_Theta, init_M, init_S); enum { THETA_ID, M_ID, S_ID }; // Names for metadata indexes @@ -250,7 +160,7 @@ Rcpp::List cpp_optimize_spherical_2( const arma::uword p = Y.n_cols; arma::mat Z = O + X * Theta.t() + M; arma::mat A = exp(Z + 0.5 * S2); - double sigma2 = arma::as_scalar(dot(w, sum(pow(M, 2) + S2, 1)) / (double(p) * w_bar) ); + double sigma2 = accu(diagmat(w) * (pow(M, 2) + S2)) / (double(p) * w_bar) ; double objective = accu(w.t() * (A - Y % Z - 0.5 * log(S2))) - 0.5 * (double(p) * w_bar) * log(sigma2) ; metadata.map(grad) = (A - Y).t() * (X.each_col() % w); @@ -269,7 +179,7 @@ Rcpp::List cpp_optimize_spherical_2( arma::mat Theta = metadata.copy(parameters.data()); // Variance parameters const arma::uword p = Y.n_cols; - const double sigma2 = arma::as_scalar(dot(w, sum(pow(M, 2) + S2, 1))) / (double(p) * w_bar); + const double sigma2 = accu(diagmat(w) * (pow(M, 2) + S2)) / (double(p) * w_bar) ; arma::sp_mat Sigma(p,p); Sigma.diag() = arma::ones(p) * sigma2; arma::sp_mat Omega(p,p); Omega.diag() = arma::ones(p) * pow(sigma2, -1); // Element-wise log-likelihood diff --git a/src/optimize_ve.cpp b/src/optimize_ve.cpp index ed6fe8bd..082d37ee 100644 --- a/src/optimize_ve.cpp +++ b/src/optimize_ve.cpp @@ -185,7 +185,7 @@ Rcpp::List cpp_optimize_vestep_spherical( ) { // Conversion from R, prepare optimization const auto init_M = Rcpp::as(init_parameters["M"]); // (n,p) - const auto init_S = Rcpp::as(init_parameters["S"]); // (n) + const auto init_S = Rcpp::as(init_parameters["S"]); // (n) const auto metadata = tuple_metadata(init_M, init_S); enum { M_ID, S_ID }; // Names for metadata indexes @@ -214,31 +214,31 @@ Rcpp::List cpp_optimize_vestep_spherical( const arma::mat M = metadata.map(params); const arma::mat S = metadata.map(params); - arma::vec S2 = S % S; + arma::mat S2 = S % S; const arma::uword p = Y.n_cols; arma::mat Z = O + X * Theta.t() + M; - arma::mat A = exp(Z + 0.5 * S2 * arma::ones(p).t()); - double n_sigma2 = dot(w, sum(pow(M, 2), 1) + double(p) * S2); + arma::mat A = exp(Z + 0.5 * S2); + double n_sigma2 = accu(diagmat(w) * (pow(M, 2) + S2)) ; double omega2 = Omega(0, 0); - double objective = accu(w.t() * (A - Y % Z)) - 0.5 * double(p) * dot(w, log(S2)) + 0.5 * n_sigma2 * omega2; + double objective = accu(w.t() * (A - Y % Z - 0.5 * log(S2))) + 0.5 * n_sigma2 * omega2; + + metadata.map(grad) = diagmat(w) * (M / omega2 + A - Y); + metadata.map(grad) = diagmat(w) * (S / omega2 + S % A - pow(S, -1)); - metadata.map(grad) = diagmat(w) * (M * omega2 + A - Y); - metadata.map(grad) = w % (double(p) * S * omega2 + S % sum(A, 1) - double(p) * pow(S, -1)); return objective; }; OptimizerResult result = minimize_objective_on_parameters(optimizer.get(), objective_and_grad, parameters); // Model and variational parameters arma::mat M = metadata.copy(parameters.data()); - arma::mat S = metadata.copy(parameters.data()); // vec(n) -> mat(n, 1) - arma::vec S2 = S % S; + arma::mat S = metadata.copy(parameters.data()); + arma::mat S2 = S % S; double omega2 = Omega(0, 0); // Element-wise log-likelihood const arma::uword p = Y.n_cols; arma::mat Z = O + X * Theta.t() + M; - arma::mat A = exp(Z + 0.5 * S2 * arma::ones(p).t()); - arma::mat loglik = sum(Y % Z - A - 0.5 * pow(M, 2) * omega2, 1) - 0.5 * double(p) * omega2 * S2 + - 0.5 * double(p) * log(S2 * omega2) + ki(Y); + arma::mat A = exp(Z + 0.5 * S2); + arma::mat loglik = sum(Y % Z - A - 0.5 * (pow(M, 2) + S2 ) * omega2 + 0.5 * log(S2 * omega2), 1) + ki(Y); return Rcpp::List::create( Rcpp::Named("status") = (int)result.status,