diff --git a/DESCRIPTION b/DESCRIPTION
index 8169301..79264b1 100644
--- a/DESCRIPTION
+++ b/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
diff --git a/NEWS.md b/NEWS.md
index a719ff7..66910b4 100644
--- a/NEWS.md
+++ b/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
 
diff --git a/R/control-lists.R b/R/control-lists.R
index eebce4a..962c29b 100644
--- a/R/control-lists.R
+++ b/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
diff --git a/R/trendfilter.R b/R/trendfilter.R
index b6ba28b..9549a5c 100644
--- a/R/trendfilter.R
+++ b/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"
+  )
 }
diff --git a/man/admm_control_list.Rd b/man/admm_control_list.Rd
index fc9c1a2..2392f9b 100644
--- a/man/admm_control_list.Rd
+++ b/man/admm_control_list.Rd
@@ -8,7 +8,7 @@ admm_control_list(
   max_iter = 10000,
   rho_scale = 1,
   tolerance = 1e-04,
-  linear_solver = c("kalman_filter", "sparse_qr"),
+  linear_solver = c("kalman_filter", "sparse_qr", "sparse_cholesky"),
   space_tolerance_ratio = sqrt(.Machine$double.eps),
   ...
 )
@@ -22,9 +22,11 @@ Lagrangian penalty parameter \eqn{\rho}.}
 \item{tolerance}{Double. The convergence tolerance for the ADMM algorithm.}
 
 \item{linear_solver}{Integer. Solver for the linear system in ADMM when
-k > 1: \code{kalman_filter} for Kalman filter or \code{sparse_qr} for sparse QR decomposition.}
+k > 1: \code{kalman_filter} for Kalman filter, \code{sparse_qr} for sparse QR
+decomposition, or \code{sparse_cholesky} for Cholesky.}
 
-\item{space_tolerance_ratio}{Double. The tolerance ratio of to detect that the signal is equally spaced.}
+\item{space_tolerance_ratio}{Double. The tolerance ratio of to detect that
+the signal is equally spaced.}
 
 \item{...}{not used}
 }
diff --git a/src/linearsystem.cpp b/src/linearsystem.cpp
index 5377d37..34cf49b 100644
--- a/src/linearsystem.cpp
+++ b/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);
 }
diff --git a/src/linearsystem.h b/src/linearsystem.h
index e7bcdfa..4eb9d49 100644
--- a/src/linearsystem.h
+++ b/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
diff --git a/tests/testthat/test-linear_solver.R b/tests/testthat/test-linear_solver.R
index 4e84239..617023d 100644
--- a/tests/testthat/test-linear_solver.R
+++ b/tests/testthat/test-linear_solver.R
@@ -6,20 +6,104 @@ k1 <- 2L
 # signal locations
 x1 <- 1:n1 / (n1 + 1) * 2 * pi
 # random standard-normally distributed signals with size n1
-y0 <- c(0.006, -0.119, 0.262, 0.043, 1.114, -1.020, -0.350, -0.889, 0.476,
-        0.865, 0.645, -0.631, 0.277, -0.480, 0.364, 0.599, 0.608, 1.755,
-        -0.041, -0.124, -0.031, 0.195, -0.093, -0.207, 2.405, -0.026,
-        -1.115, 0.474, 0.129, 0.222)
+y0 <- c(
+  0.006,
+  -0.119,
+  0.262,
+  0.043,
+  1.114,
+  -1.020,
+  -0.350,
+  -0.889,
+  0.476,
+  0.865,
+  0.645,
+  -0.631,
+  0.277,
+  -0.480,
+  0.364,
+  0.599,
+  0.608,
+  1.755,
+  -0.041,
+  -0.124,
+  -0.031,
+  0.195,
+  -0.093,
+  -0.207,
+  2.405,
+  -0.026,
+  -1.115,
+  0.474,
+  0.129,
+  0.222
+)
 # transform to sinusoidal signals with white noise
 y1 <- sin(x1) + .2 * y0
 # random Poisson signals with mean 5 for arbitrary weights on signals
-weig1 <- c(5, 5, 1, 5, 8, 3, 2, 6, 4, 5, 4, 3, 4, 6, 5, 2, 9, 4, 1, 2, 5, 4, 3,
-           5, 3, 3, 9, 4, 6, 7)
+weig1 <- c(
+  5,
+  5,
+  1,
+  5,
+  8,
+  3,
+  2,
+  6,
+  4,
+  5,
+  4,
+  3,
+  4,
+  6,
+  5,
+  2,
+  9,
+  4,
+  1,
+  2,
+  5,
+  4,
+  3,
+  5,
+  3,
+  3,
+  9,
+  4,
+  6,
+  7
+)
 # random standard-normally distributed signals as mean adjustment
-mn1 <- c(0.485, -0.047, -0.200, 1.646, -0.594, 1.167, -0.816, 0.295, 0.784,
-         1.446, 1.328, 0.852, -0.571, 1.007, -0.741, 1.114, 0.232, -1.886,
-         -0.542, 2.111, -0.222, 0.507, -0.111, 0.199, 0.181, -0.363, 0.505,
-         -0.868)
+mn1 <- c(
+  0.485,
+  -0.047,
+  -0.200,
+  1.646,
+  -0.594,
+  1.167,
+  -0.816,
+  0.295,
+  0.784,
+  1.446,
+  1.328,
+  0.852,
+  -0.571,
+  1.007,
+  -0.741,
+  1.114,
+  0.232,
+  -1.886,
+  -0.542,
+  2.111,
+  -0.222,
+  0.507,
+  -0.111,
+  0.199,
+  0.181,
+  -0.363,
+  0.505,
+  -0.868
+)
 
 # Scenario 2: n = 50, k = 3 ----
 n2 <- 50L
@@ -27,24 +111,162 @@ k2 <- 3L
 # signal locations
 x2 <- 1:n2 / (n2 + 1) * 2 * pi
 # random signals from Normal(0,1) with size n2
-y0 <- c(-0.065, 0.654, 1.397, 0.626, 1.233, 0.818, 0.586, -0.693, -0.345,
-        -0.550, 0.618, -0.216, -1.355, 0.131, -0.461, 0.581, -1.709, -0.881,
-        0.506, 1.749, -0.449, 1.856, -0.840, -2.980, -0.355, -1.221, 0.450,
-        0.006, 0.350, 0.790, 0.405, 1.467, 0.890, 0.282, -1.038, 0.208, -0.039,
-        1.166, 1.192, 0.027, 1.078, 2.363, -0.080, -0.074, -0.638, -0.539,
-        -1.721, 0.646, -0.899, 0.600)
+y0 <- c(
+  -0.065,
+  0.654,
+  1.397,
+  0.626,
+  1.233,
+  0.818,
+  0.586,
+  -0.693,
+  -0.345,
+  -0.550,
+  0.618,
+  -0.216,
+  -1.355,
+  0.131,
+  -0.461,
+  0.581,
+  -1.709,
+  -0.881,
+  0.506,
+  1.749,
+  -0.449,
+  1.856,
+  -0.840,
+  -2.980,
+  -0.355,
+  -1.221,
+  0.450,
+  0.006,
+  0.350,
+  0.790,
+  0.405,
+  1.467,
+  0.890,
+  0.282,
+  -1.038,
+  0.208,
+  -0.039,
+  1.166,
+  1.192,
+  0.027,
+  1.078,
+  2.363,
+  -0.080,
+  -0.074,
+  -0.638,
+  -0.539,
+  -1.721,
+  0.646,
+  -0.899,
+  0.600
+)
 y2 <- sin(x2) + .2 * y0
 # random signals from Poisson(5) for arbitrary weights on signals
-weig2 <- c(10, 2, 7, 9, 6, 0, 4, 5, 4, 6, 6, 3, 7, 3, 2, 4, 3, 4, 4, 7, 4, 3,
-           5, 4, 5, 7, 7, 5, 4, 5, 7, 2, 6, 3, 6, 1, 7, 6, 7, 6, 4, 7, 2, 3,
-           5, 9, 6, 3, 4, 3)
+weig2 <- c(
+  10,
+  2,
+  7,
+  9,
+  6,
+  0,
+  4,
+  5,
+  4,
+  6,
+  6,
+  3,
+  7,
+  3,
+  2,
+  4,
+  3,
+  4,
+  4,
+  7,
+  4,
+  3,
+  5,
+  4,
+  5,
+  7,
+  7,
+  5,
+  4,
+  5,
+  7,
+  2,
+  6,
+  3,
+  6,
+  1,
+  7,
+  6,
+  7,
+  6,
+  4,
+  7,
+  2,
+  3,
+  5,
+  9,
+  6,
+  3,
+  4,
+  3
+)
 # random signals from Normal(0,1) as mean adjustment
-mn2 <- c(0.194, -0.579, 0.240, -1.641, -0.009, 0.865, -0.362, -0.739, 0.463,
-        0.854, 0.103, -0.718, -0.110, 0.367, 0.131, 0.368, -1.452, 2.719,
-        -0.314, 0.846, 0.017, -0.899, -0.569, 0.116, -0.849, 0.969, -0.614,
-        1.181, 1.432, 0.369, 0.427, 0.200, 0.218, -0.678, -0.513, -1.233,
-        -1.161, -0.601, 1.064, -3.186, -1.311, 0.603, -0.092, -0.429, 0.997,
-        -0.419, -2.362)
+mn2 <- c(
+  0.194,
+  -0.579,
+  0.240,
+  -1.641,
+  -0.009,
+  0.865,
+  -0.362,
+  -0.739,
+  0.463,
+  0.854,
+  0.103,
+  -0.718,
+  -0.110,
+  0.367,
+  0.131,
+  0.368,
+  -1.452,
+  2.719,
+  -0.314,
+  0.846,
+  0.017,
+  -0.899,
+  -0.569,
+  0.116,
+  -0.849,
+  0.969,
+  -0.614,
+  1.181,
+  1.432,
+  0.369,
+  0.427,
+  0.200,
+  0.218,
+  -0.678,
+  -0.513,
+  -1.233,
+  -1.161,
+  -0.601,
+  1.064,
+  -3.186,
+  -1.311,
+  0.603,
+  -0.092,
+  -0.429,
+  0.997,
+  -0.419,
+  -2.362
+)
 
 test_that("test a single iterate of linear system solvers yields same results", {
   rho <- 2
@@ -54,44 +276,94 @@ test_that("test a single iterate of linear system solvers yields same results",
     Diagonal(n1, weig1) + rho * Matrix::crossprod(Dkx),
     weig1 * y1 + rho * Matrix::crossprod(Dkx, mn1)
   ))
-  theta_sparseQR <- linear_single_solve_test(1, y1, weig1, x1, rho, mn1)
+  theta_sparse_qr <- linear_single_solve_test(1, y1, weig1, x1, rho, mn1)
   theta_kf <- linear_single_solve_test(2, y1, weig1, x1, rho, mn1)
+  theta_cholesky <- linear_single_solve_test(3, y1, weig1, x1, rho, mn1)
   expect_equal(theta_kf, theta)
-  expect_equal(theta_sparseQR, theta)
+  expect_equal(theta_sparse_qr, theta)
+  expect_equal(theta_cholesky, theta)
   # Scenario 2
   Dkx <- dspline::d_mat(k2, x2, FALSE)
   theta <- as.vector(Matrix::solve(
     Diagonal(n2, weig2) + rho * Matrix::crossprod(Dkx),
     weig2 * y2 + rho * Matrix::crossprod(Dkx, mn2)
   ))
-  theta_sparseQR <- linear_single_solve_test(1, y2, weig2, x2, rho, mn2)
+  theta_sparse_qr <- linear_single_solve_test(1, y2, weig2, x2, rho, mn2)
   theta_kf <- linear_single_solve_test(2, y2, weig2, x2, rho, mn2)
+  theta_cholesky <- linear_single_solve_test(3, y2, weig2, x2, rho, mn2)
   expect_equal(theta_kf, theta)
-  expect_equal(theta_sparseQR, theta)
+  expect_equal(theta_sparse_qr, theta)
+  expect_equal(theta_cholesky, theta)
 })
 
 test_that("test linear solvers yield same estimates for single lambda", {
   lam <- 5
   # Scenario 1
   mod_kf <- trendfilter(
-    y1, x1, weights = weig1, k = k1, lambda = lam,
+    y1,
+    x1,
+    weights = weig1,
+    k = k1,
+    lambda = lam,
     control = trendfilter_control_list(
-      admm_control = admm_control_list(linear_solver = "kalman_filter")))
-  mod_sparseQR <- trendfilter(
-    y1, x1, weights = weig1, k = k1, lambda = lam,
+      admm_control = admm_control_list(linear_solver = "kalman_filter")
+    )
+  )
+  mod_sparse_qr <- trendfilter(
+    y1,
+    x1,
+    weights = weig1,
+    k = k1,
+    lambda = lam,
     control = trendfilter_control_list(
-      admm_control = admm_control_list(linear_solver = "sparse_qr")))
-  expect_equal(mod_kf$theta, mod_sparseQR$theta)
+      admm_control = admm_control_list(linear_solver = "sparse_qr")
+    )
+  )
+  mod_cholesky <- trendfilter(
+    y1,
+    x1,
+    weights = weig1,
+    k = k1,
+    lambda = lam,
+    control = trendfilter_control_list(
+      admm_control = admm_control_list(linear_solver = "sparse_cholesky")
+    )
+  )
+  expect_equal(mod_kf$theta, mod_sparse_qr$theta)
+  expect_equal(mod_kf$theta, mod_cholesky$theta)
   # Scenario 2
   mod_kf <- trendfilter(
-    y2, x2, weights = weig2, k = k2, lambda = lam,
+    y2,
+    x2,
+    weights = weig2,
+    k = k2,
+    lambda = lam,
     control = trendfilter_control_list(
-      admm_control = admm_control_list(linear_solver = "kalman_filter")))
+      admm_control = admm_control_list(linear_solver = "kalman_filter")
+    )
+  )
   mod_sparse_qr <- trendfilter(
-    y2, x2, weights = weig2, lambda = lam, k = k2,
+    y2,
+    x2,
+    weights = weig2,
+    lambda = lam,
+    k = k2,
+    control = trendfilter_control_list(
+      admm_control = admm_control_list(linear_solver = "sparse_qr")
+    )
+  )
+  mod_cholesky <- trendfilter(
+    y2,
+    x2,
+    weights = weig2,
+    lambda = lam,
+    k = k2,
     control = trendfilter_control_list(
-      admm_control = admm_control_list(linear_solver = "sparse_qr")))
+      admm_control = admm_control_list(linear_solver = "sparse_cholesky")
+    )
+  )
   expect_equal(mod_kf$theta, mod_sparse_qr$theta)
+  expect_equal(mod_kf$theta, mod_cholesky$theta)
 })
 
 test_that("test linear solvers run with no errors for lambda sequences", {
@@ -99,49 +371,116 @@ test_that("test linear solvers run with no errors for lambda sequences", {
   # Scenario 1
   expect_no_error(
     mod_kf <- trendfilter(
-      y1, x1, weights = weig1, k = k1, nlambda = nlam,
+      y1,
+      x1,
+      weights = weig1,
+      k = k1,
+      nlambda = nlam,
       control = trendfilter_control_list(
-        admm_control = admm_control_list(linear_solver = "kalman_filter"))
+        admm_control = admm_control_list(linear_solver = "kalman_filter")
+      )
     )
   )
   expect_no_error(
     mod_sparse_qr <- trendfilter(
-      y1, x1, weights = weig1, k = k1, nlambda = nlam,
+      y1,
+      x1,
+      weights = weig1,
+      k = k1,
+      nlambda = nlam,
       control = trendfilter_control_list(
-        admm_control = admm_control_list(linear_solver = "sparse_qr"))
+        admm_control = admm_control_list(linear_solver = "sparse_qr")
+      )
+    )
+  )
+  expect_no_error(
+    mod_cholesky <- trendfilter(
+      y1,
+      x1,
+      weights = weig1,
+      k = k1,
+      nlambda = nlam,
+      control = trendfilter_control_list(
+        admm_control = admm_control_list(linear_solver = "sparse_cholesky")
+      )
     )
   )
   # equal estimates from the first model
-  expect_equal(mod_kf$theta[,1], mod_sparse_qr$theta[,1])
+  expect_equal(mod_kf$theta[, 1], mod_sparse_qr$theta[, 1])
+  expect_equal(mod_kf$theta[, 1], mod_cholesky$theta[, 1])
   # Scenario 2
-  skip("The rest of this takes too long to run.")
+  # skip("The rest of this takes too long to run.")
   expect_no_error(
     mod_kf <- trendfilter(
-      y2, x2, weights = weig2, k = k2, nlambda = 2,
+      y2,
+      x2,
+      weights = weig2,
+      k = k2,
+      nlambda = 2,
       control = trendfilter_control_list(
-        admm_control = admm_control_list(linear_solver = "kalman_filter"))
+        admm_control = admm_control_list(linear_solver = "kalman_filter")
+      )
     )
   )
   expect_no_error(
     mod_sparse_qr <- trendfilter(
-      y2, x2, weights = weig2, k = k2, nlambda = 2,
+      y2,
+      x2,
+      weights = weig2,
+      k = k2,
+      nlambda = 2,
       control = trendfilter_control_list(
-        admm_control = admm_control_list(linear_solver = "sparse_qr"))
+        admm_control = admm_control_list(linear_solver = "sparse_qr")
+      )
     )
   )
-  # equal estimates from the first model
-  expect_equal(mod_kf$theta[,1], mod_sparse_qr$theta[,1])
+  expect_no_error(
+    mod_cholesky <- trendfilter(
+      y2,
+      x2,
+      weights = weig2,
+      k = k2,
+      nlambda = 2,
+      control = trendfilter_control_list(
+        admm_control = admm_control_list(linear_solver = "sparse_cholesky")
+      )
+    )
+  )
+  expect_equal(mod_kf$theta[, 1], mod_sparse_qr$theta[, 1])
+  expect_equal(mod_kf$theta[, 1], mod_cholesky$theta[, 1])
 })
 
 test_that("test marginal cases", {
   x <- 1:n1 / (n1 + 1) * 2 * pi
   w0 <- rep(0, n1) # zero weights on all signals
   expect_no_error(trendfilter(
-    y = rep(1, n1), x, weights = w0, k = k1, lambda = 10,
+    y = rep(1, n1),
+    x,
+    weights = w0,
+    k = k1,
+    lambda = 10,
     control = trendfilter_control_list(
-      admm_control = admm_control_list(linear_solver = "kalman_filter"))))
+      admm_control = admm_control_list(linear_solver = "kalman_filter")
+    )
+  ))
+  expect_no_error(trendfilter(
+    y = rep(1, n1),
+    x,
+    weights = w0,
+    k = k1,
+    lambda = 10,
+    control = trendfilter_control_list(
+      admm_control = admm_control_list(linear_solver = "sparse_qr")
+    )
+  ))
   expect_no_error(trendfilter(
-    y = rep(1, n1), x, weights = w0, k = k1, lambda = 10,
+    y = rep(1, n1),
+    x,
+    weights = w0,
+    k = k1,
+    lambda = 10,
     control = trendfilter_control_list(
-      admm_control = admm_control_list(linear_solver = "sparse_qr"))))
+      admm_control = admm_control_list(linear_solver = "sparse_cholesky")
+    )
+  ))
 })