From 209414c388423bf71d9680ec4f2ce00751c3b9c2 Mon Sep 17 00:00:00 2001
From: Johan Larsson <johanlarsson@outlook.com>
Date: Tue, 6 Feb 2024 13:13:20 +0100
Subject: [PATCH] fix: avoid infinite loop by bounding stopping criterion

Introduce an EPSILON constant at 1e-10 to the relative stopping criterion,
to avoid machine precision issues when the primal is very small. Also
add a new test to check sparse vs dense comparisons
---
 CMakeLists.txt        | 10 +++++++--
 src/slope/constants.h |  8 +++++++
 src/slope/slope.h     |  5 +++--
 tests/gaussian.cpp    |  2 +-
 tests/sparse.cpp      | 49 +++++++++++++++++++++++++++++++++++++++++++
 5 files changed, 69 insertions(+), 5 deletions(-)
 create mode 100644 src/slope/constants.h
 create mode 100644 tests/sparse.cpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 531d2e9..94e38ff 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -36,8 +36,14 @@ if(BUILD_TESTING)
   endif()
 
   add_executable(
-    tests tests/gaussian.cpp tests/prox.cpp tests/qnorm.cpp
-          tests/lambda_sequence.cpp tests/path.cpp tests/standardize.cpp)
+    tests
+    tests/gaussian.cpp
+    tests/prox.cpp
+    tests/qnorm.cpp
+    tests/lambda_sequence.cpp
+    tests/path.cpp
+    tests/standardize.cpp
+    tests/sparse.cpp)
   target_link_libraries(tests PRIVATE slope Catch2::Catch2WithMain
                                       Eigen3::Eigen)
   target_include_directories(tests PUBLIC tests/ "${PROJECT_BINARY_DIR}")
diff --git a/src/slope/constants.h b/src/slope/constants.h
new file mode 100644
index 0000000..23da389
--- /dev/null
+++ b/src/slope/constants.h
@@ -0,0 +1,8 @@
+/**
+ * @internal
+ * @file
+ * @brief Definitions of constants used in libslope
+ */
+#pragma once
+
+const double EPSILON = 1e-10;
diff --git a/src/slope/slope.h b/src/slope/slope.h
index 5a70d6d..29d4ff1 100644
--- a/src/slope/slope.h
+++ b/src/slope/slope.h
@@ -7,6 +7,7 @@
 
 #include "cd.h"
 #include "clusters.h"
+#include "constants.h"
 #include "helpers.h"
 #include "math.h"
 #include "objectives.h"
@@ -390,7 +391,7 @@ class Slope
 
         dual_gaps.emplace_back(dual_gap);
 
-        double tol_scaled = std::abs(primal) * this->tol;
+        double tol_scaled = (std::abs(primal) + EPSILON) * this->tol;
 
         if (this->print_level > 1) {
           std::cout << indent(1) << "IRLS iteration: " << it_outer << std::endl
@@ -431,7 +432,7 @@ class Slope
 
             double dual_gap_inner = primal_inner - dual_inner;
 
-            double tol_inner = std::abs(primal_inner) * this->tol;
+            double tol_inner = (std::abs(primal_inner) + EPSILON) * this->tol;
 
             if (this->print_level > 2) {
               std::cout << indent(2) << "iteration: " << it << std::endl
diff --git a/tests/gaussian.cpp b/tests/gaussian.cpp
index f0654bc..f0cf56a 100644
--- a/tests/gaussian.cpp
+++ b/tests/gaussian.cpp
@@ -42,7 +42,7 @@ TEST_CASE("Simple low-dimensional design", "[gaussian][basic]")
   double gap = dual_gaps.front().back();
   double primal = primals.front().back();
 
-  REQUIRE(gap <= primal * 1e-4);
+  REQUIRE(gap <= (primal + 1e-10) * 1e-4);
 }
 
 TEST_CASE("X is identity", "[gaussian][identity]")
diff --git a/tests/sparse.cpp b/tests/sparse.cpp
new file mode 100644
index 0000000..36a4799
--- /dev/null
+++ b/tests/sparse.cpp
@@ -0,0 +1,49 @@
+#include "test_helpers.hpp"
+#include <Eigen/Core>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers.hpp>
+#include <slope/slope.h>
+
+TEST_CASE("Sparse and dense methods agree", "[gaussian][sparse]")
+{
+  using namespace Catch::Matchers;
+
+  Eigen::Vector<double, 3> beta;
+  Eigen::Vector<double, 10> y;
+  Eigen::Matrix<double, 10, 3> x_dense;
+
+  x_dense << 0.0, 0.13339576, 0.49361983, 0.17769259, 0.66565742, 0.36972579,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.94280368, 0.0, 0.0, 0.3499374, 0.0, 0.22377115,
+    0.0, 0.0, 0.96893287, 0.95858229, 0.70486475, 0.60885162, 0.0, 0.0,
+    0.92902639, 0.0, 0.4978676, 0.0, 0.50022619;
+  Eigen::SparseMatrix<double> x_sparse = x_dense.sparseView();
+
+  beta << 1, 2, -0.9;
+
+  y = x_dense * beta;
+
+  Eigen::ArrayXd alpha = Eigen::ArrayXd::Zero(1);
+  Eigen::Array<double, 3, 1> lambda;
+  lambda << 0.5, 0.5, 0.2;
+
+  slope::Slope model;
+
+  model.setIntercept(false);
+  model.setStandardize(false);
+  model.setPrintLevel(3);
+  model.setPgdFreq(1);
+
+  model.fit(x_sparse, y, alpha, lambda);
+  auto coefs_sparse = model.getCoefs();
+
+  model.fit(x_dense, y, alpha, lambda);
+  auto coefs_dense = model.getCoefs();
+
+  auto dual_gaps = model.getDualGaps();
+  auto primals = model.getPrimals();
+
+  Eigen::VectorXd coef_sparse = coefs_sparse.col(0);
+  Eigen::VectorXd coef_dense = coefs_dense.col(0);
+
+  REQUIRE_THAT(coef_sparse, VectorApproxEqual(coef_dense, 1e-4));
+}