QP solver

src/IPM/HSD/HSD.jl

@@ -6,7 +6,7 @@ Solver for the homogeneous self-dual algorithm.
 mutable struct HSD{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T}
 
     # Problem data, in standard form
-    dat::IPMData{T, Tv, Tb, Ta}
+    dat::LPData{T, Tv, Tb, Ta}
 
     # =================
     #   Book-keeping
@@ -32,9 +32,9 @@ mutable struct HSD{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T}
     regG::T   # gap regularization
 
     function HSD(
-        dat::IPMData{T, Tv, Tb, Ta}, kkt_options::KKTOptions{T}
+        dat::LPData{T, Tv, Tb, Ta}, kkt_options::KKTOptions{T}
     ) where{T, Tv<:AbstractVector{T}, Tb<:AbstractVector{Bool}, Ta<:AbstractMatrix{T}}
-        
+
         m, n = dat.nrow, dat.ncol
         p = sum(dat.lflag) + sum(dat.uflag)
 
@@ -161,15 +161,15 @@ function update_solver_status!(hsd::HSD{T}, ϵp::T, ϵd::T, ϵg::T, ϵi::T) wher
     else
         hsd.dual_status = Sln_Unknown
     end
-    
+
     # Check for optimal solution
     if ρp <= ϵp && ρd <= ϵd && ρg <= ϵg
         hsd.primal_status = Sln_Optimal
         hsd.dual_status   = Sln_Optimal
         hsd.solver_status = Trm_Optimal
         return nothing
     end
-    
+
     # Check for infeasibility certificates
     if max(
         norm(dat.A * pt.x, Inf),
@@ -244,7 +244,7 @@ function ipm_optimize!(hsd::HSD{T}, params::IPMOptions{T}) where{T}
     hsd.pt.y   .= zero(T)
     hsd.pt.zl  .= one(T) .* dat.lflag
     hsd.pt.zu  .= one(T) .* dat.uflag
-    
+
     hsd.pt.τ   = one(T)
     hsd.pt.κ   = one(T)
 
@@ -268,12 +268,12 @@ function ipm_optimize!(hsd::HSD{T}, params::IPMOptions{T}) where{T}
         if params.OutputLevel > 0
             # Display log
             @printf "%4d" hsd.niter
-            
+
             # Objectives
             ϵ = dat.objsense ? one(T) : -one(T)
             @printf "  %+14.7e" ϵ * hsd.primal_objective
             @printf "  %+14.7e" ϵ * hsd.dual_objective
-            
+
             # Residuals
             @printf "  %8.2e" max(hsd.res.rp_nrm, hsd.res.ru_nrm)
             @printf " %8.2e" hsd.res.rd_nrm
@@ -306,14 +306,14 @@ function ipm_optimize!(hsd::HSD{T}, params::IPMOptions{T}) where{T}
             || hsd.solver_status == Trm_DualInfeasible
         )
             break
-        elseif hsd.niter >= params.IterationsLimit 
+        elseif hsd.niter >= params.IterationsLimit
             hsd.solver_status = Trm_IterationLimit
             break
         elseif ttot >= params.TimeLimit
             hsd.solver_status = Trm_TimeLimit
             break
         end
-        
+
 
         # TODO: step
         # For now, include the factorization in the step function
@@ -325,7 +325,7 @@ function ipm_optimize!(hsd::HSD{T}, params::IPMOptions{T}) where{T}
             if isa(err, PosDefException) || isa(err, SingularException)
                 # Numerical trouble while computing the factorization
                 hsd.solver_status = Trm_NumericalProblem
-    
+
             elseif isa(err, OutOfMemoryError)
                 # Out of memory
                 hsd.solver_status = Trm_MemoryLimit
@@ -349,4 +349,4 @@ function ipm_optimize!(hsd::HSD{T}, params::IPMOptions{T}) where{T}
 
     return nothing
 
-end
+end

src/IPM/MPC/MPC.jl

@@ -3,10 +3,10 @@
 
 Implements Mehrotra's Predictor-Corrector interior-point algorithm.
 """
-mutable struct MPC{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T}
+mutable struct MPC{T, Tv, Tb, Ta, Tq, Tk} <: AbstractIPMOptimizer{T}
 
     # Problem data, in standard form
-    dat::IPMData{T, Tv, Tb, Ta}
+    dat::IPMData{T, Tv, Tb, Ta, Tq}
 
     # =================
     #   Book-keeping
@@ -48,8 +48,8 @@ mutable struct MPC{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T}
     regD::Tv  # Dual regularization
 
     function MPC(
-        dat::IPMData{T, Tv, Tb, Ta}, kkt_options::KKTOptions{T}
-    ) where{T, Tv<:AbstractVector{T}, Tb<:AbstractVector{Bool}, Ta<:AbstractMatrix{T}}
+        dat::IPMData{T, Tv, Tb, Ta, Tq}, kkt_options::KKTOptions{T}
+    ) where{T, Tv<:AbstractVector{T}, Tb<:AbstractVector{Bool}, Ta<:AbstractMatrix{T}, Tq<:AbstractMatrix{T}}
 
         m, n = dat.nrow, dat.ncol
         p = sum(dat.lflag) + sum(dat.uflag)
@@ -76,10 +76,10 @@ mutable struct MPC{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T}
         regP = tones(Tv, n)
         regD = tones(Tv, m)
 
-        kkt = KKT.setup(kkt_options.Factory.T, dat.A; kkt_options.Factory.options...)
+        kkt = KKT.setup(kkt_options.Factory.T, dat.A, dat.Q; kkt_options.Factory.options...)
         Tk = typeof(kkt)
 
-        return new{T, Tv, Tb, Ta, Tk}(dat,
+        return new{T, Tv, Tb, Ta, Tq, Tk}(dat,
             0, Trm_Unknown, Sln_Unknown, Sln_Unknown,
             T(Inf), T(-Inf),
             TimerOutput(),
@@ -103,6 +103,10 @@ function compute_residuals!(mpc::MPC{T}) where{T}
     pt, res = mpc.pt, mpc.res
     dat = mpc.dat
 
+
+    Qx = zeros(T, pt.n)
+    mul!(Qx, dat.Q, pt.x, true, true)
+
     # Primal residual
     # rp = b - A*x
     res.rp .= dat.b
@@ -123,6 +127,7 @@ function compute_residuals!(mpc::MPC{T}) where{T}
     # Dual residual
     # rd = c - (A'y + zl - zu)
     res.rd .= dat.c
+    axpy!(one(T), Qx, res.rd)
     mul!(res.rd, transpose(dat.A), pt.y, -one(T), one(T))
     @. res.rd += pt.zu .* dat.uflag - pt.zl .* dat.lflag
 
@@ -133,9 +138,11 @@ function compute_residuals!(mpc::MPC{T}) where{T}
     res.rd_nrm = norm(res.rd, Inf)
 
     # Compute primal and dual bounds
-    mpc.primal_objective = dot(dat.c, pt.x) + dat.c0
+    xᵗQx = dot(pt.x, Qx) / 2
+    mpc.primal_objective = xᵗQx + dot(dat.c, pt.x) + dat.c0
     mpc.dual_objective   = (
-        dot(dat.b, pt.y)
+        - xᵗQx
+        + dot(dat.b, pt.y)
         + dot(dat.l .* dat.lflag, pt.zl)
         - dot(dat.u .* dat.uflag, pt.zu)
     ) + dat.c0
@@ -232,7 +239,7 @@ function ipm_optimize!(mpc::MPC{T}, params::IPMOptions{T}) where{T}
         @printf "\nOptimizer info (MPC)\n"
         @printf "Constraints  : %d\n" dat.nrow
         @printf "Variables    : %d\n" dat.ncol
-        bmin, bmax = extrema(dat.b)
+        bmin, bmax = length(dat.b) == 0 ? (zero(T), zero(T)) : extrema(dat.b)
         @printf "RHS          : [%+.2e, %+.2e]\n" bmin bmax
         lmin, lmax = extrema(dat.l .* dat.lflag)
         @printf "Lower bounds : [%+.2e, %+.2e]\n" lmin lmax

src/IPM/ipmdata.jl

@@ -1,31 +1,30 @@
 """
-    IPMData{T, Tv, Ta}
+    IPMData{T, Tv, Tb, Ta, Tq}
 
 Holds data about an interior point method.
 
 The problem is represented as
 ```
-min   c'x + c0
+min   ¹/₂ x'Qx + c'x + c0
 s.t.  A x = b
       l ≤ x ≤ u
 ```
 where `l`, `u` may take infinite values.
 """
-struct IPMData{T, Tv, Tb, Ta}
+struct IPMData{T, Tv, Tb, Ta, Tq}
 
     # Problem size
     nrow::Int
     ncol::Int
 
     # Objective
     objsense::Bool  # min (true) or max (false)
-    c0::T
+    Q::Tq
     c::Tv
+    c0::T
 
-    # Constraint matrix
+    # Linear constraints
     A::Ta
-
-    # RHS
     b::Tv
 
     # Variable bounds (may contain infinite values)
@@ -39,22 +38,30 @@ struct IPMData{T, Tv, Tb, Ta}
     uflag::Tb
 
     function IPMData(
-        A::Ta, b::Tv, objsense::Bool, c::Tv, c0::T, l::Tv, u::Tv
-    ) where{T, Tv<:AbstractVector{T}, Ta<:AbstractMatrix{T}}
+        A::Ta, b::Tv, objsense::Bool, Q::Tq, c::Tv, c0::T, l::Tv, u::Tv
+    ) where{T, Tv<:AbstractVector{T}, Ta<:AbstractMatrix{T}, Tq<:AbstractMatrix{T}}
         nrow, ncol = size(A)
 
         lflag = isfinite.(l)
         uflag = isfinite.(u)
         Tb = typeof(lflag)
 
-        return new{T, Tv, Tb, Ta}(
+        return new{T, Tv, Tb, Ta, Tq}(
             nrow, ncol,
-            objsense, c0, c,
+            objsense, Q, c, c0,
             A, b, l, u, lflag, uflag
         )
     end
 end
 
+const LPData{T, Tv, Tb, Ta} = IPMData{T, Tv, Tb, Ta, ZeroMatrix{T}}
+
+function LPData(A, b, objsense, c, c0, l, u)
+    m, n = size(A)
+    Q = ZeroMatrix{eltype(A)}(n, n)
+    return IPMData(A, b, objsense, Q, c, c0, l, u)
+end
+
 # TODO: extract IPM data from presolved problem
 """
     IPMData(pb::ProblemData, options::MatrixOptions)
@@ -107,7 +114,7 @@ function IPMData(pb::ProblemData{T}, mfact::Factory) where{T}
             # lb <= a'x <= ub
             # Two options:
             # --> a'x + s = ub, 0 <= s <= ub - lb
-            # --> a'x - s = lb, 0 <= s <= ub - lb 
+            # --> a'x - s = lb, 0 <= s <= ub - lb
             push!(sind, i)
             push!(sval, one(T))
             push!(lslack, zero(T))
@@ -168,6 +175,6 @@ function IPMData(pb::ProblemData{T}, mfact::Factory) where{T}
     # Variable bounds
     l = [pb.lvar; lslack]
     u = [pb.uvar; uslack]
-    
-    return IPMData(A, b, pb.objsense, c, c0, l, u)
-end
+
+    return LPData(A, b, pb.objsense, c, c0, l, u)
+end

src/KKT/KKT.jl

@@ -3,6 +3,7 @@ module KKT
 using LinearAlgebra
 
 import ..Tulip.Factory
+using ..TLPLinearAlgebra
 
 export AbstractKKTSolver, KKTOptions
 
@@ -37,9 +38,9 @@ function setup(Ts::Type{<:AbstractKKTSolver}, args...; kwargs...)
     return Ts(args...; kwargs...)
 end
 
-# 
+#
 # Specialized implementations should extend the functions below
-# 
+#
 
 """
     update!(kkt, θinv, regP, regD)
@@ -122,4 +123,4 @@ function default_factory(::Type{T}) where{T}
     end
 end
 
-end  # module
+end  # module