From 0c341164a1cd387acca934440028100610926573 Mon Sep 17 00:00:00 2001 From: Paul Goulart Date: Wed, 8 Mar 2023 20:45:09 +0000 Subject: [PATCH] release v0.4.1 (#119) * Constraint reduction for inequalities with infinite RHS (#115) * internal changes for sdp efficiency and rust compatibility * Constructor float type defaults (fixes #117) --- CHANGELOG.md | 21 +++- Project.toml | 2 +- README.md | 2 +- src/Clarabel.jl | 12 ++ src/cones/compositecone_type.jl | 37 +++--- src/cones/cone_types.jl | 10 +- src/cones/coneops_compositecone.jl | 4 +- src/cones/coneops_nncone.jl | 8 +- src/cones/coneops_psdtrianglecone.jl | 51 +++++--- src/cones/coneops_zerocone.jl | 2 +- src/info_print.jl | 20 +-- src/kktsolvers/direct-ldl/directldl_utils.jl | 20 +-- src/kktsolvers/kktsolver_directldl.jl | 13 +- src/presolver.jl | 123 +++++++++++++++++++ src/settings.jl | 11 +- src/solution.jl | 36 +++--- src/solver.jl | 11 +- src/types.jl | 33 +++-- src/utils/mathutils.jl | 12 ++ src/variables.jl | 35 ------ test/OptTests/presolve.jl | 90 ++++++++++++++ test/UnitTests/test_constructors.jl | 60 +++++++++ test/run_solver_tests.jl | 1 + 23 files changed, 470 insertions(+), 144 deletions(-) create mode 100644 src/presolver.jl create mode 100644 test/OptTests/presolve.jl create mode 100644 test/UnitTests/test_constructors.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index 9ff9d4bd..f7d486f1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,17 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/). Version numbering in this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). We aim to keep the core solver functionality and minor releases in sync between the Rust/Python and Julia implementations. Small fixes that affect one implementation only may result in the patch release versions differing. +## [0.4.1] - 2023-08-03 +### Changed + +Added optional feature to remove inequality constraints with very large upper bounds. This feature is enabled by default but can be turned off using the `presolve_enable` setting. + +Bug fix in equilibration for NN and zero cones. +### Julia specific changes + +Internal implementation of composite cone logic updated to more closely match the rust version. + +Internal modifications to SDP cone implementation to reduce allocations. ## [0.4.0] - 2023-25-02 ### Changed @@ -43,8 +54,8 @@ Version numbering in this project adheres to [Semantic Versioning](https://semve - Initial release - -[0.4.0]: https://github.com/pyo3/pyo3/compare/v0.4.0...v0.3.0 -[0.3.0]: https://github.com/pyo3/pyo3/compare/v0.3.0...v0.2.0 -[0.2.0]: https://github.com/pyo3/pyo3/compare/v0.2.0...v0.1.0 -[0.1.0]: https://github.com/PyO3/pyo3/tree/0.1.0 +[0.4.1]: https://github.com/oxfordcontrol/Clarabel.jl/compare/v0.4.0...v0.4.1 +[0.4.0]: https://github.com/oxfordcontrol/Clarabel.jl/compare/v0.4.0...v0.3.0 +[0.3.0]: https://github.com/oxfordcontrol/Clarabel.jl/compare/v0.3.0...v0.2.0 +[0.2.0]: https://github.com/oxfordcontrol/Clarabel.jl/compare/v0.2.0...v0.1.0 +[0.1.0]: https://github.com/oxfordcontrol/Clarabel.jl/tree/0.1.0 diff --git a/Project.toml b/Project.toml index 19616680..d06bebf6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Clarabel" uuid = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" authors = ["Paul Goulart "] -version = "0.4.0" +version = "0.4.1" [deps] AMD = "14f7f29c-3bd6-536c-9a0b-7339e30b5a3e" diff --git a/README.md b/README.md index bd2182e8..e35ea25f 100644 --- a/README.md +++ b/README.md @@ -10,7 +10,7 @@ Interior Point Conic Optimization for Julia - +

diff --git a/src/Clarabel.jl b/src/Clarabel.jl index e78b9407..ce619d0b 100644 --- a/src/Clarabel.jl +++ b/src/Clarabel.jl @@ -6,6 +6,17 @@ module Clarabel const DefaultInt = LinearAlgebra.BlasInt const IdentityMatrix = UniformScaling{Bool} + #internal constraint RHS limits. This let block + #hides the INFINITY field in the module and makes + #it accessible only via the get/set provided + let + _INFINITY_DEFAULT = 1e20 + INFINITY = _INFINITY_DEFAULT + global default_infinity() = INFINITY = _INFINITY_DEFAULT; + global set_infinity(v::Float64) = INFINITY = Float64(v) + global get_infinity() = INFINITY + end + #version / release info include("./version.jl") @@ -20,6 +31,7 @@ module Clarabel include("./settings.jl") include("./conicvector.jl") include("./statuscodes.jl") + include("./presolver.jl") include("./types.jl") include("./variables.jl") include("./residuals.jl") diff --git a/src/cones/compositecone_type.jl b/src/cones/compositecone_type.jl index 18d155cb..796ff088 100644 --- a/src/cones/compositecone_type.jl +++ b/src/cones/compositecone_type.jl @@ -8,10 +8,8 @@ struct CompositeCone{T} <: AbstractCone{T} cones::Vector{AbstractCone{T}} - #API type specs and count of each cone type - cone_specs::Vector{SupportedCone} - types::Vector{DataType} - type_counts::Dict{DataType,Int} + #count of each cone type + type_counts::Dict{Type,Int} #overall size of the composite cone numel::DefaultInt @@ -24,38 +22,31 @@ struct CompositeCone{T} <: AbstractCone{T} # the flag for symmetric cone check _is_symmetric::Bool - function CompositeCone{T}(cone_specs::Vector{<:SupportedCone}) where {T} - - #make copy to protect from user interference after setup, - #and explicitly cast as the abstract type. This prevents - #errors arising from input vectors that are all the same cone, - #and therefore more concretely typed than we want - cone_specs = Vector{SupportedCone}(cone_specs) + function CompositeCone{T}(cone_specs::Vector{SupportedCone}) where {T} ncones = length(cone_specs) cones = Vector{AbstractCone{T}}(undef,ncones) - types = Vector{DataType}(undef,ncones) #count the number of each cone type - type_counts = Dict{DataType,Int}() - for coneT in keys(ConeDict) - type_counts[coneT] = count(C->isa(C,coneT), cone_specs) + type_counts = Dict{Type,Int}() + for (key, val) in ConeDict + type_counts[val] = count(C->isa(C,key), cone_specs) end #assumed symmetric to start _is_symmetric = true #create cones with the given dims - for i in eachindex(cone_specs) - types[i] = typeof(cone_specs[i]) - if types[i] == ExponentialConeT - cones[i] = ConeDict[typeof(cone_specs[i])]{T}() + for (i, coneT) in enumerate(cone_specs) + typeT = typeof(coneT) + if typeT == ExponentialConeT + cones[i] = ConeDict[typeT]{T}() _is_symmetric = false - elseif types[i] == PowerConeT - cones[i] = ConeDict[typeof(cone_specs[i])]{T}(T(cone_specs[i].α)) + elseif typeT == PowerConeT + cones[i] = ConeDict[typeT]{T}(T(cone_specs[i].α)) _is_symmetric = false else - cones[i] = ConeDict[typeof(cone_specs[i])]{T}(cone_specs[i].dim) + cones[i] = ConeDict[typeT]{T}(cone_specs[i].dim) end end @@ -68,7 +59,7 @@ struct CompositeCone{T} <: AbstractCone{T} headidx = Vector{Int}(undef,length(cones)) _make_headidx!(headidx,cones) - return new(cones,cone_specs,types,type_counts,numel,degree,headidx,_is_symmetric) + return new(cones,type_counts,numel,degree,headidx,_is_symmetric) end end diff --git a/src/cones/cone_types.jl b/src/cones/cone_types.jl index 0fe3790b..c45eb4ef 100644 --- a/src/cones/cone_types.jl +++ b/src/cones/cone_types.jl @@ -47,7 +47,7 @@ struct NonnegativeCone{T} <: AbstractCone{T} function NonnegativeCone{T}(dim) where {T} - dim >= 1 || throw(DomainError(dim, "dimension must be positive")) + dim >= 0 || throw(DomainError(dim, "dimension must be nonnegative")) w = zeros(T,dim) λ = zeros(T,dim) return new(dim,w,λ) @@ -112,9 +112,10 @@ mutable struct PSDConeWork{T} B::Matrix{T} Hs::Matrix{T} - #workspace for various internal use + #workspace for various internal uses workmat1::Matrix{T} workmat2::Matrix{T} + workmat3::Matrix{T} workvec::Vector{T} function PSDConeWork{T}(n::Int) where {T} @@ -133,10 +134,11 @@ mutable struct PSDConeWork{T} workmat1 = zeros(T,n,n) workmat2 = zeros(T,n,n) + workmat3 = zeros(T,n,n) workvec = zeros(T,(n*(n+1))>>1) return new(cholS,cholZ,SVD,λ,Λisqrt,R,Rinv, - kronRR,B,Hs,workmat1,workmat2,workvec) + kronRR,B,Hs,workmat1,workmat2,workmat3,workvec) end end @@ -239,7 +241,7 @@ PowerCone(args...) = PowerCone{DefaultFloat}(args...) A Dict that maps the user-facing SupportedCone types to the types used internally in the solver. See [SupportedCone](@ref) """ -const ConeDict = Dict( +const ConeDict = Dict{DataType,Type}( ZeroConeT => ZeroCone, NonnegativeConeT => NonnegativeCone, SecondOrderConeT => SecondOrderCone, diff --git a/src/cones/coneops_compositecone.jl b/src/cones/coneops_compositecone.jl index 9c1cb0f8..92142af6 100644 --- a/src/cones/coneops_compositecone.jl +++ b/src/cones/coneops_compositecone.jl @@ -143,7 +143,7 @@ function update_scaling!( # update cone scalings by passing subview to each of # the appropriate cone types. - for (cone,type,si,zi) in zip(cones,cones.types,s.views,z.views) + for (cone,si,zi) in zip(cones,s.views,z.views) @conedispatch is_scaling_success = update_scaling!(cone,si,zi,μ,scaling_strategy) # YC: currently, only check whether SOC variables are in the interior; # we could extend the interior checkfor other cones @@ -247,7 +247,7 @@ function step_length( s = s.views # Force symmetric cones first. - for (cone,type,dzi,dsi,zi,si) in zip(cones,cones.types,dz,ds,z,s) + for (cone,dzi,dsi,zi,si) in zip(cones,dz,ds,z,s) if !is_symmetric(cone) continue end @conedispatch (nextαz,nextαs) = step_length(cone,dzi,dsi,zi,si,settings,α) diff --git a/src/cones/coneops_nncone.jl b/src/cones/coneops_nncone.jl index 32f26e82..7ea23a91 100644 --- a/src/cones/coneops_nncone.jl +++ b/src/cones/coneops_nncone.jl @@ -9,7 +9,7 @@ function rectify_equilibration!( ) where{T} #allow elementwise equilibration scaling - δ .= e + δ .= one(T) return false end @@ -19,7 +19,11 @@ function margins( pd::PrimalOrDualCone ) where{T} - α = minimum(z) + if(length(z) > 0) + α = minimum(z) + else #catch cone of dimension zero + α = floatmax(T) + end # total positive margin is the sum of the # nonnegative elements in the vector, since diff --git a/src/cones/coneops_psdtrianglecone.jl b/src/cones/coneops_psdtrianglecone.jl index 0c83498a..ba99689b 100644 --- a/src/cones/coneops_psdtrianglecone.jl +++ b/src/cones/coneops_psdtrianglecone.jl @@ -14,8 +14,9 @@ function margins( Z = K.work.workmat1 _svec_to_mat!(Z,z,K) - α = eigvals(Symmetric(Z),1:1)[1] #minimum eigenvalue - β = sum(eigvals(Symmetric(Z),0,floatmax(T))) #sum of positive eigenvalues + e = eigvals!(Symmetric(Z)) + α = minimum(e) #minimum eigenvalue + β = reduce((x,y) -> y > 0 ? x + y : x, e, init = 0.) # = sum(e[e.>0]) (no alloc) (α,β) end @@ -94,16 +95,21 @@ function update_scaling!( #R is the same size as L2'*L1, #so use it as temporary workspace - f.R .= L2'*L1 + mul!(f.R,L2',L1) #R = L2'*L1 + f.SVD = svd(f.R) #assemble λ (diagonal), R and Rinv. f.λ .= f.SVD.S f.Λisqrt.diag .= inv.(sqrt.(f.λ)) - # PJG : allocating - f.R .= L1*(f.SVD.V)*f.Λisqrt - f.Rinv .= f.Λisqrt*(f.SVD.U)'*L2' + #f.R = L1*(f.SVD.V)*f.Λisqrt + mul!(f.R,L1,f.SVD.V); + mul!(f.R,f.R,f.Λisqrt) #mul! can take Rinv twice because Λ is diagonal + + #f.Rinv .= f.Λisqrt*(f.SVD.U)'*L2' + mul!(f.Rinv,f.SVD.U',L2') + mul!(f.Rinv,f.Λisqrt,f.Rinv) #mul! can take Rinv twice because Λ is diagonal return is_scaling_success = true end @@ -258,14 +264,18 @@ function mul_W!( (X,Y) = (K.work.workmat1,K.work.workmat2) map((M,v)->_svec_to_mat!(M,v,K),(X,Y),(x,y)) + tmp = K.work.workmat3 R = K.work.R - # PJG : allocating if is_transpose === :T - Y .+= α*(R*X*R') #W^T*x + #Y .+= α*(R*X*R') #W^T*x + mul!(tmp,X,R') + mul!(Y,R,tmp,α,one(T)) else # :N - Y .+= α*(R'*X*R) #W*x + #Y .+= α*(R'*X*R) #W*x + mul!(tmp,R',X) + mul!(Y,tmp,R,α,one(T)) end _mat_to_svec!(y,Y,K) @@ -287,14 +297,18 @@ function mul_Winv!( (X,Y) = (K.work.workmat1,K.work.workmat2) map((M,v)->_svec_to_mat!(M,v,K),(X,Y),(x,y)) + tmp = K.work.workmat3 Rinv = K.work.Rinv - # PJG : allocating if is_transpose === :T - Y .+= α*(Rinv*X*Rinv') #W^{-T}*x + #Y .+= α*(Rinv*X*Rinv') #W^{-T}*x + mul!(tmp,X,Rinv') + mul!(Y,Rinv,tmp,α,one(T)) else # :N - Y .+= α*(Rinv'*X*Rinv) #W^{-1}*x + #Y .+= α*(Rinv'*X*Rinv) #W^{-1}*x + mul!(tmp,Rinv',X) + mul!(Y,tmp,Rinv,α,one(T)) end _mat_to_svec!(y,Y,K) @@ -338,8 +352,13 @@ function circ_op!( (Y,Z) = (K.work.workmat1,K.work.workmat2) map((M,v)->_svec_to_mat!(M,v,K),(Y,Z),(y,z)) - # PJG : allocating - Y .= (Y*Z + Z*Y)/2 + tmp = K.work.workmat3; + + #Y .= (Y*Z + Z*Y)/2 + # NB: works b/c Y and Z are both symmetric + mul!(tmp,Y,Z) + Y .= tmp + symmetric_part!(Y) _mat_to_svec!(x,Y,K) return nothing @@ -383,9 +402,7 @@ function _step_length_psd_component( # we only need to populate the upper # triangle lrscale!(Λisqrt.diag,Δ,Λisqrt.diag) - M = Symmetric(Δ) - - γ = eigvals(M,1:1)[1] #minimum eigenvalue + γ = eigvals!(Symmetric(Δ),1:1)[1] #minimum eigenvalue if γ < 0 return min(inv(-γ),αmax) else diff --git a/src/cones/coneops_zerocone.jl b/src/cones/coneops_zerocone.jl index 9aa402f7..fbc2eded 100644 --- a/src/cones/coneops_zerocone.jl +++ b/src/cones/coneops_zerocone.jl @@ -19,7 +19,7 @@ function rectify_equilibration!( ) where{T} #allow elementwise equilibration scaling - δ .= e + δ .= one(T) return false end diff --git a/src/info_print.jl b/src/info_print.jl index e3dc2fb2..005fa4eb 100644 --- a/src/info_print.jl +++ b/src/info_print.jl @@ -17,18 +17,22 @@ function info_print_configuration( if(settings.verbose == false) return end + if(is_reduced(data.presolver)) + @printf("\npresolve: removed %i constraints\n", count_reduced(data.presolver)) + end + @printf("\nproblem:\n") @printf(" variables = %i\n", data.n) @printf(" constraints = %i\n", data.m) @printf(" nnz(P) = %i\n", nnz(data.P)) @printf(" nnz(A) = %i\n", nnz(data.A)) @printf(" cones (total) = %i\n", length(cones)) - print_conedims_by_type(cones, ZeroConeT) - print_conedims_by_type(cones, NonnegativeConeT) - print_conedims_by_type(cones, SecondOrderConeT) - print_conedims_by_type(cones, PSDTriangleConeT) - print_conedims_by_type(cones, PowerConeT) - print_conedims_by_type(cones, ExponentialConeT) + print_conedims_by_type(cones, ZeroCone) + print_conedims_by_type(cones, NonnegativeCone) + print_conedims_by_type(cones, SecondOrderCone) + print_conedims_by_type(cones, PSDTriangleCone) + print_conedims_by_type(cones, PowerCone) + print_conedims_by_type(cones, ExponentialCone) print_settings(settings, T) @printf("\n") @@ -176,8 +180,8 @@ function print_conedims_by_type(cones::CompositeCone{T}, type) where {T} return #don't report if none end - nvars = map(K->Clarabel.numel(K), cones[isa.(cones.cone_specs,type)]) - name = rpad(string(type)[1:end-5],11) #drops "ConeT part" + nvars = map(K->Clarabel.numel(K), cones[isa.(cones,type)]) + name = rpad(string(nameof(type))[1:end-4],11) #drops "Cone" part @printf(" : %s = %i, ", name, count) if count == 1 diff --git a/src/kktsolvers/direct-ldl/directldl_utils.jl b/src/kktsolvers/direct-ldl/directldl_utils.jl index b2a6dbc0..16c63945 100644 --- a/src/kktsolvers/direct-ldl/directldl_utils.jl +++ b/src/kktsolvers/direct-ldl/directldl_utils.jl @@ -33,16 +33,16 @@ struct LDLDataMap Hsblocks = _allocate_kkt_Hsblocks(Int, cones) #now do the SOC expansion pieces - nsoc = cones.type_counts[Clarabel.SecondOrderConeT] + nsoc = cones.type_counts[Clarabel.SecondOrderCone] p = 2*nsoc SOC_D = zeros(Int,p) SOC_u = Vector{Vector{Int}}(undef,nsoc) SOC_v = Vector{Vector{Int}}(undef,nsoc) - count = 1; + count = 1 for (i,cone) in enumerate(cones) - if isa(cones.cone_specs[i],Clarabel.SecondOrderConeT) + if isa(cone,Clarabel.SecondOrderCone) SOC_u[count] = Vector{Int}(undef,numel(cone)) SOC_v[count] = Vector{Int}(undef,numel(cone)) count += 1 @@ -83,7 +83,7 @@ function _assemble_kkt_matrix( ) where{T} (m,n) = (size(A,1), size(P,1)) - n_socs = cones.type_counts[Clarabel.SecondOrderConeT] + n_socs = cones.type_counts[Clarabel.SecondOrderCone] p = 2*n_socs maps = LDLDataMap(P,A,cones) @@ -158,11 +158,11 @@ function _kkt_assemble_colcounts( #count dense columns for each SOC socidx = 1 #which SOC are we working on? - for i in eachindex(cones) - if isa(cones.cone_specs[i],Clarabel.SecondOrderConeT) + for (i,cone) in enumerate(cones) + if isa(cone,Clarabel.SecondOrderCone) #we will add the u and v columns for this cone - nvars = numel(cones[i]) + nvars = numel(cone) headidx = cones.headidx[i] #which column does u go into? @@ -230,10 +230,10 @@ function _kkt_assemble_fill( #fill in dense columns for each SOC socidx = 1 #which SOC are we working on? - for i in eachindex(cones) - if isa(cones.cone_specs[i],Clarabel.SecondOrderConeT) + for (i,cone) in enumerate(cones) + if isa(cone,Clarabel.SecondOrderCone) - nvars = numel(cones[i]) + nvars = numel(cone) headidx = cones.headidx[i] #which column does u go into (if triu)? diff --git a/src/kktsolvers/kktsolver_directldl.jl b/src/kktsolvers/kktsolver_directldl.jl index 4f1124db..c10650e1 100644 --- a/src/kktsolvers/kktsolver_directldl.jl +++ b/src/kktsolvers/kktsolver_directldl.jl @@ -47,7 +47,7 @@ mutable struct DirectLDLKKTSolver{T} <: AbstractKKTSolver{T} #solving in sparse format. Need this many #extra variables for SOCs - p = 2*cones.type_counts[SecondOrderConeT] + p = 2*cones.type_counts[SecondOrderCone] #LHS/RHS/work for iterative refinement x = Vector{T}(undef,n+m+p) @@ -212,14 +212,13 @@ function _kktsolver_update_inner!( #update the scaled u and v columns. cidx = 1 #which of the SOCs are we working on? - for (i,K) = enumerate(cones) - if isa(cones.cone_specs[i],SecondOrderConeT) - - η2 = K.η^2 + for (i,cone) = enumerate(cones) + if isa(cone,SecondOrderCone) + η2 = cone.η^2 #off diagonal columns (or rows) - _update_values!(ldlsolver,KKT,map.SOC_u[cidx],K.u) - _update_values!(ldlsolver,KKT,map.SOC_v[cidx],K.v) + _update_values!(ldlsolver,KKT,map.SOC_u[cidx],cone.u) + _update_values!(ldlsolver,KKT,map.SOC_v[cidx],cone.v) _scale_values!(ldlsolver,KKT,map.SOC_u[cidx],-η2) _scale_values!(ldlsolver,KKT,map.SOC_v[cidx],-η2) diff --git a/src/presolver.jl b/src/presolver.jl new file mode 100644 index 00000000..85ed1d83 --- /dev/null +++ b/src/presolver.jl @@ -0,0 +1,123 @@ + +struct PresolverRowReductionIndex + + # vector of length = original RHS. Entries are false + # for those rows that should be eliminated before solve + keep_logical::Vector{Bool} + + # vector of length = reduced RHS, taking values + # that map reduced b back to their original index + # This is just findall(keep_logical) and is held for + # efficient solution repopulation + keep_index::Vector{Int64} + +end + +struct Presolver{T} + + # possibly reduced internal copy of user cone specification + cone_specs::Vector{SupportedCone} + + # record of reduced constraints for NN cones with inf bounds + reduce_map::Union{Nothing,PresolverRowReductionIndex} + + # size of original and reduced RHS, respectively + mfull::Int64 + mreduced::Int64 + + # inf bound that was taken from the module level + # and should be applied throughout. Held here so + # that any subsequent change to the module's state + # won't mess up our solver mid-solve + infbound::Float64 + + function Presolver{T}( + A::AbstractMatrix{T}, + b::Vector{T}, + cone_specs::Vector{<:SupportedCone}, + settings::Settings{T} + ) where {T} + + infbound = Clarabel.get_infinity() + + # make copy of cone_specs to protect from user interference after + # setup and explicitly cast as the abstract type. This also prevents + # errors arising from input vectors that are all the same cone, and + # therefore more concretely typed than we want + cone_specs = Vector{SupportedCone}(cone_specs) + mfull = length(b) + + (reduce_map, mreduced) = let + if settings.presolve_enable + reduce_cones!(cone_specs,b,T(infbound)) + else + (nothing,mfull) + end + end + + return new(cone_specs,reduce_map,mfull, mreduced, infbound) + + end + +end + +Presolver(args...) = Presolver{DefaultFloat}(args...) + +is_reduced(ps::Presolver{T}) where {T} = !isnothing(ps.reduce_map) +count_reduced(ps::Presolver{T}) where {T} = ps.mfull - ps.mreduced + + +function reduce_cones!( + cone_specs::Vector{<:SupportedCone}, + b::Vector{T}, + infbound::T) where {T} + + keep_logical = trues(length(b)) + mreduced = length(b) + + # we loop through b and remove any entries that are both infinite + # and in a nonnegative cone + + is_reduced = false + bptr = 1 # index into the b vector + + for (cidx,cone) in enumerate(cone_specs) + + numel_cone = nvars(cone) + + # only try to reduce nn cones. Make a slight contraction + # so that we are firmly "less than" here + infbound *= (1-10*eps(T)) + + if isa(cone, NonnegativeConeT) + num_finite = 0 + for i in bptr:(bptr + numel_cone - 1) + if b[i] < infbound + num_finite += 1 + else + keep_logical[i] = false + mreduced -= 1 + end + end + if num_finite < numel_cone + # contract the cone to a smaller size + cone_specs[cidx] = NonnegativeConeT(num_finite) + is_reduced = true + end + end + + bptr += numel_cone + end + + outoption = + let + if is_reduced + keep_index = findall(keep_logical) + PresolverRowReductionIndex(keep_logical, keep_index) + else + nothing + end + end + + (outoption, mreduced) +end diff --git a/src/settings.jl b/src/settings.jl index 9c9c50e9..ea3d16a1 100644 --- a/src/settings.jl +++ b/src/settings.jl @@ -58,6 +58,8 @@ iterative\\_refinement\\_reltol | 1e-12 | iterative refinement relat iterative\\_refinement\\_abstol | 1e-12 | iterative refinement absolute tolerance iterative\\_refinement\\_max\\_iter | 10 | iterative refinement maximum iterations iterative\\_refinement\\_stop\\_ratio | 5.0 | iterative refinement stalling tolerance +__Preprocessing Settings +presolve_enable | true | enable presolve constraint reduction """ Base.@kwdef mutable struct Settings{T <: AbstractFloat} @@ -115,15 +117,18 @@ Base.@kwdef mutable struct Settings{T <: AbstractFloat} iterative_refinement_abstol::T = 1e-12 iterative_refinement_max_iter::Int = 10 - iterative_refinement_stop_ratio::T = 5 + iterative_refinement_stop_ratio::T = 5 + + #preprocessing + presolve_enable::Bool = true end Settings(args...) = Settings{DefaultFloat}(args...) -function Settings(d::Dict) +function Settings{T}(d::Dict) where{T} - settings = Settings() + settings = Settings{T}() settings_populate!(settings,d) return settings end diff --git a/src/solution.jl b/src/solution.jl index c3477f84..c56c0803 100644 --- a/src/solution.jl +++ b/src/solution.jl @@ -10,11 +10,6 @@ function solution_finalize!( solution.status = info.status solution.obj_val = info.cost_primal - #copy internal variables and undo homogenization - solution.x .= variables.x - solution.z .= variables.z - solution.s .= variables.s - #if we have an infeasible problem, normalize #using κ to get an infeasibility certificate. #Otherwise use τ to get a solution. @@ -25,18 +20,29 @@ function solution_finalize!( scaleinv = one(T) / variables.τ end - @. solution.x *= scaleinv - @. solution.z *= scaleinv - @. solution.s *= scaleinv + #also undo the equilibration + d = data.equilibration.d; dinv = data.equilibration.dinv + e = data.equilibration.e; einv = data.equilibration.einv + cscale = data.equilibration.c[] + + @. solution.x = variables.x * d * scaleinv + + map = data.presolver.reduce_map + if !isnothing(map) + map = data.presolver.reduce_map + @. solution.z[map.keep_index] = variables.z * e * (scaleinv / cscale) + @. solution.s[map.keep_index] = variables.s * einv * scaleinv - #undo the equilibration - d = data.equilibration.d; dinv = data.equilibration.dinv - e = data.equilibration.e; einv = data.equilibration.einv - cscale = data.equilibration.c[] + #eliminated constraints get huge slacks + #and are assumed to be nonbinding + @. solution.s[!map.keep_logical] = T(data.presolver.infbound) + @. solution.z[!map.keep_logical] = zero(T) - @. solution.x *= d - @. solution.z *= e ./ cscale - @. solution.s *= einv + else + @. solution.z = variables.z * e * (scaleinv / cscale) + @. solution.s = variables.s * einv * scaleinv + end + solution.iterations = info.iterations solution.solve_time = info.solve_time diff --git a/src/solver.jl b/src/solver.jl index 8ea1bf28..2cfc81d7 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -90,8 +90,13 @@ function setup!( @timeit s.timers "setup!" begin - s.cones = CompositeCone{T}(cones) - s.data = DefaultProblemData{T}(P,q,A,b,s.cones) + #reduce the cone sizes. (A,b) will be reduced + #within the problem data constructor. Also makes + #an internal copy of the user cone specification + presolver = Presolver{T}(A,b,cones,s.settings) + + s.cones = CompositeCone{T}(presolver.cone_specs) + s.data = DefaultProblemData{T}(P,q,A,b,s.cones,presolver) s.data.m == s.cones.numel || throw(DimensionMismatch()) s.variables = DefaultVariables{T}(s.data.n,s.cones) @@ -116,7 +121,7 @@ function setup!( s.prev_vars = DefaultVariables{T}(s.data.n,s.cones) # user facing results go here - s.solution = DefaultSolution{T}(s.data.m,s.data.n) + s.solution = DefaultSolution{T}(s.data.presolver.mfull,s.data.n) end diff --git a/src/types.jl b/src/types.jl index e5fd3f4c..3c7580e8 100644 --- a/src/types.jl +++ b/src/types.jl @@ -160,31 +160,50 @@ mutable struct DefaultProblemData{T} <: AbstractProblemData{T} normq::T normb::T + presolver::Presolver{T} + function DefaultProblemData{T}( P::AbstractMatrix{T}, q::AbstractVector{T}, A::AbstractMatrix{T}, b::AbstractVector{T}, - cones::CompositeCone{T} + cones::CompositeCone{T}, + presolver::Presolver{T}, ) where {T} # dimension checks will have already been # performed during problem setup, so skip here - (m,n) = size(A) #take an internal copy of all problem #data, since we are going to scale it P = triu(P) - A = deepcopy(A) q = deepcopy(q) - b = deepcopy(b) + + (A,b) = let + map = presolver.reduce_map; + if !isnothing(map) + ( + A[map.keep_logical,:], + b[map.keep_logical] + ) + else + (deepcopy(A),deepcopy(b)) + end + end + + #cap entries in b at INFINITY. This is important + #for inf values that were not in a reduced cone + @. b .= min(b,T(presolver.infbound)) + + #this ensures m is the *reduced* size m + (m,n) = size(A) equilibration = DefaultEquilibration{T}(n,cones) normq = norm(q, Inf) normb = norm(b, Inf) - new(P,q,A,b,n,m,equilibration,normq,normb) + new(P,q,A,b,n,m,equilibration,normq,normb,presolver) end @@ -343,8 +362,8 @@ function Solver{T}() where {T} end #partial user defined settings -function Solver(d::Dict) where {T} - Solver{T}(Settings(d)) +function Solver{T}(d::Dict) where {T} + Solver{T}(Settings{T}(d)) end Solver(args...; kwargs...) = Solver{DefaultFloat}(args...; kwargs...) diff --git a/src/utils/mathutils.jl b/src/utils/mathutils.jl index 32806760..2d6dc198 100644 --- a/src/utils/mathutils.jl +++ b/src/utils/mathutils.jl @@ -219,6 +219,18 @@ function lrscale!(L::AbstractVector{T},M::Matrix{T},R::AbstractVector{T}) where return M end +# Set A = (A + A') / 2. Assumes A is real +function symmetric_part!(A::Matrix{T}) where T <: Real + n = LinearAlgebra.checksquare(A) + @inbounds for r in 1:n + for c in 1:r-1 + val = (A[r, c] + A[c, r]) / 2 + A[r,c] = A[c,r] = val + end + end + return A +end + #Julia SparseArrays dot function is very slow for Symmtric #matrices. See https://github.com/JuliaSparse/SparseArrays.jl/issues/83 diff --git a/src/variables.jl b/src/variables.jl index 8a0bcfdf..e44c894f 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -225,41 +225,6 @@ function variables_unit_initialization!( return nothing end -function variables_finalize!( - variables::DefaultVariables{T}, - equil::DefaultEquilibration{T}, - status::SolverStatus -) where {T} - - #undo the homogenization - # - #if we have an infeasible problem, normalize - #using κ to get an infeasibility certificate. - #Otherwise use τ to get a solution. - if(status == PRIMAL_INFEASIBLE || status == DUAL_INFEASIBLE) - scaleinv = one(T) / variables.κ - else - scaleinv = one(T) / variables.τ - end - - @. variables.x *= scaleinv - @. variables.z *= scaleinv - @. variables.s *= scaleinv - variables.τ *= scaleinv - variables.κ *= scaleinv - - #undo the equilibration - d = equil.d; dinv = equil.dinv - e = equil.e; einv = equil.einv - cscale = equil.c[] - - @. variables.x *= d - @. variables.z *= e ./ cscale - @. variables.s *= einv - -end - - function variables_rescale!(variables) scale = max(variables.τ,variables.κ) diff --git a/test/OptTests/presolve.jl b/test/OptTests/presolve.jl new file mode 100644 index 00000000..22cc7317 --- /dev/null +++ b/test/OptTests/presolve.jl @@ -0,0 +1,90 @@ +using Test, LinearAlgebra, SparseArrays + +#if not run in full test setup, just do it for one float type +@isdefined(UnitTestFloats) || (UnitTestFloats = [Float64]) + +function presolver_test_data(Type::Type{T}) where {T <: AbstractFloat} + + P = sparse(I(3)*one(T)) + A = sparse(I(3)*one(T)) + A = [A;-A].*2 + c = T[3.;-2.;1.] + b = ones(T,6) + cones = [Clarabel.NonnegativeConeT(3), Clarabel.NonnegativeConeT(3)] + + return (P,c,A,b,cones) +end + + +@testset "Presolver Tests" begin + + for FloatT in UnitTestFloats + + tol = FloatT(1e-3) + + #force default infinite bounds in solver + Clarabel.default_infinity() + + @testset "Presolver tests (T = $(FloatT))" begin + + @testset "single unbounded constraint" begin + + # Just one bound to zero + P,c,A,b,cones = presolver_test_data(FloatT) + solver = Clarabel.Solver{FloatT}() + settings = Clarabel.Settings{FloatT}() + b[4] = 1e30 + Clarabel.setup!(solver,P,c,A,b,cones) + Clarabel.solve!(solver) + @test solver.solution.status == Clarabel.SOLVED + @test (length(solver.variables.z) == 5) + @test (all(solver.solution.z[4] .== zero(FloatT))) + @test (all(solver.solution.s[4] .== Clarabel.get_infinity())) + + end + + @testset "completely redundant cone" begin + + # all bounds on first cone huge + P,c,A,b,cones = presolver_test_data(FloatT) + solver = Clarabel.Solver{FloatT}() + settings = Clarabel.Settings{FloatT}() + b[1:3] .= 1e30 + Clarabel.setup!(solver,P,c,A,b,cones) + Clarabel.solve!(solver) + @test solver.solution.status == Clarabel.SOLVED + @test (length(solver.variables.z) == 3) + @test (all(solver.solution.z[1:3] .== zero(FloatT))) + @test (all(solver.solution.s[1:3] .== Clarabel.get_infinity())) + @test isapprox(norm(solver.solution.x .- FloatT[-0.5; 2; -0.5]), zero(FloatT), atol=tol) + + end + + @testset "every constraint redundant" begin + + # all bounds on first cone huge + P,c,A,b,cones = presolver_test_data(FloatT) + solver = Clarabel.Solver{FloatT}() + settings = Clarabel.Settings{FloatT}() + b[1:6] .= 1e30 + Clarabel.setup!(solver,P,c,A,b,cones) + Clarabel.solve!(solver) + @test solver.solution.status == Clarabel.SOLVED + @test (length(solver.variables.z) == 0) + @test isapprox(norm(solver.solution.x .+ c), zero(FloatT), atol=tol) + + end + + @testset "settable bound" begin + + thebound = Clarabel.get_infinity() + Clarabel.set_infinity(1e21) + @test (Clarabel.get_infinity() == 1e21) + Clarabel.set_infinity(thebound) + @test (Clarabel.get_infinity() == thebound) + + end + + end #end Presolver Tests (FloatT)" + end # UnitTestFloats +end \ No newline at end of file diff --git a/test/UnitTests/test_constructors.jl b/test/UnitTests/test_constructors.jl new file mode 100644 index 00000000..8fd90017 --- /dev/null +++ b/test/UnitTests/test_constructors.jl @@ -0,0 +1,60 @@ +using Test, Clarabel + + +@testset "Settings and Solver constructors" begin + + for FloatT in [Float32,Float64,BigFloat] + + @testset "Settings and Solver (T = $(FloatT))" begin + + @test begin + #no settings + Clarabel.Solver{FloatT}() + + #settings via struct + Clarabel.Solver{FloatT}(Clarabel.Settings{FloatT}()) + + #settings via dict + Clarabel.Solver{FloatT}(Dict("verbose" => true)) + + #settings via parameter list + Clarabel.Settings{FloatT}(verbose = true) + + #settings via dict + Clarabel.Settings{FloatT}(Dict("verbose" => true)) + + true + end + + end + end + + @testset "Settings and Solver constructors (default types)" begin + + @test begin + #no settings + Clarabel.Solver() + + #settings via struct + Clarabel.Solver(Clarabel.Settings()) + + #settings via dict + Clarabel.Solver(Dict("verbose" => true)) + + #settings via parameter list + Clarabel.Settings(verbose = true) + + #settings via dict + Clarabel.Settings(Dict("verbose" => true)) + + #fail on mixed types (assumes 64 bit default) + @test_throws MethodError Clarabel.Solver{Float32}(Clarabel.Settings{Float64}()) + @test_throws MethodError Clarabel.Solver(Clarabel.Settings{Float32}()) + @test_throws MethodError Clarabel.Solver{Float32}(Clarabel.Settings()) + + true + end + end + +end +nothing diff --git a/test/run_solver_tests.jl b/test/run_solver_tests.jl index ad867aad..badfffd6 100644 --- a/test/run_solver_tests.jl +++ b/test/run_solver_tests.jl @@ -16,6 +16,7 @@ UnitTestFloats = [Float64,BigFloat] include("./OptTests/basic_exp.jl") include("./OptTests/basic_pow.jl") include("./OptTests/basic_sdp.jl") + include("./OptTests/presolve.jl") end