add option to directly use the initial estimates for every lambda value

R/pense_regression.R

@@ -28,11 +28,14 @@
 #' @param mscale_opts options for the M-scale estimation. See [mscale_algorithm_options] for details.
 #' @param enpy_opts options for the ENPY initial estimates, created with the
 #'                  [enpy_options] function. See [enpy_initial_estimates] for details.
+#' @param zero_based also compute the PENSE regularization path starting from the 0 vector.
+#' @param share_initials use all initial estimates for every lambda value.
 #'
 #' @export
 pense <- function(x, y, alpha, lambdas, cold_lambdas, penalty_loadings, additional_initial_estimates,
                   include_intercept = TRUE, bdp = 0.25, cc, eps = 1e-6, algorithm_opts, nr_tracks = 10, explore_it = 10,
-                  sparse = TRUE, mscale_opts = mscale_algorithm_options(), enpy_opts = enpy_options()) {
+                  sparse = TRUE, mscale_opts = mscale_algorithm_options(), enpy_opts = enpy_options(),
+                  zero_based = TRUE, share_initials = FALSE) {
   optional_args <- list()
 
   # Normalize input
@@ -63,12 +66,12 @@ pense <- function(x, y, alpha, lambdas, cold_lambdas, penalty_loadings, addition
                      eps = .as(eps[[1L]], 'numeric'),
                      explore_it = .as(explore_it[[1L]], 'integer'),
                      nr_tracks = .as(nr_tracks[[1L]], 'integer'),
+                     compute_0_based = isTRUE(zero_based),
+                     same_initials_for_all = isTRUE(share_initials),
                      mscale = .full_mscale_algo_options(bdp = bdp, cc = cc, mscale_opts = mscale_opts))
 
   # Propagate `sparse` to the inner optimizer
-  if (!is.null(pense_opts$algo_opts$sparse)) {
-    pense_opts$algo_opts$sparse <- sparse
-  }
+  pense_opts$algo_opts$sparse <- sparse
 
   # If using the MM algorithm, ensure that the EN options are set.
   if (pense_opts$algorithm == 1L) {
@@ -113,7 +116,13 @@ pense <- function(x, y, alpha, lambdas, cold_lambdas, penalty_loadings, addition
   }
 
   init_ests <- if (!inherits(additional_initial_estimates, 'pense_estimate_list')) {
-    .make_initest_list(additional_initial_estimates, lambdas, sparse = sparse)
+    if (pense_opts$same_initials_for_all) {
+      additional_initial_estimates
+    } else {
+      .make_initest_list(additional_initial_estimates, lambdas, sparse = sparse)
+    }
+  } else if (pense_opts$same_initials_for_all) {
+    unlist(additional_initial_estimates, recursive = FALSE)
   } else {
     additional_initial_estimates
   }

R/utilities.R

@@ -50,7 +50,7 @@ mloc <- function (x, scale = mad(x), rho, cc, opts = mscale_algorithm_options())
   if (missing(cc)) {
     cc <- NULL
   }
-  opts <- .full_mscale_algo_options(bdp, cc, opts)
+  opts <- .full_mscale_algo_options(.5, cc, opts)
   opts$rho <- rho_function(rho)
   .Call(C_mloc, .as(x, 'numeric'), scale, opts)
 }
@@ -70,10 +70,6 @@ mloc <- function (x, scale = mad(x), rho, cc, opts = mscale_algorithm_options())
 mlocscale <- function (x, bdp = 0.25, location_rho = c('bisquare', 'huber'),
                        scale_cc = consistency_const(bdp, 'bisquare'),
                        location_cc, opts = mscale_algorithm_options()) {
-  # No checks for NA values!
-  if (missing(cc)) {
-    cc <- NULL
-  }
   opts <- .full_mscale_algo_options(bdp, scale_cc, opts)
   loc_opts <- list(rho = rho_function(location_rho))
   if (!missing(location_cc)) {

man/en_admm_options.Rd

@@ -4,7 +4,7 @@
 \alias{en_admm_options}
 \title{Options for the ADMM Elastic Net Algorithm}
 \usage{
-en_admm_options(max_it = 1000, eps = 1e-09, tau, sparse = FALSE,
+en_admm_options(max_it = 1000, eps = 1e-08, tau, sparse = FALSE,
   admm_type = c("auto", "linearized", "var-stepsize"),
   tau_lower_mult = 0.01, tau_adjustment_lower = 0.98,
   tau_adjustment_upper = 0.999)

man/mloc.Rd

@@ -4,7 +4,7 @@
 \alias{mloc}
 \title{Compute the M-estimate of Location}
 \usage{
-mloc(x, scale = mad(x), rho, cc, opts = mest_options())
+mloc(x, scale = mad(x), rho, cc, opts = mscale_algorithm_options())
 }
 \arguments{
 \item{x}{numeric values.}
@@ -17,8 +17,8 @@ mloc(x, scale = mad(x), rho, cc, opts = mest_options())
 By default, chosen to achieve 95% efficiency under the Normal
 model.}
 
-\item{opts}{a list of options for the M-estimating equation, see
-\link{mest_options} for details.}
+\item{opts}{a list of options for the M-estimating algorithm, see
+\link{mscale_algorithm_options} for details.}
 }
 \value{
 the m-scale estimate.

man/mlocscale.Rd

@@ -6,7 +6,7 @@
 \usage{
 mlocscale(x, bdp = 0.25, location_rho = c("bisquare", "huber"),
   scale_cc = consistency_const(bdp, "bisquare"), location_cc,
-  opts = mest_options())
+  opts = mscale_algorithm_options())
 }
 \arguments{
 \item{x}{numeric values.}
@@ -14,7 +14,7 @@ mlocscale(x, bdp = 0.25, location_rho = c("bisquare", "huber"),
 \item{bdp}{desired breakdown point (between 0 and 0.5).}
 
 \item{opts}{a list of options for the M-estimating equation,
-see \link{mest_options} for details.}
+see \link{mscale_algorithm_options} for details.}
 
 \item{cc}{cutoff value for the bisquare rho function. By default, chosen
 for a consistent estimate under the Normal model.}

man/mscale.Rd

@@ -5,7 +5,7 @@
 \title{Compute the M-Scale of Centered Values}
 \usage{
 mscale(x, bdp = 0.25, cc = consistency_const(bdp, "bisquare"),
-  opts = mest_options())
+  opts = mscale_algorithm_options())
 }
 \arguments{
 \item{x}{numeric values.}
@@ -15,7 +15,7 @@ mscale(x, bdp = 0.25, cc = consistency_const(bdp, "bisquare"),
 \item{cc}{cutoff value for the bisquare rho function. By default, chosen
 for a consistent estimate under the Normal model.}
 
-\item{opts}{a list of options for the M-scale equation, see \link{mest_options}
+\item{opts}{a list of options for the M-scale estimation algorithm, see \link{mscale_algorithm_options}
 for details.}
 }
 \value{

man/mscale_algorithm_options.Rd

@@ -4,7 +4,7 @@
 \alias{mscale_algorithm_options}
 \title{Options for the M-scale Estimation Algorithm}
 \usage{
-mscale_algorithm_options(max_it = 200, eps = 1e-06)
+mscale_algorithm_options(max_it = 200, eps = 1e-08)
 }
 \arguments{
 \item{max_it}{maximum number of iterations.}

man/pense.Rd

@@ -8,7 +8,8 @@ pense(x, y, alpha, lambdas, cold_lambdas, penalty_loadings,
   additional_initial_estimates, include_intercept = TRUE, bdp = 0.25,
   cc, eps = 1e-06, algorithm_opts, nr_tracks = 10, explore_it = 10,
   sparse = TRUE, mscale_opts = mscale_algorithm_options(),
-  enpy_opts = enpy_options())
+  enpy_opts = enpy_options(), zero_based = TRUE,
+  share_initials = FALSE)
 }
 \arguments{
 \item{x}{\code{n} by \code{p} matrix of numeric predictors.}
@@ -55,6 +56,10 @@ solutions.}
 
 \item{enpy_opts}{options for the ENPY initial estimates, created with the
 \link{enpy_options} function. See \link{enpy_initial_estimates} for details.}
+
+\item{zero_based}{also compute the PENSE regularization path starting from the 0 vector.}
+
+\item{share_initials}{use all initial estimates for every lambda value.}
 }
 \description{
 Compute the PENSE Regularization Path

man/pense_admm_options.Rd

@@ -5,7 +5,7 @@
 \title{Options for the S-Estimate Algorithm}
 \usage{
 pense_admm_options(max_it = 5000, tau, sparse = TRUE,
-  prox_eps = 1e-09, prox_max_it = 200, prox_oscillate_window = 4,
+  prox_eps = 1e-08, prox_max_it = 200, prox_oscillate_window = 4,
   prox_minimum_step_size = 0.01, prox_wolfe_c1 = 1e-04,
   prox_wolfe_c2 = 0.9, prox_step_size_adj = 0.9,
   prox_max_small_steps = 10)

src/r_interface.cc

@@ -36,7 +36,6 @@ const R_CallMethodDef kExportedCallMethods[] = {
   {"C_lsen_regression", (DL_FUNC) &LsEnRegression, 5},
   {"C_pense_regression", (DL_FUNC) &PenseEnRegression, 8},
   {"C_pense_max_lambda", (DL_FUNC) &PenseMaxLambda, 4},
-  {"C_s_proximal_op", (DL_FUNC) &SProximalOp, 7},
   // {"C_pensem_en_regression_dal", (DL_FUNC) &PensemEnRegressionDAL, 7},
   // {"C_pensem_ridge_regression", (DL_FUNC) &PensemRidgeRegression, 6},
   {"C_penpy", (DL_FUNC) &PenPyInitialEstimator, 6},

src/r_pense_regression.cc

@@ -8,6 +8,8 @@
 
 #include "r_pense_regression.hpp"
 
+#include <unordered_set>
+
 #include "rcpp_integration.hpp"
 #include "r_interface_utils.hpp"
 #include "alias.hpp"
@@ -43,12 +45,56 @@ using DenseCoefs = nsoptim::RegressionCoefficients<arma::vec>;
 template<typename Optimizer>
 using StartCoefficientsList = FwdList<FwdList<typename Optimizer::Coefficients>>;
 
+template<typename Optimizer>
+using CoefficientsList = FwdList<typename Optimizer::Coefficients>;
+
 template<typename Optimizer>
 using PenaltyList = FwdList<typename Optimizer::PenaltyFunction>;
 
 namespace {
 constexpr int kDefaultExploreIt = 10;
 constexpr int kDefaultTracks = 10;
+constexpr bool kDefaultCompute0Based = true;
+constexpr bool kDefaultSameInitialForAll = false;
+
+//! A hash functor which always returns 0.
+template<typename T>
+class DummyHasher {
+ public:
+  constexpr std::size_t operator()(const typename T::Optimum&) const noexcept {
+    return 0;
+  }
+};
+
+//! A functor to compare two optima for equivalence.
+template<typename T>
+class OptimumComparator {
+ public:
+  explicit OptimumComparator(const double eps) noexcept : eps_(eps) {}
+  bool operator()(const typename T::Optimum& a, const typename T::Optimum& b) const noexcept {
+    // First check if the value of the objective function is similar.
+    if (std::abs(a.objf_value - b.objf_value) < eps_) {
+      // Value of the objective function is similar. Check if the intercept is similar.
+      const double int_diff = a.coefs.intercept - b.coefs.intercept;
+      if (int_diff * int_diff < eps_) {
+        // The intercept is similar. Check if the slope is also similar.
+        const double beta_diff = arma::norm(a.coefs.beta - b.coefs.beta, 2);
+        if (beta_diff * beta_diff < eps_) {
+          // The slope is also similar. Return true.
+          return true;
+        }
+      }
+    }
+    return false;
+  }
+
+ private:
+  const double eps_;
+};
+
+//! Alias for a set of optima.
+template<typename T>
+using OptimumSet = std::unordered_set<typename T::Optimum, DummyHasher<T>, OptimumComparator<T>>;
 
 //! Expand a list of PY Results to a list of start coefficients.
 //!
@@ -231,60 +277,92 @@ SEXP PenseRegressionImpl(SOptimizer optimizer, SEXP r_x, SEXP r_y, SEXP r_penalt
   optimizer.convergence_tolerance(eps);
 
   Metrics metrics("pense");
+  FwdList<OptimumSet<SOptimizer>> reg_path;
+
+  // Initialize the list of optima with empty sets.
+  for (auto pen_it = penalties.cbegin(), pen_end = penalties.cend(); pen_it != pen_end; ++pen_it) {
+    reg_path.emplace_front(1, DummyHasher<SOptimizer>(), OptimumComparator<SOptimizer>(eps));
+  }
 
   // Compute the initial estimators
   const auto cold_starts = EnpyInitialEstimates<SOptimizer>(loss, penalties, r_penalties, r_enpy_inds, r_enpy_opts,
                                                             optional_args, &metrics);
 
   // Compute the cold-based solutions
-  FwdList<Optima<SOptimizer>> reg_path_cold;
-
   if (!cold_starts.empty()) {
-    reg_path_cold = RegularizationPath(loss, penalties, optimizer, cold_starts,
-                                       GetFallback(pense_opts, "explore_it", kDefaultExploreIt),
-                                       GetFallback(pense_opts, "nr_tracks", kDefaultTracks), std::sqrt(eps));
+    auto reg_path_cold = RegularizationPath(loss, penalties, optimizer, cold_starts,
+                                            GetFallback(pense_opts, "explore_it", kDefaultExploreIt),
+                                            GetFallback(pense_opts, "nr_tracks", kDefaultTracks), std::sqrt(eps));
+
+    auto insert_it = reg_path.begin();
+    for (auto&& optima : reg_path_cold) {
+      for (auto&& optimum : optima) {
+        insert_it->emplace(std::move(optimum));
+      }
+      ++insert_it;
+    }
   }
 
   //! Check for a user interrupt...
   Rcpp::checkUserInterrupt();
 
   // Compute the 0-based solutions
-  auto reg_path_0 = RegularizationPath(loss, penalties, optimizer);
+  if (GetFallback(pense_opts, "compute_0_based", kDefaultCompute0Based)) {
+    auto reg_path_0 = RegularizationPath(loss, penalties, optimizer);
+    auto insert_it = reg_path.begin();
+    for (auto&& optimum : reg_path_0) {
+      insert_it->emplace(std::move(optimum));
+      ++insert_it;
+    }
 
-  //! Check for a user interrupt...
-  Rcpp::checkUserInterrupt();
+    //! Check for a user interrupt...
+    Rcpp::checkUserInterrupt();
+  }
 
   // Compute the "others"-based solutions
-  auto other_initial_ests = as<StartCoefficientsList<SOptimizer>>(r_initial_ests);
-  FwdList<Optima<SOptimizer>> reg_path_others;
-  if (!other_initial_ests.empty()) {
-    reg_path_others = RegularizationPath(loss, penalties, optimizer, other_initial_ests,
-                                         GetFallback(pense_opts, "explore_it", kDefaultExploreIt),
-                                         GetFallback(pense_opts, "nr_tracks", kDefaultTracks), std::sqrt(eps));
+  // FwdList<Optima<SOptimizer>> reg_path_others;
+  if (GetFallback(pense_opts, "same_initials_for_all", kDefaultSameInitialForAll)) {
+    // For all penalties use the same initial estimates.
+    auto initial_ests = as<CoefficientsList<SOptimizer>>(r_initial_ests);
+    for (auto&& start : initial_ests) {
+      auto reg_path_others = RegularizationPath(loss, penalties, optimizer, start);
+
+      // Add optimum to the other optima.
+      auto insert_it = reg_path.begin();
+      for (auto&& optimum : reg_path_others) {
+        insert_it->emplace(std::move(optimum));
+        insert_it++;
+      }
+      //! Check for a user interrupt...
+      Rcpp::checkUserInterrupt();
+    }
+  } else {
+    // For every penalty use separate initial estimates.
+    auto other_initial_ests = as<StartCoefficientsList<SOptimizer>>(r_initial_ests);
+    if (!other_initial_ests.empty()) {
+      auto reg_path_others = RegularizationPath(loss, penalties, optimizer, other_initial_ests,
+                                                GetFallback(pense_opts, "explore_it", kDefaultExploreIt),
+                                                GetFallback(pense_opts, "nr_tracks", kDefaultTracks), std::sqrt(eps));
+      auto insert_it = reg_path.begin();
+      for (auto&& optima : reg_path_others) {
+        for (auto&& optimum : optima) {
+          insert_it->emplace(std::move(optimum));
+        }
+        ++insert_it;
+      }
+    }
   }
 
-  // Combine all solutions
-  auto reg_path_0_it = reg_path_0.begin();
-  auto reg_path_cold_it = reg_path_cold.begin();
-  auto reg_path_others_it = reg_path_others.begin();
-  auto reg_path_0_end = reg_path_0.end();
-  auto reg_path_cold_end = reg_path_cold.end();
-  auto reg_path_others_end = reg_path_others.end();
-
-  auto&& reg_path_0_metrics = metrics.CreateSubMetrics("reg_path_0");
-  auto&& reg_path_cold_metrics = metrics.CreateSubMetrics("reg_path_cold");
-  auto&& reg_path_others_metrics = metrics.CreateSubMetrics("reg_path_others");
-
   Rcpp::List combined_reg_path;
-  while ((reg_path_0_it != reg_path_0_end) || (reg_path_cold_it != reg_path_cold_end) ||
-         (reg_path_others_it != reg_path_others_end)) {
+  for (auto&& optima : reg_path) {
     Rcpp::List solutions;
-
-    reg_path_0_it = AddEstimate(reg_path_0_it, reg_path_0_end, &reg_path_0_metrics, &solutions, "0");
-    reg_path_cold_it = AddEstimates(reg_path_cold_it, reg_path_cold_end, &reg_path_cold_metrics, &solutions, "cold");
-    reg_path_others_it = AddEstimates(reg_path_others_it, reg_path_others_end, &reg_path_others_metrics, &solutions,
-                                      "others");
-
+    Metrics& sub_metrics = metrics.CreateSubMetrics("lambda");
+    for (auto&& optimum : optima) {
+      if (optimum.metrics) {
+        sub_metrics.AddSubMetrics(*optimum.metrics);
+      }
+      solutions.push_back(WrapOptimum(optimum));
+    }
     combined_reg_path.push_back(solutions);
   }
 
@@ -500,23 +578,5 @@ SEXP PenseMaxLambda(SEXP r_x, SEXP r_y, SEXP r_pense_opts, SEXP r_optional_args)
   END_RCPP
 }
 
-SEXP SProximalOp(SEXP r_u, SEXP r_v_prev, SEXP r_x, SEXP r_y, SEXP r_lambda, SEXP r_prox_opts,
-                 SEXP r_mscale_opts) noexcept {
-  BEGIN_RCPP
-  ConstRegressionDataPtr data = MakePredictorResponseData(r_x, r_y);
-  Mscale<RhoBisquare> mscale(as<Rcpp::List>(r_mscale_opts));
-  PenseProximalOperator proxop(as<Rcpp::List>(r_prox_opts));
-  SLoss loss(data, mscale, true);
-  auto u = MakeVectorView(r_u);
-  auto v_prev = MakeVectorView(r_v_prev);
-
-  proxop.loss(&loss);
-  const arma::vec v = proxop(*u, *v_prev, 0, *REAL(r_lambda));
-
-  return Rcpp::wrap(v);
-
-  END_RCPP
-}
-
 }  // namespace r_interface
 }  // namespace pense

src/r_pense_regression.hpp

@@ -37,18 +37,6 @@ SEXP PenseEnRegression(SEXP x, SEXP y, SEXP penalties, SEXP initial_ests, SEXP e
 //!                       `pen_loadings` ... optional vector of length `p` with non-negative penalty loadings.
 SEXP PenseMaxLambda(SEXP x, SEXP y, SEXP pense_opts, SEXP optional_args) noexcept;
 
-
-//! Proximal operator for the S-loss.
-//!
-//! @param u numeric argument to the proximal operator with `n` elements.
-//! @param v_prev starting point for the proximal operator with `n` elements.
-//! @param x numeric predictor matrix with `n` rows and `p` columns.
-//! @param y numeric response vector with `n` elements.
-//! @param lambda proximal operator.
-//! @param prox_opts options for the proximal operator.
-//! @param mscale_opts options for the M-scale.
-SEXP SProximalOp(SEXP u, SEXP v_prev, SEXP x, SEXP y, SEXP lambda, SEXP prox_opts, SEXP mscale_opts) noexcept;
-
 }  // namespace r_interface
 }  // namespace pense