From 295a4372ee8b362b2ddd67c20b8bb7b4da6ca9fb Mon Sep 17 00:00:00 2001 From: Geoffroy Leconte Date: Sun, 3 Dec 2023 15:54:19 -0500 Subject: [PATCH 1/2] feasibility presolve --- src/presolve/presolve.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/presolve/presolve.jl b/src/presolve/presolve.jl index b8b7fd0..c9858aa 100644 --- a/src/presolve/presolve.jl +++ b/src/presolve/presolve.jl @@ -325,6 +325,8 @@ function presolve( status = feasible ? :first_order : :infeasible, solution = qmp.xps, objective = obj(qm, qmp.xps), + primal_feas = feasible ? zero(T) : T(Inf), + dual_feas = feasible ? zero(T) : T(Inf), multipliers = zeros(T, ncon), multipliers_L = max.(s, zero(T)), multipliers_U = max.(.-s, zero(T)), From 59154fe072ac19fe24e7fec361e776c513faad6c Mon Sep 17 00:00:00 2001 From: Geoffroy Leconte Date: Mon, 4 Dec 2023 13:59:52 -0500 Subject: [PATCH 2/2] primal cstr update --- src/presolve/empty_rows.jl | 1 + src/presolve/postsolve_utils.jl | 15 +++++- src/presolve/presolve.jl | 87 ++++++++++++++++++------------ src/presolve/primal_constraints.jl | 59 ++++++++++++++++++++ src/presolve/remove_ifix.jl | 13 ++--- test/test_presolve.jl | 31 ++++++++++- 6 files changed, 165 insertions(+), 41 deletions(-) diff --git a/src/presolve/empty_rows.jl b/src/presolve/empty_rows.jl index c497e8a..a883960 100644 --- a/src/presolve/empty_rows.jl +++ b/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 diff --git a/src/presolve/postsolve_utils.jl b/src/presolve/postsolve_utils.jl index 65610ae..30e3166 100644 --- a/src/presolve/postsolve_utils.jl +++ b/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 \ No newline at end of file diff --git a/src/presolve/presolve.jl b/src/presolve/presolve.jl index c9858aa..45e6a65 100644 --- a/src/presolve/presolve.jl +++ b/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,18 +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), + solution = sol.x, + objective = obj(qm, sol.x), primal_feas = feasible ? zero(T) : T(Inf), dual_feas = feasible ? zero(T) : T(Inf), - multipliers = zeros(T, ncon), - multipliers_L = max.(s, zero(T)), - multipliers_U = max.(.-s, zero(T)), + 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), @@ -351,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, @@ -376,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) @@ -391,6 +400,7 @@ 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] @@ -398,6 +408,26 @@ function postsolve!( 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} @@ -406,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) diff --git a/src/presolve/primal_constraints.jl b/src/presolve/primal_constraints.jl index 266a5fd..906f8a1 100644 --- a/src/presolve/primal_constraints.jl +++ b/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 diff --git a/src/presolve/remove_ifix.jl b/src/presolve/remove_ifix.jl index 344824f..0591bf7 100644 --- a/src/presolve/remove_ifix.jl +++ b/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 diff --git a/test/test_presolve.jl b/test/test_presolve.jl index 3da546f..62ef503 100644 --- a/test/test_presolve.jl +++ b/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