Google OR-Tools v9.14
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
primal_dual_hybrid_gradient.cc
Go to the documentation of this file.
1// Copyright 2010-2025 Google LLC
2// Licensed under the Apache License, Version 2.0 (the "License");
3// you may not use this file except in compliance with the License.
4// You may obtain a copy of the License at
5//
6// http://www.apache.org/licenses/LICENSE-2.0
7//
8// Unless required by applicable law or agreed to in writing, software
9// distributed under the License is distributed on an "AS IS" BASIS,
10// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11// See the License for the specific language governing permissions and
12// limitations under the License.
13
14// We solve a QP, which we call the "original QP", by applying preprocessing
15// including presolve and rescaling, which produces a new QP that we call the
16// "working QP". We then solve the working QP using the Primal Dual Hybrid
17// Gradient algorithm (PDHG). The optimality criteria are evaluated using the
18// original QP. There are three main modules in this file:
19// * The free function `PrimalDualHybridGradient()`, which is the user API, and
20// is responsible for input validation that doesn't use
21// ShardedQuadraticProgram, creating a `PreprocessSolver`, and calling
22// `PreprocessSolver::PreprocessAndSolve()`.
23// * The class `PreprocessSolver`, which is responsible for everything that
24// touches the original QP, including input validation that uses
25// ShardedQuadraticProgram, the preprocessing, converting solutions to the
26// working QP back to solutions to the original QP, and termination checks. It
27// also creates a `Solver` object and calls `Solver::Solve()`.
28// * The class `Solver`, which is responsible for everything that only touches
29// the working QP. It keeps a pointer to `PreprocessSolver` and calls methods
30// on it when it needs access to the original QP, e.g. termination checks.
31// When feasibility polishing is enabled the main solve's `Solver` object
32// creates additional `Solver` objects periodically to do the feasibility
33// polishing (in `Solver::TryPrimalPolishing()` and
34// `Solver::TryDualPolishing()`).
35// The main reason for having two separate classes `PreprocessSolver` and
36// `Solver` is the fact that feasibility polishing mode uses a single
37// `PreprocessSolver` object with multiple `Solver` objects.
38
40
41#include <algorithm>
42#include <atomic>
43#include <cmath>
44#include <cstdint>
45#include <functional>
46#include <limits>
47#include <optional>
48#include <random>
49#include <string>
50#include <utility>
51#include <vector>
52
53#include "Eigen/Core"
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"
68#include "ortools/base/timer.h"
86
88
89namespace {
90
91using ::Eigen::VectorXd;
92using ::operations_research::SolverLogger;
93
94using IterationStatsCallback =
95 std::function<void(const IterationCallbackInfo&)>;
96
97// Computes a `num_threads` that is capped by the problem size and `num_shards`,
98// if specified, to avoid creating unusable threads.
99int NumThreads(const int num_threads, const int num_shards,
100 const QuadraticProgram& qp, SolverLogger& logger) {
101 int capped_num_threads = num_threads;
102 if (num_shards > 0) {
103 capped_num_threads = std::min(capped_num_threads, num_shards);
104 }
105 const int64_t problem_limit = std::max(qp.variable_lower_bounds.size(),
106 qp.constraint_lower_bounds.size());
107 capped_num_threads =
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.");
114 }
115 return capped_num_threads;
116}
117
118// If `num_shards` is positive, returns it. Otherwise returns a reasonable
119// number of shards to use with `ShardedQuadraticProgram` for the given
120// `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;
124}
125
126std::string ConvergenceInformationString(
127 const ConvergenceInformation& convergence_information,
128 const RelativeConvergenceInformation& relative_information,
129 const OptimalityNorm residual_norm) {
130 constexpr absl::string_view kFormatStr =
131 "%#12.6g %#12.6g %#12.6g | %#12.6g %#12.6g %#12.6g | %#12.6g %#12.6g | "
132 "%#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(
162 kFormatStr,
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.";
176 }
177 LOG(FATAL) << "Invalid residual norm " << residual_norm << ".";
178}
179
180std::string ConvergenceInformationShortString(
181 const ConvergenceInformation& convergence_information,
182 const RelativeConvergenceInformation& relative_information,
183 const OptimalityNorm residual_norm) {
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(
203 kFormatStr,
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.";
211 }
212 LOG(FATAL) << "Invalid residual norm " << residual_norm << ".";
213}
214
215// Logs a description of `iter_stats`, based on the
216// `iter_stats.convergence_information` entry with
217// `.candidate_type()==preferred_candidate` if one exists, otherwise based on
218// the first value, if any. `termination_criteria.optimality_norm` determines
219// which residual norms from `iter_stats.convergence_information` are used.
220void LogIterationStats(int verbosity_level, bool use_feasibility_polishing,
221 IterationType iteration_type,
222 const IterationStats& iter_stats,
223 const TerminationCriteria& termination_criteria,
224 const QuadraticProgramBoundNorms& bound_norms,
225 PointType preferred_candidate, SolverLogger& logger) {
226 std::string iteration_string =
227 verbosity_level >= 3
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 =
234 GetConvergenceInformation(iter_stats, preferred_candidate);
235 if (!convergence_information.has_value() &&
236 iter_stats.convergence_information_size() > 0) {
237 convergence_information = iter_stats.convergence_information(0);
238 }
239 const char* phase_string = [&]() {
240 if (use_feasibility_polishing) {
241 switch (iteration_type) {
243 return "n ";
245 return "p ";
247 return "d ";
251 return "t ";
252 }
253 } else {
254 return "";
255 }
256 }();
257 if (convergence_information.has_value()) {
258 const char* iterate_string = [&]() {
259 if (verbosity_level >= 4) {
260 switch (convergence_information->candidate_type()) {
262 return "C ";
264 return "A ";
266 return "D ";
268 return "F ";
269 default:
270 return "? ";
271 }
272 } else {
273 return "";
274 }
275 }();
276 const RelativeConvergenceInformation relative_information =
278 EffectiveOptimalityCriteria(termination_criteria),
279 *convergence_information, bound_norms);
280 std::string convergence_string =
281 verbosity_level >= 3
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, " | ",
289 convergence_string);
290 } else {
291 // No convergence information, just log the basic work stats.
292 SOLVER_LOG(&logger, phase_string, verbosity_level >= 4 ? "? " : "",
293 iteration_string);
294 }
295}
296
297void LogIterationStatsHeader(int verbosity_level,
298 bool use_feasibility_polishing,
299 SolverLogger& logger) {
300 std::string work_labels =
301 verbosity_level >= 3
302 ? absl::StrFormat("%6s %8s %6s", "iter#", "kkt_pass", "time")
303 : absl::StrFormat("%6s %6s", "iter#", "time");
304 std::string convergence_labels =
305 verbosity_level >= 3
306 ? absl::StrFormat(
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",
310 "dual_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, " | ",
315 convergence_labels);
316}
317
318enum class InnerStepOutcome {
319 kSuccessful,
320 kForceNumericalTermination,
321};
322
323// Makes the closing changes to `solve_log` and builds a `SolverResult`.
324// NOTE: `primal_solution`, `dual_solution`, and `solve_log` are passed by
325// value. To avoid unnecessary copying, move these objects (i.e. use
326// `std::move()`).
327SolverResult ConstructSolverResult(VectorXd primal_solution,
328 VectorXd dual_solution,
329 const IterationStats& stats,
330 TerminationReason termination_reason,
331 PointType output_type, SolveLog solve_log) {
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;
337 return SolverResult{.primal_solution = std::move(primal_solution),
338 .dual_solution = std::move(dual_solution),
339 .solve_log = std::move(solve_log)};
340}
341
342// See comment at top of file.
343class PreprocessSolver {
344 public:
345 // Assumes that `qp` and `params` are valid.
346 // Note that the `qp` is intentionally passed by value.
347 // `logger` is captured, and must outlive the class.
348 // NOTE: Many `PreprocessSolver` methods accept a `params` argument. This is
349 // passed as an argument instead of stored as a member variable to support
350 // using different `params` in different contexts with the same
351 // `PreprocessSolver` object.
352 explicit PreprocessSolver(QuadraticProgram qp,
353 const PrimalDualHybridGradientParams& params,
354 SolverLogger* logger);
355
356 // Not copyable or movable (because `glop::MainLpPreprocessor` isn't).
357 PreprocessSolver(const PreprocessSolver&) = delete;
358 PreprocessSolver& operator=(const PreprocessSolver&) = delete;
359 PreprocessSolver(PreprocessSolver&&) = delete;
360 PreprocessSolver& operator=(PreprocessSolver&&) = delete;
361
362 // Zero is used if `initial_solution` is nullopt. If `interrupt_solve` is not
363 // nullptr, then the solver will periodically check if
364 // `interrupt_solve->load()` is true, in which case the solve will terminate
365 // with `TERMINATION_REASON_INTERRUPTED_BY_USER`. Ownership is not
366 // transferred. If `iteration_stats_callback` is not nullptr, then at each
367 // termination step (when iteration stats are logged),
368 // `iteration_stats_callback` will also be called with those iteration stats.
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);
374
375 // Returns a `TerminationReasonAndPointType` when the termination criteria are
376 // satisfied, otherwise returns nothing. The pointers to working_* can be
377 // nullptr if an iterate of that type is not available. For the iterate types
378 // that are available, uses the primal and dual vectors to compute solution
379 // statistics and adds them to the stats proto.
380 // `full_stats` is used for checking work-based termination criteria,
381 // but `stats` is updated with convergence information.
382 // NOTE: The primal and dual input pairs should be scaled solutions.
383 // TODO(user): Move this long argument list into a struct.
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);
396
397 // Computes solution statistics for the primal and dual input pair, which
398 // should be a scaled solution. If `convergence_information != nullptr`,
399 // populates it with convergence information, and similarly if
400 // `infeasibility_information != nullptr`, populates it with infeasibility
401 // information.
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;
407
408 // Returns a `SolverResult` for the original problem, given a `SolverResult`
409 // from the scaled or preprocessed problem. Also computes the reduced costs.
410 // NOTE: `result` is passed by value. To avoid unnecessary copying, move this
411 // object (i.e. use `std::move()`).
412 SolverResult ConstructOriginalSolverResult(
413 const PrimalDualHybridGradientParams& params, SolverResult result,
414 SolverLogger& logger) const;
415
416 const ShardedQuadraticProgram& ShardedWorkingQp() const {
417 return sharded_qp_;
418 }
419
420 // Swaps `variable_lower_bounds` and `variable_upper_bounds` with the ones on
421 // the sharded quadratic program.
422 void SwapVariableBounds(VectorXd& variable_lower_bounds,
423 VectorXd& variable_upper_bounds) {
424 sharded_qp_.SwapVariableBounds(variable_lower_bounds,
425 variable_upper_bounds);
426 }
427
428 // Swaps `constraint_lower_bounds` and `constraint_upper_bounds` with the ones
429 // on the sharded quadratic program.
430 void SwapConstraintBounds(VectorXd& constraint_lower_bounds,
431 VectorXd& constraint_upper_bounds) {
432 sharded_qp_.SwapConstraintBounds(constraint_lower_bounds,
433 constraint_upper_bounds);
434 }
435
436 // Swaps `objective` with the objective vector on the quadratic program.
437 // Swapping `objective_matrix` is not yet supported because it hasn't been
438 // needed.
439 void SwapObjectiveVector(VectorXd& objective) {
440 sharded_qp_.SwapObjectiveVector(objective);
441 }
442
443 const QuadraticProgramBoundNorms& OriginalBoundNorms() const {
444 return original_bound_norms_;
445 }
446
447 SolverLogger& Logger() { return logger_; }
448
449 private:
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())) {}
460
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;
466 };
467
468 // TODO(user): experiment with different preprocessor types.
469 static glop::GlopParameters PreprocessorParameters(
470 const PrimalDualHybridGradientParams& params);
471
472 // If presolve is enabled, moves `sharded_qp_` to
473 // `presolve_info_.sharded_original_qp` and computes the presolved linear
474 // program and installs it in `sharded_qp_`. Clears `initial_solution` if
475 // presolve is enabled. If presolve solves the problem completely returns the
476 // appropriate `TerminationReason`. Otherwise returns nullopt. If presolve
477 // is disabled or an error occurs modifies nothing and returns nullopt.
478 std::optional<TerminationReason> ApplyPresolveIfEnabled(
479 const PrimalDualHybridGradientParams& params,
480 std::optional<PrimalAndDualSolution>* initial_solution);
481
482 void ComputeAndApplyRescaling(const PrimalDualHybridGradientParams& params,
483 VectorXd& starting_primal_solution,
484 VectorXd& starting_dual_solution);
485
486 void LogQuadraticProgramStats(const QuadraticProgramStats& stats) const;
487
488 double InitialPrimalWeight(const PrimalDualHybridGradientParams& params,
489 double l2_norm_primal_linear_objective,
490 double l2_norm_constraint_bounds) const;
491
492 PrimalAndDualSolution RecoverOriginalSolution(
493 PrimalAndDualSolution working_solution) const;
494
495 // Adds one entry of `PointMetadata` to `stats` using the input solutions.
496 void AddPointMetadata(const PrimalDualHybridGradientParams& params,
497 const VectorXd& primal_solution,
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;
502
503 const QuadraticProgram& Qp() const { return sharded_qp_.Qp(); }
504
505 const int num_threads_;
506 const int num_shards_;
507
508 // The bound norms of the original problem.
509 QuadraticProgramBoundNorms original_bound_norms_;
510
511 // This is the QP that PDHG is run on. It is modified by presolve and
512 // rescaling, if those are enabled, and then serves as the
513 // `ShardedWorkingQp()` when calling `Solver::Solve()`. The original problem
514 // is available in `presolve_info_->sharded_original_qp` if
515 // `presolve_info_.has_value()`, and otherwise can be obtained by undoing the
516 // scaling of `sharded_qp_` by `col_scaling_vec_` and `row_scaling_vec_`.
517 ShardedQuadraticProgram sharded_qp_;
518
519 // Set iff presolve is enabled.
520 std::optional<PresolveInfo> presolve_info_;
521
522 // The scaling vectors that map the original (or presolved) quadratic program
523 // to the working version. See
524 // `ShardedQuadraticProgram::RescaleQuadraticProgram()` for details.
525 VectorXd col_scaling_vec_;
526 VectorXd row_scaling_vec_;
527
528 // A counter used to trigger writing iteration headers.
529 int log_counter_ = 0;
530 absl::Time time_of_last_log_ = absl::InfinitePast();
531 SolverLogger& logger_;
532 IterationStatsCallback iteration_stats_callback_;
533};
534
535// See comment at top of file.
536class Solver {
537 public:
538 // `preprocess_solver` should not be nullptr, and the `PreprocessSolver`
539 // object must outlive this `Solver` object. Ownership is not transferred.
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);
545
546 // Not copyable or movable (because there are const members).
547 Solver(const Solver&) = delete;
548 Solver& operator=(const Solver&) = delete;
549 Solver(Solver&&) = delete;
550 Solver& operator=(Solver&&) = delete;
551
552 const QuadraticProgram& WorkingQp() const { return ShardedWorkingQp().Qp(); }
553
554 const ShardedQuadraticProgram& ShardedWorkingQp() const {
555 return preprocess_solver_->ShardedWorkingQp();
556 }
557
558 // Runs PDHG iterations on the instance that has been initialized in `Solver`.
559 // If `interrupt_solve` is not nullptr, then the solver will periodically
560 // check if `interrupt_solve->load()` is true, in which case the solve will
561 // terminate with `TERMINATION_REASON_INTERRUPTED_BY_USER`. Ownership is not
562 // transferred.
563 // `solve_log` should contain initial problem statistics.
564 // On return, `SolveResult.reduced_costs` will be empty, and the solution will
565 // be to the preprocessed/scaled problem.
566 SolverResult Solve(IterationType iteration_type,
567 const std::atomic<bool>* interrupt_solve,
568 SolveLog solve_log);
569
570 private:
571 struct NextSolutionAndDelta {
572 VectorXd value;
573 // `delta` is `value` - current_solution.
574 VectorXd delta;
575 };
576
577 struct DistanceBasedRestartInfo {
578 double distance_moved_last_restart_period;
579 int length_of_last_restart_period;
580 };
581
582 // Movement terms (weighted squared norms of primal and dual deltas) larger
583 // than this cause termination because iterates are diverging, and likely to
584 // cause infinite and NaN values.
585 constexpr static double kDivergentMovement = 1.0e100;
586
587 // The total number of iterations in feasibility polishing is at most
588 // `4 * iterations_completed_ / kFeasibilityIterationFraction`.
589 // One factor of two is because there are both primal and dual feasibility
590 // polishing phases, and the other factor of two is because
591 // `next_feasibility_polishing_iteration` increases by a factor of 2 each
592 // feasibility polishing phase, so the sum of iteration limits is at most
593 // twice the last value.
594 constexpr static int kFeasibilityIterationFraction = 8;
595
596 // Attempts to solve primal and dual feasibility subproblems starting at the
597 // average iterate, for at most `iteration_limit` iterations each. If
598 // successful, returns a `SolverResult`, otherwise nullopt. Appends
599 // information about the polishing phases to
600 // `solve_log.feasibility_polishing_details`.
601 std::optional<SolverResult> TryFeasibilityPolishing(
602 int iteration_limit, const std::atomic<bool>* interrupt_solve,
603 SolveLog& solve_log);
604
605 // Tries to find primal feasibility, adds the solve log details to
606 // `solve_log.feasibility_polishing_details`, and returns the result.
607 SolverResult TryPrimalPolishing(VectorXd starting_primal_solution,
608 int iteration_limit,
609 const std::atomic<bool>* interrupt_solve,
610 SolveLog& solve_log);
611
612 // Tries to find dual feasibility, adds the solve log details to
613 // `solve_log.feasibility_polishing_details`, and returns the result.
614 SolverResult TryDualPolishing(VectorXd starting_dual_solution,
615 int iteration_limit,
616 const std::atomic<bool>* interrupt_solve,
617 SolveLog& solve_log);
618
619 NextSolutionAndDelta ComputeNextPrimalSolution(double primal_step_size) const;
620
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;
625
626 std::pair<double, double> ComputeMovementTerms(
627 const VectorXd& delta_primal, const VectorXd& delta_dual) const;
628
629 double ComputeMovement(const VectorXd& delta_primal,
630 const VectorXd& delta_dual) const;
631
632 double ComputeNonlinearity(const VectorXd& delta_primal,
633 const VectorXd& next_dual_product) const;
634
635 // Sets current_primal_product_ and current_dual_product_ based on
636 // current_primal_solution_ and current_dual_solution_ respectively.
637 void SetCurrentPrimalAndDualProducts();
638
639 // Creates all the simple-to-compute statistics in stats.
640 IterationStats CreateSimpleIterationStats(RestartChoice restart_used) const;
641
642 // Returns the total work so far from the main iterations and feasibility
643 // polishing phases.
644 IterationStats TotalWorkSoFar(const SolveLog& solve_log) const;
645
646 RestartChoice ChooseRestartToApply(bool is_major_iteration);
647
648 VectorXd PrimalAverage() const;
649
650 VectorXd DualAverage() const;
651
652 double ComputeNewPrimalWeight() const;
653
654 // Picks the primal and dual solutions according to `output_type`, and makes
655 // the closing changes to `solve_log`. This function should only be called
656 // once when the solver is finishing its execution.
657 // NOTE: `primal_solution` and `dual_solution` are used as the output except
658 // when `output_type` is `POINT_TYPE_CURRENT_ITERATE` or
659 // `POINT_TYPE_ITERATE_DIFFERENCE`, in which case the values are computed from
660 // `Solver` data.
661 // NOTE: `primal_solution`, `dual_solution`, and `solve_log` are passed by
662 // value. To avoid unnecessary copying, move these objects (i.e. use
663 // `std::move()`).
664 SolverResult PickSolutionAndConstructSolverResult(
665 VectorXd primal_solution, VectorXd dual_solution,
666 const IterationStats& stats, TerminationReason termination_reason,
667 PointType output_type, SolveLog solve_log) const;
668
669 double DistanceTraveledFromLastStart(const VectorXd& primal_solution,
670 const VectorXd& dual_solution) const;
671
672 LocalizedLagrangianBounds ComputeLocalizedBoundsAtCurrent() const;
673
674 LocalizedLagrangianBounds ComputeLocalizedBoundsAtAverage() const;
675
676 // Applies the given `RestartChoice`. If a restart is chosen, updates the
677 // state of the algorithm accordingly and computes a new primal weight.
678 void ApplyRestartChoice(RestartChoice restart_to_apply);
679
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);
685
686 bool ShouldDoAdaptiveRestartHeuristic(double candidate_normalized_gap) const;
687
688 RestartChoice DetermineDistanceBasedRestartChoice() const;
689
690 void ResetAverageToCurrent();
691
692 void LogNumericalTermination(const Eigen::VectorXd& primal_delta,
693 const Eigen::VectorXd& dual_delta) const;
694
695 void LogInnerIterationLimitHit() const;
696
697 // Takes a step based on the Malitsky and Pock linesearch algorithm.
698 // (https://arxiv.org/pdf/1608.08883.pdf)
699 // The current implementation is provably convergent (at an optimal rate)
700 // for LP programs (provided we do not change the primal weight at every major
701 // iteration). Further, we have observed that this rule is very sensitive to
702 // the parameter choice whenever we apply the primal weight recomputation
703 // heuristic.
704 InnerStepOutcome TakeMalitskyPockStep();
705
706 // Takes a step based on the adaptive heuristic presented in Section 3.1 of
707 // https://arxiv.org/pdf/2106.04756.pdf (further generalized to QP).
708 InnerStepOutcome TakeAdaptiveStep();
709
710 // Takes a constant-size step.
711 InnerStepOutcome TakeConstantSizeStep();
712
713 const PrimalDualHybridGradientParams params_;
714
715 VectorXd current_primal_solution_;
716 VectorXd current_dual_solution_;
717 VectorXd current_primal_delta_;
718 VectorXd current_dual_delta_;
719
720 ShardedWeightedAverage primal_average_;
721 ShardedWeightedAverage dual_average_;
722
723 double step_size_;
724 double primal_weight_;
725
726 PreprocessSolver* preprocess_solver_;
727
728 // For Malitsky-Pock linesearch only: `step_size_` / previous_step_size
729 double ratio_last_two_step_sizes_;
730 // For adaptive restarts only.
731 double normalized_gap_at_last_trial_ =
732 std::numeric_limits<double>::infinity();
733 // For adaptive restarts only.
734 double normalized_gap_at_last_restart_ =
735 std::numeric_limits<double>::infinity();
736
737 // `preprocessing_time_sec_` is added to the current value of `timer_` when
738 // computing `cumulative_time_sec` in iteration stats.
739 double preprocessing_time_sec_;
740 WallTimer timer_;
741 int iterations_completed_;
742 int num_rejected_steps_;
743 // A cache of `constraint_matrix * current_primal_solution_`.
744 // Malitsky-Pock linesearch only.
745 std::optional<VectorXd> current_primal_product_;
746 // A cache of `constraint_matrix.transpose() * current_dual_solution_`.
747 VectorXd current_dual_product_;
748 // The primal point at which the algorithm was last restarted from, or
749 // the initial primal starting point if no restart has occurred.
750 VectorXd last_primal_start_point_;
751 // The dual point at which the algorithm was last restarted from, or
752 // the initial dual starting point if no restart has occurred.
753 VectorXd last_dual_start_point_;
754 // Information for deciding whether to trigger a distance-based restart.
755 // The distances are initialized to +inf to force a restart during the first
756 // major iteration check.
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,
761 };
762};
763
764PreprocessSolver::PreprocessSolver(QuadraticProgram qp,
765 const PrimalDualHybridGradientParams& params,
766 SolverLogger* logger)
767 : num_threads_(
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),
772 logger_(*logger) {}
773
774SolverResult ErrorSolverResult(const TerminationReason reason,
775 const std::string& message,
776 SolverLogger& logger) {
777 SolveLog error_log;
778 error_log.set_termination_reason(reason);
779 error_log.set_termination_string(message);
780 SOLVER_LOG(&logger,
781 "The solver did not run because of invalid input: ", message);
782 return SolverResult{.solve_log = error_log};
783}
784
785// If `check_excessively_small_values`, in addition to checking for extreme
786// values that PDLP can't handle, also check for extremely small constraint
787// bounds, variable bound gaps, and objective vector values that PDLP can handle
788// but that cause problems for other code such as glop's presolve.
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);
798 }
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, "."),
805 logger);
806 }
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.");
814 }
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, "."),
824 logger);
825 }
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, "."),
835 logger);
836 }
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);
840 }
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 "
845 "absolute value ",
846 problem_stats.combined_bounds_max(),
847 " which exceeds limit of ", kExcessiveInputValue, "."),
848 logger);
849 }
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 "
856 "absolute value ",
857 problem_stats.combined_bounds_min(),
858 " which is less than the limit of ",
859 kExcessivelySmallInputValue, "."),
860 logger);
861 }
862 if (problem_stats.combined_bounds_max() >
863 kMaxDynamicRange * problem_stats.combined_bounds_min()) {
864 SOLVER_LOG(&logger,
865 "WARNING: Combined constraint bounds vector has largest "
866 "absolute value ",
867 problem_stats.combined_bounds_max(),
868 " and smallest non-zero absolute value ",
869 problem_stats.combined_bounds_min(),
870 "; performance may suffer.");
871 }
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);
875 }
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 "
880 "absolute value ",
881 problem_stats.combined_variable_bounds_max(),
882 " which exceeds limit of ", kExcessiveInputValue, "."),
883 logger);
884 }
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 "
892 "absolute value ",
893 problem_stats.combined_variable_bounds_min(),
894 " which is less than the limit of ",
895 kExcessivelySmallInputValue, "."),
896 logger);
897 }
898 if (problem_stats.combined_variable_bounds_max() >
899 kMaxDynamicRange * problem_stats.combined_variable_bounds_min()) {
901 &logger,
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.");
907 }
908 if (problem_stats.variable_bound_gaps_max() >
909 kMaxDynamicRange * problem_stats.variable_bound_gaps_min()) {
910 SOLVER_LOG(&logger,
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.");
916 }
917 if (std::isnan(objective_offset)) {
918 return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
919 "Objective offset is NAN.", logger);
920 }
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, "."),
927 logger);
928 }
929 if (std::isnan(problem_stats.objective_vector_l2_norm())) {
930 return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
931 "Objective vector has a NAN.", logger);
932 }
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, "."),
939 logger);
940 }
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, "."),
950 logger);
951 }
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.");
959 }
960 if (std::isnan(problem_stats.objective_matrix_l2_norm())) {
961 return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
962 "Objective matrix has a NAN.", logger);
963 }
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, "."),
970 logger);
971 }
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.");
979 }
980 return std::nullopt;
981}
982
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()),
994 logger);
995 }
996 if (std::isnan(
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);
1000 }
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,
1006 absl::StrCat(
1007 "Initial primal solution has an entry with absolute value ", norm,
1008 " which exceeds limit of ", kExcessiveInputValue),
1009 logger);
1010 }
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()),
1018 logger);
1019 }
1020 if (std::isnan(
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);
1024 }
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),
1032 logger);
1033 }
1034 return std::nullopt;
1035}
1036
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) {
1042 WallTimer timer;
1043 timer.Start();
1044 SolveLog solve_log;
1045 if (params.verbosity_level() >= 1) {
1046 SOLVER_LOG(&logger_, "Solving with PDLP parameters: ", params);
1047 }
1048 if (Qp().problem_name.has_value()) {
1049 solve_log.set_instance_name(*Qp().problem_name);
1050 }
1051 *solve_log.mutable_params() = params;
1052 sharded_qp_.ReplaceLargeConstraintBoundsWithInfinity(
1053 params.infinite_constraint_bound_threshold());
1054 if (!HasValidBounds(sharded_qp_)) {
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 == "
1060 "-inf.",
1061 logger_);
1062 }
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)
1067 .all();
1068 })) {
1069 return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
1070 "The objective is not convex (i.e., the objective "
1071 "matrix contains negative or NAN entries).",
1072 logger_);
1073 }
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;
1082 }
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;
1088 }
1089 }
1090 original_bound_norms_ = BoundNormsFromProblemStats(original_problem_stats);
1091 const std::string preprocessing_string = absl::StrCat(
1092 params.presolve_options().use_glop() ? "presolving and " : "",
1093 "rescaling:");
1094 if (params.verbosity_level() >= 1) {
1095 SOLVER_LOG(&logger_, "Problem stats before ", preprocessing_string);
1096 LogQuadraticProgramStats(solve_log.original_problem_stats());
1097 }
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()) {
1102 // Glop also feeds zero primal and dual solutions when the preprocessor
1103 // has a non-INIT status. When the preprocessor status is optimal the
1104 // vectors have length 0. When the status is something else the lengths
1105 // may be non-zero, but that's OK since we don't promise to produce a
1106 // meaningful solution in that case.
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 =
1117 CheckIterateTerminationCriteria(params.termination_criteria(),
1118 iteration_stats, original_bound_norms_,
1119 /*force_numerical_termination=*/false);
1120 if (!earned_termination.has_value()) {
1121 earned_termination = CheckSimpleTerminationCriteria(
1122 params.termination_criteria(), iteration_stats, interrupt_solve);
1123 }
1124 TerminationReason final_termination_reason;
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;
1130 } else {
1131 if (*maybe_terminate == TERMINATION_REASON_OPTIMAL) {
1132 final_termination_reason = TERMINATION_REASON_NUMERICAL_ERROR;
1133 SOLVER_LOG(
1134 &logger_,
1135 "WARNING: Presolve claimed to solve the LP optimally but the "
1136 "solution doesn't satisfy the optimality criteria.");
1137 } else {
1138 final_termination_reason = *maybe_terminate;
1139 }
1140 }
1141 return ConstructOriginalSolverResult(
1142 params,
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)),
1147 logger_);
1148 }
1149
1150 VectorXd starting_primal_solution;
1151 VectorXd starting_dual_solution;
1152 // The current solution is updated by `ComputeAndApplyRescaling`.
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);
1156 } else {
1157 SetZero(sharded_qp_.PrimalSharder(), starting_primal_solution);
1158 SetZero(sharded_qp_.DualSharder(), starting_dual_solution);
1159 }
1160 // The following projections are necessary since all our checks assume that
1161 // the primal and dual variable bounds are satisfied.
1162 ProjectToPrimalVariableBounds(sharded_qp_, starting_primal_solution);
1163 ProjectToDualVariableBounds(sharded_qp_, starting_dual_solution);
1164
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());
1171 }
1172
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,
1181 /*desired_relative_error=*/0.2, /*failure_probability=*/0.0005,
1182 random);
1183 // With high probability, `lipschitz_result.singular_value` is within
1184 // +/- `lipschitz_result.estimated_relative_error
1185 // * lipschitz_result.singular_value`
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;
1191 } else {
1192 // This initial step size is designed to err on the side of being too big.
1193 // This is because
1194 // (i) too-big steps are rejected and hence don't hurt us other than
1195 // wasting
1196 // an iteration and
1197 // (ii) the step size adjustment algorithm shrinks the step size as far as
1198 // needed in a single iteration but raises it slowly.
1199 // The tiny constant is there to keep the step size finite in the case of a
1200 // trivial LP with no constraints.
1201 step_size =
1202 1.0 /
1203 std::max(
1204 1.0e-20,
1205 solve_log.preprocessed_problem_stats().constraint_matrix_abs_max());
1206 }
1207 step_size *= params.initial_step_size_scaling();
1208
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());
1212
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_);
1219}
1220
1221glop::GlopParameters PreprocessSolver::PreprocessorParameters(
1222 const PrimalDualHybridGradientParams& params) {
1223 glop::GlopParameters glop_params;
1224 // TODO(user): Test if dualization helps or hurts performance.
1225 glop_params.set_solve_dual_problem(glop::GlopParameters::NEVER_DO);
1226 // Experiments show that this preprocessing step can hurt because it relaxes
1227 // variable bounds.
1228 glop_params.set_use_implied_free_preprocessor(false);
1229 // We do our own scaling.
1230 glop_params.set_use_scaling(false);
1231 if (params.presolve_options().has_glop_parameters()) {
1232 glop_params.MergeFrom(params.presolve_options().glop_parameters());
1233 }
1234 return glop_params;
1235}
1236
1237TerminationReason GlopStatusToTerminationReason(
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:
1253 default:
1254 SOLVER_LOG(&logger, "WARNING: Unexpected preprocessor status ",
1255 glop_status);
1257 }
1258}
1259
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;
1266 }
1267 if (!IsLinearProgram(Qp())) {
1268 SOLVER_LOG(&logger_,
1269 "WARNING: Skipping presolve, which is only supported for linear "
1270 "programs");
1271 return std::nullopt;
1272 }
1273 absl::StatusOr<MPModelProto> model = QpToMpModelProto(Qp());
1274 if (!model.ok()) {
1275 SOLVER_LOG(&logger_,
1276 "WARNING: Skipping presolve because of error converting to "
1277 "MPModelProto: ",
1278 model.status());
1279 return std::nullopt;
1280 }
1281 if (initial_solution->has_value()) {
1282 SOLVER_LOG(&logger_,
1283 "WARNING: Ignoring initial solution. Initial solutions are "
1284 "ignored when presolve is on.");
1285 initial_solution->reset();
1286 }
1287 glop::LinearProgram glop_lp;
1288 glop::MPModelProtoToLinearProgram(*model, &glop_lp);
1289 // Save RAM.
1290 model->Clear();
1291 presolve_info_.emplace(std::move(sharded_qp_), params);
1292 // To simplify our code we ignore the return value indicating whether
1293 // postprocessing is required. We simply call `RecoverSolution()`
1294 // unconditionally, which may do nothing.
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);
1300 // This will only fail if given an invalid LP, which shouldn't happen.
1301 absl::StatusOr<QuadraticProgram> presolved_qp =
1302 QpFromMpModelProto(output, /*relax_integer_variables=*/false);
1303 CHECK_OK(presolved_qp.status());
1304 // `MPModelProto` doesn't support scaling factors, so if `glop_lp` has an
1305 // `objective_scaling_factor` it won't be set in output and `presolved_qp`.
1306 // The scaling factor of `presolved_qp` isn't actually used anywhere, but we
1307 // set it for completeness.
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());
1311 // A status of `INIT` means the preprocessor created a (usually) smaller
1312 // problem that needs solving. Other statuses mean the preprocessor solved
1313 // the problem completely.
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(),
1318 logger_);
1319 }
1320 return std::nullopt;
1321}
1322
1323void PreprocessSolver::ComputeAndApplyRescaling(
1324 const PrimalDualHybridGradientParams& params,
1325 VectorXd& starting_primal_solution, VectorXd& starting_dual_solution) {
1326 ScalingVectors scaling = ApplyRescaling(
1327 RescalingOptions{.l_inf_ruiz_iterations = params.l_inf_ruiz_iterations(),
1328 .l2_norm_rescaling = params.l2_norm_rescaling()},
1329 sharded_qp_);
1330 row_scaling_vec_ = std::move(scaling.row_scaling_vec);
1331 col_scaling_vec_ = std::move(scaling.col_scaling_vec);
1332
1333 CoefficientWiseQuotientInPlace(col_scaling_vec_, sharded_qp_.PrimalSharder(),
1334 starting_primal_solution);
1335 CoefficientWiseQuotientInPlace(row_scaling_vec_, sharded_qp_.DualSharder(),
1336 starting_dual_solution);
1337}
1338
1339void PreprocessSolver::LogQuadraticProgramStats(
1340 const QuadraticProgramStats& stats) const {
1341 SOLVER_LOG(&logger_,
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) {
1347 SOLVER_LOG(&logger_,
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()));
1354 SOLVER_LOG(
1355 &logger_,
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()));
1361 SOLVER_LOG(
1362 &logger_,
1363 absl::StrFormat(
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()));
1368 }
1369 if (!IsLinearProgram(Qp())) {
1370 SOLVER_LOG(&logger_,
1371 absl::StrFormat("There are %i nonzero diagonal coefficients in "
1372 "the objective matrix.",
1373 stats.objective_matrix_num_nonzeros()));
1374 SOLVER_LOG(
1375 &logger_,
1376 absl::StrFormat(
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()));
1381 }
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()));
1389 SOLVER_LOG(
1390 &logger_,
1391 absl::StrFormat(
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()));
1397}
1398
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();
1405 }
1406 if (l2_norm_primal_linear_objective > 0.0 &&
1407 l2_norm_constraint_bounds > 0.0) {
1408 // The hand-wavy motivation for this choice is that the objective vector
1409 // has units of (objective units)/(primal units) and the constraint
1410 // bounds vector has units of (objective units)/(dual units),
1411 // therefore this ratio has units (dual units)/(primal units). By
1412 // dimensional analysis, these are the same units as the primal weight.
1413 return l2_norm_primal_linear_objective / l2_norm_constraint_bounds;
1414 } else {
1415 return 1.0;
1416 }
1417}
1418
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()) {
1423 // We compute statuses relative to the working problem so we can detect when
1424 // variables are at their bounds without floating-point roundoff induced by
1425 // scaling.
1426 glop_solution = internal::ComputeStatuses(Qp(), working_solution);
1427 }
1428 CoefficientWiseProductInPlace(col_scaling_vec_, sharded_qp_.PrimalSharder(),
1429 working_solution.primal_solution);
1430 CoefficientWiseProductInPlace(row_scaling_vec_, sharded_qp_.DualSharder(),
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());
1439 // We got the working QP by calling `LinearProgramToMPModelProto()` and
1440 // `QpFromMpModelProto()`. We need to negate the duals if the LP resulting
1441 // from presolve was a max problem.
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;
1445 }
1446 }
1447 presolve_info_->preprocessor.RecoverSolution(&glop_solution);
1448 PrimalAndDualSolution solution;
1449 solution.primal_solution =
1450 Eigen::Map<Eigen::VectorXd>(glop_solution.primal_values.data(),
1451 glop_solution.primal_values.size().value());
1452 solution.dual_solution =
1453 Eigen::Map<Eigen::VectorXd>(glop_solution.dual_values.data(),
1454 glop_solution.dual_values.size().value());
1455 // We called `QpToMpModelProto()` and `MPModelProtoToLinearProgram()` to
1456 // convert our original QP into input for glop's preprocessor. The former
1457 // multiplies the objective vector by `objective_scaling_factor`, which
1458 // multiplies the duals by that factor as well. To undo this we divide by
1459 // `objective_scaling_factor`.
1460 solution.dual_solution /=
1461 presolve_info_->sharded_original_qp.Qp().objective_scaling_factor;
1462 // Glop's preprocessor sometimes violates the primal bounds constraints. To
1463 // be safe we project both primal and dual.
1464 ProjectToPrimalVariableBounds(presolve_info_->sharded_original_qp,
1465 solution.primal_solution);
1466 ProjectToDualVariableBounds(presolve_info_->sharded_original_qp,
1467 solution.dual_solution);
1468 return solution;
1469 } else {
1470 return working_solution;
1471 }
1472}
1473
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) {
1480 CHECK_EQ(primal_solution.size(), sharded_qp.PrimalSize());
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());
1484
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())
1494 .count();
1495 })));
1496
1497 // Most of the computation from the previous `ParallelSumOverShards` is
1498 // duplicated here. However the overhead shouldn't be too large, and using
1499 // `ParallelSumOverShards` is simpler than just using `ParallelForEachShard`.
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()))
1511 .count();
1512 })));
1513
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))
1524 .count();
1525 })));
1526
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)))
1541 .count();
1542 })));
1543}
1544
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());
1555 SetRandomProjections(sharded_qp_, primal_solution, dual_solution,
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,
1560 metadata);
1561 }
1562 *stats.add_point_metadata() = metadata;
1563}
1564
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);
1592 }
1593 // Undoing presolve doesn't work for iterate differences.
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);
1603 }
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_);
1613 }
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) {
1619 // If the convergence information doesn't contain an average iterate, the
1620 // previous `LogIterationStats()` will report some other iterate (usually
1621 // the current one) so don't repeat logging it.
1622 if (GetConvergenceInformation(stats, POINT_TYPE_AVERAGE_ITERATE) !=
1623 std::nullopt) {
1624 LogIterationStats(
1625 params.verbosity_level(), params.use_feasibility_polishing(),
1626 iteration_type, stats, params.termination_criteria(),
1627 original_bound_norms_, POINT_TYPE_CURRENT_ITERATE, logger_);
1628 }
1629 }
1630 time_of_last_log_ = logging_time;
1631 if (++log_counter_ >= kLogEvery) {
1632 log_counter_ = 0;
1633 }
1634 }
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_});
1641 }
1642
1643 if (const auto termination = CheckIterateTerminationCriteria(
1644 params.termination_criteria(), stats, original_bound_norms_,
1645 force_numerical_termination);
1646 termination.has_value()) {
1647 return termination;
1648 }
1649 return CheckSimpleTerminationCriteria(params.termination_criteria(),
1650 full_stats, interrupt_solve);
1651}
1652
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 =
1659 EffectiveOptimalityCriteria(params.termination_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()) {
1667 // Undoing presolve doesn't make sense for iterate differences.
1668 CHECK_NE(candidate_type, POINT_TYPE_ITERATE_DIFFERENCE);
1669
1670 PrimalAndDualSolution original = RecoverOriginalSolution(
1671 {.primal_solution = working_primal, .dual_solution = working_dual});
1672 if (convergence_information != nullptr) {
1673 *convergence_information = ComputeConvergenceInformation(
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,
1678 candidate_type);
1679 }
1680 if (infeasibility_information != nullptr) {
1681 VectorXd primal_copy =
1682 CloneVector(original.primal_solution,
1683 presolve_info_->sharded_original_qp.PrimalSharder());
1684 ProjectToPrimalVariableBounds(presolve_info_->sharded_original_qp,
1685 primal_copy,
1686 /*use_feasibility_bounds=*/true);
1687 // Only iterate differences should be able to violate the dual variable
1688 // bounds, and iterate differences aren't used when presolve is enabled.
1689 *infeasibility_information = ComputeInfeasibilityInformation(
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);
1694 }
1695 } else {
1696 if (convergence_information != nullptr) {
1697 *convergence_information = ComputeConvergenceInformation(
1698 params, sharded_qp_, col_scaling_vec_, row_scaling_vec_,
1699 working_primal, working_dual, primal_epsilon_ratio,
1700 dual_epsilon_ratio, candidate_type);
1701 }
1702 if (infeasibility_information != nullptr) {
1703 VectorXd primal_copy =
1704 CloneVector(working_primal, sharded_qp_.PrimalSharder());
1705 ProjectToPrimalVariableBounds(sharded_qp_, primal_copy,
1706 /*use_feasibility_bounds=*/true);
1707 if (candidate_type == POINT_TYPE_ITERATE_DIFFERENCE) {
1708 // Iterate differences might not satisfy the dual variable bounds.
1709 VectorXd dual_copy =
1710 CloneVector(working_dual, sharded_qp_.DualSharder());
1711 ProjectToDualVariableBounds(sharded_qp_, dual_copy);
1712 *infeasibility_information = ComputeInfeasibilityInformation(
1713 params, sharded_qp_, col_scaling_vec_, row_scaling_vec_,
1714 primal_copy, dual_copy, working_primal, candidate_type);
1715 } else {
1716 *infeasibility_information = ComputeInfeasibilityInformation(
1717 params, sharded_qp_, col_scaling_vec_, row_scaling_vec_,
1718 primal_copy, working_dual, working_primal, candidate_type);
1719 }
1720 }
1721 }
1722}
1723
1724// `result` is used both as the input and as the temporary that will be
1725// returned.
1726SolverResult PreprocessSolver::ConstructOriginalSolverResult(
1727 const PrimalDualHybridGradientParams& params, SolverResult result,
1728 SolverLogger& logger) const {
1729 const bool use_zero_primal_objective =
1730 result.solve_log.termination_reason() ==
1732 if (presolve_info_.has_value()) {
1733 // Transform the solutions so they match the original unscaled problem.
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) {
1740 ProjectToPrimalVariableBounds(presolve_info_->sharded_original_qp,
1741 result.primal_solution,
1742 /*use_feasibility_bounds=*/true);
1743 }
1744 // If presolve is enabled, no projection of the dual is performed when
1745 // checking for infeasibility.
1746 result.dual_solution = std::move(original_solution.dual_solution);
1747 // `RecoverOriginalSolution` doesn't recover reduced costs so we need to
1748 // compute them with respect to the original problem.
1749 result.reduced_costs = ReducedCosts(
1750 params, presolve_info_->sharded_original_qp, result.primal_solution,
1751 result.dual_solution, use_zero_primal_objective);
1752 } else {
1753 if (result.solve_log.termination_reason() ==
1754 TERMINATION_REASON_DUAL_INFEASIBLE) {
1755 ProjectToPrimalVariableBounds(sharded_qp_, result.primal_solution,
1756 /*use_feasibility_bounds=*/true);
1757 }
1758 if (result.solve_log.termination_reason() ==
1759 TERMINATION_REASON_PRIMAL_INFEASIBLE) {
1760 ProjectToDualVariableBounds(sharded_qp_, result.dual_solution);
1761 }
1762 result.reduced_costs =
1763 ReducedCosts(params, sharded_qp_, result.primal_solution,
1764 result.dual_solution, use_zero_primal_objective);
1765 // Transform the solutions so they match the original unscaled problem.
1766 CoefficientWiseProductInPlace(col_scaling_vec_, sharded_qp_.PrimalSharder(),
1767 result.primal_solution);
1768 CoefficientWiseProductInPlace(row_scaling_vec_, sharded_qp_.DualSharder(),
1769 result.dual_solution);
1771 col_scaling_vec_, sharded_qp_.PrimalSharder(), result.reduced_costs);
1772 }
1773 IterationType iteration_type;
1774 switch (result.solve_log.solution_type()) {
1776 iteration_type = IterationType::kFeasibilityPolishingTermination;
1777 break;
1779 iteration_type = IterationType::kPresolveTermination;
1780 break;
1781 default:
1782 iteration_type = IterationType::kNormalTermination;
1783 break;
1784 }
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_});
1791 }
1792
1793 if (params.verbosity_level() >= 1) {
1794 SOLVER_LOG(&logger, "Termination reason: ",
1795 TerminationReason_Name(result.solve_log.termination_reason()));
1796 SOLVER_LOG(&logger, "Solution point type: ",
1797 PointType_Name(result.solve_log.solution_type()));
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);
1806 const auto& convergence_info = GetConvergenceInformation(
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());
1812 }
1813 }
1814 }
1815 return result;
1816}
1817
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)
1823 : params_(params),
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) {}
1831
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),
1838 };
1839 const QuadraticProgram& qp = WorkingQp();
1840 // This computes the primal portion of the PDHG algorithm:
1841 // argmin_x[gradient(f)(`current_primal_solution_`)^T x + g(x)
1842 // + `current_dual_solution_`^T K x
1843 // + (0.5 / `primal_step_size`) * norm(x - `current_primal_solution_`)^2]
1844 // See Sections 2 - 3 of Chambolle and Pock and the comment in the header.
1845 // We omitted the constant terms from Chambolle and Pock's (7).
1846 // This minimization is easy to do in closed form since it can be separated
1847 // into independent problems for each of the primal variables.
1848 ShardedWorkingQp().PrimalSharder().ParallelForEachShard(
1849 [&](const Sharder::Shard& shard) {
1850 if (!IsLinearProgram(qp)) {
1851 // TODO(user): Does changing this to auto (so it becomes an
1852 // Eigen deferred result), or inlining it below, change performance?
1853 const VectorXd diagonal_scaling =
1854 primal_step_size *
1855 shard(qp.objective_matrix->diagonal()).array() +
1856 1.0;
1857 shard(result.value) =
1858 (shard(current_primal_solution_) -
1859 primal_step_size *
1860 (shard(qp.objective_vector) - shard(current_dual_product_)))
1861 // Scale i-th element by 1 / (1 + `primal_step_size` * Q_{ii})
1862 .cwiseQuotient(diagonal_scaling)
1863 .cwiseMin(shard(qp.variable_upper_bounds))
1864 .cwiseMax(shard(qp.variable_lower_bounds));
1865 } else {
1866 // The formula in the LP case is simplified for better performance.
1867 shard(result.value) =
1868 (shard(current_primal_solution_) -
1869 primal_step_size *
1870 (shard(qp.objective_vector) - shard(current_dual_product_)))
1871 .cwiseMin(shard(qp.variable_upper_bounds))
1872 .cwiseMax(shard(qp.variable_lower_bounds));
1873 }
1874 shard(result.delta) =
1875 shard(result.value) - shard(current_primal_solution_);
1876 });
1877 return result;
1878}
1879
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),
1888 };
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));
1898 });
1899 }
1900 ShardedWorkingQp().TransposedConstraintMatrixSharder().ParallelForEachShard(
1901 [&](const Sharder::Shard& shard) {
1902 VectorXd temp;
1903 if (next_primal_product) {
1904 CHECK(current_primal_product_.has_value());
1905 temp = shard(current_dual_solution_) -
1906 dual_step_size *
1907 (-extrapolation_factor * shard(*current_primal_product_) +
1908 (extrapolation_factor + 1) * shard(*next_primal_product));
1909 } else {
1910 temp = shard(current_dual_solution_) -
1911 dual_step_size *
1912 shard(ShardedWorkingQp().TransposedConstraintMatrix())
1913 .transpose() *
1914 extrapolated_primal.value();
1915 }
1916
1917 // Each element of the argument of `.cwiseMin()` is the critical point
1918 // of the respective 1D minimization problem if it's negative.
1919 // Likewise the argument to the `.cwiseMax()` is the critical point if
1920 // positive.
1921 shard(result.value) =
1922 VectorXd::Zero(temp.size())
1923 .cwiseMin(temp +
1924 dual_step_size * shard(qp.constraint_upper_bounds))
1925 .cwiseMax(temp +
1926 dual_step_size * shard(qp.constraint_lower_bounds));
1927 shard(result.delta) =
1928 (shard(result.value) - shard(current_dual_solution_));
1929 });
1930 return result;
1931}
1932
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())};
1937}
1938
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;
1945}
1946
1947double Solver::ComputeNonlinearity(const VectorXd& delta_primal,
1948 const VectorXd& next_dual_product) const {
1949 // Lemma 1 in Chambolle and Pock includes a term with L_f, the Lipshitz
1950 // constant of f. This is zero in our formulation.
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_));
1956 });
1957}
1958
1959void Solver::SetCurrentPrimalAndDualProducts() {
1960 if (params_.linesearch_rule() ==
1961 PrimalDualHybridGradientParams::MALITSKY_POCK_LINESEARCH_RULE) {
1962 current_primal_product_ = TransposedMatrixVectorProduct(
1963 ShardedWorkingQp().TransposedConstraintMatrix(),
1964 current_primal_solution_,
1965 ShardedWorkingQp().TransposedConstraintMatrixSharder());
1966 } else {
1967 current_primal_product_.reset();
1968 }
1969 current_dual_product_ = TransposedMatrixVectorProduct(
1970 WorkingQp().constraint_matrix, current_dual_solution_,
1971 ShardedWorkingQp().ConstraintMatrixSharder());
1972}
1973
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;
1981 }
1982 stats.set_iteration_number(iterations_completed_);
1983 stats.set_cumulative_rejected_steps(num_rejected_steps_);
1984 // TODO(user): This formula doesn't account for kkt passes in major
1985 // iterations.
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_);
1993 return stats;
1994}
1995
1996double Solver::DistanceTraveledFromLastStart(
1997 const VectorXd& primal_solution, const VectorXd& dual_solution) const {
1998 return std::sqrt((0.5 * primal_weight_) *
1999 SquaredDistance(primal_solution,
2000 last_primal_start_point_,
2001 ShardedWorkingQp().PrimalSharder()) +
2002 (0.5 / primal_weight_) *
2003 SquaredDistance(dual_solution, last_dual_start_point_,
2004 ShardedWorkingQp().DualSharder()));
2005}
2006
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 /*primal_product=*/current_primal_product_.has_value()
2015 ? &current_primal_product_.value()
2016 : nullptr,
2017 &current_dual_product_, params_.use_diagonal_qp_trust_region_solver(),
2018 params_.diagonal_qp_trust_region_solver_tolerance());
2019}
2020
2021LocalizedLagrangianBounds Solver::ComputeLocalizedBoundsAtAverage() const {
2022 // TODO(user): These vectors are recomputed again for termination checks
2023 // and again if we eventually restart to the average.
2024 VectorXd average_primal = PrimalAverage();
2025 VectorXd average_dual = DualAverage();
2026
2027 const double distance_traveled_by_average =
2028 DistanceTraveledFromLastStart(average_primal, average_dual);
2029
2031 ShardedWorkingQp(), average_primal, average_dual,
2032 PrimalDualNorm::kEuclideanNorm, primal_weight_,
2033 distance_traveled_by_average,
2034 /*primal_product=*/nullptr, /*dual_product=*/nullptr,
2035 params_.use_diagonal_qp_trust_region_solver(),
2036 params_.diagonal_qp_trust_region_solver_tolerance());
2037}
2038
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);
2046}
2047
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;
2053}
2054
2055// TODO(user): Review / cleanup adaptive heuristic.
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()) {
2061 return true;
2062 }
2063 if (gap_reduction_ratio < params_.necessary_reduction_for_restart() &&
2064 candidate_normalized_gap > normalized_gap_at_last_trial_) {
2065 // We've made the "necessary" amount of progress, and iterates appear to
2066 // be getting worse, so restart.
2067 return true;
2068 }
2069 return false;
2070}
2071
2072RestartChoice Solver::DetermineDistanceBasedRestartChoice() const {
2073 // The following checks are safeguards that normally should not be triggered.
2074 if (primal_average_.NumTerms() == 0) {
2076 } else if (distance_based_restart_info_.length_of_last_restart_period == 0) {
2078 }
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;
2085
2086 // A restart should be triggered when the normalized distance traveled by
2087 // the average is at least a constant factor smaller than the last.
2088 // TODO(user): Experiment with using `.necessary_reduction_for_restart()`
2089 // as a heuristic when deciding if a restart should be triggered.
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)) {
2094 // Restart at current solution when it yields a smaller normalized potential
2095 // function value than the average (heuristic suggested by ohinder@).
2096 if (AverageHasBetterPotential(ComputeLocalizedBoundsAtAverage(),
2097 ComputeLocalizedBoundsAtCurrent())) {
2099 } else {
2101 }
2102 } else {
2104 }
2105}
2106
2107RestartChoice Solver::ChooseRestartToApply(const bool is_major_iteration) {
2108 if (!primal_average_.HasNonzeroWeight() &&
2109 !dual_average_.HasNonzeroWeight()) {
2111 }
2112 // TODO(user): This forced restart is very important for the performance of
2113 // `ADAPTIVE_HEURISTIC`. Test if the impact comes primarily from the first
2114 // forced restart (which would unseat a good initial starting point that could
2115 // prevent restarts early in the solve) or if it's really needed for the full
2116 // duration of the solve. If it is really needed, should we then trigger major
2117 // iterations on powers of two?
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())) {
2125 } else {
2127 }
2128 }
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;
2141 RestartChoice choice;
2142 if (AverageHasBetterPotential(local_bounds_at_average,
2143 local_bounds_at_current)) {
2144 normalized_gap = NormalizedGap(local_bounds_at_average);
2146 } else {
2147 normalized_gap = NormalizedGap(local_bounds_at_current);
2149 }
2150 if (ShouldDoAdaptiveRestartHeuristic(normalized_gap)) {
2151 return choice;
2152 } else {
2153 normalized_gap_at_last_trial_ = normalized_gap;
2155 }
2156 }
2157 case PrimalDualHybridGradientParams::ADAPTIVE_DISTANCE_BASED: {
2158 return DetermineDistanceBasedRestartChoice();
2159 }
2160 default:
2161 LOG(FATAL) << "Unrecognized restart_strategy "
2162 << params_.restart_strategy();
2164 }
2165 } else {
2167 }
2168}
2169
2170VectorXd Solver::PrimalAverage() const {
2171 if (primal_average_.HasNonzeroWeight()) {
2172 return primal_average_.ComputeAverage();
2173 } else {
2174 return current_primal_solution_;
2175 }
2176}
2177
2178VectorXd Solver::DualAverage() const {
2179 if (dual_average_.HasNonzeroWeight()) {
2180 return dual_average_.ComputeAverage();
2181 } else {
2182 return current_dual_solution_;
2183 }
2184}
2185
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());
2193 // This choice of a nonzero tolerance balances performance and numerical
2194 // issues caused by very huge or very tiny weights. It was picked as the best
2195 // among {0.0, 1.0e-20, 2.0e-16, 1.0e-10, 1.0e-5} on the preprocessed MIPLIB
2196 // dataset. The effect of changing this value is relatively minor overall.
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_;
2201 }
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_);
2210 }
2211 return new_primal_weight;
2212}
2213
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(),
2221 primal_solution);
2222 AssignVector(current_dual_solution_, ShardedWorkingQp().DualSharder(),
2223 dual_solution);
2224 break;
2226 AssignVector(current_primal_delta_, ShardedWorkingQp().PrimalSharder(),
2227 primal_solution);
2228 AssignVector(current_dual_delta_, ShardedWorkingQp().DualSharder(),
2229 dual_solution);
2230 break;
2233 break;
2234 default:
2235 // Default to average whenever `output_type` is `POINT_TYPE_NONE`.
2236 output_type = POINT_TYPE_AVERAGE_ITERATE;
2237 break;
2238 }
2239 return ConstructSolverResult(
2240 std::move(primal_solution), std::move(dual_solution), stats,
2241 termination_reason, output_type, std::move(solve_log));
2242}
2243
2244void Solver::ApplyRestartChoice(const RestartChoice restart_to_apply) {
2245 switch (restart_to_apply) {
2248 return;
2250 if (params_.verbosity_level() >= 4) {
2251 SOLVER_LOG(&preprocess_solver_->Logger(),
2252 "Restarted to current on iteration ", iterations_completed_,
2253 " after ", primal_average_.NumTerms(), " iterations");
2254 }
2255 break;
2257 if (params_.verbosity_level() >= 4) {
2258 SOLVER_LOG(&preprocess_solver_->Logger(),
2259 "Restarted to average on iteration ", iterations_completed_,
2260 " after ", primal_average_.NumTerms(), " iterations");
2261 }
2262 current_primal_solution_ = primal_average_.ComputeAverage();
2263 current_dual_solution_ = dual_average_.ComputeAverage();
2264 SetCurrentPrimalAndDualProducts();
2265 break;
2266 }
2267 primal_weight_ = ComputeNewPrimalWeight();
2268 ratio_last_two_step_sizes_ = 1;
2269 if (params_.restart_strategy() ==
2270 PrimalDualHybridGradientParams::ADAPTIVE_HEURISTIC) {
2271 // It's important for the theory that the distances here are calculated
2272 // given the new primal weight.
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) {
2282 // Update parameters for distance-based restarts.
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()};
2287 }
2288 primal_average_.Clear();
2289 dual_average_.Clear();
2290 AssignVector(current_primal_solution_, ShardedWorkingQp().PrimalSharder(),
2291 /*dest=*/last_primal_start_point_);
2292 AssignVector(current_dual_solution_, ShardedWorkingQp().DualSharder(),
2293 /*dest=*/last_dual_start_point_);
2294}
2295
2296// Adds the work (`iteration_number`, `cumulative_kkt_matrix_passes`,
2297// `cumulative_rejected_steps`, and `cumulative_time_sec`) from
2298// `additional_work_stats` to `stats` and returns the result.
2299// `stats` is intentionally passed by value.
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());
2312 return stats;
2313}
2314
2315// Accumulates and returns the work (`iteration_number`,
2316// `cumulative_kkt_matrix_passes`, `cumulative_rejected_steps`, and
2317// `cumulative_time_sec`) from `solve_log.feasibility_polishing_details`.
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());
2324 }
2325 return result;
2326}
2327
2328bool TerminationReasonIsInterrupted(const TerminationReason reason) {
2330}
2331
2332bool TerminationReasonIsWorkLimitNotInterrupted(
2333 const TerminationReason reason) {
2334 return reason == TERMINATION_REASON_ITERATION_LIMIT ||
2337}
2338
2339// Note: `TERMINATION_REASON_INTERRUPTED_BY_USER` is treated as a work limit
2340// (that was determined in real-time by the user).
2341bool TerminationReasonIsWorkLimit(const TerminationReason reason) {
2342 return TerminationReasonIsWorkLimitNotInterrupted(reason) ||
2343 TerminationReasonIsInterrupted(reason);
2344}
2345
2346bool DoFeasibilityPolishingAfterLimitsReached(
2347 const PrimalDualHybridGradientParams& params,
2348 const TerminationReason reason) {
2349 if (TerminationReasonIsWorkLimitNotInterrupted(reason)) {
2350 return params.apply_feasibility_polishing_after_limits_reached();
2351 }
2352 if (TerminationReasonIsInterrupted(reason)) {
2353 return params.apply_feasibility_polishing_if_solver_is_interrupted();
2354 }
2355 return false;
2356}
2357
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;
2367 // Just decide what to do for now. The actual restart, if any, is
2368 // performed after the termination check.
2369 const RestartChoice restart = force_numerical_termination
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 =
2376 CheckSimpleTerminationCriteria(params_.termination_criteria(),
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;
2381 // We check termination on every major iteration.
2382 DCHECK(!is_major_iteration || check_termination);
2383 if (check_termination) {
2384 // Check for termination and update iteration stats with both simple and
2385 // solution statistics. The later are computationally harder to compute and
2386 // hence only computed here.
2387 VectorXd primal_average = PrimalAverage();
2388 VectorXd dual_average = DualAverage();
2389
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 ? &current_primal_delta_
2398 : nullptr,
2399 current_dual_delta_.size() > 0 ? &current_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;
2404 }
2405 // We've terminated.
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 "
2416 "limits reached";
2417 return *feasibility_result;
2418 }
2419 }
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));
2426 }
2427 } else if (params_.record_iteration_stats()) {
2428 // Record simple iteration stats only.
2429 *solve_log.add_iteration_stats() = stats;
2430 }
2431 ApplyRestartChoice(restart);
2432 return std::nullopt;
2433}
2434
2435void Solver::ResetAverageToCurrent() {
2436 primal_average_.Clear();
2437 dual_average_.Clear();
2438 primal_average_.Add(current_primal_solution_, /*weight=*/1.0);
2439 dual_average_.Add(current_dual_solution_, /*weight=*/1.0);
2440}
2441
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);
2447 SOLVER_LOG(&preprocess_solver_->Logger(),
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_);
2452 }
2453}
2454
2455void Solver::LogInnerIterationLimitHit() const {
2456 SOLVER_LOG(&preprocess_solver_->Logger(),
2457 "WARNING: Inner iteration limit reached at iteration ",
2458 iterations_completed_);
2459}
2460
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);
2466 // The theory by Malitsky and Pock holds for any new_step_size in the interval
2467 // [`step_size`, `step_size` * sqrt(1 + `ratio_last_two_step_sizes_`)].
2468 // `dilating_coeff` determines where in this interval the new step size lands.
2469 // NOTE: Malitsky and Pock use theta for `ratio_last_two_step_sizes`.
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;
2486 });
2487
2488 for (bool accepted_step = false; !accepted_step; ++inner_iterations) {
2489 if (inner_iterations >= 60) {
2490 LogInnerIterationLimitHit();
2491 ResetAverageToCurrent();
2492 outcome = InnerStepOutcome::kForceNumericalTermination;
2493 break;
2494 }
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);
2500
2501 VectorXd next_dual_product = TransposedMatrixVectorProduct(
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) {
2511 // Accept new_step_size as a good step.
2512 step_size_ = new_primal_step_size * primal_weight_;
2513 ratio_last_two_step_sizes_ = new_last_two_step_sizes_ratio;
2514 // Malitsky and Pock guarantee uses a nonsymmetric weighted average,
2515 // the primal variable average involves the initial point, while the dual
2516 // doesn't. See Theorem 2 in https://arxiv.org/pdf/1608.08883.pdf for
2517 // details.
2518 if (!primal_average_.HasNonzeroWeight()) {
2519 primal_average_.Add(
2520 current_primal_solution_,
2521 /*weight=*/new_primal_step_size * new_last_two_step_sizes_ratio);
2522 }
2523
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 /*weight=*/new_primal_step_size);
2530 dual_average_.Add(current_dual_solution_,
2531 /*weight=*/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;
2543 }
2544 current_primal_delta_ = std::move(next_primal_solution.delta);
2545 current_dual_delta_ = std::move(next_dual_solution.delta);
2546 break;
2547 } else {
2548 new_primal_step_size = step_size_downscaling * new_primal_step_size;
2549 }
2550 }
2551 // `inner_iterations` isn't incremented for the accepted step.
2552 num_rejected_steps_ += inner_iterations;
2553 return outcome;
2554}
2555
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;
2564 break;
2565 }
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, /*extrapolation_factor=*/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;
2579 break;
2580 } else if (movement > kDivergentMovement) {
2581 LogNumericalTermination(next_primal_solution.delta,
2582 next_dual_solution.delta);
2583 outcome = InnerStepOutcome::kForceNumericalTermination;
2584 break;
2585 }
2586 VectorXd next_dual_product = TransposedMatrixVectorProduct(
2587 WorkingQp().constraint_matrix, next_dual_solution.value,
2588 ShardedWorkingQp().ConstraintMatrixSharder());
2589 const double nonlinearity =
2590 ComputeNonlinearity(next_primal_solution.delta, next_dual_product);
2591
2592 // See equation (5) in https://arxiv.org/pdf/2106.04756.pdf.
2593 const double step_size_limit =
2594 nonlinearity > 0 ? movement / nonlinearity
2595 : std::numeric_limits<double>::infinity();
2596
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_, /*weight=*/step_size_);
2605 dual_average_.Add(current_dual_solution_, /*weight=*/step_size_);
2606 accepted_step = true;
2607 }
2608 const double total_steps_attempted =
2609 num_rejected_steps_ + inner_iterations + iterations_completed_ + 1;
2610 // Our step sizes are a factor 1 - (`total_steps_attempted` + 1)^(-
2611 // `step_size_reduction_exponent`) smaller than they could be as a margin to
2612 // reduce rejected steps.
2613 // The std::isinf() test protects against NAN if std::pow() == 1.0.
2614 const double first_term =
2615 std::isinf(step_size_limit)
2616 ? step_size_limit
2617 : (1 - std::pow(total_steps_attempted + 1.0,
2618 -params_.adaptive_linesearch_parameters()
2619 .step_size_reduction_exponent())) *
2620 step_size_limit;
2621 const double second_term =
2622 (1 + std::pow(total_steps_attempted + 1.0,
2623 -params_.adaptive_linesearch_parameters()
2624 .step_size_growth_exponent())) *
2625 step_size_;
2626 // From the first term when we have to reject a step, `step_size_`
2627 // decreases by a factor of at least 1 - (`total_steps_attempted` + 1)^(-
2628 // `step_size_reduction_exponent`). From the second term we increase
2629 // `step_size_` by a factor of at most 1 + (`total_steps_attempted` +
2630 // 1)^(-`step_size_growth_exponent`) Therefore if more than order
2631 // (`total_steps_attempted` + 1)^(`step_size_reduction_exponent`
2632 // - `step_size_growth_exponent`) fraction of the time we have a rejected
2633 // step, we overall decrease `step_size_`. When `step_size_` is
2634 // sufficiently small we stop having rejected steps.
2635 step_size_ = std::min(first_term, second_term);
2636 }
2637 // `inner_iterations` is incremented for the accepted step.
2638 num_rejected_steps_ += inner_iterations - 1;
2639 return outcome;
2640}
2641
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, /*extrapolation_factor=*/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;
2660 }
2661 VectorXd next_dual_product = TransposedMatrixVectorProduct(
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_, /*weight=*/step_size_);
2671 dual_average_.Add(current_dual_solution_, /*weight=*/step_size_);
2672 return InnerStepOutcome::kSuccessful;
2673}
2674
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));
2679 return full_stats;
2680}
2681
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()));
2697 return details;
2698}
2699
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 =
2704 EffectiveOptimalityCriteria(params_.termination_criteria());
2705
2706 VectorXd average_primal = PrimalAverage();
2707 VectorXd average_dual = DualAverage();
2708
2709 ConvergenceInformation first_convergence_info;
2710 preprocess_solver_->ComputeConvergenceAndInfeasibilityFromWorkingSolution(
2711 params_, average_primal, average_dual, POINT_TYPE_AVERAGE_ITERATE,
2712 &first_convergence_info, nullptr);
2713
2714 // NOTE: Even though the objective gap sometimes is reduced by feasibility
2715 // polishing, it is usually increased, and an experiment (on MIPLIB2017)
2716 // shows that this test reduces the iteration count by 3-4% on average.
2717 if (!ObjectiveGapMet(optimality_criteria, first_convergence_info)) {
2718 std::optional<TerminationReasonAndPointType> simple_termination_reason =
2719 CheckSimpleTerminationCriteria(params_.termination_criteria(),
2720 TotalWorkSoFar(solve_log),
2721 interrupt_solve);
2722 if (!(simple_termination_reason.has_value() &&
2723 DoFeasibilityPolishingAfterLimitsReached(
2724 params_, simple_termination_reason->reason))) {
2725 if (params_.verbosity_level() >= 2) {
2726 SOLVER_LOG(&preprocess_solver_->Logger(),
2727 "Skipping feasibility polishing because the objective gap "
2728 "is too large.");
2729 }
2730 return std::nullopt;
2731 }
2732 }
2733
2734 if (params_.verbosity_level() >= 2) {
2735 SOLVER_LOG(&preprocess_solver_->Logger(),
2736 "Starting primal feasibility polishing");
2737 }
2738 SolverResult primal_result = TryPrimalPolishing(
2739 std::move(average_primal), iteration_limit, interrupt_solve, solve_log);
2740
2741 if (params_.verbosity_level() >= 2) {
2742 SOLVER_LOG(
2743 &preprocess_solver_->Logger(),
2744 "Primal feasibility polishing termination reason: ",
2745 TerminationReason_Name(primal_result.solve_log.termination_reason()));
2746 }
2747 if (TerminationReasonIsWorkLimit(
2748 primal_result.solve_log.termination_reason())) {
2749 // Have we also reached the overall work limit? If so, consider finishing
2750 // the final polishing phase.
2751 std::optional<TerminationReasonAndPointType> simple_termination_reason =
2752 CheckSimpleTerminationCriteria(params_.termination_criteria(),
2753 TotalWorkSoFar(solve_log),
2754 interrupt_solve);
2755 if (!(simple_termination_reason.has_value() &&
2756 DoFeasibilityPolishingAfterLimitsReached(
2757 params_, simple_termination_reason->reason))) {
2758 return std::nullopt;
2759 }
2760 } else if (primal_result.solve_log.termination_reason() !=
2761 TERMINATION_REASON_OPTIMAL) {
2762 // Note: `TERMINATION_REASON_PRIMAL_INFEASIBLE` could happen normally, but
2763 // we haven't ensured that the correct solution is returned in that case, so
2764 // we ignore the polishing result indicating infeasibility.
2765 // `TERMINATION_REASON_NUMERICAL_ERROR` can occur, but would be surprising
2766 // and interesting. Other termination reasons are probably bugs.
2767 SOLVER_LOG(&preprocess_solver_->Logger(),
2768 "WARNING: Primal feasibility polishing terminated with error ",
2769 primal_result.solve_log.termination_reason());
2770 return std::nullopt;
2771 }
2772
2773 if (params_.verbosity_level() >= 2) {
2774 SOLVER_LOG(&preprocess_solver_->Logger(),
2775 "Starting dual feasibility polishing");
2776 }
2777 SolverResult dual_result = TryDualPolishing(
2778 std::move(average_dual), iteration_limit, interrupt_solve, solve_log);
2779
2780 if (params_.verbosity_level() >= 2) {
2781 SOLVER_LOG(
2782 &preprocess_solver_->Logger(),
2783 "Dual feasibility polishing termination reason: ",
2784 TerminationReason_Name(dual_result.solve_log.termination_reason()));
2785 }
2786
2787 IterationStats full_stats = TotalWorkSoFar(solve_log);
2788 std::optional<TerminationReasonAndPointType> simple_termination_reason =
2789 CheckSimpleTerminationCriteria(params_.termination_criteria(), full_stats,
2790 interrupt_solve);
2791 if (TerminationReasonIsWorkLimit(
2792 dual_result.solve_log.termination_reason())) {
2793 // Have we also reached the overall work limit? If so, consider falling out
2794 // of the "if" test and returning the polishing solution anyway.
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);
2807 } else {
2808 return std::nullopt;
2809 }
2810 } else if (dual_result.solve_log.termination_reason() !=
2811 TERMINATION_REASON_OPTIMAL) {
2812 // Note: The comment in the corresponding location when checking the
2813 // termination reason for primal feasibility polishing applies here
2814 // too.
2815 SOLVER_LOG(&preprocess_solver_->Logger(),
2816 "WARNING: Dual feasibility polishing terminated with error ",
2817 dual_result.solve_log.termination_reason());
2818 return std::nullopt;
2819 }
2820
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) {
2826 SOLVER_LOG(&preprocess_solver_->Logger(),
2827 "solution stats for polished solution:");
2828 LogIterationStatsHeader(params_.verbosity_level(),
2829 /*use_feasibility_polishing=*/true,
2830 preprocess_solver_->Logger());
2831 LogIterationStats(params_.verbosity_level(),
2832 /*use_feasibility_polishing=*/true,
2833 IterationType::kFeasibilityPolishingTermination,
2834 full_stats, params_.termination_criteria(),
2835 preprocess_solver_->OriginalBoundNorms(),
2836 POINT_TYPE_FEASIBILITY_POLISHING_SOLUTION,
2837 preprocess_solver_->Logger());
2838 }
2839 std::optional<TerminationReasonAndPointType> earned_termination =
2840 CheckIterateTerminationCriteria(params_.termination_criteria(),
2841 full_stats,
2842 preprocess_solver_->OriginalBoundNorms(),
2843 /*force_numerical_termination=*/false);
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);
2854 }
2855 // Note: A typical termination check would now call
2856 // `CheckSimpleTerminationCriteria`. However, there is no obvious iterate to
2857 // use for a `SolverResult`. Instead, allow the main iteration loop to notice
2858 // if the termination criteria checking has caused us to meet the simple
2859 // termination criteria. This would be a rare case anyway, since dual
2860 // feasibility polishing didn't terminate because of a work limit, and
2861 // `full_stats` is created shortly after `TryDualPolishing` returned.
2862 return std::nullopt;
2863}
2864
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());
2873 } else {
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()));
2882 }
2883 return criteria;
2884}
2885
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;
2897 }
2898
2899 // This will save the original objective after the swap.
2900 VectorXd objective;
2901 SetZero(ShardedWorkingQp().PrimalSharder(), objective);
2902 preprocess_solver_->SwapObjectiveVector(objective);
2903
2904 TerminationCriteria::DetailedOptimalityCriteria criteria =
2905 EffectiveOptimalityCriteria(params_.termination_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;
2913 // TODO(user): Evaluate disabling primal weight updates for this primal
2914 // feasibility problem.
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;
2922 // Stop `timer_`, since the time in `primal_solver.Solve()` will be recorded
2923 // by its timer.
2924 timer_.Stop();
2925 SolverResult primal_result = primal_solver.Solve(
2926 IterationType::kPrimalFeasibility, interrupt_solve, primal_solve_log);
2927 timer_.Start();
2928 // Restore the original objective.
2929 preprocess_solver_->SwapObjectiveVector(objective);
2930
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;
2936}
2937
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;
2942 };
2943 sharder.ParallelForEachShard([&](const Sharder::Shard& shard) {
2944 shard(output) = shard(input).unaryExpr(make_finite_values_zero);
2945 });
2946 return output;
2947}
2948
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;
2961 }
2962
2963 // These will initially contain the homogenous variable and constraint
2964 // bounds, but will contain the original variable and constraint bounds
2965 // after the swap.
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);
2978
2979 TerminationCriteria::DetailedOptimalityCriteria criteria =
2980 EffectiveOptimalityCriteria(params_.termination_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;
2988 // TODO(user): Evaluate disabling primal weight updates for this dual
2989 // feasibility problem.
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;
2997 // Stop `timer_`, since the time in `dual_solver.Solve()` will be recorded
2998 // by its timer.
2999 timer_.Stop();
3000 SolverResult dual_result = dual_solver.Solve(IterationType::kDualFeasibility,
3001 interrupt_solve, dual_solve_log);
3002 timer_.Start();
3003 // Restore the original bounds.
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);
3012 return dual_result;
3013}
3014
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();
3019 timer_.Start();
3020 last_primal_start_point_ =
3021 CloneVector(current_primal_solution_, ShardedWorkingQp().PrimalSharder());
3022 last_dual_start_point_ =
3023 CloneVector(current_dual_solution_, ShardedWorkingQp().DualSharder());
3024 // Note: Any cached values computed here also need to be recomputed after a
3025 // restart.
3026
3027 ratio_last_two_step_sizes_ = 1;
3028 SetCurrentPrimalAndDualProducts();
3029
3030 // This is set to true if we can't proceed any more because of numerical
3031 // issues. We may or may not have found the optimal solution.
3032 bool force_numerical_termination = false;
3033
3034 int next_feasibility_polishing_iteration = 100;
3035
3036 num_rejected_steps_ = 0;
3037
3038 IterationStats work_from_feasibility_polishing =
3039 WorkFromFeasibilityPolishing(solve_log);
3040 for (iterations_completed_ = 0;; ++iterations_completed_) {
3041 // This code performs the logic of the major iterations and termination
3042 // checks. It may modify the current solution and primal weight (e.g., when
3043 // performing a restart).
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();
3050 }
3051
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;
3061 }
3062 next_feasibility_polishing_iteration *= 2;
3063 // Update work to include new feasibility phases.
3064 work_from_feasibility_polishing = WorkFromFeasibilityPolishing(solve_log);
3065 }
3066
3067 // TODO(user): If we use a step rule that could reject many steps in a
3068 // row, we should add a termination check within this loop also. For the
3069 // Malitsky and Pock rule, we perform a termination check and declare
3070 // NUMERICAL_ERROR whenever we hit 60 inner iterations.
3071 InnerStepOutcome outcome;
3072 switch (params_.linesearch_rule()) {
3073 case PrimalDualHybridGradientParams::MALITSKY_POCK_LINESEARCH_RULE:
3074 outcome = TakeMalitskyPockStep();
3075 break;
3076 case PrimalDualHybridGradientParams::ADAPTIVE_LINESEARCH_RULE:
3077 outcome = TakeAdaptiveStep();
3078 break;
3079 case PrimalDualHybridGradientParams::CONSTANT_STEP_SIZE_RULE:
3080 outcome = TakeConstantSizeStep();
3081 break;
3082 default:
3083 LOG(FATAL) << "Unrecognized linesearch rule "
3084 << params_.linesearch_rule();
3085 }
3086 if (outcome == InnerStepOutcome::kForceNumericalTermination) {
3087 force_numerical_termination = true;
3088 }
3089 } // loop over iterations
3090}
3091
3092} // namespace
3093
3096 const std::atomic<bool>* interrupt_solve,
3097 std::function<void(const std::string&)> message_callback,
3098 IterationStatsCallback iteration_stats_callback) {
3099 return PrimalDualHybridGradient(std::move(qp), params, std::nullopt,
3100 interrupt_solve, std::move(message_callback),
3101 std::move(iteration_stats_callback));
3102}
3103
3104// See comment at top of file.
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) {
3111 SolverLogger logger;
3112 logger.EnableLogging(true);
3113 if (message_callback) {
3114 logger.AddInfoLoggingCallback(std::move(message_callback));
3115 } else {
3116 logger.SetLogToStdOut(true);
3117 }
3118 const absl::Status params_status =
3120 if (!params_status.ok()) {
3121 return ErrorSolverResult(TERMINATION_REASON_INVALID_PARAMETER,
3122 params_status.ToString(), logger);
3123 }
3124 if (!qp.constraint_matrix.isCompressed()) {
3125 return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
3126 "constraint_matrix must be in compressed format. "
3127 "Call constraint_matrix.makeCompressed()",
3128 logger);
3129 }
3130 const absl::Status dimensions_status = ValidateQuadraticProgramDimensions(qp);
3131 if (!dimensions_status.ok()) {
3132 return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
3133 dimensions_status.ToString(), logger);
3134 }
3135 if (qp.objective_scaling_factor == 0) {
3136 return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
3137 "The objective scaling factor cannot be zero.",
3138 logger);
3139 }
3140 if (params.use_feasibility_polishing() && !IsLinearProgram(qp)) {
3141 return ErrorSolverResult(
3143 "use_feasibility_polishing is only implemented for linear programs.",
3144 logger);
3145 }
3146 PreprocessSolver solver(std::move(qp), params, &logger);
3147 return solver.PreprocessAndSolve(params, std::move(initial_solution),
3148 interrupt_solve,
3149 std::move(iteration_stats_callback));
3150}
3151
3152namespace internal {
3153
3156 glop::ProblemSolution glop_solution(
3157 glop::RowIndex(solution.dual_solution.size()),
3158 glop::ColIndex(solution.primal_solution.size()));
3159 // This doesn't matter much as glop's preprocessor doesn't use this much.
3160 // We pick IMPRECISE since we are often calling this code early in the solve.
3161 glop_solution.status = glop::ProblemStatus::IMPRECISE;
3162 for (glop::RowIndex i{0}; i.value() < solution.dual_solution.size(); ++i) {
3163 if (qp.constraint_lower_bounds[i.value()] ==
3164 qp.constraint_upper_bounds[i.value()]) {
3165 glop_solution.constraint_statuses[i] =
3167 } else if (solution.dual_solution[i.value()] > 0) {
3168 glop_solution.constraint_statuses[i] =
3170 } else if (solution.dual_solution[i.value()] < 0) {
3171 glop_solution.constraint_statuses[i] =
3173 } else {
3175 }
3176 }
3177
3178 for (glop::ColIndex i{0}; i.value() < solution.primal_solution.size(); ++i) {
3179 const bool at_lb = solution.primal_solution[i.value()] <=
3180 qp.variable_lower_bounds[i.value()];
3181 const bool at_ub = solution.primal_solution[i.value()] >=
3182 qp.variable_upper_bounds[i.value()];
3183 // Note that `ShardedWeightedAverage` is designed so that variables at their
3184 // bounds will be exactly at their bounds even with floating-point roundoff.
3185 if (at_lb) {
3186 if (at_ub) {
3188 } else {
3189 glop_solution.variable_statuses[i] =
3191 }
3192 } else {
3193 if (at_ub) {
3194 glop_solution.variable_statuses[i] =
3196 } else {
3198 }
3199 }
3200 }
3201 return glop_solution;
3202}
3203
3204} // namespace internal
3205
3206} // namespace operations_research::pdlp
double Get() const
Definition timer.h:44
void Start()
When Start() is called multiple times, only the most recent is used.
Definition timer.h:30
void SetLogToStdOut(bool enable)
Should all messages be displayed on stdout ?
Definition logging.h:53
void AddInfoLoggingCallback(std::function< void(const std::string &message)> callback)
Definition logging.cc:29
Fractional SquaredNorm(const SparseColumn &v)
Definition lp_utils.cc:36
constexpr Fractional kInfinity
Infinity for type Fractional.
Definition lp_types.h:87
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)
Definition solve.cc:62
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)
Definition sharder.cc:175
InfeasibilityInformation ComputeInfeasibilityInformation(const PrimalDualHybridGradientParams &params, 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 &params)
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)
Definition sharder.cc:254
double LInfNorm(const VectorXd &vector, const Sharder &sharder)
Definition sharder.cc:232
TerminationCriteria::DetailedOptimalityCriteria EffectiveOptimalityCriteria(const TerminationCriteria &termination_criteria)
Computes the effective optimality criteria for a TerminationCriteria.
double Distance(const VectorXd &vector1, const VectorXd &vector2, const Sharder &sharder)
Definition sharder.cc:261
VectorXd ReducedCosts(const PrimalDualHybridGradientParams &params, 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)
Definition sharder.cc:160
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)
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)
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)
Definition sharder.cc:213
void CoefficientWiseQuotientInPlace(const VectorXd &scale, const Sharder &sharder, VectorXd &dest)
Definition sharder.cc:220
VectorXd ZeroVector(const Sharder &sharder)
Like VectorXd::Zero(sharder.NumElements()).
Definition sharder.cc:181
VectorXd CloneVector(const VectorXd &vec, const Sharder &sharder)
Definition sharder.cc:207
SolverResult PrimalDualHybridGradient(QuadraticProgram qp, const PrimalDualHybridGradientParams &params, 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()).
Definition sharder.cc:187
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 &params, 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)
Definition sharder.cc:250
RelativeConvergenceInformation ComputeRelativeResiduals(const TerminationCriteria::DetailedOptimalityCriteria &optimality_criteria, const ConvergenceInformation &stats, const QuadraticProgramBoundNorms &bound_norms)
void AssignVector(const VectorXd &vec, const Sharder &sharder, VectorXd &dest)
Definition sharder.cc:201
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
STL namespace.
static int input(yyscan_t yyscanner)
Contains the solution of a LinearProgram as returned by a preprocessor.
Definition lp_data.h:670
ConstraintStatusColumn constraint_statuses
Definition lp_data.h:699
ProblemStatus status
The solution status.
Definition lp_data.h:678
Eigen::SparseMatrix< double, Eigen::ColMajor, int64_t > constraint_matrix
#define SOLVER_LOG(logger,...)
Definition logging.h:110