pdlp uses SolverLogger

This commit is contained in:
Laurent Perron
2024-02-15 08:46:25 +01:00
parent 1cbaea9a7b
commit ed9b82dc91
13 changed files with 524 additions and 279 deletions

View File

@@ -117,6 +117,7 @@ cc_library(
"//ortools/lp_data",
"//ortools/lp_data:base",
"//ortools/lp_data:proto_utils",
"//ortools/util:logging",
"@com_google_absl//absl/algorithm:container",
"@com_google_absl//absl/status",
"@com_google_absl//absl/status:statusor",
@@ -252,6 +253,7 @@ cc_library(
":sharder",
"//ortools/base",
"//ortools/base:threadpool",
"//ortools/util:logging",
"@com_google_absl//absl/memory",
"@com_google_absl//absl/strings",
"@eigen//:eigen3",

View File

@@ -60,8 +60,7 @@ struct ResidualNorms {
// any corrections (so the returned `.objective_correction == 0` and
// `.objective_full_correction == 0`). `sharded_qp` is assumed to be the scaled
// problem. If `use_homogeneous_constraint_bounds` is set to true the residuals
// are computed with upper and lower bounds zeroed out (note that we only zero
// out the bounds that are finite in the original problem).
// are computed with all finite bounds mapped to zero.
// NOTE: `componentwise_residual_offset` only affects the value of
// `l_inf_componentwise_residual` in the returned `ResidualNorms`.
ResidualNorms PrimalResidualNorms(
@@ -221,13 +220,13 @@ ResidualNorms DualResidualNorms(const PrimalDualHybridGradientParams& params,
// unscaled_primal_gradient = primal_gradient / scale, so the scales
// cancel out.
if (primal_gradient_shard[i] == 0.0) continue;
const double bound_for_rc = primal_gradient_shard[i] > 0.0
? lower_bound_shard[i]
: upper_bound_shard[i];
const double upper_bound = upper_bound_shard[i];
const double lower_bound = lower_bound_shard[i];
const double bound_for_rc =
primal_gradient_shard[i] > 0.0 ? lower_bound : upper_bound;
dual_full_correction += bound_for_rc * primal_gradient_shard[i];
VariableBounds effective_bounds = EffectiveVariableBounds(
params, primal_solution_shard[i], lower_bound_shard[i],
upper_bound_shard[i]);
params, primal_solution_shard[i], lower_bound, upper_bound);
// The dual correction (using the appropriate bound) is applied even
// if the gradient is handled as a residual, so that the dual
// objective is convex.
@@ -453,13 +452,46 @@ ConvergenceInformation ComputeConvergenceInformation(
return result;
}
namespace {
double PrimalRayMaxSignViolation(const ShardedQuadraticProgram& sharded_qp,
const VectorXd& col_scaling_vec,
const VectorXd& scaled_primal_ray) {
VectorXd primal_ray_local_max_sign_violation(
sharded_qp.PrimalSharder().NumShards());
sharded_qp.PrimalSharder().ParallelForEachShard(
[&](const Sharder::Shard& shard) {
const auto lower_bound_shard =
shard(sharded_qp.Qp().variable_lower_bounds);
const auto upper_bound_shard =
shard(sharded_qp.Qp().variable_upper_bounds);
const auto ray_shard = shard(scaled_primal_ray);
const auto scale_shard = shard(col_scaling_vec);
double local_max = 0.0;
for (int64_t i = 0; i < ray_shard.size(); ++i) {
if (std::isfinite(lower_bound_shard[i])) {
local_max = std::max(local_max, -ray_shard[i] * scale_shard[i]);
}
if (std::isfinite(upper_bound_shard[i])) {
local_max = std::max(local_max, ray_shard[i] * scale_shard[i]);
}
}
primal_ray_local_max_sign_violation[shard.Index()] = local_max;
});
return primal_ray_local_max_sign_violation.lpNorm<Eigen::Infinity>();
}
} // namespace
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, PointType candidate_type) {
const Eigen::VectorXd& scaled_dual_ray,
const Eigen::VectorXd& primal_solution_for_residual_tests,
PointType candidate_type) {
const QuadraticProgram& qp = scaled_sharded_qp.Qp();
CHECK_EQ(col_scaling_vec.size(), scaled_sharded_qp.PrimalSize());
CHECK_EQ(row_scaling_vec.size(), scaled_sharded_qp.DualSize());
@@ -479,8 +511,9 @@ InfeasibilityInformation ComputeInfeasibilityInformation(
// We don't use `dual_residuals.l_inf_componentwise_residual`, so don't need
// to set `componentwise_residual_offset` to a meaningful value.
ResidualNorms dual_residuals = DualResidualNorms(
params, scaled_sharded_qp, col_scaling_vec, scaled_primal_ray,
scaled_primal_gradient, /*componentwise_residual_offset=*/0.0);
params, scaled_sharded_qp, col_scaling_vec,
primal_solution_for_residual_tests, scaled_primal_gradient,
/*componentwise_residual_offset=*/0.0);
double dual_ray_objective =
DualObjectiveBoundsTerm(scaled_sharded_qp, scaled_dual_ray) +
@@ -501,32 +534,12 @@ InfeasibilityInformation ComputeInfeasibilityInformation(
PrimalResidualNorms(scaled_sharded_qp, row_scaling_vec, scaled_primal_ray,
/*componentwise_residual_offset=*/0.0,
/*use_homogeneous_constraint_bounds=*/true);
// `primal_residuals` contains the violations of the linear constraints. The
// signs of the components are also constrained by the presence or absence
// of variable bounds.
VectorXd primal_ray_local_sign_max_violation(
scaled_sharded_qp.PrimalSharder().NumShards());
scaled_sharded_qp.PrimalSharder().ParallelForEachShard(
[&](const Sharder::Shard& shard) {
const auto lower_bound_shard =
shard(scaled_sharded_qp.Qp().variable_lower_bounds);
const auto upper_bound_shard =
shard(scaled_sharded_qp.Qp().variable_upper_bounds);
const auto ray_shard = shard(scaled_primal_ray);
const auto scale_shard = shard(col_scaling_vec);
double local_max = 0.0;
for (int64_t i = 0; i < ray_shard.size(); ++i) {
if (std::isfinite(lower_bound_shard[i])) {
local_max = std::max(local_max, -ray_shard[i] * scale_shard[i]);
}
if (std::isfinite(upper_bound_shard[i])) {
local_max = std::max(local_max, ray_shard[i] * scale_shard[i]);
}
}
primal_ray_local_sign_max_violation[shard.Index()] = local_max;
});
const double primal_ray_sign_max_violation =
primal_ray_local_sign_max_violation.lpNorm<Eigen::Infinity>();
// The primal ray should have been projected onto the feasibility bounds, so
// that it has the correct signs.
DCHECK_EQ(PrimalRayMaxSignViolation(scaled_sharded_qp, col_scaling_vec,
scaled_primal_ray),
0.0);
if (l_inf_primal > 0.0) {
VectorXd scaled_objective_product =
@@ -534,10 +547,8 @@ InfeasibilityInformation ComputeInfeasibilityInformation(
result.set_primal_ray_quadratic_norm(
LInfNorm(scaled_objective_product, scaled_sharded_qp.PrimalSharder()) /
l_inf_primal);
result.set_max_primal_ray_infeasibility(
std::max(primal_residuals.l_inf_residual,
primal_ray_sign_max_violation) /
l_inf_primal);
result.set_max_primal_ray_infeasibility(primal_residuals.l_inf_residual /
l_inf_primal);
result.set_primal_ray_linear_objective(
Dot(scaled_primal_ray, qp.objective_vector,
scaled_sharded_qp.PrimalSharder()) /

View File

@@ -65,13 +65,18 @@ ConvergenceInformation ComputeConvergenceInformation(
// row_scaling_vec)`. `scaled_primal_ray` and `scaled_dual_ray` are
// potential certificates for the scaled problem. The stats are computed with
// respect to the implicit original problem.
// `primal_solution_for_residual_tests` is used instead of `scaled_primal_ray`
// when deciding whether or not to treat a primal gradient as a dual residual
// or not.
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, PointType candidate_type);
const Eigen::VectorXd& scaled_dual_ray,
const Eigen::VectorXd& primal_solution_for_residual_tests,
PointType candidate_type);
// Computes the reduced costs vector, objective_matrix * `primal_solution` +
// objective_vector - constraint_matrix * `dual_solution`, when

View File

@@ -14,6 +14,7 @@
#include "ortools/pdlp/iteration_stats.h"
#include <cmath>
#include <limits>
#include <optional>
#include <utility>

File diff suppressed because it is too large Load Diff

View File

@@ -17,6 +17,7 @@
#include <atomic>
#include <functional>
#include <optional>
#include <string>
#include "Eigen/Core"
#include "ortools/lp_data/lp_data.h"
@@ -128,6 +129,12 @@ struct IterationCallbackInfo {
// if `interrupt_solve->load()` is true, in which case the solve will terminate
// with `TERMINATION_REASON_INTERRUPTED_BY_USER`.
//
// If `message_callback` is not nullptr, solver logging will be passed to
// `message_callback`. Otherwise solver logging will be written to stdout.
// In either case, the amount of logging is controlled by `verbosity_level`.
// In particular, if `verbosity_level == 0`, there will be no logging in either
// case.
//
// If `iteration_stats_callback` is not nullptr, then at each termination step
// (when iteration stats are logged), `iteration_stats_callback` will also
// be called with those iteration stats.
@@ -144,6 +151,7 @@ struct IterationCallbackInfo {
SolverResult PrimalDualHybridGradient(
QuadraticProgram qp, const PrimalDualHybridGradientParams& params,
const std::atomic<bool>* interrupt_solve = nullptr,
std::function<void(const std::string&)> message_callback = nullptr,
std::function<void(const IterationCallbackInfo&)> iteration_stats_callback =
nullptr);
@@ -156,6 +164,7 @@ SolverResult PrimalDualHybridGradient(
QuadraticProgram qp, const PrimalDualHybridGradientParams& params,
std::optional<PrimalAndDualSolution> initial_solution,
const std::atomic<bool>* interrupt_solve = nullptr,
std::function<void(const std::string&)> message_callback = nullptr,
std::function<void(const IterationCallbackInfo&)> iteration_stats_callback =
nullptr);

View File

@@ -47,6 +47,7 @@
#include "ortools/pdlp/solvers.pb.h"
#include "ortools/pdlp/termination.h"
#include "ortools/pdlp/test_util.h"
#include "testing/base/public/googletest.h"
namespace operations_research::pdlp {
namespace {
@@ -443,6 +444,36 @@ TEST_P(PrimalDualHybridGradientLPTest, InfeasibleDual) {
EXPECT_LE(output.solve_log.iteration_count(), iteration_upperbound);
}
TEST_P(PrimalDualHybridGradientLPTest, InfeasiblePrimalWithReducedCosts) {
// A trivial LP that has reduced costs within the dual ray.
// min x
// Constraint: 2 <= x
// Variable: 0 <= x <= 1
QuadraticProgram lp(1, 1);
lp.objective_vector = Eigen::VectorXd{{1}};
lp.constraint_lower_bounds = Eigen::VectorXd{{2}};
lp.constraint_upper_bounds =
Eigen::VectorXd{{std::numeric_limits<double>::infinity()}};
lp.variable_lower_bounds = Eigen::VectorXd{{0}};
lp.variable_upper_bounds = Eigen::VectorXd{{1}};
lp.constraint_matrix.coeffRef(0, 0) = 1.0;
lp.constraint_matrix.makeCompressed();
const int iteration_upperbound = 100;
PrimalDualHybridGradientParams params =
CreateSolverParamsForFixture(iteration_upperbound,
/*eps_optimal_absolute=*/1.0e-6);
params.set_major_iteration_frequency(5);
params.mutable_termination_criteria()->set_eps_primal_infeasible(1.0e-6);
SolverResult output = PrimalDualHybridGradient(lp, params);
EXPECT_EQ(output.solve_log.termination_reason(),
TERMINATION_REASON_PRIMAL_INFEASIBLE);
const auto& dual = output.dual_solution;
EXPECT_GT(dual[0], 0.0);
EXPECT_EQ(output.reduced_costs[0], -dual[0]);
EXPECT_LT(output.solve_log.iteration_count(), iteration_upperbound);
}
TEST_P(PrimalDualHybridGradientLPTest, InfeasiblePrimalDual) {
const int iteration_upperbound = 600;
PrimalDualHybridGradientParams params =
@@ -1425,6 +1456,7 @@ TEST(PrimalDualHybridGradientTest, CallsCallback) {
int callback_count = 0;
SolverResult output = PrimalDualHybridGradient(
TestLp(), params, /*interrupt_solve=*/nullptr,
/*message_callback=*/nullptr,
[&callback_count](const IterationCallbackInfo& callback_info) {
++callback_count;
});
@@ -1525,7 +1557,8 @@ TEST(PrimalDualHybridGradientTest, RespectsInterruptFromCallback) {
};
const SolverResult output =
PrimalDualHybridGradient(TestLp(), params, &interrupt_solve, callback);
PrimalDualHybridGradient(TestLp(), params, &interrupt_solve,
/*message_callback=*/nullptr, callback);
EXPECT_EQ(output.solve_log.termination_reason(),
TERMINATION_REASON_INTERRUPTED_BY_USER);
EXPECT_GE(output.solve_log.iteration_count(), 10);
@@ -1921,6 +1954,7 @@ TEST_F(FeasibilityPolishingPrimalTest, CallsCallbackForAllThreePhases) {
bool polishing_termination_is_last = false;
SolverResult output = PrimalDualHybridGradient(
lp_, params_, /*interrupt_solve=*/nullptr,
/*message_callback=*/nullptr,
[&](const IterationCallbackInfo& callback_info) {
++callback_count[callback_info.iteration_type];
polishing_termination_is_last =
@@ -2061,6 +2095,36 @@ TEST(PresolveTest, PresolveInfeasible) {
EXPECT_EQ(output.solve_log.iteration_count(), 0);
}
TEST(PresolveTest, WarnsInitialSolution) {
PrimalAndDualSolution initial_solution;
initial_solution.primal_solution.setZero(6);
initial_solution.dual_solution.setZero(3);
PrimalDualHybridGradientParams params;
params.mutable_presolve_options()->set_use_glop(true);
CaptureTestStdout();
SolverResult output = PrimalDualHybridGradient(CorrelationClusteringStarLp(),
params, initial_solution);
EXPECT_THAT(GetCapturedTestStdout(),
AllOf(HasSubstr("initial solution"), HasSubstr("presolve")));
}
TEST(PresolveTest, WarnsInitialSolutionViaCallback) {
PrimalAndDualSolution initial_solution;
initial_solution.primal_solution.setZero(6);
initial_solution.dual_solution.setZero(3);
PrimalDualHybridGradientParams params;
params.mutable_presolve_options()->set_use_glop(true);
std::vector<std::string> message_output;
SolverResult output = PrimalDualHybridGradient(
CorrelationClusteringStarLp(), params, initial_solution,
/*interrupt_solve=*/nullptr,
[&message_output](const std::string& message) {
message_output.push_back(message);
});
EXPECT_THAT(message_output, Contains(AllOf(HasSubstr("initial solution"),
HasSubstr("presolve"))));
}
TEST_P(PresolveDualScalingTest, Dualize) {
auto [dualize, negate_and_scale_objective] = GetParam();
PrimalDualHybridGradientParams params;

View File

@@ -723,13 +723,26 @@ bool HasValidBounds(const ShardedQuadraticProgram& sharded_qp) {
}
void ProjectToPrimalVariableBounds(const ShardedQuadraticProgram& sharded_qp,
VectorXd& primal) {
VectorXd& primal,
const bool use_feasibility_bounds) {
const auto make_finite_values_zero = [](const double x) {
return std::isfinite(x) ? 0.0 : x;
};
sharded_qp.PrimalSharder().ParallelForEachShard(
[&](const Sharder::Shard& shard) {
const QuadraticProgram& qp = sharded_qp.Qp();
shard(primal) = shard(primal)
.cwiseMin(shard(qp.variable_upper_bounds))
.cwiseMax(shard(qp.variable_lower_bounds));
if (use_feasibility_bounds) {
shard(primal) =
shard(primal)
.cwiseMin(shard(qp.variable_upper_bounds)
.unaryExpr(make_finite_values_zero))
.cwiseMax(shard(qp.variable_lower_bounds)
.unaryExpr(make_finite_values_zero));
} else {
shard(primal) = shard(primal)
.cwiseMin(shard(qp.variable_upper_bounds))
.cwiseMax(shard(qp.variable_lower_bounds));
}
});
}
@@ -742,6 +755,9 @@ void ProjectToDualVariableBounds(const ShardedQuadraticProgram& sharded_qp,
const auto upper_bound_shard = shard(qp.constraint_upper_bounds);
auto dual_shard = shard(dual);
// TODO(user): Investigate whether it is more efficient to
// use .cwiseMax() + .cwiseMin() with unaryExpr(s) that map
// upper_bound_shard and lower_bound_shard to appropriate values.
for (int64_t i = 0; i < dual_shard.size(); ++i) {
if (!std::isfinite(upper_bound_shard[i])) {
dual_shard[i] = std::max(dual_shard[i], 0.0);

View File

@@ -190,8 +190,11 @@ SingularValueAndIterations EstimateMaximumSingularValueOfConstraintMatrix(
bool HasValidBounds(const ShardedQuadraticProgram& sharded_qp);
// Projects `primal` onto the variable bounds constraints.
// If `use_feasibility_bounds == true`, all finite variable bounds are replaced
// by zero.
void ProjectToPrimalVariableBounds(const ShardedQuadraticProgram& sharded_qp,
Eigen::VectorXd& primal);
Eigen::VectorXd& primal,
bool use_feasibility_bounds = false);
// Projects `dual` onto the dual variable bounds; see
// https://developers.google.com/optimization/lp/pdlp_math#dual_variable_bounds.

View File

@@ -631,6 +631,14 @@ TEST(ProjectToPrimalVariableBoundsTest, TestLp) {
EXPECT_THAT(primal, ElementsAre(-3, -2, 5, 3.5));
}
TEST(ProjectToPrimalVariableBoundsTest, TestLpWithFeasibility) {
ShardedQuadraticProgram qp(TestLp(), /*num_threads=*/2,
/*num_shards=*/2);
Eigen::VectorXd primal{{-3, -3, 5, 5}};
ProjectToPrimalVariableBounds(qp, primal, /*use_feasibility_bounds=*/true);
EXPECT_THAT(primal, ElementsAre(-3, 0, 0, 0));
}
TEST(ProjectToDualVariableBoundsTest, TestLp) {
ShardedQuadraticProgram qp(TestLp(), /*num_threads=*/2,
/*num_shards=*/2);

View File

@@ -26,6 +26,7 @@
#include "ortools/base/threadpool.h"
#include "ortools/pdlp/quadratic_program.h"
#include "ortools/pdlp/sharder.h"
#include "ortools/util/logging.h"
namespace operations_research::pdlp {
@@ -35,7 +36,8 @@ namespace {
// a single column.
void WarnIfMatrixUnbalanced(
const Eigen::SparseMatrix<double, Eigen::ColMajor, int64_t>& matrix,
absl::string_view matrix_name, int64_t density_limit) {
absl::string_view matrix_name, int64_t density_limit,
operations_research::SolverLogger* logger) {
if (matrix.cols() == 0) return;
int64_t worst_column = 0;
for (int64_t col = 0; col < matrix.cols(); ++col) {
@@ -46,23 +48,35 @@ void WarnIfMatrixUnbalanced(
if (matrix.col(worst_column).nonZeros() > density_limit) {
// TODO(user): fix this automatically in presolve instead of asking the
// user to do it.
LOG(WARNING)
<< "The " << matrix_name << " has "
<< matrix.col(worst_column).nonZeros() << " non-zeros in its "
<< worst_column
<< "th column. For best parallelization all rows and columns should "
"have at most order "
<< density_limit
<< " non-zeros. Consider rewriting the QP to split the corresponding "
"variable or constraint.";
if (logger) {
SOLVER_LOG(
logger, "WARNING: The ", matrix_name, " has ",
matrix.col(worst_column).nonZeros(), " non-zeros in its ",
worst_column,
"th column. For best parallelization all rows and columns should "
"have at most order ",
density_limit,
" non-zeros. Consider rewriting the QP to split the corresponding "
"variable or constraint.");
} else {
LOG(WARNING)
<< "The " << matrix_name << " has "
<< matrix.col(worst_column).nonZeros() << " non-zeros in its "
<< worst_column
<< "th column. For best parallelization all rows and columns should "
"have at most order "
<< density_limit
<< " non-zeros. Consider rewriting the QP to split the corresponding "
"variable or constraint.";
}
}
}
} // namespace
ShardedQuadraticProgram::ShardedQuadraticProgram(QuadraticProgram qp,
const int num_threads,
const int num_shards)
ShardedQuadraticProgram::ShardedQuadraticProgram(
QuadraticProgram qp, const int num_threads, const int num_shards,
operations_research::SolverLogger* logger)
: qp_(std::move(qp)),
transposed_constraint_matrix_(qp_.constraint_matrix.transpose()),
thread_pool_(num_threads == 1
@@ -85,10 +99,10 @@ ShardedQuadraticProgram::ShardedQuadraticProgram(QuadraticProgram qp,
qp_.constraint_lower_bounds.size();
const int64_t column_density_limit = work_per_iteration / num_threads;
WarnIfMatrixUnbalanced(qp_.constraint_matrix, "constraint matrix",
column_density_limit);
column_density_limit, logger);
WarnIfMatrixUnbalanced(transposed_constraint_matrix_,
"transposed constraint matrix",
column_density_limit);
"transposed constraint matrix", column_density_limit,
logger);
}
}

View File

@@ -24,6 +24,7 @@
#include "ortools/base/threadpool.h"
#include "ortools/pdlp/quadratic_program.h"
#include "ortools/pdlp/sharder.h"
#include "ortools/util/logging.h"
namespace operations_research::pdlp {
@@ -37,7 +38,10 @@ class ShardedQuadraticProgram {
public:
// Requires `num_shards` >= `num_threads` >= 1.
// Note that the `qp` is intentionally passed by value.
ShardedQuadraticProgram(QuadraticProgram qp, int num_threads, int num_shards);
// If `logger` is not nullptr, warns about unbalanced matrices using it;
// otherwise warns via Google standard logging.
ShardedQuadraticProgram(QuadraticProgram qp, int num_threads, int num_shards,
operations_research::SolverLogger* logger = nullptr);
// Movable but not copyable.
ShardedQuadraticProgram(const ShardedQuadraticProgram&) = delete;

View File

@@ -204,12 +204,13 @@ message ConvergenceInformation {
// infeasibility (i.e. has no solution); see also TerminationCriteria.
message InfeasibilityInformation {
// Let x_ray be the algorithm's estimate of the primal extreme ray where x_ray
// is a vector scaled such that its infinity norm is one. A simple and typical
// is a vector that satisfies the sign constraints for a ray, scaled such that
// its infinity norm is one (the sign constraints are the variable bound
// constraints, with all finite bounds mapped to zero). A simple and typical
// choice of x_ray is x_ray = x / | x |_∞ where x is the current primal
// iterate. For this value compute the maximum absolute error in the primal
// linear program with the right hand side and finite variable bounds set to
// zero. This error refers to both the linear constraints and sign constraints
// on the ray.
// iterate projected onto the primal ray sign constraints. For this value
// compute the maximum absolute error in the primal linear program with the
// right hand side set to zero.
optional double max_primal_ray_infeasibility = 1;
// The value of the linear part of the primal objective (ignoring additive
@@ -223,12 +224,12 @@ message InfeasibilityInformation {
optional double primal_ray_quadratic_norm = 3;
// Let (y_ray, r_ray) be the algorithm's estimate of the dual and reduced cost
// extreme ray where (y_ray, r_ray) is a vector scaled such that its infinity
// norm is one. A simple and typical choice of y_ray is
// (y_ray, r_ray) = (y, r) / max(| y |_∞, | r |_∞) where y is the
// current dual iterate and r is the current dual reduced costs.
// Consider the quadratic program we are solving but with the objective (both
// quadratic and linear terms) set to zero. This forms a linear program
// extreme ray where (y_ray, r_ray) is a vector (satisfying the dual variable
// constraints) scaled such that its infinity norm is one. A simple and
// typical choice of y_ray is (y_ray, r_ray) = (y, r) / max(| y |_∞, | r |_∞)
// where y is the current dual iterate and r is the current dual reduced
// costs. Consider the quadratic program we are solving but with the objective
// (both quadratic and linear terms) set to zero. This forms a linear program
// (label this linear program (1)) with no objective. Take the dual of (1) and
// compute the maximum absolute value of the constraint error for
// (y_ray, r_ray) to obtain the value of max_dual_ray_infeasibility.