fix: #18 and #15

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: trendfilter
 Title: Univariate trend filtering
-Version: 0.0.2.9003
+Version: 0.0.2.9004
 Authors@R: c(
     person("Addison", "Hu", , "mail@huisaddison.com", role = c("aut", "cre")),
     person("Ryan J.", "Tibshirani", , "ryantibs@berkeley.edu", role = c("aut", "cph")),

R/trendfilter.R

@@ -128,6 +128,10 @@ trendfilter <- function(y,
   assert_numeric(lambda, finite = TRUE, lower = 0, null.ok = TRUE)
   assert_class(control, "trendfilter_control")
 
+  if (!is.null(lambda)) {
+    lambda_max <- max(lambda)
+    lambda_min <- min(lambda)
+  }
   lambda_min <- lambda_min %||% -1.0
   lambda_max <- lambda_max %||% -1.0
   lambda <- sort(lambda, decreasing = TRUE) %||% double(nlambda)

src/trendfilter.cpp

@@ -71,22 +71,24 @@ void admm_single_lambda(int n, const Eigen::VectorXd& y, const NumericVector& xd
   Eigen::Ref<Eigen::VectorXd> alpha, Eigen::Ref<Eigen::VectorXd> u, int& iter,
   double& obj_val, const Eigen::SparseMatrix<double>& dk_mat_sq,
   const Eigen::MatrixXd& denseD, const Eigen::VectorXd& s_seq, double lam,
-  int max_iter, double rho, double tol = 1e-5, int linear_solver = 2,
+  int max_iter, double rho_scale, double tol = 1e-5, int linear_solver = 2,
   bool equal_space = false) {
   // Initialize internals
   VectorXd tmp(n-k);
   VectorXd Dth_tmp(alpha.size());
   VectorXd alpha_old(alpha);
-  VectorXd wy = (y.array()*weights).matrix();
-  double rr, ss;
 
-  // LinearSystem linear_system;
-  // Technically, can form one SparseQR object, analyze the pattern once,
-  // and then re-use it.
-  // So call analyzePattern once, and then factorize repeatedly.
-  // https://eigen.tuxfamily.org/dox/classEigen_1_1SparseQR.html#aba8ae81fd3d4ce9139eccb6b7a0256b2
-  linear_system.construct(y, weights, k, rho, dk_mat_sq, denseD, s_seq, linear_solver);
-  linear_system.compute(linear_solver);
+  double rr, ss;
+  double rho = lam * rho_scale;
+  if (rho > Eigen::NumTraits<double>::epsilon()) {
+    // LinearSystem linear_system;
+    // Technically, can form one SparseQR object, analyze the pattern once,
+    // and then re-use it.
+    // So call analyzePattern once, and then factorize repeatedly.
+    // https://eigen.tuxfamily.org/dox/classEigen_1_1SparseQR.html#aba8ae81fd3d4ce9139eccb6b7a0256b2
+    linear_system.construct(y, weights, k, rho, dk_mat_sq, denseD, s_seq, linear_solver);
+    linear_system.compute(linear_solver);
+  }
 
   // Perform ADMM updates
   int computation_info;
@@ -97,15 +99,19 @@ void admm_single_lambda(int n, const Eigen::VectorXd& y, const NumericVector& xd
     if (iter % 1000 == 0) Rcpp::checkUserInterrupt(); // check if killed
 
     // theta update
-    std::tie(theta, computation_info) = linear_system.solve(y, weights,
-        alpha + u, k, xd, rho, denseD, s_seq, linear_solver, equal_space);
+    if (rho > Eigen::NumTraits<double>::epsilon()) {
+      std::tie(theta, computation_info) = linear_system.solve(y, weights,
+          alpha + u, k, xd, rho, denseD, s_seq, linear_solver, equal_space);
+    } else {
+      theta = y;
+    }
     // if (computation_info > 1) {
     //  std::cerr << "Eigen Sparse QR solve returned nonzero exit status.\n";
     // }
     Dth_tmp = Dkv(theta, k, xd);
     tmp = Dth_tmp - u;
     // alpha update
-    alpha = tf_dp(tmp, lam / rho);
+    alpha = tf_dp(tmp, 1 / rho_scale);
     // u update
     u += alpha - Dth_tmp;
     // double cur_objective = tf_objective(y, theta, xd, weights, lam, k);
@@ -139,7 +145,8 @@ Rcpp::List admm_lambda_seq(
 
   int n = x.size();
 
-  if (lambda[0] < tol / 100 && lambda_max <= 0) {
+
+  if (lambda_max < -tol) {
     lambda_max = get_lambda_max(x, y, weights, k);
   }
   get_lambda_seq(lambda, lambda_max, lambda_min, lambda_min_ratio, nlambda);
@@ -190,20 +197,23 @@ Rcpp::List admm_lambda_seq(
   }
 
   Eigen::MatrixXd alpha(n-k, nlambda);
+  VectorXd u(n-k);
 
   // Initialize ADMM variables
   // Project onto Legendre polynomials to initialize for largest lambda.
-  theta.col(0) = project_polynomials(x, y, weights, k);
-  alpha.col(0) = Dkv(theta.col(0), k, x);
-  VectorXd u = init_u((theta.col(0) - y)/(lambda[0]*rho_scale), x, k, weights);
 
+  if (lambda_max > 0) {
+    theta.col(0) = project_polynomials(x, y, weights, std::fmin(k, 3));
+    alpha.col(0) = Dkv(theta.col(0), k, x);
+    u = init_u((theta.col(0) - y) / (lambda_max*rho_scale), x, k, weights);
+  }
 
   for (int i = 0; i < nlambda; i++) {
     Rcpp::checkUserInterrupt();
     admm_single_lambda(n, y, x, weights, k,
       theta.col(i), alpha.col(i), u,
       iters[i], objective_val[i],
-      dk_mat_sq, denseD, s_seq, lambda[i], max_iter, lambda[i]*rho_scale,
+      dk_mat_sq, denseD, s_seq, lambda[i], max_iter, rho_scale,
       tol, linear_solver, equal_space);
     dof[i] = calc_degrees_of_freedom(alpha.col(i), k);
     if (i + 1 < nlambda) {

src/utils.cpp

@@ -224,11 +224,7 @@ void get_lambda_seq(
     double lambda_min_ratio = 1e-5,
     int n_lambda = 50) {
 
-  if (!(lambda.array() < 1e-12).all()) {
-    lambda_min = lambda.minCoeff();
-    lambda_max = lambda.maxCoeff();
-    n_lambda = lambda.size();
-  } else {
+  if (lambda_max > 0 && n_lambda > 1) {
     double lmpad = lambda_min_ratio * lambda_max;
     lambda_min = (lambda_min < 0) ? lmpad : lambda_min;
     double ns = static_cast<double>(n_lambda) - 1;

tests/testthat/test-lambda_seq.R

@@ -28,4 +28,10 @@ test_that("get_lambda_seq works", {
     get_lambda_seq_r(double(10), 10, 0, 1e-4, 10),
     c(10^seq(log10(10), log10(1e-3), length.out = 9), 0)
   )
+
+  # allow 0 for lambda glmgen/trendfilter#18
+  expect_equal(
+    get_lambda_seq_r(double(1), -1, -1, 1e-4, 1),
+    0
+  )
 })