Safeguard in ConvexifiedHessian if regularization fails

uno/ingredients/hessian_models/ConvexifiedHessian.cpp

@@ -1,12 +1,15 @@
 // Copyright (c) 2024 Charlie Vanaret
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
+#include <stdexcept>
 #include "ConvexifiedHessian.hpp"
+#include "ingredients/hessian_models/UnstableRegularization.hpp"
 #include "reformulation/OptimizationProblem.hpp"
 #include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
 #include "solvers/SymmetricIndefiniteLinearSolverFactory.hpp"
 #include "tools/Logger.hpp"
 #include "options/Options.hpp"
+#include "tools/Infinity.hpp"
 #include "tools/Statistics.hpp"
 
 namespace uno {
@@ -15,7 +18,8 @@ namespace uno {
          // inertia-based convexification needs a linear solver
          linear_solver(SymmetricIndefiniteLinearSolverFactory::create(dimension, maximum_number_nonzeros, options)),
          regularization_initial_value(options.get_double("regularization_initial_value")),
-         regularization_increase_factor(options.get_double("regularization_increase_factor")) {
+         regularization_increase_factor(options.get_double("regularization_increase_factor")),
+         regularization_failure_threshold(options.get_double("regularization_failure_threshold")) {
    }
 
    void ConvexifiedHessian::evaluate(Statistics& statistics, const OptimizationProblem& problem, const Vector<double>& primal_variables,
@@ -57,7 +61,9 @@ namespace uno {
          else {
             DEBUG << "rank: " << this->linear_solver->rank() << ", negative eigenvalues: " << this->linear_solver->number_negative_eigenvalues() << '\n';
             regularization_factor = (regularization_factor == 0.) ? this->regularization_initial_value : this->regularization_increase_factor * regularization_factor;
-            assert(is_finite(regularization_factor) && "The regularization coefficient diverged");
+            if (regularization_factor > this->regularization_failure_threshold) {
+               throw UnstableRegularization();
+            }
          }
       }
       statistics.set("regularization", regularization_factor);

uno/ingredients/hessian_models/ConvexifiedHessian.hpp

@@ -22,6 +22,7 @@ namespace uno {
       std::unique_ptr<DirectSymmetricIndefiniteLinearSolver<size_t, double>> linear_solver; /*!< Solver that computes the inertia */
       const double regularization_initial_value{};
       const double regularization_increase_factor{};
+      const double regularization_failure_threshold{};
 
       void regularize(Statistics& statistics, SymmetricMatrix<size_t, double>& hessian, size_t number_original_variables);
    };

uno/ingredients/hessian_models/UnstableRegularization.hpp

@@ -0,0 +1,18 @@
+// Copyright (c) 2024 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#ifndef UNO_UNSTABLEREGULARIZATION_H
+#define UNO_UNSTABLEREGULARIZATION_H
+
+#include <exception>
+
+namespace uno {
+   struct UnstableRegularization : public std::exception {
+
+      [[nodiscard]] const char* what() const noexcept override {
+         return "The inertia correction got unstable (delta_w > threshold)";
+      }
+   };
+} // namespace
+
+#endif // UNO_UNSTABLEREGULARIZATION_H

uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp

@@ -8,19 +8,13 @@
 #include "SymmetricMatrix.hpp"
 #include "SparseStorageFactory.hpp"
 #include "RectangularMatrix.hpp"
+#include "ingredients/hessian_models/UnstableRegularization.hpp"
 #include "model/Model.hpp"
 #include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
 #include "options/Options.hpp"
 #include "tools/Statistics.hpp"
 
 namespace uno {
-   struct UnstableRegularization : public std::exception {
-
-      [[nodiscard]] const char* what() const noexcept override {
-         return "The inertia correction got unstable (delta_w > threshold)";
-      }
-   };
-
    template <typename ElementType>
    class SymmetricIndefiniteLinearSystem {
    public:
@@ -158,7 +152,7 @@ namespace uno {
             this->previous_primal_regularization = this->primal_regularization;
          }
          else {
-            auto [number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues] = linear_solver.get_inertia();
+            std::tie(number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues) = linear_solver.get_inertia();
             DEBUG << "Expected inertia (" << size_primal_block << ", " << size_dual_block << ", 0), ";
             DEBUG << "got (" << number_pos_eigenvalues << ", " << number_neg_eigenvalues << ", " << number_zero_eigenvalues << ")\n";
             DEBUG << "Number of attempts: " << number_attempts << "\n";