diff --git a/DESCRIPTION b/DESCRIPTION index d0b3e4b..8169301 100644 --- a/DESCRIPTION +++ b/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")), diff --git a/R/trendfilter.R b/R/trendfilter.R index 751c6dc..b6ba28b 100644 --- a/R/trendfilter.R +++ b/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) diff --git a/src/trendfilter.cpp b/src/trendfilter.cpp index d02c268..21a5dfa 100644 --- a/src/trendfilter.cpp +++ b/src/trendfilter.cpp @@ -71,22 +71,24 @@ void admm_single_lambda(int n, const Eigen::VectorXd& y, const NumericVector& xd Eigen::Ref alpha, Eigen::Ref u, int& iter, double& obj_val, const Eigen::SparseMatrix& 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::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::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) { diff --git a/src/utils.cpp b/src/utils.cpp index 3dda48e..e36140f 100644 --- a/src/utils.cpp +++ b/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(n_lambda) - 1; diff --git a/tests/testthat/test-lambda_seq.R b/tests/testthat/test-lambda_seq.R index b0c0cc6..aeb5e6a 100644 --- a/tests/testthat/test-lambda_seq.R +++ b/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 + ) })