54#include "Eigen/SparseCore"
55#include "absl/algorithm/container.h"
56#include "absl/base/nullability.h"
57#include "absl/log/check.h"
58#include "absl/log/log.h"
59#include "absl/status/status.h"
60#include "absl/status/statusor.h"
61#include "absl/strings/str_cat.h"
62#include "absl/strings/str_format.h"
63#include "absl/strings/string_view.h"
64#include "absl/time/clock.h"
65#include "absl/time/time.h"
66#include "google/protobuf/repeated_ptr_field.h"
91using ::Eigen::VectorXd;
92using ::operations_research::SolverLogger;
94using IterationStatsCallback =
99int NumThreads(
const int num_threads,
const int num_shards,
101 int capped_num_threads = num_threads;
102 if (num_shards > 0) {
103 capped_num_threads = std::min(capped_num_threads, num_shards);
105 const int64_t problem_limit = std::max(qp.variable_lower_bounds.size(),
106 qp.constraint_lower_bounds.size());
108 static_cast<int>(std::min(int64_t{capped_num_threads}, problem_limit));
109 capped_num_threads = std::max(capped_num_threads, 1);
110 if (capped_num_threads != num_threads) {
111 SOLVER_LOG(&logger,
"WARNING: Reducing num_threads from ", num_threads,
112 " to ", capped_num_threads,
113 " because additional threads would be useless.");
115 return capped_num_threads;
121int NumShards(
const int num_threads,
const int num_shards) {
122 if (num_shards > 0)
return num_shards;
123 return num_threads == 1 ? 1 : 4 * num_threads;
126std::string ConvergenceInformationString(
130 constexpr absl::string_view kFormatStr =
131 "%#12.6g %#12.6g %#12.6g | %#12.6g %#12.6g %#12.6g | %#12.6g %#12.6g | "
133 switch (residual_norm) {
135 return absl::StrFormat(
136 kFormatStr, relative_information.relative_l_inf_primal_residual,
137 relative_information.relative_l_inf_dual_residual,
138 relative_information.relative_optimality_gap,
139 convergence_information.l_inf_primal_residual(),
140 convergence_information.l_inf_dual_residual(),
141 convergence_information.primal_objective() -
142 convergence_information.dual_objective(),
143 convergence_information.primal_objective(),
144 convergence_information.dual_objective(),
145 convergence_information.l2_primal_variable(),
146 convergence_information.l2_dual_variable());
148 return absl::StrFormat(kFormatStr,
149 relative_information.relative_l2_primal_residual,
150 relative_information.relative_l2_dual_residual,
151 relative_information.relative_optimality_gap,
152 convergence_information.l2_primal_residual(),
153 convergence_information.l2_dual_residual(),
154 convergence_information.primal_objective() -
155 convergence_information.dual_objective(),
156 convergence_information.primal_objective(),
157 convergence_information.dual_objective(),
158 convergence_information.l2_primal_variable(),
159 convergence_information.l2_dual_variable());
161 return absl::StrFormat(
163 convergence_information.l_inf_componentwise_primal_residual(),
164 convergence_information.l_inf_componentwise_dual_residual(),
165 relative_information.relative_optimality_gap,
166 convergence_information.l_inf_primal_residual(),
167 convergence_information.l_inf_dual_residual(),
168 convergence_information.primal_objective() -
169 convergence_information.dual_objective(),
170 convergence_information.primal_objective(),
171 convergence_information.dual_objective(),
172 convergence_information.l2_primal_variable(),
173 convergence_information.l2_dual_variable());
175 LOG(FATAL) <<
"Residual norm not specified.";
177 LOG(FATAL) <<
"Invalid residual norm " << residual_norm <<
".";
180std::string ConvergenceInformationShortString(
184 constexpr absl::string_view kFormatStr =
185 "%#10.4g %#10.4g %#10.4g | %#10.4g %#10.4g";
186 switch (residual_norm) {
188 return absl::StrFormat(
189 kFormatStr, relative_information.relative_l_inf_primal_residual,
190 relative_information.relative_l_inf_dual_residual,
191 relative_information.relative_optimality_gap,
192 convergence_information.primal_objective(),
193 convergence_information.dual_objective());
195 return absl::StrFormat(kFormatStr,
196 relative_information.relative_l2_primal_residual,
197 relative_information.relative_l2_dual_residual,
198 relative_information.relative_optimality_gap,
199 convergence_information.primal_objective(),
200 convergence_information.dual_objective());
202 return absl::StrFormat(
204 convergence_information.l_inf_componentwise_primal_residual(),
205 convergence_information.l_inf_componentwise_dual_residual(),
206 relative_information.relative_optimality_gap,
207 convergence_information.primal_objective(),
208 convergence_information.dual_objective());
210 LOG(FATAL) <<
"Residual norm not specified.";
212 LOG(FATAL) <<
"Invalid residual norm " << residual_norm <<
".";
220void LogIterationStats(
int verbosity_level,
bool use_feasibility_polishing,
226 std::string iteration_string =
228 ? absl::StrFormat(
"%6d %8.1f %6.1f", iter_stats.iteration_number(),
229 iter_stats.cumulative_kkt_matrix_passes(),
230 iter_stats.cumulative_time_sec())
231 : absl::StrFormat(
"%6d %6.1f", iter_stats.iteration_number(),
232 iter_stats.cumulative_time_sec());
233 auto convergence_information =
235 if (!convergence_information.has_value() &&
236 iter_stats.convergence_information_size() > 0) {
237 convergence_information = iter_stats.convergence_information(0);
239 const char* phase_string = [&]() {
240 if (use_feasibility_polishing) {
241 switch (iteration_type) {
257 if (convergence_information.has_value()) {
258 const char* iterate_string = [&]() {
259 if (verbosity_level >= 4) {
260 switch (convergence_information->candidate_type()) {
279 *convergence_information, bound_norms);
280 std::string convergence_string =
282 ? ConvergenceInformationString(
283 *convergence_information, relative_information,
284 termination_criteria.optimality_norm())
285 : ConvergenceInformationShortString(
286 *convergence_information, relative_information,
287 termination_criteria.optimality_norm());
288 SOLVER_LOG(&logger, phase_string, iterate_string, iteration_string,
" | ",
292 SOLVER_LOG(&logger, phase_string, verbosity_level >= 4 ?
"? " :
"",
297void LogIterationStatsHeader(
int verbosity_level,
298 bool use_feasibility_polishing,
300 std::string work_labels =
302 ? absl::StrFormat(
"%6s %8s %6s",
"iter#",
"kkt_pass",
"time")
303 : absl::StrFormat(
"%6s %6s",
"iter#",
"time");
304 std::string convergence_labels =
307 "%12s %12s %12s | %12s %12s %12s | %12s %12s | %12s %12s",
308 "rel_prim_res",
"rel_dual_res",
"rel_gap",
"prim_resid",
309 "dual_resid",
"obj_gap",
"prim_obj",
"dual_obj",
"prim_var_l2",
311 : absl::StrFormat(
"%10s %10s %10s | %10s %10s",
"rel_p_res",
312 "rel_d_res",
"rel_gap",
"prim_obj",
"dual_obj");
313 SOLVER_LOG(&logger, use_feasibility_polishing ?
"f " :
"",
314 verbosity_level >= 4 ?
"I " :
"", work_labels,
" | ",
318enum class InnerStepOutcome {
320 kForceNumericalTermination,
328 VectorXd dual_solution,
332 solve_log.set_iteration_count(stats.iteration_number());
333 solve_log.set_termination_reason(termination_reason);
334 solve_log.set_solution_type(output_type);
335 solve_log.set_solve_time_sec(stats.cumulative_time_sec());
336 *solve_log.mutable_solution_stats() = stats;
338 .dual_solution = std::move(dual_solution),
339 .solve_log = std::move(solve_log)};
343class PreprocessSolver {
352 explicit PreprocessSolver(QuadraticProgram qp,
353 const PrimalDualHybridGradientParams& params,
354 SolverLogger* logger);
357 PreprocessSolver(
const PreprocessSolver&) =
delete;
358 PreprocessSolver& operator=(
const PreprocessSolver&) =
delete;
359 PreprocessSolver(PreprocessSolver&&) =
delete;
360 PreprocessSolver& operator=(PreprocessSolver&&) =
delete;
369 SolverResult PreprocessAndSolve(
370 const PrimalDualHybridGradientParams& params,
371 std::optional<PrimalAndDualSolution> initial_solution,
372 const std::atomic<bool>* interrupt_solve,
373 IterationStatsCallback iteration_stats_callback);
384 std::optional<TerminationReasonAndPointType>
385 UpdateIterationStatsAndCheckTermination(
386 const PrimalDualHybridGradientParams& params,
387 bool force_numerical_termination,
const VectorXd& working_primal_current,
388 const VectorXd& working_dual_current,
389 const VectorXd* working_primal_average,
390 const VectorXd* working_dual_average,
391 const VectorXd* working_primal_delta,
const VectorXd* working_dual_delta,
392 const VectorXd& last_primal_start_point,
393 const VectorXd& last_dual_start_point,
394 const std::atomic<bool>* interrupt_solve,
IterationType iteration_type,
395 const IterationStats& full_stats, IterationStats& stats);
402 void ComputeConvergenceAndInfeasibilityFromWorkingSolution(
403 const PrimalDualHybridGradientParams& params,
404 const VectorXd& working_primal,
const VectorXd& working_dual,
405 PointType candidate_type, ConvergenceInformation* convergence_information,
406 InfeasibilityInformation* infeasibility_information)
const;
412 SolverResult ConstructOriginalSolverResult(
413 const PrimalDualHybridGradientParams& params, SolverResult result,
414 SolverLogger& logger)
const;
416 const ShardedQuadraticProgram& ShardedWorkingQp()
const {
422 void SwapVariableBounds(VectorXd& variable_lower_bounds,
423 VectorXd& variable_upper_bounds) {
424 sharded_qp_.SwapVariableBounds(variable_lower_bounds,
425 variable_upper_bounds);
430 void SwapConstraintBounds(VectorXd& constraint_lower_bounds,
431 VectorXd& constraint_upper_bounds) {
432 sharded_qp_.SwapConstraintBounds(constraint_lower_bounds,
433 constraint_upper_bounds);
439 void SwapObjectiveVector(VectorXd& objective) {
440 sharded_qp_.SwapObjectiveVector(objective);
443 const QuadraticProgramBoundNorms& OriginalBoundNorms()
const {
444 return original_bound_norms_;
447 SolverLogger& Logger() {
return logger_; }
450 struct PresolveInfo {
451 explicit PresolveInfo(ShardedQuadraticProgram original_qp,
452 const PrimalDualHybridGradientParams& params)
453 : preprocessor_parameters(PreprocessorParameters(params)),
454 preprocessor(&preprocessor_parameters),
455 sharded_original_qp(std::move(original_qp)),
456 trivial_col_scaling_vec(
457 OnesVector(sharded_original_qp.PrimalSharder())),
458 trivial_row_scaling_vec(
459 OnesVector(sharded_original_qp.DualSharder())) {}
461 glop::GlopParameters preprocessor_parameters;
462 glop::MainLpPreprocessor preprocessor;
463 ShardedQuadraticProgram sharded_original_qp;
464 bool presolved_problem_was_maximization =
false;
465 const VectorXd trivial_col_scaling_vec, trivial_row_scaling_vec;
469 static glop::GlopParameters PreprocessorParameters(
470 const PrimalDualHybridGradientParams& params);
478 std::optional<TerminationReason> ApplyPresolveIfEnabled(
479 const PrimalDualHybridGradientParams& params,
480 std::optional<PrimalAndDualSolution>* initial_solution);
482 void ComputeAndApplyRescaling(
const PrimalDualHybridGradientParams& params,
483 VectorXd& starting_primal_solution,
484 VectorXd& starting_dual_solution);
486 void LogQuadraticProgramStats(
const QuadraticProgramStats& stats)
const;
488 double InitialPrimalWeight(
const PrimalDualHybridGradientParams& params,
489 double l2_norm_primal_linear_objective,
490 double l2_norm_constraint_bounds)
const;
492 PrimalAndDualSolution RecoverOriginalSolution(
493 PrimalAndDualSolution working_solution)
const;
496 void AddPointMetadata(
const PrimalDualHybridGradientParams& params,
498 const VectorXd& dual_solution,
PointType point_type,
499 const VectorXd& last_primal_start_point,
500 const VectorXd& last_dual_start_point,
501 IterationStats& stats)
const;
503 const QuadraticProgram& Qp()
const {
return sharded_qp_.Qp(); }
505 const int num_threads_;
506 const int num_shards_;
509 QuadraticProgramBoundNorms original_bound_norms_;
517 ShardedQuadraticProgram sharded_qp_;
520 std::optional<PresolveInfo> presolve_info_;
525 VectorXd col_scaling_vec_;
526 VectorXd row_scaling_vec_;
529 int log_counter_ = 0;
530 absl::Time time_of_last_log_ = absl::InfinitePast();
531 SolverLogger& logger_;
532 IterationStatsCallback iteration_stats_callback_;
540 explicit Solver(
const PrimalDualHybridGradientParams& params,
541 VectorXd starting_primal_solution,
542 VectorXd starting_dual_solution,
double initial_step_size,
543 double initial_primal_weight,
544 PreprocessSolver* preprocess_solver);
547 Solver(
const Solver&) =
delete;
548 Solver& operator=(
const Solver&) =
delete;
549 Solver(Solver&&) =
delete;
550 Solver& operator=(Solver&&) =
delete;
552 const QuadraticProgram& WorkingQp()
const {
return ShardedWorkingQp().Qp(); }
554 const ShardedQuadraticProgram& ShardedWorkingQp()
const {
555 return preprocess_solver_->ShardedWorkingQp();
567 const std::atomic<bool>* interrupt_solve,
571 struct NextSolutionAndDelta {
577 struct DistanceBasedRestartInfo {
578 double distance_moved_last_restart_period;
579 int length_of_last_restart_period;
585 constexpr static double kDivergentMovement = 1.0e100;
594 constexpr static int kFeasibilityIterationFraction = 8;
601 std::optional<SolverResult> TryFeasibilityPolishing(
602 int iteration_limit,
const std::atomic<bool>* interrupt_solve,
603 SolveLog& solve_log);
607 SolverResult TryPrimalPolishing(VectorXd starting_primal_solution,
609 const std::atomic<bool>* interrupt_solve,
610 SolveLog& solve_log);
614 SolverResult TryDualPolishing(VectorXd starting_dual_solution,
616 const std::atomic<bool>* interrupt_solve,
617 SolveLog& solve_log);
619 NextSolutionAndDelta ComputeNextPrimalSolution(
double primal_step_size)
const;
621 NextSolutionAndDelta ComputeNextDualSolution(
622 double dual_step_size,
double extrapolation_factor,
623 const NextSolutionAndDelta& next_primal_solution,
624 const VectorXd* absl_nullable next_primal_product =
nullptr)
const;
626 std::pair<double, double> ComputeMovementTerms(
627 const VectorXd& delta_primal,
const VectorXd& delta_dual)
const;
629 double ComputeMovement(
const VectorXd& delta_primal,
630 const VectorXd& delta_dual)
const;
632 double ComputeNonlinearity(
const VectorXd& delta_primal,
633 const VectorXd& next_dual_product)
const;
637 void SetCurrentPrimalAndDualProducts();
640 IterationStats CreateSimpleIterationStats(
RestartChoice restart_used)
const;
644 IterationStats TotalWorkSoFar(
const SolveLog& solve_log)
const;
648 VectorXd PrimalAverage()
const;
650 VectorXd DualAverage()
const;
652 double ComputeNewPrimalWeight()
const;
664 SolverResult PickSolutionAndConstructSolverResult(
667 PointType output_type, SolveLog solve_log)
const;
670 const VectorXd& dual_solution)
const;
672 LocalizedLagrangianBounds ComputeLocalizedBoundsAtCurrent()
const;
674 LocalizedLagrangianBounds ComputeLocalizedBoundsAtAverage()
const;
680 std::optional<SolverResult> MajorIterationAndTerminationCheck(
681 IterationType iteration_type,
bool force_numerical_termination,
682 const std::atomic<bool>* interrupt_solve,
683 const IterationStats& work_from_feasibility_polishing,
684 SolveLog& solve_log);
686 bool ShouldDoAdaptiveRestartHeuristic(
double candidate_normalized_gap)
const;
690 void ResetAverageToCurrent();
692 void LogNumericalTermination(
const Eigen::VectorXd& primal_delta,
693 const Eigen::VectorXd& dual_delta)
const;
695 void LogInnerIterationLimitHit()
const;
704 InnerStepOutcome TakeMalitskyPockStep();
708 InnerStepOutcome TakeAdaptiveStep();
711 InnerStepOutcome TakeConstantSizeStep();
713 const PrimalDualHybridGradientParams params_;
715 VectorXd current_primal_solution_;
716 VectorXd current_dual_solution_;
717 VectorXd current_primal_delta_;
718 VectorXd current_dual_delta_;
720 ShardedWeightedAverage primal_average_;
721 ShardedWeightedAverage dual_average_;
724 double primal_weight_;
726 PreprocessSolver* preprocess_solver_;
729 double ratio_last_two_step_sizes_;
731 double normalized_gap_at_last_trial_ =
732 std::numeric_limits<double>::infinity();
734 double normalized_gap_at_last_restart_ =
735 std::numeric_limits<double>::infinity();
739 double preprocessing_time_sec_;
741 int iterations_completed_;
742 int num_rejected_steps_;
745 std::optional<VectorXd> current_primal_product_;
747 VectorXd current_dual_product_;
750 VectorXd last_primal_start_point_;
753 VectorXd last_dual_start_point_;
757 DistanceBasedRestartInfo distance_based_restart_info_ = {
758 .distance_moved_last_restart_period =
759 std::numeric_limits<double>::infinity(),
760 .length_of_last_restart_period = 1,
768 NumThreads(params.num_threads(), params.num_shards(), qp, *logger)),
769 num_shards_(NumShards(num_threads_, params.num_shards())),
770 sharded_qp_(std::move(qp), num_threads_, num_shards_,
771 params.scheduler_type(), nullptr),
774SolverResult ErrorSolverResult(
const TerminationReason reason,
775 const std::string& message,
778 error_log.set_termination_reason(reason);
779 error_log.set_termination_string(message);
781 "The solver did not run because of invalid input: ", message);
782 return SolverResult{.solve_log = error_log};
789std::optional<SolverResult> CheckProblemStats(
790 const QuadraticProgramStats& problem_stats,
const double objective_offset,
791 bool check_excessively_small_values,
SolverLogger& logger) {
792 const double kExcessiveInputValue = 1e50;
793 const double kExcessivelySmallInputValue = 1e-50;
794 const double kMaxDynamicRange = 1e20;
795 if (std::isnan(problem_stats.constraint_matrix_l2_norm())) {
796 return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
797 "Constraint matrix has a NAN.", logger);
799 if (problem_stats.constraint_matrix_abs_max() > kExcessiveInputValue) {
800 return ErrorSolverResult(
801 TERMINATION_REASON_INVALID_PROBLEM,
802 absl::StrCat(
"Constraint matrix has a non-zero with absolute value ",
803 problem_stats.constraint_matrix_abs_max(),
804 " which exceeds limit of ", kExcessiveInputValue,
"."),
807 if (problem_stats.constraint_matrix_abs_max() >
808 kMaxDynamicRange * problem_stats.constraint_matrix_abs_min()) {
810 &logger,
"WARNING: Constraint matrix has largest absolute value ",
811 problem_stats.constraint_matrix_abs_max(),
812 " and smallest non-zero absolute value ",
813 problem_stats.constraint_matrix_abs_min(),
" performance may suffer.");
815 if (problem_stats.constraint_matrix_col_min_l_inf_norm() > 0 &&
816 problem_stats.constraint_matrix_col_min_l_inf_norm() <
817 kExcessivelySmallInputValue) {
818 return ErrorSolverResult(
819 TERMINATION_REASON_INVALID_PROBLEM,
820 absl::StrCat(
"Constraint matrix has a column with Linf norm ",
821 problem_stats.constraint_matrix_col_min_l_inf_norm(),
822 " which is less than limit of ",
823 kExcessivelySmallInputValue,
"."),
826 if (problem_stats.constraint_matrix_row_min_l_inf_norm() > 0 &&
827 problem_stats.constraint_matrix_row_min_l_inf_norm() <
828 kExcessivelySmallInputValue) {
829 return ErrorSolverResult(
830 TERMINATION_REASON_INVALID_PROBLEM,
831 absl::StrCat(
"Constraint matrix has a row with Linf norm ",
832 problem_stats.constraint_matrix_row_min_l_inf_norm(),
833 " which is less than limit of ",
834 kExcessivelySmallInputValue,
"."),
837 if (std::isnan(problem_stats.combined_bounds_l2_norm())) {
838 return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
839 "Constraint bounds vector has a NAN.", logger);
841 if (problem_stats.combined_bounds_max() > kExcessiveInputValue) {
842 return ErrorSolverResult(
843 TERMINATION_REASON_INVALID_PROBLEM,
844 absl::StrCat(
"Combined constraint bounds vector has a non-zero with "
846 problem_stats.combined_bounds_max(),
847 " which exceeds limit of ", kExcessiveInputValue,
"."),
850 if (check_excessively_small_values &&
851 problem_stats.combined_bounds_min() > 0 &&
852 problem_stats.combined_bounds_min() < kExcessivelySmallInputValue) {
853 return ErrorSolverResult(
854 TERMINATION_REASON_INVALID_PROBLEM,
855 absl::StrCat(
"Combined constraint bounds vector has a non-zero with "
857 problem_stats.combined_bounds_min(),
858 " which is less than the limit of ",
859 kExcessivelySmallInputValue,
"."),
862 if (problem_stats.combined_bounds_max() >
863 kMaxDynamicRange * problem_stats.combined_bounds_min()) {
865 "WARNING: Combined constraint bounds vector has largest "
867 problem_stats.combined_bounds_max(),
868 " and smallest non-zero absolute value ",
869 problem_stats.combined_bounds_min(),
870 "; performance may suffer.");
872 if (std::isnan(problem_stats.combined_variable_bounds_l2_norm())) {
873 return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
874 "Variable bounds vector has a NAN.", logger);
876 if (problem_stats.combined_variable_bounds_max() > kExcessiveInputValue) {
877 return ErrorSolverResult(
878 TERMINATION_REASON_INVALID_PROBLEM,
879 absl::StrCat(
"Combined variable bounds vector has a non-zero with "
881 problem_stats.combined_variable_bounds_max(),
882 " which exceeds limit of ", kExcessiveInputValue,
"."),
885 if (check_excessively_small_values &&
886 problem_stats.combined_variable_bounds_min() > 0 &&
887 problem_stats.combined_variable_bounds_min() <
888 kExcessivelySmallInputValue) {
889 return ErrorSolverResult(
890 TERMINATION_REASON_INVALID_PROBLEM,
891 absl::StrCat(
"Combined variable bounds vector has a non-zero with "
893 problem_stats.combined_variable_bounds_min(),
894 " which is less than the limit of ",
895 kExcessivelySmallInputValue,
"."),
898 if (problem_stats.combined_variable_bounds_max() >
899 kMaxDynamicRange * problem_stats.combined_variable_bounds_min()) {
902 "WARNING: Combined variable bounds vector has largest absolute value ",
903 problem_stats.combined_variable_bounds_max(),
904 " and smallest non-zero absolute value ",
905 problem_stats.combined_variable_bounds_min(),
906 "; performance may suffer.");
908 if (problem_stats.variable_bound_gaps_max() >
909 kMaxDynamicRange * problem_stats.variable_bound_gaps_min()) {
911 "WARNING: Variable bound gap vector has largest absolute value ",
912 problem_stats.variable_bound_gaps_max(),
913 " and smallest non-zero absolute value ",
914 problem_stats.variable_bound_gaps_min(),
915 "; performance may suffer.");
917 if (std::isnan(objective_offset)) {
918 return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
919 "Objective offset is NAN.", logger);
921 if (std::abs(objective_offset) > kExcessiveInputValue) {
922 return ErrorSolverResult(
923 TERMINATION_REASON_INVALID_PROBLEM,
924 absl::StrCat(
"Objective offset ", objective_offset,
925 " has absolute value which exceeds limit of ",
926 kExcessiveInputValue,
"."),
929 if (std::isnan(problem_stats.objective_vector_l2_norm())) {
930 return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
931 "Objective vector has a NAN.", logger);
933 if (problem_stats.objective_vector_abs_max() > kExcessiveInputValue) {
934 return ErrorSolverResult(
935 TERMINATION_REASON_INVALID_PROBLEM,
936 absl::StrCat(
"Objective vector has a non-zero with absolute value ",
937 problem_stats.objective_vector_abs_max(),
938 " which exceeds limit of ", kExcessiveInputValue,
"."),
941 if (check_excessively_small_values &&
942 problem_stats.objective_vector_abs_min() > 0 &&
943 problem_stats.objective_vector_abs_min() < kExcessivelySmallInputValue) {
944 return ErrorSolverResult(
945 TERMINATION_REASON_INVALID_PROBLEM,
946 absl::StrCat(
"Objective vector has a non-zero with absolute value ",
947 problem_stats.objective_vector_abs_min(),
948 " which is less than the limit of ",
949 kExcessivelySmallInputValue,
"."),
952 if (problem_stats.objective_vector_abs_max() >
953 kMaxDynamicRange * problem_stats.objective_vector_abs_min()) {
954 SOLVER_LOG(&logger,
"WARNING: Objective vector has largest absolute value ",
955 problem_stats.objective_vector_abs_max(),
956 " and smallest non-zero absolute value ",
957 problem_stats.objective_vector_abs_min(),
958 "; performance may suffer.");
960 if (std::isnan(problem_stats.objective_matrix_l2_norm())) {
961 return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
962 "Objective matrix has a NAN.", logger);
964 if (problem_stats.objective_matrix_abs_max() > kExcessiveInputValue) {
965 return ErrorSolverResult(
966 TERMINATION_REASON_INVALID_PROBLEM,
967 absl::StrCat(
"Objective matrix has a non-zero with absolute value ",
968 problem_stats.objective_matrix_abs_max(),
969 " which exceeds limit of ", kExcessiveInputValue,
"."),
972 if (problem_stats.objective_matrix_abs_max() >
973 kMaxDynamicRange * problem_stats.objective_matrix_abs_min()) {
974 SOLVER_LOG(&logger,
"WARNING: Objective matrix has largest absolute value ",
975 problem_stats.objective_matrix_abs_max(),
976 " and smallest non-zero absolute value ",
977 problem_stats.objective_matrix_abs_min(),
978 "; performance may suffer.");
983std::optional<SolverResult> CheckInitialSolution(
984 const ShardedQuadraticProgram& sharded_qp,
985 const PrimalAndDualSolution& initial_solution,
SolverLogger& logger) {
986 const double kExcessiveInputValue = 1e50;
987 if (initial_solution.primal_solution.size() != sharded_qp.PrimalSize()) {
988 return ErrorSolverResult(
989 TERMINATION_REASON_INVALID_INITIAL_SOLUTION,
990 absl::StrCat(
"Initial primal solution has size ",
991 initial_solution.primal_solution.size(),
992 " which differs from problem primal size ",
993 sharded_qp.PrimalSize()),
997 Norm(initial_solution.primal_solution, sharded_qp.PrimalSharder()))) {
998 return ErrorSolverResult(TERMINATION_REASON_INVALID_INITIAL_SOLUTION,
999 "Initial primal solution has a NAN.", logger);
1001 if (
const double norm =
LInfNorm(initial_solution.primal_solution,
1002 sharded_qp.PrimalSharder());
1003 norm > kExcessiveInputValue) {
1004 return ErrorSolverResult(
1005 TERMINATION_REASON_INVALID_INITIAL_SOLUTION,
1007 "Initial primal solution has an entry with absolute value ", norm,
1008 " which exceeds limit of ", kExcessiveInputValue),
1011 if (initial_solution.dual_solution.size() != sharded_qp.DualSize()) {
1012 return ErrorSolverResult(
1013 TERMINATION_REASON_INVALID_INITIAL_SOLUTION,
1014 absl::StrCat(
"Initial dual solution has size ",
1015 initial_solution.dual_solution.size(),
1016 " which differs from problem dual size ",
1017 sharded_qp.DualSize()),
1021 Norm(initial_solution.dual_solution, sharded_qp.DualSharder()))) {
1022 return ErrorSolverResult(TERMINATION_REASON_INVALID_INITIAL_SOLUTION,
1023 "Initial dual solution has a NAN.", logger);
1025 if (
const double norm =
1026 LInfNorm(initial_solution.dual_solution, sharded_qp.DualSharder());
1027 norm > kExcessiveInputValue) {
1028 return ErrorSolverResult(
1029 TERMINATION_REASON_INVALID_INITIAL_SOLUTION,
1030 absl::StrCat(
"Initial dual solution has an entry with absolute value ",
1031 norm,
" which exceeds limit of ", kExcessiveInputValue),
1034 return std::nullopt;
1037SolverResult PreprocessSolver::PreprocessAndSolve(
1038 const PrimalDualHybridGradientParams& params,
1039 std::optional<PrimalAndDualSolution> initial_solution,
1040 const std::atomic<bool>* interrupt_solve,
1041 IterationStatsCallback iteration_stats_callback) {
1045 if (params.verbosity_level() >= 1) {
1046 SOLVER_LOG(&logger_,
"Solving with PDLP parameters: ", params);
1048 if (Qp().problem_name.has_value()) {
1049 solve_log.set_instance_name(*Qp().problem_name);
1051 *solve_log.mutable_params() = params;
1052 sharded_qp_.ReplaceLargeConstraintBoundsWithInfinity(
1053 params.infinite_constraint_bound_threshold());
1055 return ErrorSolverResult(
1056 TERMINATION_REASON_INVALID_PROBLEM,
1057 "The input problem has invalid bounds (after replacing large "
1058 "constraint bounds with infinity): some variable or constraint has "
1059 "lower_bound > upper_bound, lower_bound == inf, or upper_bound == "
1063 if (Qp().objective_matrix.has_value() &&
1064 !sharded_qp_.PrimalSharder().ParallelTrueForAllShards(
1065 [&](
const Sharder::Shard& shard) ->
bool {
1066 return (shard(Qp().objective_matrix->diagonal()).array() >= 0.0)
1069 return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
1070 "The objective is not convex (i.e., the objective "
1071 "matrix contains negative or NAN entries).",
1074 *solve_log.mutable_original_problem_stats() =
ComputeStats(sharded_qp_);
1075 const QuadraticProgramStats& original_problem_stats =
1076 solve_log.original_problem_stats();
1077 if (
auto maybe_result =
1078 CheckProblemStats(original_problem_stats, Qp().objective_offset,
1079 params.presolve_options().use_glop(), logger_);
1080 maybe_result.has_value()) {
1081 return *maybe_result;
1083 if (initial_solution.has_value()) {
1084 if (
auto maybe_result =
1085 CheckInitialSolution(sharded_qp_, *initial_solution, logger_);
1086 maybe_result.has_value()) {
1087 return *maybe_result;
1091 const std::string preprocessing_string = absl::StrCat(
1092 params.presolve_options().use_glop() ?
"presolving and " :
"",
1094 if (params.verbosity_level() >= 1) {
1095 SOLVER_LOG(&logger_,
"Problem stats before ", preprocessing_string);
1096 LogQuadraticProgramStats(solve_log.original_problem_stats());
1098 iteration_stats_callback_ = std::move(iteration_stats_callback);
1099 std::optional<TerminationReason> maybe_terminate =
1100 ApplyPresolveIfEnabled(params, &initial_solution);
1101 if (maybe_terminate.has_value()) {
1107 IterationStats iteration_stats;
1108 iteration_stats.set_cumulative_time_sec(timer.
Get());
1109 solve_log.set_preprocessing_time_sec(iteration_stats.cumulative_time_sec());
1110 VectorXd working_primal =
ZeroVector(sharded_qp_.PrimalSharder());
1111 VectorXd working_dual =
ZeroVector(sharded_qp_.DualSharder());
1112 ComputeConvergenceAndInfeasibilityFromWorkingSolution(
1113 params, working_primal, working_dual, POINT_TYPE_PRESOLVER_SOLUTION,
1114 iteration_stats.add_convergence_information(),
1115 iteration_stats.add_infeasibility_information());
1116 std::optional<TerminationReasonAndPointType> earned_termination =
1118 iteration_stats, original_bound_norms_,
1120 if (!earned_termination.has_value()) {
1122 params.termination_criteria(), iteration_stats, interrupt_solve);
1125 if (earned_termination.has_value() &&
1126 (earned_termination->reason == TERMINATION_REASON_OPTIMAL ||
1127 earned_termination->reason == TERMINATION_REASON_PRIMAL_INFEASIBLE ||
1128 earned_termination->reason == TERMINATION_REASON_DUAL_INFEASIBLE)) {
1129 final_termination_reason = earned_termination->reason;
1131 if (*maybe_terminate == TERMINATION_REASON_OPTIMAL) {
1135 "WARNING: Presolve claimed to solve the LP optimally but the "
1136 "solution doesn't satisfy the optimality criteria.");
1138 final_termination_reason = *maybe_terminate;
1141 return ConstructOriginalSolverResult(
1143 ConstructSolverResult(
1144 std::move(working_primal), std::move(working_dual),
1145 std::move(iteration_stats), final_termination_reason,
1146 POINT_TYPE_PRESOLVER_SOLUTION, std::move(solve_log)),
1150 VectorXd starting_primal_solution;
1151 VectorXd starting_dual_solution;
1153 if (initial_solution.has_value()) {
1154 starting_primal_solution = std::move(initial_solution->primal_solution);
1155 starting_dual_solution = std::move(initial_solution->dual_solution);
1157 SetZero(sharded_qp_.PrimalSharder(), starting_primal_solution);
1158 SetZero(sharded_qp_.DualSharder(), starting_dual_solution);
1165 ComputeAndApplyRescaling(params, starting_primal_solution,
1166 starting_dual_solution);
1167 *solve_log.mutable_preprocessed_problem_stats() =
ComputeStats(sharded_qp_);
1168 if (params.verbosity_level() >= 1) {
1169 SOLVER_LOG(&logger_,
"Problem stats after ", preprocessing_string);
1170 LogQuadraticProgramStats(solve_log.preprocessed_problem_stats());
1173 double step_size = 0.0;
1174 if (params.linesearch_rule() ==
1175 PrimalDualHybridGradientParams::CONSTANT_STEP_SIZE_RULE) {
1176 std::mt19937 random(1);
1177 double inverse_step_size;
1178 const auto lipschitz_result =
1180 sharded_qp_, std::nullopt, std::nullopt,
1186 const double lipschitz_term_upper_bound =
1187 lipschitz_result.singular_value /
1188 (1.0 - lipschitz_result.estimated_relative_error);
1189 inverse_step_size = lipschitz_term_upper_bound;
1190 step_size = inverse_step_size > 0.0 ? 1.0 / inverse_step_size : 1.0;
1205 solve_log.preprocessed_problem_stats().constraint_matrix_abs_max());
1207 step_size *= params.initial_step_size_scaling();
1209 const double primal_weight = InitialPrimalWeight(
1210 params, solve_log.preprocessed_problem_stats().objective_vector_l2_norm(),
1211 solve_log.preprocessed_problem_stats().combined_bounds_l2_norm());
1213 Solver solver(params, starting_primal_solution, starting_dual_solution,
1214 step_size, primal_weight,
this);
1215 solve_log.set_preprocessing_time_sec(timer.
Get());
1216 SolverResult result = solver.Solve(IterationType::kNormal, interrupt_solve,
1217 std::move(solve_log));
1218 return ConstructOriginalSolverResult(params, std::move(result), logger_);
1221glop::GlopParameters PreprocessSolver::PreprocessorParameters(
1222 const PrimalDualHybridGradientParams& params) {
1223 glop::GlopParameters glop_params;
1225 glop_params.set_solve_dual_problem(glop::GlopParameters::NEVER_DO);
1228 glop_params.set_use_implied_free_preprocessor(
false);
1230 glop_params.set_use_scaling(
false);
1231 if (params.presolve_options().has_glop_parameters()) {
1232 glop_params.MergeFrom(params.presolve_options().glop_parameters());
1238 const glop::ProblemStatus glop_status,
SolverLogger& logger) {
1239 switch (glop_status) {
1240 case glop::ProblemStatus::OPTIMAL:
1242 case glop::ProblemStatus::INVALID_PROBLEM:
1244 case glop::ProblemStatus::ABNORMAL:
1245 case glop::ProblemStatus::IMPRECISE:
1247 case glop::ProblemStatus::PRIMAL_INFEASIBLE:
1248 case glop::ProblemStatus::DUAL_INFEASIBLE:
1249 case glop::ProblemStatus::INFEASIBLE_OR_UNBOUNDED:
1250 case glop::ProblemStatus::DUAL_UNBOUNDED:
1251 case glop::ProblemStatus::PRIMAL_UNBOUNDED:
1254 SOLVER_LOG(&logger,
"WARNING: Unexpected preprocessor status ",
1260std::optional<TerminationReason> PreprocessSolver::ApplyPresolveIfEnabled(
1261 const PrimalDualHybridGradientParams& params,
1262 std::optional<PrimalAndDualSolution>*
const initial_solution) {
1263 const bool presolve_enabled = params.presolve_options().use_glop();
1264 if (!presolve_enabled) {
1265 return std::nullopt;
1269 "WARNING: Skipping presolve, which is only supported for linear "
1271 return std::nullopt;
1276 "WARNING: Skipping presolve because of error converting to "
1279 return std::nullopt;
1281 if (initial_solution->has_value()) {
1283 "WARNING: Ignoring initial solution. Initial solutions are "
1284 "ignored when presolve is on.");
1285 initial_solution->reset();
1287 glop::LinearProgram glop_lp;
1288 glop::MPModelProtoToLinearProgram(*model, &glop_lp);
1291 presolve_info_.emplace(std::move(sharded_qp_), params);
1295 presolve_info_->preprocessor.Run(&glop_lp);
1296 presolve_info_->presolved_problem_was_maximization =
1297 glop_lp.IsMaximizationProblem();
1298 MPModelProto output;
1299 glop::LinearProgramToMPModelProto(glop_lp, &output);
1301 absl::StatusOr<QuadraticProgram> presolved_qp =
1303 CHECK_OK(presolved_qp.status());
1308 presolved_qp->objective_scaling_factor = glop_lp.objective_scaling_factor();
1309 sharded_qp_ = ShardedQuadraticProgram(std::move(*presolved_qp), num_threads_,
1310 num_shards_, params.scheduler_type());
1314 if (presolve_info_->preprocessor.status() != glop::ProblemStatus::INIT) {
1315 col_scaling_vec_ =
OnesVector(sharded_qp_.PrimalSharder());
1316 row_scaling_vec_ =
OnesVector(sharded_qp_.DualSharder());
1317 return GlopStatusToTerminationReason(presolve_info_->preprocessor.status(),
1320 return std::nullopt;
1323void PreprocessSolver::ComputeAndApplyRescaling(
1324 const PrimalDualHybridGradientParams& params,
1325 VectorXd& starting_primal_solution, VectorXd& starting_dual_solution) {
1327 RescalingOptions{.l_inf_ruiz_iterations = params.l_inf_ruiz_iterations(),
1328 .l2_norm_rescaling = params.l2_norm_rescaling()},
1330 row_scaling_vec_ = std::move(scaling.row_scaling_vec);
1331 col_scaling_vec_ = std::move(scaling.col_scaling_vec);
1334 starting_primal_solution);
1336 starting_dual_solution);
1339void PreprocessSolver::LogQuadraticProgramStats(
1340 const QuadraticProgramStats& stats)
const {
1342 absl::StrFormat(
"There are %i variables, %i constraints, and %i "
1343 "constraint matrix nonzeros.",
1344 stats.num_variables(), stats.num_constraints(),
1345 stats.constraint_matrix_num_nonzeros()));
1346 if (Qp().constraint_matrix.nonZeros() > 0) {
1348 absl::StrFormat(
"Absolute values of nonzero constraint matrix "
1349 "elements: largest=%f, "
1350 "smallest=%f, avg=%f",
1351 stats.constraint_matrix_abs_max(),
1352 stats.constraint_matrix_abs_min(),
1353 stats.constraint_matrix_abs_avg()));
1356 absl::StrFormat(
"Constraint matrix, infinity norm: max(row & col)=%f, "
1357 "min_col=%f, min_row=%f",
1358 stats.constraint_matrix_abs_max(),
1359 stats.constraint_matrix_col_min_l_inf_norm(),
1360 stats.constraint_matrix_row_min_l_inf_norm()));
1364 "Constraint bounds statistics (max absolute value per row): "
1365 "largest=%f, smallest=%f, avg=%f, l2_norm=%f",
1366 stats.combined_bounds_max(), stats.combined_bounds_min(),
1367 stats.combined_bounds_avg(), stats.combined_bounds_l2_norm()));
1371 absl::StrFormat(
"There are %i nonzero diagonal coefficients in "
1372 "the objective matrix.",
1373 stats.objective_matrix_num_nonzeros()));
1377 "Absolute values of nonzero objective matrix elements: largest=%f, "
1378 "smallest=%f, avg=%f",
1379 stats.objective_matrix_abs_max(), stats.objective_matrix_abs_min(),
1380 stats.objective_matrix_abs_avg()));
1382 SOLVER_LOG(&logger_, absl::StrFormat(
"Absolute values of objective vector "
1383 "elements: largest=%f, smallest=%f, "
1384 "avg=%f, l2_norm=%f",
1385 stats.objective_vector_abs_max(),
1386 stats.objective_vector_abs_min(),
1387 stats.objective_vector_abs_avg(),
1388 stats.objective_vector_l2_norm()));
1392 "Gaps between variable upper and lower bounds: #finite=%i of %i, "
1393 "largest=%f, smallest=%f, avg=%f",
1394 stats.variable_bound_gaps_num_finite(), stats.num_variables(),
1395 stats.variable_bound_gaps_max(), stats.variable_bound_gaps_min(),
1396 stats.variable_bound_gaps_avg()));
1399double PreprocessSolver::InitialPrimalWeight(
1400 const PrimalDualHybridGradientParams& params,
1401 const double l2_norm_primal_linear_objective,
1402 const double l2_norm_constraint_bounds)
const {
1403 if (params.has_initial_primal_weight()) {
1404 return params.initial_primal_weight();
1406 if (l2_norm_primal_linear_objective > 0.0 &&
1407 l2_norm_constraint_bounds > 0.0) {
1413 return l2_norm_primal_linear_objective / l2_norm_constraint_bounds;
1419PrimalAndDualSolution PreprocessSolver::RecoverOriginalSolution(
1420 PrimalAndDualSolution working_solution)
const {
1421 glop::ProblemSolution glop_solution(glop::RowIndex{0}, glop::ColIndex{0});
1422 if (presolve_info_.has_value()) {
1426 glop_solution = internal::ComputeStatuses(Qp(), working_solution);
1429 working_solution.primal_solution);
1431 working_solution.dual_solution);
1432 if (presolve_info_.has_value()) {
1433 glop_solution.primal_values =
1434 glop::DenseRow(working_solution.primal_solution.begin(),
1435 working_solution.primal_solution.end());
1436 glop_solution.dual_values =
1437 glop::DenseColumn(working_solution.dual_solution.begin(),
1438 working_solution.dual_solution.end());
1442 if (presolve_info_->presolved_problem_was_maximization) {
1443 for (glop::RowIndex i{0};
i < glop_solution.dual_values.size(); ++
i) {
1444 glop_solution.dual_values[
i] *= -1;
1447 presolve_info_->preprocessor.RecoverSolution(&glop_solution);
1450 Eigen::Map<Eigen::VectorXd>(glop_solution.primal_values.data(),
1451 glop_solution.primal_values.size().value());
1453 Eigen::Map<Eigen::VectorXd>(glop_solution.dual_values.data(),
1454 glop_solution.dual_values.size().value());
1461 presolve_info_->sharded_original_qp.Qp().objective_scaling_factor;
1470 return working_solution;
1474void SetActiveSetInformation(
const ShardedQuadraticProgram& sharded_qp,
1475 const VectorXd& primal_solution,
1476 const VectorXd& dual_solution,
1477 const VectorXd& primal_start_point,
1478 const VectorXd& dual_start_point,
1479 PointMetadata& metadata) {
1481 CHECK_EQ(dual_solution.size(), sharded_qp.DualSize());
1482 CHECK_EQ(primal_start_point.size(), sharded_qp.PrimalSize());
1483 CHECK_EQ(dual_start_point.size(), sharded_qp.DualSize());
1485 const QuadraticProgram& qp = sharded_qp.Qp();
1486 metadata.set_active_primal_variable_count(
1487 static_cast<int64_t
>(sharded_qp.PrimalSharder().ParallelSumOverShards(
1488 [&](
const Sharder::Shard& shard) {
1489 const auto primal_shard = shard(primal_solution);
1490 const auto lower_bound_shard = shard(qp.variable_lower_bounds);
1491 const auto upper_bound_shard = shard(qp.variable_upper_bounds);
1492 return (primal_shard.array() > lower_bound_shard.array() &&
1493 primal_shard.array() < upper_bound_shard.array())
1500 metadata.set_active_primal_variable_change(
1501 static_cast<int64_t
>(sharded_qp.PrimalSharder().ParallelSumOverShards(
1502 [&](
const Sharder::Shard& shard) {
1503 const auto primal_shard = shard(primal_solution);
1504 const auto primal_start_shard = shard(primal_start_point);
1505 const auto lower_bound_shard = shard(qp.variable_lower_bounds);
1506 const auto upper_bound_shard = shard(qp.variable_upper_bounds);
1507 return ((primal_shard.array() > lower_bound_shard.array() &&
1508 primal_shard.array() < upper_bound_shard.array()) !=
1509 (primal_start_shard.array() > lower_bound_shard.array() &&
1510 primal_start_shard.array() < upper_bound_shard.array()))
1514 metadata.set_active_dual_variable_count(
1515 static_cast<int64_t
>(sharded_qp.DualSharder().ParallelSumOverShards(
1516 [&](
const Sharder::Shard& shard) {
1517 const auto dual_shard = shard(dual_solution);
1518 const auto lower_bound_shard = shard(qp.constraint_lower_bounds);
1519 const auto upper_bound_shard = shard(qp.constraint_upper_bounds);
1520 const double kInfinity = std::numeric_limits<double>::infinity();
1521 return (dual_shard.array() != 0.0 ||
1522 (lower_bound_shard.array() == -kInfinity &&
1523 upper_bound_shard.array() == kInfinity))
1527 metadata.set_active_dual_variable_change(
1528 static_cast<int64_t
>(sharded_qp.DualSharder().ParallelSumOverShards(
1529 [&](
const Sharder::Shard& shard) {
1530 const auto dual_shard = shard(dual_solution);
1531 const auto dual_start_shard = shard(dual_start_point);
1532 const auto lower_bound_shard = shard(qp.constraint_lower_bounds);
1533 const auto upper_bound_shard = shard(qp.constraint_upper_bounds);
1534 const double kInfinity = std::numeric_limits<double>::infinity();
1535 return ((dual_shard.array() != 0.0 ||
1536 (lower_bound_shard.array() == -kInfinity &&
1537 upper_bound_shard.array() == kInfinity)) !=
1538 (dual_start_shard.array() != 0.0 ||
1539 (lower_bound_shard.array() == -kInfinity &&
1540 upper_bound_shard.array() == kInfinity)))
1545void PreprocessSolver::AddPointMetadata(
1546 const PrimalDualHybridGradientParams& params,
1547 const VectorXd& primal_solution,
const VectorXd& dual_solution,
1548 PointType point_type,
const VectorXd& last_primal_start_point,
1549 const VectorXd& last_dual_start_point, IterationStats& stats)
const {
1550 PointMetadata metadata;
1551 metadata.set_point_type(point_type);
1552 std::vector<int> random_projection_seeds(
1553 params.random_projection_seeds().begin(),
1554 params.random_projection_seeds().end());
1556 random_projection_seeds, metadata);
1557 if (point_type != POINT_TYPE_ITERATE_DIFFERENCE) {
1558 SetActiveSetInformation(sharded_qp_, primal_solution, dual_solution,
1559 last_primal_start_point, last_dual_start_point,
1562 *stats.add_point_metadata() = metadata;
1565std::optional<TerminationReasonAndPointType>
1566PreprocessSolver::UpdateIterationStatsAndCheckTermination(
1567 const PrimalDualHybridGradientParams& params,
1568 bool force_numerical_termination,
const VectorXd& working_primal_current,
1569 const VectorXd& working_dual_current,
1570 const VectorXd* working_primal_average,
1571 const VectorXd* working_dual_average,
const VectorXd* working_primal_delta,
1572 const VectorXd* working_dual_delta,
const VectorXd& last_primal_start_point,
1573 const VectorXd& last_dual_start_point,
1574 const std::atomic<bool>* interrupt_solve,
1575 const IterationType iteration_type,
const IterationStats& full_stats,
1576 IterationStats& stats) {
1577 ComputeConvergenceAndInfeasibilityFromWorkingSolution(
1578 params, working_primal_current, working_dual_current,
1579 POINT_TYPE_CURRENT_ITERATE, stats.add_convergence_information(),
1580 stats.add_infeasibility_information());
1581 AddPointMetadata(params, working_primal_current, working_dual_current,
1582 POINT_TYPE_CURRENT_ITERATE, last_primal_start_point,
1583 last_dual_start_point, stats);
1584 if (working_primal_average !=
nullptr && working_dual_average !=
nullptr) {
1585 ComputeConvergenceAndInfeasibilityFromWorkingSolution(
1586 params, *working_primal_average, *working_dual_average,
1587 POINT_TYPE_AVERAGE_ITERATE, stats.add_convergence_information(),
1588 stats.add_infeasibility_information());
1589 AddPointMetadata(params, *working_primal_average, *working_dual_average,
1590 POINT_TYPE_AVERAGE_ITERATE, last_primal_start_point,
1591 last_dual_start_point, stats);
1594 if (!presolve_info_.has_value() && working_primal_delta !=
nullptr &&
1595 working_dual_delta !=
nullptr) {
1596 ComputeConvergenceAndInfeasibilityFromWorkingSolution(
1597 params, *working_primal_delta, *working_dual_delta,
1598 POINT_TYPE_ITERATE_DIFFERENCE,
nullptr,
1599 stats.add_infeasibility_information());
1600 AddPointMetadata(params, *working_primal_delta, *working_dual_delta,
1601 POINT_TYPE_ITERATE_DIFFERENCE, last_primal_start_point,
1602 last_dual_start_point, stats);
1604 constexpr int kLogEvery = 15;
1605 absl::Time logging_time = absl::Now();
1606 if (params.verbosity_level() >= 2 &&
1607 (params.log_interval_seconds() == 0.0 ||
1608 logging_time - time_of_last_log_ >=
1609 absl::Seconds(params.log_interval_seconds()))) {
1610 if (log_counter_ == 0) {
1611 LogIterationStatsHeader(params.verbosity_level(),
1612 params.use_feasibility_polishing(), logger_);
1614 LogIterationStats(params.verbosity_level(),
1615 params.use_feasibility_polishing(), iteration_type, stats,
1616 params.termination_criteria(), original_bound_norms_,
1617 POINT_TYPE_AVERAGE_ITERATE, logger_);
1618 if (params.verbosity_level() >= 4) {
1625 params.verbosity_level(), params.use_feasibility_polishing(),
1626 iteration_type, stats, params.termination_criteria(),
1627 original_bound_norms_, POINT_TYPE_CURRENT_ITERATE, logger_);
1630 time_of_last_log_ = logging_time;
1631 if (++log_counter_ >= kLogEvery) {
1635 if (iteration_stats_callback_ !=
nullptr) {
1636 iteration_stats_callback_(
1637 {.iteration_type = iteration_type,
1638 .termination_criteria = params.termination_criteria(),
1639 .iteration_stats = stats,
1640 .bound_norms = original_bound_norms_});
1644 params.termination_criteria(), stats, original_bound_norms_,
1645 force_numerical_termination);
1646 termination.has_value()) {
1650 full_stats, interrupt_solve);
1653void PreprocessSolver::ComputeConvergenceAndInfeasibilityFromWorkingSolution(
1654 const PrimalDualHybridGradientParams& params,
1655 const VectorXd& working_primal,
const VectorXd& working_dual,
1656 PointType candidate_type, ConvergenceInformation* convergence_information,
1657 InfeasibilityInformation* infeasibility_information)
const {
1658 const TerminationCriteria::DetailedOptimalityCriteria criteria =
1660 const double primal_epsilon_ratio =
1661 EpsilonRatio(criteria.eps_optimal_primal_residual_absolute(),
1662 criteria.eps_optimal_primal_residual_relative());
1663 const double dual_epsilon_ratio =
1664 EpsilonRatio(criteria.eps_optimal_dual_residual_absolute(),
1665 criteria.eps_optimal_dual_residual_relative());
1666 if (presolve_info_.has_value()) {
1668 CHECK_NE(candidate_type, POINT_TYPE_ITERATE_DIFFERENCE);
1670 PrimalAndDualSolution original = RecoverOriginalSolution(
1671 {.primal_solution = working_primal, .dual_solution = working_dual});
1672 if (convergence_information !=
nullptr) {
1674 params, presolve_info_->sharded_original_qp,
1675 presolve_info_->trivial_col_scaling_vec,
1676 presolve_info_->trivial_row_scaling_vec, original.primal_solution,
1677 original.dual_solution, primal_epsilon_ratio, dual_epsilon_ratio,
1680 if (infeasibility_information !=
nullptr) {
1681 VectorXd primal_copy =
1683 presolve_info_->sharded_original_qp.PrimalSharder());
1690 params, presolve_info_->sharded_original_qp,
1691 presolve_info_->trivial_col_scaling_vec,
1692 presolve_info_->trivial_row_scaling_vec, primal_copy,
1693 original.dual_solution, original.primal_solution, candidate_type);
1696 if (convergence_information !=
nullptr) {
1698 params, sharded_qp_, col_scaling_vec_, row_scaling_vec_,
1699 working_primal, working_dual, primal_epsilon_ratio,
1700 dual_epsilon_ratio, candidate_type);
1702 if (infeasibility_information !=
nullptr) {
1703 VectorXd primal_copy =
1704 CloneVector(working_primal, sharded_qp_.PrimalSharder());
1707 if (candidate_type == POINT_TYPE_ITERATE_DIFFERENCE) {
1709 VectorXd dual_copy =
1710 CloneVector(working_dual, sharded_qp_.DualSharder());
1713 params, sharded_qp_, col_scaling_vec_, row_scaling_vec_,
1714 primal_copy, dual_copy, working_primal, candidate_type);
1717 params, sharded_qp_, col_scaling_vec_, row_scaling_vec_,
1718 primal_copy, working_dual, working_primal, candidate_type);
1726SolverResult PreprocessSolver::ConstructOriginalSolverResult(
1727 const PrimalDualHybridGradientParams& params, SolverResult result,
1729 const bool use_zero_primal_objective =
1730 result.solve_log.termination_reason() ==
1732 if (presolve_info_.has_value()) {
1734 PrimalAndDualSolution original_solution = RecoverOriginalSolution(
1735 {.primal_solution = std::move(result.primal_solution),
1736 .dual_solution = std::move(result.dual_solution)});
1737 result.primal_solution = std::move(original_solution.primal_solution);
1738 if (result.solve_log.termination_reason() ==
1739 TERMINATION_REASON_DUAL_INFEASIBLE) {
1741 result.primal_solution,
1746 result.dual_solution = std::move(original_solution.dual_solution);
1750 params, presolve_info_->sharded_original_qp, result.primal_solution,
1751 result.dual_solution, use_zero_primal_objective);
1753 if (result.solve_log.termination_reason() ==
1754 TERMINATION_REASON_DUAL_INFEASIBLE) {
1758 if (result.solve_log.termination_reason() ==
1759 TERMINATION_REASON_PRIMAL_INFEASIBLE) {
1762 result.reduced_costs =
1763 ReducedCosts(params, sharded_qp_, result.primal_solution,
1764 result.dual_solution, use_zero_primal_objective);
1767 result.primal_solution);
1769 result.dual_solution);
1771 col_scaling_vec_, sharded_qp_.PrimalSharder(), result.reduced_costs);
1774 switch (result.solve_log.solution_type()) {
1776 iteration_type = IterationType::kFeasibilityPolishingTermination;
1779 iteration_type = IterationType::kPresolveTermination;
1782 iteration_type = IterationType::kNormalTermination;
1785 if (iteration_stats_callback_ !=
nullptr) {
1786 iteration_stats_callback_(
1787 {.iteration_type = iteration_type,
1788 .termination_criteria = params.termination_criteria(),
1789 .iteration_stats = result.solve_log.solution_stats(),
1790 .bound_norms = original_bound_norms_});
1793 if (params.verbosity_level() >= 1) {
1798 SOLVER_LOG(&logger,
"Final solution stats:");
1799 LogIterationStatsHeader(params.verbosity_level(),
1800 params.use_feasibility_polishing(), logger);
1801 LogIterationStats(params.verbosity_level(),
1802 params.use_feasibility_polishing(), iteration_type,
1803 result.solve_log.solution_stats(),
1804 params.termination_criteria(), original_bound_norms_,
1805 result.solve_log.solution_type(), logger);
1807 result.solve_log.solution_stats(), result.solve_log.solution_type());
1808 if (convergence_info.has_value()) {
1809 if (std::isfinite(convergence_info->corrected_dual_objective())) {
1810 SOLVER_LOG(&logger,
"Dual objective after infeasibility correction: ",
1811 convergence_info->corrected_dual_objective());
1818Solver::Solver(
const PrimalDualHybridGradientParams& params,
1819 VectorXd starting_primal_solution,
1820 VectorXd starting_dual_solution,
const double initial_step_size,
1821 const double initial_primal_weight,
1822 PreprocessSolver* preprocess_solver)
1824 current_primal_solution_(
std::move(starting_primal_solution)),
1825 current_dual_solution_(
std::move(starting_dual_solution)),
1826 primal_average_(&preprocess_solver->ShardedWorkingQp().PrimalSharder()),
1827 dual_average_(&preprocess_solver->ShardedWorkingQp().DualSharder()),
1828 step_size_(initial_step_size),
1829 primal_weight_(initial_primal_weight),
1830 preprocess_solver_(preprocess_solver) {}
1832Solver::NextSolutionAndDelta Solver::ComputeNextPrimalSolution(
1833 double primal_step_size)
const {
1834 const int64_t primal_size = ShardedWorkingQp().PrimalSize();
1835 NextSolutionAndDelta result = {
1836 .value = VectorXd(primal_size),
1837 .delta = VectorXd(primal_size),
1839 const QuadraticProgram& qp = WorkingQp();
1848 ShardedWorkingQp().PrimalSharder().ParallelForEachShard(
1849 [&](
const Sharder::Shard& shard) {
1853 const VectorXd diagonal_scaling =
1855 shard(qp.objective_matrix->diagonal()).array() +
1857 shard(result.value) =
1858 (shard(current_primal_solution_) -
1860 (shard(qp.objective_vector) - shard(current_dual_product_)))
1862 .cwiseQuotient(diagonal_scaling)
1863 .cwiseMin(shard(qp.variable_upper_bounds))
1864 .cwiseMax(shard(qp.variable_lower_bounds));
1867 shard(result.value) =
1868 (shard(current_primal_solution_) -
1870 (shard(qp.objective_vector) - shard(current_dual_product_)))
1871 .cwiseMin(shard(qp.variable_upper_bounds))
1872 .cwiseMax(shard(qp.variable_lower_bounds));
1874 shard(result.delta) =
1875 shard(result.value) - shard(current_primal_solution_);
1880Solver::NextSolutionAndDelta Solver::ComputeNextDualSolution(
1881 double dual_step_size,
double extrapolation_factor,
1882 const NextSolutionAndDelta& next_primal_solution,
1883 const VectorXd* absl_nullable next_primal_product)
const {
1884 const int64_t dual_size = ShardedWorkingQp().DualSize();
1885 NextSolutionAndDelta result = {
1886 .value = VectorXd(dual_size),
1887 .delta = VectorXd(dual_size),
1889 const QuadraticProgram& qp = WorkingQp();
1890 std::optional<VectorXd> extrapolated_primal;
1891 if (!next_primal_product) {
1892 extrapolated_primal.emplace(ShardedWorkingQp().PrimalSize());
1893 ShardedWorkingQp().PrimalSharder().ParallelForEachShard(
1894 [&](
const Sharder::Shard& shard) {
1895 shard(*extrapolated_primal) =
1896 (shard(next_primal_solution.value) +
1897 extrapolation_factor * shard(next_primal_solution.delta));
1900 ShardedWorkingQp().TransposedConstraintMatrixSharder().ParallelForEachShard(
1901 [&](
const Sharder::Shard& shard) {
1903 if (next_primal_product) {
1904 CHECK(current_primal_product_.has_value());
1905 temp = shard(current_dual_solution_) -
1907 (-extrapolation_factor * shard(*current_primal_product_) +
1908 (extrapolation_factor + 1) * shard(*next_primal_product));
1910 temp = shard(current_dual_solution_) -
1912 shard(ShardedWorkingQp().TransposedConstraintMatrix())
1914 extrapolated_primal.value();
1921 shard(result.value) =
1922 VectorXd::Zero(temp.size())
1924 dual_step_size * shard(qp.constraint_upper_bounds))
1926 dual_step_size * shard(qp.constraint_lower_bounds));
1927 shard(result.delta) =
1928 (shard(result.value) - shard(current_dual_solution_));
1933std::pair<double, double> Solver::ComputeMovementTerms(
1934 const VectorXd& delta_primal,
const VectorXd& delta_dual)
const {
1935 return {
SquaredNorm(delta_primal, ShardedWorkingQp().PrimalSharder()),
1936 SquaredNorm(delta_dual, ShardedWorkingQp().DualSharder())};
1939double Solver::ComputeMovement(
const VectorXd& delta_primal,
1940 const VectorXd& delta_dual)
const {
1941 const auto [primal_squared_norm, dual_squared_norm] =
1942 ComputeMovementTerms(delta_primal, delta_dual);
1943 return (0.5 * primal_weight_ * primal_squared_norm) +
1944 (0.5 / primal_weight_) * dual_squared_norm;
1947double Solver::ComputeNonlinearity(
const VectorXd& delta_primal,
1948 const VectorXd& next_dual_product)
const {
1951 return ShardedWorkingQp().PrimalSharder().ParallelSumOverShards(
1952 [&](
const Sharder::Shard& shard) {
1953 return -shard(delta_primal)
1954 .dot(shard(next_dual_product) -
1955 shard(current_dual_product_));
1959void Solver::SetCurrentPrimalAndDualProducts() {
1960 if (params_.linesearch_rule() ==
1961 PrimalDualHybridGradientParams::MALITSKY_POCK_LINESEARCH_RULE) {
1963 ShardedWorkingQp().TransposedConstraintMatrix(),
1964 current_primal_solution_,
1965 ShardedWorkingQp().TransposedConstraintMatrixSharder());
1967 current_primal_product_.reset();
1970 WorkingQp().constraint_matrix, current_dual_solution_,
1971 ShardedWorkingQp().ConstraintMatrixSharder());
1974IterationStats Solver::CreateSimpleIterationStats(
1975 RestartChoice restart_used)
const {
1976 IterationStats stats;
1977 double num_kkt_passes_per_rejected_step = 1.0;
1978 if (params_.linesearch_rule() ==
1979 PrimalDualHybridGradientParams::MALITSKY_POCK_LINESEARCH_RULE) {
1980 num_kkt_passes_per_rejected_step = 0.5;
1982 stats.set_iteration_number(iterations_completed_);
1983 stats.set_cumulative_rejected_steps(num_rejected_steps_);
1986 stats.set_cumulative_kkt_matrix_passes(iterations_completed_ +
1987 num_kkt_passes_per_rejected_step *
1988 num_rejected_steps_);
1989 stats.set_cumulative_time_sec(preprocessing_time_sec_ + timer_.Get());
1990 stats.set_restart_used(restart_used);
1991 stats.set_step_size(step_size_);
1992 stats.set_primal_weight(primal_weight_);
1996double Solver::DistanceTraveledFromLastStart(
1997 const VectorXd& primal_solution,
const VectorXd& dual_solution)
const {
1998 return std::sqrt((0.5 * primal_weight_) *
2000 last_primal_start_point_,
2001 ShardedWorkingQp().PrimalSharder()) +
2002 (0.5 / primal_weight_) *
2004 ShardedWorkingQp().DualSharder()));
2007LocalizedLagrangianBounds Solver::ComputeLocalizedBoundsAtCurrent()
const {
2008 const double distance_traveled_by_current = DistanceTraveledFromLastStart(
2009 current_primal_solution_, current_dual_solution_);
2011 ShardedWorkingQp(), current_primal_solution_, current_dual_solution_,
2012 PrimalDualNorm::kEuclideanNorm, primal_weight_,
2013 distance_traveled_by_current,
2014 current_primal_product_.has_value()
2015 ? ¤t_primal_product_.value()
2017 ¤t_dual_product_, params_.use_diagonal_qp_trust_region_solver(),
2018 params_.diagonal_qp_trust_region_solver_tolerance());
2021LocalizedLagrangianBounds Solver::ComputeLocalizedBoundsAtAverage()
const {
2024 VectorXd average_primal = PrimalAverage();
2025 VectorXd average_dual = DualAverage();
2027 const double distance_traveled_by_average =
2028 DistanceTraveledFromLastStart(average_primal, average_dual);
2031 ShardedWorkingQp(), average_primal, average_dual,
2032 PrimalDualNorm::kEuclideanNorm, primal_weight_,
2033 distance_traveled_by_average,
2035 params_.use_diagonal_qp_trust_region_solver(),
2036 params_.diagonal_qp_trust_region_solver_tolerance());
2039bool AverageHasBetterPotential(
2040 const LocalizedLagrangianBounds& local_bounds_at_average,
2041 const LocalizedLagrangianBounds& local_bounds_at_current) {
2042 return BoundGap(local_bounds_at_average) /
2043 MathUtil::Square(local_bounds_at_average.radius) <
2044 BoundGap(local_bounds_at_current) /
2045 MathUtil::Square(local_bounds_at_current.radius);
2048double NormalizedGap(
2049 const LocalizedLagrangianBounds& local_bounds_at_candidate) {
2050 const double distance_traveled_by_candidate =
2051 local_bounds_at_candidate.radius;
2052 return BoundGap(local_bounds_at_candidate) / distance_traveled_by_candidate;
2056bool Solver::ShouldDoAdaptiveRestartHeuristic(
2057 double candidate_normalized_gap)
const {
2058 const double gap_reduction_ratio =
2059 candidate_normalized_gap / normalized_gap_at_last_restart_;
2060 if (gap_reduction_ratio < params_.sufficient_reduction_for_restart()) {
2063 if (gap_reduction_ratio < params_.necessary_reduction_for_restart() &&
2064 candidate_normalized_gap > normalized_gap_at_last_trial_) {
2072RestartChoice Solver::DetermineDistanceBasedRestartChoice()
const {
2074 if (primal_average_.NumTerms() == 0) {
2076 }
else if (distance_based_restart_info_.length_of_last_restart_period == 0) {
2079 const int restart_period_length = primal_average_.NumTerms();
2080 const double distance_moved_this_restart_period_by_average =
2081 DistanceTraveledFromLastStart(primal_average_.ComputeAverage(),
2082 dual_average_.ComputeAverage());
2083 const double distance_moved_last_restart_period =
2084 distance_based_restart_info_.distance_moved_last_restart_period;
2090 if ((distance_moved_this_restart_period_by_average / restart_period_length) <
2091 params_.sufficient_reduction_for_restart() *
2092 (distance_moved_last_restart_period /
2093 distance_based_restart_info_.length_of_last_restart_period)) {
2096 if (AverageHasBetterPotential(ComputeLocalizedBoundsAtAverage(),
2097 ComputeLocalizedBoundsAtCurrent())) {
2107RestartChoice Solver::ChooseRestartToApply(
const bool is_major_iteration) {
2108 if (!primal_average_.HasNonzeroWeight() &&
2109 !dual_average_.HasNonzeroWeight()) {
2118 const int restart_length = primal_average_.NumTerms();
2119 if (restart_length >= iterations_completed_ / 2 &&
2120 params_.restart_strategy() ==
2121 PrimalDualHybridGradientParams::ADAPTIVE_HEURISTIC) {
2122 if (AverageHasBetterPotential(ComputeLocalizedBoundsAtAverage(),
2123 ComputeLocalizedBoundsAtCurrent())) {
2129 if (is_major_iteration) {
2130 switch (params_.restart_strategy()) {
2131 case PrimalDualHybridGradientParams::NO_RESTARTS:
2133 case PrimalDualHybridGradientParams::EVERY_MAJOR_ITERATION:
2135 case PrimalDualHybridGradientParams::ADAPTIVE_HEURISTIC: {
2136 const LocalizedLagrangianBounds local_bounds_at_average =
2137 ComputeLocalizedBoundsAtAverage();
2138 const LocalizedLagrangianBounds local_bounds_at_current =
2139 ComputeLocalizedBoundsAtCurrent();
2140 double normalized_gap;
2142 if (AverageHasBetterPotential(local_bounds_at_average,
2143 local_bounds_at_current)) {
2144 normalized_gap = NormalizedGap(local_bounds_at_average);
2147 normalized_gap = NormalizedGap(local_bounds_at_current);
2150 if (ShouldDoAdaptiveRestartHeuristic(normalized_gap)) {
2153 normalized_gap_at_last_trial_ = normalized_gap;
2157 case PrimalDualHybridGradientParams::ADAPTIVE_DISTANCE_BASED: {
2158 return DetermineDistanceBasedRestartChoice();
2161 LOG(FATAL) <<
"Unrecognized restart_strategy "
2162 << params_.restart_strategy();
2170VectorXd Solver::PrimalAverage()
const {
2171 if (primal_average_.HasNonzeroWeight()) {
2172 return primal_average_.ComputeAverage();
2174 return current_primal_solution_;
2178VectorXd Solver::DualAverage()
const {
2179 if (dual_average_.HasNonzeroWeight()) {
2180 return dual_average_.ComputeAverage();
2182 return current_dual_solution_;
2186double Solver::ComputeNewPrimalWeight()
const {
2187 const double primal_distance =
2188 Distance(current_primal_solution_, last_primal_start_point_,
2189 ShardedWorkingQp().PrimalSharder());
2190 const double dual_distance =
2191 Distance(current_dual_solution_, last_dual_start_point_,
2192 ShardedWorkingQp().DualSharder());
2197 constexpr double kNonzeroTol = 1.0e-10;
2198 if (primal_distance <= kNonzeroTol || primal_distance >= 1.0 / kNonzeroTol ||
2199 dual_distance <= kNonzeroTol || dual_distance >= 1.0 / kNonzeroTol) {
2200 return primal_weight_;
2202 const double smoothing_param = params_.primal_weight_update_smoothing();
2203 const double unsmoothed_new_primal_weight = dual_distance / primal_distance;
2204 const double new_primal_weight =
2205 std::exp(smoothing_param * std::log(unsmoothed_new_primal_weight) +
2206 (1.0 - smoothing_param) * std::log(primal_weight_));
2207 if (params_.verbosity_level() >= 4) {
2208 SOLVER_LOG(&preprocess_solver_->Logger(),
"New computed primal weight is ",
2209 new_primal_weight,
" at iteration ", iterations_completed_);
2211 return new_primal_weight;
2214SolverResult Solver::PickSolutionAndConstructSolverResult(
2215 VectorXd primal_solution, VectorXd dual_solution,
2216 const IterationStats& stats, TerminationReason termination_reason,
2217 PointType output_type, SolveLog solve_log)
const {
2218 switch (output_type) {
2220 AssignVector(current_primal_solution_, ShardedWorkingQp().PrimalSharder(),
2222 AssignVector(current_dual_solution_, ShardedWorkingQp().DualSharder(),
2226 AssignVector(current_primal_delta_, ShardedWorkingQp().PrimalSharder(),
2228 AssignVector(current_dual_delta_, ShardedWorkingQp().DualSharder(),
2239 return ConstructSolverResult(
2240 std::move(primal_solution), std::move(dual_solution), stats,
2241 termination_reason, output_type, std::move(solve_log));
2244void Solver::ApplyRestartChoice(
const RestartChoice restart_to_apply) {
2245 switch (restart_to_apply) {
2250 if (params_.verbosity_level() >= 4) {
2252 "Restarted to current on iteration ", iterations_completed_,
2253 " after ", primal_average_.NumTerms(),
" iterations");
2257 if (params_.verbosity_level() >= 4) {
2259 "Restarted to average on iteration ", iterations_completed_,
2260 " after ", primal_average_.NumTerms(),
" iterations");
2262 current_primal_solution_ = primal_average_.ComputeAverage();
2263 current_dual_solution_ = dual_average_.ComputeAverage();
2264 SetCurrentPrimalAndDualProducts();
2267 primal_weight_ = ComputeNewPrimalWeight();
2268 ratio_last_two_step_sizes_ = 1;
2269 if (params_.restart_strategy() ==
2270 PrimalDualHybridGradientParams::ADAPTIVE_HEURISTIC) {
2273 const LocalizedLagrangianBounds local_bounds_at_last_restart =
2274 ComputeLocalizedBoundsAtCurrent();
2275 const double distance_traveled_since_last_restart =
2276 local_bounds_at_last_restart.radius;
2277 normalized_gap_at_last_restart_ =
BoundGap(local_bounds_at_last_restart) /
2278 distance_traveled_since_last_restart;
2279 normalized_gap_at_last_trial_ = std::numeric_limits<double>::infinity();
2280 }
else if (params_.restart_strategy() ==
2281 PrimalDualHybridGradientParams::ADAPTIVE_DISTANCE_BASED) {
2283 distance_based_restart_info_ = {
2284 .distance_moved_last_restart_period = DistanceTraveledFromLastStart(
2285 current_primal_solution_, current_dual_solution_),
2286 .length_of_last_restart_period = primal_average_.NumTerms()};
2288 primal_average_.Clear();
2289 dual_average_.Clear();
2290 AssignVector(current_primal_solution_, ShardedWorkingQp().PrimalSharder(),
2291 last_primal_start_point_);
2292 AssignVector(current_dual_solution_, ShardedWorkingQp().DualSharder(),
2293 last_dual_start_point_);
2300IterationStats AddWorkStats(IterationStats stats,
2301 const IterationStats& additional_work_stats) {
2302 stats.set_iteration_number(stats.iteration_number() +
2303 additional_work_stats.iteration_number());
2304 stats.set_cumulative_kkt_matrix_passes(
2305 stats.cumulative_kkt_matrix_passes() +
2306 additional_work_stats.cumulative_kkt_matrix_passes());
2307 stats.set_cumulative_rejected_steps(
2308 stats.cumulative_rejected_steps() +
2309 additional_work_stats.cumulative_rejected_steps());
2310 stats.set_cumulative_time_sec(stats.cumulative_time_sec() +
2311 additional_work_stats.cumulative_time_sec());
2318IterationStats WorkFromFeasibilityPolishing(
const SolveLog& solve_log) {
2319 IterationStats result;
2320 for (
const FeasibilityPolishingDetails& feasibility_polishing_detail :
2321 solve_log.feasibility_polishing_details()) {
2322 result = AddWorkStats(std::move(result),
2323 feasibility_polishing_detail.solution_stats());
2328bool TerminationReasonIsInterrupted(
const TerminationReason reason) {
2332bool TerminationReasonIsWorkLimitNotInterrupted(
2333 const TerminationReason reason) {
2341bool TerminationReasonIsWorkLimit(
const TerminationReason reason) {
2342 return TerminationReasonIsWorkLimitNotInterrupted(reason) ||
2343 TerminationReasonIsInterrupted(reason);
2346bool DoFeasibilityPolishingAfterLimitsReached(
2347 const PrimalDualHybridGradientParams& params,
2348 const TerminationReason reason) {
2349 if (TerminationReasonIsWorkLimitNotInterrupted(reason)) {
2350 return params.apply_feasibility_polishing_after_limits_reached();
2352 if (TerminationReasonIsInterrupted(reason)) {
2353 return params.apply_feasibility_polishing_if_solver_is_interrupted();
2358std::optional<SolverResult> Solver::MajorIterationAndTerminationCheck(
2359 const IterationType iteration_type,
const bool force_numerical_termination,
2360 const std::atomic<bool>* interrupt_solve,
2361 const IterationStats& work_from_feasibility_polishing,
2362 SolveLog& solve_log) {
2363 const int major_iteration_cycle =
2364 iterations_completed_ % params_.major_iteration_frequency();
2365 const bool is_major_iteration =
2366 major_iteration_cycle == 0 && iterations_completed_ > 0;
2371 : ChooseRestartToApply(is_major_iteration);
2372 IterationStats stats = CreateSimpleIterationStats(restart);
2373 IterationStats full_work_stats =
2374 AddWorkStats(stats, work_from_feasibility_polishing);
2375 std::optional<TerminationReasonAndPointType> simple_termination_reason =
2377 full_work_stats, interrupt_solve);
2378 const bool check_termination =
2379 major_iteration_cycle % params_.termination_check_frequency() == 0 ||
2380 simple_termination_reason.has_value() || force_numerical_termination;
2382 DCHECK(!is_major_iteration || check_termination);
2383 if (check_termination) {
2387 VectorXd primal_average = PrimalAverage();
2388 VectorXd dual_average = DualAverage();
2390 const std::optional<TerminationReasonAndPointType>
2391 maybe_termination_reason =
2392 preprocess_solver_->UpdateIterationStatsAndCheckTermination(
2393 params_, force_numerical_termination, current_primal_solution_,
2394 current_dual_solution_,
2395 primal_average_.HasNonzeroWeight() ? &primal_average :
nullptr,
2396 dual_average_.HasNonzeroWeight() ? &dual_average :
nullptr,
2397 current_primal_delta_.size() > 0 ? ¤t_primal_delta_
2399 current_dual_delta_.size() > 0 ? ¤t_dual_delta_ :
nullptr,
2400 last_primal_start_point_, last_dual_start_point_,
2401 interrupt_solve, iteration_type, full_work_stats, stats);
2402 if (params_.record_iteration_stats()) {
2403 *solve_log.add_iteration_stats() = stats;
2406 if (maybe_termination_reason.has_value()) {
2407 if (iteration_type == IterationType::kNormal &&
2408 DoFeasibilityPolishingAfterLimitsReached(
2409 params_, maybe_termination_reason->reason)) {
2410 const std::optional<SolverResult> feasibility_result =
2411 TryFeasibilityPolishing(
2412 iterations_completed_ / kFeasibilityIterationFraction,
2413 interrupt_solve, solve_log);
2414 if (feasibility_result.has_value()) {
2415 LOG(INFO) <<
"Returning result from feasibility polishing after "
2417 return *feasibility_result;
2420 IterationStats terminating_full_stats =
2421 AddWorkStats(stats, work_from_feasibility_polishing);
2422 return PickSolutionAndConstructSolverResult(
2423 std::move(primal_average), std::move(dual_average),
2424 terminating_full_stats, maybe_termination_reason->reason,
2425 maybe_termination_reason->type, std::move(solve_log));
2427 }
else if (params_.record_iteration_stats()) {
2429 *solve_log.add_iteration_stats() = stats;
2431 ApplyRestartChoice(restart);
2432 return std::nullopt;
2435void Solver::ResetAverageToCurrent() {
2436 primal_average_.Clear();
2437 dual_average_.Clear();
2438 primal_average_.Add(current_primal_solution_, 1.0);
2439 dual_average_.Add(current_dual_solution_, 1.0);
2442void Solver::LogNumericalTermination(
const Eigen::VectorXd& primal_delta,
2443 const Eigen::VectorXd& dual_delta)
const {
2444 if (params_.verbosity_level() >= 2) {
2445 auto [primal_squared_norm, dual_squared_norm] =
2446 ComputeMovementTerms(primal_delta, dual_delta);
2448 "Forced numerical termination at iteration ",
2449 iterations_completed_,
" with primal delta squared norm ",
2450 primal_squared_norm,
" dual delta squared norm ",
2451 dual_squared_norm,
" primal weight ", primal_weight_);
2455void Solver::LogInnerIterationLimitHit()
const {
2457 "WARNING: Inner iteration limit reached at iteration ",
2458 iterations_completed_);
2461InnerStepOutcome Solver::TakeMalitskyPockStep() {
2462 InnerStepOutcome outcome = InnerStepOutcome::kSuccessful;
2463 const double primal_step_size = step_size_ / primal_weight_;
2464 NextSolutionAndDelta next_primal_solution =
2465 ComputeNextPrimalSolution(primal_step_size);
2470 double dilating_coeff =
2471 1 + (params_.malitsky_pock_parameters().step_size_interpolation() *
2472 (sqrt(1 + ratio_last_two_step_sizes_) - 1));
2473 double new_primal_step_size = primal_step_size * dilating_coeff;
2474 double step_size_downscaling =
2475 params_.malitsky_pock_parameters().step_size_downscaling_factor();
2476 double contraction_factor =
2477 params_.malitsky_pock_parameters().linesearch_contraction_factor();
2478 const double dual_weight = primal_weight_ * primal_weight_;
2479 int inner_iterations = 0;
2480 VectorXd next_primal_product(current_dual_solution_.size());
2481 ShardedWorkingQp().TransposedConstraintMatrixSharder().ParallelForEachShard(
2482 [&](
const Sharder::Shard& shard) {
2483 shard(next_primal_product) =
2484 shard(ShardedWorkingQp().TransposedConstraintMatrix()).transpose() *
2485 next_primal_solution.value;
2488 for (
bool accepted_step =
false; !accepted_step; ++inner_iterations) {
2489 if (inner_iterations >= 60) {
2490 LogInnerIterationLimitHit();
2491 ResetAverageToCurrent();
2492 outcome = InnerStepOutcome::kForceNumericalTermination;
2495 const double new_last_two_step_sizes_ratio =
2496 new_primal_step_size / primal_step_size;
2497 NextSolutionAndDelta next_dual_solution = ComputeNextDualSolution(
2498 dual_weight * new_primal_step_size, new_last_two_step_sizes_ratio,
2499 next_primal_solution, &next_primal_product);
2502 WorkingQp().constraint_matrix, next_dual_solution.value,
2503 ShardedWorkingQp().ConstraintMatrixSharder());
2504 double delta_dual_norm =
2505 Norm(next_dual_solution.delta, ShardedWorkingQp().DualSharder());
2506 double delta_dual_prod_norm =
2507 Distance(current_dual_product_, next_dual_product,
2508 ShardedWorkingQp().PrimalSharder());
2509 if (primal_weight_ * new_primal_step_size * delta_dual_prod_norm <=
2510 contraction_factor * delta_dual_norm) {
2512 step_size_ = new_primal_step_size * primal_weight_;
2513 ratio_last_two_step_sizes_ = new_last_two_step_sizes_ratio;
2518 if (!primal_average_.HasNonzeroWeight()) {
2519 primal_average_.Add(
2520 current_primal_solution_,
2521 new_primal_step_size * new_last_two_step_sizes_ratio);
2524 current_primal_solution_ = std::move(next_primal_solution.value);
2525 current_dual_solution_ = std::move(next_dual_solution.value);
2526 current_dual_product_ = std::move(next_dual_product);
2527 current_primal_product_ = std::move(next_primal_product);
2528 primal_average_.Add(current_primal_solution_,
2529 new_primal_step_size);
2530 dual_average_.Add(current_dual_solution_,
2531 new_primal_step_size);
2532 const double movement =
2533 ComputeMovement(next_primal_solution.delta, next_dual_solution.delta);
2534 if (movement == 0.0) {
2535 LogNumericalTermination(next_primal_solution.delta,
2536 next_dual_solution.delta);
2537 ResetAverageToCurrent();
2538 outcome = InnerStepOutcome::kForceNumericalTermination;
2539 }
else if (movement > kDivergentMovement) {
2540 LogNumericalTermination(next_primal_solution.delta,
2541 next_dual_solution.delta);
2542 outcome = InnerStepOutcome::kForceNumericalTermination;
2544 current_primal_delta_ = std::move(next_primal_solution.delta);
2545 current_dual_delta_ = std::move(next_dual_solution.delta);
2548 new_primal_step_size = step_size_downscaling * new_primal_step_size;
2552 num_rejected_steps_ += inner_iterations;
2556InnerStepOutcome Solver::TakeAdaptiveStep() {
2557 InnerStepOutcome outcome = InnerStepOutcome::kSuccessful;
2558 int inner_iterations = 0;
2559 for (
bool accepted_step =
false; !accepted_step; ++inner_iterations) {
2560 if (inner_iterations >= 60) {
2561 LogInnerIterationLimitHit();
2562 ResetAverageToCurrent();
2563 outcome = InnerStepOutcome::kForceNumericalTermination;
2566 const double primal_step_size = step_size_ / primal_weight_;
2567 const double dual_step_size = step_size_ * primal_weight_;
2568 NextSolutionAndDelta next_primal_solution =
2569 ComputeNextPrimalSolution(primal_step_size);
2570 NextSolutionAndDelta next_dual_solution = ComputeNextDualSolution(
2571 dual_step_size, 1.0, next_primal_solution);
2572 const double movement =
2573 ComputeMovement(next_primal_solution.delta, next_dual_solution.delta);
2574 if (movement == 0.0) {
2575 LogNumericalTermination(next_primal_solution.delta,
2576 next_dual_solution.delta);
2577 ResetAverageToCurrent();
2578 outcome = InnerStepOutcome::kForceNumericalTermination;
2580 }
else if (movement > kDivergentMovement) {
2581 LogNumericalTermination(next_primal_solution.delta,
2582 next_dual_solution.delta);
2583 outcome = InnerStepOutcome::kForceNumericalTermination;
2587 WorkingQp().constraint_matrix, next_dual_solution.value,
2588 ShardedWorkingQp().ConstraintMatrixSharder());
2589 const double nonlinearity =
2590 ComputeNonlinearity(next_primal_solution.delta, next_dual_product);
2593 const double step_size_limit =
2594 nonlinearity > 0 ? movement / nonlinearity
2595 : std::numeric_limits<double>::infinity();
2597 if (step_size_ <= step_size_limit) {
2598 current_primal_solution_ = std::move(next_primal_solution.value);
2599 current_dual_solution_ = std::move(next_dual_solution.value);
2600 current_dual_product_ = std::move(next_dual_product);
2601 current_primal_product_.reset();
2602 current_primal_delta_ = std::move(next_primal_solution.delta);
2603 current_dual_delta_ = std::move(next_dual_solution.delta);
2604 primal_average_.Add(current_primal_solution_, step_size_);
2605 dual_average_.Add(current_dual_solution_, step_size_);
2606 accepted_step =
true;
2608 const double total_steps_attempted =
2609 num_rejected_steps_ + inner_iterations + iterations_completed_ + 1;
2614 const double first_term =
2615 std::isinf(step_size_limit)
2617 : (1 - std::pow(total_steps_attempted + 1.0,
2618 -params_.adaptive_linesearch_parameters()
2619 .step_size_reduction_exponent())) *
2621 const double second_term =
2622 (1 + std::pow(total_steps_attempted + 1.0,
2623 -params_.adaptive_linesearch_parameters()
2624 .step_size_growth_exponent())) *
2635 step_size_ = std::min(first_term, second_term);
2638 num_rejected_steps_ += inner_iterations - 1;
2642InnerStepOutcome Solver::TakeConstantSizeStep() {
2643 const double primal_step_size = step_size_ / primal_weight_;
2644 const double dual_step_size = step_size_ * primal_weight_;
2645 NextSolutionAndDelta next_primal_solution =
2646 ComputeNextPrimalSolution(primal_step_size);
2647 NextSolutionAndDelta next_dual_solution = ComputeNextDualSolution(
2648 dual_step_size, 1.0, next_primal_solution);
2649 const double movement =
2650 ComputeMovement(next_primal_solution.delta, next_dual_solution.delta);
2651 if (movement == 0.0) {
2652 LogNumericalTermination(next_primal_solution.delta,
2653 next_dual_solution.delta);
2654 ResetAverageToCurrent();
2655 return InnerStepOutcome::kForceNumericalTermination;
2656 }
else if (movement > kDivergentMovement) {
2657 LogNumericalTermination(next_primal_solution.delta,
2658 next_dual_solution.delta);
2659 return InnerStepOutcome::kForceNumericalTermination;
2662 WorkingQp().constraint_matrix, next_dual_solution.value,
2663 ShardedWorkingQp().ConstraintMatrixSharder());
2664 current_primal_solution_ = std::move(next_primal_solution.value);
2665 current_dual_solution_ = std::move(next_dual_solution.value);
2666 current_dual_product_ = std::move(next_dual_product);
2667 current_primal_product_.reset();
2668 current_primal_delta_ = std::move(next_primal_solution.delta);
2669 current_dual_delta_ = std::move(next_dual_solution.delta);
2670 primal_average_.Add(current_primal_solution_, step_size_);
2671 dual_average_.Add(current_dual_solution_, step_size_);
2672 return InnerStepOutcome::kSuccessful;
2675IterationStats Solver::TotalWorkSoFar(
const SolveLog& solve_log)
const {
2676 IterationStats stats = CreateSimpleIterationStats(RESTART_CHOICE_NO_RESTART);
2677 IterationStats full_stats =
2678 AddWorkStats(stats, WorkFromFeasibilityPolishing(solve_log));
2682FeasibilityPolishingDetails BuildFeasibilityPolishingDetails(
2683 PolishingPhaseType phase_type,
int iteration_count,
2684 const PrimalDualHybridGradientParams& params,
const SolveLog& solve_log) {
2685 FeasibilityPolishingDetails details;
2686 details.set_polishing_phase_type(phase_type);
2687 details.set_main_iteration_count(iteration_count);
2688 *details.mutable_params() = params;
2689 details.set_termination_reason(solve_log.termination_reason());
2690 details.set_iteration_count(solve_log.iteration_count());
2691 details.set_solve_time_sec(solve_log.solve_time_sec());
2692 *details.mutable_solution_stats() = solve_log.solution_stats();
2693 details.set_solution_type(solve_log.solution_type());
2694 absl::c_copy(solve_log.iteration_stats(),
2695 google::protobuf::RepeatedPtrFieldBackInserter(
2696 details.mutable_iteration_stats()));
2700std::optional<SolverResult> Solver::TryFeasibilityPolishing(
2701 const int iteration_limit,
const std::atomic<bool>* interrupt_solve,
2702 SolveLog& solve_log) {
2703 TerminationCriteria::DetailedOptimalityCriteria optimality_criteria =
2706 VectorXd average_primal = PrimalAverage();
2707 VectorXd average_dual = DualAverage();
2709 ConvergenceInformation first_convergence_info;
2710 preprocess_solver_->ComputeConvergenceAndInfeasibilityFromWorkingSolution(
2711 params_, average_primal, average_dual, POINT_TYPE_AVERAGE_ITERATE,
2712 &first_convergence_info,
nullptr);
2718 std::optional<TerminationReasonAndPointType> simple_termination_reason =
2720 TotalWorkSoFar(solve_log),
2722 if (!(simple_termination_reason.has_value() &&
2723 DoFeasibilityPolishingAfterLimitsReached(
2724 params_, simple_termination_reason->reason))) {
2725 if (params_.verbosity_level() >= 2) {
2727 "Skipping feasibility polishing because the objective gap "
2730 return std::nullopt;
2734 if (params_.verbosity_level() >= 2) {
2736 "Starting primal feasibility polishing");
2738 SolverResult primal_result = TryPrimalPolishing(
2739 std::move(average_primal), iteration_limit, interrupt_solve, solve_log);
2741 if (params_.verbosity_level() >= 2) {
2743 &preprocess_solver_->Logger(),
2744 "Primal feasibility polishing termination reason: ",
2747 if (TerminationReasonIsWorkLimit(
2748 primal_result.solve_log.termination_reason())) {
2751 std::optional<TerminationReasonAndPointType> simple_termination_reason =
2753 TotalWorkSoFar(solve_log),
2755 if (!(simple_termination_reason.has_value() &&
2756 DoFeasibilityPolishingAfterLimitsReached(
2757 params_, simple_termination_reason->reason))) {
2758 return std::nullopt;
2760 }
else if (primal_result.solve_log.termination_reason() !=
2761 TERMINATION_REASON_OPTIMAL) {
2768 "WARNING: Primal feasibility polishing terminated with error ",
2769 primal_result.solve_log.termination_reason());
2770 return std::nullopt;
2773 if (params_.verbosity_level() >= 2) {
2775 "Starting dual feasibility polishing");
2777 SolverResult dual_result = TryDualPolishing(
2778 std::move(average_dual), iteration_limit, interrupt_solve, solve_log);
2780 if (params_.verbosity_level() >= 2) {
2782 &preprocess_solver_->Logger(),
2783 "Dual feasibility polishing termination reason: ",
2787 IterationStats full_stats = TotalWorkSoFar(solve_log);
2788 std::optional<TerminationReasonAndPointType> simple_termination_reason =
2791 if (TerminationReasonIsWorkLimit(
2792 dual_result.solve_log.termination_reason())) {
2795 if (simple_termination_reason.has_value() &&
2796 DoFeasibilityPolishingAfterLimitsReached(
2797 params_, simple_termination_reason->reason)) {
2798 preprocess_solver_->ComputeConvergenceAndInfeasibilityFromWorkingSolution(
2799 params_, primal_result.primal_solution, dual_result.dual_solution,
2800 POINT_TYPE_FEASIBILITY_POLISHING_SOLUTION,
2801 full_stats.add_convergence_information(),
nullptr);
2802 return ConstructSolverResult(
2803 std::move(primal_result.primal_solution),
2804 std::move(dual_result.dual_solution), full_stats,
2805 simple_termination_reason->reason,
2806 POINT_TYPE_FEASIBILITY_POLISHING_SOLUTION, solve_log);
2808 return std::nullopt;
2810 }
else if (dual_result.solve_log.termination_reason() !=
2811 TERMINATION_REASON_OPTIMAL) {
2816 "WARNING: Dual feasibility polishing terminated with error ",
2817 dual_result.solve_log.termination_reason());
2818 return std::nullopt;
2821 preprocess_solver_->ComputeConvergenceAndInfeasibilityFromWorkingSolution(
2822 params_, primal_result.primal_solution, dual_result.dual_solution,
2823 POINT_TYPE_FEASIBILITY_POLISHING_SOLUTION,
2824 full_stats.add_convergence_information(),
nullptr);
2825 if (params_.verbosity_level() >= 2) {
2827 "solution stats for polished solution:");
2828 LogIterationStatsHeader(params_.verbosity_level(),
2830 preprocess_solver_->Logger());
2831 LogIterationStats(params_.verbosity_level(),
2833 IterationType::kFeasibilityPolishingTermination,
2834 full_stats, params_.termination_criteria(),
2835 preprocess_solver_->OriginalBoundNorms(),
2836 POINT_TYPE_FEASIBILITY_POLISHING_SOLUTION,
2837 preprocess_solver_->Logger());
2839 std::optional<TerminationReasonAndPointType> earned_termination =
2842 preprocess_solver_->OriginalBoundNorms(),
2844 if (earned_termination.has_value() ||
2845 (simple_termination_reason.has_value() &&
2846 DoFeasibilityPolishingAfterLimitsReached(
2847 params_, simple_termination_reason->reason))) {
2848 return ConstructSolverResult(
2849 std::move(primal_result.primal_solution),
2850 std::move(dual_result.dual_solution), full_stats,
2851 earned_termination.has_value() ? earned_termination->reason
2852 : simple_termination_reason->reason,
2853 POINT_TYPE_FEASIBILITY_POLISHING_SOLUTION, solve_log);
2862 return std::nullopt;
2865TerminationCriteria ReduceWorkLimitsByPreviousWork(
2866 TerminationCriteria criteria,
const int iteration_limit,
2867 const IterationStats& previous_work,
2868 bool apply_feasibility_polishing_after_limits_reached) {
2869 if (apply_feasibility_polishing_after_limits_reached) {
2870 criteria.set_iteration_limit(iteration_limit);
2871 criteria.set_kkt_matrix_pass_limit(std::numeric_limits<double>::infinity());
2872 criteria.set_time_sec_limit(std::numeric_limits<double>::infinity());
2874 criteria.set_iteration_limit(std::max(
2875 0, std::min(iteration_limit, criteria.iteration_limit() -
2876 previous_work.iteration_number())));
2877 criteria.set_kkt_matrix_pass_limit(
2878 std::max(0.0, criteria.kkt_matrix_pass_limit() -
2879 previous_work.cumulative_kkt_matrix_passes()));
2880 criteria.set_time_sec_limit(std::max(
2881 0.0, criteria.time_sec_limit() - previous_work.cumulative_time_sec()));
2886SolverResult Solver::TryPrimalPolishing(
2887 VectorXd starting_primal_solution,
const int iteration_limit,
2888 const std::atomic<bool>* interrupt_solve, SolveLog& solve_log) {
2889 PrimalDualHybridGradientParams primal_feasibility_params = params_;
2890 *primal_feasibility_params.mutable_termination_criteria() =
2891 ReduceWorkLimitsByPreviousWork(
2892 params_.termination_criteria(), iteration_limit,
2893 TotalWorkSoFar(solve_log),
2894 params_.apply_feasibility_polishing_after_limits_reached());
2895 if (params_.apply_feasibility_polishing_if_solver_is_interrupted()) {
2896 interrupt_solve =
nullptr;
2901 SetZero(ShardedWorkingQp().PrimalSharder(), objective);
2902 preprocess_solver_->SwapObjectiveVector(objective);
2904 TerminationCriteria::DetailedOptimalityCriteria criteria =
2906 const double kInfinity = std::numeric_limits<double>::infinity();
2907 criteria.set_eps_optimal_dual_residual_absolute(kInfinity);
2908 criteria.set_eps_optimal_dual_residual_relative(kInfinity);
2909 criteria.set_eps_optimal_objective_gap_absolute(kInfinity);
2910 criteria.set_eps_optimal_objective_gap_relative(kInfinity);
2911 *primal_feasibility_params.mutable_termination_criteria()
2912 ->mutable_detailed_optimality_criteria() = criteria;
2915 VectorXd primal_feasibility_starting_dual;
2916 SetZero(ShardedWorkingQp().DualSharder(), primal_feasibility_starting_dual);
2917 Solver primal_solver(primal_feasibility_params,
2918 std::move(starting_primal_solution),
2919 std::move(primal_feasibility_starting_dual), step_size_,
2920 primal_weight_, preprocess_solver_);
2921 SolveLog primal_solve_log;
2925 SolverResult primal_result = primal_solver.Solve(
2926 IterationType::kPrimalFeasibility, interrupt_solve, primal_solve_log);
2929 preprocess_solver_->SwapObjectiveVector(objective);
2931 *solve_log.add_feasibility_polishing_details() =
2932 BuildFeasibilityPolishingDetails(
2933 POLISHING_PHASE_TYPE_PRIMAL_FEASIBILITY, iterations_completed_,
2934 primal_feasibility_params, primal_result.solve_log);
2935 return primal_result;
2938VectorXd MapFiniteValuesToZero(
const Sharder& sharder,
const VectorXd&
input) {
2939 VectorXd output(
input.size());
2940 const auto make_finite_values_zero = [](
const double x) {
2941 return std::isfinite(x) ? 0.0 :
x;
2943 sharder.ParallelForEachShard([&](
const Sharder::Shard& shard) {
2944 shard(output) = shard(
input).unaryExpr(make_finite_values_zero);
2949SolverResult Solver::TryDualPolishing(VectorXd starting_dual_solution,
2950 const int iteration_limit,
2951 const std::atomic<bool>* interrupt_solve,
2952 SolveLog& solve_log) {
2953 PrimalDualHybridGradientParams dual_feasibility_params = params_;
2954 *dual_feasibility_params.mutable_termination_criteria() =
2955 ReduceWorkLimitsByPreviousWork(
2956 params_.termination_criteria(), iteration_limit,
2957 TotalWorkSoFar(solve_log),
2958 params_.apply_feasibility_polishing_after_limits_reached());
2959 if (params_.apply_feasibility_polishing_if_solver_is_interrupted()) {
2960 interrupt_solve =
nullptr;
2966 VectorXd constraint_lower_bounds = MapFiniteValuesToZero(
2967 ShardedWorkingQp().DualSharder(), WorkingQp().constraint_lower_bounds);
2968 VectorXd constraint_upper_bounds = MapFiniteValuesToZero(
2969 ShardedWorkingQp().DualSharder(), WorkingQp().constraint_upper_bounds);
2970 VectorXd variable_lower_bounds = MapFiniteValuesToZero(
2971 ShardedWorkingQp().PrimalSharder(), WorkingQp().variable_lower_bounds);
2972 VectorXd variable_upper_bounds = MapFiniteValuesToZero(
2973 ShardedWorkingQp().PrimalSharder(), WorkingQp().variable_upper_bounds);
2974 preprocess_solver_->SwapConstraintBounds(constraint_lower_bounds,
2975 constraint_upper_bounds);
2976 preprocess_solver_->SwapVariableBounds(variable_lower_bounds,
2977 variable_upper_bounds);
2979 TerminationCriteria::DetailedOptimalityCriteria criteria =
2981 const double kInfinity = std::numeric_limits<double>::infinity();
2982 criteria.set_eps_optimal_primal_residual_absolute(kInfinity);
2983 criteria.set_eps_optimal_primal_residual_relative(kInfinity);
2984 criteria.set_eps_optimal_objective_gap_absolute(kInfinity);
2985 criteria.set_eps_optimal_objective_gap_relative(kInfinity);
2986 *dual_feasibility_params.mutable_termination_criteria()
2987 ->mutable_detailed_optimality_criteria() = criteria;
2990 VectorXd dual_feasibility_starting_primal;
2991 SetZero(ShardedWorkingQp().PrimalSharder(), dual_feasibility_starting_primal);
2992 Solver dual_solver(dual_feasibility_params,
2993 std::move(dual_feasibility_starting_primal),
2994 std::move(starting_dual_solution), step_size_,
2995 primal_weight_, preprocess_solver_);
2996 SolveLog dual_solve_log;
3000 SolverResult dual_result = dual_solver.Solve(IterationType::kDualFeasibility,
3001 interrupt_solve, dual_solve_log);
3004 preprocess_solver_->SwapConstraintBounds(constraint_lower_bounds,
3005 constraint_upper_bounds);
3006 preprocess_solver_->SwapVariableBounds(variable_lower_bounds,
3007 variable_upper_bounds);
3008 *solve_log.add_feasibility_polishing_details() =
3009 BuildFeasibilityPolishingDetails(
3010 POLISHING_PHASE_TYPE_DUAL_FEASIBILITY, iterations_completed_,
3011 dual_feasibility_params, dual_result.solve_log);
3015SolverResult Solver::Solve(
const IterationType iteration_type,
3016 const std::atomic<bool>* interrupt_solve,
3017 SolveLog solve_log) {
3018 preprocessing_time_sec_ = solve_log.preprocessing_time_sec();
3020 last_primal_start_point_ =
3021 CloneVector(current_primal_solution_, ShardedWorkingQp().PrimalSharder());
3022 last_dual_start_point_ =
3023 CloneVector(current_dual_solution_, ShardedWorkingQp().DualSharder());
3027 ratio_last_two_step_sizes_ = 1;
3028 SetCurrentPrimalAndDualProducts();
3032 bool force_numerical_termination =
false;
3034 int next_feasibility_polishing_iteration = 100;
3036 num_rejected_steps_ = 0;
3038 IterationStats work_from_feasibility_polishing =
3039 WorkFromFeasibilityPolishing(solve_log);
3040 for (iterations_completed_ = 0;; ++iterations_completed_) {
3044 const std::optional<SolverResult> maybe_result =
3045 MajorIterationAndTerminationCheck(
3046 iteration_type, force_numerical_termination, interrupt_solve,
3047 work_from_feasibility_polishing, solve_log);
3048 if (maybe_result.has_value()) {
3049 return maybe_result.value();
3052 if (params_.use_feasibility_polishing() &&
3053 iteration_type == IterationType::kNormal &&
3054 iterations_completed_ >= next_feasibility_polishing_iteration) {
3055 const std::optional<SolverResult> feasibility_result =
3056 TryFeasibilityPolishing(
3057 iterations_completed_ / kFeasibilityIterationFraction,
3058 interrupt_solve, solve_log);
3059 if (feasibility_result.has_value()) {
3060 return *feasibility_result;
3062 next_feasibility_polishing_iteration *= 2;
3064 work_from_feasibility_polishing = WorkFromFeasibilityPolishing(solve_log);
3071 InnerStepOutcome outcome;
3072 switch (params_.linesearch_rule()) {
3073 case PrimalDualHybridGradientParams::MALITSKY_POCK_LINESEARCH_RULE:
3074 outcome = TakeMalitskyPockStep();
3076 case PrimalDualHybridGradientParams::ADAPTIVE_LINESEARCH_RULE:
3077 outcome = TakeAdaptiveStep();
3079 case PrimalDualHybridGradientParams::CONSTANT_STEP_SIZE_RULE:
3080 outcome = TakeConstantSizeStep();
3083 LOG(FATAL) <<
"Unrecognized linesearch rule "
3084 << params_.linesearch_rule();
3086 if (outcome == InnerStepOutcome::kForceNumericalTermination) {
3087 force_numerical_termination =
true;
3096 const std::atomic<bool>* interrupt_solve,
3097 std::function<
void(
const std::string&)> message_callback,
3098 IterationStatsCallback iteration_stats_callback) {
3100 interrupt_solve, std::move(message_callback),
3101 std::move(iteration_stats_callback));
3107 std::optional<PrimalAndDualSolution> initial_solution,
3108 const std::atomic<bool>* interrupt_solve,
3109 std::function<
void(
const std::string&)> message_callback,
3110 IterationStatsCallback iteration_stats_callback) {
3113 if (message_callback) {
3118 const absl::Status params_status =
3120 if (!params_status.ok()) {
3122 params_status.ToString(), logger);
3126 "constraint_matrix must be in compressed format. "
3127 "Call constraint_matrix.makeCompressed()",
3131 if (!dimensions_status.ok()) {
3133 dimensions_status.ToString(), logger);
3137 "The objective scaling factor cannot be zero.",
3141 return ErrorSolverResult(
3143 "use_feasibility_polishing is only implemented for linear programs.",
3146 PreprocessSolver solver(std::move(qp), params, &logger);
3147 return solver.PreprocessAndSolve(params, std::move(initial_solution),
3149 std::move(iteration_stats_callback));
3157 glop::RowIndex(
solution.dual_solution.size()),
3158 glop::ColIndex(
solution.primal_solution.size()));
3162 for (glop::RowIndex i{0}; i.value() <
solution.dual_solution.size(); ++i) {
3167 }
else if (
solution.dual_solution[i.value()] > 0) {
3170 }
else if (
solution.dual_solution[i.value()] < 0) {
3178 for (glop::ColIndex i{0}; i.value() <
solution.primal_solution.size(); ++i) {
3179 const bool at_lb =
solution.primal_solution[i.value()] <=
3181 const bool at_ub =
solution.primal_solution[i.value()] >=
3201 return glop_solution;
void Start()
When Start() is called multiple times, only the most recent is used.
bool use_feasibility_polishing() const
void EnableLogging(bool enable)
void SetLogToStdOut(bool enable)
Should all messages be displayed on stdout ?
void AddInfoLoggingCallback(std::function< void(const std::string &message)> callback)
Fractional SquaredNorm(const SparseColumn &v)
constexpr Fractional kInfinity
Infinity for type Fractional.
@ TERMINATION_REASON_OPTIMAL
@ TERMINATION_REASON_NUMERICAL_ERROR
TerminationReason
The reason a call to Solve() terminates.
absl::StatusOr< SolveResult > Solve(const Model &model, const SolverType solver_type, const SolveArguments &solve_args, const SolverInitArguments &init_args)
glop::ProblemSolution ComputeStatuses(const QuadraticProgram &qp, const PrimalAndDualSolution &solution)
Validation utilities for solvers.proto.
std::optional< TerminationReasonAndPointType > CheckSimpleTerminationCriteria(const TerminationCriteria &criteria, const IterationStats &stats, const std::atomic< bool > *interrupt_solve)
absl::StatusOr< MPModelProto > QpToMpModelProto(const QuadraticProgram &qp)
void SetZero(const Sharder &sharder, VectorXd &dest)
InfeasibilityInformation ComputeInfeasibilityInformation(const PrimalDualHybridGradientParams ¶ms, const ShardedQuadraticProgram &scaled_sharded_qp, const Eigen::VectorXd &col_scaling_vec, const Eigen::VectorXd &row_scaling_vec, const Eigen::VectorXd &scaled_primal_ray, const Eigen::VectorXd &scaled_dual_ray, const Eigen::VectorXd &primal_solution_for_residual_tests, PointType candidate_type)
QuadraticProgramStats ComputeStats(const ShardedQuadraticProgram &qp)
Returns a QuadraticProgramStats for a ShardedQuadraticProgram.
std::optional< ConvergenceInformation > GetConvergenceInformation(const IterationStats &stats, PointType candidate_type)
absl::Status ValidateQuadraticProgramDimensions(const QuadraticProgram &qp)
absl::Status ValidatePrimalDualHybridGradientParams(const PrimalDualHybridGradientParams ¶ms)
void SetRandomProjections(const ShardedQuadraticProgram &sharded_qp, const Eigen::VectorXd &primal_solution, const Eigen::VectorXd &dual_solution, const std::vector< int > &random_projection_seeds, PointMetadata &metadata)
double SquaredDistance(const VectorXd &vector1, const VectorXd &vector2, const Sharder &sharder)
double LInfNorm(const VectorXd &vector, const Sharder &sharder)
TerminationCriteria::DetailedOptimalityCriteria EffectiveOptimalityCriteria(const TerminationCriteria &termination_criteria)
Computes the effective optimality criteria for a TerminationCriteria.
@ POINT_TYPE_FEASIBILITY_POLISHING_SOLUTION
@ POINT_TYPE_PRESOLVER_SOLUTION
@ POINT_TYPE_AVERAGE_ITERATE
@ POINT_TYPE_ITERATE_DIFFERENCE
@ POINT_TYPE_CURRENT_ITERATE
double Distance(const VectorXd &vector1, const VectorXd &vector2, const Sharder &sharder)
VectorXd ReducedCosts(const PrimalDualHybridGradientParams ¶ms, const ShardedQuadraticProgram &sharded_qp, const VectorXd &primal_solution, const VectorXd &dual_solution, bool use_zero_primal_objective)
VectorXd TransposedMatrixVectorProduct(const Eigen::SparseMatrix< double, Eigen::ColMajor, int64_t > &matrix, const VectorXd &vector, const Sharder &sharder)
const ::std::string & PointType_Name(T value)
QuadraticProgramBoundNorms BoundNormsFromProblemStats(const QuadraticProgramStats &stats)
double EpsilonRatio(const double epsilon_absolute, const double epsilon_relative)
SingularValueAndIterations EstimateMaximumSingularValueOfConstraintMatrix(const ShardedQuadraticProgram &sharded_qp, const std::optional< VectorXd > &primal_solution, const std::optional< VectorXd > &dual_solution, const double desired_relative_error, const double failure_probability, std::mt19937 &mt_generator)
@ RESTART_CHOICE_RESTART_TO_AVERAGE
@ RESTART_CHOICE_UNSPECIFIED
@ RESTART_CHOICE_NO_RESTART
@ RESTART_CHOICE_WEIGHTED_AVERAGE_RESET
void ProjectToPrimalVariableBounds(const ShardedQuadraticProgram &sharded_qp, VectorXd &primal, const bool use_feasibility_bounds)
absl::StatusOr< QuadraticProgram > QpFromMpModelProto(const MPModelProto &proto, bool relax_integer_variables, bool include_names)
const ::std::string & TerminationReason_Name(T value)
bool IsLinearProgram(const QuadraticProgram &qp)
@ OPTIMALITY_NORM_UNSPECIFIED
@ OPTIMALITY_NORM_L_INF_COMPONENTWISE
void ProjectToDualVariableBounds(const ShardedQuadraticProgram &sharded_qp, VectorXd &dual)
std::optional< TerminationReasonAndPointType > CheckIterateTerminationCriteria(const TerminationCriteria &criteria, const IterationStats &stats, const QuadraticProgramBoundNorms &bound_norms, const bool force_numerical_termination)
void CoefficientWiseProductInPlace(const VectorXd &scale, const Sharder &sharder, VectorXd &dest)
void CoefficientWiseQuotientInPlace(const VectorXd &scale, const Sharder &sharder, VectorXd &dest)
VectorXd ZeroVector(const Sharder &sharder)
Like VectorXd::Zero(sharder.NumElements()).
VectorXd CloneVector(const VectorXd &vec, const Sharder &sharder)
SolverResult PrimalDualHybridGradient(QuadraticProgram qp, const PrimalDualHybridGradientParams ¶ms, const std::atomic< bool > *interrupt_solve, std::function< void(const std::string &)> message_callback, IterationStatsCallback iteration_stats_callback)
VectorXd OnesVector(const Sharder &sharder)
Like VectorXd::Ones(sharder.NumElements()).
bool HasValidBounds(const ShardedQuadraticProgram &sharded_qp)
ScalingVectors ApplyRescaling(const RescalingOptions &rescaling_options, ShardedQuadraticProgram &sharded_qp)
bool ObjectiveGapMet(const TerminationCriteria::DetailedOptimalityCriteria &optimality_criteria, const ConvergenceInformation &stats)
ConvergenceInformation ComputeConvergenceInformation(const PrimalDualHybridGradientParams ¶ms, const ShardedQuadraticProgram &scaled_sharded_qp, const Eigen::VectorXd &col_scaling_vec, const Eigen::VectorXd &row_scaling_vec, const Eigen::VectorXd &scaled_primal_solution, const Eigen::VectorXd &scaled_dual_solution, const double componentwise_primal_residual_offset, const double componentwise_dual_residual_offset, PointType candidate_type)
LocalizedLagrangianBounds ComputeLocalizedLagrangianBounds(const ShardedQuadraticProgram &sharded_qp, const VectorXd &primal_solution, const VectorXd &dual_solution, const PrimalDualNorm primal_dual_norm, const double primal_weight, const double radius, const VectorXd *primal_product, const VectorXd *dual_product, const bool use_diagonal_qp_trust_region_solver, const double diagonal_qp_trust_region_solver_tolerance)
double Norm(const VectorXd &vector, const Sharder &sharder)
@ TERMINATION_REASON_TIME_LIMIT
@ TERMINATION_REASON_KKT_MATRIX_PASS_LIMIT
@ TERMINATION_REASON_PRIMAL_OR_DUAL_INFEASIBLE
@ TERMINATION_REASON_INVALID_PROBLEM
@ TERMINATION_REASON_PRIMAL_INFEASIBLE
@ TERMINATION_REASON_ITERATION_LIMIT
@ TERMINATION_REASON_INVALID_PARAMETER
@ TERMINATION_REASON_OTHER
@ TERMINATION_REASON_INTERRUPTED_BY_USER
RelativeConvergenceInformation ComputeRelativeResiduals(const TerminationCriteria::DetailedOptimalityCriteria &optimality_criteria, const ConvergenceInformation &stats, const QuadraticProgramBoundNorms &bound_norms)
void AssignVector(const VectorXd &vec, const Sharder &sharder, VectorXd &dest)
double BoundGap(const LocalizedLagrangianBounds &bounds)
@ kNormal
An intermediate iteration in the "main" phase.
@ kPrimalFeasibility
An intermediate iteration during a primal feasibility polishing phase.
@ kFeasibilityPolishingTermination
Terminating with a solution found by feasibility polishing.
@ kPresolveTermination
Terminating with a solution found by presolve.
@ kDualFeasibility
An intermediate iteration during a dual feasibility polishing phase.
@ kNormalTermination
Terminating with a solution found by the "main" phase.
Select next search node to expand Select next item_i to add this new search node to the search Generate a new search node where item_i is not in the knapsack Check validity of this new partial solution(using propagators) - If valid
static int input(yyscan_t yyscanner)
Contains the solution of a LinearProgram as returned by a preprocessor.
ConstraintStatusColumn constraint_statuses
VariableStatusRow variable_statuses
ProblemStatus status
The solution status.
double objective_scaling_factor
Eigen::VectorXd variable_lower_bounds
Eigen::VectorXd constraint_lower_bounds
Eigen::SparseMatrix< double, Eigen::ColMajor, int64_t > constraint_matrix
Eigen::VectorXd variable_upper_bounds
Eigen::VectorXd constraint_upper_bounds
#define SOLVER_LOG(logger,...)