diff --git a/.Rbuildignore b/.Rbuildignore index d96d2f7..42bec2f 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -4,3 +4,4 @@ ^\.github$ ^data-raw$ ^README_files/ +^\.Rprofile$ diff --git a/.gitignore b/.gitignore index 5b6a065..7234786 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,6 @@ .Rhistory .RData .Ruserdata + +.Rprofile +.DS_Store diff --git a/DESCRIPTION b/DESCRIPTION index ec868aa..d7b4663 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,10 +1,12 @@ Package: learner Type: Package Title: Latent Space-Based Transfer Learning -Version: 1.0.1 +Version: 1.1.0 Authors@R: c(person("Sean", "McGrath", role = c("aut", "cre"), email = "sean.mcgrath514@gmail.com", comment = c(ORCID = "0000-0002-7281-3516")), + person("Jodie", "Li", role = "aut", + email = "ljuanshuo@gmail.com"), person("Ryan", "ODea", role = c("aut"), email = "ryanodea@hsph.harvard.edu"), person("Cenhao", "Zhu", role = c("aut"), @@ -14,14 +16,16 @@ Authors@R: c(person("Sean", "McGrath", role = c("aut", "cre"), Maintainer: Sean McGrath Description: Implements transfer learning methods for low-rank matrix estimation. These methods leverage similarity in the latent row and column - spaces between the source and target populations to improve estimation in - the target population. The methods include the LatEnt spAce-based tRaNsfer + spaces between one or more source populations and a target population to + improve estimation in the target population. The methods include the + LatEnt spAce-based tRaNsfer lEaRning (LEARNER) method and the direct projection LEARNER (D-LEARNER) method described by McGrath et al. (2024) . + Supports weighted combinations of source spaces, separate row- and + column-space penalties, and cross-validation for penalty selection. License: GPL (>=3) Encoding: UTF-8 LazyData: true -RoxygenNote: 7.3.2 URL: https://github.com/stmcg/learner BugReports: https://github.com/stmcg/learner/issues Imports: ScreeNOT, @@ -33,3 +37,4 @@ Suggests: Config/testthat/edition: 3 Depends: R (>= 2.10) +RoxygenNote: 7.3.2 diff --git a/NEWS.md b/NEWS.md index a538a8b..e9ac59c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,12 @@ +### learner version 1.1.0 (unreleased) + +* Add support for multiple source populations and optional source weights in + `learner()`, `dlearner()`, and `cv.learner()`. +* Add `dat_multisource`, an example dataset with one target and two source + populations with independently generated noise. +* Add separate row- and column-space penalties in `learner()` and + cross-validation over these penalties in `cv.learner()`. + ### learner version 1.0.0 (2025-03-02) * Re-implemented most of `learner` and `cv.learner` in Rcpp, resulting in significant speed-ups diff --git a/R/RcppExports.R b/R/RcppExports.R index ec831f9..d1b4b70 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, max_value) { - .Call(`_learner_learner_cpp`, Y_source, Y_target, r, lambda1, lambda2, step_size, max_iter, threshold, max_value) +learner_cpp <- function(Y_source, Y_target, source_weights, r, lambda1_row, lambda1_col, lambda2, step_size, max_iter, threshold, max_value) { + .Call(`_learner_learner_cpp`, Y_source, Y_target, source_weights, r, lambda1_row, lambda1_col, 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, 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) +cv_learner_cpp <- function(Y_source, Y_target, source_weights, row_grid, col_grid, balance_grid, step_size, max_iter, threshold, n_cores, r, max_value, index_set) { + .Call(`_learner_cv_learner_cpp`, Y_source, Y_target, source_weights, row_grid, col_grid, balance_grid, step_size, max_iter, threshold, n_cores, r, max_value, index_set) } omp_max_threads <- function() { diff --git a/R/data-multisource.R b/R/data-multisource.R new file mode 100644 index 0000000..d98d4d9 --- /dev/null +++ b/R/data-multisource.R @@ -0,0 +1,43 @@ +#' Simulated data set: Multiple sources with different noise levels +#' +#' Simulated matrices for one target population and two source populations. +#' The sources share the target's latent row and column spaces, but source A +#' is observed more precisely than source B. This data set illustrates fitting +#' LEARNER and D-LEARNER with multiple sources. +#' +#' @docType data +#' @format A list with four components. All matrices have 100 rows and 50 columns, +#' with corresponding rows and columns aligned across populations. +#' \describe{ +#' \item{\code{Y_source}}{A named list of observed source matrices, \code{A} +#' and \code{B}, suitable for the \code{Y_source} argument.} +#' \item{\code{Y_target}}{The observed target matrix.} +#' \item{\code{Theta_source}}{A named list of the true rank-3 source matrices, +#' \code{A} and \code{B}. The two true matrices are identical.} +#' \item{\code{Theta_target}}{The true rank-3 target matrix.} +#' } +#' @details +#' The source signal is obtained by reversing the three nonzero singular values +#' of the target signal while keeping its singular vectors. Each observed matrix +#' is generated by adding independent, mean-zero normal noise to its own signal. +#' Noise is independent across entries and populations, with variances 0.10 for +#' the target, 0.01 for source A, and 0.50 for source B. Thus, source A has less +#' noise than the target and source B has more. Source B is not generated from +#' the observed source A matrix. +#' +#' The true signal matrices are included for illustration and estimation-error +#' calculations; they are not inputs to the fitting functions. These simulated +#' data are a usage example, not a reproduction of the paper's simulation study. +#' +#' @examples +#' data(dat_multisource) +#' fit <- dlearner( +#' Y_source = dat_multisource$Y_source, +#' Y_target = dat_multisource$Y_target, +#' r = 3, source_weights = c(0.7, 0.3) +#' ) +#' dim(fit$dlearner_estimate) +#' @seealso \code{\link{dat_highsim}}, \code{\link{dat_modsim}}, +#' \code{\link{learner}}, \code{\link{dlearner}} +#' @keywords datasets +"dat_multisource" diff --git a/R/dlearner.R b/R/dlearner.R index 0454de2..014d9ad 100644 --- a/R/dlearner.R +++ b/R/dlearner.R @@ -1,66 +1,106 @@ -#' Latent space-based transfer learning +#' Direct projection transfer learning #' -#' This function applies the Direct project LatEnt spAce-based tRaNsfer lEaRning (D-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. +#' @description +#' This function applies the Direct projection LatEnt spAce-based tRaNsfer +#' lEaRning (D-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. -#' -#' @return A list with the following components: -#' \item{dlearner_estimate}{matrix containing the D-LEARNER estimate of the target population knowledge graph.} -#' \item{r}{rank value used.} +#' @param Y_source Matrix containing the source population data, with the same rows +#' and columns as \code{Y_target}. For multiple sources, supply a list of +#' matrices. +#' @param Y_target Matrix containing the target population data. Missing entries +#' are not allowed. +#' @param r Rank of the knowledge graphs. If omitted, ScreeNOT (Donoho et al. 2023) +#' is applied to the source population data to select the rank. +#' @param source_weights Relative weights of the sources, in the order of +#' \code{Y_source}. Default is equal weights. #' #' @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}. +#' 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 rank-\eqn{r} 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}. +#' 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{\hat{U}_{1}^{\top}\hat{U}_{1} Y_0 \hat{V}_{1}^{\top}\hat{V}_{1}}. +#' This method estimates \eqn{\Theta_0} by +#' +#' \deqn{\hat U_1\hat U_1^\top Y_0 \hat V_1\hat V_1^\top.} +#' +#' That is, the target population data are projected onto the latent row and +#' column spaces of the source population, without iterative fitting or penalty +#' selection. +#' +#' \strong{Multiple sources:} +#' +#' Supply a list of source matrices to combine information from several +#' populations. Their projection operators are averaged using +#' \code{source_weights}, and the estimate becomes +#' +#' \deqn{\widehat\Theta_0=\bar P_U Y_0\bar P_V,\qquad +#' \bar P_U=\sum_{k=1}^K w_k\hat U_k\hat U_k^\top,\qquad +#' \bar P_V=\sum_{k=1}^K w_k\hat V_k\hat V_k^\top.} +#' +#' where \eqn{\hat U_k} and \eqn{\hat V_k} come from the truncated SVD of +#' the \eqn{k}th source and the weights \eqn{w_k} are given by +#' \code{source_weights}. If \code{r} is not specified, the rank is selected +#' using the first source. +#' +#' @return +#' A list with the following components: +#' \item{dlearner_estimate}{matrix containing the D-LEARNER estimate of the target population knowledge graph.} +#' \item{r}{rank value used.} #' #' @references #' 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. #' +#' @seealso +#' \code{\link{learner}} for a target fit that can adapt the source +#' spaces through penalized estimation. +#' #' @examples #' res <- dlearner(Y_source = dat_highsim$Y_source, #' Y_target = dat_highsim$Y_target) #' +#' # Multiple sources with fixed demonstration weights. +#' res_multi <- dlearner( +#' Y_source = dat_multisource$Y_source, +#' Y_target = dat_multisource$Y_target, r = 3, +#' source_weights = c(0.7, 0.3)) +#' dim(res_multi$dlearner_estimate) +#' #' @export -dlearner <- function(Y_source, Y_target, 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 (any(is.na(Y_target))){ - stop('Y_target cannot have NA values') - } - - 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) +dlearner <- function(Y_source, Y_target, r, source_weights = NULL) { + sources <- prepare_sources(Y_source, Y_target, source_weights, + allow_missing_target = FALSE) + r <- resolve_rank(if (missing(r)) NULL else r, sources$matrices[[1]]) + bases <- lapply(sources$matrices, svd, nu = r, nv = r) + if (length(bases) == 1L) { + # Preserve the original single-source multiplication order. + s <- bases[[1]] + estimate <- s$u %*% (t(s$u) %*% Y_target %*% s$v) %*% t(s$v) + } else { + # Average source projectors, not raw source matrices. + middle <- matrix(0, nrow(Y_target), ncol(Y_target)) + for (k in seq_along(bases)) { + u <- bases[[k]]$u + middle <- middle + sources$weights[k] * (u %*% crossprod(u, Y_target)) + } + estimate <- matrix(0, nrow(Y_target), ncol(Y_target)) + for (k in seq_along(bases)) { + v <- bases[[k]]$v + estimate <- estimate + sources$weights[k] * ((middle %*% v) %*% t(v)) + } } - - 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) - - colnames(dlearner_estimate) <- colnames(Y_source) - rownames(dlearner_estimate) <- rownames(Y_source) - - return(list(dlearner_estimate = dlearner_estimate, r = r)) + dimnames(estimate) <- dimnames(sources$matrices[[1]]) + list(dlearner_estimate = estimate, r = r) } - diff --git a/R/input-validation.R b/R/input-validation.R new file mode 100644 index 0000000..d942698 --- /dev/null +++ b/R/input-validation.R @@ -0,0 +1,79 @@ +# Normalize one source or a source list without changing the user's ordering. +# The first source determines the default rank and LEARNER initialization. +prepare_sources <- function(Y_source, Y_target, source_weights = NULL, + allow_missing_target = TRUE) { + if (!is.matrix(Y_target) || !is.numeric(Y_target) || + any(dim(Y_target) == 0L) || any(is.infinite(Y_target))) { + stop('Y_target must be a nonempty numeric matrix without infinite values.') + } + if (all(is.na(Y_target))) { + stop('Y_target must contain observed values.') + } + if (!allow_missing_target && anyNA(Y_target)) { + stop('dlearner requires a complete Y_target without missing values.') + } + sources <- if (is.matrix(Y_source)) list(Y_source) else Y_source + if (!is.list(sources) || length(sources) == 0L) { + stop('Y_source must be a matrix or a nonempty list of matrices.') + } + for (source in sources) { + if (!is.matrix(source) || !is.numeric(source) || any(!is.finite(source))) { + stop('Every source must be a numeric matrix with finite values and no NA values.') + } + if (!identical(dim(source), dim(Y_target))) { + stop('All source matrices and Y_target must have the same dimensions.') + } + } + if (is.null(source_weights)) source_weights <- rep(1, length(sources)) + if (!is.numeric(source_weights) || length(source_weights) != length(sources) || + any(!is.finite(source_weights)) || any(source_weights < 0) || + !any(source_weights > 0)) { + stop('source_weights must have one finite nonnegative value per source and a positive total.') + } + weights <- source_weights / max(source_weights) + list(matrices = sources, weights = weights / sum(weights)) +} + +resolve_rank <- function(r, source) { + if (is.null(r)) { + r <- max(ScreeNOT::adaptiveHardThresholding( + Y = source, k = min(dim(source)) / 3)$r, 1) + } + if (!is.numeric(r) || length(r) != 1L || !is.finite(r) || + r < 1 || r != floor(r) || r > min(dim(source))) { + stop('r must be an integer between 1 and the smaller matrix dimension.') + } + as.integer(r) +} + +validate_penalty_grid <- function(values, name) { + if (!is.numeric(values) || length(values) == 0L || + any(!is.finite(values)) || any(values < 0)) { + stop(name, ' must be a nonempty vector of finite nonnegative numbers.', call. = FALSE) + } +} + +validate_count <- function(value, name, lower = 1, upper = .Machine$integer.max) { + if (!is.numeric(value) || length(value) != 1L || !is.finite(value) || + value != floor(value) || value < lower || value > upper) { + stop(name, ' must be an integer between ', lower, ' and ', upper, '.', call. = FALSE) + } +} + +prepare_control <- function(control, step_size) { + if (!is.list(control)) stop('control must be a list.') + if (!is.numeric(step_size) || length(step_size) != 1L || + !is.finite(step_size) || step_size <= 0) { + stop('step_size must be a finite positive number.') + } + if (is.null(control$max_iter)) control$max_iter <- 100L + if (is.null(control$threshold)) control$threshold <- 0.001 + if (is.null(control$max_value)) control$max_value <- 10 + validate_count(control$max_iter, 'control$max_iter') + validate_space_penalty(control$threshold, 'control$threshold') + if (!is.numeric(control$max_value) || length(control$max_value) != 1L || + !is.finite(control$max_value) || control$max_value <= 0) { + stop('control$max_value must be a finite positive number.') + } + control +} diff --git a/R/learner.R b/R/learner.R index 5116111..00fc67d 100644 --- a/R/learner.R +++ b/R/learner.R @@ -1,29 +1,71 @@ #' 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 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 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.} +#' @description +#' Performs k-fold cross-validation to select the source-space and balance +#' penalties for \code{\link{learner}}. Returns prediction errors and the +#' penalty combination with the smallest error. +#' +#' @param Y_source Source population data, as in \code{\link{learner}}. +#' Supply a matrix for one source or a list of matrices for multiple sources. +#' @param Y_target Target population data matrix, as in \code{\link{learner}}. +#' Only observed entries are used to evaluate predictions. +#' @param r Rank used in every fold, as in \code{\link{learner}}. It is not +#' selected by cross-validation. +#' @param lambda_1_all Candidate values for the common source-space penalty +#' \eqn{\lambda_1}. All combinations with \code{lambda_2_all} are evaluated. +#' @param lambda_2_all Candidate values for the balance penalty \eqn{\lambda_2}. +#' @param step_size Step size for each fit, as in \code{\link{learner}}. +#' @param n_folds Number of cross-validation folds. Default is \code{4}. +#' @param n_cores Number of CPU cores used to evaluate the candidate penalty +#' combinations in parallel. Default is \code{1}, i.e., no parallelization. +#' @param control Optimization settings for each fit; see \code{\link{learner}}. +#' @param lambda_1_row_all,lambda_1_col_all Candidate values for separate +#' penalties on the row and column factors. Supply both instead of +#' \code{lambda_1_all}; all combinations with \code{lambda_2_all} are evaluated. +#' @param source_weights Relative source weights, as in \code{\link{learner}}. The weights +#' are not selected by cross-validation. Default is equal weights. #' #' @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. +#' 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 comparing +#' the entries in the test data sets with those imputed from the LEARNER +#' estimates. See McGrath et al. (2024) for further details. +#' +#' With separate penalties, the procedure selects the row, column, and balance +#' penalties together. Both modes can be used with multiple sources. Call +#' \code{set.seed()} first to reproduce the random partition. +#' +#' To obtain the final estimate, apply \code{\link{learner}} with the selected +#' penalties to the full data (see Examples). Fits that reach the iteration +#' limit still contribute to the held out error, so check the settings in +#' \code{control} before interpreting the selection. +#' +#' @return +#' A list containing: +#' \item{lambda_1_min}{Common source-space penalty in the selected combination. +#' \code{NA} when separate penalties are used.} +#' \item{lambda_2_min}{Balance penalty in the selected combination.} +#' \item{lambda_1_row_min, lambda_1_col_min}{Row and column penalties in the +#' selected combination. Returned only when separate penalties are used.} +#' \item{mse_all}{Sum of held-out squared errors for each candidate combination. +#' For common penalties, rows correspond to \code{lambda_1_all} and columns to +#' \code{lambda_2_all}. For separate penalties, the three dimensions correspond +#' to the row, column, and balance grids, in that order. Divide by the number +#' of observed target entries to obtain mean squared errors.} +#' \item{r}{Rank used in the fits.} #' #' @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. #' +#' @seealso +#' \code{\link{learner}} for fitting with the selected penalties. +#' #' @examples #' res <- cv.learner(Y_source = dat_highsim$Y_source, #' Y_target = dat_highsim$Y_target, @@ -31,139 +73,273 @@ #' lambda_2_all = c(1, 10, 100), #' step_size = 0.003) #' +#' # Refit using the selected penalties and all observed target entries. +#' fit <- learner(dat_highsim$Y_source, dat_highsim$Y_target, r = res$r, +#' lambda_1 = res$lambda_1_min, lambda_2 = res$lambda_2_min, +#' step_size = 0.003) +#' fit$convergence_criterion +#' +#' # Multiple sources with fixed demonstration weights. +#' sources <- dat_multisource$Y_source +#' target <- dat_multisource$Y_target +#' weights <- c(0.7, 0.3) +#' +#' set.seed(815) +#' res_separate <- cv.learner( +#' Y_source = sources, Y_target = target, r = 3, +#' source_weights = weights, +#' lambda_1_row_all = c(1, 10), lambda_1_col_all = c(1, 3), +#' lambda_2_all = c(0.1, 1), step_size = 0.0003, +#' n_folds = 2, control = list(max_iter = 100)) +#' dim(res_separate$mse_all) # row grid x column grid x balance grid #' +#' fit <- learner( +#' Y_source = sources, Y_target = target, r = res_separate$r, +#' source_weights = weights, +#' lambda_1_row = res_separate$lambda_1_row_min, +#' lambda_1_col = res_separate$lambda_1_col_min, +#' lambda_2 = res_separate$lambda_2_min, step_size = 0.0003, +#' control = list(max_iter = 100)) +#' fit$convergence_criterion +#' # This short example demonstrates the workflow; check convergence before +#' # interpreting results or comparing statistical performance. #' #' @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.') - } +cv.learner <- function(Y_source, Y_target, r, lambda_1_all = NULL, lambda_2_all, + step_size, n_folds = 4, n_cores = 1, control = list(), + lambda_1_row_all = NULL, lambda_1_col_all = NULL, + source_weights = NULL) { + sources <- prepare_sources(Y_source, Y_target, source_weights) + r <- resolve_rank(if (missing(r)) NULL else r, sources$matrices[[1]]) + control <- prepare_control(control, step_size) + validate_penalty_grid(lambda_2_all, 'lambda_2_all') - 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) - } - 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 - } - - if (n_folds < 2){ - stop('n_folds must be 2 or greater.') + # A common penalty uses a 2D grid; independent penalties use a 3D grid. + separate <- !is.null(lambda_1_row_all) || !is.null(lambda_1_col_all) + if (separate) { + if (!is.null(lambda_1_all)) { + stop('Provide either lambda_1_all or both separate grids, not both modes.') + } + if (is.null(lambda_1_row_all) || is.null(lambda_1_col_all)) { + stop('Both lambda_1_row_all and lambda_1_col_all must be provided together.') + } + validate_penalty_grid(lambda_1_row_all, 'lambda_1_row_all') + validate_penalty_grid(lambda_1_col_all, 'lambda_1_col_all') + grid <- expand.grid(row = lambda_1_row_all, col = lambda_1_col_all, + balance = lambda_2_all) + grid_dim <- c(length(lambda_1_row_all), length(lambda_1_col_all), length(lambda_2_all)) + } else { + validate_penalty_grid(lambda_1_all, 'lambda_1_all') + grid <- expand.grid(row = lambda_1_all, balance = lambda_2_all) + grid$col <- grid$row + grid_dim <- c(length(lambda_1_all), length(lambda_2_all)) } - ## Creating training and testing data sets + # Only observed target entries can be validation entries. 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)] + validate_count(n_folds, 'n_folds', lower = 2, upper = n_indices) + validate_count(n_cores, 'n_cores') + indices <- sample(available_indices, size = n_indices, replace = FALSE) - 1L + index_set <- vector('list', n_folds) + for (fold in seq_len(n_folds)) { + start <- floor((fold - 1) * n_indices / n_folds) + 1L + end <- floor(fold * n_indices / n_folds) + index_set[[fold]] <- indices[seq.int(start, end)] } - 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, - index_set) - return(result) + scores <- cv_learner_cpp(sources$matrices, Y_target, sources$weights, + grid$row, grid$col, grid$balance, step_size, + control$max_iter, control$threshold, n_cores, r, + control$max_value, index_set) + if (!any(is.finite(scores))) stop('No candidate produced finite validation errors.') + if (any(!is.finite(scores))) warning('Some candidates produced non-finite validation errors and were excluded.') + scores[!is.finite(scores)] <- Inf + best <- which.min(scores) + # Keep the historical mse_all scale: summed held-out squared errors. + mse_all <- array(scores, dim = grid_dim) + if (!separate) { + return(list(lambda_1_min = grid$row[best], lambda_2_min = grid$balance[best], + mse_all = mse_all, r = r)) + } + list(lambda_1_min = NA_real_, lambda_2_min = grid$balance[best], + lambda_1_row_min = grid$row[best], lambda_1_col_min = grid$col[best], + mse_all = mse_all, r = r) } #' 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.} +#' @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 +#' @param Y_source Matrix containing the source population data, with the same rows +#' and columns as \code{Y_target}. For multiple sources, supply a list of +#' matrices. +#' @param Y_target Matrix containing the target population data. Missing entries +#' are allowed. +#' @param r Rank of the knowledge graphs. If omitted, ScreeNOT (Donoho et al. 2023) +#' is applied to the source population data to select the rank. +#' @param lambda_1 Common source-space penalty. Larger values borrow more +#' information from the source population; zero does not use the source spaces. +#' Use \code{\link{cv.learner}} to select penalty values. +#' @param lambda_2 Balance penalty. Larger values more strongly encourage +#' \eqn{U} and \eqn{V} to have comparable scale. +#' @param step_size Step size in the optimization algorithm. Smaller values may +#' reduce oscillation but require more iterations. +#' @param control List of settings controlling when optimization stops: +#' \describe{ +#' \item{\code{max_iter}}{Maximum number of updates; default \code{100}.} +#' \item{\code{threshold}}{Stop when the absolute change between consecutive +#' objective values is smaller than this value; default \code{0.001}.} +#' \item{\code{max_value}}{Stop if the objective grows beyond this multiple of +#' its preceding value; default \code{10}.}} +#' @param lambda_1_row,lambda_1_col Separate source-space penalties for the row +#' factor \eqn{U} and column factor \eqn{V}, respectively. Use these when the two +#' sides should borrow different amounts of source information. Supply both +#' instead of \code{lambda_1}; equal values reproduce the common penalty. +#' @param source_weights Relative weights of the sources, in the order of +#' \code{Y_source}. Default is equal weights. #' +#' @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}. +#' 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 rank-\eqn{r} 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}. +#' 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 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\{\|UV^\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\}} #' -#' 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. +#' where \eqn{\mathcal{P}_\perp(\hat U_1) = I - \hat U_1\hat U_1^\top} and +#' \eqn{\mathcal{P}_\perp(\hat V_1) = I - \hat V_1\hat V_1^\top}. The four terms measure +#' target fit, departures from the two source spaces, and imbalance between the +#' target factors, respectively. If \code{Y_target} has missing entries, the +#' first term is computed over the observed entries and divided by the observed +#' fraction (the number of observed entries divided by \eqn{pq}). The penalty +#' terms are unchanged. The same adjustment applies to the extensions below. +#' +#' At each iteration, \eqn{U} and \eqn{V} are updated by a gradient +#' descent step, starting from the truncated SVD of the source population data. +#' These updates are repeated until convergence. +#' +#' \strong{Extensions:} +#' +#' Separate penalties \code{lambda_1_row} and \code{lambda_1_col} replace +#' \eqn{\lambda_1} in the second and third terms. For multiple sources, the +#' source projectors are replaced by weighted combinations: +#' +#' \deqn{\bar P_U=\sum_{k=1}^K w_k\hat U_k\hat U_k^\top,\qquad +#' \bar P_V=\sum_{k=1}^K w_k\hat V_k\hat V_k^\top.} +#' +#' where \eqn{\hat U_k} and \eqn{\hat V_k} come from the truncated SVD of +#' the \eqn{k}th source and the weights \eqn{w_k} are given by +#' \code{source_weights}. The full objective is +#' +#' \deqn{\begin{aligned} +#' L(U,V) ={}& \|UV^\top-Y_0\|_F^2 \\ +#' &+\lambda_{1,\mathrm{row}}\|(I_p-\bar P_U)U\|_F^2 \\ +#' &+\lambda_{1,\mathrm{col}}\|(I_q-\bar P_V)V\|_F^2 \\ +#' &+\lambda_2\|U^\top U-V^\top V\|_F^2. +#' \end{aligned}} +#' +#' where \eqn{I_p} and \eqn{I_q} are identity matrices. Equal space penalties +#' give the common-penalty version. With multiple sources, the first source is +#' used to initialize the algorithm and, if \code{r} is not specified, to +#' select the rank. +#' +#' @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}{stopping code: \code{1}, the absolute objective-change threshold was met; \code{2}, the iteration limit was reached; \code{3}, excessive objective growth or non-finite objective/factors were detected.} +#' \item{r}{rank value used.} #' #' @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. #' +#' @seealso +#' \code{\link{cv.learner}} for penalty selection; +#' \code{\link{dlearner}} for direct projection onto source spaces. +#' #' @examples +#' # Single source with a common source-space penalty. #' res <- learner(Y_source = dat_highsim$Y_source, #' Y_target = dat_highsim$Y_target, #' lambda_1 = 1, lambda_2 = 1, #' step_size = 0.003) #' +#' # Same source, with separate penalties for the two factors. +#' res_separate <- learner(Y_source = dat_highsim$Y_source, +#' Y_target = dat_highsim$Y_target, +#' lambda_1_row = 10, lambda_1_col = 1, +#' lambda_2 = 1, step_size = 0.003) +#' +#' # Multiple sources with fixed demonstration weights. +#' res_multi <- learner( +#' Y_source = dat_multisource$Y_source, +#' Y_target = dat_multisource$Y_target, r = 3, +#' source_weights = c(0.7, 0.3), +#' lambda_1_row = 10, lambda_1_col = 1, lambda_2 = 1, +#' step_size = 0.0003, control = list(max_iter = 100)) +#' dim(res_multi$learner_estimate) +#' res_multi$convergence_criterion +#' # Short runs demonstrate usage; they need not meet the stopping threshold. +#' #' @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.') +learner <- function(Y_source, Y_target, r, lambda_1 = NULL, lambda_2, step_size, + control = list(), lambda_1_row = NULL, lambda_1_col = NULL, + source_weights = NULL) { + sources <- prepare_sources(Y_source, Y_target, source_weights) + r <- resolve_rank(if (missing(r)) NULL else r, sources$matrices[[1]]) + control <- prepare_control(control, step_size) + # Resolve the original common penalty or the two independent penalties. + separate <- !is.null(lambda_1_row) || !is.null(lambda_1_col) + if (separate){ + if (!is.null(lambda_1)){ + stop('Provide either lambda_1 or both lambda_1_row and lambda_1_col, not both modes.') + } + if (is.null(lambda_1_row) || is.null(lambda_1_col)){ + stop('Both lambda_1_row and lambda_1_col must be provided together.') + } + validate_space_penalty(lambda_1_row, 'lambda_1_row') + validate_space_penalty(lambda_1_col, 'lambda_1_col') + } else { + if (is.null(lambda_1)){ + stop('Provide lambda_1 or both lambda_1_row and lambda_1_col.') + } + validate_space_penalty(lambda_1, 'lambda_1') + lambda_1_row <- lambda_1 + lambda_1_col <- lambda_1 } + validate_space_penalty(lambda_2, 'lambda_2') - 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) - } - 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, control$max_iter, control$threshold, control$max_value) + result <- learner_cpp(sources$matrices, Y_target, sources$weights, r, + lambda_1_row, lambda_1_col, lambda_2, step_size, + control$max_iter, control$threshold, control$max_value) return(result) } - +validate_space_penalty <- function(value, name) { + if (!is.numeric(value) || length(value) != 1L || + !is.finite(value) || value < 0){ + stop(name, ' must be a finite nonnegative numeric scalar.', call. = FALSE) + } +} diff --git a/data-raw/gen-data.R b/data-raw/gen-data.R index 32d302a..53a7863 100644 --- a/data-raw/gen-data.R +++ b/data-raw/gen-data.R @@ -55,3 +55,38 @@ dat_modsim <- list(Y_source = Y_source, Theta_target = Theta_target) usethis::use_data(dat_modsim, overwrite = TRUE) + +################################################################################ +## Multiple sources with different noise levels +################################################################################ + +# This block can be run on its own without regenerating the existing datasets. +local({ + set.seed(20261014) + p <- 100; q <- 50; r <- 3 + + Theta_target <- matrix(rnorm(p * r), p, r) %*% + matrix(rnorm(r * q), r, q) + target_svd <- svd(Theta_target, nu = r, nv = r) + # Both sources share the target's latent spaces, with reversed signal strengths. + Theta_source <- target_svd$u %*% + diag(rev(target_svd$d[seq_len(r)]), nrow = r, ncol = r) %*% + t(target_svd$v) + + # Variances, not standard deviations. Each call draws independent noise. + noise_variance <- c(target = 0.10, A = 0.01, B = 0.50) + Y_target <- Theta_target + matrix( + rnorm(p * q, sd = sqrt(noise_variance['target'])), p, q) + Y_A <- Theta_source + matrix( + rnorm(p * q, sd = sqrt(noise_variance['A'])), p, q) + Y_B <- Theta_source + matrix( + rnorm(p * q, sd = sqrt(noise_variance['B'])), p, q) + + dat_multisource <- list( + Y_source = list(A = Y_A, B = Y_B), + Y_target = Y_target, + Theta_source = list(A = Theta_source, B = Theta_source), + Theta_target = Theta_target + ) + usethis::use_data(dat_multisource, overwrite = TRUE, version = 2) +}) diff --git a/data/dat_multisource.rda b/data/dat_multisource.rda new file mode 100644 index 0000000..21395b4 Binary files /dev/null and b/data/dat_multisource.rda differ diff --git a/inst/examples/multiple-sources.R b/inst/examples/multiple-sources.R new file mode 100644 index 0000000..fd72a87 --- /dev/null +++ b/inst/examples/multiple-sources.R @@ -0,0 +1,88 @@ +# Demonstration using the built-in simulated multiple-source data. +# Each population has independently generated noise; source A is more precise. +sources <- learner::dat_multisource$Y_source +target <- learner::dat_multisource$Y_target +# Illustrative weights, not weights estimated or tuned from these data. +weights <- c(0.7, 0.3) + +# Tune the U-space, V-space, and balance penalties separately. +set.seed(815) +cv_result <- learner::cv.learner( + Y_source = sources, Y_target = target, r = 3, + source_weights = weights, + lambda_1_row_all = c(1, 10), + lambda_1_col_all = c(1, 3), + lambda_2_all = c(0.1, 1), + step_size = 0.0003, n_folds = 3, + control = list(max_iter = 2000) +) +print(c(row = cv_result$lambda_1_row_min, + col = cv_result$lambda_1_col_min, + balance = cv_result$lambda_2_min)) +print(dim(cv_result$mse_all)) # row grid x column grid x balance grid + +# Refit with the selected penalties using all observed target entries. +fit <- learner::learner( + Y_source = sources, Y_target = target, r = cv_result$r, + source_weights = weights, + lambda_1_row = cv_result$lambda_1_row_min, + lambda_1_col = cv_result$lambda_1_col_min, + lambda_2 = cv_result$lambda_2_min, + step_size = 0.0003, control = list(max_iter = 5000) +) +print(fit$learner_estimate[1:5, 1:5]) +print(fit$convergence_criterion) +# 1: objective-change threshold; 2: iteration limit; 3: excessive growth/breakdown. +# This example demonstrates the workflow, not scientific performance or convergence. + +# D-LEARNER also accepts multiple sources with the same weights. +direct_fit <- learner::dlearner(sources, target, r = 3, source_weights = weights) +print(dim(direct_fit$dlearner_estimate)) + +# Compare source choices using D-LEARNER at the same fixed rank. +# D-LEARNER avoids optimization and penalty tuning as additional differences. +# Truth is available only because this is a simulated dataset. +source_choices <- list( + A_only = sources$A, + B_only = sources$B, + Weighted_A_B = sources +) +comparison <- do.call(rbind, lapply(names(source_choices), function(label) { + estimate <- learner::dlearner( + source_choices[[label]], target, r = 3, + source_weights = if (label == 'Weighted_A_B') weights else NULL + )$dlearner_estimate + data.frame(source = label, + frobenius_error = norm(estimate - learner::dat_multisource$Theta_target, 'F')) +})) +print(comparison, row.names = FALSE) +# Lower Frobenius error is better. These are single-dataset results, not CV MSE. +# The weights 0.7/0.3 are illustrative; combining sources need not beat source A. + +# A second setting: both sources have appreciable noise. +# Keep the built-in data unchanged; variances are target 0.10, A 0.25, B 0.50. +# Reuse the true signals, with new independent noise for every population. +set.seed(20261016) +add_noise <- function(signal, variance) { + signal + matrix(rnorm(length(signal), sd = sqrt(variance)), + nrow(signal), ncol(signal)) +} +target_noisy <- add_noise(learner::dat_multisource$Theta_target, variance = 0.10) +sources_noisy <- list( + A = add_noise(learner::dat_multisource$Theta_source$A, variance = 0.25), + B = add_noise(learner::dat_multisource$Theta_source$B, variance = 0.50) +) +estimates_noisy <- list( + Target_only = learner::dlearner(target_noisy, target_noisy, r = 3)$dlearner_estimate, + A_only = learner::dlearner(sources_noisy$A, target_noisy, r = 3)$dlearner_estimate, + B_only = learner::dlearner(sources_noisy$B, target_noisy, r = 3)$dlearner_estimate, + Weighted_A_B = learner::dlearner(sources_noisy, target_noisy, r = 3, + source_weights = c(0.7, 0.3))$dlearner_estimate +) +errors_noisy <- vapply(estimates_noisy, function(estimate) { + norm(estimate - learner::dat_multisource$Theta_target, "F") +}, numeric(1)) +print(round(errors_noisy, 3)) +# Here the weighted combination improves on either source alone, but the +# target-only SVD has the smallest error. Fixed weights and one noise draw +# illustrate a possible outcome, not a general performance guarantee. diff --git a/man/cv.learner.Rd b/man/cv.learner.Rd index 8c0bd0c..609ca27 100644 --- a/man/cv.learner.Rd +++ b/man/cv.learner.Rd @@ -8,45 +8,87 @@ cv.learner( Y_source, Y_target, r, - lambda_1_all, + lambda_1_all = NULL, lambda_2_all, step_size, n_folds = 4, n_cores = 1, - control = list() + control = list(), + lambda_1_row_all = NULL, + lambda_1_col_all = NULL, + source_weights = NULL ) } \arguments{ -\item{Y_source}{matrix containing the source population data, as in \code{\link{learner}}} +\item{Y_source}{Source population data, as in \code{\link{learner}}. +Supply a matrix for one source or a list of matrices for multiple sources.} -\item{Y_target}{matrix containing the target population data, as in \code{\link{learner}}} +\item{Y_target}{Target population data matrix, as in \code{\link{learner}}. +Only observed entries are used to evaluate predictions.} -\item{r}{(optional) integer specifying the rank of the knowledge graphs, as in \code{\link{learner}}} +\item{r}{Rank used in every fold, as in \code{\link{learner}}. It is not +selected by cross-validation.} -\item{lambda_1_all}{vector of numerics specifying the candidate values of \eqn{\lambda_1} (see Details)} +\item{lambda_1_all}{Candidate values for the common source-space penalty +\eqn{\lambda_1}. All combinations with \code{lambda_2_all} are evaluated.} -\item{lambda_2_all}{vector of numerics specifying the candidate values of \eqn{\lambda_2} (see Details)} +\item{lambda_2_all}{Candidate values for the balance penalty \eqn{\lambda_2}.} -\item{step_size}{numeric scalar specifying the step size for the Newton steps in the numerical optimization algorithm, as in \code{\link{learner}}} +\item{step_size}{Step size for each fit, as in \code{\link{learner}}.} -\item{n_folds}{an integer specify the number of cross-validation folds. The default is \code{4}.} +\item{n_folds}{Number of cross-validation folds. Default is \code{4}.} -\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{n_cores}{Number of CPU cores used to evaluate the candidate penalty +combinations in parallel. 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{control}{Optimization settings for each fit; see \code{\link{learner}}.} + +\item{lambda_1_row_all, lambda_1_col_all}{Candidate values for separate +penalties on the row and column factors. Supply both instead of +\code{lambda_1_all}; all combinations with \code{lambda_2_all} are evaluated.} + +\item{source_weights}{Relative source weights, as in \code{\link{learner}}. The weights +are not selected by cross-validation. Default is equal weights.} } \value{ -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.} +A list containing: +\item{lambda_1_min}{Common source-space penalty in the selected combination. +\code{NA} when separate penalties are used.} +\item{lambda_2_min}{Balance penalty in the selected combination.} +\item{lambda_1_row_min, lambda_1_col_min}{Row and column penalties in the +selected combination. Returned only when separate penalties are used.} +\item{mse_all}{Sum of held-out squared errors for each candidate combination. +For common penalties, rows correspond to \code{lambda_1_all} and columns to +\code{lambda_2_all}. For separate penalties, the three dimensions correspond +to the row, column, and balance grids, in that order. Divide by the number +of observed target entries to obtain mean squared errors.} +\item{r}{Rank used in the fits.} } \description{ -This function performs k-fold cross-validation to select the nuisance parameters \eqn{(\lambda_1, \lambda_2)} for \code{\link{learner}}. +Performs k-fold cross-validation to select the source-space and balance +penalties for \code{\link{learner}}. Returns prediction errors and the +penalty combination with the smallest error. } \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. +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 comparing +the entries in the test data sets with those imputed from the LEARNER +estimates. See McGrath et al. (2024) for further details. + +With separate penalties, the procedure selects the row, column, and balance +penalties together. Both modes can be used with multiple sources. Call +\code{set.seed()} first to reproduce the random partition. + +To obtain the final estimate, apply \code{\link{learner}} with the selected +penalties to the full data (see Examples). Fits that reach the iteration +limit still contribute to the held out error, so check the settings in +\code{control} before interpreting the selection. } \examples{ res <- cv.learner(Y_source = dat_highsim$Y_source, @@ -55,9 +97,41 @@ res <- cv.learner(Y_source = dat_highsim$Y_source, lambda_2_all = c(1, 10, 100), step_size = 0.003) +# Refit using the selected penalties and all observed target entries. +fit <- learner(dat_highsim$Y_source, dat_highsim$Y_target, r = res$r, + lambda_1 = res$lambda_1_min, lambda_2 = res$lambda_2_min, + step_size = 0.003) +fit$convergence_criterion +# Multiple sources with fixed demonstration weights. +sources <- dat_multisource$Y_source +target <- dat_multisource$Y_target +weights <- c(0.7, 0.3) + +set.seed(815) +res_separate <- cv.learner( + Y_source = sources, Y_target = target, r = 3, + source_weights = weights, + lambda_1_row_all = c(1, 10), lambda_1_col_all = c(1, 3), + lambda_2_all = c(0.1, 1), step_size = 0.0003, + n_folds = 2, control = list(max_iter = 100)) +dim(res_separate$mse_all) # row grid x column grid x balance grid + +fit <- learner( + Y_source = sources, Y_target = target, r = res_separate$r, + source_weights = weights, + lambda_1_row = res_separate$lambda_1_row_min, + lambda_1_col = res_separate$lambda_1_col_min, + lambda_2 = res_separate$lambda_2_min, step_size = 0.0003, + control = list(max_iter = 100)) +fit$convergence_criterion +# This short example demonstrates the workflow; check convergence before +# interpreting results or comparing statistical performance. } \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. } +\seealso{ +\code{\link{learner}} for fitting with the selected penalties. +} diff --git a/man/dat_multisource.Rd b/man/dat_multisource.Rd new file mode 100644 index 0000000..32bc0f1 --- /dev/null +++ b/man/dat_multisource.Rd @@ -0,0 +1,54 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data-multisource.R +\docType{data} +\name{dat_multisource} +\alias{dat_multisource} +\title{Simulated data set: Multiple sources with different noise levels} +\format{ +A list with four components. All matrices have 100 rows and 50 columns, + with corresponding rows and columns aligned across populations. +\describe{ + \item{\code{Y_source}}{A named list of observed source matrices, \code{A} + and \code{B}, suitable for the \code{Y_source} argument.} + \item{\code{Y_target}}{The observed target matrix.} + \item{\code{Theta_source}}{A named list of the true rank-3 source matrices, + \code{A} and \code{B}. The two true matrices are identical.} + \item{\code{Theta_target}}{The true rank-3 target matrix.} +} +} +\usage{ +dat_multisource +} +\description{ +Simulated matrices for one target population and two source populations. +The sources share the target's latent row and column spaces, but source A +is observed more precisely than source B. This data set illustrates fitting +LEARNER and D-LEARNER with multiple sources. +} +\details{ +The source signal is obtained by reversing the three nonzero singular values +of the target signal while keeping its singular vectors. Each observed matrix +is generated by adding independent, mean-zero normal noise to its own signal. +Noise is independent across entries and populations, with variances 0.10 for +the target, 0.01 for source A, and 0.50 for source B. Thus, source A has less +noise than the target and source B has more. Source B is not generated from +the observed source A matrix. + +The true signal matrices are included for illustration and estimation-error +calculations; they are not inputs to the fitting functions. These simulated +data are a usage example, not a reproduction of the paper's simulation study. +} +\examples{ +data(dat_multisource) +fit <- dlearner( + Y_source = dat_multisource$Y_source, + Y_target = dat_multisource$Y_target, + r = 3, source_weights = c(0.7, 0.3) +) +dim(fit$dlearner_estimate) +} +\seealso{ +\code{\link{dat_highsim}}, \code{\link{dat_modsim}}, + \code{\link{learner}}, \code{\link{dlearner}} +} +\keyword{datasets} diff --git a/man/dlearner.Rd b/man/dlearner.Rd index f1c4563..4d9bbbb 100644 --- a/man/dlearner.Rd +++ b/man/dlearner.Rd @@ -2,16 +2,23 @@ % Please edit documentation in R/dlearner.R \name{dlearner} \alias{dlearner} -\title{Latent space-based transfer learning} +\title{Direct projection transfer learning} \usage{ -dlearner(Y_source, Y_target, r) +dlearner(Y_source, Y_target, r, source_weights = NULL) } \arguments{ -\item{Y_source}{matrix containing the source population data} +\item{Y_source}{Matrix containing the source population data, with the same rows +and columns as \code{Y_target}. For multiple sources, supply a list of +matrices.} -\item{Y_target}{matrix containing the target population data} +\item{Y_target}{Matrix containing the target population data. Missing entries +are not allowed.} -\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{r}{Rank of the knowledge graphs. If omitted, ScreeNOT (Donoho et al. 2023) +is applied to the source population data to select the rank.} + +\item{source_weights}{Relative weights of the sources, in the order of +\code{Y_source}. Default is equal weights.} } \value{ A list with the following components: @@ -19,26 +26,65 @@ A list with the following components: \item{r}{rank value used.} } \description{ -This function applies the Direct project LatEnt spAce-based tRaNsfer lEaRning (D-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 applies the Direct projection LatEnt spAce-based tRaNsfer +lEaRning (D-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}. +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 rank-\eqn{r} 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}. +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{\hat{U}_{1}^{\top}\hat{U}_{1} Y_0 \hat{V}_{1}^{\top}\hat{V}_{1}}. +This method estimates \eqn{\Theta_0} by + +\deqn{\hat U_1\hat U_1^\top Y_0 \hat V_1\hat V_1^\top.} + +That is, the target population data are projected onto the latent row and +column spaces of the source population, without iterative fitting or penalty +selection. + +\strong{Multiple sources:} + +Supply a list of source matrices to combine information from several +populations. Their projection operators are averaged using +\code{source_weights}, and the estimate becomes + +\deqn{\widehat\Theta_0=\bar P_U Y_0\bar P_V,\qquad +\bar P_U=\sum_{k=1}^K w_k\hat U_k\hat U_k^\top,\qquad +\bar P_V=\sum_{k=1}^K w_k\hat V_k\hat V_k^\top.} + +where \eqn{\hat U_k} and \eqn{\hat V_k} come from the truncated SVD of +the \eqn{k}th source and the weights \eqn{w_k} are given by +\code{source_weights}. If \code{r} is not specified, the rank is selected +using the first source. } \examples{ res <- dlearner(Y_source = dat_highsim$Y_source, Y_target = dat_highsim$Y_target) +# Multiple sources with fixed demonstration weights. +res_multi <- dlearner( + Y_source = dat_multisource$Y_source, + Y_target = dat_multisource$Y_target, r = 3, + source_weights = c(0.7, 0.3)) +dim(res_multi$dlearner_estimate) + } \references{ 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. } +\seealso{ +\code{\link{learner}} for a target fit that can adapt the source +spaces through penalized estimation. +} diff --git a/man/learner.Rd b/man/learner.Rd index e1bdce4..a56809b 100644 --- a/man/learner.Rd +++ b/man/learner.Rd @@ -4,63 +4,160 @@ \alias{learner} \title{Latent space-based transfer learning} \usage{ -learner(Y_source, Y_target, r, lambda_1, lambda_2, step_size, control = list()) +learner( + Y_source, + Y_target, + r, + lambda_1 = NULL, + lambda_2, + step_size, + control = list(), + lambda_1_row = NULL, + lambda_1_col = NULL, + source_weights = NULL +) } \arguments{ -\item{Y_source}{matrix containing the source population data} +\item{Y_source}{Matrix containing the source population data, with the same rows +and columns as \code{Y_target}. For multiple sources, supply a list of +matrices.} -\item{Y_target}{matrix containing the target population data} +\item{Y_target}{Matrix containing the target population data. Missing entries +are allowed.} -\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{r}{Rank of the knowledge graphs. If omitted, ScreeNOT (Donoho et al. 2023) +is applied to the source population data to select the rank.} -\item{lambda_1}{numeric scalar specifying the value of \eqn{\lambda_1} (see Details)} +\item{lambda_1}{Common source-space penalty. Larger values borrow more +information from the source population; zero does not use the source spaces. +Use \code{\link{cv.learner}} to select penalty values.} -\item{lambda_2}{numeric scalar specifying the value of \eqn{\lambda_2} (see Details)} +\item{lambda_2}{Balance penalty. Larger values more strongly encourage +\eqn{U} and \eqn{V} to have comparable scale.} -\item{step_size}{numeric scalar specifying the step size for the Newton steps in the numerical optimization algorithm} +\item{step_size}{Step size in the optimization algorithm. Smaller values may +reduce oscillation but require more 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{control}{List of settings controlling when optimization stops: +\describe{ +\item{\code{max_iter}}{Maximum number of updates; default \code{100}.} +\item{\code{threshold}}{Stop when the absolute change between consecutive +objective values is smaller than this value; default \code{0.001}.} +\item{\code{max_value}}{Stop if the objective grows beyond this multiple of +its preceding value; default \code{10}.}}} + +\item{lambda_1_row, lambda_1_col}{Separate source-space penalties for the row +factor \eqn{U} and column factor \eqn{V}, respectively. Use these when the two +sides should borrow different amounts of source information. Supply both +instead of \code{lambda_1}; equal values reproduce the common penalty.} + +\item{source_weights}{Relative weights of the sources, in the order of +\code{Y_source}. Default is equal weights.} } \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{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}{stopping code: \code{1}, the absolute objective-change threshold was met; \code{2}, the iteration limit was reached; \code{3}, excessive objective growth or non-finite objective/factors were detected.} \item{r}{rank value used.} } \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. +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}. +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 rank-\eqn{r} 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}. +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 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\{\|UV^\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\hat U_1^\top} and +\eqn{\mathcal{P}_\perp(\hat V_1) = I - \hat V_1\hat V_1^\top}. The four terms measure +target fit, departures from the two source spaces, and imbalance between the +target factors, respectively. If \code{Y_target} has missing entries, the +first term is computed over the observed entries and divided by the observed +fraction (the number of observed entries divided by \eqn{pq}). The penalty +terms are unchanged. The same adjustment applies to the extensions below. + +At each iteration, \eqn{U} and \eqn{V} are updated by a gradient +descent step, starting from the truncated SVD of the source population data. +These updates are repeated until convergence. + +\strong{Extensions:} + +Separate penalties \code{lambda_1_row} and \code{lambda_1_col} replace +\eqn{\lambda_1} in the second and third terms. For multiple sources, the +source projectors are replaced by weighted combinations: -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. +\deqn{\bar P_U=\sum_{k=1}^K w_k\hat U_k\hat U_k^\top,\qquad + \bar P_V=\sum_{k=1}^K w_k\hat V_k\hat V_k^\top.} + +where \eqn{\hat U_k} and \eqn{\hat V_k} come from the truncated SVD of +the \eqn{k}th source and the weights \eqn{w_k} are given by +\code{source_weights}. The full objective is + +\deqn{\begin{aligned} +L(U,V) ={}& \|UV^\top-Y_0\|_F^2 \\ +&+\lambda_{1,\mathrm{row}}\|(I_p-\bar P_U)U\|_F^2 \\ +&+\lambda_{1,\mathrm{col}}\|(I_q-\bar P_V)V\|_F^2 \\ +&+\lambda_2\|U^\top U-V^\top V\|_F^2. +\end{aligned}} + +where \eqn{I_p} and \eqn{I_q} are identity matrices. Equal space penalties +give the common-penalty version. With multiple sources, the first source is +used to initialize the algorithm and, if \code{r} is not specified, to +select the rank. } \examples{ +# Single source with a common source-space penalty. res <- learner(Y_source = dat_highsim$Y_source, Y_target = dat_highsim$Y_target, lambda_1 = 1, lambda_2 = 1, step_size = 0.003) +# Same source, with separate penalties for the two factors. +res_separate <- learner(Y_source = dat_highsim$Y_source, + Y_target = dat_highsim$Y_target, + lambda_1_row = 10, lambda_1_col = 1, + lambda_2 = 1, step_size = 0.003) + +# Multiple sources with fixed demonstration weights. +res_multi <- learner( + Y_source = dat_multisource$Y_source, + Y_target = dat_multisource$Y_target, r = 3, + source_weights = c(0.7, 0.3), + lambda_1_row = 10, lambda_1_col = 1, lambda_2 = 1, + step_size = 0.0003, control = list(max_iter = 100)) +dim(res_multi$learner_estimate) +res_multi$convergence_criterion +# Short runs demonstrate usage; they need not meet the stopping threshold. + } \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. } +\seealso{ +\code{\link{cv.learner}} for penalty selection; +\code{\link{dlearner}} for direct projection onto source spaces. +} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index f39d9a3..0574cb1 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -12,43 +12,46 @@ 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, 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) { +List learner_cpp(const Rcpp::List& Y_source, const Eigen::MatrixXd& Y_target, const std::vector& source_weights, int r, double lambda1_row, double lambda1_col, 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 source_weightsSEXP, SEXP rSEXP, SEXP lambda1_rowSEXP, SEXP lambda1_colSEXP, 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; - Rcpp::traits::input_parameter< const Eigen::MatrixXd& >::type Y_source(Y_sourceSEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::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 source_weights(source_weightsSEXP); Rcpp::traits::input_parameter< int >::type r(rSEXP); - Rcpp::traits::input_parameter< double >::type lambda1(lambda1SEXP); + Rcpp::traits::input_parameter< double >::type lambda1_row(lambda1_rowSEXP); + Rcpp::traits::input_parameter< double >::type lambda1_col(lambda1_colSEXP); 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::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)); + rcpp_result_gen = Rcpp::wrap(learner_cpp(Y_source, Y_target, source_weights, r, lambda1_row, lambda1_col, 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, 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) { +std::vector cv_learner_cpp(const Rcpp::List& Y_source, const Eigen::MatrixXd& Y_target, const std::vector& source_weights, const std::vector& row_grid, const std::vector& col_grid, const std::vector& balance_grid, double step_size, 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 source_weightsSEXP, SEXP row_gridSEXP, SEXP col_gridSEXP, SEXP balance_gridSEXP, SEXP step_sizeSEXP, 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; - Rcpp::traits::input_parameter< const Eigen::MatrixXd& >::type Y_source(Y_sourceSEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::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< const std::vector& >::type source_weights(source_weightsSEXP); + Rcpp::traits::input_parameter< const std::vector& >::type row_grid(row_gridSEXP); + Rcpp::traits::input_parameter< const std::vector& >::type col_grid(col_gridSEXP); + Rcpp::traits::input_parameter< const std::vector& >::type balance_grid(balance_gridSEXP); 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::traits::input_parameter< double >::type max_value(max_valueSEXP); 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)); + rcpp_result_gen = Rcpp::wrap(cv_learner_cpp(Y_source, Y_target, source_weights, row_grid, col_grid, balance_grid, step_size, max_iter, threshold, n_cores, r, max_value, index_set)); return rcpp_result_gen; END_RCPP } @@ -64,8 +67,8 @@ 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, 12}, + {"_learner_learner_cpp", (DL_FUNC) &_learner_learner_cpp, 11}, + {"_learner_cv_learner_cpp", (DL_FUNC) &_learner_cv_learner_cpp, 13}, {"_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 dbf2ed3..f6d09b0 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -12,31 +12,56 @@ using namespace Rcpp; using namespace Eigen; -//-------------------------------------------------------------- -// Internal worker function: computes the learner estimate - thread safe for OMP -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, - double max_value) { - - 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())); - 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; - - // Precompute reusable matrices - Eigen::MatrixXd U_trunc_T = U_trunc.transpose(); - Eigen::MatrixXd V_trunc_T = V_trunc.transpose(); +// Source decompositions are prepared once on the R thread and shared read-only. +struct SourceSpaces { + std::vector u, v; + std::vector weights; + Eigen::MatrixXd initial_u, initial_v; +}; + +SourceSpaces prepare_spaces(const Rcpp::List &sources, + const std::vector &weights, int r) { + SourceSpaces spaces; + spaces.weights = weights; + for (R_xlen_t k = 0; k < sources.size(); ++k) { + Eigen::MatrixXd source = Rcpp::as(sources[k]); + Eigen::BDCSVD svd(source, Eigen::ComputeThinU | Eigen::ComputeThinV); + spaces.u.push_back(svd.matrixU().leftCols(r)); + spaces.v.push_back(svd.matrixV().leftCols(r)); + if (k == 0) { + Eigen::MatrixXd root = svd.singularValues().head(r).array().sqrt().matrix().asDiagonal(); + spaces.initial_u = spaces.u[0] * root; + spaces.initial_v = spaces.v[0] * root; + } + } + return spaces; +} + +// Apply I - sum(w_k B_k B_k') without allocating a p-by-p projector. +Eigen::MatrixXd space_residual(const std::vector &bases, + const std::vector &weights, + const Eigen::MatrixXd &value) { + Eigen::MatrixXd projected = Eigen::MatrixXd::Zero(value.rows(), value.cols()); + for (size_t k = 0; k < bases.size(); ++k) { + projected.noalias() += weights[k] * (bases[k] * (bases[k].transpose() * value)); + } + return value - projected; +} +// Workers use no R objects or allocations, including inside the OpenMP loop. +struct FitResult { + Eigen::MatrixXd estimate; + std::vector objective; + int criterion; +}; + +FitResult learner_worker(const SourceSpaces &spaces, + const Eigen::MatrixXd &Y_target, + double lambda1_row, double lambda1_col, double lambda2, + double step_size, int max_iter, double threshold, + double max_value) { + Eigen::MatrixXd U = spaces.initial_u; + Eigen::MatrixXd V = spaces.initial_v; double perc_nonmissing = 1.0 - (static_cast((Y_target.array().isNaN()).count()) / Y_target.size()); bool missing = Y_target.hasNaN(); @@ -47,8 +72,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 * (U - U_trunc * (U_trunc_T * U)).squaredNorm() + - lambda1 * (V - V_trunc * (V_trunc_T * V)).squaredNorm() + + lambda1_row * space_residual(spaces.u, spaces.weights, U).squaredNorm() + + lambda1_col * space_residual(spaces.v, spaces.weights, V).squaredNorm() + lambda2 * (U.transpose() * U - V.transpose() * V).squaredNorm(); } @@ -71,12 +96,20 @@ List learner_worker(const Eigen::MatrixXd &Y_source, adjusted_theta = adjusted_theta.array().isNaN().select(temp, adjusted_theta); } + Eigen::MatrixXd residual_u = space_residual(spaces.u, spaces.weights, U); + Eigen::MatrixXd residual_v = space_residual(spaces.v, spaces.weights, V); + // A weighted average of projectors is generally not idempotent. + // The derivative of ||A U||^2 is 2 A' A U; here A is symmetric. + Eigen::MatrixXd penalty_u = spaces.u.size() == 1 ? residual_u : + space_residual(spaces.u, spaces.weights, residual_u); + Eigen::MatrixXd penalty_v = spaces.v.size() == 1 ? residual_v : + space_residual(spaces.v, spaces.weights, residual_v); Eigen::MatrixXd grad_U = (2.0 / perc_nonmissing) * (U * V_tilde - adjusted_theta * V) - + lambda1 * 2 * (U - U_trunc * (U_trunc_T * U)) - + lambda2 * 4 * U * (U_tilde - V_tilde); + + lambda1_row * 2 * penalty_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 * (V - V_trunc * (V_trunc_T * V)) - + lambda2 * 4 * V * (V_tilde - U_tilde); + + lambda1_col * 2 * penalty_v + + lambda2 * 4 * V * (V_tilde - U_tilde); double grad_U_norm = grad_U.norm(); double grad_V_norm = grad_V.norm(); @@ -93,11 +126,14 @@ 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 * (U - U_trunc * (U_trunc_T * U)).squaredNorm() + - lambda1 * (V - V_trunc * (V_trunc_T * V)).squaredNorm() + + lambda1_row * space_residual(spaces.u, spaces.weights, U).squaredNorm() + + lambda1_col * space_residual(spaces.v, spaces.weights, V).squaredNorm() + lambda2 * (U.transpose() * U - V.transpose() * V).squaredNorm(); } obj_values.push_back(obj); + if (!std::isfinite(obj) || !U.allFinite() || !V.allFinite()) { + return FitResult{U_best * V_best.transpose(), obj_values, 3}; + } if (obj < obj_best) { obj_best = obj; @@ -117,79 +153,63 @@ List learner_worker(const Eigen::MatrixXd &Y_source, obj_init = obj; } - return List::create( - Named("learner_estimate") = U_best * V_best.transpose(), - Named("objective_values") = obj_values, - Named("convergence_criterion") = convergence_criterion - ); + return FitResult{U_best * V_best.transpose(), obj_values, convergence_criterion}; } -//-------------------------------------------------------------- -// 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, double 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") = out["learner_estimate"], - Named("objective_values") = out["objective_values"], - Named("convergence_criterion") = out["convergence_criterion"], - Named("r") = r - ); +List learner_cpp(const Rcpp::List &Y_source, const Eigen::MatrixXd &Y_target, + const std::vector &source_weights, int r, + double lambda1_row, double lambda1_col, double lambda2, + double step_size, int max_iter, double threshold, double max_value) { + SourceSpaces spaces = prepare_spaces(Y_source, source_weights, r); + FitResult out = learner_worker(spaces, Y_target, lambda1_row, lambda1_col, + lambda2, step_size, max_iter, threshold, max_value); + return List::create(Named("learner_estimate") = out.estimate, + Named("objective_values") = out.objective, + Named("convergence_criterion") = out.criterion, + Named("r") = r); } -//-------------------------------------------------------------- +// All three vectors contain one entry per candidate combination, in R array order. // [[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, 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(); - - 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 % p; - int col = idx / p; - Y_train(row, col) = std::numeric_limits::quiet_NaN(); - } - 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 % p; - int col = idx / p; - double diff = learner_estimate(row, col) - Y_target(row, col); - fold_mse += diff * diff; - } - mse += fold_mse; +std::vector cv_learner_cpp(const Rcpp::List &Y_source, + const Eigen::MatrixXd &Y_target, + const std::vector &source_weights, + const std::vector &row_grid, + const std::vector &col_grid, + const std::vector &balance_grid, + double step_size, int max_iter, double threshold, + int n_cores, int r, double max_value, + const std::vector> &index_set) { + SourceSpaces spaces = prepare_spaces(Y_source, source_weights, r); + int p = Y_target.rows(); + int n_grid = row_grid.size(); + std::vector scores(n_grid, 0.0); + +#pragma omp parallel for schedule(dynamic) num_threads(n_cores) + for (int g = 0; g < n_grid; ++g) { + double score = 0.0; + for (const auto &fold : index_set) { + Eigen::MatrixXd training = Y_target; + for (int idx : fold) training(idx % p, idx / p) = std::numeric_limits::quiet_NaN(); + FitResult fit = learner_worker(spaces, training, row_grid[g], col_grid[g], + balance_grid[g], step_size, max_iter, threshold, max_value); + if (!fit.estimate.allFinite() || + !std::all_of(fit.objective.begin(), fit.objective.end(), + [](double x) { return std::isfinite(x); })) { + score = std::numeric_limits::infinity(); + break; + } + double fold_score = 0.0; + for (int idx : fold) { + double difference = fit.estimate(idx % p, idx / p) - Y_target(idx % p, idx / p); + fold_score += difference * difference; } - mse_all(i, j) = mse; + score += fold_score; } + scores[g] = score; } - - Eigen::Index minRow, minCol; - mse_all.minCoeff(&minRow, &minCol); - - return List::create( - Named("lambda_1_min") = lambda1_all[minRow], - Named("lambda_2_min") = lambda2_all[minCol], - Named("mse_all") = mse_all, - Named("r") = r - ); + return scores; } // -------------------------------------------------------------- @@ -201,4 +221,3 @@ int omp_max_threads() { return 1; #endif } - diff --git a/tests/testthat/fixtures/README.md b/tests/testthat/fixtures/README.md new file mode 100644 index 0000000..ea3a456 --- /dev/null +++ b/tests/testthat/fixtures/README.md @@ -0,0 +1,54 @@ +# Fixed reference results for the extensions + +`multisource-v030.rds` contains inputs, fitting arguments, and expected outputs +computed by Sean McGrath's base-R **learnerv2 0.3.0** implementation from the +`learner-paper` repository's simulations materials. The tests do not create +expected values by calling the `learner` implementation under test. + +## Coverage + +- Three simulated 9-by-7 source matrices with weights 2, 3, 5 and fixed rank 2. +- Common penalties and unequal U/V penalties, for complete and missing targets. +- Complete fitted matrices and all 15 objective values, rank, and stopping code. +- Separate-penalty CV with 2-by-2-by-2 candidates, three folds, and fixed RNG + settings, for complete and missing targets. All scores and selected values + are checked; each saved CV case has a unique minimum. +- The complete multiple-source D-LEARNER estimate, also checked at generation + time against explicit dense weighted source operators. + +The saved inputs are synthetic. They are generated with seed 405 and then stored +in full, rather than regenerated during tests. CV uses seed 514 with RNG kinds +`Mersenne-Twister`, `Inversion`, and `Rejection`; its test restores the incoming +RNG state. The fixture records the reference version, source-file MD5 hashes, +R version, and RNG kinds so the source of the expected values can be audited. + +`expect_equal()` uses tolerance `1e-8` to allow ordinary floating-point and SVD +library differences. Matrix estimates and objective histories are checked, not +raw U/V factors, whose SVD signs are arbitrary. Rank is fixed to avoid depending +on future changes in automatic rank selection. + +The 15-update examples use `threshold = 0` and a large growth limit to compare a +fixed optimization path. They do not demonstrate convergence or statistical +accuracy. They do not assert that the two implementations have identical stopping +behavior for all settings. The analytic tests in `test-stopping-behavior.R` +separately check the current package's documented stopping rules. + +Normal tests only read this fixture and require neither `learnerv2` nor `foreach`. + +## Deliberate regeneration + +Run the generator manually from the package root, using an extracted, unmodified +copy of learnerv2 0.3.0. The generator requires `foreach`; it verifies the package +name and version and sources the reference R files without loading `learner`. + +```sh +Rscript tests/testthat/fixtures/generate-reference.R \ + /path/to/learnerv2 \ + tests/testthat/fixtures/multisource-v030.rds +``` + +The generator refuses to overwrite an existing baseline unless `--overwrite` is +added. Do not regenerate the fixture merely to make a failing test pass. First +understand whether a change is an intended method/interface change, an optimization +difference, or a bug; review any new baseline independently and update this note. +Neither testthat nor R CMD check runs the generator automatically. diff --git a/tests/testthat/fixtures/generate-reference.R b/tests/testthat/fixtures/generate-reference.R new file mode 100644 index 0000000..fa1b5c0 --- /dev/null +++ b/tests/testthat/fixtures/generate-reference.R @@ -0,0 +1,87 @@ +# Developer-only script. It is never run by testthat or R CMD check. +# Usage: Rscript generate-reference.R /path/to/learnerv2-0.3.0 /path/to/output.rds +# Add --overwrite only after reviewing an intentional change to the baseline. +args <- commandArgs(trailingOnly = TRUE) +if (length(args) < 2L) stop('Supply the extracted reference package directory and output .rds path.') +reference_dir <- normalizePath(args[1], mustWork = TRUE) +output <- args[2] +if (file.exists(output) && !('--overwrite' %in% args)) { + stop('Reference already exists; review the change before using --overwrite.') +} +description <- read.dcf(file.path(reference_dir, 'DESCRIPTION')) +stopifnot(description[1, 'Package'] == 'learnerv2', description[1, 'Version'] == '0.3.0') +# The reference CV implementation uses foreach even for sequential execution. +if (!requireNamespace('foreach', quietly = TRUE)) stop('Install foreach to generate this fixture.') +reference <- new.env(parent = asNamespace('foreach')) +reference_files <- file.path('R', c('learner.R', 'cv_learner.R', 'dlearner.R')) +for (f in reference_files) sys.source(file.path(reference_dir, f), envir = reference) + +# Inputs are saved in the fixture, so tests do not regenerate them from the RNG. +RNGkind('Mersenne-Twister', 'Inversion', 'Rejection') +set.seed(405) +sources <- list(matrix(rnorm(63), 9, 7), matrix(rnorm(63), 9, 7), + matrix(rnorm(63), 9, 7)) +target <- matrix(rnorm(63), 9, 7) +masked <- target +masked[c(3, 8, 19)] <- NA_real_ +targets <- list(complete = target, missing = masked) +control <- list(max_iter = 15L, threshold = 0, max_value = 1e6) +rank <- 2L +weights <- c(2, 3, 5) +step <- 0.01 +penalties <- list(common = list(lambda_1 = 3), + separate = list(lambda_1_row = 3, lambda_1_col = 0.25)) +grids <- list(lambda_1_row_all = c(0.25, 3), lambda_1_col_all = c(0.5, 2), + lambda_2_all = c(0, 1)) +fit_cases <- list() +cv_cases <- list() +for (target_name in names(targets)) { + for (penalty_name in names(penalties)) { + fit_args <- c(list(Y_source = sources, Y_target = targets[[target_name]], + r = rank, source_weights = weights, lambda_2 = 1, + step_size = step, control = control), penalties[[penalty_name]]) + reference_args <- fit_args + reference_args$r <- NULL + reference_args$r_source <- rank + reference_args$r_target <- rank + expected <- do.call(reference$learner, reference_args) + stopifnot(expected$convergence_criterion == 2L, + length(expected$objective_values) == control$max_iter, + all(is.finite(expected$objective_values))) + fit_cases[[paste(target_name, penalty_name, sep = '_')]] <- + list(args = fit_args, expected = expected) + } + cv_args <- c(list(Y_source = sources, Y_target = targets[[target_name]], + r = rank, source_weights = weights, step_size = step, + control = control, n_folds = 3L, n_cores = 1L), grids) + reference_args <- cv_args + reference_args$r <- NULL + reference_args$r_source <- rank + reference_args$r_target <- rank + set.seed(514) + expected <- do.call(reference$cv.learner, reference_args) + stopifnot(all(is.finite(expected$mse_all))) + sorted_scores <- sort(as.vector(expected$mse_all)) + stopifnot(sorted_scores[2] - sorted_scores[1] > 1e-5) + cv_cases[[target_name]] <- list(args = cv_args, seed = 514L, expected = expected) +} +direct_args <- list(Y_source = sources, Y_target = target, r = rank, + source_weights = weights) +direct_expected <- do.call(reference$dlearner, direct_args) +# Independently verify D-LEARNER by forming the dense weighted operators. +bases <- lapply(sources, svd, nu = rank, nv = rank) +w <- weights / sum(weights) +pu <- Reduce(`+`, Map(function(b, wk) wk * tcrossprod(b$u), bases, w)) +pv <- Reduce(`+`, Map(function(b, wk) wk * tcrossprod(b$v), bases, w)) +stopifnot(isTRUE(all.equal(direct_expected$dlearner_estimate, + pu %*% target %*% pv, tolerance = 1e-12))) +fixture <- list( + provenance = list(reference_package = 'learnerv2', reference_version = '0.3.0', + reference_md5 = setNames(as.character(tools::md5sum( + file.path(reference_dir, reference_files))), reference_files), + r_version = as.character(getRversion()), + rng_kind = RNGkind(), tolerance = 1e-8), + fit = fit_cases, cv = cv_cases, + direct = list(args = direct_args, expected = direct_expected)) +saveRDS(fixture, output, version = 2) +cat('Wrote independently generated reference fixture:', output, '\n') diff --git a/tests/testthat/fixtures/multisource-v030.rds b/tests/testthat/fixtures/multisource-v030.rds new file mode 100644 index 0000000..7377dc5 Binary files /dev/null and b/tests/testthat/fixtures/multisource-v030.rds differ diff --git a/tests/testthat/test-cv-multisource.R b/tests/testthat/test-cv-multisource.R new file mode 100644 index 0000000..1351de5 --- /dev/null +++ b/tests/testthat/test-cv-multisource.R @@ -0,0 +1,135 @@ +multi_example <- function() { + set.seed(705) + list(a = matrix(rnorm(35), 7, 5), b = matrix(rnorm(35), 7, 5), + y = matrix(rnorm(35), 7, 5)) +} + +# Rebuild the held-out predictions explicitly to check grid axes and masking. +manual_cv_scores <- function(sources, target, row, col, balance, weights = NULL, + seed = 92, folds = 3) { + set.seed(seed) + observed <- which(!is.na(target)) + indices <- sample(observed, length(observed), replace = FALSE) + groups <- lapply(seq_len(folds), function(fold) { + indices[seq.int(floor((fold - 1) * length(indices) / folds) + 1, + floor(fold * length(indices) / folds))] + }) + scores <- array(NA_real_, c(length(row), length(col), length(balance))) + for (i in seq_along(row)) for (j in seq_along(col)) for (k in seq_along(balance)) { + scores[i, j, k] <- sum(vapply(groups, function(held_out) { + train <- target + train[held_out] <- NA_real_ + fit <- learner(sources, train, r = 2, lambda_1_row = row[i], + lambda_1_col = col[j], lambda_2 = balance[k], + source_weights = weights, step_size = 0.01, + control = list(max_iter = 9, threshold = 0)) + sum((fit$learner_estimate[held_out] - target[held_out])^2) + }, numeric(1))) + } + scores +} + +test_that("independent CV axes and scores match explicit held-out predictions", { + d <- multi_example() + for (multi in c(FALSE, TRUE)) for (missing in c(FALSE, TRUE)) { + sources <- if (multi) list(d$a, d$b) else d$a + weights <- if (multi) c(2, 1) else NULL + y <- d$y + if (missing) y[c(3, 10, 14)] <- NA_real_ + row <- c(0, 3) + col <- c(0.2, 2, 5) + balance <- c(0, 1) + reference <- manual_cv_scores(sources, y, row, col, balance, weights) + set.seed(92) + actual <- cv.learner(sources, y, r = 2, lambda_1_row_all = row, + lambda_1_col_all = col, lambda_2_all = balance, + source_weights = weights, step_size = 0.01, n_folds = 3, + control = list(max_iter = 9, threshold = 0)) + expect_equal(actual$mse_all, reference, tolerance = 1e-10) + best <- arrayInd(which.min(reference), dim(reference)) + expect_equal(actual$lambda_1_row_min, row[best[1]]) + expect_equal(actual$lambda_1_col_min, col[best[2]]) + expect_equal(actual$lambda_2_min, balance[best[3]]) + expect_true(is.na(actual$lambda_1_min)) + } +}) + +test_that("tied CV is the diagonal of the independent grid", { + d <- multi_example() + fit <- function(...) { + set.seed(18) + cv.learner(list(d$a, d$b), d$y, r = 2, lambda_2_all = c(0, 1), + step_size = 0.01, control = list(max_iter = 8), ...) + } + tied <- fit(lambda_1_all = c(1, 3)) + separate <- fit(lambda_1_row_all = c(1, 3), lambda_1_col_all = c(1, 3)) + expect_named(tied, c('lambda_1_min', 'lambda_2_min', 'mse_all', 'r')) + expect_equal(dim(tied$mse_all), c(2L, 2L)) + for (i in 1:2) expect_equal(tied$mse_all[i, ], separate$mse_all[i, i, ], tolerance = 1e-10) + single <- fit(lambda_1_row_all = 1, lambda_1_col_all = 3) + expect_equal(dim(single$mse_all), c(1L, 1L, 2L)) + parallel <- fit(lambda_1_row_all = c(1, 3), lambda_1_col_all = c(1, 3), n_cores = 2) + expect_equal(parallel, separate, tolerance = 1e-12) +}) + +test_that("source representations and normalized weights are consistent", { + d <- multi_example() + fit <- function(source, weights = NULL) { + learner(source, d$y, r = 2, lambda_1_row = 3, lambda_1_col = 0.5, + lambda_2 = 1, step_size = 0.01, source_weights = weights, + control = list(max_iter = 12, threshold = 0)) + } + expect_equal(fit(d$a), fit(list(d$a)), tolerance = 1e-12) + expect_equal(fit(d$a), fit(list(d$a, d$a)), tolerance = 1e-10) + expect_equal(fit(d$a), fit(list(d$a, d$b), c(1, 0)), tolerance = 1e-10) + expect_equal(fit(list(d$a, d$b), c(2, 1)), + fit(list(d$a, d$b), c(2/3, 1/3)), tolerance = 1e-12) + expect_equal(fit(list(d$a, d$b)), fit(list(d$a, d$b), c(1, 1)), tolerance = 1e-12) + expect_gt(max(abs(fit(list(d$a, d$b), c(1, 0))$learner_estimate - + fit(list(d$a, d$b), c(1, 3))$learner_estimate)), 1e-6) +}) + +test_that("D-LEARNER matches weighted projectors and preserves matrix calls", { + d <- multi_example() + ua <- svd(d$a, nu = 2, nv = 2) + ub <- svd(d$b, nu = 2, nv = 2) + pu <- 0.7 * tcrossprod(ua$u) + 0.3 * tcrossprod(ub$u) + pv <- 0.7 * tcrossprod(ua$v) + 0.3 * tcrossprod(ub$v) + result <- dlearner(list(d$a, d$b), d$y, 2, c(7, 3)) + expect_equal(result$dlearner_estimate, pu %*% d$y %*% pv, tolerance = 1e-12) + expect_equal(dlearner(d$a, d$y, 2), dlearner(list(d$a), d$y, 2), tolerance = 1e-12) + expect_equal(dlearner(d$a, d$y, 2), dlearner(list(d$a, d$b), d$y, 2, c(1, 0)), tolerance = 1e-12) + expect_equal(result, dlearner(list(d$b, d$a), d$y, 2, c(3, 7)), tolerance = 1e-12) +}) + +test_that("invalid sources, weights, ranks, folds, and grids are rejected", { + d <- multi_example() + fit <- function(source = list(d$a, d$b), weights = NULL, rank = 2) { + learner(source, d$y, r = rank, source_weights = weights, + lambda_1 = 1, lambda_2 = 1, step_size = 0.01) + } + expect_error(fit(list()), 'nonempty list') + expect_error(fit(list(d$a, 2)), 'numeric matrix') + expect_error(fit(list(d$a, d$b[1:3, ])), 'same dimensions') + bad <- d$b; bad[1] <- NA_real_ + expect_error(fit(list(d$a, bad)), 'finite values') + for (w in list(c(0, 0), c(-1, 2), c(1, NA), c(1, Inf), 1, 'x')) { + expect_error(fit(weights = w), 'source_weights') + } + for (r in list(0, 1.5, 6, NA, Inf, c(1, 2))) expect_error(fit(rank = r), 'r must') + cv <- function(...) cv.learner(d$a, d$y, r = 2, step_size = 0.01, ...) + expect_error(cv(lambda_1_row_all = 1, lambda_2_all = 1), 'provided together') + expect_error(cv(lambda_1_all = 1, lambda_1_row_all = 1, + lambda_1_col_all = 1, lambda_2_all = 1), 'not both modes') + for (grid in list(numeric(), NA, Inf, -1, '1')) { + expect_error(cv(lambda_1_all = grid, lambda_2_all = 1), 'nonempty vector') + } + expect_error(cv(lambda_1_all = 1, lambda_2_all = numeric()), 'lambda_2_all') + for (folds in c(1, 36, 2.5, NA)) { + expect_error(cv(lambda_1_all = 1, lambda_2_all = 1, n_folds = folds), 'n_folds') + } + y <- matrix(NA_real_, 7, 5) + expect_error(learner(d$a, y, 2, 1, 1, 0.01), 'observed values') + y[1] <- 1 + expect_error(cv.learner(d$a, y, 2, 1, 1, 0.01), 'n_folds') +}) diff --git a/tests/testthat/test-failure-handling.R b/tests/testthat/test-failure-handling.R new file mode 100644 index 0000000..d7cb877 --- /dev/null +++ b/tests/testthat/test-failure-handling.R @@ -0,0 +1,99 @@ +# Small, deterministic matrices keep error-path tests independent of fixtures. +failure_example <- function() { + list(source = diag(c(4, 3, 2, 1)), target = diag(c(1, 2, 3, 4))) +} + +test_that('all fitting interfaces reject malformed target matrices', { + d <- failure_example() + infinite <- d$target + infinite[1, 1] <- Inf + calls <- list( + function(y) learner(d$source, y, r = 2, lambda_1 = 1, + lambda_2 = 1, step_size = 0.001), + function(y) cv.learner(d$source, y, r = 2, lambda_1_all = 1, + lambda_2_all = 1, step_size = 0.001), + function(y) dlearner(d$source, y, r = 2)) + for (target in list(as.vector(d$target), matrix('x', 4, 4), + matrix(numeric(), 0, 4), infinite)) { + for (fit in calls) { + expect_error(fit(target), 'Y_target must be a nonempty numeric matrix', + fixed = TRUE) + } + } +}) + +test_that('missing-target messages distinguish an empty target from D-LEARNER requirements', { + d <- failure_example() + missing <- d$target + missing[] <- NA_real_ + expect_error(learner(d$source, missing, r = 2, lambda_1 = 1, + lambda_2 = 1, step_size = 0.001), + '^Y_target must contain observed values\\.$') + expect_error(cv.learner(d$source, missing, r = 2, lambda_1_all = 1, + lambda_2_all = 1, step_size = 0.001), + '^Y_target must contain observed values\\.$') + d$target[1, 1] <- NA_real_ + expect_error(dlearner(d$source, d$target, r = 2), + 'dlearner requires a complete Y_target without missing values', fixed = TRUE) +}) + +test_that('learner and CV reject invalid optimization settings', { + d <- failure_example() + calls <- list( + function(...) learner(d$source, d$target, r = 2, lambda_1 = 1, + lambda_2 = 1, ...), + function(...) cv.learner(d$source, d$target, r = 2, lambda_1_all = 1, + lambda_2_all = 1, ...)) + for (fit in calls) { + expect_error(fit(step_size = 0.001, control = 3), + 'control must be a list', fixed = TRUE) + for (bad in c(0, -1, Inf, NA_real_)) { + expect_error(fit(step_size = bad), + 'step_size must be a finite positive number', fixed = TRUE) + expect_error(fit(step_size = 0.001, control = list(max_value = bad)), + 'control$max_value must be a finite positive number', fixed = TRUE) + } + } +}) + +test_that('a non-finite update stops immediately and retains the finite initialization', { + d <- failure_example() + # A finite but extreme step overflows the objective on the first update. + fit <- learner(d$source, d$target, r = 2, lambda_1 = 1, lambda_2 = 1, + step_size = 1e200, control = list(max_iter = 3)) + expect_equal(fit$convergence_criterion, 3) + expect_length(fit$objective_values, 1) + expect_false(is.finite(fit$objective_values[1])) + expect_true(all(is.finite(fit$learner_estimate))) + expect_equal(fit$learner_estimate, diag(c(4, 3, 0, 0)), tolerance = 1e-12) +}) + +test_that('CV reports when no candidate has a finite validation score', { + d <- failure_example() + set.seed(42) + expect_error( + cv.learner(d$source, d$target, r = 2, lambda_1_all = c(0, 1), + lambda_2_all = 1, step_size = 1e200, n_folds = 2, + control = list(max_iter = 3)), + 'No candidate produced finite validation errors', fixed = TRUE) +}) + +test_that('CV excludes a failed candidate and preserves the valid candidate score', { + d <- failure_example() + cv <- function(grid) { + set.seed(42) + cv.learner(d$source, d$target, r = 2, lambda_1_all = grid, + lambda_2_all = 1, step_size = 0.001, n_folds = 2, + control = list(max_iter = 3)) + } + reference <- cv(0) + # An extreme finite penalty causes non-finite updates in just that candidate. + expect_warning(result <- cv(c(1e308, 0)), + 'Some candidates produced non-finite validation errors and were excluded', + fixed = TRUE) + expect_equal(result$mse_all[1, 1], Inf) + expect_true(is.finite(result$mse_all[2, 1])) + expect_equal(result$mse_all[2, 1], reference$mse_all[1, 1], tolerance = 1e-12) + expect_equal(result$lambda_1_min, reference$lambda_1_min) + expect_equal(result$lambda_2_min, reference$lambda_2_min) +}) diff --git a/tests/testthat/test-main.R b/tests/testthat/test-main.R index 7cf7bcb..4a1313f 100644 --- a/tests/testthat/test-main.R +++ b/tests/testthat/test-main.R @@ -27,7 +27,7 @@ test_that("cv.learner output is correct", { 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) + expect_equal(result$mse_all, expected_result, tolerance = 1e-3) }) test_that("learner does not fail", { @@ -114,5 +114,5 @@ test_that("cv.learner output is correct with missing data", { 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) + expect_equal(result$mse_all, expected_result, tolerance = 1e-3) }) diff --git a/tests/testthat/test-multisource-data.R b/tests/testthat/test-multisource-data.R new file mode 100644 index 0000000..0cecc33 --- /dev/null +++ b/tests/testthat/test-multisource-data.R @@ -0,0 +1,31 @@ +test_that('multiple-source example contains aligned observations and rank-three signals', { + d <- dat_multisource + expect_named(d, c('Y_source', 'Y_target', 'Theta_source', 'Theta_target')) + expect_named(d$Y_source, c('A', 'B')) + expect_named(d$Theta_source, c('A', 'B')) + matrices <- c(d$Y_source, list(d$Y_target), d$Theta_source, list(d$Theta_target)) + for (m in matrices) { + expect_identical(dim(m), c(100L, 50L)) + expect_true(all(is.finite(m))) + } + for (m in c(d$Theta_source, list(d$Theta_target))) { + expect_equal(sum(svd(m, nu = 0, nv = 0)$d > 1e-8), 3) + } + expect_equal(d$Theta_source$A, d$Theta_source$B) + target <- svd(d$Theta_target, nu = 3, nv = 3) + source <- svd(d$Theta_source$A, nu = 3, nv = 3) + expect_equal(tcrossprod(target$u), tcrossprod(source$u), tolerance = 1e-10) + expect_equal(tcrossprod(target$v), tcrossprod(source$v), tolerance = 1e-10) +}) + +test_that('saved example has the intended noise scales without shared source noise', { + d <- dat_multisource + noise <- cbind(target = as.vector(d$Y_target - d$Theta_target), + A = as.vector(d$Y_source$A - d$Theta_source$A), + B = as.vector(d$Y_source$B - d$Theta_source$B)) + # Fixed saved data: broad checks catch variance/SD mistakes and reused noise. + ratios <- apply(noise, 2, var) / c(target = 0.10, A = 0.01, B = 0.50) + expect_true(all(ratios > 0.8 & ratios < 1.2)) + correlations <- cor(noise) + expect_true(all(abs(correlations[upper.tri(correlations)]) < 0.05)) +}) diff --git a/tests/testthat/test-reference-regression.R b/tests/testthat/test-reference-regression.R new file mode 100644 index 0000000..bbc28cb --- /dev/null +++ b/tests/testthat/test-reference-regression.R @@ -0,0 +1,51 @@ +# These expected values come from Sean's base-R learnerv2 0.3.0 implementation, +# not from the function under test. See fixtures/README.md for provenance. +reference_fixture <- function() { + readRDS(test_path('fixtures', 'multisource-v030.rds')) +} + +with_reference_rng <- function(seed, kind, code) { + old_kind <- RNGkind() + had_seed <- exists('.Random.seed', envir = .GlobalEnv, inherits = FALSE) + if (had_seed) old_seed <- get('.Random.seed', envir = .GlobalEnv) + on.exit({ + do.call(RNGkind, as.list(old_kind)) + if (had_seed) { + assign('.Random.seed', old_seed, envir = .GlobalEnv) + } else if (exists('.Random.seed', envir = .GlobalEnv, inherits = FALSE)) { + rm('.Random.seed', envir = .GlobalEnv) + } + }) + do.call(RNGkind, as.list(kind)) + set.seed(seed) + force(code) +} + +for (case_name in c('complete_common', 'complete_separate', + 'missing_common', 'missing_separate')) { + test_that(paste('multisource fit matches the saved reference:', case_name), { + fixture <- reference_fixture() + case <- fixture$fit[[case_name]] + actual <- do.call(learner, case$args) + # Includes every estimated entry, every objective value, rank, and stop code. + expect_equal(actual, case$expected, tolerance = fixture$provenance$tolerance) + }) +} + +for (case_name in c('complete', 'missing')) { + test_that(paste('separate-penalty CV matches the saved reference:', case_name), { + fixture <- reference_fixture() + case <- fixture$cv[[case_name]] + actual <- with_reference_rng(case$seed, fixture$provenance$rng_kind, + do.call(cv.learner, case$args)) + # Includes all eight scores and the selected penalties, not only dimensions. + expect_equal(actual, case$expected, tolerance = fixture$provenance$tolerance) + }) +} + +test_that('multisource D-LEARNER matches the saved reference', { + fixture <- reference_fixture() + actual <- do.call(dlearner, fixture$direct$args) + expect_equal(actual, fixture$direct$expected, + tolerance = fixture$provenance$tolerance) +}) diff --git a/tests/testthat/test-separate-penalties.R b/tests/testthat/test-separate-penalties.R new file mode 100644 index 0000000..878bd24 --- /dev/null +++ b/tests/testthat/test-separate-penalties.R @@ -0,0 +1,144 @@ +penalty_example <- function() { + set.seed(421) + list(source = matrix(rnorm(20), 5, 4), + target = matrix(rnorm(20), 5, 4)) +} + +test_that("equal separate penalties preserve common and positional calls", { + dat <- penalty_example() + for (with_missing in c(FALSE, TRUE)) { + target <- dat$target + if (with_missing) target[c(2, 7)] <- NA_real_ + control <- list(max_iter = 12, threshold = 0) + common <- learner(dat$source, target, 1, 3, 1, 0.01, control) + separate <- learner(dat$source, target, r = 1, lambda_2 = 1, + step_size = 0.01, control = control, + lambda_1_row = 3, lambda_1_col = 3) + expect_equal(separate, common, tolerance = 1e-12) + } +}) + +test_that("transposing data swaps the two independent space penalties", { + dat <- penalty_example() + fit <- function(source, target, row, col) { + learner(source, target, r = 1, lambda_1_row = row, + lambda_1_col = col, lambda_2 = 1, step_size = 0.02, + control = list(max_iter = 12, threshold = 0)) + } + original <- fit(dat$source, dat$target, 3, 0.25) + transposed <- fit(t(dat$source), t(dat$target), 0.25, 3) + expect_equal(original$learner_estimate, t(transposed$learner_estimate), + tolerance = 1e-9) + expect_equal(original$objective_values, transposed$objective_values, + tolerance = 1e-9) +}) + +# A finite-difference oracle checks the C++ updates against the objective, +# without copying the analytic gradient formulas used by the implementation. +numeric_penalty_fit <- function(source, target, row, col, balance, + step_size, iterations, weights = NULL) { + sources <- if (is.matrix(source)) list(source) else source + bases <- lapply(sources, svd, nu = 1, nv = 1) + s <- bases[[1]] + if (is.null(weights)) weights <- rep(1, length(bases)) + weights <- weights / sum(weights) + p <- nrow(sources[[1]]) + q <- ncol(sources[[1]]) + pu <- Reduce(`+`, Map(function(b, w) w * tcrossprod(b$u), bases, weights)) + pv <- Reduce(`+`, Map(function(b, w) w * tcrossprod(b$v), bases, weights)) + u <- s$u * sqrt(s$d[1]) + v <- s$v * sqrt(s$d[1]) + theta <- c(u, v) + observed <- !is.na(target) + objective <- function(x) { + U <- matrix(x[seq_len(p)], p, 1) + V <- matrix(x[p + seq_len(q)], q, 1) + residual <- (U %*% t(V) - target)[observed] + sum(residual^2) / mean(observed) + + row * sum(((diag(p) - pu) %*% U)^2) + + col * sum(((diag(q) - pv) %*% V)^2) + + balance * (sum(U^2) - sum(V^2))^2 + } + best <- theta + best_value <- objective(theta) + values <- numeric(iterations) + for (i in seq_len(iterations)) { + gradient <- vapply(seq_along(theta), function(j) { + direction <- numeric(length(theta)) + direction[j] <- 1e-6 + (objective(theta + direction) - objective(theta - direction)) / 2e-6 + }, numeric(1)) + for (indices in list(seq_len(p), p + seq_len(q))) { + theta[indices] <- theta[indices] - step_size * + sqrt(sum(theta[indices]^2)) / + (sqrt(sum(gradient[indices]^2)) + 1e-12) * gradient[indices] + } + values[i] <- objective(theta) + if (values[i] < best_value) { + best <- theta + best_value <- values[i] + } + } + list(estimate = outer(best[seq_len(p)], best[p + seq_len(q)]), + objective = values) +} + +test_that("unequal and zero penalties follow the intended objective", { + dat <- penalty_example() + for (with_missing in c(FALSE, TRUE)) { + target <- dat$target + if (with_missing) target[c(2, 7)] <- NA_real_ + for (penalties in list(c(3, 0.25), c(0.25, 3), c(0, 3), c(3, 0))) { + reference <- numeric_penalty_fit(dat$source, target, penalties[1], + penalties[2], 1, 0.02, 12) + actual <- learner(dat$source, target, r = 1, + lambda_1_row = penalties[1], + lambda_1_col = penalties[2], lambda_2 = 1, + step_size = 0.02, + control = list(max_iter = 12, threshold = 0, + max_value = 1e6)) + expect_equal(actual$objective_values, reference$objective, + tolerance = 1e-6) + expect_equal(actual$learner_estimate, reference$estimate, + tolerance = 1e-6) + } + } +}) + +test_that("invalid or ambiguous space penalties fail clearly", { + dat <- penalty_example() + fit <- function(...) learner(dat$source, dat$target, r = 1, + lambda_2 = 1, step_size = 0.01, ...) + expect_error(fit(), "Provide lambda_1") + expect_error(fit(lambda_1_row = 1), "provided together") + expect_error(fit(lambda_1_col = 1), "provided together") + expect_error(fit(lambda_1 = 1, lambda_1_row = 1, lambda_1_col = 1), + "not both modes") + expect_error(fit(lambda_1 = 1, lambda_1_row = 1), "not both modes") + invalid <- list(numeric(0), c(1, 2), NA_real_, NaN, Inf, -1, "1", TRUE) + for (value in invalid) { + expect_error(fit(lambda_1 = value), "lambda_1 must be") + expect_error(fit(lambda_1_row = value, lambda_1_col = 1), + "lambda_1_row must be") + expect_error(fit(lambda_1_row = 1, lambda_1_col = value), + "lambda_1_col must be") + } +}) + + +test_that("multiple-source updates match finite differences of the weighted objective", { + dat <- penalty_example() + set.seed(25) + sources <- list(dat$source, matrix(rnorm(20), 5, 4), matrix(rnorm(20), 5, 4)) + for (masked in c(FALSE, TRUE)) { + y <- dat$target + if (masked) y[c(2, 7)] <- NA_real_ + reference <- numeric_penalty_fit(sources, y, 3, 0.25, 1, 0.02, 12, + weights = c(0.2, 0.3, 0.5)) + actual <- learner(sources, y, r = 1, lambda_1_row = 3, lambda_1_col = 0.25, + lambda_2 = 1, step_size = 0.02, source_weights = c(2, 3, 5), + control = list(max_iter = 12, threshold = 0, max_value = 1e6)) + expect_equal(actual$objective_values, reference$objective, tolerance = 1e-6) + expect_equal(actual$learner_estimate, reference$estimate, tolerance = 1e-6) + } +}) diff --git a/tests/testthat/test-stopping-behavior.R b/tests/testthat/test-stopping-behavior.R new file mode 100644 index 0000000..7bfd930 --- /dev/null +++ b/tests/testthat/test-stopping-behavior.R @@ -0,0 +1,51 @@ +# Analytic examples distinguish a stopping code from optimization accuracy. +# The source has distinct singular values 4, 2, 1; its rank-one approximation +# and initial target estimate are exactly diag(4, 0, 0) in a 4-by-3 matrix. +stopping_example <- function() { + source <- matrix(0, 4, 3) + source[cbind(1:3, 1:3)] <- c(4, 2, 1) + target <- matrix(0, 4, 3) + target[1, 1] <- 4 + list(source = source, target = target) +} + +test_that('unchanged consecutive objectives meet a positive threshold', { + d <- stopping_example() + fit <- learner(d$source, d$target, r = 1, lambda_1 = 1, lambda_2 = 1, + step_size = 0.01, control = list(max_iter = 5, threshold = 0.001)) + expect_equal(fit$convergence_criterion, 1) + expect_equal(fit$objective_values, c(0, 0), tolerance = 1e-12) + expect_equal(fit$learner_estimate, d$target, tolerance = 1e-12) +}) + +test_that('a zero threshold uses a strict inequality and reaches the budget', { + d <- stopping_example() + fit <- learner(d$source, d$target, r = 1, lambda_1 = 1, lambda_2 = 1, + step_size = 0.01, control = list(max_iter = 5, threshold = 0)) + expect_equal(fit$convergence_criterion, 2) + expect_length(fit$objective_values, 5) + expect_equal(fit$learner_estimate, d$target, tolerance = 1e-12) +}) + +test_that('the returned estimate can retain the initialization after a bad step', { + d <- stopping_example() + # Opposite target and an intentionally huge step make the first update worse. + fit <- learner(d$source, -d$target, r = 1, lambda_1 = 1, lambda_2 = 1, + step_size = 10, control = list(max_iter = 1, threshold = 0)) + initial_objective <- sum((d$target - (-d$target))^2) + expect_equal(fit$convergence_criterion, 2) + expect_length(fit$objective_values, 1) + expect_gt(fit$objective_values[1], initial_objective) + expect_equal(fit$learner_estimate, d$target, tolerance = 1e-12) +}) + +test_that('excessive objective growth returns code three and the best estimate', { + d <- stopping_example() + fit <- learner(d$source, -d$target, r = 1, lambda_1 = 1, lambda_2 = 1, + step_size = 10, + control = list(max_iter = 5, threshold = 0, max_value = 2)) + expect_equal(fit$convergence_criterion, 3) + expect_length(fit$objective_values, 2) + expect_gt(fit$objective_values[2], 2 * fit$objective_values[1]) + expect_equal(fit$learner_estimate, d$target, tolerance = 1e-12) +})