Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Reduce the number of symbolic analyses #104

Merged
merged 11 commits into from
Dec 5, 2024
4 changes: 2 additions & 2 deletions .github/julia/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -18,11 +18,11 @@ import Uno_jll
Create a new `AmplNLWriter.Optimizer` object that uses Uno as the backing
solver.
"""
function Optimizer(options = String["logger=SILENT"])
function Optimizer(options = String["logger=SILENT", "unbounded_objective_threshold=-1e15"])
return AmplNLWriter.Optimizer(Uno_jll.amplexe, options)
end

Optimizer_LP() = Optimizer(["logger=SILENT", "preset=filterslp", "max_iterations=10000"])
Optimizer_LP() = Optimizer(["logger=SILENT", "preset=filterslp", "max_iterations=10000", "unbounded_objective_threshold=-1e15"])

# This testset runs https://github.com/jump-dev/MINLPTests.jl
@testset "MINLPTests" begin
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -3,12 +3,13 @@

#include <cmath>
#include "PrimalDualInteriorPointSubproblem.hpp"
#include "optimization/Direction.hpp"
#include "optimization/Iterate.hpp"
#include "ingredients/hessian_models/HessianModelFactory.hpp"
#include "linear_algebra/SparseStorageFactory.hpp"
#include "linear_algebra/SymmetricIndefiniteLinearSystem.hpp"
#include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
#include "solvers/SymmetricIndefiniteLinearSolverFactory.hpp"
#include "optimization/Direction.hpp"
#include "optimization/Iterate.hpp"
#include "optimization/WarmstartInformation.hpp"
#include "preprocessing/Preprocessing.hpp"
#include "reformulation/l1RelaxedProblem.hpp"
Expand Down Expand Up @@ -123,26 +124,6 @@ namespace uno {

void PrimalDualInteriorPointSubproblem::evaluate_functions(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate,
const Multipliers& current_multipliers, const WarmstartInformation& warmstart_information) {
// barrier Lagrangian Hessian
if (warmstart_information.objective_changed || warmstart_information.constraints_changed) {
// original Lagrangian Hessian
this->hessian_model->evaluate(statistics, problem, current_iterate.primals, current_multipliers.constraints);

// diagonal barrier terms (grouped by variable)
for (size_t variable_index: Range(problem.number_variables)) {
double diagonal_barrier_term = 0.;
if (is_finite(problem.variable_lower_bound(variable_index))) { // lower bounded
const double distance_to_bound = current_iterate.primals[variable_index] - problem.variable_lower_bound(variable_index);
diagonal_barrier_term += current_multipliers.lower_bounds[variable_index] / distance_to_bound;
}
if (is_finite(problem.variable_upper_bound(variable_index))) { // upper bounded
const double distance_to_bound = current_iterate.primals[variable_index] - problem.variable_upper_bound(variable_index);
diagonal_barrier_term += current_multipliers.upper_bounds[variable_index] / distance_to_bound;
}
this->hessian_model->hessian.insert(diagonal_barrier_term, variable_index, variable_index);
}
}

// barrier objective gradient
if (warmstart_information.objective_changed) {
// original objective gradient
Expand Down Expand Up @@ -174,6 +155,26 @@ namespace uno {
problem.evaluate_constraints(current_iterate, this->constraints);
problem.evaluate_constraint_jacobian(current_iterate, this->constraint_jacobian);
}

// barrier Lagrangian Hessian
if (warmstart_information.objective_changed || warmstart_information.constraints_changed) {
// original Lagrangian Hessian
this->hessian_model->evaluate(statistics, problem, current_iterate.primals, current_multipliers.constraints);

// diagonal barrier terms (grouped by variable)
for (size_t variable_index: Range(problem.number_variables)) {
double diagonal_barrier_term = 0.;
if (is_finite(problem.variable_lower_bound(variable_index))) { // lower bounded
const double distance_to_bound = current_iterate.primals[variable_index] - problem.variable_lower_bound(variable_index);
diagonal_barrier_term += current_multipliers.lower_bounds[variable_index] / distance_to_bound;
}
if (is_finite(problem.variable_upper_bound(variable_index))) { // upper bounded
const double distance_to_bound = current_iterate.primals[variable_index] - problem.variable_upper_bound(variable_index);
diagonal_barrier_term += current_multipliers.upper_bounds[variable_index] / distance_to_bound;
}
this->hessian_model->hessian.insert(diagonal_barrier_term, variable_index, variable_index);
}
}
}

void PrimalDualInteriorPointSubproblem::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate,
Expand All @@ -199,7 +200,7 @@ namespace uno {
this->evaluate_functions(statistics, problem, current_iterate, current_multipliers, warmstart_information);

// compute the primal-dual solution
this->assemble_augmented_system(statistics, problem, current_multipliers);
this->assemble_augmented_system(statistics, problem, current_multipliers, warmstart_information);
this->augmented_system.solve(*this->linear_solver);
assert(direction.status == SubproblemStatus::OPTIMAL && "The primal-dual perturbed subproblem was not solved to optimality");
this->number_subproblems_solved++;
Expand All @@ -209,13 +210,14 @@ namespace uno {
}

void PrimalDualInteriorPointSubproblem::assemble_augmented_system(Statistics& statistics, const OptimizationProblem& problem,
const Multipliers& current_multipliers) {
const Multipliers& current_multipliers, const WarmstartInformation& warmstart_information) {
// assemble, factorize and regularize the augmented matrix
this->augmented_system.assemble_matrix(this->hessian_model->hessian, this->constraint_jacobian, problem.number_variables, problem.number_constraints);
this->augmented_system.factorize_matrix(problem.model, *this->linear_solver);
DEBUG << "Testing factorization with regularization factors (0, 0)\n";
this->augmented_system.factorize_matrix(*this->linear_solver, warmstart_information);
const double dual_regularization_parameter = std::pow(this->barrier_parameter(), this->parameters.regularization_exponent);
this->augmented_system.regularize_matrix(statistics, problem.model, *this->linear_solver, problem.number_variables, problem.number_constraints,
dual_regularization_parameter);
this->augmented_system.regularize_matrix(statistics, *this->linear_solver, problem.number_variables, problem.number_constraints,
dual_regularization_parameter, warmstart_information);

// check the inertia
[[maybe_unused]] auto [number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues] = this->linear_solver->get_inertia();
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,8 @@ namespace uno {
const Vector<double>& primal_direction, double tau);
[[nodiscard]] static double dual_fraction_to_boundary(const OptimizationProblem& problem, const Multipliers& current_multipliers,
Multipliers& direction_multipliers, double tau);
void assemble_augmented_system(Statistics& statistics, const OptimizationProblem& problem, const Multipliers& current_multipliers);
void assemble_augmented_system(Statistics& statistics, const OptimizationProblem& problem, const Multipliers& current_multipliers,
const WarmstartInformation& warmstart_information);
void assemble_augmented_rhs(const OptimizationProblem& problem, const Multipliers& current_multipliers);
void assemble_primal_dual_direction(const OptimizationProblem& problem, const Vector<double>& current_primals, const Multipliers& current_multipliers,
Vector<double>& direction_primals, Multipliers& direction_multipliers);
Expand Down
55 changes: 29 additions & 26 deletions uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,9 @@
#include "RectangularMatrix.hpp"
#include "ingredients/hessian_models/UnstableRegularization.hpp"
#include "model/Model.hpp"
#include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
#include "optimization/WarmstartInformation.hpp"
#include "options/Options.hpp"
#include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
#include "tools/Statistics.hpp"

namespace uno {
Expand All @@ -26,14 +27,13 @@ namespace uno {
const Options& options);
void assemble_matrix(const SymmetricMatrix<size_t, double>& hessian, const RectangularMatrix<double>& constraint_jacobian,
size_t number_variables, size_t number_constraints);
void factorize_matrix(const Model& model, DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver);
void regularize_matrix(Statistics& statistics, const Model& model, DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver,
size_t size_primal_block, size_t size_dual_block, ElementType dual_regularization_parameter);
void factorize_matrix(DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver, const WarmstartInformation& warmstart_information);
void regularize_matrix(Statistics& statistics, DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver,
size_t size_primal_block, size_t size_dual_block, ElementType dual_regularization_parameter, const WarmstartInformation& warmstart_information);
void solve(DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver);
// [[nodiscard]] T get_primal_regularization() const;

protected:
size_t number_factorizations{0};
ElementType primal_regularization{0.};
ElementType dual_regularization{0.};
ElementType previous_primal_regularization{0.};
Expand Down Expand Up @@ -89,34 +89,36 @@ namespace uno {
}

template <typename ElementType>
void SymmetricIndefiniteLinearSystem<ElementType>::factorize_matrix(const Model& /*model*/,
DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver) {
if (true || this->number_factorizations == 0) {
void SymmetricIndefiniteLinearSystem<ElementType>::factorize_matrix(DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver,
const WarmstartInformation& warmstart_information) {
if (warmstart_information.problem_structure_changed) {
DEBUG << "Performing symbolic analysis of the indefinite system\n";
linear_solver.do_symbolic_analysis(this->matrix);
}
DEBUG << "Performing numerical factorization of the indefinite system\n";
linear_solver.do_numerical_factorization(this->matrix);
this->number_factorizations++;
}

// the matrix has been factorized prior to calling this function
template <typename ElementType>
void SymmetricIndefiniteLinearSystem<ElementType>::regularize_matrix(Statistics& statistics, const Model& model,
void SymmetricIndefiniteLinearSystem<ElementType>::regularize_matrix(Statistics& statistics,
DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver, size_t size_primal_block, size_t size_dual_block,
ElementType dual_regularization_parameter) {
ElementType dual_regularization_parameter, const WarmstartInformation& warmstart_information) {
DEBUG2 << "Original matrix\n" << this->matrix << '\n';
this->primal_regularization = ElementType(0.);
this->dual_regularization = ElementType(0.);
size_t number_attempts = 1;
DEBUG << "Testing factorization with regularization factors (" << this->primal_regularization << ", " << this->dual_regularization << ")\n";
DEBUG << "Number of attempts: " << number_attempts << "\n\n";

auto [number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues] = linear_solver.get_inertia();
DEBUG << "Expected inertia (" << size_primal_block << ", " << size_dual_block << ", 0)\n";
DEBUG << "Estimated inertia (" << number_pos_eigenvalues << ", " << number_neg_eigenvalues << ", " << number_zero_eigenvalues << ")\n";

if (not linear_solver.matrix_is_singular() && linear_solver.number_negative_eigenvalues() == size_dual_block) {
DEBUG << "Inertia is good\n";
if (number_pos_eigenvalues == size_primal_block && number_neg_eigenvalues == size_dual_block && number_zero_eigenvalues == 0) {
DEBUG << "The inertia is correct\n";
statistics.set("regulariz", this->primal_regularization);
return;
}
auto [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\n";

// set the constraint regularization coefficient
if (linear_solver.matrix_is_singular()) {
Expand All @@ -141,19 +143,20 @@ namespace uno {
while (not good_inertia) {
DEBUG << "Testing factorization with regularization factors (" << this->primal_regularization << ", " << this->dual_regularization << ")\n";
DEBUG2 << this->matrix << '\n';
this->factorize_matrix(model, linear_solver);
this->factorize_matrix(linear_solver, warmstart_information);
number_attempts++;
DEBUG << "Number of attempts: " << number_attempts << "\n";

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)\n";
DEBUG << "Estimated inertia (" << number_pos_eigenvalues << ", " << number_neg_eigenvalues << ", " << number_zero_eigenvalues << ")\n";

if (not linear_solver.matrix_is_singular() && linear_solver.number_negative_eigenvalues() == size_dual_block) {
if (number_pos_eigenvalues == size_primal_block && number_neg_eigenvalues == size_dual_block && number_zero_eigenvalues == 0) {
good_inertia = true;
DEBUG << "Factorization was a success\n";
DEBUG << "The inertia is correct\n";
this->previous_primal_regularization = this->primal_regularization;
}
else {
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";
if (this->previous_primal_regularization == 0. || this->threshold_unsuccessful_attempts < number_attempts) {
this->primal_regularization *= this->primal_regularization_fast_increase_factor;
}
Expand Down Expand Up @@ -188,4 +191,4 @@ namespace uno {
*/
} // namespace

#endif // UNO_SYMMETRICINDEFINITELINEARSYSTEM_H
#endif // UNO_SYMMETRICINDEFINITELINEARSYSTEM_H
8 changes: 4 additions & 4 deletions uno/optimization/WarmstartInformation.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -10,15 +10,15 @@ namespace uno {
std::cout << "Constraints: " << std::boolalpha << this->constraints_changed << '\n';
std::cout << "Constraint bounds: " << std::boolalpha << this->constraint_bounds_changed << '\n';
std::cout << "Variable bounds: " << std::boolalpha << this->variable_bounds_changed << '\n';
std::cout << "Problem: " << std::boolalpha << this->problem_changed << '\n';
std::cout << "Problem structure: " << std::boolalpha << this->problem_structure_changed << '\n';
}

void WarmstartInformation::no_changes() {
this->objective_changed = false;
this->constraints_changed = false;
this->constraint_bounds_changed = false;
this->variable_bounds_changed = false;
this->problem_changed = false;
this->problem_structure_changed = false;
}

void WarmstartInformation::iterate_changed() {
Expand All @@ -33,14 +33,14 @@ namespace uno {
this->constraints_changed = true;
this->constraint_bounds_changed = true;
this->variable_bounds_changed = true;
this->problem_changed = true;
this->problem_structure_changed = true;
}

void WarmstartInformation::only_objective_changed() {
this->objective_changed = true;
this->constraints_changed = false;
this->constraint_bounds_changed = false;
this->variable_bounds_changed = false;
this->problem_changed = false;
this->problem_structure_changed = false;
}
} // namespace
2 changes: 1 addition & 1 deletion uno/optimization/WarmstartInformation.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ namespace uno {
bool constraints_changed{true};
bool constraint_bounds_changed{true};
bool variable_bounds_changed{true};
bool problem_changed{true};
bool problem_structure_changed{true};

void display() const;
void no_changes();
Expand Down
4 changes: 2 additions & 2 deletions uno/solvers/BQPD/BQPDSolver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -165,8 +165,8 @@ namespace uno {

BQPDMode BQPDSolver::determine_mode(const WarmstartInformation& warmstart_information) const {
BQPDMode mode = (this->number_calls == 0) ? BQPDMode::ACTIVE_SET_EQUALITIES : BQPDMode::USER_DEFINED;
// if problem changed, use cold start
if (warmstart_information.problem_changed) {
// if problem structure changed, use cold start
if (warmstart_information.problem_structure_changed) {
mode = BQPDMode::ACTIVE_SET_EQUALITIES;
}
// if only the variable bounds changed, reuse the active set estimate and the Jacobian information
Expand Down
2 changes: 1 addition & 1 deletion uno/solvers/MA57/MA57Solver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -165,7 +165,7 @@ namespace uno {
// build the internal matrix representation
this->row_indices.clear();
this->column_indices.clear();
for (const auto [row_index, column_index, element]: matrix) {
for (const auto [row_index, column_index, _]: matrix) {
this->row_indices.emplace_back(static_cast<int>(row_index + this->fortran_shift));
this->column_indices.emplace_back(static_cast<int>(column_index + this->fortran_shift));
}
Expand Down
Loading
Loading