PDLP declared scaled a2864 instance infeasible

[ohinder]

Version: 9.10
Language: C++
Solver: PDLP
OS: Linux

The issue

PDLP is declaring a scaled version of the a2864 instance from the Mittelmann library infeasible.

The instance should be feasible. This was verified by solving the instance with Gurobi that found an optimal solution.

However, after 64 iterations PDLP declared the problem infeasible.

Relevant files

mps instance

decoded_solve_log.json

solver-output.txt

The suspected cause

Inspection of the code showed that the that the objective gradient was being included in the computation of the max_dual_ray_infeasibility, based on the logs it appears the algorithm stumbling across a point where the dual infeasibility (including the objective gradient) is zero and the dual objective is 1.68 * 10^{-15}. I believe there is a numerical error causing the slightly positive dual objective since it is very small (~10^{-15}). Replacing

http://github.com/google/or-tools/blob/93dafb79a0261b5211f4eb17437aeebd7f45e543/ortools/pdlp/iteration_stats.cc#L522C35-L522C66

with

result.set_dual_ray_objective(dual_ray_objective / l_inf_dual - 1e-14);

Fixed the issue and the algorithm rapidly found the optimal solution ~2624 iterations, using the command line interface, i.e.,

 "$path_to_OR_tools"/build/bin/pdlp_solve --input "$path_to_instances"/"$instance_name".mps --sol_file "$path_to_solution_files"/no-polish-"$instance_name".sol --solve_log_file "$path_to_log_files"/no-polish-"$instance_name".json --params "verbosity_level: 4 num_threads: 16 termination_criteria {detailed_optimality_criteria {eps_optimal_primal_residual_absolute: 1.0e-8 eps_optimal_primal_residual_relative: 0.0 eps_optimal_dual_residual_absolute: 1.0e-8 eps_optimal_dual_residual_relative: 0.0 eps_optimal_objective_gap_absolute: 0.0 eps_optimal_objective_gap_relative: 1.0e-2} eps_primal_infeasible: 1.0e-9 eps_dual_infeasible: 1.0e-9 optimality_norm: OPTIMALITY_NORM_L_INF} use_feasibility_polishing: false handle_some_primal_gradients_on_finite_bounds_as_residuals: false"

Other relevant material

Documentation on infeasibility certificates for PDLP
https://developers.google.com/optimization/lp/pdlp_math#infeasibility_identification

Comment on previous suspected cause

Initially, I though that the primal infeasibility certificate was not being computed with the homogenous residuals but actually it is just being computed in a slightly confusing way (and differently from the primal certificate). In particular,

or-tools/ortools/pdlp/iteration_stats.cc

Lines 505 to 528 in 93dafb7

	InfeasibilityInformation result;
	// Compute primal infeasibility information.
	VectorXd scaled_primal_gradient = PrimalGradientFromObjectiveProduct(
	scaled_sharded_qp, scaled_dual_ray,
	ZeroVector(scaled_sharded_qp.PrimalSharder()),
	/*use_zero_primal_objective=*/true);
	// 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,
	primal_solution_for_residual_tests, scaled_primal_gradient,
	/*componentwise_residual_offset=*/0.0);
	double dual_ray_objective =
	DualObjectiveBoundsTerm(scaled_sharded_qp, scaled_dual_ray) +
	dual_residuals.objective_correction;
	if (l_inf_dual > 0) {
	result.set_dual_ray_objective(dual_ray_objective / l_inf_dual);
	result.set_max_dual_ray_infeasibility(dual_residuals.l_inf_residual /
	l_inf_dual);
	} else {
	result.set_dual_ray_objective(0.0);
	result.set_max_dual_ray_infeasibility(0.0);
	}

shows that the "scaled_primal_gradient" is being computed with the PrimalGradientFromObjectiveProduct function with the option use_zero_primal_objective=true. Thus this is really the scaled_homogenous_primal_gradient"

Dual infeasibility certificate computes the homogenous residuals through an option in the PrimalResidualNorms function:

or-tools/ortools/pdlp/iteration_stats.cc

Lines 66 to 70 in 93dafb7

	ResidualNorms PrimalResidualNorms(
	const ShardedQuadraticProgram& sharded_qp, const VectorXd& row_scaling_vec,
	const VectorXd& scaled_primal_solution,
	const double componentwise_residual_offset,
	bool use_homogeneous_constraint_bounds = false) {

or-tools/ortools/pdlp/iteration_stats.cc

Lines 530 to 536 in 93dafb7

	// Compute dual infeasibility information. We don't use
	// `primal_residuals.l_inf_componentwise_residual`, so don't need to set
	// `componentwise_residual_offset` to a meaningful value.
	ResidualNorms primal_residuals =
	PrimalResidualNorms(scaled_sharded_qp, row_scaling_vec, scaled_primal_ray,
	/*componentwise_residual_offset=*/0.0,
	/*use_homogeneous_constraint_bounds=*/true);

I would suggest inside the ComputeInfeasibilityInformation function renaming scaled_primal_gradient to scaled_homogeneous_primal_gradient, dual_residuals to homogeneous_dual_residuals, and primal_residuals to homogeneous_primal_residuals to improve readability