Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

WIP: preprocessing using QR factorizations #73

Closed
wants to merge 1 commit into from
Closed

WIP: preprocessing using QR factorizations #73

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
76 changes: 64 additions & 12 deletions src/nativeinterface.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -177,7 +177,7 @@ function load_data!(
G::AbstractMatrix{Float64},
h::Vector{Float64},
cone::Cone;
check::Bool=false, # check rank conditions
check::Bool=true, # check rank conditions
)

# check data consistency
Expand All @@ -201,10 +201,73 @@ function load_data!(
@assert dimension(cone.prms[k]) == length(cone.idxs[k])
end


if issparse(A)
error("not implemented for sparse A")
end

prseps = 1e-10 # presolve epsilon


# preprocess primal equality constraints
# F = qr(A')
M = A
v = b

samp = rand(size(M, 2), 3)
sol = M'\samp
keeprows = [sum(abs(sol[i,j]) for j in 1:3) > prseps for i in 1:size(M, 1)]

if norm(M*(M[keeprows,:]\v[keeprows,:]) - v, Inf) > prseps
error("some primal equality constraints are inconsistent")
end

A = A[keeprows,:]
b = b[keeprows]
p = length(b)

# preprocess dual equality constraints
M = hcat(A', G')
v = c

samp = rand(size(M, 2), 3)
F = qr(M')
@show M'
@show F
sol = F\samp
keeprows = [sum(abs(sol[i,j]) for j in 1:3) > prseps for i in 1:size(M, 1)]

if norm(M*(M[keeprows,:]\v[keeprows,:]) - v, Inf) > prseps
error("some dual equality constraints are inconsistent")
end

A = A[:,keeprows]
G = G[:,keeprows]
c = c[keeprows]


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 n < p
# error("number of equality constraints ($p) exceeds number of variables ($n)")
# end
# 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")
# 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")
end
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)
# TODO don't reorder to match original: store indices and permute optimal solution
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'))
Expand All @@ -230,17 +293,6 @@ function load_data!(
@assert norm(A'*Ri - Q1) < 1e-8 # TODO delete later
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")
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")
end
end

alf.c = c
alf.A = A
alf.b = b
Expand Down
159 changes: 127 additions & 32 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,19 +8,20 @@ using Test

# native interface tests

verbflag = false # Hypatia verbose option

# TODO interpolation tests

# load optimizer builder functions from 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"))

# verbflag = false # Hypatia verbose option
#
# # TODO interpolation tests
#
# # load optimizer builder functions from 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
# @testset "native interface tests" begin
# include(joinpath(@__DIR__, "native.jl"))
# end

# MathOptInterface tests

Expand Down Expand Up @@ -53,34 +54,128 @@ MOIU.@model(HypatiaModelData,
optimizer = MOIU.CachingOptimizer(HypatiaModelData{Float64}(), Hypatia.Optimizer())

config = MOIT.TestConfig(
atol=1e-4,
rtol=1e-4,
atol=1e-3,
rtol=1e-3,
solve=true,
query=true,
modify_lhs=true,
duals=true,
infeas_certificates=true,
)

@testset "Continuous linear problems" begin
MOIT.contlineartest(
MOIB.SplitInterval{Float64}(
optimizer
),
config)
#
# function linear15test(model::MOI.ModelLike, config::MOIT.TestConfig)
# atol = config.atol
# rtol = config.rtol
# # minimize 0
# # s.t. 0 == 0
# # x == 1
# @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
# @test MOI.supports(model, MOI.ObjectiveSense())
# @test MOI.supports_constraint(model, MOI.VectorAffineFunction{Float64}, MOI.Zeros)
#
# MOI.empty!(model)
# @test MOI.is_empty(model)
#
# x = MOI.add_variables(model, 1)
# # Create a VectorAffineFunction with two rows, but only
# # one term, belonging to the second row. The first row,
# # which is empty, is essentially a constraint that 0 == 0.
# c = MOI.add_constraint(model,
# MOI.VectorAffineFunction(
# MOI.VectorAffineTerm.(2, MOI.ScalarAffineTerm.([1.0], x)),
# zeros(2)
# ),
# MOI.Zeros(2)
# )
#
# MOI.set(model,
# MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
# MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([0.0], x), 0.0))
# MOI.set(model, MOI.ObjectiveSense(), MOI.MinSense)
#
# if config.solve
# MOI.optimize!(model)
#
# @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
#
# @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint
#
# @test MOI.get(model, MOI.ObjectiveValue()) ≈ 0 atol=atol rtol=rtol
#
# @test MOI.get(model, MOI.VariablePrimal(), x[1]) ≈ 0 atol=atol rtol=rtol
# end
# end

function linear11test(model::MOI.ModelLike, config::MOIT.TestConfig)
atol = config.atol
rtol = config.rtol
# simple 2 variable, 1 constraint problem
# min x + y
# st x + y >= 1
# x + y >= 2

@test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
@test MOI.supports(model, MOI.ObjectiveSense())
@test MOI.supports_constraint(model, MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64})
@test MOI.supports_constraint(model, MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})

MOI.empty!(model)
@test MOI.is_empty(model)

v = MOI.add_variables(model, 2)

c1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,1.0], v), 0.0), MOI.GreaterThan(1.0))
c2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,1.0], v), 0.0), MOI.GreaterThan(2.0))

MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,1.0], v), 0.0))
MOI.set(model, MOI.ObjectiveSense(), MOI.MinSense)

if config.solve
MOI.optimize!(model)

@test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
@test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint
@test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.0 atol=atol rtol=rtol
end

c3 = MOI.transform(model, c2, MOI.LessThan(2.0))

@test isa(c3, MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}})
@test MOI.is_valid(model, c2) == false
@test MOI.is_valid(model, c3) == true

if config.solve
MOI.optimize!(model)

@test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
@test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint
@test MOI.get(model, MOI.ObjectiveValue()) ≈ 1.0 atol=atol rtol=rtol
end
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)
# linear15test(optimizer, config)
linear11test(optimizer, config)

# @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 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)
# end
end


return nothing