From e76ecdf237c29031becf075f4bf0b92254ce53ca Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Tue, 17 Jun 2025 12:48:42 +0200 Subject: [PATCH] minor cleaning --- ortools/pdlp/BUILD.bazel | 1 + ortools/pdlp/primal_dual_hybrid_gradient.cc | 97 +++++++++++++++------ 2 files changed, 71 insertions(+), 27 deletions(-) diff --git a/ortools/pdlp/BUILD.bazel b/ortools/pdlp/BUILD.bazel index 7a19ebc39a..0059d6d122 100644 --- a/ortools/pdlp/BUILD.bazel +++ b/ortools/pdlp/BUILD.bazel @@ -144,6 +144,7 @@ cc_library( "//ortools/lp_data:proto_utils", "//ortools/util:logging", "@abseil-cpp//absl/algorithm:container", + "@abseil-cpp//absl/base:nullability", "@abseil-cpp//absl/status", "@abseil-cpp//absl/status:statusor", "@abseil-cpp//absl/strings", diff --git a/ortools/pdlp/primal_dual_hybrid_gradient.cc b/ortools/pdlp/primal_dual_hybrid_gradient.cc index d166ef7700..b86f3e9f1c 100644 --- a/ortools/pdlp/primal_dual_hybrid_gradient.cc +++ b/ortools/pdlp/primal_dual_hybrid_gradient.cc @@ -53,6 +53,7 @@ #include "Eigen/Core" #include "Eigen/SparseCore" #include "absl/algorithm/container.h" +#include "absl/base/nullability.h" #include "absl/log/check.h" #include "absl/log/log.h" #include "absl/status/status.h" @@ -619,7 +620,8 @@ class Solver { NextSolutionAndDelta ComputeNextDualSolution( double dual_step_size, double extrapolation_factor, - const NextSolutionAndDelta& next_primal) const; + const NextSolutionAndDelta& next_primal_solution, + const VectorXd* absl_nullable next_primal_product = nullptr) const; std::pair ComputeMovementTerms( const VectorXd& delta_primal, const VectorXd& delta_dual) const; @@ -630,6 +632,10 @@ class Solver { double ComputeNonlinearity(const VectorXd& delta_primal, const VectorXd& next_dual_product) const; + // Sets current_primal_product_ and current_dual_product_ based on + // current_primal_solution_ and current_dual_solution_ respectively. + void SetCurrentPrimalAndDualProducts(); + // Creates all the simple-to-compute statistics in stats. IterationStats CreateSimpleIterationStats(RestartChoice restart_used) const; @@ -734,6 +740,9 @@ class Solver { WallTimer timer_; int iterations_completed_; int num_rejected_steps_; + // A cache of `constraint_matrix * current_primal_solution_`. + // Malitsky-Pock linesearch only. + std::optional current_primal_product_; // A cache of `constraint_matrix.transpose() * current_dual_solution_`. VectorXd current_dual_product_; // The primal point at which the algorithm was last restarted from, or @@ -1870,31 +1879,41 @@ Solver::NextSolutionAndDelta Solver::ComputeNextPrimalSolution( Solver::NextSolutionAndDelta Solver::ComputeNextDualSolution( double dual_step_size, double extrapolation_factor, - const NextSolutionAndDelta& next_primal_solution) const { + const NextSolutionAndDelta& next_primal_solution, + const VectorXd* absl_nullable next_primal_product) const { const int64_t dual_size = ShardedWorkingQp().DualSize(); NextSolutionAndDelta result = { .value = VectorXd(dual_size), .delta = VectorXd(dual_size), }; const QuadraticProgram& qp = WorkingQp(); - VectorXd extrapolated_primal(ShardedWorkingQp().PrimalSize()); - ShardedWorkingQp().PrimalSharder().ParallelForEachShard( - [&](const Sharder::Shard& shard) { - shard(extrapolated_primal) = - (shard(next_primal_solution.value) + - extrapolation_factor * shard(next_primal_solution.delta)); - }); - // TODO(user): Refactor this multiplication so that we only do one matrix - // vector multiply for the primal variable. This only applies to Malitsky and - // Pock and not to the adaptive step size rule. + std::optional extrapolated_primal; + if (!next_primal_product) { + extrapolated_primal.emplace(ShardedWorkingQp().PrimalSize()); + ShardedWorkingQp().PrimalSharder().ParallelForEachShard( + [&](const Sharder::Shard& shard) { + shard(*extrapolated_primal) = + (shard(next_primal_solution.value) + + extrapolation_factor * shard(next_primal_solution.delta)); + }); + } ShardedWorkingQp().TransposedConstraintMatrixSharder().ParallelForEachShard( [&](const Sharder::Shard& shard) { - VectorXd temp = - shard(current_dual_solution_) - - dual_step_size * - shard(ShardedWorkingQp().TransposedConstraintMatrix()) - .transpose() * - extrapolated_primal; + VectorXd temp; + if (next_primal_product) { + CHECK(current_primal_product_.has_value()); + temp = shard(current_dual_solution_) - + dual_step_size * + (-extrapolation_factor * shard(*current_primal_product_) + + (extrapolation_factor + 1) * shard(*next_primal_product)); + } else { + temp = shard(current_dual_solution_) - + dual_step_size * + shard(ShardedWorkingQp().TransposedConstraintMatrix()) + .transpose() * + extrapolated_primal.value(); + } + // Each element of the argument of `.cwiseMin()` is the critical point // of the respective 1D minimization problem if it's negative. // Likewise the argument to the `.cwiseMax()` is the critical point if @@ -1937,6 +1956,21 @@ double Solver::ComputeNonlinearity(const VectorXd& delta_primal, }); } +void Solver::SetCurrentPrimalAndDualProducts() { + if (params_.linesearch_rule() == + PrimalDualHybridGradientParams::MALITSKY_POCK_LINESEARCH_RULE) { + current_primal_product_ = TransposedMatrixVectorProduct( + ShardedWorkingQp().TransposedConstraintMatrix(), + current_primal_solution_, + ShardedWorkingQp().TransposedConstraintMatrixSharder()); + } else { + current_primal_product_.reset(); + } + current_dual_product_ = TransposedMatrixVectorProduct( + WorkingQp().constraint_matrix, current_dual_solution_, + ShardedWorkingQp().ConstraintMatrixSharder()); +} + IterationStats Solver::CreateSimpleIterationStats( RestartChoice restart_used) const { IterationStats stats; @@ -1977,8 +2011,10 @@ LocalizedLagrangianBounds Solver::ComputeLocalizedBoundsAtCurrent() const { ShardedWorkingQp(), current_primal_solution_, current_dual_solution_, PrimalDualNorm::kEuclideanNorm, primal_weight_, distance_traveled_by_current, - /*primal_product=*/nullptr, ¤t_dual_product_, - params_.use_diagonal_qp_trust_region_solver(), + /*primal_product=*/current_primal_product_.has_value() + ? ¤t_primal_product_.value() + : nullptr, + ¤t_dual_product_, params_.use_diagonal_qp_trust_region_solver(), params_.diagonal_qp_trust_region_solver_tolerance()); } @@ -2225,9 +2261,7 @@ void Solver::ApplyRestartChoice(const RestartChoice restart_to_apply) { } current_primal_solution_ = primal_average_.ComputeAverage(); current_dual_solution_ = dual_average_.ComputeAverage(); - current_dual_product_ = TransposedMatrixVectorProduct( - WorkingQp().constraint_matrix, current_dual_solution_, - ShardedWorkingQp().ConstraintMatrixSharder()); + SetCurrentPrimalAndDualProducts(); break; } primal_weight_ = ComputeNewPrimalWeight(); @@ -2443,6 +2477,14 @@ InnerStepOutcome Solver::TakeMalitskyPockStep() { params_.malitsky_pock_parameters().linesearch_contraction_factor(); const double dual_weight = primal_weight_ * primal_weight_; int inner_iterations = 0; + VectorXd next_primal_product(current_dual_solution_.size()); + ShardedWorkingQp().TransposedConstraintMatrixSharder().ParallelForEachShard( + [&](const Sharder::Shard& shard) { + shard(next_primal_product) = + shard(ShardedWorkingQp().TransposedConstraintMatrix()).transpose() * + next_primal_solution.value; + }); + for (bool accepted_step = false; !accepted_step; ++inner_iterations) { if (inner_iterations >= 60) { LogInnerIterationLimitHit(); @@ -2454,7 +2496,7 @@ InnerStepOutcome Solver::TakeMalitskyPockStep() { new_primal_step_size / primal_step_size; NextSolutionAndDelta next_dual_solution = ComputeNextDualSolution( dual_weight * new_primal_step_size, new_last_two_step_sizes_ratio, - next_primal_solution); + next_primal_solution, &next_primal_product); VectorXd next_dual_product = TransposedMatrixVectorProduct( WorkingQp().constraint_matrix, next_dual_solution.value, @@ -2482,6 +2524,7 @@ InnerStepOutcome Solver::TakeMalitskyPockStep() { current_primal_solution_ = std::move(next_primal_solution.value); current_dual_solution_ = std::move(next_dual_solution.value); current_dual_product_ = std::move(next_dual_product); + current_primal_product_ = std::move(next_primal_product); primal_average_.Add(current_primal_solution_, /*weight=*/new_primal_step_size); dual_average_.Add(current_dual_solution_, @@ -2555,6 +2598,7 @@ InnerStepOutcome Solver::TakeAdaptiveStep() { current_primal_solution_ = std::move(next_primal_solution.value); current_dual_solution_ = std::move(next_dual_solution.value); current_dual_product_ = std::move(next_dual_product); + current_primal_product_.reset(); current_primal_delta_ = std::move(next_primal_solution.delta); current_dual_delta_ = std::move(next_dual_solution.delta); primal_average_.Add(current_primal_solution_, /*weight=*/step_size_); @@ -2620,6 +2664,7 @@ InnerStepOutcome Solver::TakeConstantSizeStep() { current_primal_solution_ = std::move(next_primal_solution.value); current_dual_solution_ = std::move(next_dual_solution.value); current_dual_product_ = std::move(next_dual_product); + current_primal_product_.reset(); current_primal_delta_ = std::move(next_primal_solution.delta); current_dual_delta_ = std::move(next_dual_solution.delta); primal_average_.Add(current_primal_solution_, /*weight=*/step_size_); @@ -2980,9 +3025,7 @@ SolverResult Solver::Solve(const IterationType iteration_type, // restart. ratio_last_two_step_sizes_ = 1; - current_dual_product_ = TransposedMatrixVectorProduct( - WorkingQp().constraint_matrix, current_dual_solution_, - ShardedWorkingQp().ConstraintMatrixSharder()); + SetCurrentPrimalAndDualProducts(); // This is set to true if we can't proceed any more because of numerical // issues. We may or may not have found the optimal solution.