fix postsolve, check primal constraints (#90)

Author: geoffroyleconte

src/presolve/free_rows.jl

@@ -0,0 +1,24 @@
+struct FreeRow{T, S} <: PresolveOperation{T, S}
+  i::Int # idx of free row
+end
+
+function free_rows!(
+  operations::Vector{PresolveOperation{T, S}},
+  lcon::S,
+  ucon::S,
+  ncon,
+  row_cnt,
+  kept_rows::Vector{Bool},
+) where {T, S}
+  free_row_pass = false
+  for i = 1:ncon
+    (kept_rows[i] && lcon[i] == -T(Inf) && ucon[i] == T(Inf)) || continue
+    free_row_pass = true
+    row_cnt[i] = -1
+    kept_rows[i] = false
+    push!(operations, FreeRow{T, S}(i))
+  end
+  return free_row_pass
+end
+
+postsolve!(pt::OutputPoint, operation::FreeRow) = nothing

src/presolve/postsolve_utils.jl

@@ -32,11 +32,11 @@ function restore_ilow_iupp!(ilow, iupp, kept_cols)
       offset += 1
     end
     if c_low ≤ nlow && ilow[c_low] + offset == i
-      ilow[c_low] += 1
+      ilow[c_low] += offset
       c_low += 1
     end
     if c_upp ≤ nupp && iupp[c_upp] + offset == i
-      iupp[c_upp] += 1
+      iupp[c_upp] += offset
       c_upp += 1
     end
   end

src/presolve/presolve.jl

@@ -70,6 +70,8 @@ include("empty_rows.jl")
 include("singleton_rows.jl")
 include("unconstrained_reductions.jl")
 include("linear_singleton_columns.jl")
+include("primal_constraints.jl")
+include("free_rows.jl")
 include("postsolve_utils.jl")
 
 mutable struct PresolvedData{T, S}
@@ -210,6 +212,23 @@ function presolve(
       kept_cols,
     )
 
+    infeasible_cst = primal_constraints!(
+      operations,
+      arows,
+      lcon,
+      ucon,
+      lvar,
+      uvar,
+      nvar,
+      ncon,
+      kept_rows,
+      kept_cols,
+      row_cnt,
+      col_cnt,
+    )
+
+    free_rows_pass = free_rows!(operations, lcon, ucon, ncon, row_cnt, kept_rows)
+
     psdata.c0, ifix_pass = remove_ifix!(
       operations,
       hcols,
@@ -228,11 +247,11 @@ function presolve(
       xps,
     )
 
-    infeasible = check_bounds(lvar, uvar, lcon, ucon, nvar, ncon, kept_rows, kept_cols)
+    infeasible_bnd = check_bounds(lvar, uvar, lcon, ucon, nvar, ncon, kept_rows, kept_cols)
 
-    keep_iterating =
-      (empty_row_pass || singl_row_pass || ifix_pass || free_lsc_pass) &&
-      (!unbounded || !infeasible)
+    infeasible = infeasible_bnd || infeasible_cst
+    reduction_pass = empty_row_pass || singl_row_pass || ifix_pass || free_lsc_pass || free_rows_pass
+    keep_iterating = reduction_pass && !unbounded && !infeasible
     keep_iterating && (nb_pass += 1)
   end
 
@@ -295,7 +314,8 @@ function presolve(
       multipliers_U = s_u,
       iter = 0,
       elapsed_time = time() - start_time,
-      solver_specific = Dict(:presolvedQM => nothing),
+      solver_specific = Dict(:presolvedQM => nothing,
+                             :psoperations => operations),
     )
   else
     psmeta = NLPModelMeta{T, S}(
@@ -319,7 +339,7 @@ function presolve(
       ps,
       iter = 0,
       elapsed_time = time() - start_time,
-      solver_specific = Dict(:presolvedQM => ps),
+      solver_specific = Dict(:presolvedQM => ps, :psoperations => operations),
     )
   end
 end
@@ -343,12 +363,12 @@ function postsolve!(
 end
 
 """
-    x, y, s_l, s_u = postsolve(qm::QuadraticModel{T, S}, psqm::PresolvedQuadraticModel{T, S}, 
-                               x_in::S, y_in::S,
-                               s_l_in::SparseVector{T, Int},
-                               s_u_in::SparseVector{T, Int}) where {T, S}
+    pt = postsolve(qm::QuadraticModel{T, S}, psqm::PresolvedQuadraticModel{T, S}, 
+                   x_in::S, y_in::S,
+                   s_l_in::SparseVector{T, Int},
+                   s_u_in::SparseVector{T, Int}) where {T, S}
 
-Retrieve the solution `x, y, s_l, s_u` of the original QP `qm` given the solution of the presolved QP (`psqm`)
+Retrieve the solution `(x, y, s_l, s_u)` of the original QP `qm` given the solution of the presolved QP (`psqm`)
 `x_in, y_in, s_l_in, s_u_in`.
 """
 function postsolve(

src/presolve/primal_constraints.jl

@@ -0,0 +1,62 @@
+function primal_constraints!(
+  operations::Vector{PresolveOperation{T, S}},
+  arows::Vector{Row{T}},
+  lcon::S,
+  ucon::S,
+  lvar::S,
+  uvar::S,
+  nvar,
+  ncon,
+  kept_rows,
+  kept_cols,
+  row_cnt,
+  col_cnt,
+) where {T, S}
+
+  for i=1:ncon
+    (kept_rows[i] && !(lcon[i] == -T(Inf) && ucon[i] == T(Inf))) || continue
+    rowi = arows[i]
+    uconi2 = zero(T)
+    lconi2 = zero(T)
+    for (j, aij) in zip(rowi.nzind, rowi.nzval)
+      kept_cols[j] || continue
+      uconi2 += (aij > zero(T)) ? aij * uvar[j] : aij * lvar[j]
+      lconi2 += (aij > zero(T)) ? aij * lvar[j] : aij * uvar[j]
+    end
+    if uconi2 < lcon[i] || lconi2 > ucon[i]
+      @warn "implied bounds for row $i lead to a primal infeasible constraint"
+      return true
+    elseif uconi2 == lcon[i]
+      for (j, aij) in zip(rowi.nzind, rowi.nzval)
+        kept_cols[j] || continue
+        if aij > zero(T)
+          lvar[j] = uvar[j]
+        else
+          uvar[j] = lvar[j]
+        end
+      end
+    elseif lconi2 == ucon[i]
+      for (j, aij) in zip(rowi.nzind, rowi.nzval)
+        kept_cols[j] || continue
+        if aij > zero(T)
+          uvar[j] = lvar[j]
+        else
+          lvar[j] = uvar[j]
+        end
+      end
+    elseif lcon[i] < lconi2 && uconi2 < ucon[i]
+      kept_rows[i] = false
+      row_cnt[i] = -1
+      for j in rowi.nzind
+        kept_cols[j] || continue
+        col_cnt[j] -= 1
+      end
+      push!(operations, EmptyRow{T, S}(i))
+    elseif lcon[i] < lconi2
+      lcon[i] = -T(Inf)
+    elseif uconi2 < ucon[i]
+      ucon[i] = T(Inf)
+    end
+  end
+  return false
+end

src/presolve/remove_ifix.jl

@@ -1,6 +1,9 @@
 struct RemoveIfix{T, S} <: PresolveOperation{T, S}
   j::Int
   xj::T
+  cj::T
+  acolj::Col{T}
+  hcolj::Col{T}
 end
 
 # ̃xᵀ̃Hx̃ + ̃ĉᵀx̃ + lⱼ²Hⱼⱼ + cⱼxⱼ + c₀
@@ -41,23 +44,42 @@ function remove_ifix!(
     end
 
     # remove ifix in A cols
-    for k = 1:length(acols[j].nzind)
-      i = acols[j].nzind[k]
+    c_acolj = 0
+    for (i, aij) in zip(acols[j].nzind, acols[j].nzval)
       if kept_rows[i]
         row_cnt[i] -= 1
-        con_offset = acols[j].nzval[k] * xj
+        con_offset = aij * xj
         lcon[i] -= con_offset
         ucon[i] -= con_offset
+        c_acolj += 1 # count number of rows in acols[j]
       end
     end
 
+    # store acolj for postsolve
+    acolj = Col(zeros(Int, c_acolj), zeros(T, c_acolj))
+    c_acolj = 1
+    for k in 1:length(acols[j].nzind)
+      i = acols[j].nzind[k]
+      kept_rows[i] || continue
+      aij = acols[j].nzval[k]
+      acolj.nzind[c_acolj] = i
+      acolj.nzval[c_acolj] = aij
+      c_acolj += 1
+    end
+
+    # store hcolj for postsolve
+    hcolj = Col{T}(
+      [i for i in hcols[j].nzind if kept_cols[i]],
+      [hij for (i, hij) in zip(hcols[j].nzind, hcols[j].nzval) if kept_cols[i]],
+    )
+
     # update c0 with c[j] coeff
     c0_offset += c[j] * xj
     xps[j] = xj
     kept_cols[j] = false
     col_cnt[j] = -1
 
-    push!(operations, RemoveIfix{T, S}(j, xj))
+    push!(operations, RemoveIfix{T, S}(j, xj, c[j], acolj, hcolj))
   end
   # update c0
   c0ps = c0 + c0_offset
@@ -66,5 +88,14 @@ function remove_ifix!(
 end
 
 function postsolve!(pt::OutputPoint{T, S}, operation::RemoveIfix{T, S}) where {T, S}
-  pt.x[operation.j] = operation.xj
+  j = operation.j
+  pt.x[j] = operation.xj
+  ATyj = @views dot(operation.acolj.nzval, pt.y[operation.acolj.nzind])
+  Hxj = @views dot(operation.hcolj.nzval, pt.x[operation.hcolj.nzind])
+  s = operation.cj + Hxj - ATyj
+  if s > zero(T)
+    pt.s_l[j] = s
+  else
+    pt.s_u[j] = -s
+  end
 end

src/presolve/singleton_rows.jl

@@ -37,7 +37,6 @@ function singleton_rows!(
     end
     Ax = rowi.nzval[k]
 
-    col_cnt[j] -= 1
     if Ax > zero(T)
       lvar2 = lcon[i] / Ax
       uvar2 = ucon[i] / Ax
@@ -47,32 +46,33 @@ function singleton_rows!(
     else
       error("remove explicit zeros in A")
     end
-    if lvar[j] < lvar2
+    if lvar[j] ≤ lvar2
       lvar[j] = lvar2
       tightened_lvar = true
     end
-    if uvar[j] > uvar2
+    if uvar[j] ≥ uvar2
       uvar[j] = uvar2
       tightened_uvar = true
     end
 
     row_cnt[i] = -1
     kept_rows[i] = false
+    col_cnt[j] -= 1
     push!(operations, SingletonRow{T, S}(i, j, Ax, tightened_lvar, tightened_uvar))
   end
   return singl_row_pass
 end
 
 function postsolve!(pt::OutputPoint{T, S}, operation::SingletonRow{T, S}) where {T, S}
   i, j = operation.i, operation.j
-  dual_slack = -pt.s_l[j] + pt.s_u[j]
   Aij = operation.Aij
-  if operation.tightened_lvar
-    pt.y[i] = dual_slack / Aij
+  pt.y[i] = zero(T)
+  if operation.tightened_lvar 
+    pt.y[i] += pt.s_l[j] / Aij 
     pt.s_l[j] = zero(T)
   end
   if operation.tightened_uvar
-    pt.y[i] = dual_slack / Aij
+    pt.y[i] -= pt.s_u[j] / Aij
     pt.s_u[j] = zero(T)
   end
   if !operation.tightened_lvar && !operation.tightened_uvar

test/test_presolve.jl

@@ -18,7 +18,7 @@
     SparseMatrixCOO(tril(H)),
     A = SparseMatrixCOO(A),
     lcon = [-3.0; -4.0],
-    ucon = [-2.0; Inf],
+    ucon = [20.0; Inf],
     lvar = l,
     uvar = u,
     c0 = 0.0,
@@ -83,7 +83,7 @@ end
     0.0 2.0 1.0
   ]
   b = [0.0; 0.0; 0.0; 4.0; 0.0; 3]
-  l = [0.0; 0.0; 0]
+  l = [-2.0; -4.0; 0]
   u = [Inf; 2.0; Inf]
   T = eltype(c)
   qp = QuadraticModel(
@@ -170,7 +170,7 @@ end
   s_l = sparse([3.0; 4.0; 2.0])
   s_u = sparse([0.0; 3.0; 0.0])
   pt_out = postsolve(qp, psqp, x_in, y_in, s_l, s_u)
-  @test pt_out.y == [2.0; 0.0; 0.0; 0.0; 0.0; 4.0]
+  @test pt_out.y == [2.0; 0.0; 0.0; 1.0; 0.0; 4.0]
 end
 
 @testset "presolve singleton rows and ifix" begin