From ae6812c1732f5b0b96ce88eaa62d528e892b8cd9 Mon Sep 17 00:00:00 2001 From: mtanneau Date: Sat, 26 Dec 2020 12:34:29 -0500 Subject: [PATCH 1/2] Initial QP implementation (MPC only) --- src/IPM/HSD/HSD.jl | 24 ++++++++-------- src/IPM/MPC/MPC.jl | 25 +++++++++++------ src/IPM/ipmdata.jl | 37 ++++++++++++++---------- src/KKT/KKT.jl | 7 +++-- src/KKT/cholmod.jl | 45 +++++++++++++++++++++++------- src/LinearAlgebra/LinearAlgebra.jl | 5 +++- src/LinearAlgebra/zero.jl | 20 +++++++++++++ src/Tulip.jl | 2 +- test/IPM/HSD.jl | 10 +++---- test/IPM/MPC.jl | 10 +++---- 10 files changed, 124 insertions(+), 61 deletions(-) create mode 100644 src/LinearAlgebra/zero.jl diff --git a/src/IPM/HSD/HSD.jl b/src/IPM/HSD/HSD.jl index 1c073932..a65aa72e 100644 --- a/src/IPM/HSD/HSD.jl +++ b/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,7 +161,7 @@ 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 @@ -169,7 +169,7 @@ function update_solver_status!(hsd::HSD{T}, ϵp::T, ϵd::T, ϵg::T, ϵi::T) wher 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 \ No newline at end of file +end diff --git a/src/IPM/MPC/MPC.jl b/src/IPM/MPC/MPC.jl index eec253e6..a3f72cf5 100644 --- a/src/IPM/MPC/MPC.jl +++ b/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 diff --git a/src/IPM/ipmdata.jl b/src/IPM/ipmdata.jl index 94999778..bed13568 100644 --- a/src/IPM/ipmdata.jl +++ b/src/IPM/ipmdata.jl @@ -1,17 +1,17 @@ """ - 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 @@ -19,13 +19,12 @@ struct IPMData{T, Tv, Tb, Ta} # 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 \ No newline at end of file + + return LPData(A, b, pb.objsense, c, c0, l, u) +end diff --git a/src/KKT/KKT.jl b/src/KKT/KKT.jl index f8ea8044..d48db63d 100644 --- a/src/KKT/KKT.jl +++ b/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 \ No newline at end of file +end # module diff --git a/src/KKT/cholmod.jl b/src/KKT/cholmod.jl index 0a3e9957..1d1482a5 100644 --- a/src/KKT/cholmod.jl +++ b/src/KKT/cholmod.jl @@ -24,14 +24,14 @@ model.params.KKT.Factory = Tulip.Factory(CholmodSolver, normal_equations=true) abstract type CholmodSolver <: AbstractKKTSolver{Float64} end function CholmodSolver( - A::AbstractMatrix{Float64}; + args...; normal_equations::Bool=false, kwargs... ) if normal_equations - return CholmodSPD(A; kwargs...) + return CholmodSPD(args...; kwargs...) else - return CholmodSQD(A; kwargs...) + return CholmodSQD(args...; kwargs...) end end @@ -55,6 +55,8 @@ mutable struct CholmodSQD <: CholmodSolver # and add flag (default to false) to know whether user ordering should be used # Problem data + Q::SparseMatrixCSC{Float64} + dQ::Vector{Float64} # diagonal of Q A::SparseMatrixCSC{Float64} θ::Vector{Float64} regP::Vector{Float64} # primal-dual regularization @@ -71,17 +73,40 @@ mutable struct CholmodSQD <: CholmodSolver m, n = size(A) θ = ones(Float64, n) + dQ = zeros(Float64, n) + S = [ - spdiagm(0 => -θ) A'; + -spdiagm(0 => θ) A'; spzeros(Float64, m, n) spdiagm(0 => ones(m)) ] # TODO: Symbolic factorization only F = ldlt(Symmetric(S)) - return new(m, n, A, θ, ones(Float64, n), ones(Float64, m), S, F) + return new(m, n, spzeros(Float64, n, n), dQ, A, θ, ones(Float64, n), ones(Float64, m), S, F) end + function CholmodSQD(A::AbstractMatrix{Float64}, Q::SparseMatrixCSC{Float64}) + m, n = size(A) + θ = ones(Float64, n) + + dQ = Vector(diag(Q, 0)) + + # TODO: improve this + U = sparse(UpperTriangular(Q)) + S = [ + -U - spdiagm(0 => θ) A'; + spzeros(Float64, m, n) spdiagm(0 => ones(m)) + ] + + # TODO: Symbolic factorization only + F = ldlt(Symmetric(S)) + + return new(m, n, Q, dQ, A, θ, ones(Float64, n), ones(Float64, m), S, F) + end + + CholmodSQD(A::AbstractMatrix{Float64}, ::ZeroMatrix{Float64}) = CholmodSQD(A, spzeros(Float64, size(A, 2), size(A, 2))) + end backend(::CholmodSQD) = "CHOLMOD" @@ -122,7 +147,7 @@ function update!( # Update S. S is stored as upper-triangular and only its diagonal changes. @inbounds for j in 1:kkt.n k = kkt.S.colptr[1+j] - 1 - kkt.S.nzval[k] = -kkt.θ[j] - regP[j] + kkt.S.nzval[k] = -kkt.dQ[j] - kkt.θ[j] - regP[j] end @inbounds for i in 1:kkt.m k = kkt.S.colptr[1+kkt.n+i] - 1 @@ -147,7 +172,7 @@ function solve!( ξp::Vector{Float64}, ξd::Vector{Float64} ) m, n = kkt.m, kkt.n - + # Set-up right-hand side ξ = [ξd; ξp] @@ -244,7 +269,7 @@ function update!( regP::AbstractVector{Float64}, regD::AbstractVector{Float64} ) - + # Sanity checks length(θ) == kkt.n || throw(DimensionMismatch( "θ has length $(length(θ)) but linear solver is for n=$(kkt.n)." @@ -287,7 +312,7 @@ function solve!( d = one(Float64) ./ (kkt.θ .+ kkt.regP) D = Diagonal(d) - + # Set-up right-hand side ξ_ = ξp .+ kkt.A * (D * ξd) @@ -305,6 +330,6 @@ function solve!( # resD = - D \ dx + kkt.A' * dy - ξd # println("\n|resP| = $(norm(resP, Inf))\n|resD| = $(norm(resD, Inf))") - + return nothing end diff --git a/src/LinearAlgebra/LinearAlgebra.jl b/src/LinearAlgebra/LinearAlgebra.jl index e10ab928..40ba4544 100644 --- a/src/LinearAlgebra/LinearAlgebra.jl +++ b/src/LinearAlgebra/LinearAlgebra.jl @@ -3,6 +3,7 @@ module TLPLinearAlgebra using LinearAlgebra using SparseArrays export construct_matrix +export ZeroMatrix import ..Tulip.Factory @@ -31,4 +32,6 @@ construct_matrix( aI::Vector{Int}, aJ::Vector{Int}, aV::Vector{T} ) where{T} = sparse(aI, aJ, aV, m, n) -end # module \ No newline at end of file +include("zero.jl") + +end # module diff --git a/src/LinearAlgebra/zero.jl b/src/LinearAlgebra/zero.jl new file mode 100644 index 00000000..0517ea98 --- /dev/null +++ b/src/LinearAlgebra/zero.jl @@ -0,0 +1,20 @@ +import Base: + getindex, size +import LinearAlgebra.mul! + +""" + ZeroMatrix{T} <: AbstractMatrix{T} + +A type for representing matrices whose coefficients are all zero. +""" +struct ZeroMatrix{T} <: AbstractMatrix{T} + m::Int + n::Int + + ZeroMatrix{T}(m, n) where{T} = new{T}(m, n) +end + +size(A::ZeroMatrix) = (A.m, A.n) +getindex(::ZeroMatrix{T}, ::Integer, ::Integer) where{T} = zero(T) + +mul!(x::AbstractVector{T}, ::ZeroMatrix{T}, y::AbstractVector{T}, α::Number, β::Number) where{T} = rmul!(x, β) diff --git a/src/Tulip.jl b/src/Tulip.jl index 0fc89b56..a744cdee 100644 --- a/src/Tulip.jl +++ b/src/Tulip.jl @@ -13,7 +13,7 @@ include("utils.jl") # Linear algebra include("LinearAlgebra/LinearAlgebra.jl") -import .TLPLinearAlgebra.construct_matrix +using .TLPLinearAlgebra const TLA = TLPLinearAlgebra # KKT solvers diff --git a/test/IPM/HSD.jl b/test/IPM/HSD.jl index 062ec070..577f75ef 100644 --- a/test/IPM/HSD.jl +++ b/test/IPM/HSD.jl @@ -60,7 +60,7 @@ function run_tests_hsd(T::Type) c0 = zero(T) l = Vector{T}([0, 0]) u = Vector{T}([2, 2]) - dat = Tulip.IPMData(A, b, true, c, c0, l, u) + dat = Tulip.LPData(A, b, true, c, c0, l, u) hsd = TLP.HSD(dat, kkt_options) @@ -73,7 +73,7 @@ function run_tests_hsd(T::Type) hsd.pt.y .= T.([0, 1]) hsd.pt.zl .= T.([0, 0]) hsd.pt.zu .= T.([0, 0]) - + hsd.pt.τ = 1 hsd.pt.κ = 0 hsd.pt.μ = 0 @@ -83,12 +83,12 @@ function run_tests_hsd(T::Type) @testset "Residuals" begin @inferred TLP.compute_residuals!(hsd) TLP.compute_residuals!(hsd) - + @test isapprox(hsd.res.rp_nrm, zero(T); atol=ϵ, rtol=ϵ) @test isapprox(hsd.res.ru_nrm, zero(T); atol=ϵ, rtol=ϵ) @test isapprox(hsd.res.rd_nrm, zero(T); atol=ϵ, rtol=ϵ) @test isapprox(hsd.res.rg_nrm, zero(T); atol=ϵ, rtol=ϵ) - + end @testset "Convergence" begin @@ -110,4 +110,4 @@ end for T in TvTYPES @testset "$T" begin run_tests_hsd(T) end end -end \ No newline at end of file +end diff --git a/test/IPM/MPC.jl b/test/IPM/MPC.jl index 2818a271..a1d8a120 100644 --- a/test/IPM/MPC.jl +++ b/test/IPM/MPC.jl @@ -60,7 +60,7 @@ function run_tests_mpc(T::Type) c0 = zero(T) l = Vector{T}([0, 0]) u = Vector{T}([2, 2]) - dat = Tulip.IPMData(A, b, true, c, c0, l, u) + dat = Tulip.LPData(A, b, true, c, c0, l, u) ipm = TLP.MPC(dat, kkt_options) @@ -73,7 +73,7 @@ function run_tests_mpc(T::Type) ipm.pt.y .= T.([0, 1]) ipm.pt.zl .= T.([0, 0]) ipm.pt.zu .= T.([0, 0]) - + ipm.pt.τ = 1 ipm.pt.κ = 0 ipm.pt.μ = 0 @@ -83,12 +83,12 @@ function run_tests_mpc(T::Type) @testset "Residuals" begin @inferred TLP.compute_residuals!(ipm) TLP.compute_residuals!(ipm) - + @test isapprox(ipm.res.rp_nrm, zero(T); atol=ϵ, rtol=ϵ) @test isapprox(ipm.res.ru_nrm, zero(T); atol=ϵ, rtol=ϵ) @test isapprox(ipm.res.rd_nrm, zero(T); atol=ϵ, rtol=ϵ) @test isapprox(ipm.res.rg_nrm, zero(T); atol=ϵ, rtol=ϵ) - + end @testset "Convergence" begin @@ -110,4 +110,4 @@ end for T in TvTYPES @testset "$T" begin run_tests_mpc(T) end end -end \ No newline at end of file +end From 62b4cd3bb39be3cbfa82a3987513f352ff4089a6 Mon Sep 17 00:00:00 2001 From: mtanneau Date: Thu, 18 Mar 2021 22:41:27 -0400 Subject: [PATCH 2/2] Fix LDLF to support QP --- src/KKT/cholmod.jl | 15 +++++++-------- src/KKT/ldlfact.jl | 32 +++++++++++++++++++++++++------- 2 files changed, 32 insertions(+), 15 deletions(-) diff --git a/src/KKT/cholmod.jl b/src/KKT/cholmod.jl index 1d1482a5..90424b6f 100644 --- a/src/KKT/cholmod.jl +++ b/src/KKT/cholmod.jl @@ -55,21 +55,21 @@ mutable struct CholmodSQD <: CholmodSolver # and add flag (default to false) to know whether user ordering should be used # Problem data - Q::SparseMatrixCSC{Float64} + Q::SparseMatrixCSC{Float64,Int} dQ::Vector{Float64} # diagonal of Q - A::SparseMatrixCSC{Float64} + A::SparseMatrixCSC{Float64,Int} θ::Vector{Float64} regP::Vector{Float64} # primal-dual regularization regD::Vector{Float64} # dual regularization # Left-hand side matrix - S::SparseMatrixCSC{Float64, Int} # TODO: use `CHOLMOD.Sparse` instead + S::SparseMatrixCSC{Float64,Int} # TODO: use `CHOLMOD.Sparse` instead # Factorization F::CHOLMOD.Factor{Float64} # TODO: constructor with initial memory allocation - function CholmodSQD(A::AbstractMatrix{Float64}) + function CholmodSQD(A::SparseMatrixCSC{Float64,Int}) m, n = size(A) θ = ones(Float64, n) @@ -85,8 +85,7 @@ mutable struct CholmodSQD <: CholmodSolver return new(m, n, spzeros(Float64, n, n), dQ, A, θ, ones(Float64, n), ones(Float64, m), S, F) end - - function CholmodSQD(A::AbstractMatrix{Float64}, Q::SparseMatrixCSC{Float64}) + function CholmodSQD(A::SparseMatrixCSC{Float64,Int}, Q::SparseMatrixCSC{Float64,Int}) m, n = size(A) θ = ones(Float64, n) @@ -104,8 +103,8 @@ mutable struct CholmodSQD <: CholmodSolver return new(m, n, Q, dQ, A, θ, ones(Float64, n), ones(Float64, m), S, F) end - - CholmodSQD(A::AbstractMatrix{Float64}, ::ZeroMatrix{Float64}) = CholmodSQD(A, spzeros(Float64, size(A, 2), size(A, 2))) + CholmodSQD(A, ::ZeroMatrix) = CholmodSQD(A) + CholmodSQD(A) = CholmodSQD(sparse(A)) end diff --git a/src/KKT/ldlfact.jl b/src/KKT/ldlfact.jl index 56f55951..0b31d20f 100644 --- a/src/KKT/ldlfact.jl +++ b/src/KKT/ldlfact.jl @@ -24,7 +24,9 @@ mutable struct LDLFactSQD{T<:Real} <: AbstractKKTSolver{T} # and add flag (default to false) to know whether user ordering should be used # Problem data - A::SparseMatrixCSC{T, Int} + Q::SparseMatrixCSC{T,Int} + dQ::Vector{T} # diagonal of Q + A::SparseMatrixCSC{T,Int} θ::Vector{T} regP::Vector{T} # primal regularization regD::Vector{T} # dual regularization @@ -35,11 +37,11 @@ mutable struct LDLFactSQD{T<:Real} <: AbstractKKTSolver{T} # Factorization F::LDLF.LDLFactorization{T} - # TODO: constructor with initial memory allocation - function LDLFactSQD(A::AbstractMatrix{T}) where{T<:Real} + function LDLFactSQD(A::SparseMatrixCSC{T,Int}) where{T<:Real} m, n = size(A) θ = ones(T, n) + dQ = zeros(T, n) S = [ spdiagm(0 => -θ) A'; spzeros(T, m, n) spdiagm(0 => ones(T, m)) @@ -49,9 +51,26 @@ mutable struct LDLFactSQD{T<:Real} <: AbstractKKTSolver{T} # TODO: symbolic factorization only F = LDLF.ldl_analyze(Symmetric(S)) - return new{T}(m, n, A, θ, ones(T, n), ones(T, m), S, F) + return new{T}(m, n, spzeros(T, n, n), dQ, A, θ, ones(T, n), ones(T, m), S, F) end + function LDLFactSQD(A::SparseMatrixCSC{T,Int}, Q::SparseMatrixCSC{T,Int}) where{T<:Real} + m, n = size(A) + θ = ones(T, n) + + dQ = Vector(diag(Q, 0)) + U = sparse(UpperTriangular(Q)) + S = [ + -U - spdiagm(0 => -θ) A'; + spzeros(T, m, n) spdiagm(0 => ones(T, m)) + ] + + F = LDLF.ldl_analyze(Symmetric(S)) + + return new{T}(m, n, Q, dQ, A, θ, ones(T, n), ones(T, m), S, F) + end + LDLFactSQD(A) = LDLFactSQD(sparse(A)) + LDLFactSQD(A, ::ZeroMatrix) = LDLFactSQD(A) end setup(::Type{LDLFactSQD}, A) = LDLFactSQD(A) @@ -95,14 +114,13 @@ function update!( # S is stored as upper-triangular and only its diagonal changes. @inbounds for j in 1:kkt.n k = kkt.S.colptr[1+j] - 1 - kkt.S.nzval[k] = -kkt.θ[j] - regP[j] + kkt.S.nzval[k] = -kkt.dQ[j] -kkt.θ[j] - regP[j] end @inbounds for i in 1:kkt.m k = kkt.S.colptr[1+kkt.n+i] - 1 kkt.S.nzval[k] = regD[i] end - # TODO: PSD-ness checks try LDLF.ldl_factorize!(Symmetric(kkt.S), kkt.F) catch err @@ -124,7 +142,7 @@ function solve!( ξp::Vector{T}, ξd::Vector{T} ) where{T<:Real} m, n = kkt.m, kkt.n - + # Set-up right-hand side ξ = [ξd; ξp]