From 7129cc0e7c7b6049017a4b2c7bcf6fe9a4471911 Mon Sep 17 00:00:00 2001 From: Ryan O'Dea <70209371+ryan-odea@users.noreply.github.com> Date: Wed, 19 Feb 2025 01:04:20 -0500 Subject: [PATCH 01/16] kernel implementation --- DESCRIPTION | 9 +- NAMESPACE | 5 +- R/RcppExports.R | 11 ++ R/cv_learner.R | 156 ---------------------------- R/learner-package.R | 6 ++ R/learner.R | 223 ++++++++-------------------------------- man/cv.learner.Rd | 16 +-- man/learner.Rd | 70 ++++--------- src/.gitignore | 3 + src/Makevars | 3 + src/Makevars.win | 3 + src/RcppExports.cpp | 62 +++++++++++ src/kernels.cpp | 191 ++++++++++++++++++++++++++++++++++ src/learner-package.cpp | 2 + 14 files changed, 367 insertions(+), 393 deletions(-) create mode 100644 R/RcppExports.R delete mode 100644 R/cv_learner.R create mode 100644 R/learner-package.R create mode 100644 src/.gitignore create mode 100644 src/Makevars create mode 100644 src/Makevars.win create mode 100644 src/RcppExports.cpp create mode 100644 src/kernels.cpp create mode 100644 src/learner-package.cpp diff --git a/DESCRIPTION b/DESCRIPTION index f31d101..02bca85 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -8,7 +8,9 @@ Authors@R: c(person("Sean", "McGrath", role = c("aut", "cre"), person("Cenhao", "Zhu", role = c("aut"), email = "cenhaoz@mit.edu"), person("Rui", "Duan", role = c("aut"), - email = "rduan@hsph.harvard.edu")) + email = "rduan@hsph.harvard.edu"), + person("Ryan", "ODea", role = c("ctb"), + email = "ryanodea@hsph.harvard.edu")) Maintainer: Sean McGrath Description: Implements transfer learning methods for low-rank matrix estimation. These methods leverage similarity in the latent row and column @@ -22,7 +24,10 @@ LazyData: true RoxygenNote: 7.3.2 URL: https://github.com/stmcg/learner BugReports: https://github.com/stmcg/learner/issues -Imports: doParallel, foreach, ScreeNOT +Imports: ScreeNOT, + Rcpp (>= 1.0.11), + RcppEigen +LinkingTo: Rcpp, RcppEigen Suggests: testthat (>= 3.0.0) Config/testthat/edition: 3 diff --git a/NAMESPACE b/NAMESPACE index adc0742..1ab3843 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,5 +3,6 @@ export(cv.learner) export(dlearner) export(learner) -import(doParallel) -import(foreach) +import(RcppEigen) +importFrom(Rcpp,sourceCpp) +useDynLib(learner, .registration = TRUE) diff --git a/R/RcppExports.R b/R/RcppExports.R new file mode 100644 index 0000000..0c3b1d7 --- /dev/null +++ b/R/RcppExports.R @@ -0,0 +1,11 @@ +# Generated by using Rcpp::compileAttributes() -> do not edit by hand +# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +learner_cpp <- function(Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold) { + .Call(`_learner_learner_cpp`, Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold) +} + +cv_learner_cpp <- function(Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r) { + .Call(`_learner_cv_learner_cpp`, Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r) +} + diff --git a/R/cv_learner.R b/R/cv_learner.R deleted file mode 100644 index a16e93a..0000000 --- a/R/cv_learner.R +++ /dev/null @@ -1,156 +0,0 @@ -#' Cross-validation for LEARNER -#' -#' This function performs k-fold cross-validation to select the nuisance parameters \eqn{(\lambda_1, \lambda_2)} for \code{\link{learner}}. -#' -#' @param Y_target matrix containing the target population data, as in \code{\link{learner}} -#' @param Y_source matrix containing the source population data, as in \code{\link{learner}} -#' @param r (optional) integer specifying the rank of the knowledge graphs, as in \code{\link{learner}} -#' @param lambda_1_all vector of numerics specifying the candidate values of \eqn{\lambda_1} (see Details) -#' @param lambda_2_all vector of numerics specifying the candidate values of \eqn{\lambda_2} (see Details) -#' @param step_size numeric scalar specifying the step size for the Newton steps in the numerical optimization algorithm, as in \code{\link{learner}} -#' @param n_folds an integer specify the number of cross-validation folds. The default is \code{4}. -#' @param n_cores an integer specifying the number of CPU cores in parallelization. Parallelization is performed across the different candidate \eqn{(\lambda_1, \lambda_2)} pairs. The default is \code{1}, i.e., no parallelization. -#' @param control a list of parameters for controlling the stopping criteria for the numerical optimization algorithm, as in \code{\link{learner}}. -#' -#' @return A list with the following elements: -#' \item{lambda_1_min}{value of \eqn{\lambda_1} with the smallest MSE} -#' \item{lambda_2_min}{value of \eqn{\lambda_2} with the smallest MSE} -#' \item{mse_all}{matrix containing MSE value for each \eqn{(\lambda_1, \lambda_2)} pair. The rows correspond to the \eqn{\lambda_1} values, and the columns correspond to the \eqn{\lambda_2} values.} -#' \item{r}{rank value used.} -#' -#' @details -#' Given sets of candidate values of \eqn{\lambda_1} and \eqn{\lambda_2}, this function performs k-fold cross-validation to select the pair \eqn{(\lambda_1, \lambda_2)} with the smallest held out error. This function randomly partitions the entries of \code{Y_target} into \eqn{k} (approximately) equally sized subsamples. The training data sets are obtained by removing one of the \eqn{k} subsamples and the corresponding test data sets are based on the held out subsamples. The \code{\link{learner}} function is applied to each training data set. The held out error is computed by the mean squared error comparing the entries in the test data sets with those imputed from the LEARNER estimates. See McGrath et al. (2024) for further details. -#' -#' @references -#' McGrath, S., Zhu, C,. Guo, M. and Duan, R. (2024). \emph{LEARNER: A transfer learning method for low-rank matrix estimation}. arXiv preprint arXiv:2412.20605. -#' -#' @examples -#' res <- cv.learner(Y_source = dat_highsim$Y_source, -#' Y_target = dat_highsim$Y_target, -#' lambda_1_all = c(1, 10, 100), -#' lambda_2_all = c(1, 10, 100), -#' step_size = 0.003) -#' -#' -#' @import foreach -#' @import doParallel -#' -#' @export - - - -cv.learner <- function(Y_source, Y_target, r, lambda_1_all, lambda_2_all, - step_size, n_folds = 4, n_cores = 1, control = list()){ - - ## Error catching - if (!identical(dim(Y_source), dim(Y_target))){ - stop('Y_source and Y_target must have the same dimensions') - } - if (any(is.na(Y_source))){ - stop('Y_source cannot have NA values.') - } - - ## Setting the rank and other parameters - p <- nrow(Y_source) - q <- ncol(Y_source) - if (missing(r)){ - max_rank <- min(p, q) / 3 - r <- max(ScreeNOT::adaptiveHardThresholding(Y = Y_source, k = max_rank)$r, 1) - } - if (!('max_iter' %in% names(control))){ - control$max_iter <- 100 - } - if (!('threshold' %in% names(control))){ - control$threshold <- 0.001 - } - if (!('max_value' %in% names(control))){ - control$max_value <- 10 - } - - svd_source <- svd(Y_source, nu = r, nv = r) - P_U <- diag(rep(1, p)) - svd_source$u %*% t(svd_source$u) - P_V <- diag(rep(1, q)) - svd_source$v %*% t(svd_source$v) - U_init <- svd_source$u %*% diag(sqrt(svd_source$d[1:r]), nrow = r, ncol = r) - V_init <- svd_source$v %*% diag(sqrt(svd_source$d[1:r]), nrow = r, ncol = r) - - if (n_folds < 2){ - stop('n_folds must be 2 or greater.') - } - n_lambda_1 <- length(lambda_1_all) - n_lambda_2 <- length(lambda_2_all) - - ## Creating training and testing data sets - available_indices <- 1:(p * q) - available_indices <- available_indices[!(available_indices %in% which(is.na(Y_target)))] - n_indices <- length(available_indices) - indices <- sample(available_indices, size = n_indices, replace = FALSE) - - index_set <- vector(mode = "list", length = n_folds) - for (fold in 1:n_folds){ - index_set[[fold]] <- indices[floor((fold - 1) * n_indices / n_folds + 1):floor(fold * n_indices / n_folds)] - } - - ## Applying LEARNER for each candidate lambda_1 and lambda_2 values - if (n_cores > 1){ - ## Parallelization - registerDoParallel(cores = n_cores) - mse_all <- foreach(lambda_1_ind = 1:n_lambda_1, .combine = 'cbind') %:% - foreach(lambda_2_ind = 1:n_lambda_2, .combine = 'c') %dopar% { - cv.learner_helper(Y_target = Y_target, index_set = index_set, n_folds = n_folds, - lambda1 = lambda_1_all[lambda_1_ind], lambda2 = lambda_2_all[lambda_2_ind], - U_init = U_init, V_init = V_init, - P_U = P_U, P_V = P_V, step_size = step_size, - control = control) - } - } else { - ## Sequential - mse_all <- foreach(lambda_1_ind = 1:n_lambda_1, .combine = 'cbind') %:% - foreach(lambda_2_ind = 1:n_lambda_2, .combine = 'c') %do% { - cv.learner_helper(Y_target = Y_target, index_set = index_set, n_folds = n_folds, - lambda1 = lambda_1_all[lambda_1_ind], lambda2 = lambda_2_all[lambda_2_ind], - U_init = U_init, V_init = V_init, - P_U = P_U, P_V = P_V, step_size = step_size, - control = control) - } - } - mse_all <- t(mse_all) - - min_ind <- which(mse_all == min(mse_all, na.rm = TRUE), arr.ind = TRUE) - lambda_1_min <- lambda_1_all[min_ind[1]] - lambda_2_min <- lambda_2_all[min_ind[2]] - - return(list( - lambda_1_min = lambda_1_min, - lambda_2_min = lambda_2_min, - mse_all = mse_all, - r = r - )) -} - -cv.learner_helper <- function(Y_target, index_set, n_folds, - lambda1, lambda2, - U_init, V_init, P_U, P_V, - step_size, control){ - norm_temp <- 0 - p <- nrow(Y_target); q <- ncol(Y_target) - for (fold in 1:n_folds){ - Y_training <- Y_target - myset <- index_set[[fold]] - Y_training[myset] <- NA - perc_nonmissing <- sum(!is.na(Y_training)) / (p * q) - - fit <- try(learner_helper(theta_hat = Y_training, - lambda1 = lambda1, - lambda2 = lambda2, - U_init = U_init, V_init = V_init, - P_U = P_U, P_V = P_V, missing = TRUE, - perc_nonmissing = perc_nonmissing, - step_size = step_size, control = control)) - if ('try-error' %in% class(fit)){ - return(NA) - } else { - norm_temp <- norm_temp + sum((fit$learner_estimate[myset] - Y_target[myset])^2) - } - } - return(norm_temp) -} diff --git a/R/learner-package.R b/R/learner-package.R new file mode 100644 index 0000000..c9c0ad8 --- /dev/null +++ b/R/learner-package.R @@ -0,0 +1,6 @@ +## usethis namespace: start +#' @import RcppEigen +#' @importFrom Rcpp sourceCpp +#' @useDynLib learner, .registration = TRUE +## usethis namespace: end +NULL diff --git a/R/learner.R b/R/learner.R index 0df2bc7..4346c25 100644 --- a/R/learner.R +++ b/R/learner.R @@ -1,198 +1,63 @@ -#' Latent space-based transfer learning +#' Cross-validation for LEARNER #' -#' This function applies the LatEnt spAce-based tRaNsfer lEaRning (LEARNER) method (McGrath et al. 2024) to leverage data from a source population to improve -#' estimation of a low rank matrix in an underrepresented target population. +#' This function performs k-fold cross-validation to select the nuisance parameters \eqn{(\lambda_1, \lambda_2)} for \code{\link{learner}}. #' -#' @param Y_target matrix containing the target population data -#' @param Y_source matrix containing the source population data -#' @param r (optional) integer specifying the rank of the knowledge graphs. By default, ScreeNOT (Donoho et al. 2023) is applied to the source population knowledge graph to select the rank. -#' @param lambda_1 numeric scalar specifying the value of \eqn{\lambda_1} (see Details) -#' @param lambda_2 numeric scalar specifying the value of \eqn{\lambda_2} (see Details) -#' @param step_size numeric scalar specifying the step size for the Newton steps in the numerical optimization algorithm -#' @param control a list of parameters for controlling the stopping criteria for the numerical optimization algorithm. The list may include the following components: -#' \tabular{ll}{ -#' \code{max_iter} \tab integer specifying the maximum number of iterations \cr -#' \code{threshold} \tab numeric scalar specifying a convergence threshold. The algorithm converges when \eqn{|\epsilon_t - \epsilon_{t-1}| < }\code{threshold}, where \eqn{\epsilon_t} denotes the value of the objective function at iteration \eqn{t}. \cr -#' \code{max_value} \tab numeric scalar used to specify the maximum value of the objective function allowed before terminating the algorithm. Specifically, the algorithm will terminate if the value of the objective function exceeds \code{max_value}\eqn{\times \epsilon_0}, where \eqn{\epsilon_0} denotes the value of the objective function at the initial point. This is used to prevent unnecessary computation time after the optimization algorithm diverges. \cr} +#' @param Y_target matrix containing the target population data, as in \code{\link{learner}} +#' @param Y_source matrix containing the source population data, as in \code{\link{learner}} +#' @param r (optional) integer specifying the rank of the knowledge graphs, as in \code{\link{learner}} +#' @param lambda_1_all vector of numerics specifying the candidate values of \eqn{\lambda_1} (see Details) +#' @param lambda_2_all vector of numerics specifying the candidate values of \eqn{\lambda_2} (see Details) +#' @param step_size numeric scalar specifying the step size for the Newton steps in the numerical optimization algorithm, as in \code{\link{learner}} +#' @param n_folds an integer specify the number of cross-validation folds. The default is \code{4}. +#' @param n_cores an integer specifying the number of CPU cores in parallelization. Parallelization is performed across the different candidate \eqn{(\lambda_1, \lambda_2)} pairs. The default is \code{1}, i.e., no parallelization. +#' @param threshold Convergence threshold. +#' @param max_iter Maximum number of iterations. #' #' @return A list with the following elements: -#' \item{learner_estimate}{matrix containing the LEARNER estimate of the target population knowledge graph} -#' \item{objective_values}{numeric vector containing the values of the objective function at each iteration} -#' \item{convergence_criterion}{integer specifying the criterion that was satisfied for terminating the numerical optimization algorithm. A value of 1 indicates the convergence threshold was satisfied; A value of 2 indicates that the maximum number of iterations was satisfied; A value of 3 indicates that the maximum value of the objective function was satisfied.} +#' \item{lambda_1_min}{value of \eqn{\lambda_1} with the smallest MSE} +#' \item{lambda_2_min}{value of \eqn{\lambda_2} with the smallest MSE} +#' \item{mse_all}{matrix containing MSE value for each \eqn{(\lambda_1, \lambda_2)} pair. The rows correspond to the \eqn{\lambda_1} values, and the columns correspond to the \eqn{\lambda_2} values.} #' \item{r}{rank value used.} #' #' @details -#' -#' \strong{Data and notation:} -#' -#' The data consists of a matrix in the target population \eqn{Y_0 \in \mathbb{R}^{p \times q}} and the source population \eqn{Y_1 \in \mathbb{R}^{p \times q}}. -#' Let \eqn{\hat{U}_{k} \hat{\Lambda}_{k} \hat{V}_{k}^{\top}} denote the truncated singular value decomposition (SVD) of \eqn{Y_k}, \eqn{k = 0, 1}. -#' -#' For \eqn{k = 0, 1}, one can view \eqn{Y_k} as a noisy version of \eqn{\Theta_k}, referred to as the knowledge graph. The target of inference is the target population knowledge graph, \eqn{\Theta_0}. -#' -#' \strong{Estimation:} -#' -#' This method estimates \eqn{\Theta_0} by \eqn{\tilde{U}\tilde{V}^{\top}}, where \eqn{(\tilde{U}, \tilde{V})} is the solution to the following optimization problem -#' \deqn{\mathrm{arg\,min}_{U \in \mathbb{R}^{p \times r}, V \in \mathbb{R}^{q \times r}} \big\{ \| U V^{\top} - Y_0 \|_F^2 + \lambda_1\| \mathcal{P}_{\perp}(\hat{U}_{1})U \|_F^2 + \lambda_1\| \mathcal{P}_{\perp}(\hat{V}_{1})V \|_F^2 + \lambda_2 \| U^{\top} U - V^{\top} V \|_F^2 \big\}} -#' where \eqn{\mathcal{P}_{\perp}(\hat{U}_{1}) = I - \hat{U}_{1}^{\top}\hat{U}_{1}} and \eqn{\mathcal{P}_{\perp}(\hat{V}_{1}) = I - \hat{V}_{1}^{\top}\hat{V}_{1}}. -#' -#' This function uses an alternating minimization strategy to solve the optimization problem. That is, this approach updates \eqn{U} by minimizing the objective function (via a gradient descent step) treating \eqn{V} as fixed. Then, \eqn{V} is updated treating \eqn{U} as fixed. These updates of \eqn{U} and \eqn{V} are repeated until convergence. +#' Given sets of candidate values of \eqn{\lambda_1} and \eqn{\lambda_2}, this function performs k-fold cross-validation to select the pair \eqn{(\lambda_1, \lambda_2)} with the smallest held out error. This function randomly partitions the entries of \code{Y_target} into \eqn{k} (approximately) equally sized subsamples. The training data sets are obtained by removing one of the \eqn{k} subsamples and the corresponding test data sets are based on the held out subsamples. The \code{\link{learner}} function is applied to each training data set. The held out error is computed by the mean squared error comparing the entries in the test data sets with those imputed from the LEARNER estimates. See McGrath et al. (2024) for further details. #' #' @references #' McGrath, S., Zhu, C,. Guo, M. and Duan, R. (2024). \emph{LEARNER: A transfer learning method for low-rank matrix estimation}. arXiv preprint arXiv:2412.20605. #' -#' Donoho, D., Gavish, M. and Romanov, E. (2023). \emph{ScreeNOT: Exact MSE-optimal singular value thresholding in correlated noise}. The Annals of Statistics, 51(1), pp.122-148. -#' #' @examples -#' res <- learner(Y_source = dat_highsim$Y_source, -#' Y_target = dat_highsim$Y_target, -#' lambda_1 = 1, lambda_2 = 1, -#' step_size = 0.003) +#' res <- cv.learner(Y_source = dat_highsim$Y_source, +#' Y_target = dat_highsim$Y_target, +#' lambda_1_all = c(1, 10, 100), +#' lambda_2_all = c(1, 10, 100), +#' step_size = 0.003) +#' +#' #' #' @export -learner <- function(Y_source, Y_target, r, lambda_1, lambda_2, step_size, - control = list()){ - ## Error catching - if (!identical(dim(Y_source), dim(Y_target))){ - stop('Y_source and Y_target must have the same dimensions') - } - if (any(is.na(Y_source))){ - stop('Y_source cannot have NA values.') - } - if (length(lambda_1) > 1){ - stop('lambda_1 must be a scalar. See the function cv.learner for selecting a suitable lambda_1 value.') - } - if (length(lambda_2) > 1){ - stop('lambda_2 must be a scalar. See the function cv.learner for selecting a suitable lambda_2 value.') - } - - ## Setting the rank and other parameters - missing <- any(is.na(Y_target)) - p <- nrow(Y_source) - q <- ncol(Y_source) - perc_nonmissing <- sum(!is.na(Y_target)) / (p * q) - if (missing(r)){ - max_rank <- min(p, q) / 3 - r <- max(ScreeNOT::adaptiveHardThresholding(Y = Y_source, k = max_rank)$r, 1) - } - if (!('max_iter' %in% names(control))){ - control$max_iter <- 100 - } - if (!('threshold' %in% names(control))){ - control$threshold <- 0.001 - } - if (!('max_value' %in% names(control))){ - control$max_value <- 10 - } - - svd_source <- svd(Y_source, nu = r, nv = r) - P_U <- diag(rep(1, p)) - svd_source$u %*% t(svd_source$u) - P_V <- diag(rep(1, q)) - svd_source$v %*% t(svd_source$v) - U_init <- svd_source$u %*% diag(sqrt(svd_source$d[1:r]), nrow = r, ncol = r) - V_init <- svd_source$v %*% diag(sqrt(svd_source$d[1:r]), nrow = r, ncol = r) - - out <- learner_helper(control = control, U_init = U_init, V_init = V_init, lambda1 = lambda_1, - lambda2 = lambda_2, theta_hat = Y_target, P_U = P_U, P_V = P_V, missing = missing, - perc_nonmissing = perc_nonmissing, step_size = step_size) - - return(out) +cv.learner <- function(Y_source, Y_target, lambda_1_all, lambda_2_all, step_size, + n_folds = 4, max_iter = 100, threshold = 1e-3, n_cores = 1, r = 1) { + result <- cv_learner_cpp(Y_source, Y_target, lambda_1_all, lambda_2_all, step_size, + n_folds, max_iter, threshold, n_cores, r) + return(result) } -learner_helper <- function(control, U_init, V_init, lambda1, lambda2, theta_hat, P_U, P_V, missing, - perc_nonmissing, step_size){ - U <- U_init - V <- V_init - max_iter <- control$max_iter - threshold <- control$threshold - max_value <- control$max_value - - obj_func_init <- f(U = U, V = V, lambda1 = lambda1, lambda2 = lambda2, theta_hat = theta_hat, P_U = P_U, P_V = P_V, missing = missing, - perc_nonmissing = perc_nonmissing) - obj_func_best <- obj_func_init; U_best <- U; V_best <- V - obj_func <- rep(NA, times = max_iter) - U_norm <- norm(U, type = 'F') - V_norm <- norm(V, type = 'F') - - for (iter in 1:max_iter){ - f_delta_U <- f_prime_u(U = U, V = V, lambda1 = lambda1, lambda2 = lambda2, theta_hat = theta_hat, P_U = P_U, missing = missing, perc_nonmissing = perc_nonmissing) - f_delta_V <- f_prime_v(U = U, V = V, lambda1 = lambda1, lambda2 = lambda2, theta_hat = theta_hat, P_V = P_V, missing = missing, perc_nonmissing = perc_nonmissing) - - f_delta_U_norm <- norm(f_delta_U, type = 'F') - f_delta_V_norm <- norm(f_delta_V, type = 'F') - - U <- U - (step_size * U_norm / f_delta_U_norm) * f_delta_U - V <- V - (step_size * V_norm / f_delta_V_norm) * f_delta_V - - obj_func[iter] <- f(U = U, V = V, lambda1 = lambda1, lambda2 = lambda2, theta_hat = theta_hat, P_U = P_U, P_V = P_V, missing = missing, perc_nonmissing = perc_nonmissing) - if (obj_func[iter] < obj_func_best){ - obj_func_best <- obj_func[iter]; U_best <- U; V_best <- V - } - if (iter > 1){ - if (abs(obj_func[iter] - obj_func[iter - 1]) < threshold){ - convergence_criterion <- 1 - break() - } - if (obj_func[iter] > max_value * obj_func_init){ - convergence_criterion <- 3 - break() - } - } - U_norm <- norm(U, type = 'F') - V_norm <- norm(V, type = 'F') - - if (iter == max_iter){ - convergence_criterion <- 2 - } - } - theta_hat <- U_best %*% t(V_best) - return(list(learner_estimate = theta_hat, objective_values = obj_func, - convergence_criterion = convergence_criterion, - r = ncol(U))) -} - - -norm_2_sq <- function(x){ - return(sum(x^2)) -} - -f <- function(U, V, lambda1, lambda2, theta_hat, P_U, P_V, missing, perc_nonmissing){ - if (missing){ - t1 <- sum((U %*% t(V) - theta_hat)^2, na.rm = TRUE) / perc_nonmissing - } else { - t1 <- norm(U %*% t(V) - theta_hat, type = 'F')^2 / perc_nonmissing - } - t2 <- lambda1 * norm(P_U %*% U, type = 'F')^2 - t3 <- lambda1 * norm(P_V %*% V, type = 'F')^2 - t4 <- lambda2 * norm(t(U) %*% U - t(V) %*% V, type = 'F')^2 - return(t1 + t2 + t3 + t4) -} - -f_prime_u <- function(U, V, lambda1, lambda2, theta_hat, P_U, missing, perc_nonmissing){ - U_tilde <- t(U) %*% U - V_tilde <- t(V) %*% V - - if (missing){ - temp <- U %*% t(V) - missing_entries <- which(is.na(theta_hat), arr.ind = T) - theta_hat[missing_entries] <- temp[missing_entries] - } - t1 <- (2 / perc_nonmissing) * (U %*% V_tilde - theta_hat %*% V) - t2 <- lambda1 * 2 * P_U %*% U - t4 <- lambda2 * 4 * U %*% (U_tilde - V_tilde) - return(t1 + t2 + t4) +#' Learner Optimization +#' +#' @param Y_source Matrix containing the source population data. +#' @param Y_target Matrix containing the target population data. +#' @param lambda_1 Regularization parameter lambda1. +#' @param lambda_2 Regularization parameter lambda2. +#' @param step_size Step size for optimization. +#' @param max_iter Maximum number of iterations. +#' @param threshold Convergence threshold. +#' @param r Rank used for the truncated SVD. +#' @return A list containing the learner estimate, objective values, convergence criterion, and rank. +#' @export +learner <- function(Y_source, Y_target, lambda_1, lambda_2, step_size, + max_iter = 100, threshold = 1e-3, r = 1) { + result <- learner_cpp(Y_source, Y_target, r, lambda_1, lambda_2, step_size, max_iter, threshold) + return(result) } -f_prime_v <- function(U, V, lambda1, lambda2, theta_hat, P_V, missing, perc_nonmissing){ - U_tilde <- t(U) %*% U - V_tilde <- t(V) %*% V - if (missing){ - temp <- U %*% t(V) - missing_entries <- which(is.na(theta_hat), arr.ind = T) - theta_hat[missing_entries] <- temp[missing_entries] - } - t1 <- (2 / perc_nonmissing) * (V %*% U_tilde - t(theta_hat) %*% U) - t3 <- lambda1 * 2 * P_V %*% V - t4 <- lambda2 * 4 * V %*% (V_tilde - U_tilde) - return(t1 + t3 + t4) -} diff --git a/man/cv.learner.Rd b/man/cv.learner.Rd index b0950c4..94d7956 100644 --- a/man/cv.learner.Rd +++ b/man/cv.learner.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/cv_learner.R +% Please edit documentation in R/learner.R \name{cv.learner} \alias{cv.learner} \title{Cross-validation for LEARNER} @@ -7,13 +7,14 @@ cv.learner( Y_source, Y_target, - r, lambda_1_all, lambda_2_all, step_size, n_folds = 4, + max_iter = 100, + threshold = 0.001, n_cores = 1, - control = list() + r = 1 ) } \arguments{ @@ -21,8 +22,6 @@ cv.learner( \item{Y_target}{matrix containing the target population data, as in \code{\link{learner}}} -\item{r}{(optional) integer specifying the rank of the knowledge graphs, as in \code{\link{learner}}} - \item{lambda_1_all}{vector of numerics specifying the candidate values of \eqn{\lambda_1} (see Details)} \item{lambda_2_all}{vector of numerics specifying the candidate values of \eqn{\lambda_2} (see Details)} @@ -31,9 +30,13 @@ cv.learner( \item{n_folds}{an integer specify the number of cross-validation folds. The default is \code{4}.} +\item{max_iter}{Maximum number of iterations.} + +\item{threshold}{Convergence threshold.} + \item{n_cores}{an integer specifying the number of CPU cores in parallelization. Parallelization is performed across the different candidate \eqn{(\lambda_1, \lambda_2)} pairs. The default is \code{1}, i.e., no parallelization.} -\item{control}{a list of parameters for controlling the stopping criteria for the numerical optimization algorithm, as in \code{\link{learner}}.} +\item{r}{(optional) integer specifying the rank of the knowledge graphs, as in \code{\link{learner}}} } \value{ A list with the following elements: @@ -56,6 +59,7 @@ res <- cv.learner(Y_source = dat_highsim$Y_source, step_size = 0.003) + } \references{ McGrath, S., Zhu, C,. Guo, M. and Duan, R. (2024). \emph{LEARNER: A transfer learning method for low-rank matrix estimation}. arXiv preprint arXiv:2412.20605. diff --git a/man/learner.Rd b/man/learner.Rd index e1bdce4..e05e175 100644 --- a/man/learner.Rd +++ b/man/learner.Rd @@ -2,65 +2,39 @@ % Please edit documentation in R/learner.R \name{learner} \alias{learner} -\title{Latent space-based transfer learning} +\title{Learner Optimization} \usage{ -learner(Y_source, Y_target, r, lambda_1, lambda_2, step_size, control = list()) +learner( + Y_source, + Y_target, + lambda_1, + lambda_2, + step_size, + max_iter = 100, + threshold = 0.001, + r = 1 +) } \arguments{ -\item{Y_source}{matrix containing the source population data} +\item{Y_source}{Matrix containing the source population data.} -\item{Y_target}{matrix containing the target population data} +\item{Y_target}{Matrix containing the target population data.} -\item{r}{(optional) integer specifying the rank of the knowledge graphs. By default, ScreeNOT (Donoho et al. 2023) is applied to the source population knowledge graph to select the rank.} +\item{lambda_1}{Regularization parameter lambda1.} -\item{lambda_1}{numeric scalar specifying the value of \eqn{\lambda_1} (see Details)} +\item{lambda_2}{Regularization parameter lambda2.} -\item{lambda_2}{numeric scalar specifying the value of \eqn{\lambda_2} (see Details)} +\item{step_size}{Step size for optimization.} -\item{step_size}{numeric scalar specifying the step size for the Newton steps in the numerical optimization algorithm} +\item{max_iter}{Maximum number of iterations.} -\item{control}{a list of parameters for controlling the stopping criteria for the numerical optimization algorithm. The list may include the following components: -\tabular{ll}{ -\code{max_iter} \tab integer specifying the maximum number of iterations \cr -\code{threshold} \tab numeric scalar specifying a convergence threshold. The algorithm converges when \eqn{|\epsilon_t - \epsilon_{t-1}| < }\code{threshold}, where \eqn{\epsilon_t} denotes the value of the objective function at iteration \eqn{t}. \cr -\code{max_value} \tab numeric scalar used to specify the maximum value of the objective function allowed before terminating the algorithm. Specifically, the algorithm will terminate if the value of the objective function exceeds \code{max_value}\eqn{\times \epsilon_0}, where \eqn{\epsilon_0} denotes the value of the objective function at the initial point. This is used to prevent unnecessary computation time after the optimization algorithm diverges. \cr}} +\item{threshold}{Convergence threshold.} + +\item{r}{Rank used for the truncated SVD.} } \value{ -A list with the following elements: -\item{learner_estimate}{matrix containing the LEARNER estimate of the target population knowledge graph} -\item{objective_values}{numeric vector containing the values of the objective function at each iteration} -\item{convergence_criterion}{integer specifying the criterion that was satisfied for terminating the numerical optimization algorithm. A value of 1 indicates the convergence threshold was satisfied; A value of 2 indicates that the maximum number of iterations was satisfied; A value of 3 indicates that the maximum value of the objective function was satisfied.} -\item{r}{rank value used.} +A list containing the learner estimate, objective values, convergence criterion, and rank. } \description{ -This function applies the LatEnt spAce-based tRaNsfer lEaRning (LEARNER) method (McGrath et al. 2024) to leverage data from a source population to improve -estimation of a low rank matrix in an underrepresented target population. -} -\details{ -\strong{Data and notation:} - -The data consists of a matrix in the target population \eqn{Y_0 \in \mathbb{R}^{p \times q}} and the source population \eqn{Y_1 \in \mathbb{R}^{p \times q}}. -Let \eqn{\hat{U}_{k} \hat{\Lambda}_{k} \hat{V}_{k}^{\top}} denote the truncated singular value decomposition (SVD) of \eqn{Y_k}, \eqn{k = 0, 1}. - -For \eqn{k = 0, 1}, one can view \eqn{Y_k} as a noisy version of \eqn{\Theta_k}, referred to as the knowledge graph. The target of inference is the target population knowledge graph, \eqn{\Theta_0}. - -\strong{Estimation:} - -This method estimates \eqn{\Theta_0} by \eqn{\tilde{U}\tilde{V}^{\top}}, where \eqn{(\tilde{U}, \tilde{V})} is the solution to the following optimization problem -\deqn{\mathrm{arg\,min}_{U \in \mathbb{R}^{p \times r}, V \in \mathbb{R}^{q \times r}} \big\{ \| U V^{\top} - Y_0 \|_F^2 + \lambda_1\| \mathcal{P}_{\perp}(\hat{U}_{1})U \|_F^2 + \lambda_1\| \mathcal{P}_{\perp}(\hat{V}_{1})V \|_F^2 + \lambda_2 \| U^{\top} U - V^{\top} V \|_F^2 \big\}} -where \eqn{\mathcal{P}_{\perp}(\hat{U}_{1}) = I - \hat{U}_{1}^{\top}\hat{U}_{1}} and \eqn{\mathcal{P}_{\perp}(\hat{V}_{1}) = I - \hat{V}_{1}^{\top}\hat{V}_{1}}. - -This function uses an alternating minimization strategy to solve the optimization problem. That is, this approach updates \eqn{U} by minimizing the objective function (via a gradient descent step) treating \eqn{V} as fixed. Then, \eqn{V} is updated treating \eqn{U} as fixed. These updates of \eqn{U} and \eqn{V} are repeated until convergence. -} -\examples{ -res <- learner(Y_source = dat_highsim$Y_source, - Y_target = dat_highsim$Y_target, - lambda_1 = 1, lambda_2 = 1, - step_size = 0.003) - -} -\references{ -McGrath, S., Zhu, C,. Guo, M. and Duan, R. (2024). \emph{LEARNER: A transfer learning method for low-rank matrix estimation}. arXiv preprint arXiv:2412.20605. - -Donoho, D., Gavish, M. and Romanov, E. (2023). \emph{ScreeNOT: Exact MSE-optimal singular value thresholding in correlated noise}. The Annals of Statistics, 51(1), pp.122-148. +Learner Optimization } diff --git a/src/.gitignore b/src/.gitignore new file mode 100644 index 0000000..22034c4 --- /dev/null +++ b/src/.gitignore @@ -0,0 +1,3 @@ +*.o +*.so +*.dll diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..a90b0d4 --- /dev/null +++ b/src/Makevars @@ -0,0 +1,3 @@ +PKG_CXXFLAGS = $(SHLIB_CXXFLAGS) $(SHLIB_OPENMP_CXXFLAGS) +PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) + diff --git a/src/Makevars.win b/src/Makevars.win new file mode 100644 index 0000000..a90b0d4 --- /dev/null +++ b/src/Makevars.win @@ -0,0 +1,3 @@ +PKG_CXXFLAGS = $(SHLIB_CXXFLAGS) $(SHLIB_OPENMP_CXXFLAGS) +PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) + diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp new file mode 100644 index 0000000..8a575a9 --- /dev/null +++ b/src/RcppExports.cpp @@ -0,0 +1,62 @@ +// Generated by using Rcpp::compileAttributes() -> do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#include +#include + +using namespace Rcpp; + +#ifdef RCPP_USE_GLOBAL_ROSTREAM +Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); +Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); +#endif + +// learner_cpp +List learner_cpp(const Eigen::MatrixXd& Y_source, const Eigen::MatrixXd& Y_target, int r, double lambda1, double lambda2, double step_size, int max_iter, double threshold); +RcppExport SEXP _learner_learner_cpp(SEXP Y_sourceSEXP, SEXP Y_targetSEXP, SEXP rSEXP, SEXP lambda1SEXP, SEXP lambda2SEXP, SEXP step_sizeSEXP, SEXP max_iterSEXP, SEXP thresholdSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const Eigen::MatrixXd& >::type Y_source(Y_sourceSEXP); + Rcpp::traits::input_parameter< const Eigen::MatrixXd& >::type Y_target(Y_targetSEXP); + Rcpp::traits::input_parameter< int >::type r(rSEXP); + Rcpp::traits::input_parameter< double >::type lambda1(lambda1SEXP); + Rcpp::traits::input_parameter< double >::type lambda2(lambda2SEXP); + Rcpp::traits::input_parameter< double >::type step_size(step_sizeSEXP); + Rcpp::traits::input_parameter< int >::type max_iter(max_iterSEXP); + Rcpp::traits::input_parameter< double >::type threshold(thresholdSEXP); + rcpp_result_gen = Rcpp::wrap(learner_cpp(Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold)); + return rcpp_result_gen; +END_RCPP +} +// cv_learner_cpp +List cv_learner_cpp(const Eigen::MatrixXd& Y_source, const Eigen::MatrixXd& Y_target, const std::vector& lambda1_all, const std::vector& lambda2_all, double step_size, int n_folds, int max_iter, double threshold, int n_cores, int r); +RcppExport SEXP _learner_cv_learner_cpp(SEXP Y_sourceSEXP, SEXP Y_targetSEXP, SEXP lambda1_allSEXP, SEXP lambda2_allSEXP, SEXP step_sizeSEXP, SEXP n_foldsSEXP, SEXP max_iterSEXP, SEXP thresholdSEXP, SEXP n_coresSEXP, SEXP rSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const Eigen::MatrixXd& >::type Y_source(Y_sourceSEXP); + Rcpp::traits::input_parameter< const Eigen::MatrixXd& >::type Y_target(Y_targetSEXP); + Rcpp::traits::input_parameter< const std::vector& >::type lambda1_all(lambda1_allSEXP); + Rcpp::traits::input_parameter< const std::vector& >::type lambda2_all(lambda2_allSEXP); + Rcpp::traits::input_parameter< double >::type step_size(step_sizeSEXP); + Rcpp::traits::input_parameter< int >::type n_folds(n_foldsSEXP); + Rcpp::traits::input_parameter< int >::type max_iter(max_iterSEXP); + Rcpp::traits::input_parameter< double >::type threshold(thresholdSEXP); + Rcpp::traits::input_parameter< int >::type n_cores(n_coresSEXP); + Rcpp::traits::input_parameter< int >::type r(rSEXP); + rcpp_result_gen = Rcpp::wrap(cv_learner_cpp(Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r)); + return rcpp_result_gen; +END_RCPP +} + +static const R_CallMethodDef CallEntries[] = { + {"_learner_learner_cpp", (DL_FUNC) &_learner_learner_cpp, 8}, + {"_learner_cv_learner_cpp", (DL_FUNC) &_learner_cv_learner_cpp, 10}, + {NULL, NULL, 0} +}; + +RcppExport void R_init_learner(DllInfo *dll) { + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} diff --git a/src/kernels.cpp b/src/kernels.cpp new file mode 100644 index 0000000..53c91eb --- /dev/null +++ b/src/kernels.cpp @@ -0,0 +1,191 @@ +#include +#include +#include +#include +#include +#include + +#ifdef _OPENMP +#include +#endif + +using namespace Rcpp; +using namespace Eigen; + +//-------------------------------------------------------------- +// Internal worker function: computes the learner estimate - thread safe for OMP +Eigen::MatrixXd learner_worker(const Eigen::MatrixXd &Y_source, + const Eigen::MatrixXd &Y_target, + int r, double lambda1, double lambda2, + double step_size, int max_iter, double threshold) { + int p = Y_source.rows(); + int q = Y_source.cols(); + + Eigen::JacobiSVD svd(Y_source, Eigen::ComputeThinU | Eigen::ComputeThinV); + Eigen::MatrixXd U_full = svd.matrixU(); + Eigen::MatrixXd V_full = svd.matrixV(); + int r_use = std::min(r, static_cast(U_full.cols())); + Eigen::MatrixXd U_trunc = U_full.leftCols(r_use); + Eigen::MatrixXd V_trunc = V_full.leftCols(r_use); + Eigen::VectorXd singular_vals = svd.singularValues().head(r_use); + + // Initialize as in R: U_init = U_trunc * sqrt(D), V_init = V_trunc * sqrt(D) + Eigen::MatrixXd sqrtD = singular_vals.array().sqrt().matrix().asDiagonal(); + Eigen::MatrixXd U = U_trunc * sqrtD; + Eigen::MatrixXd V = V_trunc * sqrtD; + + // Projection matrices (using truncated basis) + Eigen::MatrixXd P_U = Eigen::MatrixXd::Identity(p, p) - U_trunc * U_trunc.transpose(); + Eigen::MatrixXd P_V = Eigen::MatrixXd::Identity(q, q) - V_trunc * V_trunc.transpose(); + + // For now, we assume no missing data (perc_nonmissing = 1). + // double perc_nonmissing = (Y_target.array().isNaN().count() > 0) ? ( (p*q - num_missing) / double(p*q)) : 1.0; + double perc_nonmissing = 1.0; + bool missing = Y_target.hasNaN(); + + double obj_init = 0.0; + { + Eigen::MatrixXd diff = U * V.transpose() - Y_target; + if (missing) { + diff = diff.array().isNaN().select(Eigen::MatrixXd::Zero(diff.rows(), diff.cols()), diff); + } + obj_init = diff.squaredNorm() / perc_nonmissing + + lambda1 * (P_U * U).squaredNorm() + + lambda1 * (P_V * V).squaredNorm() + + lambda2 * (U.transpose() * U - V.transpose() * V).squaredNorm(); + } + + double obj_best = obj_init; + Eigen::MatrixXd U_best = U; + Eigen::MatrixXd V_best = V; + double U_norm = U.norm(); + double V_norm = V.norm(); + + for (int iter = 0; iter < max_iter; ++iter) { + Eigen::MatrixXd U_tilde = U.transpose() * U; + Eigen::MatrixXd V_tilde = V.transpose() * V; + + Eigen::MatrixXd adjusted_theta = Y_target; + if (missing) { + Eigen::MatrixXd temp = U * V.transpose(); + adjusted_theta = adjusted_theta.array().isNaN().select(temp, adjusted_theta); + } + + Eigen::MatrixXd grad_U = (2.0 / perc_nonmissing) * (U * V_tilde - adjusted_theta * V) + + lambda1 * 2 * P_U * U + + lambda2 * 4 * U * (U_tilde - V_tilde); + Eigen::MatrixXd grad_V = (2.0 / perc_nonmissing) * (V * U_tilde - adjusted_theta.transpose() * U) + + lambda1 * 2 * P_V * V + + lambda2 * 4 * V * (V_tilde - U_tilde); + + double grad_U_norm = grad_U.norm(); + double grad_V_norm = grad_V.norm(); + + U = U - (step_size * U_norm / (grad_U_norm + 1e-12)) * grad_U; + V = V - (step_size * V_norm / (grad_V_norm + 1e-12)) * grad_V; + U_norm = U.norm(); + V_norm = V.norm(); + + double obj = 0.0; + { + Eigen::MatrixXd diff = U * V.transpose() - Y_target; + if (missing) { + diff = diff.array().isNaN().select(Eigen::MatrixXd::Zero(diff.rows(), diff.cols()), diff); + } + obj = diff.squaredNorm() / perc_nonmissing + + lambda1 * (P_U * U).squaredNorm() + + lambda1 * (P_V * V).squaredNorm() + + lambda2 * (U.transpose() * U - V.transpose() * V).squaredNorm(); + } + + if (obj < obj_best) { + obj_best = obj; + U_best = U; + V_best = V; + } + + // Simple convergence check (you might refine this). + if (iter > 0 && std::abs(obj - obj_init) < threshold) { + break; + } + obj_init = obj; + } + + return U_best * V_best.transpose(); +} + +//-------------------------------------------------------------- +// Exported learner function. +// [[Rcpp::export]] +List learner_cpp(const Eigen::MatrixXd &Y_source, const Eigen::MatrixXd &Y_target, + int r, double lambda1, double lambda2, double step_size, int max_iter, double threshold) { + Eigen::MatrixXd learner_estimate = learner_worker(Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold); + return List::create( + Named("learner_estimate") = learner_estimate + ); +} + +//-------------------------------------------------------------- +// [[Rcpp::export]] +List cv_learner_cpp(const Eigen::MatrixXd &Y_source, const Eigen::MatrixXd &Y_target, + const std::vector &lambda1_all, const std::vector &lambda2_all, + double step_size, int n_folds, int max_iter, double threshold, + int n_cores, int r) { + int p = Y_source.rows(); + int q = Y_source.cols(); + int n_lambda1 = lambda1_all.size(); + int n_lambda2 = lambda2_all.size(); + + std::vector indices(p * q); + std::iota(indices.begin(), indices.end(), 0); + std::shuffle(indices.begin(), indices.end(), std::mt19937{std::random_device{}()}); + + std::vector> index_set(n_folds); + for (int fold = 0; fold < n_folds; ++fold) { + int start = fold * (indices.size() / n_folds); + int end = (fold + 1) * (indices.size() / n_folds); + index_set[fold] = std::vector(indices.begin() + start, indices.begin() + end); + } + + Eigen::MatrixXd mse_all = Eigen::MatrixXd::Zero(n_lambda1, n_lambda2); + +#pragma omp parallel for collapse(2) schedule(dynamic) num_threads(n_cores) + for (int i = 0; i < n_lambda1; ++i) { + for (int j = 0; j < n_lambda2; ++j) { + double mse = 0.0; + + for (int fold = 0; fold < n_folds; ++fold) { + Eigen::MatrixXd Y_train = Y_target; + for (int idx : index_set[fold]) { + int row = idx / q; + int col = idx % q; + Y_train(row, col) = std::numeric_limits::quiet_NaN(); + } + Eigen::MatrixXd learner_estimate = learner_worker(Y_source, Y_train, r, + lambda1_all[i], lambda2_all[j], + step_size, max_iter, threshold); + double fold_mse = 0.0; + for (int idx : index_set[fold]) { + int row = idx / q; + int col = idx % q; + double diff = learner_estimate(row, col) - Y_target(row, col); + fold_mse += diff * diff; + } + mse += fold_mse; + } + mse_all(i, j) = mse; + } + } + + Eigen::Index minRow, minCol; + mse_all.minCoeff(&minRow, &minCol); + + return List::create( + Named("lambda1_min") = lambda1_all[minRow], + Named("lambda2_min") = lambda2_all[minCol], + Named("mse_all") = mse_all, + Named("r") = r + ); +} + + diff --git a/src/learner-package.cpp b/src/learner-package.cpp new file mode 100644 index 0000000..540847e --- /dev/null +++ b/src/learner-package.cpp @@ -0,0 +1,2 @@ +#include +using namespace Rcpp; From 88d367b0de854984b70cf824538cb34f551084fb Mon Sep 17 00:00:00 2001 From: Ryan O'Dea <70209371+ryan-odea@users.noreply.github.com> Date: Wed, 19 Feb 2025 01:36:36 -0500 Subject: [PATCH 02/16] checking for OMP installation --- R/RcppExports.R | 4 ++++ R/globals.R | 11 +++++++++++ R/learner.R | 2 +- man/cv.learner.Rd | 2 +- src/RcppExports.cpp | 11 +++++++++++ src/kernels.cpp | 9 +++++++++ 6 files changed, 37 insertions(+), 2 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 0c3b1d7..3bccd12 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -9,3 +9,7 @@ cv_learner_cpp <- function(Y_source, Y_target, lambda1_all, lambda2_all, step_si .Call(`_learner_cv_learner_cpp`, Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r) } +omp_max_threads <- function() { + .Call(`_learner_omp_max_threads`) +} + diff --git a/R/globals.R b/R/globals.R index 5af0ca7..2a56d31 100644 --- a/R/globals.R +++ b/R/globals.R @@ -2,3 +2,14 @@ if(getRversion() >= "2.15.1"){ # To remove 'no visible binding for global variable ...' notes utils::globalVariables(c('lambda_1_ind', 'lambda_2_ind')) } + +.onAttach <- function(libname, pkgname) { + max_threads <- tryCatch(omp_max_threads(), error = function(e) 1) + if (max_threads < 2) { + packageStartupMessage("*******\n + This installation of learner has not detected OpenMP support\n + It will still work but will not support multithreading via the `ncores` argument + If you plan to use multithreading, please ensure you are using R>=3.4.0 and have OpenMP installed\n + *******") + } +} diff --git a/R/learner.R b/R/learner.R index 4346c25..82aad19 100644 --- a/R/learner.R +++ b/R/learner.R @@ -9,7 +9,7 @@ #' @param lambda_2_all vector of numerics specifying the candidate values of \eqn{\lambda_2} (see Details) #' @param step_size numeric scalar specifying the step size for the Newton steps in the numerical optimization algorithm, as in \code{\link{learner}} #' @param n_folds an integer specify the number of cross-validation folds. The default is \code{4}. -#' @param n_cores an integer specifying the number of CPU cores in parallelization. Parallelization is performed across the different candidate \eqn{(\lambda_1, \lambda_2)} pairs. The default is \code{1}, i.e., no parallelization. +#' @param n_cores an integer specifying the number of CPU cores in OpenMP parallelization. Parallelization is performed across the different candidate \eqn{(\lambda_1, \lambda_2)} pairs. The default is \code{1}, i.e., no parallelization. #' @param threshold Convergence threshold. #' @param max_iter Maximum number of iterations. #' diff --git a/man/cv.learner.Rd b/man/cv.learner.Rd index 94d7956..ea3f5b3 100644 --- a/man/cv.learner.Rd +++ b/man/cv.learner.Rd @@ -34,7 +34,7 @@ cv.learner( \item{threshold}{Convergence threshold.} -\item{n_cores}{an integer specifying the number of CPU cores in parallelization. Parallelization is performed across the different candidate \eqn{(\lambda_1, \lambda_2)} pairs. The default is \code{1}, i.e., no parallelization.} +\item{n_cores}{an integer specifying the number of CPU cores in OpenMP parallelization. Parallelization is performed across the different candidate \eqn{(\lambda_1, \lambda_2)} pairs. The default is \code{1}, i.e., no parallelization.} \item{r}{(optional) integer specifying the rank of the knowledge graphs, as in \code{\link{learner}}} } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 8a575a9..0eb47c2 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -49,10 +49,21 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// omp_max_threads +int omp_max_threads(); +RcppExport SEXP _learner_omp_max_threads() { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + rcpp_result_gen = Rcpp::wrap(omp_max_threads()); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_learner_learner_cpp", (DL_FUNC) &_learner_learner_cpp, 8}, {"_learner_cv_learner_cpp", (DL_FUNC) &_learner_cv_learner_cpp, 10}, + {"_learner_omp_max_threads", (DL_FUNC) &_learner_omp_max_threads, 0}, {NULL, NULL, 0} }; diff --git a/src/kernels.cpp b/src/kernels.cpp index 53c91eb..d0cb3af 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -188,4 +188,13 @@ List cv_learner_cpp(const Eigen::MatrixXd &Y_source, const Eigen::MatrixXd &Y_ta ); } +// -------------------------------------------------------------- +// [[Rcpp::export]] +int omp_max_threads() { +#ifdef _OPENMP + return omp_get_max_threads(); +#else + return 1; +#endif +} From 7b9a8d786c0145c02d302a7e11a0cdcfaaad3f01 Mon Sep 17 00:00:00 2001 From: Sean McGrath Date: Sat, 22 Feb 2025 15:58:46 -0500 Subject: [PATCH 03/16] Added some error catching back in --- R/learner.R | 38 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 36 insertions(+), 2 deletions(-) diff --git a/R/learner.R b/R/learner.R index 82aad19..dbfc28f 100644 --- a/R/learner.R +++ b/R/learner.R @@ -36,7 +36,21 @@ #' #' @export cv.learner <- function(Y_source, Y_target, lambda_1_all, lambda_2_all, step_size, - n_folds = 4, max_iter = 100, threshold = 1e-3, n_cores = 1, r = 1) { + n_folds = 4, max_iter = 100, threshold = 1e-3, n_cores = 1, r) { + ## Error catching + if (!identical(dim(Y_source), dim(Y_target))){ + stop('Y_source and Y_target must have the same dimensions') + } + if (any(is.na(Y_source))){ + stop('Y_source cannot have NA values.') + } + + if (missing(r)){ + p <- nrow(Y_source); q <- ncol(Y_source) + max_rank <- min(p, q) / 3 + r <- max(ScreeNOT::adaptiveHardThresholding(Y = Y_source, k = max_rank)$r, 1) + } + result <- cv_learner_cpp(Y_source, Y_target, lambda_1_all, lambda_2_all, step_size, n_folds, max_iter, threshold, n_cores, r) return(result) @@ -55,7 +69,27 @@ cv.learner <- function(Y_source, Y_target, lambda_1_all, lambda_2_all, step_size #' @return A list containing the learner estimate, objective values, convergence criterion, and rank. #' @export learner <- function(Y_source, Y_target, lambda_1, lambda_2, step_size, - max_iter = 100, threshold = 1e-3, r = 1) { + max_iter = 100, threshold = 1e-3, r) { + ## Error catching + if (!identical(dim(Y_source), dim(Y_target))){ + stop('Y_source and Y_target must have the same dimensions') + } + if (any(is.na(Y_source))){ + stop('Y_source cannot have NA values.') + } + if (length(lambda_1) > 1){ + stop('lambda_1 must be a scalar. See the function cv.learner for selecting a suitable lambda_1 value.') + } + if (length(lambda_2) > 1){ + stop('lambda_2 must be a scalar. See the function cv.learner for selecting a suitable lambda_2 value.') + } + + if (missing(r)){ + p <- nrow(Y_source); q <- ncol(Y_source) + max_rank <- min(p, q) / 3 + r <- max(ScreeNOT::adaptiveHardThresholding(Y = Y_source, k = max_rank)$r, 1) + } + result <- learner_cpp(Y_source, Y_target, r, lambda_1, lambda_2, step_size, max_iter, threshold) return(result) } From 8d5d7bdabee8580285ddfc28a6d52765ef553b2d Mon Sep 17 00:00:00 2001 From: Sean McGrath Date: Sat, 22 Feb 2025 16:26:01 -0500 Subject: [PATCH 04/16] Added back in the max_value based condition for terminating the optimization algorithm --- R/RcppExports.R | 8 ++++---- R/learner.R | 10 ++++++---- man/cv.learner.Rd | 5 ++++- man/learner.Rd | 5 ++++- src/RcppExports.cpp | 18 ++++++++++-------- src/kernels.cpp | 14 +++++++++----- 6 files changed, 37 insertions(+), 23 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 3bccd12..e267889 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,12 +1,12 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -learner_cpp <- function(Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold) { - .Call(`_learner_learner_cpp`, Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold) +learner_cpp <- function(Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold, max_value) { + .Call(`_learner_learner_cpp`, Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold, max_value) } -cv_learner_cpp <- function(Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r) { - .Call(`_learner_cv_learner_cpp`, Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r) +cv_learner_cpp <- function(Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r, max_value) { + .Call(`_learner_cv_learner_cpp`, Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r, max_value) } omp_max_threads <- function() { diff --git a/R/learner.R b/R/learner.R index dbfc28f..ec11040 100644 --- a/R/learner.R +++ b/R/learner.R @@ -12,6 +12,7 @@ #' @param n_cores an integer specifying the number of CPU cores in OpenMP parallelization. Parallelization is performed across the different candidate \eqn{(\lambda_1, \lambda_2)} pairs. The default is \code{1}, i.e., no parallelization. #' @param threshold Convergence threshold. #' @param max_iter Maximum number of iterations. +#' @param max_value Numeric scalar used to specify the maximum value of the objective function allowed before terminating the algorithm. Specifically, the algorithm will terminate if the value of the objective function exceeds \code{max_value}\eqn{\times \epsilon_0}, where \eqn{\epsilon_0} denotes the value of the objective function at the initial point. This is used to prevent unnecessary computation time after the optimization algorithm diverges. #' #' @return A list with the following elements: #' \item{lambda_1_min}{value of \eqn{\lambda_1} with the smallest MSE} @@ -36,7 +37,7 @@ #' #' @export cv.learner <- function(Y_source, Y_target, lambda_1_all, lambda_2_all, step_size, - n_folds = 4, max_iter = 100, threshold = 1e-3, n_cores = 1, r) { + n_folds = 4, max_iter = 100, threshold = 1e-3, n_cores = 1, r, max_value = 10) { ## Error catching if (!identical(dim(Y_source), dim(Y_target))){ stop('Y_source and Y_target must have the same dimensions') @@ -52,7 +53,7 @@ cv.learner <- function(Y_source, Y_target, lambda_1_all, lambda_2_all, step_size } result <- cv_learner_cpp(Y_source, Y_target, lambda_1_all, lambda_2_all, step_size, - n_folds, max_iter, threshold, n_cores, r) + n_folds, max_iter, threshold, n_cores, r, max_value) return(result) } @@ -66,10 +67,11 @@ cv.learner <- function(Y_source, Y_target, lambda_1_all, lambda_2_all, step_size #' @param max_iter Maximum number of iterations. #' @param threshold Convergence threshold. #' @param r Rank used for the truncated SVD. +#' @param max_value Numeric scalar used to specify the maximum value of the objective function allowed before terminating the algorithm. Specifically, the algorithm will terminate if the value of the objective function exceeds \code{max_value}\eqn{\times \epsilon_0}, where \eqn{\epsilon_0} denotes the value of the objective function at the initial point. This is used to prevent unnecessary computation time after the optimization algorithm diverges. #' @return A list containing the learner estimate, objective values, convergence criterion, and rank. #' @export learner <- function(Y_source, Y_target, lambda_1, lambda_2, step_size, - max_iter = 100, threshold = 1e-3, r) { + max_iter = 100, threshold = 1e-3, r, max_value = 10) { ## Error catching if (!identical(dim(Y_source), dim(Y_target))){ stop('Y_source and Y_target must have the same dimensions') @@ -90,7 +92,7 @@ learner <- function(Y_source, Y_target, lambda_1, lambda_2, step_size, r <- max(ScreeNOT::adaptiveHardThresholding(Y = Y_source, k = max_rank)$r, 1) } - result <- learner_cpp(Y_source, Y_target, r, lambda_1, lambda_2, step_size, max_iter, threshold) + result <- learner_cpp(Y_source, Y_target, r, lambda_1, lambda_2, step_size, max_iter, threshold, max_value) return(result) } diff --git a/man/cv.learner.Rd b/man/cv.learner.Rd index ea3f5b3..6882414 100644 --- a/man/cv.learner.Rd +++ b/man/cv.learner.Rd @@ -14,7 +14,8 @@ cv.learner( max_iter = 100, threshold = 0.001, n_cores = 1, - r = 1 + r, + max_value = 10 ) } \arguments{ @@ -37,6 +38,8 @@ cv.learner( \item{n_cores}{an integer specifying the number of CPU cores in OpenMP parallelization. Parallelization is performed across the different candidate \eqn{(\lambda_1, \lambda_2)} pairs. The default is \code{1}, i.e., no parallelization.} \item{r}{(optional) integer specifying the rank of the knowledge graphs, as in \code{\link{learner}}} + +\item{max_value}{Numeric scalar used to specify the maximum value of the objective function allowed before terminating the algorithm. Specifically, the algorithm will terminate if the value of the objective function exceeds \code{max_value}\eqn{\times \epsilon_0}, where \eqn{\epsilon_0} denotes the value of the objective function at the initial point. This is used to prevent unnecessary computation time after the optimization algorithm diverges.} } \value{ A list with the following elements: diff --git a/man/learner.Rd b/man/learner.Rd index e05e175..98bb1ca 100644 --- a/man/learner.Rd +++ b/man/learner.Rd @@ -12,7 +12,8 @@ learner( step_size, max_iter = 100, threshold = 0.001, - r = 1 + r, + max_value = 10 ) } \arguments{ @@ -31,6 +32,8 @@ learner( \item{threshold}{Convergence threshold.} \item{r}{Rank used for the truncated SVD.} + +\item{max_value}{Numeric scalar used to specify the maximum value of the objective function allowed before terminating the algorithm. Specifically, the algorithm will terminate if the value of the objective function exceeds \code{max_value}\eqn{\times \epsilon_0}, where \eqn{\epsilon_0} denotes the value of the objective function at the initial point. This is used to prevent unnecessary computation time after the optimization algorithm diverges.} } \value{ A list containing the learner estimate, objective values, convergence criterion, and rank. diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 0eb47c2..43de911 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -12,8 +12,8 @@ Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif // learner_cpp -List learner_cpp(const Eigen::MatrixXd& Y_source, const Eigen::MatrixXd& Y_target, int r, double lambda1, double lambda2, double step_size, int max_iter, double threshold); -RcppExport SEXP _learner_learner_cpp(SEXP Y_sourceSEXP, SEXP Y_targetSEXP, SEXP rSEXP, SEXP lambda1SEXP, SEXP lambda2SEXP, SEXP step_sizeSEXP, SEXP max_iterSEXP, SEXP thresholdSEXP) { +List learner_cpp(const Eigen::MatrixXd& Y_source, const Eigen::MatrixXd& Y_target, int r, double lambda1, double lambda2, double step_size, int max_iter, double threshold, double max_value); +RcppExport SEXP _learner_learner_cpp(SEXP Y_sourceSEXP, SEXP Y_targetSEXP, SEXP rSEXP, SEXP lambda1SEXP, SEXP lambda2SEXP, SEXP step_sizeSEXP, SEXP max_iterSEXP, SEXP thresholdSEXP, SEXP max_valueSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -25,13 +25,14 @@ BEGIN_RCPP Rcpp::traits::input_parameter< double >::type step_size(step_sizeSEXP); Rcpp::traits::input_parameter< int >::type max_iter(max_iterSEXP); Rcpp::traits::input_parameter< double >::type threshold(thresholdSEXP); - rcpp_result_gen = Rcpp::wrap(learner_cpp(Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold)); + Rcpp::traits::input_parameter< double >::type max_value(max_valueSEXP); + rcpp_result_gen = Rcpp::wrap(learner_cpp(Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold, max_value)); return rcpp_result_gen; END_RCPP } // cv_learner_cpp -List cv_learner_cpp(const Eigen::MatrixXd& Y_source, const Eigen::MatrixXd& Y_target, const std::vector& lambda1_all, const std::vector& lambda2_all, double step_size, int n_folds, int max_iter, double threshold, int n_cores, int r); -RcppExport SEXP _learner_cv_learner_cpp(SEXP Y_sourceSEXP, SEXP Y_targetSEXP, SEXP lambda1_allSEXP, SEXP lambda2_allSEXP, SEXP step_sizeSEXP, SEXP n_foldsSEXP, SEXP max_iterSEXP, SEXP thresholdSEXP, SEXP n_coresSEXP, SEXP rSEXP) { +List cv_learner_cpp(const Eigen::MatrixXd& Y_source, const Eigen::MatrixXd& Y_target, const std::vector& lambda1_all, const std::vector& lambda2_all, double step_size, int n_folds, int max_iter, double threshold, int n_cores, int r, double max_value); +RcppExport SEXP _learner_cv_learner_cpp(SEXP Y_sourceSEXP, SEXP Y_targetSEXP, SEXP lambda1_allSEXP, SEXP lambda2_allSEXP, SEXP step_sizeSEXP, SEXP n_foldsSEXP, SEXP max_iterSEXP, SEXP thresholdSEXP, SEXP n_coresSEXP, SEXP rSEXP, SEXP max_valueSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -45,7 +46,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< double >::type threshold(thresholdSEXP); Rcpp::traits::input_parameter< int >::type n_cores(n_coresSEXP); Rcpp::traits::input_parameter< int >::type r(rSEXP); - rcpp_result_gen = Rcpp::wrap(cv_learner_cpp(Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r)); + Rcpp::traits::input_parameter< double >::type max_value(max_valueSEXP); + rcpp_result_gen = Rcpp::wrap(cv_learner_cpp(Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r, max_value)); return rcpp_result_gen; END_RCPP } @@ -61,8 +63,8 @@ END_RCPP } static const R_CallMethodDef CallEntries[] = { - {"_learner_learner_cpp", (DL_FUNC) &_learner_learner_cpp, 8}, - {"_learner_cv_learner_cpp", (DL_FUNC) &_learner_cv_learner_cpp, 10}, + {"_learner_learner_cpp", (DL_FUNC) &_learner_learner_cpp, 9}, + {"_learner_cv_learner_cpp", (DL_FUNC) &_learner_cv_learner_cpp, 11}, {"_learner_omp_max_threads", (DL_FUNC) &_learner_omp_max_threads, 0}, {NULL, NULL, 0} }; diff --git a/src/kernels.cpp b/src/kernels.cpp index d0cb3af..c0fd6eb 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -17,7 +17,8 @@ using namespace Eigen; Eigen::MatrixXd learner_worker(const Eigen::MatrixXd &Y_source, const Eigen::MatrixXd &Y_target, int r, double lambda1, double lambda2, - double step_size, int max_iter, double threshold) { + double step_size, int max_iter, double threshold, + double max_value) { int p = Y_source.rows(); int q = Y_source.cols(); @@ -108,6 +109,9 @@ Eigen::MatrixXd learner_worker(const Eigen::MatrixXd &Y_source, if (iter > 0 && std::abs(obj - obj_init) < threshold) { break; } + if (iter > 0 && obj > max_value * obj_init) { + break; + } obj_init = obj; } @@ -118,8 +122,8 @@ Eigen::MatrixXd learner_worker(const Eigen::MatrixXd &Y_source, // Exported learner function. // [[Rcpp::export]] List learner_cpp(const Eigen::MatrixXd &Y_source, const Eigen::MatrixXd &Y_target, - int r, double lambda1, double lambda2, double step_size, int max_iter, double threshold) { - Eigen::MatrixXd learner_estimate = learner_worker(Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold); + int r, double lambda1, double lambda2, double step_size, int max_iter, double threshold, double max_value) { + Eigen::MatrixXd learner_estimate = learner_worker(Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold, max_value); return List::create( Named("learner_estimate") = learner_estimate ); @@ -130,7 +134,7 @@ List learner_cpp(const Eigen::MatrixXd &Y_source, const Eigen::MatrixXd &Y_targe List cv_learner_cpp(const Eigen::MatrixXd &Y_source, const Eigen::MatrixXd &Y_target, const std::vector &lambda1_all, const std::vector &lambda2_all, double step_size, int n_folds, int max_iter, double threshold, - int n_cores, int r) { + int n_cores, int r, double max_value) { int p = Y_source.rows(); int q = Y_source.cols(); int n_lambda1 = lambda1_all.size(); @@ -163,7 +167,7 @@ List cv_learner_cpp(const Eigen::MatrixXd &Y_source, const Eigen::MatrixXd &Y_ta } Eigen::MatrixXd learner_estimate = learner_worker(Y_source, Y_train, r, lambda1_all[i], lambda2_all[j], - step_size, max_iter, threshold); + step_size, max_iter, threshold, max_value); double fold_mse = 0.0; for (int idx : index_set[fold]) { int row = idx / q; From 8e8ee20505f9d896cfb6276fb0505375cc627a72 Mon Sep 17 00:00:00 2001 From: Sean McGrath Date: Sat, 22 Feb 2025 17:27:46 -0500 Subject: [PATCH 05/16] Included the objective function values and applicable convergence criterion in the output of the learner function --- src/kernels.cpp | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/src/kernels.cpp b/src/kernels.cpp index c0fd6eb..9cc1cb5 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -14,7 +14,7 @@ using namespace Eigen; //-------------------------------------------------------------- // Internal worker function: computes the learner estimate - thread safe for OMP -Eigen::MatrixXd learner_worker(const Eigen::MatrixXd &Y_source, +List learner_worker(const Eigen::MatrixXd &Y_source, const Eigen::MatrixXd &Y_target, int r, double lambda1, double lambda2, double step_size, int max_iter, double threshold, @@ -62,6 +62,8 @@ Eigen::MatrixXd learner_worker(const Eigen::MatrixXd &Y_source, double U_norm = U.norm(); double V_norm = V.norm(); + int convergence_criterion = 2; + std::vector obj_values; for (int iter = 0; iter < max_iter; ++iter) { Eigen::MatrixXd U_tilde = U.transpose() * U; Eigen::MatrixXd V_tilde = V.transpose() * V; @@ -98,6 +100,7 @@ Eigen::MatrixXd learner_worker(const Eigen::MatrixXd &Y_source, lambda1 * (P_V * V).squaredNorm() + lambda2 * (U.transpose() * U - V.transpose() * V).squaredNorm(); } + obj_values.push_back(obj); if (obj < obj_best) { obj_best = obj; @@ -107,15 +110,21 @@ Eigen::MatrixXd learner_worker(const Eigen::MatrixXd &Y_source, // Simple convergence check (you might refine this). if (iter > 0 && std::abs(obj - obj_init) < threshold) { + convergence_criterion = 1; break; } if (iter > 0 && obj > max_value * obj_init) { + convergence_criterion = 3; break; } obj_init = obj; } - return U_best * V_best.transpose(); + return List::create( + Named("learner_estimate") = U_best * V_best.transpose(), + Named("objective_values") = obj_values, + Named("convergence_criterion") = convergence_criterion + ); } //-------------------------------------------------------------- @@ -123,9 +132,12 @@ Eigen::MatrixXd learner_worker(const Eigen::MatrixXd &Y_source, // [[Rcpp::export]] List learner_cpp(const Eigen::MatrixXd &Y_source, const Eigen::MatrixXd &Y_target, int r, double lambda1, double lambda2, double step_size, int max_iter, double threshold, double max_value) { - Eigen::MatrixXd learner_estimate = learner_worker(Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold, max_value); + List out = learner_worker(Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold, max_value); return List::create( - Named("learner_estimate") = learner_estimate + Named("learner_estimate") = out["learner_estimate"], + Named("objective_values") = out["objective_values"], + Named("convergence_criterion") = out["convergence_criterion"], + Named("r") = r ); } @@ -165,9 +177,10 @@ List cv_learner_cpp(const Eigen::MatrixXd &Y_source, const Eigen::MatrixXd &Y_ta int col = idx % q; Y_train(row, col) = std::numeric_limits::quiet_NaN(); } - Eigen::MatrixXd learner_estimate = learner_worker(Y_source, Y_train, r, + List temp = learner_worker(Y_source, Y_train, r, lambda1_all[i], lambda2_all[j], step_size, max_iter, threshold, max_value); + Eigen::MatrixXd learner_estimate = temp["learner_estimate"]; double fold_mse = 0.0; for (int idx : index_set[fold]) { int row = idx / q; From 9f57c267182f1355a28200698a1780d45b306d8f Mon Sep 17 00:00:00 2001 From: Sean McGrath Date: Sat, 22 Feb 2025 17:41:37 -0500 Subject: [PATCH 06/16] Fixed a small eror in the convergence check; pre-allocated space for obj_values --- src/kernels.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/kernels.cpp b/src/kernels.cpp index 9cc1cb5..6793429 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -64,6 +64,7 @@ List learner_worker(const Eigen::MatrixXd &Y_source, int convergence_criterion = 2; std::vector obj_values; + obj_values.reserve(max_iter); for (int iter = 0; iter < max_iter; ++iter) { Eigen::MatrixXd U_tilde = U.transpose() * U; Eigen::MatrixXd V_tilde = V.transpose() * V; @@ -109,7 +110,7 @@ List learner_worker(const Eigen::MatrixXd &Y_source, } // Simple convergence check (you might refine this). - if (iter > 0 && std::abs(obj - obj_init) < threshold) { + if (iter > 0 && std::abs(obj - obj_values[iter - 1]) < threshold) { convergence_criterion = 1; break; } From aed237c59fc0c307ce1248912f2e6ab79c241e99 Mon Sep 17 00:00:00 2001 From: Sean McGrath Date: Tue, 25 Feb 2025 00:01:32 -0500 Subject: [PATCH 07/16] Revised computations involved P_U and P_V and changed the SVD method --- src/kernels.cpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/kernels.cpp b/src/kernels.cpp index 6793429..074b624 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -22,7 +22,7 @@ List learner_worker(const Eigen::MatrixXd &Y_source, int p = Y_source.rows(); int q = Y_source.cols(); - Eigen::JacobiSVD svd(Y_source, Eigen::ComputeThinU | Eigen::ComputeThinV); + Eigen::BDCSVD svd(Y_source, Eigen::ComputeThinU | Eigen::ComputeThinV); Eigen::MatrixXd U_full = svd.matrixU(); Eigen::MatrixXd V_full = svd.matrixV(); int r_use = std::min(r, static_cast(U_full.cols())); @@ -35,9 +35,9 @@ List learner_worker(const Eigen::MatrixXd &Y_source, Eigen::MatrixXd U = U_trunc * sqrtD; Eigen::MatrixXd V = V_trunc * sqrtD; - // Projection matrices (using truncated basis) - Eigen::MatrixXd P_U = Eigen::MatrixXd::Identity(p, p) - U_trunc * U_trunc.transpose(); - Eigen::MatrixXd P_V = Eigen::MatrixXd::Identity(q, q) - V_trunc * V_trunc.transpose(); + // Precompute reusable matrices + Eigen::MatrixXd U_trunc_T = U_trunc.transpose(); + Eigen::MatrixXd V_trunc_T = V_trunc.transpose(); // For now, we assume no missing data (perc_nonmissing = 1). // double perc_nonmissing = (Y_target.array().isNaN().count() > 0) ? ( (p*q - num_missing) / double(p*q)) : 1.0; @@ -51,8 +51,8 @@ List learner_worker(const Eigen::MatrixXd &Y_source, diff = diff.array().isNaN().select(Eigen::MatrixXd::Zero(diff.rows(), diff.cols()), diff); } obj_init = diff.squaredNorm() / perc_nonmissing + - lambda1 * (P_U * U).squaredNorm() + - lambda1 * (P_V * V).squaredNorm() + + lambda1 * (U - U_trunc * (U_trunc_T * U)).squaredNorm() + + lambda1 * (V - V_trunc * (V_trunc_T * V)).squaredNorm() + lambda2 * (U.transpose() * U - V.transpose() * V).squaredNorm(); } @@ -76,10 +76,10 @@ List learner_worker(const Eigen::MatrixXd &Y_source, } Eigen::MatrixXd grad_U = (2.0 / perc_nonmissing) * (U * V_tilde - adjusted_theta * V) - + lambda1 * 2 * P_U * U + + lambda1 * 2 * (U - U_trunc * (U_trunc_T * U)) + lambda2 * 4 * U * (U_tilde - V_tilde); Eigen::MatrixXd grad_V = (2.0 / perc_nonmissing) * (V * U_tilde - adjusted_theta.transpose() * U) - + lambda1 * 2 * P_V * V + + lambda1 * 2 * (V - V_trunc * (V_trunc_T * V)) + lambda2 * 4 * V * (V_tilde - U_tilde); double grad_U_norm = grad_U.norm(); @@ -97,8 +97,8 @@ List learner_worker(const Eigen::MatrixXd &Y_source, diff = diff.array().isNaN().select(Eigen::MatrixXd::Zero(diff.rows(), diff.cols()), diff); } obj = diff.squaredNorm() / perc_nonmissing + - lambda1 * (P_U * U).squaredNorm() + - lambda1 * (P_V * V).squaredNorm() + + lambda1 * (U - U_trunc * (U_trunc_T * U)).squaredNorm() + + lambda1 * (V - V_trunc * (V_trunc_T * V)).squaredNorm() + lambda2 * (U.transpose() * U - V.transpose() * V).squaredNorm(); } obj_values.push_back(obj); From 2fbaad7b061269f5b8f519f2b1402ddebfb974be Mon Sep 17 00:00:00 2001 From: Sean McGrath Date: Tue, 25 Feb 2025 00:29:58 -0500 Subject: [PATCH 08/16] Added back some of the original documentation and parameter structure --- R/learner.R | 97 +++++++++++++++++++++++++++++++++++++---------- man/cv.learner.Rd | 16 +++----- man/learner.Rd | 73 +++++++++++++++++++++++------------ 3 files changed, 129 insertions(+), 57 deletions(-) diff --git a/R/learner.R b/R/learner.R index ec11040..362fa08 100644 --- a/R/learner.R +++ b/R/learner.R @@ -10,9 +10,7 @@ #' @param step_size numeric scalar specifying the step size for the Newton steps in the numerical optimization algorithm, as in \code{\link{learner}} #' @param n_folds an integer specify the number of cross-validation folds. The default is \code{4}. #' @param n_cores an integer specifying the number of CPU cores in OpenMP parallelization. Parallelization is performed across the different candidate \eqn{(\lambda_1, \lambda_2)} pairs. The default is \code{1}, i.e., no parallelization. -#' @param threshold Convergence threshold. -#' @param max_iter Maximum number of iterations. -#' @param max_value Numeric scalar used to specify the maximum value of the objective function allowed before terminating the algorithm. Specifically, the algorithm will terminate if the value of the objective function exceeds \code{max_value}\eqn{\times \epsilon_0}, where \eqn{\epsilon_0} denotes the value of the objective function at the initial point. This is used to prevent unnecessary computation time after the optimization algorithm diverges. +#' @param control a list of parameters for controlling the stopping criteria for the numerical optimization algorithm, as in \code{\link{learner}}. #' #' @return A list with the following elements: #' \item{lambda_1_min}{value of \eqn{\lambda_1} with the smallest MSE} @@ -36,8 +34,8 @@ #' #' #' @export -cv.learner <- function(Y_source, Y_target, lambda_1_all, lambda_2_all, step_size, - n_folds = 4, max_iter = 100, threshold = 1e-3, n_cores = 1, r, max_value = 10) { +cv.learner <- function(Y_source, Y_target, r, lambda_1_all, lambda_2_all, + step_size, n_folds = 4, n_cores = 1, control = list()) { ## Error catching if (!identical(dim(Y_source), dim(Y_target))){ stop('Y_source and Y_target must have the same dimensions') @@ -51,27 +49,75 @@ cv.learner <- function(Y_source, Y_target, lambda_1_all, lambda_2_all, step_size max_rank <- min(p, q) / 3 r <- max(ScreeNOT::adaptiveHardThresholding(Y = Y_source, k = max_rank)$r, 1) } + if (!('max_iter' %in% names(control))){ + control$max_iter <- 100 + } + if (!('threshold' %in% names(control))){ + control$threshold <- 0.001 + } + if (!('max_value' %in% names(control))){ + control$max_value <- 10 + } result <- cv_learner_cpp(Y_source, Y_target, lambda_1_all, lambda_2_all, step_size, - n_folds, max_iter, threshold, n_cores, r, max_value) + n_folds, control$max_iter, control$threshold, n_cores, r, control$max_value) return(result) } -#' Learner Optimization -#' -#' @param Y_source Matrix containing the source population data. -#' @param Y_target Matrix containing the target population data. -#' @param lambda_1 Regularization parameter lambda1. -#' @param lambda_2 Regularization parameter lambda2. -#' @param step_size Step size for optimization. -#' @param max_iter Maximum number of iterations. -#' @param threshold Convergence threshold. -#' @param r Rank used for the truncated SVD. -#' @param max_value Numeric scalar used to specify the maximum value of the objective function allowed before terminating the algorithm. Specifically, the algorithm will terminate if the value of the objective function exceeds \code{max_value}\eqn{\times \epsilon_0}, where \eqn{\epsilon_0} denotes the value of the objective function at the initial point. This is used to prevent unnecessary computation time after the optimization algorithm diverges. -#' @return A list containing the learner estimate, objective values, convergence criterion, and rank. +#' Latent space-based transfer learning +#' +#' This function applies the LatEnt spAce-based tRaNsfer lEaRning (LEARNER) method (McGrath et al. 2024) to leverage data from a source population to improve +#' estimation of a low rank matrix in an underrepresented target population. +#' +#' @param Y_target matrix containing the target population data +#' @param Y_source matrix containing the source population data +#' @param r (optional) integer specifying the rank of the knowledge graphs. By default, ScreeNOT (Donoho et al. 2023) is applied to the source population knowledge graph to select the rank. +#' @param lambda_1 numeric scalar specifying the value of \eqn{\lambda_1} (see Details) +#' @param lambda_2 numeric scalar specifying the value of \eqn{\lambda_2} (see Details) +#' @param step_size numeric scalar specifying the step size for the Newton steps in the numerical optimization algorithm +#' @param control a list of parameters for controlling the stopping criteria for the numerical optimization algorithm. The list may include the following components: +#' \tabular{ll}{ +#' \code{max_iter} \tab integer specifying the maximum number of iterations \cr +#' \code{threshold} \tab numeric scalar specifying a convergence threshold. The algorithm converges when \eqn{|\epsilon_t - \epsilon_{t-1}| < }\code{threshold}, where \eqn{\epsilon_t} denotes the value of the objective function at iteration \eqn{t}. \cr +#' \code{max_value} \tab numeric scalar used to specify the maximum value of the objective function allowed before terminating the algorithm. Specifically, the algorithm will terminate if the value of the objective function exceeds \code{max_value}\eqn{\times \epsilon_0}, where \eqn{\epsilon_0} denotes the value of the objective function at the initial point. This is used to prevent unnecessary computation time after the optimization algorithm diverges. \cr} +#' +#' @return A list with the following elements: +#' \item{learner_estimate}{matrix containing the LEARNER estimate of the target population knowledge graph} +#' \item{objective_values}{numeric vector containing the values of the objective function at each iteration} +#' \item{convergence_criterion}{integer specifying the criterion that was satisfied for terminating the numerical optimization algorithm. A value of 1 indicates the convergence threshold was satisfied; A value of 2 indicates that the maximum number of iterations was satisfied; A value of 3 indicates that the maximum value of the objective function was satisfied.} +#' \item{r}{rank value used.} +#' +#' @details +#' +#' \strong{Data and notation:} +#' +#' The data consists of a matrix in the target population \eqn{Y_0 \in \mathbb{R}^{p \times q}} and the source population \eqn{Y_1 \in \mathbb{R}^{p \times q}}. +#' Let \eqn{\hat{U}_{k} \hat{\Lambda}_{k} \hat{V}_{k}^{\top}} denote the truncated singular value decomposition (SVD) of \eqn{Y_k}, \eqn{k = 0, 1}. +#' +#' For \eqn{k = 0, 1}, one can view \eqn{Y_k} as a noisy version of \eqn{\Theta_k}, referred to as the knowledge graph. The target of inference is the target population knowledge graph, \eqn{\Theta_0}. +#' +#' \strong{Estimation:} +#' +#' This method estimates \eqn{\Theta_0} by \eqn{\tilde{U}\tilde{V}^{\top}}, where \eqn{(\tilde{U}, \tilde{V})} is the solution to the following optimization problem +#' \deqn{\mathrm{arg\,min}_{U \in \mathbb{R}^{p \times r}, V \in \mathbb{R}^{q \times r}} \big\{ \| U V^{\top} - Y_0 \|_F^2 + \lambda_1\| \mathcal{P}_{\perp}(\hat{U}_{1})U \|_F^2 + \lambda_1\| \mathcal{P}_{\perp}(\hat{V}_{1})V \|_F^2 + \lambda_2 \| U^{\top} U - V^{\top} V \|_F^2 \big\}} +#' where \eqn{\mathcal{P}_{\perp}(\hat{U}_{1}) = I - \hat{U}_{1}^{\top}\hat{U}_{1}} and \eqn{\mathcal{P}_{\perp}(\hat{V}_{1}) = I - \hat{V}_{1}^{\top}\hat{V}_{1}}. +#' +#' This function uses an alternating minimization strategy to solve the optimization problem. That is, this approach updates \eqn{U} by minimizing the objective function (via a gradient descent step) treating \eqn{V} as fixed. Then, \eqn{V} is updated treating \eqn{U} as fixed. These updates of \eqn{U} and \eqn{V} are repeated until convergence. +#' +#' @references +#' McGrath, S., Zhu, C,. Guo, M. and Duan, R. (2024). \emph{LEARNER: A transfer learning method for low-rank matrix estimation}. arXiv preprint arXiv:2412.20605. +#' +#' Donoho, D., Gavish, M. and Romanov, E. (2023). \emph{ScreeNOT: Exact MSE-optimal singular value thresholding in correlated noise}. The Annals of Statistics, 51(1), pp.122-148. +#' +#' @examples +#' res <- learner(Y_source = dat_highsim$Y_source, +#' Y_target = dat_highsim$Y_target, +#' lambda_1 = 1, lambda_2 = 1, +#' step_size = 0.003) +#' #' @export -learner <- function(Y_source, Y_target, lambda_1, lambda_2, step_size, - max_iter = 100, threshold = 1e-3, r, max_value = 10) { +learner <- function(Y_source, Y_target, r, lambda_1, lambda_2, step_size, + control = list()) { ## Error catching if (!identical(dim(Y_source), dim(Y_target))){ stop('Y_source and Y_target must have the same dimensions') @@ -91,8 +137,17 @@ learner <- function(Y_source, Y_target, lambda_1, lambda_2, step_size, max_rank <- min(p, q) / 3 r <- max(ScreeNOT::adaptiveHardThresholding(Y = Y_source, k = max_rank)$r, 1) } + if (!('max_iter' %in% names(control))){ + control$max_iter <- 100 + } + if (!('threshold' %in% names(control))){ + control$threshold <- 0.001 + } + if (!('max_value' %in% names(control))){ + control$max_value <- 10 + } - result <- learner_cpp(Y_source, Y_target, r, lambda_1, lambda_2, step_size, max_iter, threshold, max_value) + result <- learner_cpp(Y_source, Y_target, r, lambda_1, lambda_2, step_size, control$max_iter, control$threshold, control$max_value) return(result) } diff --git a/man/cv.learner.Rd b/man/cv.learner.Rd index 6882414..8c0bd0c 100644 --- a/man/cv.learner.Rd +++ b/man/cv.learner.Rd @@ -7,15 +7,13 @@ cv.learner( Y_source, Y_target, + r, lambda_1_all, lambda_2_all, step_size, n_folds = 4, - max_iter = 100, - threshold = 0.001, n_cores = 1, - r, - max_value = 10 + control = list() ) } \arguments{ @@ -23,6 +21,8 @@ cv.learner( \item{Y_target}{matrix containing the target population data, as in \code{\link{learner}}} +\item{r}{(optional) integer specifying the rank of the knowledge graphs, as in \code{\link{learner}}} + \item{lambda_1_all}{vector of numerics specifying the candidate values of \eqn{\lambda_1} (see Details)} \item{lambda_2_all}{vector of numerics specifying the candidate values of \eqn{\lambda_2} (see Details)} @@ -31,15 +31,9 @@ cv.learner( \item{n_folds}{an integer specify the number of cross-validation folds. The default is \code{4}.} -\item{max_iter}{Maximum number of iterations.} - -\item{threshold}{Convergence threshold.} - \item{n_cores}{an integer specifying the number of CPU cores in OpenMP parallelization. Parallelization is performed across the different candidate \eqn{(\lambda_1, \lambda_2)} pairs. The default is \code{1}, i.e., no parallelization.} -\item{r}{(optional) integer specifying the rank of the knowledge graphs, as in \code{\link{learner}}} - -\item{max_value}{Numeric scalar used to specify the maximum value of the objective function allowed before terminating the algorithm. Specifically, the algorithm will terminate if the value of the objective function exceeds \code{max_value}\eqn{\times \epsilon_0}, where \eqn{\epsilon_0} denotes the value of the objective function at the initial point. This is used to prevent unnecessary computation time after the optimization algorithm diverges.} +\item{control}{a list of parameters for controlling the stopping criteria for the numerical optimization algorithm, as in \code{\link{learner}}.} } \value{ A list with the following elements: diff --git a/man/learner.Rd b/man/learner.Rd index 98bb1ca..e1bdce4 100644 --- a/man/learner.Rd +++ b/man/learner.Rd @@ -2,42 +2,65 @@ % Please edit documentation in R/learner.R \name{learner} \alias{learner} -\title{Learner Optimization} +\title{Latent space-based transfer learning} \usage{ -learner( - Y_source, - Y_target, - lambda_1, - lambda_2, - step_size, - max_iter = 100, - threshold = 0.001, - r, - max_value = 10 -) +learner(Y_source, Y_target, r, lambda_1, lambda_2, step_size, control = list()) } \arguments{ -\item{Y_source}{Matrix containing the source population data.} +\item{Y_source}{matrix containing the source population data} -\item{Y_target}{Matrix containing the target population data.} +\item{Y_target}{matrix containing the target population data} -\item{lambda_1}{Regularization parameter lambda1.} +\item{r}{(optional) integer specifying the rank of the knowledge graphs. By default, ScreeNOT (Donoho et al. 2023) is applied to the source population knowledge graph to select the rank.} -\item{lambda_2}{Regularization parameter lambda2.} +\item{lambda_1}{numeric scalar specifying the value of \eqn{\lambda_1} (see Details)} -\item{step_size}{Step size for optimization.} +\item{lambda_2}{numeric scalar specifying the value of \eqn{\lambda_2} (see Details)} -\item{max_iter}{Maximum number of iterations.} +\item{step_size}{numeric scalar specifying the step size for the Newton steps in the numerical optimization algorithm} -\item{threshold}{Convergence threshold.} - -\item{r}{Rank used for the truncated SVD.} - -\item{max_value}{Numeric scalar used to specify the maximum value of the objective function allowed before terminating the algorithm. Specifically, the algorithm will terminate if the value of the objective function exceeds \code{max_value}\eqn{\times \epsilon_0}, where \eqn{\epsilon_0} denotes the value of the objective function at the initial point. This is used to prevent unnecessary computation time after the optimization algorithm diverges.} +\item{control}{a list of parameters for controlling the stopping criteria for the numerical optimization algorithm. The list may include the following components: +\tabular{ll}{ +\code{max_iter} \tab integer specifying the maximum number of iterations \cr +\code{threshold} \tab numeric scalar specifying a convergence threshold. The algorithm converges when \eqn{|\epsilon_t - \epsilon_{t-1}| < }\code{threshold}, where \eqn{\epsilon_t} denotes the value of the objective function at iteration \eqn{t}. \cr +\code{max_value} \tab numeric scalar used to specify the maximum value of the objective function allowed before terminating the algorithm. Specifically, the algorithm will terminate if the value of the objective function exceeds \code{max_value}\eqn{\times \epsilon_0}, where \eqn{\epsilon_0} denotes the value of the objective function at the initial point. This is used to prevent unnecessary computation time after the optimization algorithm diverges. \cr}} } \value{ -A list containing the learner estimate, objective values, convergence criterion, and rank. +A list with the following elements: +\item{learner_estimate}{matrix containing the LEARNER estimate of the target population knowledge graph} +\item{objective_values}{numeric vector containing the values of the objective function at each iteration} +\item{convergence_criterion}{integer specifying the criterion that was satisfied for terminating the numerical optimization algorithm. A value of 1 indicates the convergence threshold was satisfied; A value of 2 indicates that the maximum number of iterations was satisfied; A value of 3 indicates that the maximum value of the objective function was satisfied.} +\item{r}{rank value used.} } \description{ -Learner Optimization +This function applies the LatEnt spAce-based tRaNsfer lEaRning (LEARNER) method (McGrath et al. 2024) to leverage data from a source population to improve +estimation of a low rank matrix in an underrepresented target population. +} +\details{ +\strong{Data and notation:} + +The data consists of a matrix in the target population \eqn{Y_0 \in \mathbb{R}^{p \times q}} and the source population \eqn{Y_1 \in \mathbb{R}^{p \times q}}. +Let \eqn{\hat{U}_{k} \hat{\Lambda}_{k} \hat{V}_{k}^{\top}} denote the truncated singular value decomposition (SVD) of \eqn{Y_k}, \eqn{k = 0, 1}. + +For \eqn{k = 0, 1}, one can view \eqn{Y_k} as a noisy version of \eqn{\Theta_k}, referred to as the knowledge graph. The target of inference is the target population knowledge graph, \eqn{\Theta_0}. + +\strong{Estimation:} + +This method estimates \eqn{\Theta_0} by \eqn{\tilde{U}\tilde{V}^{\top}}, where \eqn{(\tilde{U}, \tilde{V})} is the solution to the following optimization problem +\deqn{\mathrm{arg\,min}_{U \in \mathbb{R}^{p \times r}, V \in \mathbb{R}^{q \times r}} \big\{ \| U V^{\top} - Y_0 \|_F^2 + \lambda_1\| \mathcal{P}_{\perp}(\hat{U}_{1})U \|_F^2 + \lambda_1\| \mathcal{P}_{\perp}(\hat{V}_{1})V \|_F^2 + \lambda_2 \| U^{\top} U - V^{\top} V \|_F^2 \big\}} +where \eqn{\mathcal{P}_{\perp}(\hat{U}_{1}) = I - \hat{U}_{1}^{\top}\hat{U}_{1}} and \eqn{\mathcal{P}_{\perp}(\hat{V}_{1}) = I - \hat{V}_{1}^{\top}\hat{V}_{1}}. + +This function uses an alternating minimization strategy to solve the optimization problem. That is, this approach updates \eqn{U} by minimizing the objective function (via a gradient descent step) treating \eqn{V} as fixed. Then, \eqn{V} is updated treating \eqn{U} as fixed. These updates of \eqn{U} and \eqn{V} are repeated until convergence. +} +\examples{ +res <- learner(Y_source = dat_highsim$Y_source, + Y_target = dat_highsim$Y_target, + lambda_1 = 1, lambda_2 = 1, + step_size = 0.003) + +} +\references{ +McGrath, S., Zhu, C,. Guo, M. and Duan, R. (2024). \emph{LEARNER: A transfer learning method for low-rank matrix estimation}. arXiv preprint arXiv:2412.20605. + +Donoho, D., Gavish, M. and Romanov, E. (2023). \emph{ScreeNOT: Exact MSE-optimal singular value thresholding in correlated noise}. The Annals of Statistics, 51(1), pp.122-148. } From 047c7d8331749c4a324e165eb36ed5a9ec630205 Mon Sep 17 00:00:00 2001 From: Sean McGrath Date: Tue, 25 Feb 2025 00:42:50 -0500 Subject: [PATCH 09/16] Fixed name issue in the output of cv.learner. Edited the README file so that it does not display the package loading message --- README.Rmd | 2 +- README_files/figure-gfm/unnamed-chunk-7-1.png | Bin 0 -> 26109 bytes src/kernels.cpp | 4 ++-- 3 files changed, 3 insertions(+), 3 deletions(-) create mode 100644 README_files/figure-gfm/unnamed-chunk-7-1.png diff --git a/README.Rmd b/README.Rmd index 9f89d8b..595e3ff 100644 --- a/README.Rmd +++ b/README.Rmd @@ -39,7 +39,7 @@ devtools::install_github("stmcg/learner") ## Example We illustrate an example of how `learner` can be used. We first load the package. -```{r} +```{r, results="hide", message=FALSE} library(learner) ``` diff --git a/README_files/figure-gfm/unnamed-chunk-7-1.png b/README_files/figure-gfm/unnamed-chunk-7-1.png new file mode 100644 index 0000000000000000000000000000000000000000..360e568a546e5c5141abe5a4a0b6b87627aea0d3 GIT binary patch literal 26109 zcmeGEWmr`27e0!UQc4OcEh$n4Aw4RR(ka7$N)Fv!0}2)$qmt5%QX`F^q}0$20@5j+ z=NSP%-}8@m=gqmU-wVBTX3ySFtY@uz-RoZawW_l0C1P4)EG(=`@^aGZSXekoSXkH- zgcrak9`D8lu(0q*;SV3GT0WG0XlrHbpkepS$VAq}#>4?`q%Qja3+rBRuoldmUhA@C zY*`UETYErHdU^!uv0zlOm3r5Up}Av|BI^;JNUP%eVa-dU@o?#z?}x z@<3z8*l~Dh^Lb~we-$6>@wCGSpN?vUu4(Hz)xQ*|{%%d5-{WBuExofEUg zd{tTOSEDZna__aU{4&8hrv6a~^TWG;&_vY#e3$6Em{lC}FUtB7ES;+fWNU>TAHwgv zGJM2yBm3R)_u*Frv#i^Ccm$4X!y)2=nR#3+t&f_8ay}}vc2@-)_?A8zwfsIDpb}bi z%lq;;?$>SPO^2AtAo4`5=@F+_1$^+LALJ({WpJW(RxDUW= z!dG(5-_TzpO-jR&E!c@QiL^7S zkxY@x`jX$N-v>8CzP)f?qQ3Q%N;-_04DJPenvy|GJ7@lBF3tSYZtuj>cpcZ?F&0YG z;<#yXuItu=^b^*#&g{BXb1l}&i$3n=InJF#rgp47x8?ceBLa5mSt6%zUAvyb9KqSF zVf@axOH_eWurY?k(fChUSC2b3t9qjF+se0duJ;>`&5Ykq6nFekNE9X?=osj*Y>Do- zF%b2((Ro`Tqn%jM=Q!`gW=3nl>(n&V6YtP*QfD~Z|JiU>SnuL7fe)_l;;=SKa8cl= zHnDfPE%dwW9|>h+$uwRqt7kn7vAoONDhSOkKo6&Pnx!Hvn|C6ZIDDtZ zrX?ra9TOcA?Z^F^U=w@}h7N}IOpb}N7Cr8(0+uTUDPM?x!;xwui4XJ!c*SoRzm%w} zuQOIunWGMTLG+B@!$mlUD)T)!(>j!QPVMSYV}nldpLIK9%RR`22} z+f2P+m6y84^3uA(c8oz(01f7 zEp9>GY}m)J2E3{aP7aPdy#L_#jAMOGI>9|*cO&^>)Yl3e#61&Tc~d1NtefE9gjm=? za4bCVFKqCK7W~2Z5}fl-aFUX6&;N}*f%%fsaCig@3xXvteP7cJduf#DP_z4_aou&O zJDJ5=E{C1BxRpPg-7+RuO+n3=Q*3Y}GINQNiz|v3%}>qp;FEEN1eiK)jnf?C$<)?yEv8PG; zRtFnL)v%1~;`t3Beppef|1*RF9K57z2aT)$28D%vWBlL6SlEX6gb?qU^jl}+0fWL( zt;Ri{0Txa*F}dG z{Zs=Nw7+)8Pw((>HqkfRyem0dGb_Tj^T!)%jlGST-fS)I@s?EPb8p2Q13B8?GFpz< z+`$#O;51lhvd~|UXwjRUM#-*_sStfLTfeNZpAO!2T+(ll{HlrJvjKQvPt~5xZi{PpXybgB6X&5C%bqX}PV{aj6b(H;Vs$A%cUbz&q*H zF`Ta>;r}VTE70)k$J-nz%I2d0T-`yFi~7Gy4H32YN!CNfUGJ{l-LO`j4u*H9XUPQ9 zhG&eh{`01(uOJhA-x6KN+zXkVzzT_;9R8`gogjn_M~%AtY4GFI|EiiI8`}5%HI)f+ ze(*mxt**#yhix9op(NnCq7QRuW?vX8$tkrRPf#p=7)1I}+`Vw9#1gqEVep@?mdl4s z;HXYzLs9v{_1=Gex9L^7eaY>2Y>nZ~HmEAM=*yi?GgjbRY>&42QoPt?@-L{N^_au3 z&8Oe#lBy~e8Mi*$kAEep9^^K+N;6{d=tty4 zbEj3TC{*zi2|hUYn!mWnzcO|@DM|ab@}f!jC_-}{>cT`ZyVEU5 zQH398qZ}PbhFqze>SjAML5E~>!d90d$8BFRyzu5mBN-8){!GG5*3^SdP}N@o!> ztQj*3s`ytvTMNAgHeI@g#KX329|ga-aBm*;X3HaC91^=pPK$ZyVcRSFj=f4)T7OAm5(Bh>7p)6R~^UeR-P73Bzh>=h^Wm+s-*Vik$B>5 z%bwpsT?e!3FS;qaptmegm*|A%Dpsri^aYF?g-TR@-wZz@cJy7y%4o6qQtB|%;-16z zW60KQY}`?0K07D1VwH5qUpc<~uI%L9^VkouAjY$yL?_)3a2mEonRL_Ndg^ymD?8P+ z?PHeD(Y{4z5&}izRT_`8Jw0TZ_J%qQ-d%_C(@I=mIXBlmVmv-9aLDpaI}+l1lV%e< zm4z!lh*=MpI*SMdElvh8ZTI=F_GCr!m}sa9qslY$Yb5nOW8eRWs#4zq8|;i*s}H=S zk)xd_VA68gd39_QR$i@;=(C?idT#>j5~KLN*%gF`&kt(!m}hH=tjAO^L@U1Ke@7vS zRq|mDR_@f?Tbm%^Y0(}hQf5GgCaGBE6CPpKH-ct&Zv+(n!z}l&5E^5Z+fT_@4;KC? zLJ>U#7T@a3xOC?zK`54UBKB>HR!+=jn=pE=3)}W$@4p?0!;F>X(VMxWopEH_1KY0T zm%@xrbDf{|H(L3tsLGw@)wfngp(!%K=qG0J)Vq;ZijLc>Og4cgLHj3P#nzUSvCplj z3}Z!LhH(A?uZ?C-n9t9?p<;OC8*1)!L!Z4LVsC9LelBk(!2DOcQe^ewp3ptBTlnQwk{1d~yJSfToF;`G?ywIRvqoNzq{wjK_R&ulRp+ORzJ8$?- zK@DHL%SX+UZ2ix?(Oi*D=zBtEx53j&0&bdC#36M_QVUhP<@cKaUX|o~xE|Dk79m&D)s zN;>)aPjCrHHz~W2yin~%Mo+6D{7GO43nE;%7B>^NbuyQ^rGG}s?= z^0U=2_Jd4#K+M`crcH%PNiO#2vUs%OTNpu7P=DIR_&XKtWl-${fomDcRAa5$XyM)m5l?)QHB?;PwtCOBkKNs~i(?x~Hjf}+y!x(4wMlQx zF;u^nS7Ik_^JwSg;jM~)w~+^N3~Q*sP|EvgtKUY(O);23OdCXw^*(Rg{`~krGNyc_ zq0f`}ktqAv_fA{AI-DPbZ~y4e;Y1$G^UcNBlyluSvD5tL($sJU3!t3YQTG^KYmT2KqS-yBDauQ{jDO7gZ<>_`S?Tnh6$ z=LV|GY-a?)<&3K(u6UFPJnJNm6;CD<5q0TbcyY(Z>5P(kB$oelb<80LBLC;Ne2UDn zo8j0gA^{#${k|WT0CjsEU(YA$*K7Aw#lBmkn$t3dUIwhIoi9`^0Z5jZ!l+>1_huDs ze$~YAga=-fcY@P1oj+(m5IhL>NAb1FFggyGATZ~-qI16f2g=@L%72EK05UZ9LcG-< zDWNgYWKD;`Os9U>w;U-!nNvrO_({QxrIg}$REvwYS$U^t0S>%;N?_OJvth__gJ%ec z#s|@}I4^Qo57tq}!=3ct#HXtP=Y(K=G|lMYt3mSdR&Bi5oX>#w-TN4DmLkHhgRf~a z*2iUfGARs_yA8^Y_k9g#DbALXIqW7FU6upKFI_SyolikkouX6qx(q$p&rkot<2^!h zzXt5vy`*aSaZZg4lX_2Xl{;Hz_lE7C?H+`KgxoKEesHx{(V`QwKS?m(%-LNWEl)Am=( zWYfP@e1=B|AxWXG^(MGP%$dO?s`*h$%I{KA%ent_0rP;nDijH4!ZiT&N^6TvqqF^a z%-Hd79>!oE91lll@2r-|-qhfTi#whQ_2!HF_3r|Da0?L#dOeFhBa|(z<~J?ZWrap< zG{@OK5Wc~}Hdk%$;p^+h>9pvrYjq4^sCY{-HFZA5T0ZcE6&!;w$?}7h2-q9ca(Zm4 z?=g|oZ^M5UWeY5-=0~Ki0m?=_k^`OmD55RI>>rMb!0gt=Ew<>gSh3*%f*=y6Taf*& zem1Qa&nZq}8gmNB=kG3UqAWS%?E{N_s+|6`UF&oE`1WkSxiEL2%=k90aFPB-G!IuW zJ(u@!S&fMX?fF`R?-ApOqa9>^wP;)jfC6{GvYdGuh&s4?WJ_XY-SUQKZ*~NGIw!Bj zMmsM1>Gxr;!3uqF>x0F%Mn|qx36Jdh9Hj=+-M#E@dUgZ9Rd7XhWQH;#CT70#VcL(k zX#`A<<2}@<&W#x?Kzl zI;=48vV%TsJ6`}o@PzyK6p8XvU_SgXJ=(i6`*SIHR~jbIFPOrt&@ngjaz%kLRk+pS zu+C7m#Nm72v+uDG(ilT{An4avIe7tIbCg<89T9km1(sK%8LZl3CYNUq| zJaauigv<&cGn%VK;8rDcaA~Ax8aGF%Dedh{dhd*8gfNI%OgDuu*owRyaFL&gP)s#T z2(R9K;86TnVTfHx^4;0s{_ja%u?EPL>_Pge1oEw47aQ+0dptk~kSR!=&<6!L5S~4_ zN%wyezUK*_&(giyJ4QkNM8y{DS9%mGNZdis3n<=JEQi11b@58Y?uIY}NM`3>{` zbV!ZUO~I+cTpUN8gf2Ku*1e>a@F?a`j=NFuxPwLtadNEVapYW`>K6BF6u@-yovxP0gk{@DFuV z5iK#JdsFN3=ks`snTN~Ag}%JSHoNg9O3k6%T@lCMZV2;$3zsyDOqKZAHB1xQzoW`^ z_~*-JQ&ps*SH3=^6}Yi(qK02vc`JMvFHJS@z}97qD$- z&0sT#I%OJo%(D+pVZbbY48P%6>~&j-B%y;Uo>518DGI)W=mfV17GwG=^e-A;GCIhO zyq)wr4+sr7w{bh!ZEmkk!8T+SH^xBBdJ=~G87!XE+(U*OEzWWpG?~OIcD0?aC#zR< z;ixg~IZv$onaS^Bl^CxGeQRuhX&1BGM^uIr2sCp}oV=|$u$hn6nZ85$g1^~yHGy_3 zndba~`#86|Zw=}#3gzFS_E(u=BQv+~%F1Ou{1jB@Y@d^CG1JnMd8Onm+eY+Ud3DVX zq~ivEW;^{Eno*Ve{Yh_3A?4+j& zh#o8zRgHB5YQt1)`;zcT*s?oIViZ+M|8tLr zg;~$35HV%T{(PmZbm>`a7+?ob`&}`-P|YfjEm%A5{-9YRd}-jz?cs{%Jo&N+L})^{ zl5q0}`Eeg^i8Vi`R^Hot1K8*0sLO1MMQd%u^}*+O&kQ&fr7c)LvMogAw|25Lku8_9 zB85pk4oE+XTYJBZd;;@T`8jpjp6K>lYVOO1KaDOp2Su>gm)#E)RVbfjEc--TstZwt zKCC()Zx;^!na%9iW<1D8nhxsb$`+#!n3wez$?f}P|6f4N-wZ`%16 zRMQ)c>H30M=jjRC6XxH{?M%LhYJhUka9X>T<%qHw-E3s?HUH+1AC9E-QDNDCLA9!m zy~6#8H8kqLKu*$hb|_8s&$F{xt?7V$&`XiIWV>}e)-TG@vlV7r(kulwq-{^TE+O)j z36`O}&ekaY5tfCNr{;UeUH2=XYwVn)l zfi%i$boU`Ch{R~N^X)DM9!ga#d*()8emIF=$GatS~h3ur!L?y!U1t z)n=FM%tn(!q8AIw?Y@F$u)=v{rNrO32+7NE1bO_zU{OY2?vuz88ZV!%%C3CW{?{i* zE9HxP!PoA}%(O&zgA}<3Wo!5(L^bwu?^k;xz;Wox(%ZiXElJL-#9OIQ@tC?)Hs%PQ zQw1J*W9V)F`H7Ci+qeWdIT%0jpn6Bcr+>qI1q+vDbYyZzf+ z`D;cNu^EJI;zA_0by2h(#{QiQ@%Na<9DSL-+Sj3*ITcZ3H(m2nkGB7Odiv(vPeoBT z0bdAdw<%o)?9%@Bux-`WxMS&aK%0(z?Tjnn+$h`fNN_^Sd;Dn=u6mw;VK6R5;ird} z4CV>hDVKPlW@xQ>k%t6Pj> z1ubKQIh+6$>|PoB9+~^|Z{j^6Bt6liq0&Yy$!e=x&Y`H+5krgPvR9?fJ|FLzs+OXj zqTGXX=2{|H!_3I$q)rY~Ncbx|p(oyeUEbfDO&l$JV*~q<*`|Jb@SU4tw`N!477d-C z#dNUz@>o@cFx@!Eqrr{1hp-n$Q&QGN+<;bLXoR?`C2^M7J%vY;4 zQf&S`?U1^uh@c0}^MpU8?GCDSsPuP*L`e0ZBp|(d0kK@QbCm5;<8xFQrOFE}?HTxD zWWvqLU{{N60Wx=}yV+L*pFI=jg`Vcq%>6)%w}s2k?;y6pyMS%+R?Si&j+YnQ8t7_B@^8XwC0tU8!NV(D{z z;n>uD3&OA0D@1pnmPzT?q3*WJ1$(HGoI3#(<_^ZrxnI14hhf(rZvJ}9Z*hFM_r!K< zzPqskFLBF1l#CSk$^yw&B!54kG4lHjePJC^CzWH~JFo&}L$CFh!X#sjH`O$;@wW$M zZRGpc+WSZrrq(8L`oc!fN8bfPh-!{snm9lXu|6k9t_6}p4GtmDiHc|SfdMram;yw19iqugm z@Usv=SpR8&>Iqss&V;S_%zX+10S9FcK7F0VAgiQ^y|ZOM<(!62tH1lp+Eo19_>4n< zS9bzF-GJW+)EXSKoaWJw=!C2~0dMaJz!2e*T#_+L5)kX@apN){UD-5 z-!DR}1Iv`8^$LyU^IdAT`k)32DYU?*JVLLKzhZyd>QHWcVWol{=4qW<<$lf^pdSH0 zHw0=OwrQV8_`cRcnS6kUG%!$whRIHB++2YX2y67QnxpCfTf zwQFInFQvTJf1rm-#N3%B&nwOdiY8$)b0-p>xj4Yv=I(UIpo*!3me7aC2dhj~bfFXj zTnND07 zV_Okwo%18_F=i8cWV2YOW7h>)?90=t_S$<-U_yZBA~Ug&Q;+;LIxc45?-85yk4B#7 z#H@>WQbpN(i_nGQ=Re-eb3)##%P`;J9^?Bl%wpm5nxA@zxQsZvZ#D7vN<3S42=A$ zLa!t>h$&|2&Y=nNrKfa`&iG5*jdT7=?i%KF*f%oFYG*ui{}?4=bWJmyQ6jb95IN?d zSlGB)8F@cABkL4lrhb9^yqyFearsDtg?Lt<+7g{6ADNnwKn!`A#d0>3@0VW<;)-_X zg}O}z3$Ls}R%=d*ZJ0ej`P6vNxBX&x zo~Q$_-s&<8GA$FL1*yGhsuHWgJJdt>8PK{(j%H=A?a=#M=$BU%hD+1zjo@dT#WiN~ zB7$O&Ez#aJeauJPN(QQlhP(!q_xe36{UDY2jNY0bC1HhH5f?6o*;EORvD=BCyER2- zQ!;ZOYOY@nGy1bME2Y;qfcdoA3LD5W^Wl^d_W$Wnj7iS&^_~CIG3B63dTtEz6chU2bOWQ5YXHCS@dLvFMB62huuINtVE7^>+J6@uBEd5kF|cI1i63@0HNi_ z`-_exEgUR@aJ2Npz>CdWjQrqt*utF306n(FDym3w#k(bgX#4tR@YGqo1Hyv00ZfdW z4rbcg@ih;)(lNGFrC1@ZeN%~~?sTQ@>p}DpuRQD@<84rXE{NNDF=R@#aI(w#A0@(k z703NTl2DfU_R^pMNwrJA=<3_UQZkv?Op^UC>t3@&@0eB3N(2~~mC)fU3iTA3$uv=3rS2IdXEIyJIoGpJ z1CYh?(jDKLO|;_2+eU%aC&!+)b!w~!N*#d^Vxpd0E{Q+yoEok(w8=^R$a^OkyWiHH*IM(moSBzHpqM*)j zW)w{vVY62|=OjX$Mir(SLqmODHxi~s)lPhKn6FKB|1(V^E0>k6cJ2J?hv2Y(Y&3Ht zLPSyC(2(KN2+&M^u{mIbenG^gC2{kdXcYncq@&-Ozv5E z!PDdOff-P5L}nrnhunI#(Dse;zQxrDe( z^>vmp%h?jX{(L_9t$MD5=4rhbaS@a2D%%IlTXD6MNgb-vSkZ9Pw8MM@p+3K(bquM*eb^p&hylDmF01DS_MF}2CpEasF+jGWyAuIgpB}u_R{&#j zowxDnu6s%HH&nAhp8`H_vcA6LbDC();hjOjWUsS*yoK?FAf&D45!HZjp#R*cF=Wd_ z+7cX2Obwiw8qX@L9i)65ZSrP`yhlf_H+$ZQp3Q#{GyjwjQ3Vj`dWbmVCivhX_GQ$l zR*y%v70W`aE(x(SYIQM5&mJ;A!e{i5E0o?^v0TifOHA%f)jM^s#_uq>qiJ=`_K@We zC~C#3c_kqIF)2=YW7~!~ptk}qUW>S`D*z0LT~UF*h9{bvvuab_t*rfN=?lFx{*@7D z^4g=@LugP?0vJiL$gE@ALn^PayWu6H+XUjaWb-g#Y9OS_V{$#oEFNHvyq@f6r-5f# zYQh$$@d%QRK`a=@TLsEj7)5{;I($?eUUY5Ox4h&EIKSCIRFS6Mo2mGQ$Z3TK=2Nyg zml7dy{Oj!+&)(%6j^6u!vsdP2d<+|5Q27-D_N!J(J&`nDGkIvD}S_Lf2vzr7USu; zgmYl9!uHV2G9dpvP&F_{I+8U5ha!;>$$j(X*$CdDT(_E11TKMxL}CncO==rovk^#z z4l2t-oZkBH00l^>gd0~RzJ6obVlz_IIl`v$wm4(4QtV@zPs^vfzsn7>WUuxGtHHuY zrQ#WFZ&2&VJ$!_S)Rlr+L+`C>GF{rKXP2$Zilyv$gYHDzx$oEel|i5Fv7#1AMQFIt~V}%vhfCA znbgjRqv$9vzectu*Z9qAu&^1HlK-C?l%c|;|COLLsnqtpMO7vRG2X1&NdN^3wC_pm z_*iAg86;ZAggd?a{a>HZy2yKud+!A6l!W!eY25Hp1%=@`+m9#r1E(@w&glD&TPB<+>H1B}Y|GChc)0=c37id20$BaA^ z#(8C~+55&yZ`yOt_$fc#wJV!3VUN}6Pa$AC-5v~h5ch2RGUwz&nV-t zr%Tyb$CnsD6}i*6G4dk{P11A>qxFiz8bU4=P}^pGJf~Ket1%#&j>Myp!POkxpiH5( z-^#rm2)Oo3NyT2zImTZwc$HF%P0*CL%PxG#7MD7v{n%lV_w#wGLysXdcvmynJyBad zI-m@Y=zIKxGtTZo>2qEpQ04b6*qCX}0iJbmOmUh@QF^?@vLA;%*>^fzR+oq?oRd{4 z)LMSQwo0oUFsMsKGsn9TP;YJ}I>n+>j40<1Nih0b^gk|^m;9k$C=E3mWt1JaN-VY6 z|D0E5vwBbx{%VYW{wH6fww=U0TiMu=ooyzab=`R)V^4|6!$U@^JY82!jl>PzCh;#Z zhzck_S{ZsKCtYLd<5m6?XInV>V90aoJ}jexuUgjz5#a9~b6(|m+KA)6?I+`#`Z6~B zws8~r6UNTHKeJ|4y%{P@tCU`SB$5)Z*pgB7NiJMt9@5sOJ`%{UUqC!viI^W8`uq&{ z^_;gT)4cw>)|Iq=gUFWX$b(}=GKr}6yfo{xL&h0q3sNm=9_sFc)tc2&zyR>lorsvV z(<47ud}8!kyAN~j)#|F?JhW$_>H_Ko8SjJb;K2)YOqO>y2GLMZXHT6?M7H01SU!1S zej;E-Ju=NT7S3?$$6nLi_~qZ+8A7}_pOIJz`2$X|Hb`iTmMxGow&n92{8E|0WJ<U|RUU-}+WPAFhAqRT295;dB!Km!OsUiU$e%;|s25X=N2ZEVH{O4?mu-I;X1 ze7SVEkXra@Ef)VJ&C~=Jd7q8x=D}^QxyT51#g>?jcg{txb;RhVX^Q!lkfJzzR3G|D z7Sw+&fHH82rACHL&*||-Bw^wfncZs&YZf&4r9;Vb9%^jE|MTBen_rFpXS|pW2#4*W zhmIW*8L_t2gR$XKBfKVm*6PVTzZqpH#~X9Sf*du-Pf4LraecdSnoLR#5bf6EjG?TK z7f(UL9Sr@5Nz1kJ^t7(%+aY+kjzDFisxnDPdqW$^^9A6on?O{Ao({Z3$B{D?#iddC zHg>0K|6^Fc04GSwPU{ z0*bYa`dLVy|bhZrQ1#Q_C_=hagP9K|A5d>b#w9KFi?&*~Mogz9&bK$)=W_wm{W{ zkNv|r8(L}#B$Qz%H+%OgM`Wa2xUCAg=z#PBXE z&7S~k<0?Tdj;0hNR<8@f6d`9hJ|I+cubWU za>8?SPHmvRND2p^=JWSW&N${gT?&KJ$PerI=|-}cXJJ67JJFanz37jb_`!1L9O3%I z7^uy9Q|^@>oXVQxeMOqf`z24R)u5M}Z9p7tz1g>5Hr&jam`kN1t_FbLPmfRnLbeCB z>TMwm+$dV4fNjFaP5p?2gT_m0^3rg6s=rPiu}J}PWQ(F!>vuB%CC}x#Y3AzQGnJ_4 zz4WJryP!$^_Kzl$)wePDEiNKJQn!_@o%fO7M4b?)^l`=kg=+Lo^vN+N%$=R%%M3R`#Qr?A5-pA5~JAxjKX%@=Rk?%3g6Flf0EZW1nwQ5ZpYL3gNDP`uq zG}ci?p-b|LPbC-#H$n_gm3{#5z{ove;ZEL0p+1;zF{2t~*R0ywqhFSIsX1o3I65L-1$n^8Q7lJ#iUT0b_<$1_ z`nl6@;HflRY8{Cm3v%;6lr`hv1dnYXLHa9iQ7}`?$j0j z_Aj_9XIwu7SZ5mmiZUNDv~hFe|KMTy>ueO>>H}QGemFMHx+wDi0cVSQN|AKmel@6G zS*u?<4IGdQgdmgD`wVX$<82F1d*T7)Q;uo_-&tGB;IVn{jB`jC9Dm*(D_<{j37I14 zdOgQ`YYZc>OnYFhw4bgN>yaZk>bY8Fz_KrET}CT}=SB0`MAW)W8OkspYs-eb^eLv zISU2OBYn~iRdGEGiP4)6g96Q2IzY5@ZuEt!ji04AXJOe?f2avCa&~r^Y`5|DU%9e z6SnOz2_f`~G(SUz6+D0<#1yYK`5K4aTI)&09s_1v4bb-{M!c+wPxX{804+A`!uX7y zJr*8X4&Y3naJ358ED+Ru|_AO(t8{GU8ag zH^WE1SH(-vCY}m?n7`uO=KnUVw;EL8bw}vijak#k;8Lr>kLRu~ zU;?wX(CdvsxM*>gDfUx`36$2RSRMdWEr9U8ji#T8`b-^x(Zhq-&42!5{TBV;YgE;T zHaTnzo>up+wi)RdoJ#hXt?>bM!z2%mt26LE&P~^YG6BZO@Bi%<%?BHby>#;PZpR$3 zi_oP82^_-9klITbeJ_aAD_k6WM_%eE7HUPhj8|bXnO->?4urV_OCaKvWVv_gP!}jc zR6(bm>2^e%=2UHkyru1-xz#bD%58~3o@S$)1 zfSAygksoXNA+~Nf_sQK7cy~i9LCen_b1p_!RX1#;djBF(N$Yw+B?~GHw|Rl2py{+0 z?$7bL9&3{eDKV@qJhW17!ZT3Z8v=}bHqXBY0COd0+I|{Hmhf2Z{ergPv z6R5_*-#)A-UsJ3-y5TO71cLa^IKOcIP(cLOnJoSU1F?0=*j}Q_*rAx>AdqxF4yF}! zuzV6EruepEJgkx4@tx_<_t#_J8y^5CzoS4EF%~Go;pjE-JlplQmU3Oxrh#fil_|Y3V9T+?6(Aq1cRqYIgvm?ZDI{}&> zhaNTc`Hpw>6~i@F3oDg0pFhMx>=J&}vJo_Xb(vVc%0u3$v0?Tr{b{oVLX!fJd3maz z{QL4@bVVWD+PP~HXdNKa`CYvP#Nqf?=7#$lo&|N^f!QH88|d&)csx%TudoFS_ZU1{ ztRL&P2>`^>=VL(1_4^+GVOn)KB6Xcs!Xv*?T&pB3FFekJA*Y;UHn)dgO5!n{3F6o> ztEM8kXWW|y)w%K&2w{DX1Uuefcyi1+0sNua^qwR?3ud(jtad>uL&!Gf0OJFOwu~I` zPVG2?H_VJtmF_^KMX6B&=*dj>Btb1rA0QnSmn?4}K>8?T{ZVpv0vE~m+h=QO@a~{> zk{2FY?`da}9e*as=h?-q8TLpS(B;hnq!_DaZu?vOZ3|-b$n-f}mmLL^IiJ)dw_cHH zu&i48<7*}6T@U!HcayKYPcI{Dd&ta=XQd7_gT9CM^HCag$f;qJ%Ij=3qSN6K^Md@- zLKYOotMqB#PmV1E5-w~MXyN^G!tnaNcB)gbl@d45vLrV5$%};A*gQH zVbOk#QG14mR;~3{UBoWY|FqbzSn7f9$$}54e*CKd2bcgV=LG4C>f%jHT^8&0- z;Os8ydO*0Bra~6|FzG1{nFQO`BB%8r3o=oM_IC2qL?iyQ&Hk&hbZ(Jm#>TOEc4yQZ zRcsyGQ&nN?{_w1SM(zGgBWNkI)hy7pW41xoMb5yTdKT*g9 z)<)|n)CbDP_9t#dgN~*R$YL+{>5p@)Gkak3cR(1eaj0tSblt3sPBR<125dblc zPLZnliFSnBgL^ER$%gC*2KNTE^DE3$#VvcTA?^|DJG|c9Uqr7|D0WEkGdR90w$P#{nNU5Q>@--qHtaeLygO>O zlzqOtpEOef1-#1vhq^fx`r*d?*uRkTketv9i78$Cn@GZ|iOkuyX9aZY9bt0E=o5h^wCw|{ z#AcW#kKKvUbQk0lCdidTCVb*A(% z=C0|_SA(VS6_m-PSad5`-wL*;cYjHp=mQxbkc3wT9$FSGpO;S|fh7Nb^^^SDABArnlT0mra4#|QNiDMpP7pgGd9 zaa!ZA#eD!>Zdg6a8W`I4+DfdjZ6DAaR2g&tok(VR>vfHDC8e@nFyt40ezl(H(yjU2 zqG=z}dh9Gk@kv1Ec+E%s?Y#o09dQxDx)y_rQKB4oB(I1jr5Ry@_IHxQSt;}t!xsUn zbP?Vl)$jb|pu<6FD8&93E8`h5+eL=Gi<92`cJiyse}Y68y$@u@``ghYKM46@k#xy< zH7BG#!{*o+9|U=hkJFEK!eO4Tt4heqe# zl70SH^p1FpX|4l5c)ag{?n@$!^IWN!dSKnlrZ~E}RzGRbl~j!I_=^N~Imy2&uCXQW zthU~lC%i0^Z842CnA`5^p&OEjb+gks6Qg(p zHOf9t8IECLUDw0B3%~)b$p4c6hwrAEDpdMZ7)W{j^eHOe`~8eq$PE#|#TrJqT0WX9 z>PQ9pIGgGxTL1bCoU49A@`}bR`yd}NP!V63{<$i1 z<1!F@4?kaNBqcO1^xJY7^DJ;0u~$d(P2NtnxqbH%BaicHLi3$sjMSsa3IY1CWD@hG z0sm+szD#3jxxCRn5Db(U&q;KDd7%9Bk(tLc$TU#h2mEi~rO0fpw@2>Ve<6|_KcZIY z!5tsUu8e+USf6oVEu-fPPWD0eUC_Y?=o4C@iz>P;7qRvJ7G}vQK+Ein*K@2HA)a{t zTp(x-rl!8-j~R29Eu>HGqubw|LMM1$C{pU7#e}XTJ4~;ZytLt0s-I(osmXr4QzG`s zn?LV@8p{hEl{m#k?j0NYufr$S}PyP?8TB?^I+Qhr*K6X2eF79%$7V}>0VUsvV3!q z0+5slOkN07tP55sytnMH63<#=m1p0SGKX#BXyqtOHN50xf8voIDn(k5>9kyeE_a5H zd(pna4D=nQKWzCY+Z_t;hB$ZHHxOcyIG2x3n_XyYlzgZJhs37py_b>fpQZl%yqtee z>#s%K2*dm+1!@A=gsUar*ZYRytw#BJC2y$TB!|1bgFx_GkYpWiP7WQ*a_C}CbQVTL zDznk9;E=vYJ(_+X^ei-!U6mvVKpS2MZLfPts*5~<*6tn<};^S_GLUmNqwl4ch@=rhF$kGe?SXlaa(}@-WMzv!E|{p6#rgE^?}Hs=!#g5^6i|5vEbypZfM$x!LNAgQ1K-yv z%}2g|kqM#as!-|WwUbENyRWwfhk|BjDn_s0Is^uQujrxpd|8RVk~|^F*yzFNLY2V+ zM{MylUznRjT9yXX0;Ab5j)4q;<{{@Hb6&tJXx29P0~)pLa$O;$$01Pi&pZ8Wq4_w+ z=?(lVn{rEkTR_P9v2i-xwz#eHUS!7S!*!x zUIteOVreXdCi&;rase0Q6=URZ2^0=ixR2Hzx`7g0z1z_ix(otqMn z7pN>(A1(#Som@dRt+GwcYyw)IQ(i^l2kqv4sqr45f23pu8NgMYkT+O3*O*mHvU^Yk zo0CD4Po0{`=hy1flwwmsNaCx~w1GdJ|6HFcwnAvwFkUW*`_vqxLL+|&@{9}QwirGV zpCZ!$R{o0^lav6Vjmq=F^-rL`+k5>b=&y{}iqg(zi}xxt{Tn|X{T1#A7n50+gs;lI{)v(jA6%_=j^y)ZV->!TpI*IR+i_XmmP z_XF%m)u9DrHUtc8K+PG=MK$==6Elc=V8gzN==1<;60HjKE^HO?t_EN2eNTA3x3t_w zGJ_dPVqv8xzhl>UIT5fGz$&igY4gGpy;I*%v$yKdmJK?jcrH_~z}oiQ{|Yaff!Qbp zike_x*dAdNHz$L=^e)Yy$yCD{WfILRp9AZA%q5LAhb+8?jMZ&%*9OH_4_F@6Q6*XZ zu4eBl@kR`Qmb)>#^hJ7mUCo?-@0gW;mmEg?~KgD7T`VY6Jx z@IuzccYIRz9Py7rt_^k!y?3l(h-TpGm!#?dSG5%^H64VPH{>kM zO*IrPHQ76UjID^Ka+g9U_(xeMeu<4=JZyK}v`=XbS+Xs7)beg*xFeN!Csv*g+ngib zdI;Dh;__#~`BYSzlF7ZHT(?G`+CdE*64A6Zg^k*(N3vH2_UJ_Ezy&)Tj_ub zge+i^)V4$><>4675(O(IKnyIlb^+*zctNB~rWs5-8Jp%d)MZ-5=X-q3pt;#!dr_pM zw#UC$uGV6QPUwZFsP%R=dq2u^*vcVlD^O;IwaT!o%AJN>cIWu}0`Hr#eh>d-&$ld% zC55b2(k`}%UkX`QWj&&kTSMskSN<$y6(;1*8PTLJHQi=4k!>mQ0E@Rom6cJ&n3XX) zY~~OHCeB-EIQ_M&$~b(Aqi{*COgr?jC1k8|n_I4eHE~d_g7x5g(bAmJZylFS$u1Mw z%4PB2hY7=GOQr|Xt7f5Xth2XISL}&+R~ z?7z*d5nL0;A;za{HwI=vYNdc@KFQ@e??nTjU|y6Zc-CcramIi__<#6CS)c=Fk6~as zhQNKbYDd6xd+7nwBxs=P$As%_f@`Hq_dlf_o{S}`|m}^)wmM2LK z7_>NY)P_!;9yQ8{eF1`(B*t~HI~8jA{S}3c2&#@s533xo2B5m^%&9nQ^qTVX310Ql zHq?5oQbg^M-3W~Qvte@DXx&`&O|7AqDQjKQgojgo%ld`Jua}Tq&;=qM-Osgo{v}9I zDvkAsZe63rykDUZ3VQ9QIP8HuI2(vR?ZpgONHxD7)j4zH7431^frcbYZ~-mQ{v`fT z*T%*ZvKg7%dTjd{|I47da*|GgfyG#5naB3h2M6cYO=54*#(`-r6SN(xe4HWIyv+ov zUznHffF8lHh^^LXg(@I&_W=gQ-uT>HKzPNWFQp(3JAAxfc;rK2n* z$|nru3bOfh9#z16F*8N9pgr z1HN#kS^*Xu!QT77ER{e=#dr{AgU6?@4~u`EnM7=51Br^?(TZo`CBFSONCuIgW?I-inXoM)$p?O3wRoS- zOO4qD?|`2vHf8`Iir$wn>yQ6&*|30u7kF&=v@REnUU6S&D5JS0x9XN93R=RZvp2k7A9>oL3sAmsW=DZE6+d( z*IxBu+jeKy_jh-CwH9RqR`fM9j@+>SG<3@`Y|D&7>;~(=$yX~>JRz_M$(kT1SK>PN zAwv%}vNj01eSPeGW03#2QapC2u=^ z5Rk7%pj=ZyRjq#TW&ImPOhz`m;Ud4&|74`p#%T}*cxc`xW@0x;;Lm4n*h&l`1ku(`{(fx(jkE=Pu~qOn+PklfwAuq3cVJ9< zxFCGw48VhGATW-LsQ5ONZ%^DGY5{{Zsf-zY`&YSStAFJwvIkpu0rU^Po*#0pvj^(Z zBcbtAbOr$MMBu`RyMy||D=?;AC4VIhix)DL^ek~+Xs z32-NO1K-YH#Y#=8GbMn(b4Q6#Fptq3&r3y53hvw#-3NEj*gEh$p%(68#wgv%a^QWM zoG&9M8+rfZ8 z*5Q&H)DyLpB;SWR#-$E#8z0;2!U_n%_XM@LzfngNxB1wX;cYE%9)k4KgW-?p-)7!k z?0Fk?s&HbH00j!CPY~Rr(UPCyH!f#qXdbPg(+sTI=xVpi_*H}8Y^)9HQqF+YFe4qS zHi7R!ZS*O(=@aT|o&zld+Wh&{8yW{w>nb*tKI1%ctW}`Ukgj^Z}Zb*F14Ja^tnT{B<17NAdUuxH&#l z_O8_%4DErPYWQOuaX&G5ZMBjX54uW1!6dG{C>z(SMTUu66AlXlj07w_}Eees)-%fF%(_$y)Z2nOahHa!;~@zi8et@vFBU?=`|hY`Hn;Izi+0*vTPFGBMBT zi5*EL3Ioum+h>z5fHs&%b|ltaC4zwg$Cx<*d~Hpw(m5F9eq+-FP6qMyO$D6GFtbQl zrPX@%PJG|H1GPtN64P}bf*|9a(XXaP5qd}Z*kXGXz6yeGf>W}q=&>7A5P#&KEQ}4Z zMEJrZsN?P-b`*oa#7l&@{okw(ZS+gRF5HjVs0CWI6KNTRS}VJhv8_d~hG$B4-e|sa z(c~lUe#Oa(yu>sEpDcd4bEy~gddkS^`M@4b^=2EUDR1iB*HG+l4&c@&)fClMP*_ga zV05)PLcatKFOc$9ygce4Fc@5hY4<+=v5~sHMeL&x(GW6wsnGuVHCeYxy+hhR2X?BC z71hw1s?F(t#hhY_-pw@BQejyWfTCy}cXQ;TDKd=2bv!N7x z_J~K8wD_UKRy{?_aMwg$e55-+b1mdhsnFTWY(HQbzwn4x>tm*pkE>5KXvY`ZMl=~1 znNnE&=+y1C!C1~aSwq?ak$q(z0}m0KgNTxs7ChJL{?WKH=7^W(tN!$9feYD|!IS(L zzG=QvJBX~gDkI%9Vb+uq<)Uw`6=apVnBZ|zR_)H*EGw#VzFxT)ma|d3w1%f5{O0FD<-X@?Eq_1xrK z$jDMwD~Zl=>5eJ5kCP-m9II z@g+RwlfP^^5oBtblD~&EL6%Dx=T;y>#>p_}uf?#YRjo!Kd?UEb%NDq6p3EaD;?!CP z)w10unKSCMPRnmWo{bAb7~kpgERfy))G8rnNRQ`VsC74;yWO}fpv_c)>UqozQ;$8e74B>?i-TN22n6; z8E}s+z%R}P=j!43g^|9z#Eh2nnbB7lO@0@-m%&}1_p@+R)O@L)d)s7Y_Y&E=RPG_s zFTNfBgh4W$pp^aZORWBXj?1w#-LXbv_6cg1Yga~)y!KrgI1EqBp;HZVbKOvve!3~M znRF`vs!ztZ_e4|E;XTSogJKTJ$*YIb9pv zdp_iSz8K7jc6@$@TF?@JG%=9B=VY&Gx5QVOiWO|#{lZ{((*>Wdiq7GQ+eud%%~yN{ zmqUzkzcMrq;w2&`_X+z$9pPYMO>5D)zhJfRm;YRe$NP=8?a6nNTqz9L+Cg(3GFYeB(F-4+pEWv?Wqb;F(a(ssxf4q2 zeUPlEL6YP{TV!?&`j|Vf?Xm&RQzsxft{@deC&SL>={MHF3N{8>13Cj+=n%l4Y04E# zdz2OAMKK#kmv}lrpxjGI{i-;X!jE=FM?fWoY|?Z$@HQ&|PN_%l&4LKUr)Oyt zL>PzU9!Tm6zYXiuWtp5#*Y}>}QNjD<2T6LtL`4YCfne4gh)d|qk1YFZo&y{;$YD&} zt19LO(mxlw9Xgd6xk^!%6|sc9?M00x$}TBn%%lkIK9!vwDiIxdm$;eWLp29NFK z&VdlsLVmlWIF)VB4eV1<$hZOgl0u*oq!1luibblnxTVe98~(D2>Z&_TH7N(JB%`+w z&1ef|9Dih|+M;={G;X+hnk|I+VPTu1B_Xea?wCK^EH(BD>Mw<-aDSOYUZ;`r9qK}- z6YEDxVG+Lp?YXNoLV6H#B^sCC%P1!v17M~b)l4hDHtG!Be8+HF3EP^+(wKmhNo)OR z6Vk@ZZ)MuN%kV;4sWLN4ZflTF;Th^N$Rg$MchhO@7dWLSd#5((ly`??Y3P725#Io8 z=nwx?&;p@BL~pGplCc*UPmq8bvt43$|4up_T60i@@N%WrnexVb+F7-RKgoPhnIGk9 zug5DgsO;n1#B#uhCNlMbdR7DPX zKsgndK%y(;2#RYLD6GSCgLPSuRoR?QU%R+34!Pw2{a*?c$dEmSZp1KB^iw`~yg9A( z5VWsnaFg-aQZ?Ib_te+kU3C3)2cNG=S|*7xEhRgNaNl<6g(Be_Ra)^2P@}q=6eP5$ zRgo3B@K(B3UgEQQcfDTiH5A3Mbqfc>` Date: Tue, 25 Feb 2025 00:57:42 -0500 Subject: [PATCH 10/16] Updated the DESCRIPTION and NEWS files to reflect the new package version --- DESCRIPTION | 8 ++++---- NEWS.md | 3 ++- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 02bca85..6efb7dd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,16 +1,16 @@ Package: learner Type: Package Title: Latent Space-Based Transfer Learning -Version: 0.2.0 +Version: 1.0.0 Authors@R: c(person("Sean", "McGrath", role = c("aut", "cre"), email = "sean.mcgrath514@gmail.com", comment = c(ORCID = "0000-0002-7281-3516")), + person("Ryan", "ODea", role = c("aut"), + email = "ryanodea@hsph.harvard.edu") person("Cenhao", "Zhu", role = c("aut"), email = "cenhaoz@mit.edu"), person("Rui", "Duan", role = c("aut"), - email = "rduan@hsph.harvard.edu"), - person("Ryan", "ODea", role = c("ctb"), - email = "ryanodea@hsph.harvard.edu")) + email = "rduan@hsph.harvard.edu")) Maintainer: Sean McGrath Description: Implements transfer learning methods for low-rank matrix estimation. These methods leverage similarity in the latent row and column diff --git a/NEWS.md b/NEWS.md index 1eb9ec0..84fe99e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,6 @@ -### learner version 0.2.0 (TBD) +### learner version 1.0.0 (TBD) +* Implemented the optimization algorithm in Rcpp, resulting in significant computational speed-ups. * Fixed a bug in `cv.learner` when `Y_target` has missing entries. ### learner version 0.1.0 (2025-01-08) From 8bfb463b4e05be1ca851c3bf85e03926e9f133ef Mon Sep 17 00:00:00 2001 From: Sean McGrath Date: Wed, 26 Feb 2025 01:07:57 -0500 Subject: [PATCH 11/16] Included percent non-missing and added corresponding tests --- DESCRIPTION | 2 +- src/kernels.cpp | 3 +-- tests/testthat/test-main.R | 26 ++++++++++++++++++++++++++ 3 files changed, 28 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 6efb7dd..b9ea3e5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -6,7 +6,7 @@ Authors@R: c(person("Sean", "McGrath", role = c("aut", "cre"), email = "sean.mcgrath514@gmail.com", comment = c(ORCID = "0000-0002-7281-3516")), person("Ryan", "ODea", role = c("aut"), - email = "ryanodea@hsph.harvard.edu") + email = "ryanodea@hsph.harvard.edu"), person("Cenhao", "Zhu", role = c("aut"), email = "cenhaoz@mit.edu"), person("Rui", "Duan", role = c("aut"), diff --git a/src/kernels.cpp b/src/kernels.cpp index 77d011f..2a706ef 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -40,8 +40,7 @@ List learner_worker(const Eigen::MatrixXd &Y_source, Eigen::MatrixXd V_trunc_T = V_trunc.transpose(); // For now, we assume no missing data (perc_nonmissing = 1). - // double perc_nonmissing = (Y_target.array().isNaN().count() > 0) ? ( (p*q - num_missing) / double(p*q)) : 1.0; - double perc_nonmissing = 1.0; + double perc_nonmissing = 1.0 - (static_cast((Y_target.array().isNaN()).count()) / Y_target.size()); bool missing = Y_target.hasNaN(); double obj_init = 0.0; diff --git a/tests/testthat/test-main.R b/tests/testthat/test-main.R index 2de022f..a6f31f1 100644 --- a/tests/testthat/test-main.R +++ b/tests/testthat/test-main.R @@ -25,6 +25,19 @@ test_that("learner does not fail", { step_size = 0.003)) }) +test_that("learner output is correct", { + expected_result <- t(matrix(c(0.1578405, -1.54883716, 1.2205437, 0.03377010, 0.0767873, + 0.1189457, -0.12793612, -0.9755801, 0.01372606, -0.3902793, + -0.4384613, 1.21518704, -1.0718671, 0.60388502, -0.4104101, + 1.4426139, -3.13260948, 0.9175863, -1.79041482, 0.5015381, + 0.5602868, -0.05068081, -1.3878535, -0.78502749, -0.1315952), nrow = 5, ncol = 5)) + result <- learner(Y_source = dat_highsim$Y_source, + Y_target = dat_highsim$Y_target, + lambda_1 = 1, lambda_2 = 1, + step_size = 0.003) + expect_equal(result$learner_estimate[1:5, 1:5], expected_result, tolerance = 1e-5) +}) + test_that("dlearner does not fail", { expect_no_error( dlearner(Y_source = dat_highsim$Y_source, @@ -45,6 +58,19 @@ test_that("learner does not fail with missing data", { step_size = 0.003)) }) +test_that("learner output is correct with missing data", { + expected_result <- t(matrix(c(0.1402368, -1.53788299, 1.294196, 0.06241446, 0.1034029, + 0.2234374, -0.16286860, -1.017956, -0.13353267, -0.3327750, + -0.4183411, 1.20040593, -1.115950, 0.59002710, -0.4299424, + 1.4181535, -3.09632714, 1.011979, -1.75238013, 0.5554817, + 0.5603679, -0.06539629, -1.381598, -0.75764518, -0.1373639), nrow = 5, ncol = 5)) + result <- learner(Y_source = dat_highsim$Y_source, + Y_target = Y_target_highsim_missing, + lambda_1 = 1, lambda_2 = 1, + step_size = 0.003) + expect_equal(result$learner_estimate[1:5, 1:5], expected_result, tolerance = 1e-5) +}) + test_that("cv.learner does not fail with missing data", { expect_no_error( cv.learner(Y_source = dat_highsim$Y_source, From d6bf8a06dc6a1a3254a6c56f280e6590f53ad042 Mon Sep 17 00:00:00 2001 From: Sean McGrath Date: Thu, 27 Feb 2025 18:34:08 -0500 Subject: [PATCH 12/16] Revised fold index sets to respect the missing values in Y_target --- R/RcppExports.R | 4 ++-- R/learner.R | 17 ++++++++++++++++- src/RcppExports.cpp | 9 +++++---- src/kernels.cpp | 22 ++++++---------------- tests/testthat/test-main.R | 26 ++++++++++++++++++++++++++ 5 files changed, 55 insertions(+), 23 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index e267889..ec831f9 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -5,8 +5,8 @@ learner_cpp <- function(Y_source, Y_target, r, lambda1, lambda2, step_size, max_ .Call(`_learner_learner_cpp`, Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold, max_value) } -cv_learner_cpp <- function(Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r, max_value) { - .Call(`_learner_cv_learner_cpp`, Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r, max_value) +cv_learner_cpp <- function(Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r, max_value, index_set) { + .Call(`_learner_cv_learner_cpp`, Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r, max_value, index_set) } omp_max_threads <- function() { diff --git a/R/learner.R b/R/learner.R index 362fa08..5116111 100644 --- a/R/learner.R +++ b/R/learner.R @@ -59,8 +59,23 @@ cv.learner <- function(Y_source, Y_target, r, lambda_1_all, lambda_2_all, control$max_value <- 10 } + if (n_folds < 2){ + stop('n_folds must be 2 or greater.') + } + + ## Creating training and testing data sets + available_indices <- which(!is.na(Y_target)) + n_indices <- length(available_indices) + indices <- sample(available_indices, size = n_indices, replace = FALSE) - 1 + + index_set <- vector(mode = "list", length = n_folds) + for (fold in 1:n_folds){ + index_set[[fold]] <- indices[floor((fold - 1) * n_indices / n_folds + 1):floor(fold * n_indices / n_folds)] + } + result <- cv_learner_cpp(Y_source, Y_target, lambda_1_all, lambda_2_all, step_size, - n_folds, control$max_iter, control$threshold, n_cores, r, control$max_value) + n_folds, control$max_iter, control$threshold, n_cores, r, control$max_value, + index_set) return(result) } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 43de911..f39d9a3 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -31,8 +31,8 @@ BEGIN_RCPP END_RCPP } // cv_learner_cpp -List cv_learner_cpp(const Eigen::MatrixXd& Y_source, const Eigen::MatrixXd& Y_target, const std::vector& lambda1_all, const std::vector& lambda2_all, double step_size, int n_folds, int max_iter, double threshold, int n_cores, int r, double max_value); -RcppExport SEXP _learner_cv_learner_cpp(SEXP Y_sourceSEXP, SEXP Y_targetSEXP, SEXP lambda1_allSEXP, SEXP lambda2_allSEXP, SEXP step_sizeSEXP, SEXP n_foldsSEXP, SEXP max_iterSEXP, SEXP thresholdSEXP, SEXP n_coresSEXP, SEXP rSEXP, SEXP max_valueSEXP) { +List cv_learner_cpp(const Eigen::MatrixXd& Y_source, const Eigen::MatrixXd& Y_target, const std::vector& lambda1_all, const std::vector& lambda2_all, double step_size, int n_folds, int max_iter, double threshold, int n_cores, int r, double max_value, const std::vector>& index_set); +RcppExport SEXP _learner_cv_learner_cpp(SEXP Y_sourceSEXP, SEXP Y_targetSEXP, SEXP lambda1_allSEXP, SEXP lambda2_allSEXP, SEXP step_sizeSEXP, SEXP n_foldsSEXP, SEXP max_iterSEXP, SEXP thresholdSEXP, SEXP n_coresSEXP, SEXP rSEXP, SEXP max_valueSEXP, SEXP index_setSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -47,7 +47,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< int >::type n_cores(n_coresSEXP); Rcpp::traits::input_parameter< int >::type r(rSEXP); Rcpp::traits::input_parameter< double >::type max_value(max_valueSEXP); - rcpp_result_gen = Rcpp::wrap(cv_learner_cpp(Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r, max_value)); + Rcpp::traits::input_parameter< const std::vector>& >::type index_set(index_setSEXP); + rcpp_result_gen = Rcpp::wrap(cv_learner_cpp(Y_source, Y_target, lambda1_all, lambda2_all, step_size, n_folds, max_iter, threshold, n_cores, r, max_value, index_set)); return rcpp_result_gen; END_RCPP } @@ -64,7 +65,7 @@ END_RCPP static const R_CallMethodDef CallEntries[] = { {"_learner_learner_cpp", (DL_FUNC) &_learner_learner_cpp, 9}, - {"_learner_cv_learner_cpp", (DL_FUNC) &_learner_cv_learner_cpp, 11}, + {"_learner_cv_learner_cpp", (DL_FUNC) &_learner_cv_learner_cpp, 12}, {"_learner_omp_max_threads", (DL_FUNC) &_learner_omp_max_threads, 0}, {NULL, NULL, 0} }; diff --git a/src/kernels.cpp b/src/kernels.cpp index 2a706ef..7358469 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -146,23 +146,13 @@ List learner_cpp(const Eigen::MatrixXd &Y_source, const Eigen::MatrixXd &Y_targe List cv_learner_cpp(const Eigen::MatrixXd &Y_source, const Eigen::MatrixXd &Y_target, const std::vector &lambda1_all, const std::vector &lambda2_all, double step_size, int n_folds, int max_iter, double threshold, - int n_cores, int r, double max_value) { + int n_cores, int r, double max_value, + const std::vector> &index_set) { int p = Y_source.rows(); int q = Y_source.cols(); int n_lambda1 = lambda1_all.size(); int n_lambda2 = lambda2_all.size(); - std::vector indices(p * q); - std::iota(indices.begin(), indices.end(), 0); - std::shuffle(indices.begin(), indices.end(), std::mt19937{std::random_device{}()}); - - std::vector> index_set(n_folds); - for (int fold = 0; fold < n_folds; ++fold) { - int start = fold * (indices.size() / n_folds); - int end = (fold + 1) * (indices.size() / n_folds); - index_set[fold] = std::vector(indices.begin() + start, indices.begin() + end); - } - Eigen::MatrixXd mse_all = Eigen::MatrixXd::Zero(n_lambda1, n_lambda2); #pragma omp parallel for collapse(2) schedule(dynamic) num_threads(n_cores) @@ -173,8 +163,8 @@ List cv_learner_cpp(const Eigen::MatrixXd &Y_source, const Eigen::MatrixXd &Y_ta for (int fold = 0; fold < n_folds; ++fold) { Eigen::MatrixXd Y_train = Y_target; for (int idx : index_set[fold]) { - int row = idx / q; - int col = idx % q; + int row = idx % p; + int col = idx / p; Y_train(row, col) = std::numeric_limits::quiet_NaN(); } List temp = learner_worker(Y_source, Y_train, r, @@ -183,8 +173,8 @@ List cv_learner_cpp(const Eigen::MatrixXd &Y_source, const Eigen::MatrixXd &Y_ta Eigen::MatrixXd learner_estimate = temp["learner_estimate"]; double fold_mse = 0.0; for (int idx : index_set[fold]) { - int row = idx / q; - int col = idx % q; + int row = idx % p; + int col = idx / p; double diff = learner_estimate(row, col) - Y_target(row, col); fold_mse += diff * diff; } diff --git a/tests/testthat/test-main.R b/tests/testthat/test-main.R index a6f31f1..b1b14dd 100644 --- a/tests/testthat/test-main.R +++ b/tests/testthat/test-main.R @@ -17,6 +17,19 @@ test_that("cv.learner (parallel) does not fail", { n_cores = 2)) }) +test_that("cv.learner output is correct", { + expected_result <- t(matrix(c(5494.916, 5312.724, + 5341.662, 5233.828, + 5009.216, 5013.827), nrow = 2, ncol = 3)) + set.seed(1234) + result <- cv.learner(Y_source = dat_highsim$Y_source, + Y_target = dat_highsim$Y_target, + lambda_1_all = c(1, 10, 100), + lambda_2_all = c(1, 10), + step_size = 0.003) + expect_equal(result$mse, expected_result, tolerance = 1e-3) +}) + test_that("learner does not fail", { expect_no_error( learner(Y_source = dat_highsim$Y_source, @@ -79,3 +92,16 @@ test_that("cv.learner does not fail with missing data", { lambda_2_all = c(1, 10), step_size = 0.003)) }) + +test_that("cv.learner output is correct with missing data", { + expected_result <- t(matrix(c(5457.513, 5254.740, + 5291.389, 5175.707, + 4933.222, 4947.356), nrow = 2, ncol = 3)) + set.seed(1234) + result <- cv.learner(Y_source = dat_highsim$Y_source, + Y_target = Y_target_highsim_missing, + lambda_1_all = c(1, 10, 100), + lambda_2_all = c(1, 10), + step_size = 0.003) + expect_equal(result$mse, expected_result, tolerance = 1e-3) +}) From 0ab83886201e2be3161b7a4c92f3b33ba5b7ac48 Mon Sep 17 00:00:00 2001 From: Sean McGrath Date: Thu, 27 Feb 2025 18:54:41 -0500 Subject: [PATCH 13/16] Revised order of computations in dlearner to improve speed --- R/dlearner.R | 5 +++-- tests/testthat/test-main.R | 11 +++++++++++ 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/R/dlearner.R b/R/dlearner.R index 1885b82..0454de2 100644 --- a/R/dlearner.R +++ b/R/dlearner.R @@ -54,8 +54,9 @@ dlearner <- function(Y_source, Y_target, r){ } svd_source <- svd(Y_source, nu = r, nv = r) - dlearner_estimate <- (svd_source$u %*% t(svd_source$u)) %*% - Y_target %*% (svd_source$v %*% t(svd_source$v)) + dlearner_estimate <- svd_source$u %*% + (t(svd_source$u) %*% Y_target %*% svd_source$v) %*% + t(svd_source$v) colnames(dlearner_estimate) <- colnames(Y_source) rownames(dlearner_estimate) <- rownames(Y_source) diff --git a/tests/testthat/test-main.R b/tests/testthat/test-main.R index b1b14dd..7cf7bcb 100644 --- a/tests/testthat/test-main.R +++ b/tests/testthat/test-main.R @@ -57,6 +57,17 @@ test_that("dlearner does not fail", { Y_target = dat_highsim$Y_target)) }) +test_that("dlearner output is correct", { + expected_result <- t(matrix(c(0.0959171, -1.72143637, 1.15500149, 0.005478273, 0.3258085, + 0.1077137, -0.32243480, -1.03557177, -0.108195636, -0.4093049, + -0.7237275, 0.86634922, -0.43528206, 1.105802194, -0.5961351, + 1.6347696, -2.91985925, 0.04876288, -2.357504089, 0.8127403, + 0.4757450, 0.05665543, -1.50373470, -0.744076401, -0.2876033), nrow = 5, ncol = 5)) + result <- dlearner(Y_source = dat_highsim$Y_source, + Y_target = dat_highsim$Y_target) + expect_equal(result$dlearner_estimate[1:5, 1:5], expected_result, tolerance = 1e-5) +}) + # Missing data cases set.seed(1234) From 02a5a041b68194456eaac1f068786823ac89a1d9 Mon Sep 17 00:00:00 2001 From: ryan-odea <70209371+ryan-odea@users.noreply.github.com> Date: Fri, 28 Feb 2025 11:43:06 -0500 Subject: [PATCH 14/16] Clean up completed comments --- src/kernels.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/kernels.cpp b/src/kernels.cpp index 2a706ef..0daaf01 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -39,7 +39,6 @@ List learner_worker(const Eigen::MatrixXd &Y_source, Eigen::MatrixXd U_trunc_T = U_trunc.transpose(); Eigen::MatrixXd V_trunc_T = V_trunc.transpose(); - // For now, we assume no missing data (perc_nonmissing = 1). double perc_nonmissing = 1.0 - (static_cast((Y_target.array().isNaN()).count()) / Y_target.size()); bool missing = Y_target.hasNaN(); @@ -108,7 +107,7 @@ List learner_worker(const Eigen::MatrixXd &Y_source, V_best = V; } - // Simple convergence check (you might refine this). + // Checking for convergence conditions if (iter > 0 && std::abs(obj - obj_values[iter - 1]) < threshold) { convergence_criterion = 1; break; From 62a7e97963f0287b83c0eef498a2e8be11b73b48 Mon Sep 17 00:00:00 2001 From: Sean McGrath Date: Fri, 28 Feb 2025 11:56:46 -0500 Subject: [PATCH 15/16] Revised the onAttach message --- R/globals.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/globals.R b/R/globals.R index 2a56d31..7012eb0 100644 --- a/R/globals.R +++ b/R/globals.R @@ -8,7 +8,7 @@ if(getRversion() >= "2.15.1"){ if (max_threads < 2) { packageStartupMessage("*******\n This installation of learner has not detected OpenMP support\n - It will still work but will not support multithreading via the `ncores` argument + It will still work but will not support multithreading via the `n_cores` argument If you plan to use multithreading, please ensure you are using R>=3.4.0 and have OpenMP installed\n *******") } From 5c99c849b09c5f07b43949da0832f2c10bdecb3f Mon Sep 17 00:00:00 2001 From: Sean McGrath Date: Fri, 28 Feb 2025 12:03:15 -0500 Subject: [PATCH 16/16] Cleanup learner_worker --- src/kernels.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/kernels.cpp b/src/kernels.cpp index 7358469..2c99b4a 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -19,8 +19,6 @@ List learner_worker(const Eigen::MatrixXd &Y_source, int r, double lambda1, double lambda2, double step_size, int max_iter, double threshold, double max_value) { - int p = Y_source.rows(); - int q = Y_source.cols(); Eigen::BDCSVD svd(Y_source, Eigen::ComputeThinU | Eigen::ComputeThinV); Eigen::MatrixXd U_full = svd.matrixU();