diff --git a/CMakeLists.txt b/CMakeLists.txt index e05f4cab6..9d1cc0fe2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -68,6 +68,7 @@ file(GLOB UNO_SOURCE_FILES uno/ingredients/inertia_correction_strategies/*.cpp uno/ingredients/subproblem/*.cpp uno/ingredients/subproblem_solvers/*.cpp + uno/ingredients/subproblem_solvers/dogleg/*.cpp uno/model/*.cpp uno/optimization/*.cpp uno/options/*.cpp diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp index 09ef10000..5b45c08e9 100644 --- a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp @@ -179,7 +179,6 @@ namespace uno { subproblem_solver.solve(statistics, subproblem, trust_region_radius, this->initial_point, direction, current_evaluations, warmstart_information); ++this->number_subproblems_solved; - direction.norm = norm_inf(view(direction.primals, 0, subproblem.problem.get_number_original_variables())); this->initial_point.fill(0.); DEBUG3 << direction << '\n'; } diff --git a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp index bb497c304..6b8a12738 100644 --- a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp @@ -76,7 +76,6 @@ namespace uno { this->initial_point.fill(0.); this->subproblem_solver->solve(statistics, subproblem, trust_region_radius, this->initial_point, direction, current_evaluations, warmstart_information); - direction.norm = norm_inf(view(direction.primals, 0, this->original_problem.get_number_original_variables())); DEBUG3 << direction << '\n'; warmstart_information.no_changes(); } diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp index ebe8a4443..ee2efe760 100644 --- a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp +++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp @@ -226,6 +226,7 @@ namespace uno { direction.primals[variable_index] = std::min(std::max(direction.primals[variable_index], this->lower_bounds[variable_index]), this->upper_bounds[variable_index]); } + direction.norm = norm_inf(view(direction.primals, 0, subproblem.problem.get_number_original_variables())); // gather the multipliers this->set_multipliers(subproblem.number_variables, direction.multipliers); LPSolver::compute_dual_displacements(subproblem, direction.multipliers); diff --git a/uno/ingredients/subproblem_solvers/COOLinearSystem.cpp b/uno/ingredients/subproblem_solvers/COOLinearSystem.cpp index 3a5398b6d..327c55a8e 100644 --- a/uno/ingredients/subproblem_solvers/COOLinearSystem.cpp +++ b/uno/ingredients/subproblem_solvers/COOLinearSystem.cpp @@ -3,7 +3,6 @@ #include "COOLinearSystem.hpp" #include "ingredients/subproblem/Subproblem.hpp" -#include "linear_algebra/Indexing.hpp" namespace uno { COOLinearSystem::COOLinearSystem(int solver_indexing): solver_indexing(solver_indexing) { diff --git a/uno/ingredients/subproblem_solvers/EQPSolver.cpp b/uno/ingredients/subproblem_solvers/EQPSolver.cpp index 7cdeb98bf..a510aafce 100644 --- a/uno/ingredients/subproblem_solvers/EQPSolver.cpp +++ b/uno/ingredients/subproblem_solvers/EQPSolver.cpp @@ -17,6 +17,8 @@ namespace uno { linear_solver(SymmetricIndefiniteLinearSolverFactory::create(options.get_string("linear_solver"))) { } + EQPSolver::~EQPSolver() = default; + void EQPSolver::initialize_memory(const Subproblem& subproblem) { if (!subproblem.has_hessian_matrix()) { throw std::runtime_error("The subproblem does not have an explicit Hessian matrix and cannot be solved with a direct linear solver"); @@ -67,6 +69,7 @@ namespace uno { } // assemble the full primal-dual direction subproblem.assemble_primal_dual_direction(linear_system.solution, direction); + direction.norm = norm_inf(view(direction.primals, 0, subproblem.problem.get_number_original_variables())); } SolverWorkspace& EQPSolver::get_workspace() { diff --git a/uno/ingredients/subproblem_solvers/EQPSolver.hpp b/uno/ingredients/subproblem_solvers/EQPSolver.hpp index bf621c47e..1f2013445 100644 --- a/uno/ingredients/subproblem_solvers/EQPSolver.hpp +++ b/uno/ingredients/subproblem_solvers/EQPSolver.hpp @@ -16,7 +16,7 @@ namespace uno { class EQPSolver: public SubproblemSolver { public: explicit EQPSolver(const Options& options); - ~EQPSolver() override = default; + ~EQPSolver() override; void initialize_memory(const Subproblem& subproblem) override; diff --git a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp index 34766ed31..dbfc5be4d 100644 --- a/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp +++ b/uno/ingredients/subproblem_solvers/HiGHS/HiGHSSolver.cpp @@ -107,6 +107,7 @@ namespace uno { direction.multipliers.upper_bounds[variable_index] = bound_multiplier; } } + direction.norm = norm_inf(view(direction.primals, 0, subproblem.problem.get_number_original_variables())); // gather the multipliers for (size_t constraint_index = 0; constraint_index < subproblem.number_constraints; constraint_index++) { direction.multipliers.constraints[constraint_index] = solution.row_dual[constraint_index]; diff --git a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp b/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp index 8cb9c6611..c17423925 100644 --- a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp +++ b/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp @@ -12,6 +12,7 @@ #include "QPSolver.hpp" #include "QPSolverFactory.hpp" #include "WoodburyEQPSolver.hpp" +#include "dogleg/DoglegMethod.hpp" #include "ingredients/subproblem/Subproblem.hpp" #include "options/Options.hpp" #include "tools/Logger.hpp" @@ -38,7 +39,16 @@ namespace uno { const Subproblem subproblem(problem, current_iterate, hessian_model, inertia_correction_strategy); // if no inequality constraint and no trust region, allocate EQP solver // temporary fix: this is set only in interior-point methods + /*if (!subproblem.has_inequality_constraints() && uses_trust_region && subproblem.is_hessian_positive_definite()) { + // use the dogleg method + DEBUG << "Trust-region and no inequality constraints in the subproblem, allocating a dogleg solver\n"; + auto subproblem_solver = std::make_unique(options); + subproblem_solver->initialize_memory(subproblem); + return subproblem_solver; + } + else*/ if (!subproblem.has_inequality_constraints() && !uses_trust_region && options.get_string("inequality_handling_method") == "interior_point") { + // no trust region if constexpr (std::is_same_v) { DEBUG << "No inequality constraints in the subproblem, allocating an EQP solver with L-BFGS Hessian\n"; // the hessian_model we pass has type LBFGSHessian @@ -54,7 +64,7 @@ namespace uno { } } // otherwise, allocate LP/QP solver, depending on the presence of curvature in the subproblem - if (!subproblem.has_curvature()) { + else if (!subproblem.has_curvature()) { if (subproblem.number_constraints == 0) { DEBUG << "No curvature and only bound constraints in the subproblem, allocating a box LP solver\n"; auto subproblem_solver = std::make_unique(); diff --git a/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp b/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp index 639f2d638..7c5f7d118 100644 --- a/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp +++ b/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp @@ -74,6 +74,7 @@ namespace uno { // assemble the full primal-dual direction subproblem.assemble_primal_dual_direction(solution_diagonal_part, direction); + direction.norm = norm_inf(view(direction.primals, 0, subproblem.problem.get_number_original_variables())); } SolverWorkspace& WoodburyEQPSolver::get_workspace() { diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp new file mode 100644 index 000000000..c8325cead --- /dev/null +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.cpp @@ -0,0 +1,36 @@ +// Copyright (c) 2025 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "DoglegMethod.hpp" +#include "ingredients/subproblem/Subproblem.hpp" +#include "optimization/Direction.hpp" +#include "tools/Logger.hpp" + +namespace uno { + DoglegMethod::DoglegMethod(const Options& options): + workspace(options) { + } + + void DoglegMethod::initialize_memory(const Subproblem& subproblem) { + this->workspace.initialize_memory(subproblem); + } + + void DoglegMethod::solve(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius, + const Vector& /*initial_point*/, Direction& direction, Evaluations& current_evaluations, + const WarmstartInformation& warmstart_information) { + // first try the Newton step. This is the solution if within the trust region + this->workspace.compute_newton_step(statistics, subproblem, direction, current_evaluations, warmstart_information); + if (this->workspace.newton_step_squared_norm <= std::pow(trust_region_radius, 2.)) { + DEBUG << "The Newton step is within the trust region\n"; + direction.norm = norm_2(view(direction.primals, 0, subproblem.problem.get_number_original_variables())); + return; + } + // if the trust region constraint is violated, compute the dogleg path: the broken path between the Cauchy step + // and the Newton step + this->workspace.compute_dogleg(subproblem, trust_region_radius, direction, current_evaluations, warmstart_information); + } + + [[nodiscard]] SolverWorkspace& DoglegMethod::get_workspace() { + return this->workspace; + } +} // namespace \ No newline at end of file diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp new file mode 100644 index 000000000..29e812120 --- /dev/null +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegMethod.hpp @@ -0,0 +1,32 @@ +// Copyright (c) 2025 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNO_DOGLEGMETHOD_H +#define UNO_DOGLEGMETHOD_H + +#include "../SubproblemSolver.hpp" +#include "DoglegWorkspace.hpp" + +namespace uno { + // forward declaration + class Options; + + class DoglegMethod: public SubproblemSolver { + public: + explicit DoglegMethod(const Options& options); + ~DoglegMethod() override = default; + + void initialize_memory(const Subproblem& subproblem) override; + + void solve(Statistics& statistics, const Subproblem& subproblem, double trust_region_radius, + const Vector& initial_point, Direction& direction, Evaluations& current_evaluations, + const WarmstartInformation& warmstart_information) override; + + [[nodiscard]] SolverWorkspace& get_workspace() override; + + protected: + DoglegWorkspace workspace; + }; +} // namespace + +#endif // UNO_DOGLEGMETHOD_H \ No newline at end of file diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp new file mode 100644 index 000000000..4d649128b --- /dev/null +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.cpp @@ -0,0 +1,116 @@ +// Copyright (c) 2025 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "DoglegWorkspace.hpp" +#include "ingredients/subproblem/Subproblem.hpp" +#include "optimization/Direction.hpp" +#include "optimization/Iterate.hpp" +#include "optimization/OptimizationProblem.hpp" +#include "optimization/WarmstartInformation.hpp" +#include "symbolic/Sum.hpp" +#include "tools/Logger.hpp" + +namespace uno { + DoglegWorkspace::DoglegWorkspace(const Options& options): + newton_solver(options) { + } + + void DoglegWorkspace::initialize_memory(const Subproblem& subproblem) { + // Newton solver + this->newton_solver.initialize_memory(subproblem); + this->initial_point.resize(subproblem.number_variables); + // Newton step + this->g.resize(subproblem.number_variables); + this->newton_step.resize(subproblem.number_variables); + // Cauchy step + this->Hg.resize(subproblem.number_variables); + this->cauchy_step.resize(subproblem.number_variables); + } + + double DoglegWorkspace::compute_hessian_quadratic_form(const Subproblem& subproblem, const Vector& vector) const { + return this->newton_solver.get_workspace().compute_hessian_quadratic_form(subproblem, vector); + } + + void DoglegWorkspace::compute_newton_step(Statistics& statistics, const Subproblem& subproblem, Direction& direction, + Evaluations& current_evaluations, const WarmstartInformation& warmstart_information) { + if (warmstart_information.new_iterate) { + this->newton_solver.solve(statistics, subproblem, INF, this->initial_point, direction, current_evaluations, + warmstart_information); + this->newton_step = direction.primals; + this->newton_step_squared_norm = dot(this->newton_step, this->newton_step); + } + DEBUG << "Newton step: " << this->newton_step << '\n'; + } + + void DoglegWorkspace::compute_dogleg(const Subproblem& subproblem, double trust_region_radius, Direction& direction, + Evaluations& current_evaluations, const WarmstartInformation& warmstart_information) { + this->compute_cauchy_step(subproblem, current_evaluations, warmstart_information); + DEBUG << "Cauchy step: " << this->cauchy_step << '\n'; + double squared_norm_cauchy = dot(this->cauchy_step, this->cauchy_step); + if (trust_region_radius*trust_region_radius <= squared_norm_cauchy) { + DEBUG << "The Cauchy step is outside the trust region. Returning a clipped direction\n"; + direction.primals = this->cauchy_step; + direction.primals.scale(trust_region_radius/std::sqrt(squared_norm_cauchy)); + direction.norm = trust_region_radius; + DEBUG << "Scaled Cauchy direction: " << direction.primals << '\n'; + return; + } + + DEBUG << "Computing the dogleg step\n"; + // define temporary vectors + Vector u(subproblem.number_variables); + u = this->newton_step - this->cauchy_step; + Vector v(subproblem.number_variables); + v = 2.*this->cauchy_step - this->newton_step; + const double a = dot(u, u); + const double b = 2.*dot(u, v); + const double c = dot(v, v) - trust_region_radius*trust_region_radius; + DEBUG << "(a, b, c) = " << a << ", " << b << ", " << c << '\n'; + const double tau = DoglegWorkspace::compute_positive_root_quadratic_equation(a, b, c); + DEBUG << "tau = " << tau << '\n'; + if (1. <= tau && tau <= 2.) { + direction.primals = this->cauchy_step + (1. - tau)*(this->newton_step - this->cauchy_step); + direction.norm = norm_2(view(direction.primals, 0, subproblem.problem.get_number_original_variables())); + } + else { + throw std::runtime_error("The dogleg step could not be computed"); + } + } + + // private member functions + + void DoglegWorkspace::compute_cauchy_step(const Subproblem& subproblem, Evaluations& current_evaluations, + const WarmstartInformation& warmstart_information) { + if (warmstart_information.new_iterate) { + subproblem.problem.evaluate_objective_gradient(subproblem.current_iterate, this->g.data(), current_evaluations); + DEBUG << "g = " << this->g << '\n'; + // gᵀ g + this->g_squared_norm = dot(this->g, this->g); + // H g + subproblem.compute_hessian_vector_product(subproblem.current_iterate.primals.data(), this->g.data(), + this->Hg.data()); + DEBUG << "H g = " << this->Hg << '\n'; + // gᵀ H g + this->gHg = dot(this->g, this->Hg); + DEBUG << "g^T H g = " << this->gHg << '\n'; + if (this->gHg <= 0.) { + throw std::runtime_error("The Hessian is not positive definite"); + } + // Cauchy step: d_C = - (gᵀ g)/(gᵀ B g) g + this->cauchy_step = this->g; + const double scaling_factor = -this->g_squared_norm / this->gHg; + this->cauchy_step.scale(scaling_factor); + } + } + + // find the positive real root to ax^2 + bx + c = 0 + double DoglegWorkspace::compute_positive_root_quadratic_equation(double a, double b, double c) { + // avoid catastrophic cancellation + const double delta = b*b - 4.*a*c; + if (delta < 0.) { + throw std::runtime_error("No real root"); + } + // TODO check denominator + return (2.*c)/(-b - std::sqrt(delta)); + } +} // namespace \ No newline at end of file diff --git a/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp new file mode 100644 index 000000000..e2f74f03d --- /dev/null +++ b/uno/ingredients/subproblem_solvers/dogleg/DoglegWorkspace.hpp @@ -0,0 +1,54 @@ +// Copyright (c) 2025 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNO_DOGLEGWORKSPACE_H +#define UNO_DOGLEGWORKSPACE_H + +#include "../SolverWorkspace.hpp" +#include "ingredients/subproblem_solvers/EQPSolver.hpp" +#include "linear_algebra/Vector.hpp" +#include "tools/Infinity.hpp" + +namespace uno { + // forward declaration + class Direction; + class Evaluations; + class Options; + class Statistics; + class WarmstartInformation; + + class DoglegWorkspace: public SolverWorkspace { + public: + explicit DoglegWorkspace(const Options& options); + ~DoglegWorkspace() override = default; + + // Newton step + Vector g{}; + Vector newton_step{}; + double newton_step_squared_norm{INF}; + // Cauchy step + double g_squared_norm{INF}; + Vector Hg{}; + double gHg{INF}; + Vector cauchy_step{}; + + void initialize_memory(const Subproblem& subproblem); + + [[nodiscard]] double compute_hessian_quadratic_form(const Subproblem& subproblem, const Vector& vector) const override; + + void compute_newton_step(Statistics& statistics, const Subproblem& subproblem, Direction& direction, + Evaluations& current_evaluations, const WarmstartInformation& warmstart_information); + void compute_dogleg(const Subproblem& subproblem, double trust_region_radius, Direction& direction, + Evaluations& current_evaluations, const WarmstartInformation& warmstart_information); + + private: + mutable EQPSolver newton_solver; + Vector initial_point; + + void compute_cauchy_step(const Subproblem& subproblem, Evaluations& current_evaluations, + const WarmstartInformation& warmstart_information); + [[nodiscard]] static double compute_positive_root_quadratic_equation(double a, double b, double c); + }; +} // namespace + +#endif // UNO_DOGLEGWORKSPACE_H \ No newline at end of file