Feasibility presolve (#118)

* feasibility presolve

* primal cstr update

Author: geoffroyleconte

src/presolve/empty_rows.jl

@@ -20,4 +20,5 @@ end
 
 function postsolve!(sol::QMSolution, operation::EmptyRow, psd::PresolvedData)
   psd.kept_rows[operation.i] = true
+  sol.y[operation.i] = 0
 end

src/presolve/postsolve_utils.jl

@@ -1,5 +1,5 @@
 function restore_x!(kept_cols, x_in::S, x::S, nvar) where {S}
-  # put x and xps inside xout according to kept_cols
+  # put x_in inside x according to kept_cols
   cx = 0
   for i = 1:nvar
     if kept_cols[i]
@@ -66,3 +66,16 @@ function restore_s!(
   restore_x!(kept_cols, s_l_in, s_l, nvar)
   restore_x!(kept_cols, s_u_in, s_u, nvar)
 end
+
+function add_Hx!(z::Vector{T}, hcols::Vector{Col{T}}, kept_cols::Vector{Bool}, x::Vector{T}) where {T}
+  for j in 1:length(hcols)
+    hcolj = hcols[j]
+    Hxj = zero(T)
+    for (i, hij) in zip(hcolj.nzind, hcolj.nzval)
+      kept_cols[i] || continue
+      Hxj += hij * x[i]
+    end
+    z[j] += Hxj
+  end
+  return z
+end

src/presolve/presolve.jl

@@ -113,6 +113,8 @@ end
 
 mutable struct PresolvedData{T, S}
   xps::S
+  c::S # copy of c
+  z::S # c + Qx for postsolve PrimalConstraints
   arows::Vector{Row{T}}
   acols::Vector{Col{T}}
   hcols::Vector{Col{T}}
@@ -289,6 +291,26 @@ function presolve(
     psdata.c0 = qmp.c0
   end
 
+  # c_ps is the new c padded with zeros to have the same size as the init c 
+  c_ps = fill!(S(undef, nvar), zero(T))
+  restore_x!(qmp.kept_cols, qmp.c, c_ps, nvar)
+  # group all presolve info into a struct
+  psd = PresolvedData{T, S}(
+    qmp.xps,
+    c_ps,
+    copy(c_ps),
+    qmp.arows,
+    qmp.acols,
+    qmp.hcols,
+    qmp.kept_rows,
+    qmp.kept_cols,
+    nvarps,
+    nconps,
+    nvar,
+    ncon,
+    operations,
+  )
+
   # form meta
   nnzh = length(psdata.H.vals)
   if !(nnzh == length(psdata.H.rows) == length(psdata.H.cols))
@@ -318,16 +340,20 @@ function presolve(
       solver_specific = Dict(:presolvedQM => nothing),
     )
   elseif nvarps == 0
+    sol_in = QMSolution(S(undef, 0), fill!(S(undef, nconps), zero(T)), S(undef, 0), S(undef, 0))
+    sol = postsolve(psd, sol_in)
     feasible = all(qm.meta.lcon .<= qm.data.A * qmp.xps .<= qm.meta.ucon)
-    s = qm.data.c .+ Symmetric(qm.data.H, :L) * qmp.xps
+    # s = qm.data.c .+ Symmetric(qm.data.H, :L) * qmp.xps
     return GenericExecutionStats(
       qm,
       status = feasible ? :first_order : :infeasible,
-      solution = qmp.xps,
-      objective = obj(qm, qmp.xps),
-      multipliers = zeros(T, ncon),
-      multipliers_L = max.(s, zero(T)),
-      multipliers_U = max.(.-s, zero(T)),
+      solution = sol.x,
+      objective = obj(qm, sol.x),
+      primal_feas = feasible ? zero(T) : T(Inf),
+      dual_feas = feasible ? zero(T) : T(Inf),
+      multipliers = sol.y,
+      multipliers_L = sol.s_l,
+      multipliers_U = sol.s_u,
       iter = 0,
       elapsed_time = time() - start_time,
       solver_specific = Dict(:presolvedQM => nothing, :psoperations => operations),
@@ -349,19 +375,6 @@ function presolve(
       minimize = qm.meta.minimize,
       kwargs...,
     )
-    psd = PresolvedData{T, S}(
-      qmp.xps,
-      qmp.arows,
-      qmp.acols,
-      qmp.hcols,
-      qmp.kept_rows,
-      qmp.kept_cols,
-      nvarps,
-      nconps,
-      nvar,
-      ncon,
-      operations,
-    )
     ps = PresolvedQuadraticModel(psmeta, Counters(), psdata, psd)
     return GenericExecutionStats(
       ps,
@@ -374,13 +387,11 @@ function presolve(
 end
 
 function postsolve!(
-  qm::QuadraticModel{T, S},
-  psqm::PresolvedQuadraticModel{T, S},
+  psd::PresolvedData{T, S},
   sol::QMSolution{S},
   sol_in::QMSolution{S},
 ) where {T, S}
   x_in, y_in, s_l_in, s_u_in = sol_in.x, sol_in.y, sol_in.s_l, sol_in.s_u
-  psd = psqm.psd
   n_operations = length(psd.operations)
   nvar = psd.nvar
   @assert nvar == length(sol.x)
@@ -389,13 +400,34 @@ function postsolve!(
   @assert ncon == length(sol.y)
   restore_y!(psd.kept_rows, y_in, sol.y, ncon)
   restore_s!(sol.s_l, sol.s_u, s_l_in, s_u_in, psd.kept_cols)
+  # add_Hx!(psd.z, psd.hcols, psd.kept_cols) # z = c + Hx
 
   for i = n_operations:-1:1
     operation_i = psd.operations[i]
     postsolve!(sol, operation_i, psd)
   end
 end
 
+postsolve!(
+  psqm::PresolvedQuadraticModel{T, S},
+  sol::QMSolution{S},
+  sol_in::QMSolution{S},
+) where {T, S} = postsolve!(psqm.psd, sol, sol_in)
+
+function postsolve(
+  psd::PresolvedData{T, S},
+  sol_in::QMSolution{S},
+) where {T, S}
+  x = fill!(S(undef, psd.nvar), zero(T))
+  y = fill!(S(undef, psd.ncon), zero(T))
+  s_l = fill!(S(undef, psd.nvar), zero(T))
+  s_u = fill!(S(undef, psd.nvar), zero(T))
+
+  sol = QMSolution(x, y, s_l, s_u)
+  postsolve!(psd, sol, sol_in)
+  return sol
+end
+
 """
     sol = postsolve(qm::QuadraticModel{T, S}, psqm::PresolvedQuadraticModel{T, S}, 
                     sol_in::QMSolution{S}) where {T, S}
@@ -404,17 +436,8 @@ Retrieve the solution `sol = (x, y, s_l, s_u)` of the original QP `qm` given the
 `sol_in` of type [`QMSolution`](@ref).
 `sol_in.s_l` and `sol_in.s_u` can be sparse or dense vectors, but the output `sol.s_l` and `sol.s_u` are dense vectors. 
 """
-function postsolve(
+postsolve(
   qm::QuadraticModel{T, S},
   psqm::PresolvedQuadraticModel{T, S},
   sol_in::QMSolution{S},
-) where {T, S}
-  x = fill!(S(undef, psqm.psd.nvar), zero(T))
-  y = fill!(S(undef, psqm.psd.ncon), zero(T))
-  s_l = fill!(S(undef, psqm.psd.nvar), zero(T))
-  s_u = fill!(S(undef, psqm.psd.nvar), zero(T))
-
-  sol = QMSolution(x, y, s_l, s_u)
-  postsolve!(qm, psqm, sol, sol_in)
-  return sol
-end
+) where {T, S} = postsolve(psqm.psd, sol_in)

src/presolve/primal_constraints.jl

@@ -1,3 +1,8 @@
+struct PrimalConstraint{T, S} <: PresolveOperation{T, S}
+  i::Int
+  forced_lcon::Bool # row i forced to lcon[i]
+end
+
 function primal_constraints!(
   qmp::QuadraticModelPresolveData{T, S},
   operations::Vector{PresolveOperation{T, S}},
@@ -28,6 +33,7 @@ function primal_constraints!(
           uvar[j] = lvar[j]
         end
       end
+      push!(operations, PrimalConstraint{T, S}(i, true))
     elseif lconi2 == ucon[i]
       for (j, aij) in zip(rowi.nzind, rowi.nzval)
         kept_cols[j] || continue
@@ -37,6 +43,7 @@ function primal_constraints!(
           lvar[j] = uvar[j]
         end
       end
+      push!(operations, PrimalConstraint{T, S}(i, false))
     elseif lcon[i] < lconi2 && uconi2 < ucon[i]
       kept_rows[i] = false
       row_cnt[i] = -1
@@ -52,3 +59,55 @@ function primal_constraints!(
     end
   end
 end
+
+function postsolve!(
+  sol::QMSolution,
+  operation::PrimalConstraint{T, S},
+  psd::PresolvedData{T, S},
+) where {T, S}
+  i = operation.i
+  forced_lcon = operation.forced_lcon
+  x = sol.x
+  y = sol.y
+  s_l = sol.s_l
+  s_u = sol.s_u
+  arowi = psd.arows[i]
+  acols = psd.acols
+  kept_cols = psd.kept_cols
+  z = psd.z
+  z .= psd.c
+  n = length(z)
+
+  add_Hx!(z, psd.hcols, kept_cols, x) # z = c + Hx
+  for l = 1:n
+    kept_cols[l] || continue
+    for (k, akl) in zip(acols[l].nzind, acols[l].nzval) 
+      (psd.kept_rows[i] && k != i) || continue
+      z[l] -= akl * y[k]
+    end
+  end
+
+  if forced_lcon
+    yi = T(-Inf) 
+    for (l, ail) in zip(arowi.nzind, arowi.nzval)
+      (kept_cols[l]) || continue
+      trial = z[l] / ail
+      (trial > yi) && (yi = trial)
+    end
+  else
+    yi = T(Inf) 
+    for (l, ail) in zip(arowi.nzind, arowi.nzval)
+      (kept_cols[l]) || continue
+      trial = z[l] / ail
+      (trial < yi) && (yi = trial)
+    end
+  end
+  y[i] = yi
+
+  for (l, ail) in zip(arowi.nzind, arowi.nzval)
+    (kept_cols[l]) || continue
+    s = z[l] - ail * yi
+    s_l[l] = s > 0 ? s : zero(T)
+    s_u[l] = s < 0 ? -s : zero(T)
+  end
+end

src/presolve/remove_ifix.jl

@@ -68,7 +68,10 @@ function postsolve!(
   hcolj = psd.hcols[j]
   x = sol.x
   x[j] = operation.xj
+  cj = operation.cj
   y = sol.y
+  c = psd.c
+  c[j] = cj
 
   ATyj = zero(T)
   for (i, aij) in zip(acolj.nzind, acolj.nzval)
@@ -80,12 +83,10 @@ function postsolve!(
   for (i, hij) in zip(hcolj.nzind, hcolj.nzval)
     psd.kept_cols[i] || continue
     Hxj += hij * x[i]
+    (i != j) && (c[i] -= hij * x[j])
   end
 
-  s = operation.cj + Hxj - ATyj
-  if s > zero(T)
-    sol.s_l[j] = s
-  else
-    sol.s_u[j] = -s
-  end
+  s = cj + Hxj - ATyj
+  sol.s_l[j] = s > zero(T) ? s : zero(T)
+  sol.s_u[j] = s < zero(T) ? -s : zero(T)
 end

test/test_presolve.jl

@@ -165,6 +165,36 @@ end
   @test sol.y == [2.0; 0.0; 0.0; 1.0; 0.0; 4.0]
 end
 
+@testset "presolve solves problem" begin
+  Q = [6. 2. 1.
+       2. 5. 2.
+       1. 2. 4.]
+  c = [-8.; -3; -3]
+  A = [1. 0. 1.
+       0. 2. 1.]
+  b = [0.; 3]
+  l = [0.;0;0]
+  u = [Inf; Inf; Inf]
+  qp = QuadraticModel(
+    c,
+    SparseMatrixCOO(tril(Q)),
+    A=SparseMatrixCOO(A),
+    lcon=b,
+    ucon=b,
+    lvar=l,
+    uvar=u,
+    c0=0.,
+    name="QM"
+  )
+  statsps = presolve(qp)
+  psqp = statsps.solver_specific[:presolvedQM]
+  atol = eps()
+  @test isapprox(statsps.solution, [0.0; 1.5; 0.0], atol = atol)
+  @test isapprox(statsps.objective, 1.125, atol = atol)
+  @test isapprox(statsps.multipliers_L[2], 0.0, atol = atol)
+  @test isapprox(statsps.multipliers_U, [0.0; 0.0; 0.0], atol = atol)
+end
+
 @testset "presolve singleton rows and ifix" begin
   H = [
     6.0 2.0 1.0
@@ -201,7 +231,6 @@ end
 
   @test statsps.status == :first_order
   @test statsps.solution ≈ [4.0 / 3.0; 7.0 / 6.0; 2.0 / 3.0]
-  @test statsps.multipliers == zeros(length(b))
 end
 
 @testset "presolve free col singleton" begin