diff --git a/src/Bridges/Bridges.jl b/src/Bridges/Bridges.jl index 83900237b7..2e21fe5a33 100644 --- a/src/Bridges/Bridges.jl +++ b/src/Bridges/Bridges.jl @@ -43,6 +43,7 @@ Returns a `LazyBridgeOptimizer` bridging `model` for every bridge defined in thi function fullbridgeoptimizer(model::MOI.ModelLike, ::Type{T}) where T bridgedmodel = MOIB.LazyBridgeOptimizer(model, AllBridgedConstraints{T}()) add_bridge(bridgedmodel, MOIB.SplitIntervalBridge{T}) + add_bridge(bridgedmodel, MOIB.QuadtoSOCBridge{T}) add_bridge(bridgedmodel, MOIB.GeoMeanBridge{T}) add_bridge(bridgedmodel, MOIB.SquarePSDBridge{T}) add_bridge(bridgedmodel, MOIB.LogDetBridge{T}) @@ -57,6 +58,8 @@ include("intervalbridge.jl") @bridge SplitInterval SplitIntervalBridge () (MOI.Interval,) () () (MOI.SingleVariable,) (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction) () () include("rsocbridge.jl") @bridge RSOC RSOCBridge () () (MOI.RotatedSecondOrderCone,) () () () (MOI.VectorOfVariables,) (MOI.VectorAffineFunction,) +include("quadtosocbridge.jl") +@bridge QuadtoSOC QuadtoSOCBridge () (MOI.LessThan, MOI.GreaterThan) () () () (MOI.ScalarQuadraticFunction,) () () include("geomeanbridge.jl") @bridge GeoMean GeoMeanBridge () () (MOI.GeometricMeanCone,) () () () (MOI.VectorOfVariables,) (MOI.VectorAffineFunction,) include("squarepsdbridge.jl") diff --git a/src/Bridges/quadtosocbridge.jl b/src/Bridges/quadtosocbridge.jl new file mode 100644 index 0000000000..35278db383 --- /dev/null +++ b/src/Bridges/quadtosocbridge.jl @@ -0,0 +1,202 @@ +using Compat.LinearAlgebra, Compat.SparseArrays + +""" + QuadtoSOCBridge{T} + +The set of points `x` satisfying the constraint +```math +\\frac{1}{2}x^T Q x + a^T x + b \\le 0 +``` +is a convex set if `Q` is positive semidefinite and is the union of two convex +cones if `a` and `b` are zero (i.e. *homogeneous* case) and `Q` has only one +negative eigenvalue. Currently, only the non-homogeneous transformation +is implemented, see the [Note](@ref) section for more details. + +## Non-homogeneous case + +If `Q` is positive semidefinite, there exists `U` such that ``Q = U^T U``, the +inequality can then be rewritten as +```math +\\|U x\\|_2 \\le 2 (a^T x + b) +``` +which is equivalent to the membership of `(1, a^T x + b, Ux)` to the rotated +second-order cone. + +## Homogeneous case + +If `Q` has only one negative eigenvalue, the set of `x` such that ``x^T Q x \\le +0`` is the union of a convex cone and its opposite. We can choose which one to +model by checking the existence of bounds on variables as shown below. + +### Second-order cone + +If `Q` is diagonal and has eigenvalues `(1, 1, -1)`, the inequality +``x^2 + x^2 \\le z^2`` combined with ``z \\ge 0`` defines the Lorenz cone (i.e. +the second-order cone) but when combined with ``z \\le 0``, it gives the +opposite of the second order cone. Therefore, we need to check if the variable +`z` has a lower bound 0 or an upper bound 0 in order to determine which cone is + +### Rotated second-order cone + +The matrix `Q` corresponding to the inequality ``x^2 \\le 2yz`` has one +eigenvalue 1 with eigenvectors `(1, 0, 0)` and `(0, 1, -1)` and one eigenvalue +`-1` corresponding to the eigenvector `(0, 1, 1)`. Hence if we intersect this +union of two convex cone with the halfspace ``x + y \\ge 0``, we get the rotated +second-order cone and if we intersect it with the halfspace ``x + y \\le 0`` we +get the opposite of the rotated second-order cone. Note that `y` and `z` have +the same sign since `yz` is nonnegative hence ``x + y \\ge 0`` is equivalent to +``x \\ge 0`` and ``y \\ge 0``. + +### Note + +The check for existence of bound can be implemented (but inefficiently) with the +current interface but if bound is removed or transformed (e.g. `≤ 0` transformed +into `≥ 0`) then the bridge is no longer valid. For this reason the homogeneous +version of the bridge is not implemented yet. +""" +struct QuadtoSOCBridge{T} <: AbstractBridge + soc::CI{MOI.VectorAffineFunction{T}, MOI.RotatedSecondOrderCone} + dimension::Int # dimension of the SOC constraint + less_than::Bool # whether the constraint was ≤ or ≥ + set_constant::T # the constant that was on the set +end +function QuadtoSOCBridge{T}(model, func::MOI.ScalarQuadraticFunction{T}, + set::Union{MOI.LessThan{T}, + MOI.GreaterThan{T}}) where T + set_constant = MOIU.getconstant(set) + less_than = set isa MOI.LessThan + if !less_than + set_constant = -set_constant + end + Q, index_to_variable_map = matrix_from_quadratic_terms(func.quadratic_terms) + if !less_than + Compat.rmul!(Q, -1) + end + # We have L × L' ≈ Q[p, p] + L, p = try + @static if VERSION >= v"0.7-" + F = cholesky(Symmetric(Q)) + sparse(F.L), F.p + else + F = cholfact(Symmetric(Q)) + sparse(F[:L]), F[:p] + end + catch err + if err isa PosDefException + error("The optimizer supports second-order cone constraints and", + " not quadratic constraints but you entered a quadratic", + " constraint of type: `$(typeof(func))`-in-`$(typeof(set))`.", + " A bridge attempted to transform the quadratic constraint", + " to a second order cone constraint but the constraint is", + " not strongly convex, i.e., the symmetric matrix of", + " quadratic coefficients is not positive definite. Convex", + " constraints that are not strongly convex, i.e. the matrix is", + " positive semidefinite but not positive definite, are not", + " supported yet.") + else + rethrow(err) + end + end + Ux_terms = matrix_to_vector_affine_terms(L, p, index_to_variable_map) + Ux = MOI.VectorAffineFunction(Ux_terms, zeros(T, size(L, 2))) + t = MOI.ScalarAffineFunction(less_than ? MOIU.operate_terms(-, func.affine_terms) : func.affine_terms, + less_than ? set_constant - func.constant : func.constant - set_constant) + f = MOIU.operate(vcat, T, one(T), t, Ux) + dimension = MOI.output_dimension(f) + soc = MOI.add_constraint(model, f, + MOI.RotatedSecondOrderCone(dimension)) + return QuadtoSOCBridge(soc, dimension, less_than, set_constant) +end + +function matrix_from_quadratic_terms(terms::Vector{MOI.ScalarQuadraticTerm{T}}) where T + variable_to_index_map = Dict{VI, Int}() + index_to_variable_map = Dict{Int, VI}() + n = 0 + for term in terms + for variable in (term.variable_index_1, term.variable_index_2) + if !(variable in keys(variable_to_index_map)) + n += 1 + variable_to_index_map[variable] = n + index_to_variable_map[n] = variable + end + end + end + I = Int[] + J = Int[] + V = T[] + for term in terms + i = variable_to_index_map[term.variable_index_1] + j = variable_to_index_map[term.variable_index_2] + push!(I, i) + push!(J, j) + push!(V, term.coefficient) + if i != j + push!(I, i) + push!(J, j) + push!(V, term.coefficient) + end + end + # Duplicate terms are summed together in `sparse` + Q = sparse(I, J, V, n, n) + return Q, index_to_variable_map +end + +function matrix_to_vector_affine_terms(L::SparseMatrixCSC{T}, + p::Vector, + index_to_variable_map::Dict{Int, VI}) where T + # We know that L × L' ≈ Q[p, p] hence (L × L')[i, :] ≈ Q[p[i], p] + # We precompute the map to avoid having to do a dictionary lookup for every + # term + variable = map(i -> index_to_variable_map[p[i]], 1:size(L, 1)) + function term(i::Integer, j::Integer, v::T) + return MOI.VectorAffineTerm(j, MOI.ScalarAffineTerm(v, variable[i])) + end + return map(ijv -> term(ijv...), zip(findnz(L)...)) +end + +function MOI.supports_constraint(::Type{QuadtoSOCBridge{T}}, + ::Type{MOI.ScalarQuadraticFunction{T}}, + ::Type{<:Union{MOI.LessThan{T}, + MOI.GreaterThan{T}}}) where T + return true +end +function added_constraint_types(::Type{QuadtoSOCBridge{T}}) where T + return [(MOI.VectorAffineFunction{T}, MOI.RotatedSecondOrderCone)] +end +function concrete_bridge_type(::Type{<:QuadtoSOCBridge{T}}, + ::Type{MOI.ScalarQuadraticFunction{T}}, + ::Type{<:Union{MOI.LessThan{T}, + MOI.GreaterThan{T}}}) where T + return QuadtoSOCBridge{T} +end + +# Attributes, Bridge acting as an model +function MOI.get(::QuadtoSOCBridge{T}, + ::MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, + MOI.RotatedSecondOrderCone}) where T + return 1 +end +function MOI.get(bridge::QuadtoSOCBridge{T}, + ::MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, + MOI.RotatedSecondOrderCone}) where T + return [bridge.soc] +end + +# References +function MOI.delete(model::MOI.ModelLike, bridge::QuadtoSOCBridge) + MOI.delete(model, bridge.soc) +end + +# Attributes, Bridge acting as a constraint +function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal, + bridge::QuadtoSOCBridge) + soc = MOI.get(model, attr, bridge.soc) + output = sum(soc[i]^2 for i in 3:bridge.dimension) + output /= 2 + output -= soc[1] * soc[2] + if !bridge.less_than + output = -output + end + output += bridge.set_constant + return output +end diff --git a/test/bridge.jl b/test/bridge.jl index 90da0bff4d..b39acfa656 100644 --- a/test/bridge.jl +++ b/test/bridge.jl @@ -255,6 +255,40 @@ end test_delete_bridge(bridgedmock, ci, 2, ((MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone, 0),)) end + @testset "QuadtoSOC" begin + bridgedmock = MOIB.QuadtoSOC{Float64}(mock) + @testset "Error for non-convex quadratic constraints" begin + x = MOI.add_variable(bridgedmock) + @test_throws ErrorException begin + MOI.add_constraint(bridgedmock, + MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[], + [MOI.ScalarQuadraticTerm(1.0, x, x)], + 0.0), + MOI.GreaterThan(0.0)) + end + @test_throws ErrorException begin + MOI.add_constraint(bridgedmock, + MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[], + [MOI.ScalarQuadraticTerm(-1.0, x, x)], + 0.0), + MOI.LessThan(0.0)) + end + end + MOIU.set_mock_optimize!(mock, + (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/2, 7/4], MOI.FeasiblePoint)) + MOIT.qcp1test(bridgedmock, config) + MOIU.set_mock_optimize!(mock, + (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [√2], MOI.FeasiblePoint)) + MOIT.qcp2test(bridgedmock, config) + MOIU.set_mock_optimize!(mock, + (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [√2], MOI.FeasiblePoint)) + MOIT.qcp3test(bridgedmock, config) + ci = first(MOI.get(bridgedmock, + MOI.ListOfConstraintIndices{MOI.ScalarQuadraticFunction{Float64}, + MOI.LessThan{Float64}}())) + test_delete_bridge(bridgedmock, ci, 1, ((MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),)) + end + @testset "SquarePSD" begin bridgedmock = MOIB.SquarePSD{Float64}(mock) mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, ones(4),