Add multiple sources, separate penalties, and example data

.Rbuildignore

@@ -4,3 +4,4 @@
 ^\.github$
 ^data-raw$
 ^README_files/
+^\.Rprofile$

.gitignore

@@ -2,3 +2,6 @@
 .Rhistory
 .RData
 .Ruserdata
+
+.Rprofile
+.DS_Store

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 <sean.mcgrath514@gmail.com>
 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) <doi:10.48550/arXiv.2412.20605>.
+    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

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

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() {

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"

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)
 }
-

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
+}