add option to use sparse cholesky as a solver

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: trendfilter
 Title: Univariate trend filtering
-Version: 0.0.2.9004
+Version: 0.0.2.9005
 Authors@R: c(
     person("Addison", "Hu", , "mail@huisaddison.com", role = c("aut", "cre")),
     person("Ryan J.", "Tibshirani", , "ryantibs@berkeley.edu", role = c("aut", "cph")),
@@ -30,4 +30,4 @@ LinkingTo:
 Config/testthat/edition: 3
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3

NEWS.md

@@ -1,6 +1,7 @@
 # trendfilter (development version)
 
 * Initial CRAN submission.
+* add option to ues `sparse_cholesky` as a solver.
 
 # trendfilter 0.0.2
 

R/control-lists.R

@@ -8,8 +8,10 @@
 #'   Lagrangian penalty parameter \eqn{\rho}.
 #' @param tolerance Double. The convergence tolerance for the ADMM algorithm.
 #' @param linear_solver Integer. Solver for the linear system in ADMM when
-#'   k > 1: `kalman_filter` for Kalman filter or `sparse_qr` for sparse QR decomposition.
-#' @param space_tolerance_ratio Double. The tolerance ratio of to detect that the signal is equally spaced.
+#'   k > 1: `kalman_filter` for Kalman filter, `sparse_qr` for sparse QR
+#'   decomposition, or `sparse_cholesky` for Cholesky.
+#' @param space_tolerance_ratio Double. The tolerance ratio of to detect that
+#'   the signal is equally spaced.
 #' @param ... not used
 #'
 #' @return an object of class `admm_control`
@@ -20,17 +22,39 @@
 #' admm_control_list(max_iter = 10L)
 #' admm_control_list(tolerance = 1e-8)
 admm_control_list <- function(
-    max_iter = 1e4, rho_scale = 1.0, tolerance = 1e-4,
-    linear_solver = c("kalman_filter", "sparse_qr"),
-    space_tolerance_ratio = sqrt(.Machine$double.eps), ...) {
+  max_iter = 1e4,
+  rho_scale = 1.0,
+  tolerance = 1e-4,
+  linear_solver = c("kalman_filter", "sparse_qr", "sparse_cholesky"),
+  space_tolerance_ratio = sqrt(.Machine$double.eps),
+  ...
+) {
   rlang::check_dots_empty()
   assert_integerish(max_iter, lower = 1L, len = 1L)
-  assert_numeric(rho_scale, lower = .Machine$double.eps, finite = TRUE, len = 1L)
-  assert_numeric(tolerance, lower = .Machine$double.eps, finite = TRUE, len = 1L)
+  assert_numeric(
+    rho_scale,
+    lower = .Machine$double.eps,
+    finite = TRUE,
+    len = 1L
+  )
+  assert_numeric(
+    tolerance,
+    lower = .Machine$double.eps,
+    finite = TRUE,
+    len = 1L
+  )
   linear_solver <- rlang::arg_match(linear_solver)
   assert_numeric(space_tolerance_ratio, lower = 0, finite = TRUE, len = 1L)
-  structure(enlist(max_iter, rho_scale, tolerance, linear_solver, space_tolerance_ratio),
-            class = "admm_control")
+  structure(
+    enlist(
+      max_iter,
+      rho_scale,
+      tolerance,
+      linear_solver,
+      space_tolerance_ratio
+    ),
+    class = "admm_control"
+  )
 }
 
 #' @export
@@ -61,12 +85,19 @@ print.admm_control <- function(x, prefix = "An", ...) {
 #' trendfilter_control_list(obj_tol = 1e-12)
 #' trendfilter_control_list(admm_control = admm_control_list(tolerance = 1e-5))
 trendfilter_control_list <- function(
-    obj_tol = 1e-6, x_cond = 1e11, admm_control = admm_control_list(), ...) {
+  obj_tol = 1e-6,
+  x_cond = 1e11,
+  admm_control = admm_control_list(),
+  ...
+) {
   rlang::check_dots_empty()
   assert_numeric(obj_tol, lower = .Machine$double.eps, finite = TRUE, len = 1L)
   assert_numeric(x_cond, lower = 1, finite = TRUE, len = 1L)
   assert_class(admm_control, "admm_control")
-  structure(enlist(obj_tol, x_cond, admm_control), class = "trendfilter_control")
+  structure(
+    enlist(obj_tol, x_cond, admm_control),
+    class = "trendfilter_control"
+  )
 }
 
 #' @export

R/trendfilter.R

@@ -89,19 +89,21 @@
 #' out <- trendfilter(y, x)
 #'
 #' plot(out)
-trendfilter <- function(y,
-                        x = seq_along(y),
-                        weights = rep(1, n),
-                        k = 3L,
-                        family = c("gaussian", "logistic", "poisson"),
-                        method = c("admm", "pdip", "hybrid"),
-                        lambda = NULL,
-                        nlambda = 50L,
-                        lambda_max = NULL,
-                        lambda_min = NULL,
-                        lambda_min_ratio = 1e-5,
-                        standardize = TRUE,
-                        control = trendfilter_control_list()) {
+trendfilter <- function(
+  y,
+  x = seq_along(y),
+  weights = rep(1, n),
+  k = 3L,
+  family = c("gaussian", "logistic", "poisson"),
+  method = c("admm", "pdip", "hybrid"),
+  lambda = NULL,
+  nlambda = 50L,
+  lambda_max = NULL,
+  lambda_min = NULL,
+  lambda_min_ratio = 1e-5,
+  standardize = TRUE,
+  control = trendfilter_control_list()
+) {
   family <- arg_match(family)
   if (family != "gaussian") {
     cli_abort("Data family {.val {family}} is not yet implemented.")
@@ -116,13 +118,20 @@ trendfilter <- function(y,
   assert_numeric(weights, lower = 0, finite = TRUE, len = n)
   assert_integerish(k, lower = 0L, upper = n - 1L, len = 1L)
   assert_integerish(nlambda, lower = 1L, len = 1L)
-  assert_numeric(lambda_max,
-    len = 1L, lower = lambda_min %||% 0, finite = TRUE,
+  assert_numeric(
+    lambda_max,
+    len = 1L,
+    lower = lambda_min %||% 0,
+    finite = TRUE,
     null.ok = TRUE
   )
-  assert_numeric(lambda_min,
-    len = 1L, lower = 0, upper = lambda_max %||% Inf,
-    finite = TRUE, null.ok = TRUE
+  assert_numeric(
+    lambda_min,
+    len = 1L,
+    lower = 0,
+    upper = lambda_max %||% Inf,
+    finite = TRUE,
+    null.ok = TRUE
   )
   assert_numeric(lambda_min_ratio, lower = 0, upper = 1, len = 1L)
   assert_numeric(lambda, finite = TRUE, lower = 0, null.ok = TRUE)
@@ -154,29 +163,46 @@ trendfilter <- function(y,
     y <- (y - ym) / ys
   }
 
+  linear_solver_int <- match(
+    control$admm_control$linear_solver,
+    c("sparse_qr", "kalman_filter", "sparse_cholesky")
+  )
   out <- admm_lambda_seq(
-    xsc, y, wsc, k,
-    lambda, nlambda, lambda_max, lambda_min, lambda_min_ratio,
-    control$admm_control$max_iter, control$admm_control$rho_scale,
+    xsc,
+    y,
+    wsc,
+    k,
+    lambda,
+    nlambda,
+    lambda_max,
+    lambda_min,
+    lambda_min_ratio,
+    control$admm_control$max_iter,
+    control$admm_control$rho_scale,
     control$admm_control$tolerance,
-    if (k == 1L) 0L else match(control$admm_control$linear_solver, c("sparse_qr", "kalman_filter")),
+    if (k == 1L) 0L else linear_solver_int,
     control$admm_control$space_tolerance_ratio
   )
 
   alpha <- NULL
-  if (!is.null(out$alpha)) alpha <- drop(out$alpha) * ys
+  if (!is.null(out$alpha)) {
+    alpha <- drop(out$alpha) * ys
+  }
 
-  structure(enlist(
-    y = y * ys + ym,
-    x,
-    weights,
-    k,
-    theta = drop(out$theta) * ys + ym,
-    alpha = alpha,
-    lambda = out$lambda,
-    iters = out$iters,
-    objective = out$tf_objective,
-    dof = out$dof,
-    call = match.call()
-  ), class = "trendfilter")
+  structure(
+    enlist(
+      y = y * ys + ym,
+      x,
+      weights,
+      k,
+      theta = drop(out$theta) * ys + ym,
+      alpha = alpha,
+      lambda = out$lambda,
+      iters = out$iters,
+      objective = out$tf_objective,
+      dof = out$dof,
+      call = match.call()
+    ),
+    class = "trendfilter"
+  )
 }

src/linearsystem.cpp

@@ -26,8 +26,9 @@ void LinearSystem::construct(const Eigen::VectorXd& y, const Eigen::ArrayXd& wei
     int k, double rho, const Eigen::SparseMatrix<double>& dk_mat_sq, 
     const Eigen::MatrixXd& Dseq, const Eigen::VectorXd& s_seq, int solver) {
   switch(solver) {
-    case 0: 
-    case 1: {
+    case 0: [[fallthrough]];
+    case 1: [[fallthrough]];
+    case 3: {
       wy = (y.array() * weights).matrix();
       // Form Gram matrix and set up linear system for theta update
       A = rho * dk_mat_sq;
@@ -61,6 +62,10 @@ void LinearSystem::compute(int solver) {
     }
     case 2:
       break;
+    case 3: {
+      cholesky.compute(A);
+      break;
+    }
   }
 }
 
@@ -85,6 +90,12 @@ std::tuple<VectorXd,int> LinearSystem::solve(const Eigen::VectorXd& y,
       LinearSystem::kf_iter(y, weights, adj_mean, Dseq, s_seq, equal_space);
       break;
     }
+    case 3: {
+      VectorXd v = wy + rho * Dktv(adj_mean, k, x);
+      sol = cholesky.solve(v);
+      info = int(cholesky.info());
+      break;
+    }
   }
   return std::make_tuple(sol, info);
 }

src/linearsystem.h

@@ -38,6 +38,8 @@ class LinearSystem {
     double vt_b, Ft_b, Finf_b;
     Eigen::VectorXd RQR, a1, vt, Ft, Finf, Kt_b, Kinf_b, r,  r1, rtmp, sol;
     Eigen::MatrixXd T, at, P1, Pt, P1inf, Pinf, Kt, Kinf, L0, L1, Ptemp;
+    // cholesky
+    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> cholesky;
 };
 
 #endif