feat!: support generation of BH-type lambda sequence

Also add tests for the new sequence and exceptions to check for reasonable input.

CMakeLists.txt

@@ -34,7 +34,8 @@ if(BUILD_TESTING)
     FetchContent_MakeAvailable(Catch2)
   endif()
 
-  add_executable(tests tests/gaussian.cpp tests/prox.cpp tests/qnorm.cpp)
+  add_executable(tests tests/gaussian.cpp tests/prox.cpp tests/qnorm.cpp
+                       tests/lambda_sequence.cpp)
   target_link_libraries(tests PRIVATE libslope Catch2::Catch2WithMain)
   target_include_directories(tests PUBLIC tests/ "${PROJECT_BINARY_DIR}")
 

src/slope/regularization_sequence.cpp

@@ -1,19 +1,27 @@
 #include "regularization_sequence.h"
 #include "qnorm.h"
-#include <Eigen/Core>
+#include <cassert>
 
 namespace slope {
 
 Eigen::ArrayXd
-lambdaSequence(const int p, const double q)
+lambdaSequence(const int p, const double q, const std::string& type)
 {
-  Eigen::ArrayXd lambda_sequence(p);
+  Eigen::ArrayXd lambda(p);
 
-  for (int j = 0; j < p; ++j) {
-    lambda_sequence[j] = normalQuantile(1.0 - (j + 1.0) * q / (2.0 * p));
+  if (type == "bh") {
+    if (q < 0 || q > 1) {
+      throw std::invalid_argument("q must be between 0 and 1");
+    }
+
+    for (int j = 0; j < p; ++j) {
+      lambda[j] = normalQuantile(1.0 - (j + 1.0) * q / (2.0 * p));
+    }
   }
 
-  return lambda_sequence;
+  assert(lambda.minCoeff() > 0);
+
+  return lambda;
 }
 
 } // namespace slope

src/slope/regularization_sequence.h

@@ -5,6 +5,7 @@
 #include <Eigen/Core>
 #include <Eigen/SparseCore>
 #include <memory>
+#include <string>
 
 namespace slope {
 
@@ -19,7 +20,7 @@ namespace slope {
  * @return An Eigen::ArrayXd containing the generated lambda values.
  */
 Eigen::ArrayXd
-lambdaSequence(const int p, const double q);
+lambdaSequence(const int p, const double q, const std::string& type);
 
 template<typename T>
 Eigen::ArrayXd

src/slope/slope.h

@@ -11,6 +11,7 @@
 #include "sorted_l1_norm.h"
 #include "standardize.h"
 #include <Eigen/Sparse>
+#include <cassert>
 #include <iostream>
 #include <memory>
 #include <vector>
@@ -79,6 +80,9 @@ computeMaxDelta(const T& x,
  *   "gaussian".
  * @param intercept Whether to fit an intercept term. Default is true.
  * @param standardize Whether to standardize the predictors. Default is true.
+ * @param lambdaSequence The type of regularization parameter sequence to use.
+ * @param q The q parameter for computation of the lambda sequence. Default is
+ *   0.1.
  * @param path_length The number of steps in the regularization path. Default is
  *   100.
  * @param alpha_min_ratio The minimum ratio of alpha to the maximum alpha value.
@@ -106,6 +110,8 @@ slope(const T& x,
       const std::string objective_choice = "gaussian",
       bool intercept = true,
       bool standardize = true,
+      std::string lambda_type = "bh",
+      double q = 0.1,
       int path_length = 100,
       double alpha_min_ratio = 1e-4,
       int pgd_freq = 10,
@@ -120,8 +126,6 @@ slope(const T& x,
   const int n = x.rows();
   const int p = x.cols();
 
-  SortedL1Norm sl1_norm{ lambda };
-
   // standardize
   VectorXd x_centers(p);
   VectorXd x_scales(p);
@@ -146,8 +150,22 @@ slope(const T& x,
 
   VectorXd residual = z - eta;
 
+  // Setup the regularization sequence and path
+  SortedL1Norm sl1_norm{ lambda };
+
+  if (lambda.size() == 0) {
+    lambda = lambdaSequence(p, q, lambda_type);
+  } else {
+    if (lambda.size() != p) {
+      throw std::invalid_argument(
+        "lambda must be the same length as the number of predictors");
+    }
+    if (lambda.minCoeff() < 0) {
+      throw std::invalid_argument("lambda must be non-negative");
+    }
+  }
+
   if (alpha.size() == 0) {
-    // No user-supplied alpha sequence, so generate one.
     alpha = regularizationPath(x,
                                w,
                                z,

tests/lambda_sequence.cpp

@@ -0,0 +1,29 @@
+#include "test_helpers.hpp"
+#include <Eigen/Core>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers.hpp>
+#include <slope/regularization_sequence.h>
+
+TEST_CASE("Test that regularization sequence generation works",
+          "[regularization sequence]")
+{
+  SECTION("Test that BH sequence works")
+  {
+    Eigen::ArrayXd lambda1 = slope::lambdaSequence(4, 0.1, "bh");
+
+    std::vector<double> lambda1_expected = {
+      2.24140272760495, 1.95996398454005, 1.78046434169203, 1.64485362695147
+    };
+
+    REQUIRE_THAT(lambda1, VectorApproxEqual(lambda1_expected, 1e-6));
+
+    Eigen::ArrayXd lambda2 = slope::lambdaSequence(6, 0.5, "bh");
+
+    std::vector<double> lambda2_expected = {
+      1.73166439612225,  1.38299412710064,  1.15034938037601,
+      0.967421566101701, 0.812217801499913, 0.674489750196082
+    };
+
+    REQUIRE_THAT(lambda2, VectorApproxEqual(lambda2_expected, 1e-6));
+  }
+}