Fix prestogp_fit stability and beta.hat validation

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: PrestoGP
 Type: Package
 Title: Penalized Regression for Spatio-Temporal Outcomes via Gaussian Processes 
-Version: 0.2.0.9050
+Version: 0.2.0.9053
 Authors@R: c(
     person(given = "Eric",
            family = "Bair",
@@ -21,7 +21,7 @@ Authors@R: c(
 	   role = "aut"))
 Description: Simultaneous variable seletion and estimation of LUR models with spatiotemporally correlated errors that is scalable for big data.
 Depends:
-    R (>= 3.5.0)
+    R (>= 4.1.0)
 LinkingTo:
     Rcpp, RcppArmadillo, BH
 Imports:
@@ -67,3 +67,4 @@ Suggests:
     testthat (>= 3.0.0)
 Config/testthat/edition: 3
 URL: https://niehs.github.io/PrestoGP/, https://github.com/NIEHS/PrestoGP
+Config/roxygen2/version: 8.0.0

NEWS.md

@@ -1,3 +1,28 @@
+# PrestoGP 0.2.0.9053 (2026-5-26)
+
+## BUG FIXES
+
+* Fixed a bug where specifying the beta.hat parameter to `prestogp_fit`
+  caused an error due to a faulty length check.
+
+# PrestoGP 0.2.0.9052 (2026-1-20)
+
+## BUG FIXES
+
+* Created local versions of several GPvecchia utility functions (e.g.,
+  `vecchia_likelihood`, `vecchia_prediction`) that are robust to numerically
+  singular U/V matrices.
+
+* Modified `vecchia_prediction` and `vecchia_Mprediction` to avoid crashes
+  when the estimated V matrix is numerically singular.
+
+# PrestoGP 0.2.0.9051 (2026-1-14)
+
+## NEW FEATURES
+
+* Added a parameter `omp.cores` to `prestogp_fit` to specify the number of
+  cores used by OMP when computing the U matrix.
+
 # PrestoGP 0.2.0.9050 (2025-11-13)
 
 ## BUG FIXES

R/Log_Likelihood.R

@@ -59,6 +59,20 @@ negloglik.full <- function(logparams, d, y, param.seq) {
   -1 * mvtnorm::dmvnorm(y, rep(0, N), cov.mat, log = TRUE)
 }
 
+vecchia_likelihood <- function(z, vecchia.approx, covparms, nuggets,
+  covmodel = "matern") {
+  if (vecchia.approx$cond.yz == "zy")
+    warning("cond.yz='zy' will produce a poor likelihood approximation. Use 'SGV' instead.")
+  removeNAs <- getFromNamespace("removeNAs", "GPvecchia")
+  removeNAs()
+  U.obj <- createU(vecchia.approx, covparms, nuggets, covmodel)
+  res <- try(junk <- vecchia_likelihood_U(z, U.obj), silent = TRUE)
+  if (inherits(res, "try-error")) {
+    junk <- -1 * Inf
+  }
+  junk
+}
+
 #' Evaluation of the multivariate Vecchia likelihood
 #'
 #' This function is used to evaluate the multivariate Vecchia likelihood.
@@ -118,6 +132,30 @@ vecchia_Mlikelihood <- function(z, vecchia.approx, covparams) {
 }
 
 
+vecchia_likelihood_U <- function(z, U.obj) {
+  U <- U.obj$U
+  latent <- U.obj$latent
+  zord <- z[U.obj$ord.z]
+  const <- sum(!latent) * log(2 * pi)
+  z1 <- Matrix::crossprod(U[!latent, ], zord)
+  quadform.num <- sum(z1^2)
+  logdet.num <- -2 * sum(log(Matrix::diag(U)))
+  if (sum(latent) == 0) {
+    logdet.denom <- quadform.denom <- 0
+  } else {
+    U.y <- U[latent, ]
+    z2 <- as.numeric(U.y %*% z1)
+    V.ord <- U2V(U.obj)
+    z3 <- Matrix::solve(V.ord, rev(z2), system = "L")
+    quadform.denom <- sum(z3^2)
+    logdet.denom <- -2 * sum(log(Matrix::diag(V.ord)))
+  }
+  neg2loglik <- logdet.num - logdet.denom + quadform.num - quadform.denom +
+    const
+  loglik <- -neg2loglik / 2
+  loglik
+}
+
 ##############################################################################
 ### Flexible Multivariate Matern Negative Loglikelihood Function ###########
 

R/PrestoGP_CreateU_Multivariate.R

@@ -103,7 +103,7 @@ create_param_sequence <- function(P, ns = 1) {
   }
 
   param.sequence.begin <- c(1, P + 1, seq(P * (ns + 1) + 1, length = 3, by = P))
-  param.sequence.end <- c(P, ns * P, P, P, nk) %>% cumsum()
+  param.sequence.end <- c(P, ns * P, P, P, nk) |> cumsum()
   param.sequence <- cbind(param.sequence.begin, param.sequence.end)
 
   param.sequence
@@ -470,7 +470,7 @@ vecchia_Mspecify <- function(locs.list, m, locs.list.pred = NULL,
     ord.pred = ordering.pred, cond.yz = "SGV", conditioning = "NN",
     P = P, ondx = ondx, dist.func = dist.func,
     dist.func.code = dist.func.code, q.list = q.list,
-    n.neighbors = m
+    n.neighbors = m, n.cores = -1
   )
 }
 
@@ -609,7 +609,8 @@ createUMultivariate <- function(vec.approx, params, cov_func = NULL) {
     cur.qys <- do.call(cbind, q.list$q.y)
     cur.qzs <- do.call(cbind, q.list$q.z)
     # browser()
-    U <- createU_helper_mat(olocs, ondx, cur.qys, cur.qzs, vijs, aijs, full_const, nugget, sig2, uvec)
+    U <- createU_helper_mat(olocs, ondx, cur.qys, cur.qzs, vijs, aijs,
+      full_const, nugget, sig2, uvec, vec.approx$n.cores)
     # U <- sparseMatrix(i=U1[1,], j=U1[2,], x = U1[3,], triangular = TRUE)
   }
 

R/PrestoGP_Model.R

@@ -49,6 +49,8 @@
 #' @slot nscale The number of scale parameters in the model.
 #' @slot common_scale Do all columns of locs have the same scale parameter?
 #' @slot param_sequence Records the indices of the various Matern parameters.
+#' @slot omp_cores Number of cores used by OMP to compute the U matrix.
+#' @slot logparams Transformed version of the Matern parameters.
 #' See \code{\link{create_param_sequence}}.
 #'
 #' @seealso \code{\link{VecchiaModel-class}}, \code{\link{FullModel-class}},
@@ -97,6 +99,7 @@ PrestoGPModel <- setClass("PrestoGPModel",
     nscale = "numeric", # the number of scale parameters
     common_scale = "logical", # is there a common scale parameter?
     param_sequence = "matrix", # maps the indices of the various Matern parameters
+    omp_cores = "numeric", # Number of cores used by OMP to compute the U matrix
     logparams = "numeric" # transformed version of the Matern parameters
   )
 )
@@ -171,7 +174,7 @@ setGeneric(
     optim.control = list(trace = 0, reltol = 1e-3, maxit = 5000),
     penalty = c("lasso", "relaxed", "MCP", "SCAD"), alpha = 1,
     family = c("gaussian", "binomial"), nfolds = 10, foldid = NULL,
-    parallel = FALSE, cluster = NULL, adaptive = FALSE) {
+    omp.cores = -1, parallel = FALSE, cluster = NULL, adaptive = FALSE) {
     standardGeneric("prestogp_fit")
   }
 )
@@ -827,13 +830,19 @@ setMethod(
 #' what fold each observation should be assigned to in the cv.glmnet
 #' cross-validation. See \code{\link[glmnet]{cv.glmnet}} and
 #' \code{\link[ncvreg]{cv.ncvreg}}.
+#' @param omp.cores Number of parallel cores to be used by OMP when computing
+#' the U matrix in the Vecchia approximation. Defaults to -1 (which causes
+#' all available cores to be used).
 #' @param parallel Should parallel "foreach" be used to speed up the model
 #' fitting procedure where possible? Defaults to FALSE. Specifically,
 #' parallelization will be used for imputation and fitting the cv.glmnet
 #' object. See \code{\link[glmnet]{cv.glmnet}}. Note that this
 #' only applies to glmnet models (where penalty="lasso" or
 #' penalty="relaxed"). Models using ncvreg (where penalty="MCP" or
 #' penalty="SCAD") require a cluster argument for parallelization (see below).
+#' Note that parallelization will be used to compute the Vecchia
+#' approximation even if this parameter is set to FALSE (assuming that OMP is
+#' available on the system). To disable this, set omp.cores=1.
 #' @param cluster A cluster for running cv.ncvreg in parallel. See
 #' \code{\link[ncvreg]{cv.ncvreg}} and \code{\link[parallel]{makeCluster}}.
 #' This must be specified to run cv.ncvreg in parallel. It is ignored for
@@ -947,9 +956,19 @@ setMethod(
     quiet = FALSE, verbose = FALSE, optim.method = "Nelder-Mead",
     optim.control = list(trace = 0, reltol = 1e-3, maxit = 5000),
     penalty = c("lasso", "relaxed", "MCP", "SCAD"), alpha = 1,
-    family = c("gaussian", "binomial"),
-    nfolds = 10, foldid = NULL, parallel = FALSE, cluster = NULL,
+    family = c("gaussian", "binomial"), nfolds = 10, foldid = NULL,
+    omp.cores = -1, parallel = FALSE, cluster = NULL,
     adaptive = FALSE) {
+    if (!is.numeric(omp.cores)) {
+      stop("omp.cores must be a positive integer >= 1")
+    }
+    if (length(omp.cores) != 1) {
+      stop("omp.cores must be a positive integer >= 1")
+    }
+    if (omp.cores < -1 | omp.cores == 0 | omp.cores %% 1 != 0) {
+      stop("omp.cores must be a positive integer >= 1")
+    }
+    model@omp_cores <- omp.cores
     penalty <- match.arg(penalty)
     model@penalty <- penalty
     family <- match.arg(family)
@@ -994,8 +1013,7 @@ setMethod(
       if (!is.vector(beta.hat) | !is.numeric(beta.hat)) {
         stop("beta.hat parameter must be a numeric vector")
       }
-      if (length(beta.hat) != (ncol(model@X_train) +
-            length(model@locs_train))) {
+      if (length(beta.hat) != (ncol(model@X_train) + 1)) {
         stop("Length of beta.hat must match the number of predictors")
       }
       beta.hat <- as.matrix(beta.hat)

R/PrestoGP_Multivariate_Vecchia.R

@@ -542,6 +542,7 @@ setMethod("specify", "MultivariateVecchiaModel", function(model) {
   locs <- model@locs_train
   locs.scaled <- scale_locs(model, locs)
   model@vecchia_approx <- vecchia_Mspecify(locs.scaled, model@n_neighbors)
+  model@vecchia_approx$n.cores <- model@omp_cores
   if (!model@common_scale) {
     olocs.scaled <- model@vecchia_approx$locsord
     for (i in seq_along(locs)) {

R/PrestoGP_Util_Functions.R

@@ -226,6 +226,44 @@ revMat <- function(mat) {
   mat.out
 }
 
+vecchia_prediction <- function(z, vecchia.approx, covparms, nuggets,
+  var.exact, covmodel = "matern", return.values = "all") {
+  removeNAs <- getFromNamespace("removeNAs", "GPvecchia")
+  removeNAs()
+  U.obj <- createU(vecchia.approx, covparms, nuggets, covmodel)
+  V.ord <- U2V(U.obj)
+  if (length(U.obj$zero.nugg) > 0)
+    warning("Rows/cols of V have been removed for data with zero noise")
+  vecchia_mean <- getFromNamespace("vecchia_mean", "GPvecchia")
+  V.singular <- FALSE
+  res <- try(vecchia.mean <- vecchia_mean(z, U.obj, V.ord), silent = TRUE)
+  if (inherits(res, "try-error")) {
+    warning("V is numerically singular. Predicted means are unreliable.")
+    V.ord.pd <- Matrix::nearPD(V.ord)$mat
+    V.singular <- TRUE
+    vecchia.mean <- vecchia_mean(z, U.obj, V.ord.pd)
+  }
+  return.list <- list(mu.pred = vecchia.mean$mu.pred,
+    mu.obs = vecchia.mean$mu.obs, var.pred = NULL, var.obs = NULL,
+    V.ord = NULL, U.obj = NULL)
+  if (return.values == "meanmat" || return.values == "all") {
+    return.list$V.ord <- V.ord
+    return.list$U.obj <- U.obj
+  }
+  if (return.values == "meanvar" || return.values == "all") {
+    if (V.singular) {
+      stop("V is numerically singular. Prediction variance cannot be computed.")
+    }
+    if (missing(var.exact))
+      var.exact <- (sum(!vecchia.approx$obs) < 4 * 10000)
+    vecchia_var <- getFromNamespace("vecchia_var", "GPvecchia")
+    vars.vecchia <- vecchia_var(U.obj, V.ord, exact = var.exact)
+    return.list$var.pred <- vars.vecchia$vars.pred
+    return.list$var.obs <- vars.vecchia$vars.obs
+  }
+  return.list
+}
+
 #' Multivariate Vecchia prediction
 #'
 #' This function is used to make predictions based on multivariate Vecchia
@@ -296,7 +334,14 @@ vecchia_Mprediction <- function(z, vecchia.approx, covparms, var.exact = NULL, r
   #    if (length(U.obj$zero.nugg) > 0)
   #        warning("Rows/cols of V have been removed for data with zero noise")
   vecchia_mean <- getFromNamespace("vecchia_mean", "GPvecchia")
-  vecchia.mean <- vecchia_mean(z, U.obj, V.ord)
+  V.singular <- FALSE
+  res <- try(vecchia.mean <- vecchia_mean(z, U.obj, V.ord), silent = TRUE)
+  if (inherits(res, "try-error")) {
+    warning("V is numerically singular. Predicted means are unreliable.")
+    V.ord.pd <- Matrix::nearPD(V.ord)$mat
+    V.singular <- TRUE
+    vecchia.mean <- vecchia_mean(z, U.obj, V.ord.pd)
+  }
   return.list <- list(
     mu.pred = vecchia.mean$mu.pred, mu.obs = vecchia.mean$mu.obs,
     var.pred = NULL, var.obs = NULL, V.ord = NULL, U.obj = NULL
@@ -306,6 +351,9 @@ vecchia_Mprediction <- function(z, vecchia.approx, covparms, var.exact = NULL, r
     return.list$U.obj <- U.obj
   }
   if (return.values == "meanvar" || return.values == "all") {
+    if (V.singular) {
+      stop("V is numerically singular. Prediction variance cannot be computed.")
+    }
     if (is.null(var.exact)) {
       var.exact <- (sum(!vecchia.approx$obs) < 4 * 10000)
     }

R/PrestoGP_Vecchia.R

@@ -383,6 +383,9 @@ setMethod("impute_y_lod", "VecchiaModel", function(model, lodu, lodl, n.mi = 10,
 setMethod("specify", "VecchiaModel", function(model) {
   locs.scaled <- scale_locs(model, model@locs_train)
   model@vecchia_approx <- vecchia_specify(locs.scaled[[1]], model@n_neighbors)
+  if (model@omp_cores > -1) {
+    model@vecchia_approx$U.prep$n.cores <- model@omp_cores
+  }
   if (!model@common_scale) {
     olocs.scaled <- model@vecchia_approx$locsord
     for (j in 1:model@nscale) {

R/RcppExports.R

@@ -5,7 +5,7 @@ na_omit_c <- function(x) {
     .Call('_PrestoGP_na_omit_c', PACKAGE = 'PrestoGP', x)
 }
 
-createU_helper_mat <- function(olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning) {
-    .Call('_PrestoGP_createU_helper_mat', PACKAGE = 'PrestoGP', olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning)
+createU_helper_mat <- function(olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning, n_cores) {
+    .Call('_PrestoGP_createU_helper_mat', PACKAGE = 'PrestoGP', olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning, n_cores)
 }
 

src/RcppExports.cpp

@@ -23,8 +23,8 @@ BEGIN_RCPP
 END_RCPP
 }
 // createU_helper_mat
-arma::sp_mat createU_helper_mat(const arma::mat& olocs, const arma::vec& ondx, const arma::mat& curqys, const arma::mat& curqzs, const arma::mat& vijs, const arma::mat& aijs, const arma::mat& full_const, const arma::vec& nugget, const arma::vec& sig2, const arma::vec& U_beginning);
-RcppExport SEXP _PrestoGP_createU_helper_mat(SEXP olocsSEXP, SEXP ondxSEXP, SEXP curqysSEXP, SEXP curqzsSEXP, SEXP vijsSEXP, SEXP aijsSEXP, SEXP full_constSEXP, SEXP nuggetSEXP, SEXP sig2SEXP, SEXP U_beginningSEXP) {
+arma::sp_mat createU_helper_mat(const arma::mat& olocs, const arma::vec& ondx, const arma::mat& curqys, const arma::mat& curqzs, const arma::mat& vijs, const arma::mat& aijs, const arma::mat& full_const, const arma::vec& nugget, const arma::vec& sig2, const arma::vec& U_beginning, const int n_cores);
+RcppExport SEXP _PrestoGP_createU_helper_mat(SEXP olocsSEXP, SEXP ondxSEXP, SEXP curqysSEXP, SEXP curqzsSEXP, SEXP vijsSEXP, SEXP aijsSEXP, SEXP full_constSEXP, SEXP nuggetSEXP, SEXP sig2SEXP, SEXP U_beginningSEXP, SEXP n_coresSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -38,14 +38,15 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< const arma::vec& >::type nugget(nuggetSEXP);
     Rcpp::traits::input_parameter< const arma::vec& >::type sig2(sig2SEXP);
     Rcpp::traits::input_parameter< const arma::vec& >::type U_beginning(U_beginningSEXP);
-    rcpp_result_gen = Rcpp::wrap(createU_helper_mat(olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning));
+    Rcpp::traits::input_parameter< const int >::type n_cores(n_coresSEXP);
+    rcpp_result_gen = Rcpp::wrap(createU_helper_mat(olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning, n_cores));
     return rcpp_result_gen;
 END_RCPP
 }
 
 static const R_CallMethodDef CallEntries[] = {
     {"_PrestoGP_na_omit_c", (DL_FUNC) &_PrestoGP_na_omit_c, 1},
-    {"_PrestoGP_createU_helper_mat", (DL_FUNC) &_PrestoGP_createU_helper_mat, 10},
+    {"_PrestoGP_createU_helper_mat", (DL_FUNC) &_PrestoGP_createU_helper_mat, 11},
     {NULL, NULL, 0}
 };