From bdd8be444e3ce742d3ab7f2efb152ba9ccfdcc5d Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Fri, 17 Apr 2026 16:37:02 +0200 Subject: [PATCH 1/8] First draft with code converted from JSOSolvers.jl by Claude --- CMakeLists.txt | 1 + .../InequalityHandlingMethodFactory.cpp | 5 + .../SubproblemSolverFactory.hpp | 7 + .../subproblem_solvers/TRON/TRONSolver.cpp | 486 ++++++++++++++++++ .../subproblem_solvers/TRON/TRONSolver.hpp | 114 ++++ 5 files changed, 613 insertions(+) create mode 100644 uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp create mode 100644 uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index e3e10d42e4..02e4f3b62f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -90,6 +90,7 @@ file(GLOB UNO_SOURCE_FILES CONFIGURE_DEPENDS uno/ingredients/inertia_correction_strategies/*.cpp uno/ingredients/subproblem/*.cpp uno/ingredients/subproblem_solvers/*.cpp + uno/ingredients/subproblem_solvers/TRON/*.cpp uno/model/*.cpp uno/optimization/*.cpp uno/options/*.cpp diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp index c4532d5833..9d9fb6dc64 100644 --- a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp +++ b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp @@ -28,6 +28,11 @@ namespace uno { } } // from now on, the problem has inequalities + // bound-constrained problem: do not reformulate + if (!problem.has_inequality_constraints() && problem.has_bound_constraints()) { + return std::make_unique("bound-constrained method"); + } + const std::string inequality_handling_method = options.get_string("inequality_handling_method"); // inequality-constrained methods if (inequality_handling_method == "inequality_constrained") { diff --git a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp b/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp index 03aeed312b..332554be46 100644 --- a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp +++ b/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp @@ -11,6 +11,7 @@ #include "InverseNewtonSolver.hpp" #include "QPSolverFactory.hpp" #include "WoodburyEQPSolver.hpp" +#include "TRON/TRONSolver.hpp" #include "ingredients/subproblem/Subproblem.hpp" #include "options/Options.hpp" #include "tools/Logger.hpp" @@ -52,6 +53,12 @@ namespace uno { return subproblem_solver; } } + // if only bound constraints, allocate bound-constrained solver + else if (!subproblem.has_inequality_constraints() && subproblem.has_bound_constraints()) { + auto subproblem_solver = std::make_unique(); + subproblem_solver->initialize_memory(subproblem); + return subproblem_solver; + } // if no inequality constraint and no trust region, allocate EQP solver else if (!subproblem.has_inequality_constraints() && !subproblem.has_bound_constraints() && !uses_trust_region) { if constexpr (std::is_same_v) { // unconstrained diff --git a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp new file mode 100644 index 0000000000..79c77e0391 --- /dev/null +++ b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp @@ -0,0 +1,486 @@ +// Copyright (c) 2026 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include +#include +#include +#include "TRONSolver.hpp" +#include "ingredients/subproblem/Subproblem.hpp" +#include "linear_algebra/BLAS.hpp" +#include "linear_algebra/Vector.hpp" + +namespace uno { + void TRONSolver::initialize_memory(const Subproblem& subproblem) { + this->lower_bounds.resize(subproblem.number_variables); + this->upper_bounds.resize(subproblem.number_variables); + this->workspace.objective_gradient.resize(subproblem.number_variables); + } + + void TRONSolver::solve(Statistics& /*statistics*/, const Subproblem& subproblem, double trust_region_radius, + const Vector& /*initial_point*/, Direction& direction, Evaluations& current_evaluations, + const WarmstartInformation& /*warmstart_information*/) { + if (0 < subproblem.number_constraints) { + throw std::runtime_error("TRONSolver cannot solve problems with general constraints"); + } + std::cout << "TRON solver created\n"; + throw std::runtime_error("TRONSolver created"); + } + + SolverWorkspace& TRONSolver::get_workspace() { + return this->workspace; + } + + TRONStats TRONSolver::solve(const HessianOperator& hessian_operator, const std::vector& x0) { + /* TODO + const double atol = (atol > 0) ? atol : std::sqrt(std::numeric_limits::epsilon()); + const double rtol = (rtol > 0) ? rtol : std::sqrt(std::numeric_limits::epsilon()); + */ + + TRONStats stats; + + // ── Initialise x within bounds ─────────────────────────────────── + x_ = x0; + project_bounds(x_); + + auto [fx, gx_new] = fg_(x_); + gx_ = gx_new; + stats.objective = fx; + + // Projected-gradient norm at x0 + project_step(gpx_, x_, gx_, -1.0); // gpx = P(x - g) - x + double pi0 = norm_2(gpx_); + double eps_tol = atol + rtol * pi0; + double fmin = std::min(-1.0, fx) / std::numeric_limits::epsilon(); + + stats.dual_residual = pi0; + + // Trust-region radius + double radius = std::min(std::max(1.0, pi0 / 10.0), max_radius_); + double initial_radius = radius; + + // Cauchy step length (persistent across iterations) + double alpha_c = 1.0; + int num_success = 0; + + bool optimal = (pi0 <= eps_tol); + bool unbounded = (fx < fmin); + + // ── Main loop ──────────────────────────────────────────────────── + while (stats.status == SolveStatus::Unknown) { + if (optimal) { + stats.status = SolveStatus::Optimal; + break; + } + if (unbounded) { + stats.status = SolveStatus::Unbounded; + break; + } + if (stats.iter >= max_iter) { + stats.status = SolveStatus::MaxIter; + break; + } + + // Save current point + xc_ = x_; + double fc = fx; + + // ── Cauchy step ────────────────────────────────────────────── + auto cauchy_status = cauchy_step(hessian_operator, alpha_c, radius); + if (cauchy_status != SolveStatus::Unknown) { + stats.status = cauchy_status; + break; + } + + // ── Projected Newton refinement (CG) ──────────────────────── + std::string cg_info = projected_newton(hessian_operator, radius); + + // ── Ratio test ─────────────────────────────────────────────── + double slope = dot(gx_, s_); + double qs = dot(s_, Hs_) / 2.0 + slope; + + auto [fx_new, dummy_g] = fg_(x_); + (void) dummy_g; + fx = fx_new; + + double ared = fc - fx; + double pred = -qs; + if (pred >= 0.0) { + stats.status = SolveStatus::NegPred; + break; + } + + double ratio = ared / pred; + + if (ratio >= eta1_) { + // accepted + num_success++; + auto [fx2, gx2] = fg_(x_); + gx_ = gx2; + project_step(gpx_, x_, gx_, -1.0); + pi0 = norm_2(gpx_); + } + else { + // rejected + fx = fc; + x_ = xc_; + } + + // ── Update trust-region radius ─────────────────────────────── + double s_norm = norm_2(s_); + if (num_success == 0) + radius = std::min(radius, s_norm); + + if (ratio < eta1_) + radius = std::max(min_radius_, radius / 4.0); + else if (ratio >= eta2_) + radius = std::min(max_radius_, radius * 4.0); + + // ── Update stats ───────────────────────────────────────────── + stats.iter++; + stats.objective = fx; + stats.dual_residual = pi0; + + optimal = (pi0 <= eps_tol); + unbounded = (fx < fmin); + + if (verbose > 0 && stats.iter % verbose == 0) { + std::printf("iter %4d f=%-14.6e π=%-10.3e Δ=%-10.3e %s\n", stats.iter, fx, pi0, radius, cg_info.c_str()); + } + } + if (stats.status == SolveStatus::Unknown) + stats.status = SolveStatus::MaxIter; + + // TODO + // stats.solution = x_; + return stats; + } + + // protected member functions + + // Project v component-wise onto [ℓ, u] + void TRONSolver::project_bounds(Vector &v) const { + for (size_t i = 0; i < v.size(); ++i) + v[i] = std::max(this->lower_bounds[i], std::min(v[i], this->upper_bounds[i])); + } + + /// Active-set indicator: ifix[i] = true if x[i] at a bound AND gradient pushes into it + void TRONSolver::active_set(std::vector& ifix, const Vector& x) const { + for (size_t i = 0; i < x.size(); ++i) { + ifix[i] = (x[i] <= this->lower_bounds[i] || x[i] >= this->upper_bounds[i]); + } + } + + // s = P(x + alpha*d) - x + void TRONSolver::project_step(Vector& s, const Vector& x, const Vector& d, double alpha) const { + for (size_t i = 0; i < x.size(); ++i) { + // TODO dinstiguish the cases + s[i] = std::max(this->lower_bounds[i], std::min(x[i] + alpha * d[i], this->upper_bounds[i])) - x[i]; + } + } + + /// Hs = H*s, slope = gᵀs, qs = ½sᵀHs + gᵀs + std::pair TRONSolver::compute_Hs_slope_qs(const HessianOperator& hessian_operator, const Vector &s, + const Vector &g) { + hessian_operator(s, Hs_); // at xc_ + double slope = dot(g, s); + double qs = dot(s, Hs_) / 2.0 + slope; + return {slope, qs}; + } + + // This subroutine computes the number of break-points, and the minimal and maximal break-points of the projection of + // x + alpha*w on the n-dimensional interval [xl,xu]. + BreakPoints TRONSolver::compute_break_points(const Vector& x, const Vector& w) const { + BreakPoints break_points{0, 0., 0.}; + + for (size_t i = 0; i < x.size(); ++i) { + std::optional break_point = std::nullopt; + if (x[i] < this->upper_bounds[i] && w[i] > 0.) { + break_point = (this->upper_bounds[i] - x[i]) / w[i]; + } + else if (x[i] > this->lower_bounds[i] && w[i] < 0.) { + break_point = (this->lower_bounds[i] - x[i]) / w[i]; + } + if (break_point.has_value()) { + break_points.number++; + if (break_points.number == 1) { + break_points.min = *break_point; + break_points.max = *break_point; + } + else { + break_points.min = std::min(*break_point, break_points.min); + break_points.max = std::max(*break_point, break_points.max); + } + } + } + + // handle the exceptional case. + if (break_points.number == 0) { + break_points.min = 0.; + break_points.max = 0.; + } + return break_points; + } + + /** + * Backtracking projected line search: find smallest t = 2^{-k} s.t. q(s) ≤ μ₀ gᵀs + * where s = P(x + t*d) - x. + * x is updated in-place. + */ + void TRONSolver::projected_line_search(const HessianOperator& hessian_operator, Vector &x, const Vector &d, + const Vector &g) const { + const BreakPoints break_points = compute_break_points(x, d); + double alpha = 1.0; + Vector projected_step(x.size()), Hs(x.size()); + + bool search = true; + while (search && alpha > break_points.min) { + project_step(projected_step, x, d, alpha); + hessian_operator(projected_step, Hs); + double slope = dot(g, projected_step); + double qs = dot(projected_step, Hs) / 2. + slope; + if (qs <= mu0 * slope) { + search = false; + } + else { + alpha /= 2.0; + } + } + if (alpha < std::min(1., break_points.min)) { + alpha = break_points.min; + project_step(projected_step, x, d, alpha); + } + project_step(projected_step, x, d, alpha); + x += projected_step; + project_bounds(x); + + // Update Hs_ for the full step s_ + // (caller recomputes via hv_ after returning) + } + + /** + * Computes s = P(x - α g) - x satisfying sufficient decrease. + * Updates x_ = xc_ (caller is responsible for xc_ being current). + * Modifies s_, Hs_, alpha_c in-place. + */ + SolveStatus TRONSolver::cauchy_step(const HessianOperator& hessian_operator, double& alpha, double radius) { + // Negative gradient direction for breakpoints + for (size_t i = 0; i < xc_.size(); ++i) { + temp_[i] = -gx_[i]; + } + const BreakPoints break_points = compute_break_points(xc_, temp_); + + std::fill(s_.begin(), s_.end(), 0.0); + std::fill(Hs_.begin(), Hs_.end(), 0.0); + + project_step(s_, xc_, gx_, -alpha); + double s_norm = norm_2(s_); + + bool interp; + if (s_norm > mu1 * radius) { + interp = true; + } + else { + auto [slope, qs] = compute_Hs_slope_qs(hessian_operator, s_, gx_); + interp = (qs >= mu0 * slope); + } + + if (interp) { + bool search = true; + while (search) { + alpha /= sigma; + project_step(s_, xc_, gx_, -alpha); + s_norm = norm_2(s_); + if (s_norm <= mu1 * radius) { + auto [slope, qs] = compute_Hs_slope_qs(hessian_operator, s_, gx_); + search = (qs >= mu0 * slope); + } + if (alpha < std::sqrt(std::numeric_limits::min())) + return SolveStatus::SmallStep; + } + } + else { + double alpha_s = alpha; + bool search = true; + while (search && alpha <= break_points.max) { + alpha *= sigma; + project_step(s_, xc_, gx_, -alpha); + s_norm = norm_2(s_); + if (s_norm <= mu1 * radius) { + auto [slope, qs] = compute_Hs_slope_qs(hessian_operator, s_, gx_); + if (qs <= mu0 * slope) + alpha_s = alpha; + } else { + search = false; + } + } + alpha = alpha_s; + project_step(s_, xc_, gx_, -alpha); + } + + // Apply Cauchy step to x_ + for (size_t i = 0; i < x_.size(); ++i) { + x_[i] = xc_[i] + s_[i]; + } + project_bounds(x_); + return SolveStatus::Unknown; + } + + // Scale d in-place so that ‖d + t*p‖ = Delta (solve quadratic for t ≥ 0) + double TRONSolver::compute_distance_to_trust_region(const Vector& d, const Vector& p, double radius) { + const double dd = dot(d, d); + const double dp = dot(d, p); + const double pp = dot(p, p); + double discriminant = dp * dp - pp * (dd - radius * radius); + if (discriminant < 0.0) { + discriminant = 0.0; + } + return (-dp + std::sqrt(discriminant)) / pp; + // axpy(t, p, d); + } + + std::string TRONSolver::projected_newton(const HessianOperator& hessian_operator, double Delta) { + const size_t n = this->x_.size(); + + std::vector ifix(n, false); + Vector rhs(n), r(n), p(n), Hp(n), d(n); + + // Hessian-vector product with free-variable masking + // ZHZ: (Hp)_i = ifix[i] ? 0 : (H * (ifix-masked d))_i + auto masked_hv = [&](const Vector& v, Vector& Hv) { + // zero out fixed components, then apply H + for (size_t i = 0; i < n; ++i) { + d[i] = ifix[i] ? 0.0 : v[i]; + } + hessian_operator(d, Hv); + for (size_t i = 0; i < n; ++i) { + if (ifix[i]) { + Hv[i] = 0.0; + } + } + }; + + // Update Hs_ = H * s_ + hessian_operator(s_, Hs_); + + // x_ = xc_ + s_ projected + for (size_t i = 0; i < n; ++i) { + x_[i] = xc_[i] + s_[i]; + } + project_bounds(x_); + + std::string exit_status = "maximum number of iterations"; + int iters = 0; + + bool exit_optimal = false, exit_pcg = false, exit_itmax = false; + + while (!(exit_optimal || exit_pcg || exit_itmax)) { + active_set(ifix, x_); + int n_free = 0; + for (size_t i = 0; i < n; ++i) { + if (!ifix[i]) { + n_free++; + } + } + if (n_free == 0) { + exit_optimal = true; + continue; + } + + // Build RHS = -(g + Hs) for free variables + double gfnorm = 0.0; + for (size_t i = 0; i < n; ++i) { + rhs[i] = ifix[i] ? 0.0 : -(gx_[i] + Hs_[i]); + gfnorm += rhs[i] * rhs[i]; + } + double gfnorm_sqrt = std::sqrt(gfnorm); + + // ── Conjugate Gradient (Steihaug-Toint style) ──────────────── + d.fill(0.); + r = rhs; + p = r; + double norm_r = norm_2(r); + bool cg_done = false; + std::string cg_flag = "maximum number of iterations"; + + for (int cg_it = 0; cg_it < max_cgiter && !cg_done; ++cg_it) { + masked_hv(p, Hp); + double pHp = dot(p, Hp); + if (pHp <= 0.0) { + // Negative curvature: go to boundary + const double t = compute_distance_to_trust_region(d, p, Delta); + d += t*p; + cg_flag = "on trust-region boundary"; + cg_done = true; + break; + } + double alpha_cg = norm_r / pHp; + + // Trial step + for (size_t i = 0; i < n; ++i) { + w_[i] = d[i] + alpha_cg * p[i]; + } + if (norm_2(w_) >= Delta) { + const double t = compute_distance_to_trust_region(d, p, Delta); + d += t*p; + cg_flag = "on trust-region boundary"; + cg_done = true; + break; + } + d = w_; + + // Update residual r = r - alpha * Hp + // TODO + // r -= alpha_cg*Hp; + double rr_new = dot(r, r); + + if (std::sqrt(rr_new) <= cgtol * gfnorm_sqrt) { + cg_flag = "converged"; + cg_done = true; + break; + } + double beta = rr_new / norm_r; + // p = r + beta*p + for (size_t i = 0; i < n; ++i) { + p[i] = r[i] + beta * p[i]; + } + norm_r = rr_new; + } + + // Projected line search along d from current x_ + // First negate rhs (= g + Hs) for the line-search gradient arg + for (size_t i = 0; i < n; ++i) { + rhs[i] = -rhs[i]; // now = (g+Hs) restricted + } + projected_line_search(hessian_operator, x_, d, rhs); + + // w_ = x_ - (xc_ + s_) → delta_s + for (size_t i = 0; i < n; ++i) { + w_[i] = x_[i] - xc_[i] - s_[i]; + } + s_ += w_; + + hessian_operator(s_, Hs_); + + // Check optimality: ‖(g + Hs) restricted‖ ≤ cgtol * gfnorm_sqrt + double newnorm = 0.0; + for (size_t i = 0; i < n; ++i) { + double ri = ifix[i] ? 0.0 : (gx_[i] + Hs_[i]); + newnorm += ri * ri; + } + if (std::sqrt(newnorm) <= cgtol * gfnorm_sqrt) + exit_optimal = true; + else if (cg_flag == "on trust-region boundary") + exit_pcg = true; + + iters++; + if (iters >= max_cgiter) { + exit_itmax = true; + } + } + + if (exit_optimal) return "stationary point found"; + if (exit_pcg) return "on trust-region boundary"; + if (exit_itmax) return "maximum number of iterations"; + return exit_status; + } +} // namespace \ No newline at end of file diff --git a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp new file mode 100644 index 0000000000..d5f56ee2e7 --- /dev/null +++ b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp @@ -0,0 +1,114 @@ +// Copyright (c) 2026 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNO_TRONSOLVER_H +#define UNO_TRONSOLVER_H + +#include +#include +#include +#include +#include "../SubproblemSolver.hpp" +#include "../SolverWorkspace.hpp" +#include "linear_algebra/Vector.hpp" + +namespace uno { + class TRONSolverWorkspace: public SolverWorkspace { + public: + TRONSolverWorkspace() = default; + + [[nodiscard]] double compute_hessian_quadratic_form(const Subproblem& /*subproblem*/, const Vector& /*vector*/) const override { + return 0.; + } + + Vector objective_gradient; + }; + + enum class SolveStatus { + Unknown, + Optimal, + Unbounded, + MaxIter, + SmallStep, + NegPred, + Error + }; + + struct TRONStats { + SolveStatus status = SolveStatus::Unknown; + int iter = 0; + double objective = 0.0; + double dual_residual = 0.0; ///< ‖P(x - g) - x‖ (projected gradient) + std::vector solution; + }; + + struct BreakPoints { + size_t number; // number of break points + double min; // minimal break-point + double max; // maximal break-point + }; + + class TRONSolver: public SubproblemSolver { + public: + using HessianOperator = std::function& x, Vector& Hv)>; + + TRONSolver() = default; + ~TRONSolver() 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; + /** + * @param x0 Initial guess (length n); need not be feasible. + * @param opts Solver options. + * @return Execution statistics including the solution vector. + */ + TRONStats solve(const HessianOperator& hessian_operator, const std::vector& x0); + + [[nodiscard]] SolverWorkspace& get_workspace() override; + + protected: + std::vector lower_bounds; + std::vector upper_bounds; + TRONSolverWorkspace workspace{}; + + // Workspace vectors (all length n_) + Vector x_, xc_, gx_, gpx_, s_, Hs_, temp_, w_; + + // parameters + double mu0 = 1.0 / 100.0; ///< sufficient-decrease parameter ∈ (0, 0.5) + double mu1 = 1.0; ///< trust-region scaling ∈ (0, ∞) + double sigma = 10.0; ///< step-size update factor ∈ (1, ∞) + // trust-region constants + double eta1_ = 0.1; // acceptance threshold (ratio) + double eta2_ = 0.75; // "very successful" threshold + double min_radius_ = 1e-10; + double max_radius_ = std::min(1.0 / std::sqrt(2.0 * std::numeric_limits::epsilon()), 100.0); + // options + int max_iter = 100000; + int max_cgiter = 50; + double max_time = 30.0; ///< wall-clock seconds + double atol = 0.0; ///< absolute gradient tolerance (set in solve) + double rtol = 0.0; ///< relative gradient tolerance (set in solve) + double cgtol = 0.1; ///< CG sub-problem tolerance + int verbose = 0; + + std::function> (const Vector& x)> fg_; + + void project_bounds(Vector &v) const; + void active_set(std::vector& ifix, const Vector& x) const; + void project_step(Vector& s, const Vector& x, const Vector& d, double alpha) const; + [[nodiscard]] std::pair compute_Hs_slope_qs(const HessianOperator& hessian_operator, + const Vector& s, const Vector &g); + [[nodiscard]] BreakPoints compute_break_points(const Vector& x, const Vector& w) const; + void projected_line_search(const HessianOperator& hessian_operator, Vector &x, const Vector &d, + const Vector &g) const; + [[nodiscard]] SolveStatus cauchy_step(const HessianOperator& hessian_operator, double& alpha, double radius); + [[nodiscard]] static double compute_distance_to_trust_region(const Vector& d, const Vector& p, + double radius) ; + std::string projected_newton(const HessianOperator& hessian_operator, double Delta); + }; +} // namespace + +#endif // UNO_TRONSOLVER_H \ No newline at end of file From 4ca0ccc7672a7401689fae9a971f9374e5e29c6d Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Sat, 18 Apr 2026 17:21:37 +0200 Subject: [PATCH 2/8] More progress --- uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp | 9 +++------ uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp | 1 - 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp index 79c77e0391..4fd043b76f 100644 --- a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp +++ b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp @@ -56,17 +56,17 @@ namespace uno { // Trust-region radius double radius = std::min(std::max(1.0, pi0 / 10.0), max_radius_); - double initial_radius = radius; // Cauchy step length (persistent across iterations) double alpha_c = 1.0; int num_success = 0; - bool optimal = (pi0 <= eps_tol); - bool unbounded = (fx < fmin); + // ── Main loop ──────────────────────────────────────────────────── while (stats.status == SolveStatus::Unknown) { + bool optimal = (pi0 <= eps_tol); + bool unbounded = (fx < fmin); if (optimal) { stats.status = SolveStatus::Optimal; break; @@ -140,9 +140,6 @@ namespace uno { stats.objective = fx; stats.dual_residual = pi0; - optimal = (pi0 <= eps_tol); - unbounded = (fx < fmin); - if (verbose > 0 && stats.iter % verbose == 0) { std::printf("iter %4d f=%-14.6e π=%-10.3e Δ=%-10.3e %s\n", stats.iter, fx, pi0, radius, cg_info.c_str()); } diff --git a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp index d5f56ee2e7..97e36b4a26 100644 --- a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp +++ b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp @@ -6,7 +6,6 @@ #include #include -#include #include #include "../SubproblemSolver.hpp" #include "../SolverWorkspace.hpp" From a804da63c680e302faeb73ac1012f3d453e009b2 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 20 Apr 2026 22:01:31 +0200 Subject: [PATCH 3/8] Good progress: stripped down TRON to single outer iteration --- .../subproblem_solvers/TRON/TRONSolver.cpp | 582 +++++++----------- .../subproblem_solvers/TRON/TRONSolver.hpp | 60 +- uno/symbolic/UnaryNegation.hpp | 4 + 3 files changed, 272 insertions(+), 374 deletions(-) diff --git a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp index 4fd043b76f..d7ecf74698 100644 --- a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp +++ b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp @@ -2,6 +2,7 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include +#include #include #include #include "TRONSolver.hpp" @@ -30,454 +31,335 @@ namespace uno { return this->workspace; } - TRONStats TRONSolver::solve(const HessianOperator& hessian_operator, const std::vector& x0) { - /* TODO - const double atol = (atol > 0) ? atol : std::sqrt(std::numeric_limits::epsilon()); - const double rtol = (rtol > 0) ? rtol : std::sqrt(std::numeric_limits::epsilon()); - */ + void TRONSolver::solve(const ObjectiveOperator& objective_operator, const GradientOperator& gradient_operator, + const MatrixOperator& hessian_operator, const Vector& initial_point) { + SolveStatus status = SolveStatus::Unknown; - TRONStats stats; + // initialize x within bounds + this->x = initial_point; + project_onto_bounds(x); + gradient_operator(this->x, gx); - // ── Initialise x within bounds ─────────────────────────────────── - x_ = x0; - project_bounds(x_); + // projected-gradient norm at x0 + project_step(x, -1., gx, this->gpx); // projected_gradient = P(x - g) - x + double pix = norm_2(this->gpx); + double epsilon = atol + rtol * pix; - auto [fx, gx_new] = fg_(x_); - gx_ = gx_new; - stats.objective = fx; + // trust-region radius + double radius = std::min(std::max(1., pix / 10.), this->max_radius); - // Projected-gradient norm at x0 - project_step(gpx_, x_, gx_, -1.0); // gpx = P(x - g) - x - double pi0 = norm_2(gpx_); - double eps_tol = atol + rtol * pi0; - double fmin = std::min(-1.0, fx) / std::numeric_limits::epsilon(); - - stats.dual_residual = pi0; - - // Trust-region radius - double radius = std::min(std::max(1.0, pi0 / 10.0), max_radius_); - - // Cauchy step length (persistent across iterations) - double alpha_c = 1.0; - int num_success = 0; - - - - // ── Main loop ──────────────────────────────────────────────────── - while (stats.status == SolveStatus::Unknown) { - bool optimal = (pi0 <= eps_tol); - bool unbounded = (fx < fmin); - if (optimal) { - stats.status = SolveStatus::Optimal; - break; - } - if (unbounded) { - stats.status = SolveStatus::Unbounded; - break; - } - if (stats.iter >= max_iter) { - stats.status = SolveStatus::MaxIter; - break; - } - - // Save current point - xc_ = x_; - double fc = fx; - - // ── Cauchy step ────────────────────────────────────────────── - auto cauchy_status = cauchy_step(hessian_operator, alpha_c, radius); - if (cauchy_status != SolveStatus::Unknown) { - stats.status = cauchy_status; - break; - } - - // ── Projected Newton refinement (CG) ──────────────────────── - std::string cg_info = projected_newton(hessian_operator, radius); - - // ── Ratio test ─────────────────────────────────────────────── - double slope = dot(gx_, s_); - double qs = dot(s_, Hs_) / 2.0 + slope; - - auto [fx_new, dummy_g] = fg_(x_); - (void) dummy_g; - fx = fx_new; - - double ared = fc - fx; - double pred = -qs; - if (pred >= 0.0) { - stats.status = SolveStatus::NegPred; - break; - } - - double ratio = ared / pred; - - if (ratio >= eta1_) { - // accepted - num_success++; - auto [fx2, gx2] = fg_(x_); - gx_ = gx2; - project_step(gpx_, x_, gx_, -1.0); - pi0 = norm_2(gpx_); - } - else { - // rejected - fx = fc; - x_ = xc_; - } - - // ── Update trust-region radius ─────────────────────────────── - double s_norm = norm_2(s_); - if (num_success == 0) - radius = std::min(radius, s_norm); - - if (ratio < eta1_) - radius = std::max(min_radius_, radius / 4.0); - else if (ratio >= eta2_) - radius = std::min(max_radius_, radius * 4.0); + if (pix <= epsilon) { + status = SolveStatus::Optimal; + return; + } - // ── Update stats ───────────────────────────────────────────── - stats.iter++; - stats.objective = fx; - stats.dual_residual = pi0; + // TODO change point at which Hessian operator is evaluated - if (verbose > 0 && stats.iter % verbose == 0) { - std::printf("iter %4d f=%-14.6e π=%-10.3e Δ=%-10.3e %s\n", stats.iter, fx, pi0, radius, cg_info.c_str()); - } + // compute Cauchy step and store it in s + double alpha_c = 1.; + const bool cauchy_success = compute_cauchy_step(hessian_operator, this->gx, alpha_c, radius); + if (!cauchy_success) { + status = SolveStatus::Error; + return; } - if (stats.status == SolveStatus::Unknown) - stats.status = SolveStatus::MaxIter; - // TODO - // stats.solution = x_; - return stats; + // projected Newton refinement (CG) + std::string cg_info = projected_newton(hessian_operator, gx, radius); } // protected member functions - // Project v component-wise onto [ℓ, u] - void TRONSolver::project_bounds(Vector &v) const { + // project v component-wise onto [ℓ, u] + void TRONSolver::project_onto_bounds(Vector& v) const { for (size_t i = 0; i < v.size(); ++i) v[i] = std::max(this->lower_bounds[i], std::min(v[i], this->upper_bounds[i])); } - /// Active-set indicator: ifix[i] = true if x[i] at a bound AND gradient pushes into it - void TRONSolver::active_set(std::vector& ifix, const Vector& x) const { + // active-set indicator: active[i] = true if x[i] at a bound + void TRONSolver::compute_active_set(std::vector& active, const Vector& x) const { for (size_t i = 0; i < x.size(); ++i) { - ifix[i] = (x[i] <= this->lower_bounds[i] || x[i] >= this->upper_bounds[i]); + const double delta = (-INF < this->lower_bounds[i] && this->lower_bounds[i] < this->upper_bounds[i] && + this->upper_bounds[i] < INF) ? std::min(atol, rtol * (this->upper_bounds[i] - this->lower_bounds[i])) : atol; + if (x[i] == this->lower_bounds[i] && x[i] == this->upper_bounds[i]) { + active[i] = true; + } + else if (x[i] <= this->lower_bounds[i] + delta) { + active[i] = true; + } + else if (x[i] >= this->upper_bounds[i] - delta) { + active[i] = true; + } + else { + active[i] = false; + } } } // s = P(x + alpha*d) - x - void TRONSolver::project_step(Vector& s, const Vector& x, const Vector& d, double alpha) const { + void TRONSolver::project_step(const Vector& x, double alpha, const Vector& d, Vector& s) const { for (size_t i = 0; i < x.size(); ++i) { - // TODO dinstiguish the cases - s[i] = std::max(this->lower_bounds[i], std::min(x[i] + alpha * d[i], this->upper_bounds[i])) - x[i]; + if (x[i] + alpha * d[i] < this->lower_bounds[i]) { + s[i] = this->lower_bounds[i] - x[i]; + } + else if (this->upper_bounds[i] < x[i] + alpha * d[i]) { + s[i] = this->upper_bounds[i] - x[i]; + } + else { + s[i] = alpha * d[i]; + } } } - /// Hs = H*s, slope = gᵀs, qs = ½sᵀHs + gᵀs - std::pair TRONSolver::compute_Hs_slope_qs(const HessianOperator& hessian_operator, const Vector &s, - const Vector &g) { - hessian_operator(s, Hs_); // at xc_ + // Hs = H*s, slope = gᵀs, qs = ½sᵀHs + gᵀs + std::pair TRONSolver::compute_Hs_slope_qs(const MatrixOperator& hessian_operator, const Vector& s, + const Vector& g) { + hessian_operator(s, Hs); // at xc_ double slope = dot(g, s); - double qs = dot(s, Hs_) / 2.0 + slope; + double qs = dot(s, Hs) / 2. + slope; return {slope, qs}; } - // This subroutine computes the number of break-points, and the minimal and maximal break-points of the projection of - // x + alpha*w on the n-dimensional interval [xl,xu]. - BreakPoints TRONSolver::compute_break_points(const Vector& x, const Vector& w) const { - BreakPoints break_points{0, 0., 0.}; + // compute the minimal and maximal break-points of the projection of x + alpha*d on the n-dimensional interval [xl,xu]. + BreakPoints TRONSolver::compute_break_points(const Vector& x, const Vector& d) const { + BreakPoints break_points{INF, 0.}; for (size_t i = 0; i < x.size(); ++i) { std::optional break_point = std::nullopt; - if (x[i] < this->upper_bounds[i] && w[i] > 0.) { - break_point = (this->upper_bounds[i] - x[i]) / w[i]; + if (x[i] < this->upper_bounds[i] && d[i] > 0.) { + break_point = (this->upper_bounds[i] - x[i]) / d[i]; } - else if (x[i] > this->lower_bounds[i] && w[i] < 0.) { - break_point = (this->lower_bounds[i] - x[i]) / w[i]; + else if (x[i] > this->lower_bounds[i] && d[i] < 0.) { + break_point = (this->lower_bounds[i] - x[i]) / d[i]; } if (break_point.has_value()) { - break_points.number++; - if (break_points.number == 1) { - break_points.min = *break_point; - break_points.max = *break_point; - } - else { - break_points.min = std::min(*break_point, break_points.min); - break_points.max = std::max(*break_point, break_points.max); - } + break_points.min = std::min(*break_point, break_points.min); + break_points.max = std::max(*break_point, break_points.max); } } - - // handle the exceptional case. - if (break_points.number == 0) { - break_points.min = 0.; - break_points.max = 0.; - } return break_points; } - /** - * Backtracking projected line search: find smallest t = 2^{-k} s.t. q(s) ≤ μ₀ gᵀs - * where s = P(x + t*d) - x. - * x is updated in-place. - */ - void TRONSolver::projected_line_search(const HessianOperator& hessian_operator, Vector &x, const Vector &d, - const Vector &g) const { + // backtracking projected line search: find smallest t = 2^{-k} s.t. q(s) ≤ μ₀ gᵀs, where s = P(x + t*d) - x. + void TRONSolver::projected_line_search(const MatrixOperator& hessian_operator, Vector& x, const Vector& d, + const Vector& g, Vector& s) { + double alpha = 1.; const BreakPoints break_points = compute_break_points(x, d); - double alpha = 1.0; - Vector projected_step(x.size()), Hs(x.size()); bool search = true; while (search && alpha > break_points.min) { - project_step(projected_step, x, d, alpha); - hessian_operator(projected_step, Hs); - double slope = dot(g, projected_step); - double qs = dot(projected_step, Hs) / 2. + slope; + project_step(x, alpha, d, s); + const auto [slope, qs] = compute_Hs_slope_qs(hessian_operator, s, g); if (qs <= mu0 * slope) { search = false; } else { - alpha /= 2.0; + alpha /= 2.; } } if (alpha < std::min(1., break_points.min)) { alpha = break_points.min; - project_step(projected_step, x, d, alpha); + project_step(x, alpha, d, s); + hessian_operator(s, this->Hs); } - project_step(projected_step, x, d, alpha); - x += projected_step; - project_bounds(x); - - // Update Hs_ for the full step s_ - // (caller recomputes via hv_ after returning) + project_step(x, alpha, d, s); + x += s; } - /** - * Computes s = P(x - α g) - x satisfying sufficient decrease. - * Updates x_ = xc_ (caller is responsible for xc_ being current). - * Modifies s_, Hs_, alpha_c in-place. - */ - SolveStatus TRONSolver::cauchy_step(const HessianOperator& hessian_operator, double& alpha, double radius) { - // Negative gradient direction for breakpoints - for (size_t i = 0; i < xc_.size(); ++i) { - temp_[i] = -gx_[i]; - } - const BreakPoints break_points = compute_break_points(xc_, temp_); + // compute Cauchy step s = P(x - α g) - x satisfying sufficient decrease. + // returns true upon success, false upon failure + bool TRONSolver::compute_cauchy_step(const MatrixOperator& hessian_operator, const Vector& g, double& alpha, + double radius) { + // compute breakpoints along negative gradient direction + s = -g; + const BreakPoints break_points = compute_break_points(x, s); - std::fill(s_.begin(), s_.end(), 0.0); - std::fill(Hs_.begin(), Hs_.end(), 0.0); + s.fill(0.); + Hs.fill(0.); - project_step(s_, xc_, gx_, -alpha); - double s_norm = norm_2(s_); + project_step(x, -alpha, g, s); - bool interp; - if (s_norm > mu1 * radius) { - interp = true; - } - else { - auto [slope, qs] = compute_Hs_slope_qs(hessian_operator, s_, gx_); - interp = (qs >= mu0 * slope); + // interpolate or extrapolate + bool interpolate = true; + if (norm_2(s) <= mu1 * radius) { + auto [slope, qs] = compute_Hs_slope_qs(hessian_operator, s, g); + interpolate = (qs >= mu0 * slope); } - if (interp) { + if (interpolate) { bool search = true; while (search) { alpha /= sigma; - project_step(s_, xc_, gx_, -alpha); - s_norm = norm_2(s_); - if (s_norm <= mu1 * radius) { - auto [slope, qs] = compute_Hs_slope_qs(hessian_operator, s_, gx_); + project_step(x, -alpha, g, s); + if (norm_2(s) <= mu1 * radius) { + auto [slope, qs] = compute_Hs_slope_qs(hessian_operator, s, g); search = (qs >= mu0 * slope); } - if (alpha < std::sqrt(std::numeric_limits::min())) - return SolveStatus::SmallStep; + // TODO: correctly assess why this fails + if (alpha < std::sqrt(std::nextafter(0., 1.))) { + return false; + } } } - else { - double alpha_s = alpha; + else { // extrapolation + double alpha_success = alpha; bool search = true; while (search && alpha <= break_points.max) { alpha *= sigma; - project_step(s_, xc_, gx_, -alpha); - s_norm = norm_2(s_); - if (s_norm <= mu1 * radius) { - auto [slope, qs] = compute_Hs_slope_qs(hessian_operator, s_, gx_); - if (qs <= mu0 * slope) - alpha_s = alpha; - } else { + project_step(x, -alpha, g, s); + if (norm_2(s) <= mu1 * radius) { + auto [slope, qs] = compute_Hs_slope_qs(hessian_operator, s, g); + if (qs <= mu0 * slope) { + alpha_success = alpha; + } + } + else { search = false; } } - alpha = alpha_s; - project_step(s_, xc_, gx_, -alpha); + // recover the last successful step + alpha = alpha_success; + project_step(x, -alpha, g, s); } - - // Apply Cauchy step to x_ - for (size_t i = 0; i < x_.size(); ++i) { - x_[i] = xc_[i] + s_[i]; - } - project_bounds(x_); - return SolveStatus::Unknown; + return true; } - // Scale d in-place so that ‖d + t*p‖ = Delta (solve quadratic for t ≥ 0) + // find t ≥ 0 so that ‖d + t*p‖ = Delta (solve quadratic for t ≥ 0) double TRONSolver::compute_distance_to_trust_region(const Vector& d, const Vector& p, double radius) { - const double dd = dot(d, d); - const double dp = dot(d, p); - const double pp = dot(p, p); - double discriminant = dp * dp - pp * (dd - radius * radius); + const double dTd = dot(d, d); + const double dTp = dot(d, p); + const double pTp = dot(p, p); + double discriminant = dTp * dTp - pTp * (dTd - radius * radius); if (discriminant < 0.0) { discriminant = 0.0; } - return (-dp + std::sqrt(discriminant)) / pp; - // axpy(t, p, d); + return (-dTp + std::sqrt(discriminant)) / pTp; } - std::string TRONSolver::projected_newton(const HessianOperator& hessian_operator, double Delta) { - const size_t n = this->x_.size(); - - std::vector ifix(n, false); - Vector rhs(n), r(n), p(n), Hp(n), d(n); + // Conjugate Gradient (Steihaug-Toint style) + CGStatus TRONSolver::CG(Vector& d, const MatrixOperator& matrix_operator, const Vector& rhs, double radius, + double gfnorm_sqrt) { + const size_t n = d.size(); + Vector r(n), p(n), Hp(n); + d.fill(0.); + r = rhs; + p = r; + double norm_r = norm_2(r); + + for (size_t cg_it = 0; cg_it < max_cgiter; ++cg_it) { + matrix_operator(p, Hp); + double pHp = dot(p, Hp); + if (pHp <= 0.0) { + // Negative curvature: go to boundary + const double t = compute_distance_to_trust_region(d, p, radius); + d += t*p; + return CGStatus::ON_TR_BOUNDARY; + } + double alpha_cg = norm_r / pHp; - // Hessian-vector product with free-variable masking - // ZHZ: (Hp)_i = ifix[i] ? 0 : (H * (ifix-masked d))_i - auto masked_hv = [&](const Vector& v, Vector& Hv) { - // zero out fixed components, then apply H - for (size_t i = 0; i < n; ++i) { - d[i] = ifix[i] ? 0.0 : v[i]; - } - hessian_operator(d, Hv); - for (size_t i = 0; i < n; ++i) { - if (ifix[i]) { - Hv[i] = 0.0; - } - } - }; + // trial step + for (size_t i = 0; i < n; ++i) { + w[i] = d[i] + alpha_cg * p[i]; + } + if (norm_2(w) >= radius) { + const double t = compute_distance_to_trust_region(d, p, radius); + d += t*p; + return CGStatus::ON_TR_BOUNDARY; + } + d = w; - // Update Hs_ = H * s_ - hessian_operator(s_, Hs_); + // Update residual r = r - alpha * Hp + // TODO + // r -= alpha_cg*Hp; + double rr_new = dot(r, r); - // x_ = xc_ + s_ projected + if (std::sqrt(rr_new) <= cgtol * gfnorm_sqrt) { + return CGStatus::SUCCESS; + } + double beta = rr_new / norm_r; + // p = r + beta*p for (size_t i = 0; i < n; ++i) { - x_[i] = xc_[i] + s_[i]; + p[i] = r[i] + beta * p[i]; } - project_bounds(x_); + norm_r = rr_new; + } + return CGStatus::MAX_ITERATIONS; + } - std::string exit_status = "maximum number of iterations"; - int iters = 0; + std::string TRONSolver::projected_newton(const MatrixOperator& hessian_operator, const Vector& g, double radius) { + const size_t n = this->x.size(); + + // Update Hs = H * s + hessian_operator(s, Hs); + + // projected Newton step + bool exit_optimal = false, exit_pcg = false, exit_itmax = false; + std::string exit_status = "maximum number of iterations"; + size_t iters = 0; + x += s; + project_onto_bounds(x); + std::vector active(n, false); + while (!exit_optimal && !exit_pcg && !exit_itmax) { + compute_active_set(active, x); + // stop if all bounds are active + if (std::all_of(active.begin(), active.end(), [](bool active) { return active; })) { + exit_optimal = true; + continue; + } - bool exit_optimal = false, exit_pcg = false, exit_itmax = false; + // build RHS = -(g + Hs) for free variables + double gfnorm = 0.0; + for (size_t i = 0; i < n; ++i) { + this->quadratic_gradient[i] = active[i] ? 0. : -g[i]; + gfnorm += this->quadratic_gradient[i] * this->quadratic_gradient[i]; + this->quadratic_gradient[i] -= active[i] ? 0. : Hs[i]; + } + double gfnorm_sqrt = std::sqrt(gfnorm); - while (!(exit_optimal || exit_pcg || exit_itmax)) { - active_set(ifix, x_); - int n_free = 0; + // define ZHZ, the Hessian-vector product with free-variable masking + // (Hp)_i = ifix[i] ? 0 : (H * (ifix-masked d))_i + const auto ZHZ = [&](const Vector& d, Vector& result) { + // zero out fixed components of d, then apply H for (size_t i = 0; i < n; ++i) { - if (!ifix[i]) { - n_free++; - } - } - if (n_free == 0) { - exit_optimal = true; - continue; + this->d_masked[i] = active[i] ? 0. : d[i]; } - - // Build RHS = -(g + Hs) for free variables - double gfnorm = 0.0; + hessian_operator(this->d_masked, result); + // zero out fixed components of the result for (size_t i = 0; i < n; ++i) { - rhs[i] = ifix[i] ? 0.0 : -(gx_[i] + Hs_[i]); - gfnorm += rhs[i] * rhs[i]; - } - double gfnorm_sqrt = std::sqrt(gfnorm); - - // ── Conjugate Gradient (Steihaug-Toint style) ──────────────── - d.fill(0.); - r = rhs; - p = r; - double norm_r = norm_2(r); - bool cg_done = false; - std::string cg_flag = "maximum number of iterations"; - - for (int cg_it = 0; cg_it < max_cgiter && !cg_done; ++cg_it) { - masked_hv(p, Hp); - double pHp = dot(p, Hp); - if (pHp <= 0.0) { - // Negative curvature: go to boundary - const double t = compute_distance_to_trust_region(d, p, Delta); - d += t*p; - cg_flag = "on trust-region boundary"; - cg_done = true; - break; + if (active[i]) { + result[i] = 0.; } - double alpha_cg = norm_r / pHp; - - // Trial step - for (size_t i = 0; i < n; ++i) { - w_[i] = d[i] + alpha_cg * p[i]; - } - if (norm_2(w_) >= Delta) { - const double t = compute_distance_to_trust_region(d, p, Delta); - d += t*p; - cg_flag = "on trust-region boundary"; - cg_done = true; - break; - } - d = w_; - - // Update residual r = r - alpha * Hp - // TODO - // r -= alpha_cg*Hp; - double rr_new = dot(r, r); - - if (std::sqrt(rr_new) <= cgtol * gfnorm_sqrt) { - cg_flag = "converged"; - cg_done = true; - break; - } - double beta = rr_new / norm_r; - // p = r + beta*p - for (size_t i = 0; i < n; ++i) { - p[i] = r[i] + beta * p[i]; - } - norm_r = rr_new; } + }; - // Projected line search along d from current x_ - // First negate rhs (= g + Hs) for the line-search gradient arg - for (size_t i = 0; i < n; ++i) { - rhs[i] = -rhs[i]; // now = (g+Hs) restricted - } - projected_line_search(hessian_operator, x_, d, rhs); + const CGStatus cg_status = CG(this->d_cg, ZHZ, this->quadratic_gradient, radius, gfnorm_sqrt); + iters++; - // w_ = x_ - (xc_ + s_) → delta_s - for (size_t i = 0; i < n; ++i) { - w_[i] = x_[i] - xc_[i] - s_[i]; - } - s_ += w_; + // projected line search + this->quadratic_gradient.scale(-1.); + projected_line_search(ZHZ, x, d_cg, this->quadratic_gradient, w); + s += w; + hessian_operator(s, Hs); - hessian_operator(s_, Hs_); - - // Check optimality: ‖(g + Hs) restricted‖ ≤ cgtol * gfnorm_sqrt - double newnorm = 0.0; - for (size_t i = 0; i < n; ++i) { - double ri = ifix[i] ? 0.0 : (gx_[i] + Hs_[i]); - newnorm += ri * ri; - } - if (std::sqrt(newnorm) <= cgtol * gfnorm_sqrt) - exit_optimal = true; - else if (cg_flag == "on trust-region boundary") - exit_pcg = true; - - iters++; - if (iters >= max_cgiter) { - exit_itmax = true; - } + // Check optimality: ‖(g + Hs) restricted‖ ≤ cgtol * gfnorm_sqrt + double new_norm = 0.0; + for (size_t i = 0; i < n; ++i) { + const double ri = active[i] ? 0. : Hs[i] + g[i]; + new_norm += ri * ri; + } + if (std::sqrt(new_norm) <= cgtol * gfnorm_sqrt) { + exit_optimal = true; + } + else if (cg_status == CGStatus::ON_TR_BOUNDARY) { + exit_pcg = true; + } + else if (iters >= max_cgiter) { + exit_itmax = true; } - - if (exit_optimal) return "stationary point found"; - if (exit_pcg) return "on trust-region boundary"; - if (exit_itmax) return "maximum number of iterations"; - return exit_status; } + + if (exit_optimal) return "stationary point found"; + if (exit_pcg) return "on trust-region boundary"; + if (exit_itmax) return "maximum number of iterations"; + return exit_status; + } } // namespace \ No newline at end of file diff --git a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp index 97e36b4a26..408746848d 100644 --- a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp +++ b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.hpp @@ -26,13 +26,18 @@ namespace uno { enum class SolveStatus { Unknown, Optimal, - Unbounded, MaxIter, SmallStep, NegPred, Error }; + enum class CGStatus { + SUCCESS, + ON_TR_BOUNDARY, + MAX_ITERATIONS + }; + struct TRONStats { SolveStatus status = SolveStatus::Unknown; int iter = 0; @@ -42,14 +47,15 @@ namespace uno { }; struct BreakPoints { - size_t number; // number of break points double min; // minimal break-point double max; // maximal break-point }; class TRONSolver: public SubproblemSolver { public: - using HessianOperator = std::function& x, Vector& Hv)>; + using ObjectiveOperator = std::function& x)>; + using GradientOperator = std::function& x, Vector& gradient)>; + using MatrixOperator = std::function& x, Vector& Hv)>; TRONSolver() = default; ~TRONSolver() override = default; @@ -59,11 +65,12 @@ namespace uno { 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; /** - * @param x0 Initial guess (length n); need not be feasible. + * @param x Initial guess (length n); need not be feasible. * @param opts Solver options. * @return Execution statistics including the solution vector. */ - TRONStats solve(const HessianOperator& hessian_operator, const std::vector& x0); + void solve(const ObjectiveOperator& objective_operator, const GradientOperator& gradient_operator, + const MatrixOperator& hessian_operator, const Vector& initial_point); [[nodiscard]] SolverWorkspace& get_workspace() override; @@ -73,7 +80,11 @@ namespace uno { TRONSolverWorkspace workspace{}; // Workspace vectors (all length n_) - Vector x_, xc_, gx_, gpx_, s_, Hs_, temp_, w_; + Vector x, xc, x_copy, gpx, gx, s, Hs, temp_, w; + Vector d_masked; + Vector quadratic_gradient; + Vector d_cg; // CG point + double fc; // parameters double mu0 = 1.0 / 100.0; ///< sufficient-decrease parameter ∈ (0, 0.5) @@ -83,30 +94,31 @@ namespace uno { double eta1_ = 0.1; // acceptance threshold (ratio) double eta2_ = 0.75; // "very successful" threshold double min_radius_ = 1e-10; - double max_radius_ = std::min(1.0 / std::sqrt(2.0 * std::numeric_limits::epsilon()), 100.0); + double max_radius = std::min(1.0 / std::sqrt(2.0 * std::numeric_limits::epsilon()), 100.0); // options - int max_iter = 100000; - int max_cgiter = 50; - double max_time = 30.0; ///< wall-clock seconds - double atol = 0.0; ///< absolute gradient tolerance (set in solve) - double rtol = 0.0; ///< relative gradient tolerance (set in solve) + size_t max_iter = 100000; + size_t max_cgiter = 50; double cgtol = 0.1; ///< CG sub-problem tolerance - int verbose = 0; + const double atol = std::sqrt(std::numeric_limits::epsilon()); + const double rtol = std::sqrt(std::numeric_limits::epsilon()); std::function> (const Vector& x)> fg_; - void project_bounds(Vector &v) const; - void active_set(std::vector& ifix, const Vector& x) const; - void project_step(Vector& s, const Vector& x, const Vector& d, double alpha) const; - [[nodiscard]] std::pair compute_Hs_slope_qs(const HessianOperator& hessian_operator, - const Vector& s, const Vector &g); - [[nodiscard]] BreakPoints compute_break_points(const Vector& x, const Vector& w) const; - void projected_line_search(const HessianOperator& hessian_operator, Vector &x, const Vector &d, - const Vector &g) const; - [[nodiscard]] SolveStatus cauchy_step(const HessianOperator& hessian_operator, double& alpha, double radius); + void project_onto_bounds(Vector& v) const; + void compute_active_set(std::vector& active, const Vector& x) const; + void project_step(const Vector& x, double alpha, const Vector& d, Vector& s) const; + [[nodiscard]] std::pair compute_Hs_slope_qs(const MatrixOperator& hessian_operator, + const Vector& s, const Vector& g); + [[nodiscard]] BreakPoints compute_break_points(const Vector& x, const Vector& d) const; + void projected_line_search(const MatrixOperator& hessian_operator, Vector& x, const Vector& d, + const Vector& g, Vector& s); + [[nodiscard]] bool compute_cauchy_step(const MatrixOperator& hessian_operator, const Vector& g, double& alpha, + double radius); [[nodiscard]] static double compute_distance_to_trust_region(const Vector& d, const Vector& p, - double radius) ; - std::string projected_newton(const HessianOperator& hessian_operator, double Delta); + double radius); + [[nodiscard]] CGStatus CG(Vector& d, const MatrixOperator& matrix_operator, const Vector& rhs, + double radius, double gfnorm_sqrt); + std::string projected_newton(const MatrixOperator& hessian_operator, const Vector& g, double radius); }; } // namespace diff --git a/uno/symbolic/UnaryNegation.hpp b/uno/symbolic/UnaryNegation.hpp index 09994fde69..4697fcf95e 100644 --- a/uno/symbolic/UnaryNegation.hpp +++ b/uno/symbolic/UnaryNegation.hpp @@ -26,6 +26,10 @@ namespace uno { return -this->expression[index]; } + [[nodiscard]] const Expression& get_expression() const { + return this->expression; + } + protected: storage_t expression; }; From 1580eb4b6e0fee03a437a4d2c6fc13d864220c42 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Wed, 29 Apr 2026 15:05:59 +0200 Subject: [PATCH 4/8] Fixed detection of bound-constrained problems --- .../InequalityHandlingMethodFactory.cpp | 2 +- uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp index 9d9fb6dc64..c53c7399ec 100644 --- a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp +++ b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp @@ -29,7 +29,7 @@ namespace uno { } // from now on, the problem has inequalities // bound-constrained problem: do not reformulate - if (!problem.has_inequality_constraints() && problem.has_bound_constraints()) { + if (problem.number_constraints == 0 && problem.has_bound_constraints()) { return std::make_unique("bound-constrained method"); } diff --git a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp b/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp index 332554be46..79889ee6c0 100644 --- a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp +++ b/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp @@ -54,7 +54,7 @@ namespace uno { } } // if only bound constraints, allocate bound-constrained solver - else if (!subproblem.has_inequality_constraints() && subproblem.has_bound_constraints()) { + else if (subproblem.number_constraints == 0 && subproblem.has_bound_constraints()) { auto subproblem_solver = std::make_unique(); subproblem_solver->initialize_memory(subproblem); return subproblem_solver; From 632da9d3c9823d7529d8e4ef5a7fdf58ae2ddced Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Tue, 26 May 2026 19:18:20 +0200 Subject: [PATCH 5/8] Fixed includes --- uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp | 2 ++ uno/symbolic/UnaryNegation.hpp | 4 +--- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp index d7ecf74698..d191abfdc8 100644 --- a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp +++ b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp @@ -2,6 +2,7 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include +#include #include #include #include @@ -9,6 +10,7 @@ #include "ingredients/subproblem/Subproblem.hpp" #include "linear_algebra/BLAS.hpp" #include "linear_algebra/Vector.hpp" +#include "symbolic/UnaryNegation.hpp" namespace uno { void TRONSolver::initialize_memory(const Subproblem& subproblem) { diff --git a/uno/symbolic/UnaryNegation.hpp b/uno/symbolic/UnaryNegation.hpp index 4697fcf95e..7f6baba051 100644 --- a/uno/symbolic/UnaryNegation.hpp +++ b/uno/symbolic/UnaryNegation.hpp @@ -26,9 +26,7 @@ namespace uno { return -this->expression[index]; } - [[nodiscard]] const Expression& get_expression() const { - return this->expression; - } + UNO_FORWARD_ACCESSOR(get_expression, this->expression) protected: storage_t expression; From f976e2c6e9d72006e095c1689c44107d54cf2101 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Tue, 26 May 2026 19:18:29 +0200 Subject: [PATCH 6/8] Disabled TRON in SubproblemSolverFactory --- uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp b/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp index 79889ee6c0..ae0d314944 100644 --- a/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp +++ b/uno/ingredients/subproblem_solvers/SubproblemSolverFactory.hpp @@ -53,12 +53,14 @@ namespace uno { return subproblem_solver; } } + /* // if only bound constraints, allocate bound-constrained solver else if (subproblem.number_constraints == 0 && subproblem.has_bound_constraints()) { auto subproblem_solver = std::make_unique(); subproblem_solver->initialize_memory(subproblem); return subproblem_solver; } + */ // if no inequality constraint and no trust region, allocate EQP solver else if (!subproblem.has_inequality_constraints() && !subproblem.has_bound_constraints() && !uses_trust_region) { if constexpr (std::is_same_v) { // unconstrained From 6cb3cda2c14c624cc2ea9d9cfb8637999c3ce109 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Tue, 26 May 2026 21:54:42 +0200 Subject: [PATCH 7/8] Disable TRON-related decision in InequalityHandlingMethodFactory.cpp --- .../InequalityHandlingMethodFactory.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp index c53c7399ec..a55ade23c9 100644 --- a/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp +++ b/uno/ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.cpp @@ -28,10 +28,12 @@ namespace uno { } } // from now on, the problem has inequalities + /* // bound-constrained problem: do not reformulate if (problem.number_constraints == 0 && problem.has_bound_constraints()) { return std::make_unique("bound-constrained method"); } + */ const std::string inequality_handling_method = options.get_string("inequality_handling_method"); // inequality-constrained methods From 4516d1407b181ba12e3607425ce18882bca5f981 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Tue, 26 May 2026 22:02:48 +0200 Subject: [PATCH 8/8] Missing include --- uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp index d191abfdc8..8ac9c41b68 100644 --- a/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp +++ b/uno/ingredients/subproblem_solvers/TRON/TRONSolver.cpp @@ -1,6 +1,7 @@ // Copyright (c) 2026 Charlie Vanaret // Licensed under the MIT license. See LICENSE file in the project directory for details. +#include #include #include #include