Add docs and tests for infeasibility certificates

docs/make.jl

@@ -37,6 +37,7 @@ const _PAGES = [
     ],
     "Background" => [
         "background/duality.md",
+        "background/infeasibility_certificates.md",
         "background/naming_conventions.md",
     ],
     "API Reference" => [

docs/src/background/infeasibility_certificates.md

@@ -0,0 +1,160 @@
+```@meta
+CurrentModule = MathOptInterface
+DocTestSetup = quote
+    using MathOptInterface
+    const MOI = MathOptInterface
+end
+DocTestFilters = [r"MathOptInterface|MOI"]
+```
+
+# Infeasibility certificates
+
+When given a conic problem that is infeasible or unbounded, some solvers can
+produce a certificate of infeasibility. This page explains what a certificate of
+infeasibility is, and the related conventions that MathOptInterface adopts.
+
+## Conic duality
+
+MathOptInterface uses conic duality to define infeasibility certificates. A full
+explanation is given in the section [Duality](@ref), but here is a brief
+overview.
+
+### Minimization problems
+
+For a minimization problem in geometric conic form, the primal is:
+```math
+\begin{align}
+& \min_{x \in \mathbb{R}^n} & a_0^\top x + b_0
+\\
+& \;\;\text{s.t.} & A_i x + b_i & \in \mathcal{C}_i & i = 1 \ldots m,
+\end{align}
+```
+and the dual is a maximization problem in standard conic form:
+```math
+\begin{align}
+& \max_{y_1, \ldots, y_m} & -\sum_{i=1}^m b_i^\top y_i + b_0
+\\
+& \;\;\text{s.t.} & a_0 - \sum_{i=1}^m A_i^\top y_i & = 0
+\\
+& & y_i & \in \mathcal{C}_i^* & i = 1 \ldots m,
+\end{align}
+```
+where each ``\mathcal{C}_i`` is a closed convex cone and ``\mathcal{C}_i^*`` is
+its dual cone.
+
+### Maximization problems
+
+For a maximization problem in geometric conic form, the primal is:
+```math
+\begin{align}
+& \max_{x \in \mathbb{R}^n} & a_0^\top x + b_0
+\\
+& \;\;\text{s.t.} & A_i x + b_i & \in \mathcal{C}_i & i = 1 \ldots m,
+\end{align}
+```
+and the dual is a minimization problem in standard conic form:
+```math
+\begin{align}
+& \min_{y_1, \ldots, y_m} & \sum_{i=1}^m b_i^\top y_i + b_0
+\\
+& \;\;\text{s.t.} & a_0 + \sum_{i=1}^m A_i^\top y_i & = 0
+\\
+& & y_i & \in \mathcal{C}_i^* & i = 1 \ldots m.
+\end{align}
+```
+
+## Unbounded problems
+
+A problem is unbounded if and only if:
+ 1. there exists a feasible primal solution
+ 2. the dual is infeasible.
+
+A feasible primal solution—if one exists—can be obtained by setting
+[`ObjectiveSense`](@ref) to `FEASIBILITY_SENSE` before optimizing. Therefore,
+most solvers terminate after they prove the dual is infeasible via a certificate
+of dual infeasibility, but _before_ they have found a feasible primal solution.
+This is also the reason that MathOptInterface defines the `DUAL_INFEASIBLE`
+status instead of `UNBOUNDED`.
+
+A certificate of dual infeasibility is an improving ray of the primal problem.
+That is, there exists some vector ``d`` such that for all ``\eta > 0``:
+```math
+A_i (x + \eta d) + b_i \in \mathcal{C}_i,\ \ i = 1 \ldots m,
+```
+and (for minimization problems):
+```math
+a_0^\top (x + \eta d) + b_0 < a_0^\top x + b_0,
+```
+for any feasible point ``x``. The latter simplifies to ``a_0^\top d < 0``. For
+maximization problems, the inequality is reversed, so that ``a_0^\top d > 0``.
+
+If the solver has found a certificate of dual infeasibility:
+
+ * [`TerminationStatus`](@ref) must be `DUAL_INFEASIBLE`
+ * [`PrimalStatus`](@ref) must be `INFEASIBILITY_CERTIFICATE`
+ * [`VariablePrimal`](@ref) must be the corresponding value of ``d``
+ * [`ConstraintPrimal`](@ref) must be the corresponding value of ``A_i d``
+ * [`ObjectiveValue`](@ref) must be the value ``a_0^\top d``. Note that this is
+   the value of the objective function at `d`, ignoring the constant `b_0`.
+
+!!! note
+    The choice of whether to scale the ray ``d`` to have magnitude `1` is left
+    to the solver.
+
+## Infeasible problems
+
+A certificate of primal infeasibility is an improving ray of the dual problem.
+However, because infeasibility is independent of the objective function, we
+first homogenize the primal problem by removing its objective.
+
+For a minimization problem, a dual improving ray is some vector ``d`` such that
+for all ``\eta > 0``:
+```math
+\begin{align}
+-\sum_{i=1}^m A_i^\top (y_i + \eta d_i) & = 0 \\
+(y_i + \eta d_i) & \in \mathcal{C}_i^* & i = 1 \ldots m,
+\end{align}
+```
+and:
+```math
+-\sum_{i=1}^m b_i^\top (y_i + \eta d_i) > -\sum_{i=1}^m b_i^\top y_i,
+```
+for any feasible dual solution ``y``. The latter simplifies to
+``-\sum_{i=1}^m b_i^\top d_i > 0``. For a maximization problem, the inequality
+is ``\sum_{i=1}^m b_i^\top d_i < 0``. (Note that these are the same inequality,
+modulo a `-` sign.)
+
+If the solver has found a certificate of primal infeasibility:
+
+ * [`TerminationStatus`](@ref) must be `INFEASIBLE`
+ * [`DualStatus`](@ref) must be `INFEASIBILITY_CERTIFICATE`
+ * [`ConstraintDual`](@ref) must be the corresponding value of ``d``
+ * [`DualObjectiveValue`](@ref) must be the value
+   ``-\sum_{i=1}^m b_i^\top d_i`` for minimization problems and
+   ``\sum_{i=1}^m b_i^\top d_i`` for maximization problems.
+
+!!! note
+    The choice of whether to scale the ray ``d`` to have magnitude `1` is left
+    to the solver.
+
+### Infeasibility certificates of variable bounds
+
+Many linear solvers (e.g., Gurobi) do not provide explicit access to the primal
+infeasibility certificate of a variable bound. However, given a set of linear
+constraints:
+```math
+\begin{align}
+l_A \le A x \le u_A \\
+l_x \le x \le u_x,
+\end{align}
+```
+the primal certificate of the variable bounds can be computed using the primal
+certificate associated with the affine constraints, ``d``. (Note that ``d`` will
+have one element for each row of the ``A`` matrix, and that some or all of the
+elements in the vectors ``l_A`` and ``u_A`` may be ``\pm \infty``. If both
+``l_A`` and ``u_A`` are finite for some row, the corresponding element in ``d`
+ must be `0`.)
+
+Given ``d``, compute ``\bar{d} = d^\top A``. If the bound is finite, a
+certificate for the lower variable bound of ``x_i`` is ``\max\{\bar{d}_i, 0\}``,
+and a certificate for the upper variable bound is ``\min\{\bar{d}_i, 0\}``.

src/Bridges/Constraint/scalarize.jl

@@ -175,7 +175,6 @@ function MOI.set(
     bridge::ScalarizeBridge,
     value,
 )
-    # TODO do no add constant if the primal status is a ray like in Vectorize
     MOI.set.(model, attr, bridge.scalar_constraints, value .- bridge.constants)
     return
 end

src/Test/test_infeasibility_certificates.jl

@@ -0,0 +1,187 @@
+for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2)
+    o = iszero(offset) ? "" : "_offset"
+    f_unbd_name = Symbol("test_unbounded_$(Symbol(sense))$(o)")
+    f_infeas_name = Symbol("test_infeasible_$(Symbol(sense))$(o)")
+    f_infeas_affine_name = Symbol("test_infeasible_affine_$(Symbol(sense))$(o)")
+    @eval begin
+        """
+        Test the infeasibility certificates of an unbounded model.
+        """
+        function $(f_unbd_name)(model::MOI.ModelLike, config::Config)
+            @requires _supports(config, MOI.optimize!)
+            x = MOI.add_variable(model)
+            MOI.set(model, MOI.ObjectiveSense(), $sense)
+            f = 2.2 * x + $offset
+            MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+            c = if $sense == MOI.MIN_SENSE
+                MOI.add_constraint(model, 1.3 * x, MOI.LessThan(1.1))
+            else
+                MOI.add_constraint(model, 1.3 * x, MOI.GreaterThan(1.1))
+            end
+            MOI.optimize!(model)
+            @requires(
+                MOI.get(model, MOI.TerminationStatus()) == MOI.DUAL_INFEASIBLE,
+            )
+            @requires(
+                MOI.get(model, MOI.PrimalStatus()) ==
+                MOI.INFEASIBILITY_CERTIFICATE,
+            )
+            d = MOI.get(model, MOI.VariablePrimal(), x)
+            obj = MOI.get(model, MOI.ObjectiveValue())
+            if $sense == MOI.MIN_SENSE
+                @test d < -config.atol
+                @test obj < -config.atol
+            else
+                @test d > config.atol
+                @test obj > config.atol
+            end
+            @test isapprox(2.2 * d, obj, config)
+            Ad = MOI.get(model, MOI.ConstraintPrimal(), c)
+            @test isapprox(1.3 * d, Ad, config)
+            return
+        end
+
+        function setup_test(
+            ::typeof($(f_unbd_name)),
+            mock::MOIU.MockOptimizer,
+            config::Config,
+        )
+            MOIU.set_mock_optimize!(
+                mock,
+                (mock::MOIU.MockOptimizer) -> begin
+                    MOIU.mock_optimize!(
+                        mock,
+                        MOI.DUAL_INFEASIBLE,
+                        (
+                            MOI.INFEASIBILITY_CERTIFICATE,
+                            [$sense == MOI.MIN_SENSE ? -1.0 : 1.0],
+                        ),
+                        MOI.NO_SOLUTION,
+                    )
+                end,
+            )
+            return
+        end
+
+        """
+        Test the infeasibility certificates of an infeasible model.
+        """
+        function $(f_infeas_name)(model::MOI.ModelLike, config::Config)
+            @requires _supports(config, MOI.optimize!)
+            x = MOI.add_variable(model)
+            MOI.set(model, MOI.ObjectiveSense(), $sense)
+            f = 2.2 * x + $offset
+            MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+            cl = MOI.add_constraint(model, x, MOI.GreaterThan(2.5))
+            cu = MOI.add_constraint(model, 1.0 * x, MOI.LessThan(1.4))
+            MOI.optimize!(model)
+            @requires(
+                MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE,
+            )
+            @requires(
+                MOI.get(model, MOI.DualStatus()) ==
+                MOI.INFEASIBILITY_CERTIFICATE,
+            )
+            dl = MOI.get(model, MOI.ConstraintDual(), cl)
+            du = MOI.get(model, MOI.ConstraintDual(), cu)
+            obj = MOI.get(model, MOI.DualObjectiveValue())
+            if $sense == MOI.MIN_SENSE
+                @test obj > config.atol
+                @test isapprox(-(1.4 * dl + 2.5 * du), obj, config)
+            else
+                @test obj < -config.atol
+                @test isapprox(1.4 * dl + 2.5 * du, obj, config)
+            end
+            @test dl > config.atol
+            @test du < -config.atol
+            @test isapprox(dl + du, 0.0, config)
+            return
+        end
+
+        function setup_test(
+            ::typeof($(f_infeas_name)),
+            mock::MOIU.MockOptimizer,
+            config::Config,
+        )
+            MOIU.set_mock_optimize!(
+                mock,
+                (mock::MOIU.MockOptimizer) -> begin
+                    MOIU.mock_optimize!(
+                        mock,
+                        MOI.INFEASIBLE,
+                        MOI.NO_SOLUTION,
+                        MOI.INFEASIBILITY_CERTIFICATE,
+                        (MOI.VariableIndex, MOI.GreaterThan{Float64}) =>
+                            [1.0],
+                        (
+                            MOI.ScalarAffineFunction{Float64},
+                            MOI.LessThan{Float64},
+                        ) => [-1.0],
+                    )
+                end,
+            )
+            return
+        end
+
+        """
+        Test the infeasibility certificates of an infeasible model with affine
+        constraints.
+        """
+        function $(f_infeas_affine_name)(model::MOI.ModelLike, config::Config)
+            @requires _supports(config, MOI.optimize!)
+            x = MOI.add_variable(model)
+            MOI.set(model, MOI.ObjectiveSense(), $sense)
+            f = 2.2 * x + $offset
+            MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+            cl = MOI.add_constraint(model, 1.0 * x, MOI.GreaterThan(2.5))
+            cu = MOI.add_constraint(model, x, MOI.LessThan(1.4))
+            MOI.optimize!(model)
+            @requires(
+                MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE,
+            )
+            @requires(
+                MOI.get(model, MOI.DualStatus()) ==
+                MOI.INFEASIBILITY_CERTIFICATE,
+            )
+            dl = MOI.get(model, MOI.ConstraintDual(), cl)
+            du = MOI.get(model, MOI.ConstraintDual(), cu)
+            obj = MOI.get(model, MOI.DualObjectiveValue())
+            if $sense == MOI.MIN_SENSE
+                @test obj > config.atol
+                @test isapprox(-(1.4 * dl + 2.5 * du), obj, config)
+            else
+                @test obj < -config.atol
+                @test isapprox(1.4 * dl + 2.5 * du, obj, config)
+            end
+            @test dl > config.atol
+            @test du < -config.atol
+            @test isapprox(dl + du, 0.0, config)
+            return
+        end
+
+        function setup_test(
+            ::typeof($(f_infeas_affine_name)),
+            mock::MOIU.MockOptimizer,
+            config::Config,
+        )
+            MOIU.set_mock_optimize!(
+                mock,
+                (mock::MOIU.MockOptimizer) -> begin
+                    MOIU.mock_optimize!(
+                        mock,
+                        MOI.INFEASIBLE,
+                        MOI.NO_SOLUTION,
+                        MOI.INFEASIBILITY_CERTIFICATE,
+                        (
+                            MOI.ScalarAffineFunction{Float64},
+                            MOI.GreaterThan{Float64},
+                        ) => [1.0],
+                        (MOI.VariableIndex, MOI.LessThan{Float64}) =>
+                            [-1.0],
+                    )
+                end,
+            )
+            return
+        end
+    end
+end