diff --git a/.gitignore b/.gitignore index 57bff62ff..3a3307812 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ *.jl.mem deps/deps.jl scratch.jl +*/scratch.jl diff --git a/Manifest.toml b/Manifest.toml index 97e16297d..a4ca2cbca 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -15,9 +15,9 @@ version = "0.8.10" [[BinaryProvider]] deps = ["Libdl", "Pkg", "SHA", "Test"] -git-tree-sha1 = "48c147e63431adbcee69bc40b04c3f0fec0a4982" +git-tree-sha1 = "38be61810f280c3c478f5c38aaf387f8f9199275" uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232" -version = "0.5.0" +version = "0.5.1" [[Combinatorics]] deps = ["LinearAlgebra", "Polynomials", "Test"] @@ -83,6 +83,12 @@ version = "0.9.0" deps = ["LinearAlgebra", "Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +[[IterativeSolvers]] +deps = ["LinearAlgebra", "Pkg", "Printf", "Random", "RecipesBase", "SparseArrays", "Test"] +git-tree-sha1 = "a7b8056ecd43aa0d112d4b5d04180194d12c3c62" +uuid = "42fd0dbc-a981-5370-80f2-aaf504508153" +version = "0.7.1" + [[JSON]] deps = ["Dates", "Distributed", "Mmap", "Pkg", "Sockets", "Test", "Unicode"] git-tree-sha1 = "fec8e4d433072731466d37ed0061b3ba7f70eeb9" @@ -108,7 +114,7 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" [[MathOptInterface]] deps = ["Compat", "Unicode"] -git-tree-sha1 = "ba12e7ce825c1458c03f88aae84fa630d882d303" +git-tree-sha1 = "14bbbdbc746e70bbae6cc96d1fa58f14e02fae4d" repo-rev = "master" repo-url = "https://github.com/JuliaOpt/MathOptInterface.jl.git" uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" @@ -145,6 +151,12 @@ uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" deps = ["Serialization"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +[[RecipesBase]] +deps = ["Pkg", "Random", "Test"] +git-tree-sha1 = "0b3cb370ee4dc00f47f1193101600949f3dcf884" +uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" +version = "0.6.0" + [[Reexport]] deps = ["Pkg"] git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0" diff --git a/Project.toml b/Project.toml index f1a950eec..a14d08412 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" diff --git a/examples/envelope/envelope.jl b/examples/envelope/envelope.jl index 62c31bed8..d899bd2f0 100644 --- a/examples/envelope/envelope.jl +++ b/examples/envelope/envelope.jl @@ -9,15 +9,23 @@ available at https://arxiv.org/abs/1712.01792 =# using Hypatia -using SparseArrays using LinearAlgebra +using SparseArrays using DelimitedFiles using Random -function build_envelope!(alf::Hypatia.HypatiaOpt, npoly::Int, deg::Int, n::Int, d::Int; use_data::Bool=false, dense::Bool=false, rseed::Int=1) - @assert deg <= d +function build_envelope!( + npoly::Int, + deg::Int, + n::Int, + d::Int; + use_data::Bool = false, + dense::Bool = false, + rseed::Int = 1, + ) # generate interpolation + @assert deg <= d (L, U, pts, P0, P, w) = Hypatia.interpolate(n, d, calc_w=true) LWts = fill(binomial(n+d-1, n), n) wtVals = 1.0 .- pts.^2 @@ -27,7 +35,7 @@ function build_envelope!(alf::Hypatia.HypatiaOpt, npoly::Int, deg::Int, n::Int, if dense A = repeat(Array(1.0I, U, U), outer=(1, npoly)) else - A = repeat(sparse(1.0I, U, U), outer=(1, npoly)) # TODO maybe construct without repeat + A = repeat(sparse(1.0I, U, U), outer=(1, npoly)) end G = Diagonal(-1.0I, npoly*U) # TODO uniformscaling b = w @@ -44,17 +52,42 @@ function build_envelope!(alf::Hypatia.HypatiaOpt, npoly::Int, deg::Int, n::Int, cone = Hypatia.Cone([Hypatia.DualSumOfSquaresCone(U, [P, PWts...]) for k in 1:npoly], [1+(k-1)*U:k*U for k in 1:npoly]) - return Hypatia.load_data!(alf, c, A, b, G, h, cone) + return (c, A, b, G, h, cone) end -# alf = Hypatia.HypatiaOpt(maxiter=100, verbose=true) +function run_envelope() + # optionally use fixed data in folder + # select number of polynomials and degrees for the envelope + # select dimension and SOS degree (to be squared) + (c, A, b, G, h, cone) = + # build_envelope!(2, 5, 1, 5, use_data=true) + # build_envelope!(2, 5, 2, 8) + # build_envelope!(3, 5, 3, 5) + build_envelope!(2, 3, 1, 4, dense=false) -# optionally use fixed data in folder -# select number of polynomials and degrees for the envelope -# select dimension and SOS degree (to be squared) -# build_envelope!(alf, 2, 5, 1, 5, use_data=true) -# build_envelope!(alf, 2, 5, 2, 8) -# build_envelope!(alf, 3, 5, 3, 5) -# build_envelope!(alf, 2, 3, 3, 5, dense=false) + Hypatia.check_data(c, A, b, G, h, cone) + (c1, A1, b1, G1, prkeep, dukeep, Q2, RiQ1) = Hypatia.preprocess_data(c, A, b, G, useQR=true) + L = Hypatia.QRSymmCache(c1, A1, b1, G1, h, cone, Q2, RiQ1) -# @time Hypatia.solve!(alf) + opt = Hypatia.Optimizer(maxiter=100, verbose=false) + Hypatia.load_data!(opt, c1, A1, b1, G1, h, cone, L) + Hypatia.solve!(opt) + + x = zeros(length(c)) + x[dukeep] = Hypatia.get_x(opt) + y = zeros(length(b)) + y[prkeep] = Hypatia.get_y(opt) + s = Hypatia.get_s(opt) + z = Hypatia.get_z(opt) + + status = Hypatia.get_status(opt) + solvetime = Hypatia.get_solvetime(opt) + pobj = Hypatia.get_pobj(opt) + dobj = Hypatia.get_dobj(opt) + + # @show status + # @show x + # @show pobj + # @show dobj + return nothing +end diff --git a/examples/lp/lp.jl b/examples/lp/lp.jl index ba53ce232..47c8c6bc9 100644 --- a/examples/lp/lp.jl +++ b/examples/lp/lp.jl @@ -11,7 +11,16 @@ using SparseArrays using DelimitedFiles using Random -function build_lp!(alf::Hypatia.HypatiaOpt, m::Int, n::Int; use_data::Bool=false, dense::Bool=false, nzfrac::Float64=1/sqrt(n), tosparse::Bool=false, rseed::Int=1) +function build_lp!( + m::Int, + n::Int; + use_data::Bool = false, + dense::Bool = false, + nzfrac::Float64 = 1/sqrt(n), + tosparse::Bool = false, + rseed::Int = 1, + ) + # set up problem data if use_data # use provided data in data folder @@ -35,17 +44,43 @@ function build_lp!(alf::Hypatia.HypatiaOpt, m::Int, n::Int; use_data::Bool=false end G = Diagonal(-1.0I, n) # TODO uniformscaling h = zeros(n) + cone = Hypatia.Cone([Hypatia.NonnegativeCone(n)], [1:n]) - return Hypatia.load_data!(alf, c, A, b, G, h, cone) + return (c, A, b, G, h, cone) end -# alf = Hypatia.HypatiaOpt(maxiter=100, verbose=true) +function run_lp() + # optionally use fixed data in folder + # select the random matrix size, dense/sparse, sparsity fraction + (c, A, b, G, h, cone) = + # build_lp!(500, 1000, use_data=true) + # build_lp!(500, 1000) + build_lp!(15, 20) + + Hypatia.check_data(c, A, b, G, h, cone) + (c1, A1, b1, G1, prkeep, dukeep, Q2, RiQ1) = Hypatia.preprocess_data(c, A, b, G, useQR=true) + L = Hypatia.QRSymmCache(c1, A1, b1, G1, h, cone, Q2, RiQ1) -# optionally use fixed data in folder -# select the random matrix size, dense/sparse, sparsity fraction -# build_lp!(alf, 500, 1000) -# build_lp!(alf, 500, 1000, use_data=true) + opt = Hypatia.Optimizer(maxiter=100, verbose=false) + Hypatia.load_data!(opt, c1, A1, b1, G1, h, cone, L) + Hypatia.solve!(opt) -# solve it -# @time Hypatia.solve!(alf) + x = zeros(length(c)) + x[dukeep] = Hypatia.get_x(opt) + y = zeros(length(b)) + y[prkeep] = Hypatia.get_y(opt) + s = Hypatia.get_s(opt) + z = Hypatia.get_z(opt) + + status = Hypatia.get_status(opt) + solvetime = Hypatia.get_solvetime(opt) + pobj = Hypatia.get_pobj(opt) + dobj = Hypatia.get_dobj(opt) + + # @show status + # @show x + # @show pobj + # @show dobj + return nothing +end diff --git a/examples/namedpoly/namedpoly.jl b/examples/namedpoly/namedpoly.jl index 37b69a031..bf368fa32 100644 --- a/examples/namedpoly/namedpoly.jl +++ b/examples/namedpoly/namedpoly.jl @@ -12,7 +12,7 @@ using Hypatia using LinearAlgebra # list of currently available named polynomials, see https://people.sc.fsu.edu/~jburkardt/py_src/polynomials/polynomials.html -const polys = Dict{Symbol,NamedTuple}( +polys = Dict{Symbol,NamedTuple}( :butcher => (n=6, lbs=[-1,-0.1,-0.1,-1,-0.1,-0.1], ubs=[0,0.9,0.5,-0.1,-0.05,-0.03], deg=3, fn=((u,v,w,x,y,z) -> z*v^2+y*w^2-u*x^2+x^3+x^2-(1/3)*u+(4/3)*x) ), @@ -48,7 +48,11 @@ const polys = Dict{Symbol,NamedTuple}( ), ) -function build_namedpoly!(alf::Hypatia.HypatiaOpt, polyname::Symbol, d::Int) +function build_namedpoly!( + polyname::Symbol, + d::Int, + ) + # get data for named polynomial (n, lbs, ubs, deg, fn) = polys[polyname] if d < ceil(Int, deg/2) @@ -63,34 +67,58 @@ function build_namedpoly!(alf::Hypatia.HypatiaOpt, polyname::Symbol, d::Int) pts .+= (ubs + lbs)'/2 wtVals = (pts .- lbs') .* (ubs' .- pts) LWts = fill(binomial(n+d-1, n), n) - PWts = [Diagonal(sqrt.(wtVals[:,j]))*P0[:,1:LWts[j]] for j in 1:n] + PWts = [Diagonal(sqrt.(wtVals[:, j])) * P0[:, 1:LWts[j]] for j in 1:n] # set up problem data A = ones(1, U) b = [1.0,] - c = [fn(pts[j,:]...) for j in 1:U] + c = [fn(pts[j, :]...) for j in 1:U] G = Diagonal(-1.0I, U) # TODO uniformscaling? h = zeros(U) cone = Hypatia.Cone([Hypatia.DualSumOfSquaresCone(U, [P0, PWts...])], [1:U]) - return Hypatia.load_data!(alf, c, A, b, G, h, cone) + return (c, A, b, G, h, cone) end -# alf = Hypatia.HypatiaOpt(maxiter=100, verbose=false) - -# select the named polynomial to minimize and the SOS degree (to be squared) -# build_namedpoly!(alf, :butcher, 2) -# build_namedpoly!(alf, :caprasse, 4) -# build_namedpoly!(alf, :goldsteinprice, 7) -# build_namedpoly!(alf, :heart, 2) -# build_namedpoly!(alf, :lotkavolterra, 3) -# build_namedpoly!(alf, :magnetism7, 2) -# build_namedpoly!(alf, :motzkin, 7) -# build_namedpoly!(alf, :reactiondiffusion, 4) -# build_namedpoly!(alf, :robinson, 8) -# build_namedpoly!(alf, :rosenbrock, 4) -# build_namedpoly!(alf, :schwefel, 3) - -# solve it -# @time Hypatia.solve!(alf) +function run_namedpoly() + # select the named polynomial to minimize and the SOS degree (to be squared) + (c, A, b, G, h, cone) = + # build_namedpoly!(:butcher, 2) + # build_namedpoly!(:caprasse, 4) + # build_namedpoly!(:goldsteinprice, 7) + # build_namedpoly!(:heart, 2) + # build_namedpoly!(:lotkavolterra, 3) + # build_namedpoly!(:magnetism7, 2) + # build_namedpoly!(:motzkin, 7) + build_namedpoly!(:reactiondiffusion, 3) + # build_namedpoly!(:robinson, 8) + # build_namedpoly!(:rosenbrock, 4) + # build_namedpoly!(:schwefel, 3) + + Hypatia.check_data(c, A, b, G, h, cone) + (c1, A1, b1, G1, prkeep, dukeep, Q2, RiQ1) = Hypatia.preprocess_data(c, A, b, G, useQR=true) + L = Hypatia.QRSymmCache(c1, A1, b1, G1, h, cone, Q2, RiQ1) + + opt = Hypatia.Optimizer(maxiter=100, verbose=false) + Hypatia.load_data!(opt, c1, A1, b1, G1, h, cone, L) + Hypatia.solve!(opt) + + x = zeros(length(c)) + x[dukeep] = Hypatia.get_x(opt) + y = zeros(length(b)) + y[prkeep] = Hypatia.get_y(opt) + s = Hypatia.get_s(opt) + z = Hypatia.get_z(opt) + + status = Hypatia.get_status(opt) + solvetime = Hypatia.get_solvetime(opt) + pobj = Hypatia.get_pobj(opt) + dobj = Hypatia.get_dobj(opt) + + # @show status + # @show x + # @show pobj + # @show dobj + return nothing +end diff --git a/src/Hypatia.jl b/src/Hypatia.jl index cb21c4db9..677dfc684 100644 --- a/src/Hypatia.jl +++ b/src/Hypatia.jl @@ -4,14 +4,17 @@ Copyright 2018, Chris Coey and contributors module Hypatia using Printf - using SparseArrays using LinearAlgebra - using ForwardDiff - using DiffResults + using SparseArrays + import FFTW + import Combinatorics include("interpolation.jl") + + import ForwardDiff + import DiffResults include("cone.jl") - for primcone in [ + for primitivecone in [ "orthant", "dualsumofsquares", "secondorder", @@ -21,11 +24,21 @@ module Hypatia "positivesemidefinite", "ellinfinity", ] - include(joinpath(@__DIR__, "primitivecones", primcone * ".jl")) + include(joinpath(@__DIR__, "primitivecones", primitivecone * ".jl")) end + + import IterativeSolvers + include("linearsystem.jl") + for linsyssolver in [ + "qrsymm", + "naive", + ] + include(joinpath(@__DIR__, "linsyssolvers", linsyssolver * ".jl")) + end + include("nativeinterface.jl") import MathOptInterface - MOI = MathOptInterface + const MOI = MathOptInterface include("mathoptinterface.jl") end diff --git a/src/cone.jl b/src/cone.jl index d98d435f6..bed5ba6af 100644 --- a/src/cone.jl +++ b/src/cone.jl @@ -1,3 +1,6 @@ +#= +Copyright 2018, Chris Coey and contributors +=# # cone object abstract type PrimitiveCone end @@ -5,11 +8,11 @@ abstract type PrimitiveCone end # TODO reorder primitive cones so easiest ones to check incone are first mutable struct Cone prms::Vector{PrimitiveCone} - idxs::Vector{AbstractVector{Int}} + idxs::Vector{UnitRange{Int}} end -Cone() = Cone(PrimitiveCone[], AbstractVector{Int}[]) +Cone() = Cone(PrimitiveCone[], UnitRange{Int}[]) -function addprimitivecone!(cone::Cone, prm::PrimitiveCone, idx::AbstractVector{Int}) +function addprimitivecone!(cone::Cone, prm::PrimitiveCone, idx::UnitRange{Int}) @assert dimension(prm) == length(idx) push!(cone.prms, prm) push!(cone.idxs, idx) diff --git a/src/interpolation.jl b/src/interpolation.jl index 6258289b2..d3c5dc23e 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -11,9 +11,6 @@ and Matlab files in the packages - Padua2DM by M. Caliari, S. De Marchi, A. Sommariva, and M. Vianello http://profs.sci.univr.it/~caliari/software.htm =# -import FFTW -import Combinatorics - function interpolate(n, d; calc_w::Bool=false) if n == 1 return cheb2_data(d, calc_w) @@ -170,7 +167,11 @@ function approxfekete_data(n::Int, d::Int, calc_w::Bool) for xp in Combinatorics.multiexponents(n, t) col += 1 for j in 1:n - m[col] *= iseven(xp[j]) ? 2/(1 - xp[j]^2) : 0 + if iseven(xp[j]) + m[col] *= 2/(1 - xp[j]^2) + else + m[col] = 0.0 + end M[:,col] .*= u[j][:,xp[j]+1] end end diff --git a/src/linearsystem.jl b/src/linearsystem.jl new file mode 100644 index 000000000..3796340f0 --- /dev/null +++ b/src/linearsystem.jl @@ -0,0 +1,8 @@ +#= +Copyright 2018, Chris Coey and contributors + +caches for various linear system solvers, for precomputation and memory allocation +# TODO put functions for each cache type into separate files in a new folder +=# + +abstract type LinSysCache end diff --git a/src/linsyssolvers/naive.jl b/src/linsyssolvers/naive.jl new file mode 100644 index 000000000..7d9d751d5 --- /dev/null +++ b/src/linsyssolvers/naive.jl @@ -0,0 +1,148 @@ +#= +naive method that simply performs one high-dimensional linear system solve +TODO currently only does dense operations, needs to work for sparse +=# +mutable struct NaiveCache <: LinSysCache + cone + c + A + b + G + h + LHS3 + LHS3copy + rhs3 + LHS6 + LHS6copy + rhs6 + tyk + tzk + tkk + tsk + ttk + + function NaiveCache( + c::Vector{Float64}, + A::AbstractMatrix{Float64}, + b::Vector{Float64}, + G::AbstractMatrix{Float64}, + h::Vector{Float64}, + cone::Cone, + ) + + (n, p, q) = (length(c), length(b), length(h)) + L = new() + L.cone = cone + L.c = c + L.A = A + L.b = b + L.G = G + L.h = h + L.tyk = n+1 + L.tzk = L.tyk + p + L.tkk = L.tzk + q + L.tsk = L.tkk + 1 + L.ttk = L.tsk + q + # tx ty tz + L.LHS3 = [ + zeros(n,n) A' G'; + A zeros(p,p) zeros(p,q); + G zeros(q,p) Matrix(-1.0I,q,q); + ] + L.LHS3copy = similar(L.LHS3) + L.rhs3 = zeros(L.tkk-1) + # tx ty tz kap ts tau + L.LHS6 = [ + zeros(n,n) A' G' zeros(n) zeros(n,q) c; + -A zeros(p,p) zeros(p,q) zeros(p) zeros(p,q) b; + zeros(q,n) zeros(q,p) Matrix(1.0I,q,q) zeros(q) Matrix(1.0I,q,q) zeros(q); + zeros(1,n) zeros(1,p) zeros(1,q) 1.0 zeros(1,q) 1.0; + -G zeros(q,p) zeros(q,q) zeros(q) Matrix(-1.0I,q,q) h; + -c' -b' -h' -1.0 zeros(1,q) 0.0; + ] + L.LHS6copy = similar(L.LHS6) + L.rhs6 = zeros(L.ttk) + + return L + end +end + +# solve system for x, y, z +function solvelinsys3!( + rhs_tx::Vector{Float64}, + rhs_ty::Vector{Float64}, + rhs_tz::Vector{Float64}, + mu::Float64, + L::NaiveCache; + identityH::Bool = false, + ) + + rhs = L.rhs3 + rhs[1:L.tyk-1] = rhs_tx + @. rhs[L.tyk:L.tzk-1] = -rhs_ty + @. rhs[L.tzk:L.tkk-1] = -rhs_tz + + @. L.LHS3copy = L.LHS3 + if !identityH + for k in eachindex(L.cone.prms) + idxs = L.tzk - 1 .+ L.cone.idxs[k] + dim = dimension(L.cone.prms[k]) + calcHarr_prm!(view(L.LHS3copy, idxs, idxs), Matrix(-mu*I, dim, dim), L.cone.prms[k]) + end + end + + F = bunchkaufman!(Symmetric(L.LHS3copy)) + ldiv!(F, rhs) + + @. @views begin + rhs_tx = rhs[1:L.tyk-1] + rhs_ty = rhs[L.tyk:L.tzk-1] + rhs_tz = rhs[L.tzk:L.tkk-1] + end + + return nothing +end + +# solve system for x, y, z, s, kap, tau +function solvelinsys6!( + rhs_tx::Vector{Float64}, + rhs_ty::Vector{Float64}, + rhs_tz::Vector{Float64}, + rhs_ts::Vector{Float64}, + rhs_kap::Float64, + rhs_tau::Float64, + mu::Float64, + tau::Float64, + L::NaiveCache, + ) + + rhs = L.rhs6 + rhs[1:L.tyk-1] = rhs_tx + rhs[L.tyk:L.tzk-1] = rhs_ty + rhs[L.tzk:L.tkk-1] = rhs_tz + rhs[L.tkk] = rhs_kap + rhs[L.tsk:L.ttk-1] = rhs_ts + rhs[end] = rhs_tau + + @. L.LHS6copy = L.LHS6 + for k in eachindex(L.cone.prms) + dim = dimension(L.cone.prms[k]) + calcHarr_prm!(view(L.LHS6copy, L.tzk - 1 .+ L.cone.idxs[k], L.tsk - 1 .+ L.cone.idxs[k]), Matrix(mu*I, dim, dim), L.cone.prms[k]) + end + + L.LHS6copy[L.tkk, end] = mu/tau/tau + + F = qr!(L.LHS6copy) + ldiv!(F, rhs) + + @. @views begin + rhs_tx = rhs[1:L.tyk-1] + rhs_ty = rhs[L.tyk:L.tzk-1] + rhs_tz = rhs[L.tzk:L.tkk-1] + rhs_ts = rhs[L.tsk:L.ttk-1] + end + dir_kap = rhs[L.tkk] + dir_tau = rhs[end] + + return (dir_kap, dir_tau) +end diff --git a/src/linsyssolvers/qrsymm.jl b/src/linsyssolvers/qrsymm.jl new file mode 100644 index 000000000..b82b0f631 --- /dev/null +++ b/src/linsyssolvers/qrsymm.jl @@ -0,0 +1,234 @@ +#= +CVXOPT method: solve two symmetric linear systems and combine solutions +QR plus either Cholesky factorization or iterative conjugate gradients method +(1) eliminate equality constraints via QR of A' +(2) solve reduced system by Cholesky or iterative method +|0 A' G'| * |ux| = |bx| +|A 0 0 | |uy| |by| +|G 0 M | |uz| |bz| +where M = -I (for initial iterate only) or M = -Hi (Hi is Hessian inverse, pre-scaled by 1/mu) +=# +mutable struct QRSymmCache <: LinSysCache + # TODO can remove some of the intermediary prealloced arrays after github.com/JuliaLang/julia/issues/23919 is resolved + useiterative + cone + c + b + G + h + Q2 + RiQ1 + HG + GHG + GHGQ2 + Q2GHGQ2 + bxGHbz + Q1x + rhs + Q2div + Q2x + GHGxi + HGxi + x1 + y1 + z1 + # for iterative only + cgstate + lprecond + Q2sol + + function QRSymmCache( + c::Vector{Float64}, + A::AbstractMatrix{Float64}, + b::Vector{Float64}, + G::AbstractMatrix{Float64}, + h::Vector{Float64}, + cone::Cone, + Q2::AbstractMatrix{Float64}, + RiQ1::AbstractMatrix{Float64}; + useiterative::Bool = false, + ) + + L = new() + (n, p, q) = (length(c), length(b), length(h)) + nmp = n - p + L.useiterative = useiterative + L.cone = cone + L.c = c + L.b = b + L.G = G + L.h = h + L.Q2 = Q2 + L.RiQ1 = RiQ1 + L.HG = Matrix{Float64}(undef, q, n) # TODO don't enforce dense on some + L.GHG = Matrix{Float64}(undef, n, n) + L.GHGQ2 = Matrix{Float64}(undef, n, nmp) + L.Q2GHGQ2 = Matrix{Float64}(undef, nmp, nmp) + L.bxGHbz = Vector{Float64}(undef, n) + L.Q1x = Vector{Float64}(undef, n) + L.rhs = Vector{Float64}(undef, n) + L.Q2div = Vector{Float64}(undef, nmp) + L.Q2x = Vector{Float64}(undef, n) + L.GHGxi = Vector{Float64}(undef, n) + L.HGxi = Vector{Float64}(undef, q) + L.x1 = Vector{Float64}(undef, n) + L.y1 = Vector{Float64}(undef, p) + L.z1 = Vector{Float64}(undef, q) + if useiterative + cgu = zeros(nmp) + L.cgstate = IterativeSolvers.CGStateVariables{Float64, Vector{Float64}}(cgu, similar(cgu), similar(cgu)) + L.lprecond = IterativeSolvers.Identity() + L.Q2sol = zeros(nmp) + end + + return L + end +end + +QRSymmCache( + c::Vector{Float64}, + A::AbstractMatrix{Float64}, + b::Vector{Float64}, + G::AbstractMatrix{Float64}, + h::Vector{Float64}, + cone::Cone; + useiterative::Bool = false, + ) = error("to use a QRSymmCache for linear system solves, the data must be preprocessed and Q2 and RiQ1 must be passed into the QRSymmCache constructor") + +# solve system for x, y, z +function solvelinsys3!( + rhs_tx::Vector{Float64}, + rhs_ty::Vector{Float64}, + rhs_tz::Vector{Float64}, + mu::Float64, + L::QRSymmCache; + identityH::Bool = false, + ) + + F = helplhs!(mu, L, identityH=identityH) + @. rhs_ty *= -1.0 + @. rhs_tz *= -1.0 + helplinsys!(rhs_tx, rhs_ty, rhs_tz, F, L) + + return nothing +end + +# solve system for x, y, z, s, kap, tau +function solvelinsys6!( + rhs_tx::Vector{Float64}, + rhs_ty::Vector{Float64}, + rhs_tz::Vector{Float64}, + rhs_ts::Vector{Float64}, + rhs_kap::Float64, + rhs_tau::Float64, + mu::Float64, + tau::Float64, + L::QRSymmCache, + ) + + # solve two symmetric systems and combine the solutions + F = helplhs!(mu, L) + + # (x2, y2, z2) = (rhs_tx, -rhs_ty, -H*rhs_ts - rhs_tz) + @. rhs_ty *= -1.0 + @. rhs_tz *= -1.0 + if !iszero(rhs_ts) + calcHarr!(L.z1, rhs_ts, L.cone) + @. rhs_tz -= mu*L.z1 + end + helplinsys!(rhs_tx, rhs_ty, rhs_tz, F, L) + + # (x1, y1, z1) = (-c, b, H*h) + @. L.x1 = -L.c + @. L.y1 = L.b + calcHarr!(L.z1, L.h, L.cone) + @. L.z1 *= mu + helplinsys!(L.x1, L.y1, L.z1, F, L) + + # combine + dir_tau = (rhs_tau + rhs_kap + dot(L.c, rhs_tx) + dot(L.b, rhs_ty) + dot(L.h, rhs_tz))/(mu/tau/tau - dot(L.c, L.x1) - dot(L.b, L.y1) - dot(L.h, L.z1)) + @. rhs_tx += dir_tau*L.x1 + @. rhs_ty += dir_tau*L.y1 + @. rhs_tz += dir_tau*L.z1 + mul!(L.z1, L.G, rhs_tx) + @. rhs_ts = -L.z1 + L.h*dir_tau - rhs_ts + dir_kap = -dot(L.c, rhs_tx) - dot(L.b, rhs_ty) - dot(L.h, rhs_tz) - rhs_tau + + return (dir_kap, dir_tau) +end + +# calculate solution to reduced symmetric linear system +function helplinsys!( + xi::Vector{Float64}, + yi::Vector{Float64}, + zi::Vector{Float64}, + F, + L::QRSymmCache, + ) + + # bxGHbz = bx + G'*Hbz + mul!(L.bxGHbz, L.G', zi) + @. L.bxGHbz += xi + # Q1x = Q1*Ri'*by + mul!(L.Q1x, L.RiQ1', yi) + # Q2x = Q2*(K22_F\(Q2'*(bxGHbz - GHG*Q1x))) + mul!(L.rhs, L.GHG, L.Q1x) + @. L.rhs = L.bxGHbz - L.rhs + mul!(L.Q2div, L.Q2', L.rhs) + + if L.useiterative + # TODO use previous solution from same pred/corr (requires knowing if in pred or corr step, and having two Q2sol objects) + # TODO allow tol to be loose in early iterations and tighter later (maybe depend on mu) + IterativeSolvers.cg!(L.Q2sol, F, L.Q2div, maxiter=10000, tol=1e-13, statevars=L.cgstate, Pl=L.lprecond, verbose=false) + @. L.Q2div = L.Q2sol + else + ldiv!(F, L.Q2div) + end + + mul!(L.Q2x, L.Q2, L.Q2div) + # xi = Q1x + Q2x + @. xi = L.Q1x + L.Q2x + # yi = Ri*Q1'*(bxGHbz - GHG*xi) + mul!(L.GHGxi, L.GHG, xi) + @. L.bxGHbz -= L.GHGxi + mul!(yi, L.RiQ1, L.bxGHbz) + # zi = HG*xi - Hbz + mul!(L.HGxi, L.HG, xi) + @. zi = L.HGxi - zi + + return nothing +end + +# calculate or factorize LHS of symmetric positive definite linear system +function helplhs!( + mu::Float64, + L::QRSymmCache; + identityH::Bool = false, + ) + + # Q2' * G' * H * G * Q2 + if identityH + @. L.HG = L.G + else + calcHarr!(L.HG, L.G, L.cone) + @. L.HG *= mu + end + mul!(L.GHG, L.G', L.HG) + mul!(L.GHGQ2, L.GHG, L.Q2) + mul!(L.Q2GHGQ2, L.Q2', L.GHGQ2) + + if L.useiterative + return Symmetric(L.Q2GHGQ2) + else + # bunchkaufman allocates more than cholesky, but doesn't fail when approximately quasidefinite (TODO could try LDL instead) + # TODO does it matter that F could be either type? + F = cholesky!(Symmetric(L.Q2GHGQ2), check=false) + if !issuccess(F) + # verbose && # TODO pass this in + println("linear system matrix was not positive definite") + mul!(L.Q2GHGQ2, L.Q2', L.GHGQ2) + F = bunchkaufman!(Symmetric(L.Q2GHGQ2)) + end + return F + end +end diff --git a/src/mathoptinterface.jl b/src/mathoptinterface.jl index 1ba476648..baf3f190b 100644 --- a/src/mathoptinterface.jl +++ b/src/mathoptinterface.jl @@ -2,8 +2,9 @@ Copyright 2018, Chris Coey and contributors =# -mutable struct Optimizer <: MOI.AbstractOptimizer - alf::HypatiaOpt +mutable struct HypatiaOptimizer <: MOI.AbstractOptimizer + opt::Optimizer + usedense::Bool c::Vector{Float64} # linear cost vector, size n A::AbstractMatrix{Float64} # equality constraint matrix, size p*n b::Vector{Float64} # equality constraint vector, size p @@ -26,22 +27,24 @@ mutable struct Optimizer <: MOI.AbstractOptimizer pobj::Float64 dobj::Float64 - function Optimizer(alf::HypatiaOpt) - opt = new() - opt.alf = alf - return opt + function HypatiaOptimizer(opt::Optimizer, verbose::Bool, usedense::Bool) + moiopt = new() + opt.verbose = verbose + moiopt.opt = opt + moiopt.usedense = usedense + return moiopt end end -Optimizer() = Optimizer(HypatiaOpt()) +HypatiaOptimizer(; verbose::Bool=false, usedense::Bool=true) = HypatiaOptimizer(Optimizer(), verbose, usedense) -MOI.get(::Optimizer, ::MOI.SolverName) = "Hypatia" +MOI.get(::HypatiaOptimizer, ::MOI.SolverName) = "Hypatia" -MOI.is_empty(opt::Optimizer) = (get_status(opt.alf) == :NotLoaded) -MOI.empty!(opt::Optimizer) = Optimizer() # TODO empty the data and results, or just create a new one? keep options? +MOI.is_empty(moiopt::HypatiaOptimizer) = (get_status(moiopt.opt) == :NotLoaded) +MOI.empty!(moiopt::HypatiaOptimizer) = HypatiaOptimizer() # TODO empty the data and results, or just create a new one? keep options? -MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true -MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true +MOI.supports(::HypatiaOptimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true +MOI.supports(::HypatiaOptimizer, ::MOI.ObjectiveSense) = true # TODO don't restrict to Float64 type const SupportedFuns = Union{ @@ -65,7 +68,7 @@ const SupportedSets = Union{ # MOI.PowerCone, } -MOI.supports_constraint(::Optimizer, ::Type{<:SupportedFuns}, ::Type{<:SupportedSets}) = true +MOI.supports_constraint(::HypatiaOptimizer, ::Type{<:SupportedFuns}, ::Type{<:SupportedSets}) = true conefrommoi(s::MOI.SecondOrderCone) = SecondOrderCone(MOI.dimension(s)) conefrommoi(s::MOI.RotatedSecondOrderCone) = RotatedSecondOrderCone(MOI.dimension(s)) @@ -75,7 +78,13 @@ conefrommoi(s::MOI.ExponentialCone) = ExponentialCone() # build representation as min c'x s.t. Ax = b, x in K # TODO what if some variables x are in multiple cone constraints? -function MOI.copy_to(opt::Optimizer, src::MOI.ModelLike; copy_names=false, warn_attributes=true) +function MOI.copy_to( + moiopt::HypatiaOptimizer, + src::MOI.ModelLike; + copy_names::Bool = false, + warn_attributes::Bool = true, + ) + @assert !copy_names idxmap = Dict{MOI.Index, MOI.Index}() @@ -99,9 +108,9 @@ function MOI.copy_to(opt::Optimizer, src::MOI.ModelLike; copy_names=false, warn_ if MOI.get(src, MOI.ObjectiveSense()) == MOI.MaxSense Vc .*= -1.0 end - opt.objconst = obj.constant - opt.objsense = MOI.get(src, MOI.ObjectiveSense()) - opt.c = Vector(sparsevec(Jc, Vc, n)) + moiopt.objconst = obj.constant + moiopt.objsense = MOI.get(src, MOI.ObjectiveSense()) + moiopt.c = Vector(sparsevec(Jc, Vc, n)) # constraints getsrccons(F, S) = MOI.get(src, MOI.ListOfConstraintIndices{F,S}()) @@ -113,7 +122,7 @@ function MOI.copy_to(opt::Optimizer, src::MOI.ModelLike; copy_names=false, warn_ p = 0 # rows of A (equality constraint matrix) (IA, JA, VA) = (Int[], Int[], Float64[]) (Ib, Vb) = (Int[], Float64[]) - (Icpe, Vcpe) = (Int[], Float64[]) # constraint set constants for opt.constrprimeq + (Icpe, Vcpe) = (Int[], Float64[]) # constraint set constants for moiopt.constrprimeq constroffseteq = Vector{Int}() for S in (MOI.EqualTo{Float64}, MOI.Zeros) @@ -166,6 +175,7 @@ function MOI.copy_to(opt::Optimizer, src::MOI.ModelLike; copy_names=false, warn_ append!(IA, p+1:p+dim) append!(JA, idxmap[vi].value for vi in fi.variables) append!(VA, -ones(dim)) + p += dim end for ci in getsrccons(MOI.VectorAffineFunction{Float64}, MOI.Zeros) @@ -185,116 +195,153 @@ function MOI.copy_to(opt::Optimizer, src::MOI.ModelLike; copy_names=false, warn_ end push!(constroffseteq, p) - opt.A = Matrix(sparse(IA, JA, VA, p, n)) - # opt.A = dropzeros!(sparse(IA, JA, VA, p, n)) - opt.b = Vector(sparsevec(Ib, Vb, p)) - opt.numeqconstrs = i - opt.constrprimeq = Vector(sparsevec(Icpe, Vcpe, p)) - opt.constroffseteq = constroffseteq + if moiopt.usedense + moiopt.A = Matrix(sparse(IA, JA, VA, p, n)) + else + moiopt.A = dropzeros!(sparse(IA, JA, VA, p, n)) + end + moiopt.b = Vector(sparsevec(Ib, Vb, p)) # TODO if type less strongly, this can be sparse too + moiopt.numeqconstrs = i + moiopt.constrprimeq = Vector(sparsevec(Icpe, Vcpe, p)) + moiopt.constroffseteq = constroffseteq # conic constraints q = 0 # rows of G (cone constraint matrix) (IG, JG, VG) = (Int[], Int[], Float64[]) (Ih, Vh) = (Int[], Float64[]) - (Icpc, Vcpc) = (Int[], Float64[]) # constraint set constants for opt.constrprimeq + (Icpc, Vcpc) = (Int[], Float64[]) # constraint set constants for moiopt.constrprimeq constroffsetcone = Vector{Int}() cone = Cone() - # LP constraints: build up one nonnegative cone and/or one nonpositive cone - nonnegrows = Int[] - nonposrows = Int[] + # build up one nonnegative cone + nonnegstart = q - # TODO also use variable bounds to build one L_inf cone + for ci in getsrccons(MOI.SingleVariable, MOI.GreaterThan{Float64}) + i += 1 + idxmap[ci] = MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}(i) + push!(constroffsetcone, q) + q += 1 + push!(IG, q) + push!(JG, idxmap[getconfun(ci).variable].value) + push!(VG, -1.0) + push!(Ih, q) + push!(Vh, -getconset(ci).lower) + push!(Vcpc, getconset(ci).lower) + push!(Icpc, q) + end - for S in (MOI.GreaterThan{Float64}, MOI.LessThan{Float64}) - for ci in getsrccons(MOI.SingleVariable, S) - i += 1 - idxmap[ci] = MOI.ConstraintIndex{MOI.SingleVariable, S}(i) - push!(constroffsetcone, q) + for ci in getsrccons(MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) + i += 1 + idxmap[ci] = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}(i) + push!(constroffsetcone, q) + q += 1 + fi = getconfun(ci) + for vt in fi.terms + push!(IG, q) + push!(JG, idxmap[vt.variable_index].value) + push!(VG, -vt.coefficient) + end + push!(Ih, q) + push!(Vh, fi.constant - getconset(ci).lower) + push!(Vcpc, getconset(ci).lower) + push!(Icpc, q) + end + + for ci in getsrccons(MOI.VectorOfVariables, MOI.Nonnegatives) + i += 1 + idxmap[ci] = MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}(i) + push!(constroffsetcone, q) + for vj in getconfun(ci).variables q += 1 push!(IG, q) - push!(JG, idxmap[getconfun(ci).variable].value) + push!(JG, idxmap[vj].value) push!(VG, -1.0) - push!(Ih, q) - if S == MOI.GreaterThan{Float64} - push!(Vh, -getconset(ci).lower) - push!(nonnegrows, q) - push!(Vcpc, getconset(ci).lower) - else - push!(Vh, -getconset(ci).upper) - push!(nonposrows, q) - push!(Vcpc, getconset(ci).upper) - end - push!(Icpc, q) end + end - for ci in getsrccons(MOI.ScalarAffineFunction{Float64}, S) - i += 1 - idxmap[ci] = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, S}(i) - push!(constroffsetcone, q) - q += 1 - fi = getconfun(ci) - for vt in fi.terms - push!(IG, q) - push!(JG, idxmap[vt.variable_index].value) - push!(VG, -vt.coefficient) - end - push!(Ih, q) - if S == MOI.GreaterThan{Float64} - push!(Vh, fi.constant - getconset(ci).lower) - push!(nonnegrows, q) - push!(Vcpc, getconset(ci).lower) - else - push!(Vh, fi.constant - getconset(ci).upper) - push!(nonposrows, q) - push!(Vcpc, getconset(ci).upper) - end - push!(Icpc, q) + for ci in getsrccons(MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) + i += 1 + idxmap[ci] = MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}(i) + push!(constroffsetcone, q) + fi = getconfun(ci) + dim = MOI.output_dimension(fi) + for vt in fi.terms + push!(IG, q + vt.output_index) + push!(JG, idxmap[vt.scalar_term.variable_index].value) + push!(VG, -vt.scalar_term.coefficient) end + append!(Ih, q+1:q+dim) + append!(Vh, fi.constants) + q += dim end - for S in (MOI.Nonnegatives, MOI.Nonpositives) - for ci in getsrccons(MOI.VectorOfVariables, S) - i += 1 - idxmap[ci] = MOI.ConstraintIndex{MOI.VectorOfVariables, S}(i) - push!(constroffsetcone, q) - for vj in getconfun(ci).variables - q += 1 - push!(IG, q) - push!(JG, idxmap[vj].value) - push!(VG, -1.0) - if S == MOI.Nonnegatives - push!(nonnegrows, q) - else - push!(nonposrows, q) - end - end + addprimitivecone!(cone, NonnegativeCone(q - nonnegstart), nonnegstart+1:q) + + # build up one nonpositive cone + nonposstart = q + + for ci in getsrccons(MOI.SingleVariable, MOI.LessThan{Float64}) + i += 1 + idxmap[ci] = MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}(i) + push!(constroffsetcone, q) + q += 1 + push!(IG, q) + push!(JG, idxmap[getconfun(ci).variable].value) + push!(VG, -1.0) + push!(Ih, q) + push!(Vh, -getconset(ci).upper) + push!(Vcpc, getconset(ci).upper) + push!(Icpc, q) + end + + for ci in getsrccons(MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) + i += 1 + idxmap[ci] = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}(i) + push!(constroffsetcone, q) + q += 1 + fi = getconfun(ci) + for vt in fi.terms + push!(IG, q) + push!(JG, idxmap[vt.variable_index].value) + push!(VG, -vt.coefficient) end + push!(Ih, q) + push!(Vh, fi.constant - getconset(ci).upper) + push!(Vcpc, getconset(ci).upper) + push!(Icpc, q) + end - for ci in getsrccons(MOI.VectorAffineFunction{Float64}, S) - i += 1 - idxmap[ci] = MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}, S}(i) - push!(constroffsetcone, q) - fi = getconfun(ci) - dim = MOI.output_dimension(fi) - for vt in fi.terms - push!(IG, q + vt.output_index) - push!(JG, idxmap[vt.scalar_term.variable_index].value) - push!(VG, -vt.scalar_term.coefficient) - end - append!(Ih, q+1:q+dim) - append!(Vh, fi.constants) - if S == MOI.Nonnegatives - append!(nonnegrows, q+1:q+dim) - else - append!(nonposrows, q+1:q+dim) - end - q += dim + for ci in getsrccons(MOI.VectorOfVariables, MOI.Nonpositives) + i += 1 + idxmap[ci] = MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonpositives}(i) + push!(constroffsetcone, q) + for vj in getconfun(ci).variables + q += 1 + push!(IG, q) + push!(JG, idxmap[vj].value) + push!(VG, -1.0) end end - addprimitivecone!(cone, NonnegativeCone(length(nonnegrows)), nonnegrows) - addprimitivecone!(cone, NonpositiveCone(length(nonposrows)), nonposrows) + for ci in getsrccons(MOI.VectorAffineFunction{Float64}, MOI.Nonpositives) + i += 1 + idxmap[ci] = MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}, MOI.Nonpositives}(i) + push!(constroffsetcone, q) + fi = getconfun(ci) + dim = MOI.output_dimension(fi) + for vt in fi.terms + push!(IG, q + vt.output_index) + push!(JG, idxmap[vt.scalar_term.variable_index].value) + push!(VG, -vt.scalar_term.coefficient) + end + append!(Ih, q+1:q+dim) + append!(Vh, fi.constants) + q += dim + end + + addprimitivecone!(cone, NonpositiveCone(q - nonposstart), nonposstart+1:q) + + # TODO also use variable bounds to build one L_inf cone # non-LP conic constraints @@ -331,160 +378,176 @@ function MOI.copy_to(opt::Optimizer, src::MOI.ModelLike; copy_names=false, warn_ end push!(constroffsetcone, q) - opt.G = Matrix(sparse(IG, JG, VG, q, n)) - # opt.G = dropzeros!(sparse(IG, JG, VG, q, n)) - opt.h = Vector(sparsevec(Ih, Vh, q)) - opt.cone = cone - opt.constroffsetcone = constroffsetcone - opt.constrprimcone = Vector(sparsevec(Icpc, Vcpc, q)) + if moiopt.usedense + moiopt.G = Matrix(sparse(IG, JG, VG, q, n)) + else + moiopt.G = dropzeros!(sparse(IG, JG, VG, q, n)) + end + moiopt.h = Vector(sparsevec(Ih, Vh, q)) # TODO if type less strongly, this can be sparse too + moiopt.cone = cone + moiopt.constroffsetcone = constroffsetcone + moiopt.constrprimcone = Vector(sparsevec(Icpc, Vcpc, q)) return idxmap end -function MOI.optimize!(opt::Optimizer) - Hypatia.load_data!(opt.alf, opt.c, opt.A, opt.b, opt.G, opt.h, opt.cone) # dense - Hypatia.solve!(opt.alf) +function MOI.optimize!(moiopt::HypatiaOptimizer) + opt = moiopt.opt + (c, A, b, G, h, cone) = (moiopt.c, moiopt.A, moiopt.b, moiopt.G, moiopt.h, moiopt.cone) + + check_data(c, A, b, G, h, cone) + + # TODO make it optional + (c1, A1, b1, G1, prkeep, dukeep, Q2, RiQ1) = preprocess_data(c, A, b, G, useQR=true) + + L = QRSymmCache(c1, A1, b1, G1, h, cone, Q2, RiQ1) + load_data!(opt, c1, A1, b1, G1, h, cone, L) + + solve!(opt) + + moiopt.x = zeros(length(c)) + moiopt.x[dukeep] = get_x(opt) + moiopt.constrprimeq += moiopt.b - moiopt.A*moiopt.x + moiopt.y = zeros(length(b)) + moiopt.y[prkeep] = get_y(opt) - opt.x = Hypatia.get_x(opt.alf) - opt.constrprimeq += opt.b - opt.A*opt.x - opt.s = Hypatia.get_s(opt.alf) - opt.constrprimcone += opt.s - opt.y = Hypatia.get_y(opt.alf) - opt.z = Hypatia.get_z(opt.alf) + moiopt.s = get_s(opt) + moiopt.constrprimcone += moiopt.s + moiopt.z = get_z(opt) - opt.status = Hypatia.get_status(opt.alf) - opt.solvetime = Hypatia.get_solvetime(opt.alf) - opt.pobj = Hypatia.get_pobj(opt.alf) - opt.dobj = Hypatia.get_dobj(opt.alf) + moiopt.status = get_status(opt) + moiopt.solvetime = get_solvetime(opt) + moiopt.pobj = get_pobj(opt) + moiopt.dobj = get_dobj(opt) return nothing end -# function MOI.free!(opt::Optimizer) # TODO call gc on opt.alf? +# function MOI.free!(moiopt::HypatiaOptimizer) # TODO call gc on moiopt.opt? -function MOI.get(opt::Optimizer, ::MOI.TerminationStatus) +function MOI.get(moiopt::HypatiaOptimizer, ::MOI.TerminationStatus) # TODO time limit etc - if opt.status in (:Optimal, :PrimalInfeasible, :DualInfeasible, :IllPosed) + if moiopt.status in (:Optimal, :PrimalInfeasible, :DualInfeasible, :IllPosed) return MOI.Success - elseif opt.status in (:PredictorFail, :CorrectorFail) + elseif moiopt.status in (:PredictorFail, :CorrectorFail) return MOI.NumericalError - elseif opt.status == :IterationLimit + elseif moiopt.status == :IterationLimit return MOI.IterationLimit else return MOI.OtherError end end -function MOI.get(opt::Optimizer, ::MOI.ObjectiveValue) - if opt.objsense == MOI.MinSense - return opt.pobj + opt.objconst - elseif opt.objsense == MOI.MaxSense - return -opt.pobj + opt.objconst +function MOI.get(moiopt::HypatiaOptimizer, ::MOI.ObjectiveValue) + if moiopt.objsense == MOI.MinSense + return moiopt.pobj + moiopt.objconst + elseif moiopt.objsense == MOI.MaxSense + return -moiopt.pobj + moiopt.objconst else error("no objective sense is set") end end -function MOI.get(opt::Optimizer, ::MOI.ObjectiveBound) - if opt.objsense == MOI.MinSense - return opt.dobj + opt.objconst - elseif opt.objsense == MOI.MaxSense - return -opt.dobj + opt.objconst +function MOI.get(moiopt::HypatiaOptimizer, ::MOI.ObjectiveBound) + if moiopt.objsense == MOI.MinSense + return moiopt.dobj + moiopt.objconst + elseif moiopt.objsense == MOI.MaxSense + return -moiopt.dobj + moiopt.objconst else error("no objective sense is set") end end -function MOI.get(opt::Optimizer, ::MOI.ResultCount) - if opt.status in (:Optimal, :PrimalInfeasible, :DualInfeasible, :IllPosed) +function MOI.get(moiopt::HypatiaOptimizer, ::MOI.ResultCount) + if moiopt.status in (:Optimal, :PrimalInfeasible, :DualInfeasible, :IllPosed) return 1 end return 0 end -function MOI.get(opt::Optimizer, ::MOI.PrimalStatus) - if opt.status == :Optimal +function MOI.get(moiopt::HypatiaOptimizer, ::MOI.PrimalStatus) + if moiopt.status == :Optimal return MOI.FeasiblePoint - elseif opt.status == :PrimalInfeasible + elseif moiopt.status == :PrimalInfeasible return MOI.InfeasiblePoint - elseif opt.status == :DualInfeasible + elseif moiopt.status == :DualInfeasible return MOI.InfeasibilityCertificate - elseif opt.status == :IllPosed + elseif moiopt.status == :IllPosed return MOI.UnknownResultStatus # TODO later distinguish primal/dual ill posed certificates else return MOI.UnknownResultStatus end end -function MOI.get(opt::Optimizer, ::MOI.DualStatus) - if opt.status == :Optimal +function MOI.get(moiopt::HypatiaOptimizer, ::MOI.DualStatus) + if moiopt.status == :Optimal return MOI.FeasiblePoint - elseif opt.status == :PrimalInfeasible + elseif moiopt.status == :PrimalInfeasible return MOI.InfeasibilityCertificate - elseif opt.status == :DualInfeasible + elseif moiopt.status == :DualInfeasible return MOI.InfeasiblePoint - elseif opt.status == :IllPosed + elseif moiopt.status == :IllPosed return MOI.UnknownResultStatus # TODO later distinguish primal/dual ill posed certificates else return MOI.UnknownResultStatus end end -MOI.get(opt::Optimizer, ::MOI.VariablePrimal, vi::MOI.VariableIndex) = opt.x[vi.value] -MOI.get(opt::Optimizer, a::MOI.VariablePrimal, vi::Vector{MOI.VariableIndex}) = MOI.get.(opt, a, vi) +MOI.get(moiopt::HypatiaOptimizer, ::MOI.VariablePrimal, vi::MOI.VariableIndex) = moiopt.x[vi.value] +MOI.get(moiopt::HypatiaOptimizer, a::MOI.VariablePrimal, vi::Vector{MOI.VariableIndex}) = MOI.get.(moiopt, a, vi) -function MOI.get(opt::Optimizer, ::MOI.ConstraintDual, ci::MOI.ConstraintIndex{F, S}) where {F <: MOI.AbstractFunction, S <: MOI.AbstractScalarSet} +function MOI.get(moiopt::HypatiaOptimizer, ::MOI.ConstraintDual, ci::MOI.ConstraintIndex{F, S}) where {F <: MOI.AbstractFunction, S <: MOI.AbstractScalarSet} # scalar set i = ci.value - if i <= opt.numeqconstrs + if i <= moiopt.numeqconstrs # constraint is an equality - return opt.y[opt.constroffseteq[i]+1] + return moiopt.y[moiopt.constroffseteq[i]+1] else # constraint is conic - i -= opt.numeqconstrs - return opt.z[opt.constroffsetcone[i]+1] + i -= moiopt.numeqconstrs + return moiopt.z[moiopt.constroffsetcone[i]+1] end end -function MOI.get(opt::Optimizer, ::MOI.ConstraintDual, ci::MOI.ConstraintIndex{F, S}) where {F <: MOI.AbstractFunction, S <: MOI.AbstractVectorSet} +function MOI.get(moiopt::HypatiaOptimizer, ::MOI.ConstraintDual, ci::MOI.ConstraintIndex{F, S}) where {F <: MOI.AbstractFunction, S <: MOI.AbstractVectorSet} # vector set i = ci.value - if i <= opt.numeqconstrs + if i <= moiopt.numeqconstrs # constraint is an equality - os = opt.constroffseteq - return opt.y[os[i]+1:os[i+1]] + os = moiopt.constroffseteq + return moiopt.y[os[i]+1:os[i+1]] else # constraint is conic - i -= opt.numeqconstrs - os = opt.constroffsetcone - return opt.z[os[i]+1:os[i+1]] + i -= moiopt.numeqconstrs + os = moiopt.constroffsetcone + return moiopt.z[os[i]+1:os[i+1]] end end -MOI.get(opt::Optimizer, a::MOI.ConstraintDual, ci::Vector{MOI.ConstraintIndex}) = MOI.get.(opt, a, ci) +MOI.get(moiopt::HypatiaOptimizer, a::MOI.ConstraintDual, ci::Vector{MOI.ConstraintIndex}) = MOI.get.(moiopt, a, ci) -function MOI.get(opt::Optimizer, ::MOI.ConstraintPrimal, ci::MOI.ConstraintIndex{F, S}) where {F <: MOI.AbstractFunction, S <: MOI.AbstractScalarSet} +function MOI.get(moiopt::HypatiaOptimizer, ::MOI.ConstraintPrimal, ci::MOI.ConstraintIndex{F, S}) where {F <: MOI.AbstractFunction, S <: MOI.AbstractScalarSet} # scalar set i = ci.value - if i <= opt.numeqconstrs + if i <= moiopt.numeqconstrs # constraint is an equality - return opt.constrprimeq[opt.constroffseteq[i]+1] + return moiopt.constrprimeq[moiopt.constroffseteq[i]+1] else # constraint is conic - i -= opt.numeqconstrs - return opt.constrprimcone[opt.constroffsetcone[i]+1] + i -= moiopt.numeqconstrs + return moiopt.constrprimcone[moiopt.constroffsetcone[i]+1] end end -function MOI.get(opt::Optimizer, ::MOI.ConstraintPrimal, ci::MOI.ConstraintIndex{F, S}) where {F <: MOI.AbstractFunction, S <: MOI.AbstractVectorSet} +function MOI.get(moiopt::HypatiaOptimizer, ::MOI.ConstraintPrimal, ci::MOI.ConstraintIndex{F, S}) where {F <: MOI.AbstractFunction, S <: MOI.AbstractVectorSet} # vector set i = ci.value - if i <= opt.numeqconstrs + if i <= moiopt.numeqconstrs # constraint is an equality - os = opt.constroffseteq - return opt.constrprimeq[os[i]+1:os[i+1]] + os = moiopt.constroffseteq + return moiopt.constrprimeq[os[i]+1:os[i+1]] else # constraint is conic - i -= opt.numeqconstrs - os = opt.constroffsetcone - return opt.constrprimcone[os[i]+1:os[i+1]] + i -= moiopt.numeqconstrs + os = moiopt.constroffsetcone + return moiopt.constrprimcone[os[i]+1:os[i+1]] end end -MOI.get(opt::Optimizer, a::MOI.ConstraintPrimal, ci::Vector{MOI.ConstraintIndex}) = MOI.get.(opt, a, ci) +MOI.get(moiopt::HypatiaOptimizer, a::MOI.ConstraintPrimal, ci::Vector{MOI.ConstraintIndex}) = MOI.get.(moiopt, a, ci) diff --git a/src/nativeinterface.jl b/src/nativeinterface.jl index 456064592..ad4d2a70e 100644 --- a/src/nativeinterface.jl +++ b/src/nativeinterface.jl @@ -3,59 +3,8 @@ Copyright 2018, Chris Coey and contributors Copyright 2018, David Papp, Sercan Yildiz =# - -# cache for linear system solves (to avoid allocations) -# use QR + cholesky method from CVXOPT -# (1) eliminate equality constraints via QR of A' -# (2) solve reduced system by cholesky -# |0 A' G' | * |ux| = |bx| -# |A 0 0 | |uy| |by| -# |G 0 -Hi/mu| |uz| |bz| -mutable struct LinSysCache - Q2 - RiQ1 - GQ2 - HG - GHG - GHGQ2 - Q2GHGQ2 - bxGHbz - Q1x - rhs - Q2div - Q2x - GHGxi - HGxi - x1 - y1 - z1 - - function LinSysCache(Q1::AbstractMatrix{Float64}, Q2::AbstractMatrix{Float64}, Ri::AbstractMatrix{Float64}, G::AbstractMatrix{Float64}, n::Int, p::Int, q::Int) - L = new() - nmp = n - p - L.Q2 = Q2 - L.RiQ1 = Ri*Q1' - L.GQ2 = G*Q2 - L.HG = Matrix{Float64}(undef, q, n) # TODO don't enforce dense on some - L.GHG = Matrix{Float64}(undef, n, n) - L.GHGQ2 = Matrix{Float64}(undef, n, nmp) - L.Q2GHGQ2 = Matrix{Float64}(undef, nmp, nmp) - L.bxGHbz = Vector{Float64}(undef, n) - L.Q1x = Vector{Float64}(undef, n) - L.rhs = Vector{Float64}(undef, n) - L.Q2div = Vector{Float64}(undef, nmp) - L.Q2x = Vector{Float64}(undef, n) - L.GHGxi = Vector{Float64}(undef, n) - L.HGxi = Vector{Float64}(undef, q) - L.x1 = Vector{Float64}(undef, n) - L.y1 = Vector{Float64}(undef, p) - L.z1 = Vector{Float64}(undef, q) - return L - end -end - # model object containing options, problem data, linear system cache, and solution -mutable struct HypatiaOpt +mutable struct Optimizer # options verbose::Bool # if true, prints progress at each iteration tolrelopt::Float64 # relative optimality gap tolerance @@ -96,31 +45,31 @@ mutable struct HypatiaOpt pobj::Float64 # final primal objective value dobj::Float64 # final dual objective value - function HypatiaOpt(verbose, tolrelopt, tolabsopt, tolfeas, maxiter, predlinesearch, maxpredsmallsteps, predlsmulti, corrcheck, maxcorrsteps, alphacorr, maxcorrlsiters, corrlsmulti) - alf = new() - - alf.verbose = verbose - alf.tolrelopt = tolrelopt - alf.tolabsopt = tolabsopt - alf.tolfeas = tolfeas - alf.maxiter = maxiter - alf.predlinesearch = predlinesearch - alf.maxpredsmallsteps = maxpredsmallsteps - alf.predlsmulti = predlsmulti - alf.corrcheck = corrcheck - alf.maxcorrsteps = maxcorrsteps - alf.alphacorr = alphacorr - alf.maxcorrlsiters = maxcorrlsiters - alf.corrlsmulti = corrlsmulti - - alf.status = :NotLoaded - - return alf + function Optimizer(verbose, tolrelopt, tolabsopt, tolfeas, maxiter, predlinesearch, maxpredsmallsteps, predlsmulti, corrcheck, maxcorrsteps, alphacorr, maxcorrlsiters, corrlsmulti) + opt = new() + + opt.verbose = verbose + opt.tolrelopt = tolrelopt + opt.tolabsopt = tolabsopt + opt.tolfeas = tolfeas + opt.maxiter = maxiter + opt.predlinesearch = predlinesearch + opt.maxpredsmallsteps = maxpredsmallsteps + opt.predlsmulti = predlsmulti + opt.corrcheck = corrcheck + opt.maxcorrsteps = maxcorrsteps + opt.alphacorr = alphacorr + opt.maxcorrlsiters = maxcorrlsiters + opt.corrlsmulti = corrlsmulti + + opt.status = :NotLoaded + + return opt end end # initialize a model object -function HypatiaOpt(; +function Optimizer(; verbose = false, tolrelopt = 1e-6, tolabsopt = 1e-7, @@ -149,43 +98,41 @@ function HypatiaOpt(; error("maxcorrsteps must be at least 1") end - return HypatiaOpt(verbose, tolrelopt, tolabsopt, tolfeas, maxiter, predlinesearch, maxpredsmallsteps, predlsmulti, corrcheck, maxcorrsteps, alphacorr, maxcorrlsiters, corrlsmulti) + return Optimizer(verbose, tolrelopt, tolabsopt, tolfeas, maxiter, predlinesearch, maxpredsmallsteps, predlsmulti, corrcheck, maxcorrsteps, alphacorr, maxcorrlsiters, corrlsmulti) end -get_status(alf::HypatiaOpt) = alf.status -get_solvetime(alf::HypatiaOpt) = alf.solvetime -get_niters(alf::HypatiaOpt) = alf.niters +get_status(opt::Optimizer) = opt.status +get_solvetime(opt::Optimizer) = opt.solvetime +get_niters(opt::Optimizer) = opt.niters -get_x(alf::HypatiaOpt) = copy(alf.x) -get_s(alf::HypatiaOpt) = copy(alf.s) -get_y(alf::HypatiaOpt) = copy(alf.y) -get_z(alf::HypatiaOpt) = copy(alf.z) +get_x(opt::Optimizer) = copy(opt.x) +get_s(opt::Optimizer) = copy(opt.s) +get_y(opt::Optimizer) = copy(opt.y) +get_z(opt::Optimizer) = copy(opt.z) -get_tau(alf::HypatiaOpt) = alf.tau -get_kappa(alf::HypatiaOpt) = alf.kappa -get_mu(alf::HypatiaOpt) = alf.mu +get_tau(opt::Optimizer) = opt.tau +get_kappa(opt::Optimizer) = opt.kappa +get_mu(opt::Optimizer) = opt.mu -get_pobj(alf::HypatiaOpt) = dot(alf.c, alf.x) -get_dobj(alf::HypatiaOpt) = -dot(alf.b, alf.y) - dot(alf.h, alf.z) +get_pobj(opt::Optimizer) = dot(opt.c, opt.x) +get_dobj(opt::Optimizer) = -dot(opt.b, opt.y) - dot(opt.h, opt.z) -# verify problem data and load into model object -function load_data!( - alf::HypatiaOpt, +# check data for consistency +function check_data( c::Vector{Float64}, A::AbstractMatrix{Float64}, b::Vector{Float64}, G::AbstractMatrix{Float64}, h::Vector{Float64}, - cone::Cone; - check::Bool=false, # check rank conditions + cone::Cone, ) - # check data consistency - n = length(c) - p = length(b) - q = length(h) + (n, p, q) = (length(c), length(b), length(h)) @assert n > 0 @assert p + q > 0 + if n < p + println("number of equality constraints ($p) exceeds number of variables ($n)") + end if n != size(A, 2) || n != size(G, 2) error("number of variables is not consistent in A, G, and c") end @@ -201,63 +148,158 @@ function load_data!( @assert dimension(cone.prms[k]) == length(cone.idxs[k]) end - # perform QR decomposition of A' for use in linear system solves - # TODO reduce allocs, improve efficiency - # A' = [Q1 Q2] * [R1; 0] - if issparse(A) - alf.verbose && println("\nJulia is currently missing some sparse matrix methods that could improve performance; Hypatia may perform better if A is loaded as a dense matrix") - # TODO currently using dense Q1, Q2, R - probably some should be sparse - F = qr(sparse(A')) - @assert length(F.prow) == n - @assert length(F.pcol) == p - @assert istriu(F.R) - - Q = F.Q*Matrix(1.0I, n, n) - Q1 = zeros(n, p) - Q1[F.prow, F.pcol] = Q[:, 1:p] - Q2 = zeros(n, n-p) - Q2[F.prow, :] = Q[:, p+1:n] - Ri = zeros(p, p) - Ri[F.pcol, F.pcol] = inv(UpperTriangular(F.R)) + return nothing +end + +# preprocess data (optional) +function preprocess_data( + c::Vector{Float64}, + A::AbstractMatrix{Float64}, + b::Vector{Float64}, + G::AbstractMatrix{Float64}; + tol::Float64 = 1e-13, # presolve tolerance + useQR::Bool = false, # returns QR fact of A' for use in a QR-based linear system solver + ) + + (n, p) = (length(c), length(b)) + + # NOTE (pivoted) QR factorizations are usually rank-revealing but may be unreliable, see http://www.math.sjsu.edu/~foster/rankrevealingcode.html + # rank of a matrix is number of nonzero diagonal elements of R + + # preprocess dual equality constraints + dukeep = 1:n + AG = vcat(A, G) + + # get pivoted QR # TODO when Julia has a unified QR interface, replace this + if issparse(AG) + AGF = qr(AG, tol=tol) + AGR = AGF.R + AGrank = rank(AGF) # cheap else - F = qr(A') - @assert istriu(F.R) - - Q = F.Q*Matrix(1.0I, n, n) - Q1 = Q[:, 1:p] - Q2 = Q[:, p+1:n] - Ri = inv(UpperTriangular(F.R)) - @assert norm(A'*Ri - Q1) < 1e-8 # TODO delete later + AGF = qr(AG, Val(true)) + AGR = AGF.R + AGrank = 0 + for i in 1:size(AGR, 1) # TODO could replace this with rank(AF) when available for both dense and sparse + if abs(AGR[i,i]) > tol + AGrank += 1 + end + end end - if check - # check rank conditions - # TODO rank for qr decomp should be implemented in Julia - see https://github.com/JuliaLang/julia/blob/f8b52dab77415a22d28497f48407aca92fbbd4c3/stdlib/LinearAlgebra/src/qr.jl#L895 - if rank(A) < p # TODO change to rank(F) - error("A matrix is not full-row-rank; some primal equalities may be redundant or inconsistent") + if AGrank < n + if issparse(AG) + dukeep = AGF.pcol[1:AGrank] + AGQ1 = Matrix{Float64}(undef, n, AGrank) + AGQ1[AGF.prow,:] = AGF.Q*Matrix{Float64}(I, n, AGrank) # TODO could eliminate this allocation + else + dukeep = AGF.p[1:AGrank] + AGQ1 = AGF.Q*Matrix{Float64}(I, n, AGrank) # TODO could eliminate this allocation end - if rank(vcat(A, G)) < n - error("[A' G'] is not full-row-rank; some dual equalities may be redundant (i.e. primal variables can be removed) or inconsistent") + AGRiQ1 = UpperTriangular(AGR[1:AGrank,1:AGrank])\AGQ1' + + A1 = A[:,dukeep] + G1 = G[:,dukeep] + c1 = c[dukeep] + + if norm(AG'*AGRiQ1'*c1 - c, Inf) > tol + error("some dual equality constraints are inconsistent") end + + A = A1 + G = G1 + c = c1 + println("removed $(n - AGrank) out of $n dual equality constraints") + n = AGrank end - alf.c = c - alf.A = A - alf.b = b - alf.G = G - alf.h = h - alf.cone = cone - alf.L = LinSysCache(Q1, Q2, Ri, G, n, p, q) - alf.status = :Loaded + if p == 0 + # no primal equality constraints to preprocess + # TODO use I instead of dense for Q2 + return (c, A, b, G, 1:0, dukeep, Matrix{Float64}(I, n, n), Matrix{Float64}(I, 0, n)) + end + + # preprocess primal equality constraints + # get pivoted QR # TODO when Julia has a unified QR interface, replace this + # TODO LQ decomposition is the QR decomposition of A' + if issparse(A) + AF = qr(sparse(A'), tol=tol) + AR = AF.R + Arank = rank(AF) # cheap + else + AF = qr(A', Val(true)) + AR = AF.R + Arank = 0 + for i in 1:size(AR, 1) # TODO could replace this with rank(AF) when available for both dense and sparse + if abs(AR[i,i]) > tol + Arank += 1 + end + end + end + + if !useQR && Arank == p + # no primal equalities to remove and QR of A' not needed + return (c, A, b, G, 1:p, dukeep, Matrix{Float64}(undef, 0, 0), Matrix{Float64}(undef, 0, 0)) + end + + # using QR of A' (requires reordering rows) and/or some primal equalities are dependent + if issparse(A) + prkeep = AF.pcol[1:Arank] + AQ = Matrix{Float64}(undef, n, n) + AQ[AF.prow,:] = AF.Q*Matrix{Float64}(I, n, n) # TODO could eliminate this allocation + else + prkeep = AF.p[1:Arank] + AQ = AF.Q*Matrix{Float64}(I, n, n) # TODO could eliminate this allocation + end + AQ2 = AQ[:,Arank+1:n] + ARiQ1 = UpperTriangular(AR[1:Arank,1:Arank])\AQ[:,1:Arank]' + + A1 = A[prkeep,:] + b1 = b[prkeep] + + if Arank < p + # some dependent primal equalities, so check if they are consistent + x1 = ARiQ1'*b1 + if norm(A*x1 - b, Inf) > tol + error("some primal equality constraints are inconsistent") + end + println("removed $(p - Arank) out of $p primal equality constraints") + end + + A = A1 + b = b1 + + return (c, A, b, G, prkeep, dukeep, AQ2, ARiQ1) +end + +# verify problem data and load into model object +function load_data!( + opt::Optimizer, + c::Vector{Float64}, + A::AbstractMatrix{Float64}, + b::Vector{Float64}, + G::AbstractMatrix{Float64}, + h::Vector{Float64}, + cone::Cone, + L::LinSysCache, # linear system solver cache (see linsyssolvers folder) + ) + + opt.c = c + opt.A = A + opt.b = b + opt.G = G + opt.h = h + opt.cone = cone + opt.L = L + opt.status = :Loaded - return alf + return opt end # solve using predictor-corrector algorithm based on homogeneous self-dual embedding -function solve!(alf::HypatiaOpt) +function solve!(opt::Optimizer) starttime = time() - (c, A, b, G, h, cone) = (alf.c, alf.A, alf.b, alf.G, alf.h, alf.cone) + (c, A, b, G, h, cone) = (opt.c, opt.A, opt.b, opt.G, opt.h, opt.cone) (n, p, q) = (length(c), length(b), length(h)) bnu = 1.0 + barrierpar(cone) # complexity parameter nu-bar of the augmented barrier (sum of the primitive cone barrier parameters plus 1) @@ -282,7 +324,47 @@ function solve!(alf::HypatiaOpt) tmp_ts = similar(ts) # find initial primal-dual iterate - (tau, kap, mu) = findinitialiterate!(tx, ty, tz, ts, ls_ts, bnu, alf) + opt.verbose && println("\nfinding initial iterate") + + # solve linear equation + # TODO check equations + # |0 A' G'| * |tx| = |0| + # |A 0 0 | |ty| |b| + # |G 0 -I | |ts| |h| + @. tx = 0.0 + @. ty = -b + @. tz = -h + solvelinsys3!(tx, ty, tz, 1.0, opt.L, identityH=true) + @. ts = -tz + + # from ts, step along interior direction of cone + getintdir!(tmp_ts, cone) + alpha = (norm(ts) + 1e-3)#/norm(tmp_ts) # TODO tune this + @. ls_ts = ts + alpha*tmp_ts + + # continue stepping until ts is inside cone + steps = 1 + while !incone(cone) + steps += 1 + if steps > 25 + error("cannot find initial iterate") + end + alpha *= 1.5 + @. ls_ts = ts + alpha*tmp_ts + end + opt.verbose && println("$steps steps taken for initial iterate") + @. ts = ls_ts + + calcg!(tz, cone) + @. tz *= -1.0 + tau = 1.0 + kap = 1.0 + mu = (dot(tz, ts) + tau*kap)/bnu + @assert !isnan(mu) + @assert abs(1.0 - mu) < 1e-10 + # @assert calcnbhd(tau*kap, mu, copy(tz), copy(tz), cone) < 1e-6 + + opt.verbose && println("initial iterate found") # calculate tolerances for convergence tol_res_tx = inv(max(1.0, norm(c))) @@ -290,19 +372,19 @@ function solve!(alf::HypatiaOpt) tol_res_tz = inv(max(1.0, norm(h))) # calculate prediction and correction step parameters - (beta, eta, cpredfix) = getbetaeta(alf.maxcorrsteps, bnu) # beta: large neighborhood parameter, eta: small neighborhood parameter + (beta, eta, cpredfix) = getbetaeta(opt.maxcorrsteps, bnu) # beta: large neighborhood parameter, eta: small neighborhood parameter alphapredfix = cpredfix/(eta + sqrt(2*eta^2 + bnu)) # fixed predictor step size - alphapredthres = (alf.predlsmulti^alf.maxpredsmallsteps)*alphapredfix # minimum predictor step size - alphapredinit = (alf.predlinesearch ? min(100*alphapredfix, 0.9999) : alphapredfix) # predictor step size + alphapredthres = (opt.predlsmulti^opt.maxpredsmallsteps)*alphapredfix # minimum predictor step size + alphapredinit = (opt.predlinesearch ? min(100*alphapredfix, 0.9999) : alphapredfix) # predictor step size # main loop - if alf.verbose + if opt.verbose println("starting iteration") @printf("\n%5s %12s %12s %9s %9s %9s %9s %9s %9s\n", "iter", "p_obj", "d_obj", "rel_gap", "p_inf", "d_inf", "tau", "kap", "mu") flush(stdout) end - alf.status = :StartedIterating + opt.status = :StartedIterating alphapred = alphapredinit iter = 0 while true @@ -333,7 +415,7 @@ function solve!(alf::HypatiaOpt) (cx, by, hz) = (dot(c, tx), dot(b, ty), dot(h, tz)) obj_pr = cx*invtau obj_du = -(by + hz)*invtau - gap = dot(tz, ts) # TODO maybe should adapt original alfonso condition instead of using this CVXOPT condition + gap = dot(tz, ts) # TODO maybe should adapt original Alfonso condition instead of using this CVXOPT condition # TODO maybe add small epsilon to denominators that are zero to avoid NaNs, and get rid of isnans further down if obj_pr < 0.0 @@ -358,7 +440,7 @@ function solve!(alf::HypatiaOpt) infres_du = NaN end - if alf.verbose + if opt.verbose # print iteration statistics @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e\n", iter, obj_pr, obj_du, relgap, nres_pr, nres_du, tau, kap, mu) flush(stdout) @@ -366,51 +448,51 @@ function solve!(alf::HypatiaOpt) # check convergence criteria # TODO nearly primal or dual infeasible or nearly optimal cases? - if nres_pr <= alf.tolfeas && nres_du <= alf.tolfeas && (gap <= alf.tolabsopt || (!isnan(relgap) && relgap <= alf.tolrelopt)) - alf.verbose && println("optimal solution found; terminating") - alf.status = :Optimal - - alf.x = tx .*= invtau - alf.s = ts .*= invtau - alf.y = ty .*= invtau - alf.z = tz .*= invtau + if nres_pr <= opt.tolfeas && nres_du <= opt.tolfeas && (gap <= opt.tolabsopt || (!isnan(relgap) && relgap <= opt.tolrelopt)) + opt.verbose && println("optimal solution found; terminating") + opt.status = :Optimal + opt.x = tx .*= invtau + opt.s = ts .*= invtau + opt.y = ty .*= invtau + opt.z = tz .*= invtau break - elseif !isnan(infres_pr) && infres_pr <= alf.tolfeas - alf.verbose && println("primal infeasibility detected; terminating") - alf.status = :PrimalInfeasible - + elseif !isnan(infres_pr) && infres_pr <= opt.tolfeas + opt.verbose && println("primal infeasibility detected; terminating") + opt.status = :PrimalInfeasible invobj = inv(-by - hz) - alf.x = tx .= NaN - alf.s = ts .= NaN - alf.y = ty .*= invobj - alf.z = tz .*= invobj + opt.x = tx .= NaN + opt.s = ts .= NaN + opt.y = ty .*= invobj + opt.z = tz .*= invobj break - elseif !isnan(infres_du) && infres_du <= alf.tolfeas - alf.verbose && println("dual infeasibility detected; terminating") - alf.status = :DualInfeasible - + elseif !isnan(infres_du) && infres_du <= opt.tolfeas + opt.verbose && println("dual infeasibility detected; terminating") + opt.status = :DualInfeasible invobj = inv(-cx) - alf.x = tx .*= invobj - alf.s = ts .*= invobj - alf.y = ty .= NaN - alf.z = tz .= NaN + opt.x = tx .*= invobj + opt.s = ts .*= invobj + opt.y = ty .= NaN + opt.z = tz .= NaN break - elseif mu <= alf.tolfeas*1e-2 && tau <= alf.tolfeas*1e-2*min(1.0, kap) - alf.verbose && println("ill-posedness detected; terminating") - alf.status = :IllPosed - - alf.x = tx .= NaN - alf.s = ts .= NaN - alf.y = ty .= NaN - alf.z = tz .= NaN + elseif mu <= opt.tolfeas*1e-2 && tau <= opt.tolfeas*1e-2*min(1.0, kap) + opt.verbose && println("ill-posedness detected; terminating") + opt.status = :IllPosed + opt.x = tx .= NaN + opt.s = ts .= NaN + opt.y = ty .= NaN + opt.z = tz .= NaN break end # check iteration limit iter += 1 - if iter >= alf.maxiter - alf.verbose && println("iteration limit reached; terminating") - alf.status = :IterationLimit + if iter >= opt.maxiter + opt.verbose && println("iteration limit reached; terminating") + opt.status = :IterationLimit + opt.x = tx .*= invtau + opt.s = ts .*= invtau + opt.y = ty .*= invtau + opt.z = tz .*= invtau break end @@ -418,7 +500,7 @@ function solve!(alf::HypatiaOpt) # calculate prediction direction @. tmp_ts = tmp_tz @. tmp_tz = -tz - (tmp_kap, tmp_tau) = finddirection!(tmp_tx, tmp_ty, tmp_tz, tmp_ts, -kap, kap + cx + by + hz, mu, tau, alf) + (tmp_kap, tmp_tau) = solvelinsys6!(tmp_tx, tmp_ty, tmp_tz, tmp_ts, -kap, kap + cx + by + hz, mu, tau, opt.L) # determine step length alpha by line search alpha = alphapred @@ -444,16 +526,16 @@ function solve!(alf::HypatiaOpt) if nbhd < abs2(beta*ls_mu) # iterate is inside the beta-neighborhood - if !alphaprevok || (alpha > alf.predlsmulti) + if !alphaprevok || (alpha > opt.predlsmulti) # either the previous iterate was outside the beta-neighborhood or increasing alpha again will make it > 1 - if alf.predlinesearch + if opt.predlinesearch alphapred = alpha end break end alphaprevok = true - alpha = alpha/alf.predlsmulti # increase alpha + alpha = alpha/opt.predlsmulti # increase alpha continue end end @@ -464,15 +546,19 @@ function solve!(alf::HypatiaOpt) if alpha < alphapredthres # alpha is very small, so predictor has failed predfail = true - alf.verbose && println("predictor could not improve the solution ($nprediters line search steps); terminating") - alf.status = :PredictorFail + opt.verbose && println("predictor could not improve the solution ($nprediters line search steps); terminating") break end alphaprevok = false - alpha = alf.predlsmulti*alpha # decrease alpha + alpha = opt.predlsmulti*alpha # decrease alpha end if predfail + opt.status = :PredictorFail + opt.x = tx .*= invtau + opt.s = ts .*= invtau + opt.y = ty .*= invtau + opt.z = tz .*= invtau break end @@ -486,7 +572,7 @@ function solve!(alf::HypatiaOpt) mu = (dot(ts, tz) + tau*kap)/bnu # skip correction phase if allowed and current iterate is in the eta-neighborhood - if alf.corrcheck && (nbhd < abs2(eta*mu)) + if opt.corrcheck && (nbhd < abs2(eta*mu)) continue end @@ -502,12 +588,12 @@ function solve!(alf::HypatiaOpt) calcg!(g, cone) @. tmp_tz = -tz - mu*g @. tmp_ts = 0.0 - (tmp_kap, tmp_tau) = finddirection!(tmp_tx, tmp_ty, tmp_tz, tmp_ts, -kap + mu/tau, 0.0, mu, tau, alf) + (tmp_kap, tmp_tau) = solvelinsys6!(tmp_tx, tmp_ty, tmp_tz, tmp_ts, -kap + mu/tau, 0.0, mu, tau, opt.L) # determine step length alpha by line search - alpha = alf.alphacorr + alpha = opt.alphacorr ncorrlsiters = 0 - while ncorrlsiters <= alf.maxcorrlsiters + while ncorrlsiters <= opt.maxcorrlsiters ncorrlsiters += 1 @. ls_ts = ts + alpha*tmp_ts @@ -517,15 +603,14 @@ function solve!(alf::HypatiaOpt) end # primal iterate tx is outside the cone - if ncorrlsiters == alf.maxcorrlsiters + if ncorrlsiters == opt.maxcorrlsiters # corrector failed corrfail = true - alf.verbose && println("corrector could not improve the solution ($ncorrlsiters line search steps); terminating") - alf.status = :CorrectorFail + opt.verbose && println("corrector could not improve the solution ($ncorrlsiters line search steps); terminating") break end - alpha = alf.corrlsmulti*alpha # decrease alpha + alpha = opt.corrlsmulti*alpha # decrease alpha end if corrfail break @@ -541,33 +626,38 @@ function solve!(alf::HypatiaOpt) mu = (dot(ts, tz) + tau*kap)/bnu # finish if allowed and current iterate is in the eta-neighborhood, or if taken max steps - if (ncorrsteps == alf.maxcorrsteps) || alf.corrcheck + if (ncorrsteps == opt.maxcorrsteps) || opt.corrcheck @. ls_tz = tz nbhd = calcnbhd(tau*kap, mu, ls_tz, g, cone) if nbhd <= abs2(eta*mu) break - elseif ncorrsteps == alf.maxcorrsteps + elseif ncorrsteps == opt.maxcorrsteps # outside eta neighborhood, so corrector failed corrfail = true - alf.verbose && println("corrector phase finished outside the eta-neighborhood ($ncorrsteps correction steps); terminating") - alf.status = :CorrectorFail + opt.verbose && println("corrector phase finished outside the eta-neighborhood ($ncorrsteps correction steps); terminating") break end end end if corrfail + opt.status = :CorrectorFail + invtau = inv(tau) + opt.x = tx .*= invtau + opt.s = ts .*= invtau + opt.y = ty .*= invtau + opt.z = tz .*= invtau break end end - alf.verbose && println("\nfinished in $iter iterations; internal status is $(alf.status)\n") + opt.verbose && println("\nfinished in $iter iterations; internal status is $(opt.status)\n") # calculate solution and iteration statistics - alf.niters = iter - alf.tau = tau - alf.kap = kap - alf.mu = mu - alf.solvetime = time() - starttime + opt.niters = iter + opt.tau = tau + opt.kap = kap + opt.mu = mu + opt.solvetime = time() - starttime return nothing end @@ -580,171 +670,6 @@ function calcnbhd(tk, mu, ls_tz, g, cone) return (tk - mu)^2 + dot(ls_tz, g) end -# calculate initial central primal-dual iterate -function findinitialiterate!(tx::Vector{Float64}, ty::Vector{Float64}, tz::Vector{Float64}, ts::Vector{Float64}, ls_ts::Vector{Float64}, bnu::Float64, alf::HypatiaOpt) - (b, G, h, cone) = (alf.b, alf.G, alf.h, alf.cone) - L = alf.L - (Q2, RiQ1, GQ2, GHGQ2, Q2GHGQ2, bxGHbz, Q1x, rhs, Q2div, Q2x, HGxi) = (L.Q2, L.RiQ1, L.GQ2, L.GHGQ2, L.Q2GHGQ2, L.bxGHbz, L.Q1x, L.rhs, L.Q2div, L.Q2x, L.HGxi) - - alf.verbose && println("\nfinding initial iterate") - - # solve linear equation - # |0 A' G'| * |tx| = |0| - # |A 0 0 | |ty| |b| - # |G 0 -I| |ts| |h| - - # Q1x = Q1*Ri'*b - # Q2x = Q2*(Symmetric(GQ2'*GQ2)\(GQ2'*(h - G*Q1x))) - # tx = Q1x + Q2x - # ts = h - G*tx - # ty = Ri*Q1'*G'*ts - - # bunchkaufman allocates more than cholesky, but doesn't fail when approximately quasidefinite (TODO could try LDL instead) - # TODO does it matter that F could be either type? - mul!(Q2GHGQ2, GQ2', GQ2) - F = cholesky!(Symmetric(Q2GHGQ2), check=false) - if !issuccess(F) - alf.verbose && println("linear system matrix was not positive definite") - mul!(Q2GHGQ2, Q2', GHGQ2) - F = bunchkaufman!(Symmetric(Q2GHGQ2)) - end - - # Q1x = Q1*Ri'*b - mul!(Q1x, RiQ1', b) - # Q2x = Q2*(F\(GQ2'*(h - G*Q1x))) - mul!(HGxi, G, Q1x) - @. HGxi = h - HGxi - mul!(Q2div, GQ2', HGxi) - ldiv!(F, Q2div) - mul!(Q2x, Q2, Q2div) - # tx = Q1x + Q2x - @. tx = Q1x + Q2x - # ts = h - G*tx - mul!(ts, G, tx) - @. ts = h - ts - # ty = Ri*Q1'*G'*ts - mul!(bxGHbz, G', ts) - mul!(ty, RiQ1, bxGHbz) - - # from ts, step along interior direction of cone until ts is inside cone - @. ls_ts = ts - if !incone(cone) - tmp_ts = getintdir!(tz, cone) - alpha = 1.0 # TODO starting alpha maybe should depend on ls_ts (eg norm like in Hypatia) in case 1.0 is too large/small - steps = 0 - @. ls_ts = ts + alpha*tmp_ts - while !incone(cone) - steps += 1 - if steps > 25 - error("cannot find initial iterate") - end - alpha *= 1.5 - @. ls_ts = ts + alpha*tmp_ts - end - alf.verbose && println("$steps steps taken for initial iterate") - @. ts = ls_ts - end - - calcg!(tz, cone) - @. tz *= -1.0 - tau = 1.0 - kap = 1.0 - mu = (dot(tz, ts) + tau*kap)/bnu - @assert abs(1.0 - mu) < 1e-10 - # @assert calcnbhd(tau*kap, mu, copy(tz), copy(tz), cone) < 1e-6 - - alf.verbose && println("initial iterate found") - - return (tau, kap, mu) -end - -# calculate new prediction or correction direction given rhs of KKT linear system -function finddirection!(rhs_tx::Vector{Float64}, rhs_ty::Vector{Float64}, rhs_tz::Vector{Float64}, rhs_ts::Vector{Float64}, rhs_kap::Float64, rhs_tau::Float64, mu::Float64, tau::Float64, alf::HypatiaOpt) - (c, b, G, h, cone) = (alf.c, alf.b, alf.G, alf.h, alf.cone) - L = alf.L - (Q2, HG, GHG, GHGQ2, Q2GHGQ2, x1, y1, z1) = (L.Q2, L.HG, L.GHG, L.GHGQ2, L.Q2GHGQ2, L.x1, L.y1, L.z1) - - # solve two symmetric systems and combine the solutions - # use QR + cholesky method from CVXOPT - # (1) eliminate equality constraints via QR of A' - # (2) solve reduced system by cholesky - # |0 A' G' | * |ux| = |bx| - # |A 0 0 | |uy| |by| - # |G 0 -Hi/mu| |uz| |bz| - - # A' = [Q1 Q2] * [R1; 0] - # factorize Q2' * G' * mu*H * G * Q2 - calcHarr!(HG, G, cone) - @. HG *= mu - mul!(GHG, G', HG) - mul!(GHGQ2, GHG, Q2) - mul!(Q2GHGQ2, Q2', GHGQ2) - - # bunchkaufman allocates more than cholesky, but doesn't fail when approximately quasidefinite (TODO could try LDL instead) - # TODO does it matter that F could be either type? - F = cholesky!(Symmetric(Q2GHGQ2), check=false) - if !issuccess(F) - alf.verbose && println("linear system matrix was not positive definite") - mul!(Q2GHGQ2, Q2', GHGQ2) - F = bunchkaufman!(Symmetric(Q2GHGQ2)) - end - - # (x2, y2, z2) = (rhs_tx, -rhs_ty, -mu*H*rhs_ts - rhs_tz) - @. rhs_ty *= -1.0 - @. rhs_tz *= -1.0 - if !iszero(rhs_ts) - calcHarr!(z1, rhs_ts, cone) - @. rhs_tz -= mu*z1 - end - calcxyz!(rhs_tx, rhs_ty, rhs_tz, F, G, L) - - # (x1, y1, z1) = (-c, b, mu*H*h) - @. x1 = -c - @. y1 = b - calcHarr!(z1, h, cone) - @. z1 *= mu - calcxyz!(x1, y1, z1, F, G, L) - - # combine - dir_tau = (rhs_tau + rhs_kap + dot(c, rhs_tx) + dot(b, rhs_ty) + dot(h, rhs_tz))/(mu/tau/tau - dot(c, x1) - dot(b, y1) - dot(h, z1)) - @. rhs_tx += dir_tau*x1 - @. rhs_ty += dir_tau*y1 - @. rhs_tz += dir_tau*z1 - mul!(z1, G, rhs_tx) - @. rhs_ts = -z1 + h*dir_tau - rhs_ts - dir_kap = -dot(c, rhs_tx) - dot(b, rhs_ty) - dot(h, rhs_tz) - rhs_tau - - return (dir_kap, dir_tau) -end - -# calculate solution to reduced symmetric linear systems -function calcxyz!(xi::Vector{Float64}, yi::Vector{Float64}, zi::Vector{Float64}, F, G::AbstractMatrix{Float64}, L::LinSysCache) - (Q2, RiQ1, HG, GHG, bxGHbz, Q1x, rhs, Q2div, Q2x, GHGxi, HGxi) = (L.Q2, L.RiQ1, L.HG, L.GHG, L.bxGHbz, L.Q1x, L.rhs, L.Q2div, L.Q2x, L.GHGxi, L.HGxi) - - # bxGHbz = bx + G'*Hbz - mul!(bxGHbz, G', zi) - @. bxGHbz += xi - # Q1x = Q1*Ri'*by - mul!(Q1x, RiQ1', yi) - # Q2x = Q2*(K22_F\(Q2'*(bxGHbz - GHG*Q1x))) - mul!(rhs, GHG, Q1x) - @. rhs = bxGHbz - rhs - mul!(Q2div, Q2', rhs) - ldiv!(F, Q2div) - mul!(Q2x, Q2, Q2div) - # xi = Q1x + Q2x - @. xi = Q1x + Q2x - # yi = Ri*Q1'*(bxGHbz - GHG*xi) - mul!(GHGxi, GHG, xi) - @. bxGHbz -= GHGxi - mul!(yi, RiQ1, bxGHbz) - # zi = HG*xi - Hbz - mul!(HGxi, HG, xi) - @. zi = HGxi - zi - - return (xi, yi, zi) -end - # get neighborhood parameters depending on magnitude of barrier parameter and maximum number of correction steps function getbetaeta(maxcorrsteps::Int, bnu::Float64) if maxcorrsteps <= 2 diff --git a/test/examples.jl b/test/examples.jl new file mode 100644 index 000000000..a6b0d4688 --- /dev/null +++ b/test/examples.jl @@ -0,0 +1,210 @@ +#= +Copyright 2018, Chris Coey and contributors +=# + +function _envelope1(verbose::Bool, lscachetype) + # dense methods + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_envelope!(2, 5, 1, 5, use_data=true, dense=true) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ -25.502777 atol=1e-4 rtol=1e-4 + @test r.niters <= 35 + + # sparse methods + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_envelope!(2, 5, 1, 5, use_data=true, dense=false) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ -25.502777 atol=1e-4 rtol=1e-4 + @test r.niters <= 35 +end + +function _envelope2(verbose::Bool, lscachetype) + # dense methods + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_envelope!(2, 4, 2, 7, dense=true) + rd = fullsolve(opt, c, A, b, G, h, cone) + @test rd.status == :Optimal + @test rd.niters <= 60 + @test rd.pobj ≈ rd.dobj atol=1e-4 rtol=1e-4 + + # sparse methods + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_envelope!(2, 4, 2, 7, dense=false) + rs = fullsolve(opt, c, A, b, G, h, cone) + @test rs.status == :Optimal + @test rs.niters <= 60 + @test rs.pobj ≈ rs.dobj atol=1e-4 rtol=1e-4 + + @test rs.pobj ≈ rd.pobj atol=1e-4 rtol=1e-4 +end + +function _envelope3(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_envelope!(2, 3, 3, 5, dense=false) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 +end + +function _envelope4(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose, tolrelopt=1e-5, tolabsopt=1e-6, tolfeas=1e-6) + (c, A, b, G, h, cone) = build_envelope!(2, 3, 4, 4, dense=false) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 +end + +function _lp1(verbose::Bool, lscachetype) + # dense methods + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_lp!(50, 100, dense=true, tosparse=false) + rd = fullsolve(opt, c, A, b, G, h, cone) + @test rd.status == :Optimal + @test rd.niters <= 45 + @test rd.pobj ≈ rd.dobj atol=1e-4 rtol=1e-4 + + # sparse methods + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_lp!(50, 100, dense=true, tosparse=true) + rs = fullsolve(opt, c, A, b, G, h, cone) + @test rs.status == :Optimal + @test rs.niters <= 45 + @test rs.pobj ≈ rs.dobj atol=1e-4 rtol=1e-4 + + @test rs.pobj ≈ rd.pobj atol=1e-4 rtol=1e-4 +end + +function _lp2(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_lp!(500, 1000, use_data=true, dense=true) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 75 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ 2055.807 atol=1e-4 rtol=1e-4 +end + +function _lp3(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_lp!(500, 1000, dense=false, nzfrac=10/1000) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 70 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 +end + +# for namedpoly tests, most optimal values are taken from https://people.sc.fsu.edu/~jburkardt/py_src/polynomials/polynomials.html + +function _namedpoly1(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_namedpoly!(:butcher, 2) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 40 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ -1.4393333333 atol=1e-4 rtol=1e-4 +end + +function _namedpoly2(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose, tolfeas=5e-7) + (c, A, b, G, h, cone) = build_namedpoly!(:caprasse, 4) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 45 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ -3.1800966258 atol=1e-4 rtol=1e-4 +end + +function _namedpoly3(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose, tolfeas=1e-10) + (c, A, b, G, h, cone) = build_namedpoly!(:goldsteinprice, 7) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 60 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ 3 atol=1e-4 rtol=1e-4 +end + +function _namedpoly4(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_namedpoly!(:heart, 2) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + # @test r.niters <= 40 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ -1.36775 atol=1e-4 rtol=1e-4 +end + +function _namedpoly5(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_namedpoly!(:lotkavolterra, 3) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 35 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ -20.8 atol=1e-4 rtol=1e-4 +end + +function _namedpoly6(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_namedpoly!(:magnetism7, 2) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + # @test r.niters <= 40 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ -0.25 atol=1e-4 rtol=1e-4 +end + +function _namedpoly7(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose, tolrelopt=1e-5, tolabsopt=1e-6, tolfeas=1e-6) + (c, A, b, G, h, cone) = build_namedpoly!(:motzkin, 7) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 40 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ 0 atol=1e-4 rtol=1e-4 +end + +function _namedpoly8(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_namedpoly!(:reactiondiffusion, 4) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 35 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ -36.71269068 atol=1e-4 rtol=1e-4 +end + +function _namedpoly9(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_namedpoly!(:robinson, 8) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 40 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ 0.814814 atol=1e-4 rtol=1e-4 +end + +function _namedpoly10(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose, tolfeas=1.1e-8) + (c, A, b, G, h, cone) = build_namedpoly!(:rosenbrock, 3) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 65 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ 0 atol=1e-4 rtol=1e-4 +end + +function _namedpoly11(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) + (c, A, b, G, h, cone) = build_namedpoly!(:schwefel, 4) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 50 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ 0 atol=1e-4 rtol=1e-4 +end diff --git a/test/moi.jl b/test/moi.jl new file mode 100644 index 000000000..a71c169e8 --- /dev/null +++ b/test/moi.jl @@ -0,0 +1,66 @@ +#= +Copyright 2018, Chris Coey and contributors +=# + +using MathOptInterface +MOI = MathOptInterface +MOIT = MOI.Test +MOIB = MOI.Bridges +MOIU = MOI.Utilities + +MOIU.@model(HypatiaModelData, + (), + ( + MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, + # MOI.Interval, + ), + ( + MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, + MOI.SecondOrderCone, MOI.RotatedSecondOrderCone, + MOI.PositiveSemidefiniteConeTriangle, + MOI.ExponentialCone, + # MOI.PowerCone, + ), + (), + (MOI.SingleVariable,), + (MOI.ScalarAffineFunction,), + (MOI.VectorOfVariables,), + (MOI.VectorAffineFunction,), + ) + +function testmoi(verbose::Bool, usedense::Bool) + optimizer = MOIU.CachingOptimizer(HypatiaModelData{Float64}(), Hypatia.HypatiaOptimizer(verbose=verbose, usedense=usedense)) + + config = MOIT.TestConfig( + atol = 1e-3, + rtol = 1e-3, + solve = true, + query = true, + modify_lhs = true, + duals = true, + infeas_certificates = true, + ) + + @testset "MathOptInterface tests" begin + @testset "Continuous linear problems" begin + MOIT.contlineartest( + MOIB.SplitInterval{Float64}( + optimizer + ), + config) + end + @testset "Continuous conic problems" begin + exclude = ["rootdet", "logdet", "sdp"] # TODO MOI does not yet support scaled PSD triangle + MOIT.contconictest( + MOIB.GeoMean{Float64}( + # MOIB.SquarePSD{Float64}( + # MOIB.LogDet{Float64}( + # MOIB.RootDet{Float64}( + optimizer + ),#))), + config, exclude) + end + end + + return nothing +end diff --git a/test/native.jl b/test/native.jl index 10dcf7e5a..1bc051dc2 100644 --- a/test/native.jl +++ b/test/native.jl @@ -2,34 +2,93 @@ Copyright 2018, Chris Coey and contributors =# -@testset "small lp 1: nonnegative vs nonpositive orthant" begin +function _consistent1(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) Random.seed!(1) - (n, p, q) = (10, 8, 10) + (n, p, q) = (30, 15, 30) + A = rand(-9.0:9.0, p, n) + G = Matrix(1.0I, q, n) + c = rand(0.0:9.0, n) + rnd1 = rand() + rnd2 = rand() + A[11:15,:] = rnd1*A[1:5,:] - rnd2*A[6:10,:] + b = A*ones(n) + rnd1 = rand() + rnd2 = rand() + A[:,11:15] = rnd1*A[:,1:5] - rnd2*A[:,6:10] + G[:,11:15] = rnd1*G[:,1:5] - rnd2*G[:,6:10] + c[11:15] = rnd1*c[1:5] - rnd2*c[6:10] + h = zeros(q) + cone = Hypatia.Cone([Hypatia.NonpositiveCone(q)], [1:q]) + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 +end + +function _inconsistent1(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) + Random.seed!(1) + (n, p, q) = (30, 15, 30) + A = rand(-9.0:9.0, p, n) + G = Matrix(-1.0I, q, n) + c = rand(0.0:9.0, n) + b = rand(p) + rnd1 = rand() + rnd2 = rand() + A[11:15,:] = rnd1*A[1:5,:] - rnd2*A[6:10,:] + b[11:15] = 2*(rnd1*b[1:5] - rnd2*b[6:10]) + h = zeros(q) + cone = Hypatia.Cone([Hypatia.NonnegativeCone(q)], [1:q]) + @test_throws ErrorException("some primal equality constraints are inconsistent") fullsolve(opt, c, A, b, G, h, cone) +end + +function _inconsistent2(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) + Random.seed!(1) + (n, p, q) = (30, 15, 30) + A = rand(-9.0:9.0, p, n) + G = Matrix(-1.0I, q, n) + c = rand(0.0:9.0, n) + b = rand(p) + rnd1 = rand() + rnd2 = rand() + A[:,11:15] = rnd1*A[:,1:5] - rnd2*A[:,6:10] + G[:,11:15] = rnd1*G[:,1:5] - rnd2*G[:,6:10] + c[11:15] = 2*(rnd1*c[1:5] - rnd2*c[6:10]) + h = zeros(q) + cone = Hypatia.Cone([Hypatia.NonnegativeCone(q)], [1:q]) + @test_throws ErrorException("some dual equality constraints are inconsistent") fullsolve(opt, c, A, b, G, h, cone) +end + +function _orthant1(verbose::Bool, lscachetype) + Random.seed!(1) + (n, p, q) = (40, 20, 40) c = rand(0.0:9.0, n) A = rand(-9.0:9.0, p, n) b = A*ones(n) h = zeros(q) - alf1 = Hypatia.HypatiaOpt(verbose=verbflag) + # nonnegative cone + opt = Hypatia.Optimizer(verbose=verbose) G = SparseMatrixCSC(-1.0I, q, n) cone = Hypatia.Cone([Hypatia.NonnegativeCone(q)], [1:q]) - Hypatia.load_data!(alf1, c, A, b, G, h, cone) - @time Hypatia.solve!(alf1) - @test Hypatia.get_status(alf1) == :Optimal + rnn = fullsolve(opt, c, A, b, G, h, cone) + @test rnn.status == :Optimal + @test rnn.pobj ≈ rnn.dobj atol=1e-4 rtol=1e-4 - alf2 = Hypatia.HypatiaOpt(verbose=verbflag) + # nonpositive cone + opt = Hypatia.Optimizer(verbose=verbose) G = SparseMatrixCSC(1.0I, q, n) cone = Hypatia.Cone([Hypatia.NonpositiveCone(q)], [1:q]) - Hypatia.load_data!(alf2, c, A, b, G, h, cone) - @time Hypatia.solve!(alf2) - @test Hypatia.get_status(alf2) == :Optimal + rnp = fullsolve(opt, c, A, b, G, h, cone) + @test rnp.status == :Optimal + @test rnp.pobj ≈ rnp.dobj atol=1e-4 rtol=1e-4 - @test Hypatia.get_pobj(alf1) ≈ Hypatia.get_pobj(alf2) atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf1) ≈ Hypatia.get_dobj(alf2) atol=1e-4 rtol=1e-4 + @test rnp.pobj ≈ rnn.pobj atol=1e-4 rtol=1e-4 end -@testset "small lp 2" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) +function _orthant2(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) Random.seed!(1) (n, p, q) = (5, 2, 10) c = rand(0.0:9.0, n) @@ -38,47 +97,45 @@ end G = rand(q, n) - Matrix(2.0I, q, n) h = G*ones(n) cone = Hypatia.Cone([Hypatia.NonnegativeCone(q)], [1:q]) - Hypatia.load_data!(alf, c, A, b, G, h, cone) - @time Hypatia.solve!(alf) - @test Hypatia.get_status(alf) == :Optimal + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 end -@testset "small lp 3" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) +function _orthant3(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) Random.seed!(1) (n, p, q) = (30, 12, 30) c = rand(0.0:9.0, n) A = rand(-9.0:9.0, p, n) b = A*ones(n) - @assert n == q G = Diagonal(1.0I, n) h = zeros(q) cone = Hypatia.Cone([Hypatia.NonpositiveCone(q)], [1:q]) - Hypatia.load_data!(alf, c, A, b, G, h, cone) - @time Hypatia.solve!(alf) - @test Hypatia.get_status(alf) == :Optimal + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 end -@testset "small L_infinity cone problem" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) +function _ellinf1(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) c = Float64[0, -1, -1] A = Float64[1 0 0; 0 1 0] b = Float64[1, 1/sqrt(2)] G = SparseMatrixCSC(-1.0I, 3, 3) h = zeros(3) cone = Hypatia.Cone([Hypatia.EllInfinityCone(3)], [1:3]) - Hypatia.load_data!(alf, c, A, b, G, h, cone) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 20 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ -1 - 1/sqrt(2) atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ Hypatia.get_pobj(alf) atol=1e-4 rtol=1e-4 - @test Hypatia.get_y(alf) ≈ [1, 1] atol=1e-4 rtol=1e-4 - @test Hypatia.get_x(alf) ≈ [1, 1/sqrt(2), 1] atol=1e-4 rtol=1e-4 + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 20 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ -1 - 1/sqrt(2) atol=1e-4 rtol=1e-4 + @test r.x ≈ [1, 1/sqrt(2), 1] atol=1e-4 rtol=1e-4 + @test r.y ≈ [1, 1] atol=1e-4 rtol=1e-4 end -@testset "small L_infinity cone problem 2" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) +function _ellinf2(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) Random.seed!(1) c = Float64[1, 0, 0, 0, 0, 0] A = rand(-9.0:9.0, 3, 6) @@ -86,337 +143,124 @@ end G = rand(6, 6) h = G*ones(6) cone = Hypatia.Cone([Hypatia.EllInfinityCone(6)], [1:6]) - Hypatia.load_data!(alf, c, A, b, G, h, cone) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 20 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ 1 atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ Hypatia.get_pobj(alf) atol=1e-4 rtol=1e-4 + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 25 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ 1 atol=1e-4 rtol=1e-4 end -@testset "small second-order cone problem" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) +function _soc1(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) c = Float64[0, -1, -1] A = Float64[1 0 0; 0 1 0] b = Float64[1, 1/sqrt(2)] G = SparseMatrixCSC(-1.0I, 3, 3) h = zeros(3) cone = Hypatia.Cone([Hypatia.SecondOrderCone(3)], [1:3]) - Hypatia.load_data!(alf, c, A, b, G, h, cone) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 15 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ -sqrt(2) atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ Hypatia.get_pobj(alf) atol=1e-4 rtol=1e-4 - @test Hypatia.get_y(alf) ≈ [sqrt(2), 0] atol=1e-4 rtol=1e-4 - @test Hypatia.get_x(alf) ≈ [1, 1/sqrt(2), 1/sqrt(2)] atol=1e-4 rtol=1e-4 + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 20 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ -sqrt(2) atol=1e-4 rtol=1e-4 + @test r.x ≈ [1, 1/sqrt(2), 1/sqrt(2)] atol=1e-4 rtol=1e-4 + @test r.y ≈ [sqrt(2), 0] atol=1e-4 rtol=1e-4 end -@testset "small rotated second-order cone problem" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) +function _rsoc1(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) c = Float64[0, 0, -1, -1] A = Float64[1 0 0 0; 0 1 0 0] b = Float64[1/2, 1] G = SparseMatrixCSC(-1.0I, 4, 4) h = zeros(4) cone = Hypatia.Cone([Hypatia.RotatedSecondOrderCone(4)], [1:4]) - Hypatia.load_data!(alf, c, A, b, G, h, cone) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 15 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ -sqrt(2) atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ Hypatia.get_pobj(alf) atol=1e-4 rtol=1e-4 - @test Hypatia.get_x(alf)[3:4] ≈ [1, 1]/sqrt(2) atol=1e-4 rtol=1e-4 + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 20 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ -sqrt(2) atol=1e-4 rtol=1e-4 + @test r.x[3:4] ≈ [1, 1]/sqrt(2) atol=1e-4 rtol=1e-4 end -@testset "small rotated second-order cone problem 2" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) +function _rsoc2(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) c = Float64[0, 0, -1] A = Float64[1 0 0; 0 1 0] b = Float64[1/2, 1]/sqrt(2) G = SparseMatrixCSC(-1.0I, 3, 3) h = zeros(3) cone = Hypatia.Cone([Hypatia.RotatedSecondOrderCone(3)], [1:3]) - Hypatia.load_data!(alf, c, A, b, G, h, cone) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 20 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ -1/sqrt(2) atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ Hypatia.get_pobj(alf) atol=1e-4 rtol=1e-4 - @test Hypatia.get_x(alf)[2] ≈ 1/sqrt(2) atol=1e-4 rtol=1e-4 + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 20 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ -1/sqrt(2) atol=1e-4 rtol=1e-4 + @test r.x[2] ≈ 1/sqrt(2) atol=1e-4 rtol=1e-4 end -@testset "small positive semidefinite cone problem" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) +function _psd1(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) c = Float64[0, -1, 0] A = Float64[1 0 0; 0 0 1] b = Float64[1/2, 1] G = SparseMatrixCSC(-1.0I, 3, 3) h = zeros(3) cone = Hypatia.Cone([Hypatia.PositiveSemidefiniteCone(3)], [1:3]) - Hypatia.load_data!(alf, c, A, b, G, h, cone) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 15 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ -1 atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ Hypatia.get_pobj(alf) atol=1e-4 rtol=1e-4 - @test Hypatia.get_x(alf)[2] ≈ 1 atol=1e-4 rtol=1e-4 + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 20 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ -1 atol=1e-4 rtol=1e-4 + @test r.x[2] ≈ 1 atol=1e-4 rtol=1e-4 end -@testset "small positive semidefinite cone problem 2" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) +function _psd2(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) c = Float64[1, 0, 1, 0, 0, 1] A = Float64[1 2 3 4 5 6; 1 1 1 1 1 1] b = Float64[10, 3] G = SparseMatrixCSC(-1.0I, 6, 6) h = zeros(6) cone = Hypatia.Cone([Hypatia.PositiveSemidefiniteCone(6)], [1:6]) - Hypatia.load_data!(alf, c, A, b, G, h, cone) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 20 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ 1.249632 atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ Hypatia.get_pobj(alf) atol=1e-4 rtol=1e-4 - @test Hypatia.get_x(alf) ≈ [0.491545, 0.647333, 0.426249, 0.571161, 0.531874, 0.331838] atol=1e-4 rtol=1e-4 + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 20 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ 1.249632 atol=1e-4 rtol=1e-4 + @test r.x ≈ [0.491545, 0.647333, 0.426249, 0.571161, 0.531874, 0.331838] atol=1e-4 rtol=1e-4 end -@testset "small exponential cone problem" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) +function _exp1(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) c = Float64[1, 1, 1] A = Float64[0 1 0; 1 0 0] b = Float64[2, 1] G = SparseMatrixCSC(-1.0I, 3, 3) h = zeros(3) cone = Hypatia.Cone([Hypatia.ExponentialCone()], [1:3]) - Hypatia.load_data!(alf, c, A, b, G, h, cone) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 20 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ (2*exp(1/2)+3) atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ Hypatia.get_pobj(alf) atol=1e-4 rtol=1e-4 - @test dot(Hypatia.get_y(alf), b) ≈ -Hypatia.get_pobj(alf) atol=1e-4 rtol=1e-4 - @test Hypatia.get_x(alf) ≈ [1, 2, 2*exp(1/2)] atol=1e-4 rtol=1e-4 - @test Hypatia.get_y(alf) ≈ -[1+exp(1/2)/2, 1+exp(1/2)] atol=1e-4 rtol=1e-4 - @test Hypatia.get_z(alf) ≈ (c + A'*Hypatia.get_y(alf)) atol=1e-4 rtol=1e-4 + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 20 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ 2*exp(1/2)+3 atol=1e-4 rtol=1e-4 + @test r.x ≈ [1, 2, 2*exp(1/2)] atol=1e-4 rtol=1e-4 + @test r.y ≈ -[1+exp(1/2)/2, 1+exp(1/2)] atol=1e-4 rtol=1e-4 + @test r.z ≈ c+A'*r.y atol=1e-4 rtol=1e-4 end -@testset "small power cone problem" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) +function _power1(verbose::Bool, lscachetype) + opt = Hypatia.Optimizer(verbose=verbose) c = Float64[1, 0, 0, -1, -1, 0] A = Float64[1 1 1/2 0 0 0; 0 0 0 0 0 1] b = Float64[2, 1] G = SparseMatrixCSC(-1.0I, 6, 6)[[4, 1, 2, 5, 3, 6], :] h = zeros(6) cone = Hypatia.Cone([Hypatia.PowerCone([0.2, 0.8]), Hypatia.PowerCone([0.4, 0.6])], [1:3, 4:6]) - Hypatia.load_data!(alf, c, A, b, G, h, cone) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 25 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ -1.80734 atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ Hypatia.get_pobj(alf) atol=1e-4 rtol=1e-4 - @test Hypatia.get_x(alf)[1:3] ≈ [0.0639314, 0.783361, 2.30542] atol=1e-4 rtol=1e-4 -end - -@testset "small dense lp example (dense vs sparse A)" begin - # dense methods - alf2 = Hypatia.HypatiaOpt(verbose=verbflag) - build_lp!(alf2, 50, 100, dense=true, tosparse=false) - @time Hypatia.solve!(alf2) - @test Hypatia.get_niters(alf2) <= 40 - @test Hypatia.get_status(alf2) == :Optimal - - # sparse methods - alf1 = Hypatia.HypatiaOpt(verbose=verbflag) - build_lp!(alf1, 50, 100, dense=true, tosparse=true) - @time Hypatia.solve!(alf1) - @test Hypatia.get_niters(alf1) <= 40 - @test Hypatia.get_status(alf1) == :Optimal - - @test Hypatia.get_pobj(alf2) ≈ Hypatia.get_pobj(alf1) atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf2) ≈ Hypatia.get_dobj(alf1) atol=1e-4 rtol=1e-4 -end - -@testset "1D poly envelope example (dense vs sparse A)" begin - # dense methods - alf2 = Hypatia.HypatiaOpt(verbose=verbflag) - build_envelope!(alf2, 2, 5, 1, 5, use_data=true, dense=true) - @time Hypatia.solve!(alf2) - @test Hypatia.get_niters(alf2) <= 30 - @test Hypatia.get_status(alf2) == :Optimal - @test Hypatia.get_pobj(alf2) ≈ -25.502777 atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf2) ≈ -25.502777 atol=1e-4 rtol=1e-4 - - # sparse methods - alf1 = Hypatia.HypatiaOpt(verbose=verbflag) - build_envelope!(alf1, 2, 5, 1, 5, use_data=true, dense=false) - @time Hypatia.solve!(alf1) - @test Hypatia.get_niters(alf1) <= 30 - @test Hypatia.get_status(alf1) == :Optimal - @test Hypatia.get_pobj(alf1) ≈ -25.502777 atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf1) ≈ -25.502777 atol=1e-4 rtol=1e-4 -end - -# most values taken from https://people.sc.fsu.edu/~jburkardt/py_src/polynomials/polynomials.html -@testset "Butcher" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) - build_namedpoly!(alf, :butcher, 2) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 40 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ -1.4393333333 atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ -1.4393333333 atol=1e-4 rtol=1e-4 -end - -@testset "Caprasse" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag, tolfeas=5e-7) - build_namedpoly!(alf, :caprasse, 4) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 45 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ -3.1800966258 atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ -3.1800966258 atol=1e-4 rtol=1e-4 -end - -@testset "Goldstein-Price" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag, tolfeas=1e-10) - build_namedpoly!(alf, :goldsteinprice, 7) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 60 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ 3 atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ 3 atol=1e-4 rtol=1e-4 -end - -# out of memory during interpolation calculations -# @testset "Heart" begin -# alf = Hypatia.HypatiaOpt(verbose=verbflag) -# build_namedpoly!(alf, :heart, 2) -# @time Hypatia.solve!(alf) -# @test Hypatia.get_status(alf) == :Optimal -# @test Hypatia.get_pobj(alf) ≈ -1.36775 atol=1e-4 rtol=1e-4 -# @test Hypatia.get_dobj(alf) ≈ -1.36775 atol=1e-4 rtol=1e-4 -# end - -@testset "Lotka-Volterra" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) - build_namedpoly!(alf, :lotkavolterra, 3) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 35 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ -20.8 atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ -20.8 atol=1e-4 rtol=1e-4 -end - -# out of memory during interpolation calculations -# @testset "Magnetism-7" begin -# alf = Hypatia.HypatiaOpt(verbose=verbflag) -# build_namedpoly!(alf, :magnetism7, 2) -# @time Hypatia.solve!(alf) -# @test Hypatia.get_status(alf) == :Optimal -# @test Hypatia.get_pobj(alf) ≈ -0.25 atol=1e-4 rtol=1e-4 -# @test Hypatia.get_dobj(alf) ≈ -0.25 atol=1e-4 rtol=1e-4 -# end - -@testset "Motzkin" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag, tolrelopt=1e-5, tolabsopt=1e-6, tolfeas=1e-6) - build_namedpoly!(alf, :motzkin, 7) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 35 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ 0 atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ 0 atol=1e-4 rtol=1e-4 -end - -@testset "Reaction-diffusion" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) - build_namedpoly!(alf, :reactiondiffusion, 4) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 35 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ -36.71269068 atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ -36.71269068 atol=1e-4 rtol=1e-4 -end - -@testset "Robinson" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) - build_namedpoly!(alf, :robinson, 8) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 40 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ 0.814814 atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ 0.814814 atol=1e-4 rtol=1e-4 -end - -@testset "Rosenbrock" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag, tolfeas=1.1e-8) - build_namedpoly!(alf, :rosenbrock, 3) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 65 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ 0 atol=1e-2 rtol=1e-2 - @test Hypatia.get_dobj(alf) ≈ 0 atol=1e-3 rtol=1e-3 -end - -@testset "Schwefel" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) - build_namedpoly!(alf, :schwefel, 4) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 50 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ 0 atol=1e-3 rtol=1e-3 - @test Hypatia.get_dobj(alf) ≈ 0 atol=1e-3 rtol=1e-3 -end - -@testset "large dense lp example (dense A)" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) - build_lp!(alf, 500, 1000, use_data=true, dense=true) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 75 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ 2055.807 atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf) ≈ 2055.807 atol=1e-4 rtol=1e-4 -end - -@testset "large sparse lp example (sparse A)" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) - build_lp!(alf, 500, 1000, dense=false, nzfrac=10/1000) - @time Hypatia.solve!(alf) - @test Hypatia.get_niters(alf) <= 70 - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ Hypatia.get_dobj(alf) atol=1e-4 rtol=1e-4 -end - -@testset "2D poly envelope example (dense vs sparse A)" begin - # dense methods - alf2 = Hypatia.HypatiaOpt(verbose=verbflag) - build_envelope!(alf2, 2, 4, 2, 7, dense=true) - @time Hypatia.solve!(alf2) - @test Hypatia.get_niters(alf2) <= 55 - @test Hypatia.get_status(alf2) == :Optimal - - # sparse methods - alf1 = Hypatia.HypatiaOpt(verbose=verbflag) - build_envelope!(alf1, 2, 4, 2, 7, dense=false) - @time Hypatia.solve!(alf1) - @test Hypatia.get_niters(alf1) <= 55 - @test Hypatia.get_status(alf1) == :Optimal - - @test Hypatia.get_pobj(alf2) ≈ Hypatia.get_pobj(alf1) atol=1e-4 rtol=1e-4 - @test Hypatia.get_dobj(alf2) ≈ Hypatia.get_dobj(alf1) atol=1e-4 rtol=1e-4 -end - -@testset "3D poly envelope example (sparse A)" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag) - build_envelope!(alf, 2, 3, 3, 5, dense=false) - @time Hypatia.solve!(alf) - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ Hypatia.get_dobj(alf) atol=1e-4 rtol=1e-4 -end - -@testset "4D poly envelope example (sparse A)" begin - alf = Hypatia.HypatiaOpt(verbose=verbflag, tolrelopt=1e-5, tolabsopt=1e-6, tolfeas=1e-6) - build_envelope!(alf, 2, 3, 4, 4, dense=false) - @time Hypatia.solve!(alf) - @test Hypatia.get_status(alf) == :Optimal - @test Hypatia.get_pobj(alf) ≈ Hypatia.get_dobj(alf) atol=1e-4 rtol=1e-4 + r = fullsolve(opt, c, A, b, G, h, cone) + @test r.status == :Optimal + @test r.niters <= 25 + @test r.pobj ≈ r.dobj atol=1e-4 rtol=1e-4 + @test r.pobj ≈ -1.80734 atol=1e-4 rtol=1e-4 + @test r.x[1:3] ≈ [0.0639314, 0.783361, 2.30542] atol=1e-4 rtol=1e-4 end diff --git a/test/runtests.jl b/test/runtests.jl index 6eac55054..352ad81e0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,83 +4,116 @@ Copyright 2018, Chris Coey and contributors using Hypatia using Test +using Random +using LinearAlgebra +using SparseArrays -# native interface tests +# TODO make it a native interface function eventually +# TODO maybe build a new high-level optimizer struct. the current optimizer struct is low-level +function fullsolve(opt::Hypatia.Optimizer, c, A, b, G, h, cone) # TODO handle lscachetype + Hypatia.check_data(c, A, b, G, h, cone) + (c1, A1, b1, G1, prkeep, dukeep, Q2, RiQ1) = Hypatia.preprocess_data(c, A, b, G, useQR=true) + + L = Hypatia.QRSymmCache(c1, A1, b1, G1, h, cone, Q2, RiQ1) + # L = Hypatia.NaiveCache(c1, A1, b1, G1, h, cone) + + Hypatia.load_data!(opt, c1, A1, b1, G1, h, cone, L) + + Hypatia.solve!(opt) + + x = zeros(length(c)) + x[dukeep] = Hypatia.get_x(opt) + y = zeros(length(b)) + y[prkeep] = Hypatia.get_y(opt) + s = Hypatia.get_s(opt) + z = Hypatia.get_z(opt) + + pobj = Hypatia.get_pobj(opt) + dobj = Hypatia.get_dobj(opt) -verbflag = false # Hypatia verbose option + status = Hypatia.get_status(opt) + stime = Hypatia.get_solvetime(opt) + niters = Hypatia.get_niters(opt) + + return (x=x, y=y, s=s, z=z, pobj=pobj, dobj=dobj, status=status, stime=stime, niters=niters) +end -# TODO interpolation tests -# load optimizer builder functions from examples folder +# native interface tests +include(joinpath(@__DIR__, "native.jl")) +verbose = false # test verbosity +lscachetype = Hypatia.QRSymmCache # linear system cache type + +@testset "native interface tests" begin + for testfun in ( + _consistent1, + _inconsistent1, + _inconsistent2, + _orthant1, + _orthant2, + _orthant3, + _ellinf1, + _ellinf2, + _soc1, + _rsoc1, + _rsoc2, + _psd1, + _psd2, + _exp1, + _power1, + ) + testfun(verbose, lscachetype) + end +end + +# examples in src/examples/ folder egs_dir = joinpath(@__DIR__, "../examples") include(joinpath(egs_dir, "envelope/envelope.jl")) include(joinpath(egs_dir, "lp/lp.jl")) include(joinpath(egs_dir, "namedpoly/namedpoly.jl")) -# run native and MOI interfaces on examples -include(joinpath(@__DIR__, "native.jl")) - - -# MathOptInterface tests - -import MathOptInterface -MOI = MathOptInterface -MOIT = MOI.Test -MOIB = MOI.Bridges -MOIU = MOI.Utilities - -MOIU.@model(HypatiaModelData, - (), - ( - MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, - # MOI.Interval, - ), - ( - MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, - MOI.SecondOrderCone, MOI.RotatedSecondOrderCone, - MOI.PositiveSemidefiniteConeTriangle, - MOI.ExponentialCone, - # MOI.PowerCone, - ), - (), - (MOI.SingleVariable,), - (MOI.ScalarAffineFunction,), - (MOI.VectorOfVariables,), - (MOI.VectorAffineFunction,), - ) - -optimizer = MOIU.CachingOptimizer(HypatiaModelData{Float64}(), Hypatia.Optimizer()) - -config = MOIT.TestConfig( - atol=1e-4, - rtol=1e-4, - solve=true, - query=true, - modify_lhs=true, - duals=true, - infeas_certificates=true, - ) - -@testset "Continuous linear problems" begin - MOIT.contlineartest( - MOIB.SplitInterval{Float64}( - optimizer - ), - config) +@testset "default examples" begin + run_envelope() + run_lp() + run_namedpoly() end -@testset "Continuous conic problems" begin - exclude = ["rootdet", "logdet", "sdp"] # TODO bridges not working? should not need to exclude in future - MOIT.contconictest( - # MOIB.SquarePSD{Float64}( - MOIB.GeoMean{Float64}( - # MOIB.LogDet{Float64}( - # MOIB.RootDet{Float64}( - optimizer - ),#))), - config, exclude) +include(joinpath(@__DIR__, "examples.jl")) +verbose = false # test verbosity +lscachetype = Hypatia.QRSymmCache # linear system cache type + +@testset "varied examples" begin + for testfun in ( + _envelope1, + _envelope2, + _envelope3, + # _envelope4, + _lp1, + # _lp2, + # _lp3, + _namedpoly1, + _namedpoly2, + _namedpoly3, + # _namedpoly4, + _namedpoly5, + # _namedpoly6, + _namedpoly7, + _namedpoly8, + _namedpoly9, + # _namedpoly10, + # _namedpoly11, + ) + testfun(verbose, lscachetype) + end end +# MathOptInterface tests +verbose = false # test verbosity +include(joinpath(@__DIR__, "moi.jl")) +testmoi(verbose, false) +testmoi(verbose, true) + + return nothing