From 607b75411912b198ba11b8b0a7b6b4f92fef4424 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 4 Nov 2024 12:34:30 +0100 Subject: [PATCH] Throw exception UnstableRegularization in ConvexifiedHessian if regularization factor becomes too large + moved UnstableRegularization to its own file --- .../hessian_models/ConvexifiedHessian.cpp | 10 ++++++++-- .../hessian_models/ConvexifiedHessian.hpp | 1 + .../hessian_models/UnstableRegularization.hpp | 18 ++++++++++++++++++ .../SymmetricIndefiniteLinearSystem.hpp | 10 ++-------- 4 files changed, 29 insertions(+), 10 deletions(-) create mode 100644 uno/ingredients/hessian_models/UnstableRegularization.hpp diff --git a/uno/ingredients/hessian_models/ConvexifiedHessian.cpp b/uno/ingredients/hessian_models/ConvexifiedHessian.cpp index e0f16d4f..5e76b1fb 100644 --- a/uno/ingredients/hessian_models/ConvexifiedHessian.cpp +++ b/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 #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& 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); diff --git a/uno/ingredients/hessian_models/ConvexifiedHessian.hpp b/uno/ingredients/hessian_models/ConvexifiedHessian.hpp index 1a090eff..795446b3 100644 --- a/uno/ingredients/hessian_models/ConvexifiedHessian.hpp +++ b/uno/ingredients/hessian_models/ConvexifiedHessian.hpp @@ -22,6 +22,7 @@ namespace uno { std::unique_ptr> 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& hessian, size_t number_original_variables); }; diff --git a/uno/ingredients/hessian_models/UnstableRegularization.hpp b/uno/ingredients/hessian_models/UnstableRegularization.hpp new file mode 100644 index 00000000..3439437c --- /dev/null +++ b/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 + +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 \ No newline at end of file diff --git a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp index e920e919..75091517 100644 --- a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp +++ b/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 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";