Add primal and dual tolerances in presolve

Project.toml

@@ -1,7 +1,7 @@
 name = "Tulip"
 uuid = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa"
 authors = ["Mathieu Tanneau <mathieu.tanneau@gmail.com>"]
-version = "0.7.4"
+version = "0.7.5"
 
 [deps]
 CodecBzip2 = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"

src/Presolve/Presolve.jl

@@ -1,5 +1,8 @@
 Base.@kwdef mutable struct PresolveOptions{T}
     Level::Int = 1  # Presolve level
+
+    TolerancePFeas::T = sqrt(eps(T))  # Primal feasibility tolerance
+    ToleranceDFeas::T = sqrt(eps(T))  # Dual   feasibility tolerance
 end
 
 """
@@ -28,6 +31,7 @@ whose dual writes
 mutable struct PresolveData{T}
     updated::Bool
     status::TerminationStatus
+    options::PresolveOptions{T}
 
     # Original problem
     pb0::ProblemData{T}
@@ -86,11 +90,15 @@ mutable struct PresolveData{T}
     # TODO: set of transformations for pre-post crush
     ops::Vector{PresolveTransformation{T}}
 
-    function PresolveData(pb::ProblemData{T}) where{T}
+    function PresolveData(
+        pb::ProblemData{T},
+        options::PresolveOptions{T}=PresolveOptions{T}()
+    ) where{T}
         ps = new{T}()
 
         ps.updated = false
         ps.status = Trm_Unknown
+        ps.options = options
 
         ps.pb0 = pb
         ps.pb_red = nothing
@@ -167,7 +175,7 @@ end
 
 # Extract pre-solved problem data, to be passed to the IPM solver
 function extract_reduced_problem!(ps::PresolveData{T}) where{T}
-    
+
     pb = ProblemData{T}()
 
     pb.ncon = sum(ps.rowflag)
@@ -193,7 +201,7 @@ function extract_reduced_problem!(ps::PresolveData{T}) where{T}
     inew = 0
     for (iold, row) in enumerate(ps.pb0.arows)
         ps.rowflag[iold] || continue
-        
+
         inew += 1
         # Compute new row
         rind = Vector{Int}(undef, ps.nzrow[iold])
@@ -219,11 +227,11 @@ function extract_reduced_problem!(ps::PresolveData{T}) where{T}
     jnew = 0
     for (jold, col) in enumerate(ps.pb0.acols)
         ps.colflag[jold] || continue
-        
+
         jnew += 1
         # Compute new row
         cind = Vector{Int}(undef, ps.nzcol[jold])
-        cval = Vector{T}(undef, ps.nzcol[jold]) 
+        cval = Vector{T}(undef, ps.nzcol[jold])
 
         k = 0
         for (iold, aij) in zip(col.nzind, col.nzval)
@@ -379,7 +387,7 @@ function presolve!(ps::PresolveData{T}) where{T}
 
     # Identify row singletons
     ps.row_singletons = [i for (i, nz) in enumerate(ps.nzrow) if ps.rowflag[i] && nz == 1]
-    
+
     # II. Passes
     ps.updated = true
     npasses = 0  # TODO: maximum number of passes
@@ -392,7 +400,7 @@ function presolve!(ps::PresolveData{T}) where{T}
         ps.status == Trm_Unknown || return ps.status
         remove_empty_columns!(ps)
         ps.status == Trm_Unknown || return ps.status
-        
+
 
         # Remove all fixed variables
         # TODO: remove empty variables as well
@@ -444,7 +452,7 @@ function presolve!(ps::PresolveData{T}) where{T}
         ps.solution.z_primal = ps.obj0
         ps.solution.z_dual = ps.obj0
     end
-    
+
     # Old <-> new index mapping
     compute_index_mapping!(ps)
 
@@ -500,7 +508,7 @@ function bounds_consistency_checks!(ps::PresolveData{T}) where{T}
             ps.status = Trm_PrimalInfeasible
             ps.updated = true
 
-            # Resize problem            
+            # Resize problem
             compute_index_mapping!(ps)
             resize!(ps.solution, ps.nrow, ps.ncol)
             ps.solution.x .= zero(T)
@@ -518,7 +526,7 @@ function bounds_consistency_checks!(ps::PresolveData{T}) where{T}
             i_ = ps.new_con_idx[i]
             ps.solution.y_lower[i_] = one(T)
             ps.solution.y_upper[i_] = one(T)
-            
+
             return
         end
     end
@@ -529,7 +537,7 @@ function bounds_consistency_checks!(ps::PresolveData{T}) where{T}
             ps.status = Trm_PrimalInfeasible
             ps.updated = true
 
-            # Resize problem            
+            # Resize problem
             compute_index_mapping!(ps)
             resize!(ps.solution, ps.nrow, ps.ncol)
             ps.solution.x .= zero(T)
@@ -547,7 +555,7 @@ function bounds_consistency_checks!(ps::PresolveData{T}) where{T}
             j_ = ps.new_var_idx[j]
             ps.solution.s_lower[j_] = one(T)
             ps.solution.s_upper[j_] = one(T)
-            
+
             return
         end
     end
@@ -677,7 +685,7 @@ function remove_dominated_columns!(ps::PresolveData{T}) where{T}
                 @debug "Col $j forces y$i >= $y_"
                 ps.ly[i] = max(ps.ly[i],  y_)
             end
-        
+
         elseif !isfinite(l) && isfinite(u)
             # Upper-bounded variable: `aij * yi ≥ cj`
             if aij > zero(T)
@@ -699,4 +707,4 @@ function remove_dominated_columns!(ps::PresolveData{T}) where{T}
         ps.status == Trm_Unknown || break
     end
     return nothing
-end
+end

src/Presolve/empty_column.jl

@@ -13,7 +13,9 @@ function remove_empty_column!(ps::PresolveData{T}, j::Int) where{T}
     cj = ps.obj[j]
     @debug "Removing empty column $j" cj lb ub
 
-    if cj > zero(T)
+    ϵ = ps.options.ToleranceDFeas
+
+    if cj > ϵ
         if isfinite(lb)
             # Set variable to lower bound
             # Update objective constant
@@ -45,7 +47,7 @@ function remove_empty_column!(ps::PresolveData{T}, j::Int) where{T}
 
             return nothing
         end
-    elseif cj < zero(T)
+    elseif cj < -ϵ
         if isfinite(ub)
             # Set variable to upper bound
             # Update objective constant
@@ -74,7 +76,7 @@ function remove_empty_column!(ps::PresolveData{T}, j::Int) where{T}
             ps.solution.z_primal = ps.solution.z_dual = -T(Inf)
             j_ = ps.new_var_idx[j]
             ps.solution.x[j_] = one(T)
-            
+
             return
         end
     else
@@ -87,7 +89,6 @@ function remove_empty_column!(ps::PresolveData{T}, j::Int) where{T}
             # Free variable with zero coefficient and empty column
             push!(ps.ops, EmptyColumn(j, zero(T), zero(T)))
         end
-
     end
 
     # Book=keeping
@@ -102,4 +103,4 @@ function postsolve!(sol::Solution{T}, op::EmptyColumn{T}) where{T}
     sol.s_lower[op.j] = pos_part(op.s)
     sol.s_upper[op.j] = neg_part(op.s)
     return nothing
-end
+end

src/Presolve/empty_row.jl

@@ -12,8 +12,9 @@ function remove_empty_row!(ps::PresolveData{T}, i::Int) where{T}
 
     # Check bounds
     lb, ub = ps.lrow[i], ps.urow[i]
+    ϵ = ps.options.TolerancePFeas
 
-    if ub < zero(T)
+    if ub < -ϵ
         # Infeasible
         @debug "Row $i is primal infeasible"
         ps.status = Trm_PrimalInfeasible
@@ -37,7 +38,7 @@ function remove_empty_row!(ps::PresolveData{T}, i::Int) where{T}
         i_ = ps.new_con_idx[i]
         ps.solution.y_upper[i_] = one(T)
         return
-    elseif lb > zero(T)
+    elseif lb > ϵ
         @debug "Row $i is primal infeasible"
         ps.status = Trm_PrimalInfeasible
         ps.updated = true
@@ -50,7 +51,7 @@ function remove_empty_row!(ps::PresolveData{T}, i::Int) where{T}
         ps.solution.y_upper .= zero(T)
         ps.solution.s_lower .= zero(T)
         ps.solution.s_upper .= zero(T)
-        
+
         # Farkas ray: y⁺_i = 1 (any > 0 value works)
         ps.solution.primal_status = Sln_Unknown
         ps.solution.dual_status = Sln_InfeasibilityCertificate

src/Tulip.jl

@@ -7,7 +7,7 @@ using SparseArrays
 
 using TimerOutputs
 
-version() = v"0.7.4"
+version() = v"0.7.5"
 
 include("utils.jl")
 

src/model.jl

@@ -55,7 +55,7 @@ function Base.empty!(m::Model{T}) where{T}
     m.presolve_data = nothing
     m.solver = nothing
     m.solution = nothing
-    
+
     return nothing
 end
 
@@ -71,7 +71,7 @@ function optimize!(model::Model{T}) where{T}
         "Number of threads must be > 0 (is $(model.params.Threads))"
     )
     BLAS.set_num_threads(model.params.Threads)
-    
+
     # Print initial stats
     if model.params.OutputLevel > 0
         @printf "\nProblem info\n"
@@ -80,13 +80,13 @@ function optimize!(model::Model{T}) where{T}
         @printf "  Variables   : %d\n" model.pbdata.nvar
         @printf "  Non-zeros   : %d\n" sum(length.([col.nzind for col in model.pbdata.acols]))
     end
-    
+
     pb_ = model.pbdata
     # Presolve
     # TODO: improve the if-else
     ps_options = model.params.Presolve
     if ps_options.Level > 0
-        model.presolve_data = PresolveData(model.pbdata)
+        model.presolve_data = PresolveData(model.pbdata, ps_options)
         t_ = @elapsed st = presolve!(model.presolve_data)
         model.status = st
 
@@ -104,7 +104,7 @@ function optimize!(model::Model{T}) where{T}
         # Check presolve status
         if st == Trm_Optimal || st == Trm_PrimalInfeasible || st == Trm_DualInfeasible || st == Trm_PrimalDualInfeasible
             model.params.OutputLevel > 0 && println("Presolve solved the problem.")
-            
+
             # Perform post-solve
             sol0 = Solution{T}(model.pbdata.ncon, model.pbdata.nvar)
             postsolve!(sol0, model.presolve_data.solution, model.presolve_data)
@@ -169,7 +169,7 @@ function _extract_solution!(sol::Solution{T},
     sol.is_primal_ray = is_primal_ray
     sol.is_dual_ray = is_dual_ray
     τ_ = (is_primal_ray || is_dual_ray) ? one(T) : inv(ipm.pt.τ)
-    
+
     @. sol.x = ipm.pt.x[1:n] * τ_
     @. sol.s_lower = ipm.pt.zl[1:n] * τ_
     @. sol.s_upper = ipm.pt.zu[1:n] * τ_
@@ -212,4 +212,4 @@ function _extract_solution!(sol::Solution{T},
     end
 
     return nothing
-end
+end

test/Presolve/empty_row.jl

@@ -57,9 +57,6 @@ function empty_row_tests(T::Type)
     #   (current problem only has 1 row)
     @test sol.y_lower[1] >  zero(T)
 
-    test_empty_row_1(T)
-    test_empty_row_2(T)
-
     return
 end
 
@@ -95,7 +92,7 @@ function test_empty_row_1(T::Type)
     sol = ps.solution
     @test sol.dual_status == Tulip.Sln_InfeasibilityCertificate
     @test sol.z_primal == sol.z_dual == T(Inf)
-    
+
     # Check Farkas ray
     #   (current problem only has 1 row)
     @test sol.y_lower[1] >  zero(T)
@@ -134,16 +131,70 @@ function test_empty_row_2(T::Type)
     sol = ps.solution
     @test sol.dual_status == Tulip.Sln_InfeasibilityCertificate
     @test sol.z_primal == sol.z_dual == T(Inf)
-    
+
     # Check Farkas ray
     #   (current problem only has 1 row)
     @test sol.y_upper[1] >  zero(T)
 
     return nothing
 end
 
+function test_empty_row_tolerances(T::Type)
+    # Adapted from https://github.com/ds4dm/Tulip.jl/issues/98
+    #=
+        min     x + y + z
+        s.t.    x + y + z == 1
+                x == ¹/₃
+                y == ¹/₃
+                z == ¹/₃
+                x, y, z, ≥ 0
+
+        In the absence of numerical tolerances, x, y, and z get eliminated,
+        but rouding errors cause the first constraint to be 0 == ϵ ≈ 1e-16,
+        thereby rendering the problem infeasible.
+    =#
+    pb = Tulip.ProblemData{T}()
+
+    m, n = 4, 3
+    A = sparse(
+        [1, 1, 1, 2, 3, 4],
+        [1, 2, 3, 1, 2, 3],
+        T[1, 1, 1, 1, 1, 1],
+        m, n
+    )
+    c = ones(T, n)
+
+    Tulip.load_problem!(pb, "test",
+        true, c, zero(T),
+        A, T[1, 1//3, 1//3, 1//3], T[1, 1//3, 1//3, 1//3],
+        zeros(T, n), fill(T(Inf), n),
+        ["row1", "row2", "row3", "row4"], ["x", "y", "z"]
+    )
+
+    ps = Tulip.PresolveData(pb)
+    Tulip.presolve!(ps)
+
+    @test ps.status == Tulip.Trm_Optimal
+    @test ps.nrow == 0
+    @test ps.ncol == 0
+
+    # Check solution status & objective value
+    sol = ps.solution
+    @test sol.primal_status == Tulip.Sln_Optimal
+    @test sol.dual_status == Tulip.Sln_Optimal
+    @test sol.z_primal ≈ 1
+    @test sol.z_dual   ≈ 1
+
+    return nothing
+end
+
 @testset "Empty row" begin
     for T in TvTYPES
-        @testset "$T" begin empty_row_tests(T) end
+        @testset "$T" begin
+            empty_row_tests(T)
+            test_empty_row_1(T)
+            test_empty_row_2(T)
+            test_empty_row_tolerances(T)
+        end
     end
-end
+end

test/runtests.jl

@@ -10,9 +10,9 @@ import Convex
 
 const TvTYPES = [Float32, Float64, BigFloat]
 
-# Check That Tulip.version() matches what's in the Project.toml
 @testset "Tulip" begin
 
+# Check That Tulip.version() matches what's in the Project.toml
 tlp_ver = Tulip.version()
 toml_ver = TOML.parsefile("../Project.toml")["version"]
 @test tlp_ver == VersionNumber(toml_ver)