diff --git a/docs/src/manual/standard_form.md b/docs/src/manual/standard_form.md index f72a40c659..0111ec8da9 100644 --- a/docs/src/manual/standard_form.md +++ b/docs/src/manual/standard_form.md @@ -93,7 +93,9 @@ The matrix-valued set types implemented in MathOptInterface.jl are: | [`LogDetConeTriangle(d)`](@ref MathOptInterface.LogDetConeTriangle) | ``\{ (t,u,X) \in \mathbb{R}^{2+d(1+d)/2} : t \le u\log(\det(X/u)), X \mbox{ is the upper triangle of a PSD matrix}, u > 0 \}`` | | [`LogDetConeSquare(d)`](@ref MathOptInterface.LogDetConeSquare) | ``\{ (t,u,X) \in \mathbb{R}^{2+d^2} : t \le u \log(\det(X/u)), X \mbox{ is a PSD matrix}, u > 0 \}`` | | [`NormSpectralCone(r, c)`](@ref MathOptInterface.NormSpectralCone) | ``\{ (t, X) \in \mathbb{R}^{1 + r \times c} : t \ge \sigma_1(X), X \mbox{ is a } r\times c\mbox{ matrix} \}`` -| [`NormNuclearCone(r, c)`](@ref MathOptInterface.NormNuclearCone) | ``\{ (t, X) \in \mathbb{R}^{1 + r \times c} : t \ge \sum_i \sigma_i(X), X \mbox{ is a } r\times c\mbox{ matrix} \}`` +| [`NormNuclearCone(r, c)`](@ref MathOptInterface.NormNuclearCone) | ``\{ (t, X) \in \mathbb{R}^{1 + r \times c} : t \ge \sum_i \sigma_i(X), X \mbox{ is a } r\times c\mbox{ matrix} \}`` | +| [`HermitianPositiveSemidefiniteConeTriangle(d)`](@ref MathOptInterface.HermitianPositiveSemidefiniteConeTriangle) | The cone of Hermitian positive semidefinite matrices, with +`side_dimension` rows and columns. | Some of these cones can take two forms: `XXXConeTriangle` and `XXXConeSquare`. diff --git a/docs/src/reference/standard_form.md b/docs/src/reference/standard_form.md index 041138dec7..3a26183f33 100644 --- a/docs/src/reference/standard_form.md +++ b/docs/src/reference/standard_form.md @@ -141,6 +141,7 @@ List of recognized matrix sets. ```@docs PositiveSemidefiniteConeTriangle PositiveSemidefiniteConeSquare +HermitianPositiveSemidefiniteConeTriangle LogDetConeTriangle LogDetConeSquare RootDetConeTriangle diff --git a/docs/src/submodules/Bridges/list_of_bridges.md b/docs/src/submodules/Bridges/list_of_bridges.md index 26b24453eb..168951e672 100644 --- a/docs/src/submodules/Bridges/list_of_bridges.md +++ b/docs/src/submodules/Bridges/list_of_bridges.md @@ -84,4 +84,5 @@ Bridges.Variable.RSOCtoSOCBridge Bridges.Variable.SOCtoRSOCBridge Bridges.Variable.VectorizeBridge Bridges.Variable.ZerosBridge +Bridges.Variable.HermitianToSymmetricPSDBridge ``` diff --git a/src/Bridges/Variable/Variable.jl b/src/Bridges/Variable/Variable.jl index d91145788b..9db20bc32d 100644 --- a/src/Bridges/Variable/Variable.jl +++ b/src/Bridges/Variable/Variable.jl @@ -22,6 +22,7 @@ include("bridges/rsoc_soc.jl") include("bridges/soc_rsoc.jl") include("bridges/vectorize.jl") include("bridges/zeros.jl") +include("bridges/hermitian.jl") """ add_all_bridges(model, ::Type{T}) where {T} @@ -38,6 +39,7 @@ function add_all_bridges(model, ::Type{T}) where {T} MOI.Bridges.add_bridge(model, SOCtoRSOCBridge{T}) MOI.Bridges.add_bridge(model, RSOCtoSOCBridge{T}) MOI.Bridges.add_bridge(model, RSOCtoPSDBridge{T}) + MOI.Bridges.add_bridge(model, HermitianToSymmetricPSDBridge{T}) return end diff --git a/src/Bridges/Variable/bridges/hermitian.jl b/src/Bridges/Variable/bridges/hermitian.jl new file mode 100644 index 0000000000..f5566b4380 --- /dev/null +++ b/src/Bridges/Variable/bridges/hermitian.jl @@ -0,0 +1,334 @@ +# Copyright (c) 2017: Miles Lubin and contributors +# Copyright (c) 2017: Google Inc. +# +# Use of this source code is governed by an MIT-style license that can be found +# in the LICENSE.md file or at https://opensource.org/licenses/MIT. + +""" + HermitianToSymmetricPSDBridge{T} <: Bridges.Variable.AbstractBridge + +`HermitianToSymmetricPSDBridge` implements the following reformulation: + + * Hermitian positive semidefinite `n x n` complex matrix to a symmetric + positive semidefinite `2n x 2n` real matrix satisfying equality constraints + described below. + +## Source node + +`HermitianToSymmetricPSDBridge` supports: + + * [`MOI.VectorOfVariables`](@ref) in + [`MOI.HermitianPositiveSemidefiniteConeTriangle`](@ref) + +## Target node + +`HermitianToSymmetricPSDBridge` creates: + + * [`MOI.VectorOfVariables`](@ref) in [`MOI.PositiveSemidefiniteConeTriangle`](@ref) + * [`MOI.ScalarAffineFunction{T}`](@ref) in [`MOI.EqualTo{T}`](@ref) + +## Reformulation + +The reformulation is best described by example. + +The Hermitian matrix: +```math +\\begin{bmatrix} + x_{11} & x_{12} + y_{12}im & x_{13} + y_{13}im\\\\ + x_{12} - y_{12}im & x_{22} & x_{23} + y_{23}im\\\\ + x_{13} - y_{13}im & x_{23} - y_{23}im & x_{33} +\\end{bmatrix} +``` +is positive semidefinite if and only if the symmetric matrix: +```math +\\begin{bmatrix} + x_{11} & x_{12} & x_{13} & 0 & y_{12} & y_{13} \\\\ + & x_{22} & x_{23} & -y_{12} & 0 & y_{23} \\\\ + & & x_{33} & -y_{13} & -y_{23} & 0 \\\\ + & & & x_{11} & x_{12} & x_{13} \\\\ + & & & & x_{22} & x_{23} \\\\ + & & & & & x_{33} +\\end{bmatrix} +``` +is positive semidefinite. + +The bridge achieves this reformulation by adding a new set of variables in +`MOI.PositiveSemidefiniteConeTriangle(6)`, and then adding three groups of +equality constraints to: + + * constrain the two `x` blocks to be equal + * force the diagonal of the `y` blocks to be `0` + * force the lower triangular of the `y` block to be the negative of the upper + triangle. +""" +struct HermitianToSymmetricPSDBridge{T} <: AbstractBridge + variables::Vector{MOI.VariableIndex} + psd::MOI.ConstraintIndex{ + MOI.VectorOfVariables, + MOI.PositiveSemidefiniteConeTriangle, + } + n::Int + ceq::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}}} +end + +const HermitianToSymmetricPSD{T,OT<:MOI.ModelLike} = + SingleBridgeOptimizer{HermitianToSymmetricPSDBridge{T},OT} + +function bridge_constrained_variable( + ::Type{HermitianToSymmetricPSDBridge{T}}, + model::MOI.ModelLike, + set::MOI.HermitianPositiveSemidefiniteConeTriangle, +) where {T} + n = set.side_dimension + variables, psd_ci = MOI.add_constrained_variables( + model, + MOI.PositiveSemidefiniteConeTriangle(2n), + ) + ceq = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}}[] + k11 = 0 + k12 = k22 = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n)) + function X21(i, j) + I, J = j, n + i + k21 = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(J - 1)) + I + return variables[k21] + end + for j in 1:n + k22 += n + for i in 1:j + k11 += 1 + k12 += 1 + k22 += 1 + f_x = MOI.Utilities.operate(-, T, variables[k11], variables[k22]) + push!(ceq, MOI.add_constraint(model, f_x, MOI.EqualTo(zero(T)))) + if i == j # y_{ii} = 0 + f_0 = convert(MOI.ScalarAffineFunction{T}, variables[k12]) + push!(ceq, MOI.add_constraint(model, f_0, MOI.EqualTo(zero(T)))) + else # y_{ij} = -y_{ji} + f_y = MOI.Utilities.operate(+, T, X21(i, j), variables[k12]) + push!(ceq, MOI.add_constraint(model, f_y, MOI.EqualTo(zero(T)))) + end + end + k12 += n + end + return HermitianToSymmetricPSDBridge(variables, psd_ci, n, ceq) +end + +function supports_constrained_variable( + ::Type{<:HermitianToSymmetricPSDBridge}, + ::Type{MOI.HermitianPositiveSemidefiniteConeTriangle}, +) + return true +end + +function MOI.Bridges.added_constrained_variable_types( + ::Type{<:HermitianToSymmetricPSDBridge}, +) + return Tuple{Type}[(MOI.PositiveSemidefiniteConeTriangle,)] +end + +function MOI.Bridges.added_constraint_types( + ::Type{HermitianToSymmetricPSDBridge{T}}, +) where {T} + return Tuple{Type,Type}[(MOI.ScalarAffineFunction{T}, MOI.EqualTo{T})] +end + +function MOI.get(bridge::HermitianToSymmetricPSDBridge, ::MOI.NumberOfVariables) + return length(bridge.variables) +end + +function MOI.get( + bridge::HermitianToSymmetricPSDBridge, + ::MOI.ListOfVariableIndices, +) + return copy(bridge.variables) +end + +function MOI.get( + ::HermitianToSymmetricPSDBridge, + ::MOI.NumberOfConstraints{ + MOI.VectorOfVariables, + MOI.PositiveSemidefiniteConeTriangle, + }, +)::Int64 + return 1 +end + +function MOI.get( + bridge::HermitianToSymmetricPSDBridge, + ::MOI.ListOfConstraintIndices{ + MOI.VectorOfVariables, + MOI.PositiveSemidefiniteConeTriangle, + }, +) + return [bridge.psd] +end + +function MOI.get( + bridge::HermitianToSymmetricPSDBridge{T}, + ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}}, +) where {T} + return length(bridge.ceq) +end + +function MOI.get( + bridge::HermitianToSymmetricPSDBridge{T}, + ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}}, +) where {T} + return copy(bridge.ceq) +end + +function MOI.delete(model::MOI.ModelLike, bridge::HermitianToSymmetricPSDBridge) + MOI.delete(model, bridge.ceq) + MOI.delete(model, bridge.variables) + return +end + +function MOI.get( + ::MOI.ModelLike, + ::MOI.ConstraintSet, + bridge::HermitianToSymmetricPSDBridge, +) + return MOI.HermitianPositiveSemidefiniteConeTriangle(bridge.n) +end + +function _matrix_indices(k) + # If `k` is a diagonal index, `s(k)` is odd and 1 + 8k is a perfect square. + n = 1 + 8k + s = isqrt(n) + j = if s^2 == n + div(s, 2) + else + # Otherwise, if it is after the diagonal index `k` but before the + # diagonal index `k'` with `s(k') = s(k) + 2`, we have + # `s(k) <= s < s(k) + 2`. + # By shifting by `+1` before `div`, we make sure to have the right + # column. + div(s + 1, 2) + end + i = k - MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(j - 1)) + return i, j +end + +function _variable_map(idx::MOI.Bridges.IndexInVector, n) + N = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n)) + if idx.value <= N + return idx.value + end + i, j = _matrix_indices(idx.value - N) + d = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(j)) + return N + j * n + d + i +end + +function _variable( + bridge::HermitianToSymmetricPSDBridge, + i::MOI.Bridges.IndexInVector, +) + return bridge.variables[_variable_map(i, bridge.n)] +end + +function MOI.get( + model::MOI.ModelLike, + attr::MOI.ConstraintPrimal, + bridge::HermitianToSymmetricPSDBridge{T}, +) where {T} + values = MOI.get(model, attr, bridge.psd) + M = MOI.dimension(MOI.get(model, MOI.ConstraintSet(), bridge)) + n = bridge.n + return [values[_variable_map(MOI.Bridges.IndexInVector(i), n)] for i in 1:M] +end + +# We don't need to take into account the equality constraints. We just need to +# sum up (with appropriate +/-) each dual variable associated with the original +# x or y element. +# The reason for this is as follows: +# Suppose for simplicity that the elements of a `2n x 2n` matrix are ordered as: +# ``` +# \\ 1 |\\ 2 +# \\ | 3 +# \\ |4_\\ +# \\ 5 +# \\ +# \\ +# ``` +# Let `H = HermitianToSymmetricPSDBridge(n)`, +# `S = PositiveSemidefiniteConeTriangle(2n)` and `ceq` be the linear space of +# `2n x 2n` symmetric matrices such that the block `1` and `5` are equal, `2` and `4` are opposite and `3` is zero. +# We consider the cone `P = S ∩ ceq`. +# We have `P = A * H` where +# ``` +# [I 0] +# [0 I] +# A = [0 0] +# [0 -I] +# [I 0] +# ``` +# Therefore, `H* = A* * P*` where +# ``` +# [I 0 0 0 I] +# A* = [0 I 0 -I 0] +# ``` +# Moreover, as `(S ∩ T)* = S* + T*` for cones `S` and `T`, we have +# ``` +# P* = S* + ceq* +# ``` +# the dual vector of `P*` is the dual vector of `S*` for which we add in the corresponding +# entries the dual of the three constraints, multiplied by the coefficients for the `EqualTo` constraints. +# Note that these contributions cancel out when we multiply them by `A*`: +# A* * (S* + ceq*) = A* * S* +# so we can just ignore them. +function MOI.get( + model::MOI.ModelLike, + attr::MOI.ConstraintDual, + bridge::HermitianToSymmetricPSDBridge{T}, +) where {T} + dual = MOI.get(model, attr, bridge.psd) + M = MOI.dimension(MOI.get(model, MOI.ConstraintSet(), bridge)) + result = zeros(T, M) + n = bridge.n + N = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n)) + k11, k12, k22 = 0, N, N + k21 = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(2n)) + 1 + k = 0 + for j in 1:n + k21 -= n + 1 - j + k22 += n + for i in 1:j + k11 += 1 + k12 += 1 + k21 -= 1 + k22 += 1 + result[k11] += dual[k11] + dual[k22] + if i != j + k += 1 + result[N+k] += dual[k12] - dual[k21] + end + end + k12 += n + k21 -= n - j + end + return result +end + +function MOI.get( + model::MOI.ModelLike, + attr::MOI.VariablePrimal, + bridge::HermitianToSymmetricPSDBridge{T}, + i::MOI.Bridges.IndexInVector, +) where {T} + return MOI.get(model, attr, _variable(bridge, i)) +end + +function MOI.Bridges.bridged_function( + bridge::HermitianToSymmetricPSDBridge{T}, + i::MOI.Bridges.IndexInVector, +) where {T} + return convert(MOI.ScalarAffineFunction{T}, _variable(bridge, i)) +end + +function unbridged_map( + bridge::HermitianToSymmetricPSDBridge{T}, + vi::MOI.VariableIndex, + i::MOI.Bridges.IndexInVector, +) where {T} + return (_variable(bridge, i) => convert(MOI.ScalarAffineFunction{T}, vi),) +end diff --git a/src/Test/Test.jl b/src/Test/Test.jl index a632366b42..638820c887 100644 --- a/src/Test/Test.jl +++ b/src/Test/Test.jl @@ -6,6 +6,8 @@ module Test +import LinearAlgebra + using MathOptInterface const MOI = MathOptInterface const MOIU = MOI.Utilities diff --git a/src/Test/test_conic.jl b/src/Test/test_conic.jl index bee4e59eae..be9186ce45 100644 --- a/src/Test/test_conic.jl +++ b/src/Test/test_conic.jl @@ -6831,3 +6831,206 @@ function setup_test( ) return end + +""" + test_conic_HermitianPositiveSemidefiniteConeTriangle_1( + model::MOI.ModelLike, + config::Config, + ) + +Test the computation of the projection of a Hermitian matrix to the cone +of Hermitian positive semidefinite matrices. +The matrix is: +``` +[ 1 -1+im] [ 1-im ?] [√3 ] [ 1+im -1+√3] +[-1-im -1 ] = [-1+√3 ?] * [ -√3] * [ ? ? ] / (6 - 2√3) +``` +so it's projection to the Hermitian PSD cone is: +``` + [1-im] [ 1+im 1-√3] +√3 / (6 - 2√3) * [1-√3] +``` +which is: +``` + [1-im] [ 1+im 1-√3] +(1 + √3) / 4 * [1-√3] +``` +which is +``` +[(1+√3)/2 -1/2+im/2] +[-1/2+im/2 (-1+√3)/2] +``` +""" +function test_conic_HermitianPositiveSemidefiniteConeTriangle_1( + model::MOI.ModelLike, + config::Config{T}, +) where {T} + @requires _supports(config, MOI.optimize!) + @requires MOI.supports_add_constrained_variables( + model, + MOI.HermitianPositiveSemidefiniteConeTriangle, + ) + @requires MOI.supports_constraint( + model, + MOI.VectorAffineFunction{T}, + MOI.SecondOrderCone, + ) + set = MOI.HermitianPositiveSemidefiniteConeTriangle(2) + x, cx = MOI.add_constrained_variables(model, set) + t = MOI.add_variable(model) + soc_set = MOI.SecondOrderCone(5) + con_soc = MOI.add_constraint( + model, + MOI.Utilities.operate( + vcat, + T, + t, + x[1] - one(T), + √T(2) * (x[2] + one(T)), + x[3] + one(T), + √T(2) * (x[4] - one(T)), + ), + soc_set, + ) + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.set(model, MOI.ObjectiveFunction{typeof(t)}(), t) + MOI.optimize!(model) + primal = [ + (T(1) + √T(3)) / T(2), + -T(1) / T(2), + (-T(1) + √T(3)) / T(2), + T(1) / T(2), + ] + soc_primal = [ + primal[1] - one(T), + √T(2) * (primal[2] + one(T)), + primal[3] + one(T), + √T(2) * (primal[4] - one(T)), + ] + t_value = LinearAlgebra.norm(soc_primal) + soc_primal = [t_value; soc_primal] + dual = [ + (T(3) - √T(3)) / T(6), + √T(3) / T(6), + (T(3) + √T(3)) / T(6), + -√T(3) / T(6), + ] + @test ≈(MOI.Utilities.set_dot(primal, dual, set), T(0), config) + soc_dual = [one(T), -dual[1], -dual[2] * √T(2), -dual[3], -dual[4] * √T(2)] + @test ≈(MOI.Utilities.set_dot(soc_primal, soc_dual, soc_set), T(0), config) + @test ≈(MOI.get(model, MOI.VariablePrimal(), x), primal, config) + @test ≈(MOI.get(model, MOI.VariablePrimal(), t), t_value, config) + @test ≈(MOI.get(model, MOI.ConstraintPrimal(), cx), primal, config) + @test ≈(MOI.get(model, MOI.ConstraintPrimal(), con_soc), soc_primal, config) + if _supports(config, MOI.ConstraintDual) + @test ≈(MOI.get(model, MOI.ConstraintDual(), cx), dual, config) + @test ≈(MOI.get(model, MOI.ConstraintDual(), con_soc), soc_dual, config) + end + return +end + +function setup_test( + ::typeof(test_conic_HermitianPositiveSemidefiniteConeTriangle_1), + model::MOIU.MockOptimizer, + ::Config{T}, +) where {T} + primal = [ + (T(1) + √T(3)) / T(2), + -T(1) / T(2), + (-T(1) + √T(3)) / T(2), + T(1) / T(2), + ] + soc_primal = [ + primal[1] - one(T), + √T(2) * (primal[2] + one(T)), + primal[3] + one(T), + √T(2) * (primal[4] - one(T)), + ] + t_value = LinearAlgebra.norm(soc_primal) + push!(primal, t_value) + soc_primal = [t_value; soc_primal] + dual = [ + (T(3) - √T(3)) / T(6), + √T(3) / T(6), + (T(3) + √T(3)) / T(6), + -√T(3) / T(6), + ] + soc_dual = [one(T), -dual[1], -dual[2] * √T(2), -dual[3], -dual[4] * √T(2)] + MOIU.set_mock_optimize!( + model, + (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!( + mock, + primal, + (MOI.VectorAffineFunction{T}, MOI.SecondOrderCone) => + [soc_dual], + ), + ) + return +end + +function version_added( + ::typeof(test_conic_HermitianPositiveSemidefiniteConeTriangle_1), +) + return v"1.7.0" +end + +""" + test_conic_HermitianPositiveSemidefiniteConeTriangle_2( + model::MOI.ModelLike, + config::Config, + ) + +Test adding [`MathOptInterface.VariableIndex`](@ref)-in-[`MathOptInterface.EqualTo](@ref)` +on [`MathOptInterface.HermitianPositiveSemidefiniteConeTriangle`](@ref) +variables. +``` +""" +function test_conic_HermitianPositiveSemidefiniteConeTriangle_2( + model::MOI.ModelLike, + config::Config{T}, +) where {T} + @requires _supports(config, MOI.optimize!) + @requires MOI.supports_add_constrained_variables( + model, + MOI.HermitianPositiveSemidefiniteConeTriangle, + ) + @requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{T}) + set = MOI.HermitianPositiveSemidefiniteConeTriangle(3) + x, cx = MOI.add_constrained_variables(model, set) + primal = T[1, 0, 1, 0, 0, 0, -1, 0, 0] + MOI.add_constraints(model, x, MOI.EqualTo.(primal)) + MOI.optimize!(model) + @test ≈(MOI.get(model, MOI.VariablePrimal(), x), primal, config) + @test ≈(MOI.get(model, MOI.ConstraintPrimal(), cx), primal, config) + if _supports(config, MOI.ConstraintDual) + dual = zeros(T, 9) + @test ≈(MOI.get(model, MOI.ConstraintDual(), cx), dual, config) + end + return +end + +function setup_test( + ::typeof(test_conic_HermitianPositiveSemidefiniteConeTriangle_2), + model::MOIU.MockOptimizer, + ::Config{T}, +) where {T} + MOIU.set_mock_optimize!( + model, + (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!( + mock, + T[1, 0, 1, 0, 0, 0, -1, 0, 0], + ( + MOI.VectorOfVariables, + MOI.HermitianPositiveSemidefiniteConeTriangle, + ) => [zeros(T, 9)], + ), + ) + model.eval_variable_constraint_dual = false + return () -> model.eval_variable_constraint_dual = true +end + +function version_added( + ::typeof(test_conic_HermitianPositiveSemidefiniteConeTriangle_2), +) + return v"1.7.0" +end diff --git a/src/Utilities/results.jl b/src/Utilities/results.jl index 10eb39b271..1646199381 100644 --- a/src/Utilities/results.jl +++ b/src/Utilities/results.jl @@ -546,6 +546,19 @@ function set_dot( return triangle_dot(x, y, MOI.side_dimension(set), 0) end +function MOI.Utilities.set_dot( + x::AbstractVector, + y::AbstractVector, + set::MOI.HermitianPositiveSemidefiniteConeTriangle, +) + sym = MOI.PositiveSemidefiniteConeTriangle(set.side_dimension) + result = set_dot(x, y, sym) + for k in (MOI.dimension(sym)+1):MOI.dimension(set) + result = MA.add_mul!!(result, 2, x[k], y[k]) + end + return result +end + function set_dot( x::AbstractVector, y::AbstractVector, @@ -597,6 +610,18 @@ function dot_coefficients( return b end +function dot_coefficients( + a::AbstractVector, + set::MOI.HermitianPositiveSemidefiniteConeTriangle, +) + sym = MOI.PositiveSemidefiniteConeTriangle(set.side_dimension) + b = dot_coefficients(a, sym) + for k in (MOI.dimension(sym)+1):MOI.dimension(set) + b[k] /= 2 + end + return b +end + function dot_coefficients(a::AbstractVector, set::MOI.RootDetConeTriangle) b = copy(a) triangle_coefficients!(b, set.side_dimension, 1) diff --git a/src/sets.jl b/src/sets.jl index 14a213006c..680d34646d 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -795,6 +795,46 @@ function triangular_form(::Type{PositiveSemidefiniteConeSquare}) return PositiveSemidefiniteConeTriangle end +""" + HermitianPositiveSemidefiniteConeTriangle(side_dimension) <: AbstractVectorSet + +The (vectorized) cone of Hermitian positive semidefinite matrices, with +`side_dimension` rows and columns. + +Becaue the matrix is Hermitian, the diagonal elements are real, and the +complex-valued lower triangular entries are obtained as the conjugate of +corresponding upper triangular entries. + +### Vectorization format + +The vectorized form starts with real part of the entries of the upper triangular +part of the matrix, given column by column as explained in +[`AbstractSymmetricMatrixSetSquare`](@ref). + +It is then followed by the imaginary part of the off-diagonal entries of the +upper triangular part, also given column by column. + +For example, the matrix +```math +\\begin{bmatrix} + 1 & 2 + 7im & 4 + 8im\\\\ + 2 - 7im & 3 & 5 + 9im\\\\ + 4 - 8im & 5 - 9im & 6 +\\end{bmatrix} +``` +has [`side_dimension`](@ref) 3 and is represented as the vector +``[1, 2, 3, 4, 5, 6, 7, 8, 9]``. +""" +struct HermitianPositiveSemidefiniteConeTriangle <: AbstractVectorSet + side_dimension::Int +end + +function dimension(set::HermitianPositiveSemidefiniteConeTriangle) + real_nnz = div(set.side_dimension * (set.side_dimension + 1), 2) + imag_nnz = div((set.side_dimension - 1) * set.side_dimension, 2) + return real_nnz + imag_nnz +end + """ LogDetConeTriangle(side_dimension) @@ -1559,6 +1599,7 @@ function Base.copy( NormNuclearCone, PositiveSemidefiniteConeTriangle, PositiveSemidefiniteConeSquare, + HermitianPositiveSemidefiniteConeTriangle, LogDetConeTriangle, LogDetConeSquare, RootDetConeTriangle, diff --git a/test/Bridges/Variable/hermitian.jl b/test/Bridges/Variable/hermitian.jl new file mode 100644 index 0000000000..b171b429e9 --- /dev/null +++ b/test/Bridges/Variable/hermitian.jl @@ -0,0 +1,119 @@ +# Copyright (c) 2017: Miles Lubin and contributors +# Copyright (c) 2017: Google Inc. +# +# Use of this source code is governed by an MIT-style license that can be found +# in the LICENSE.md file or at https://opensource.org/licenses/MIT. + +module TestVariableHermitianToSymmetricPSD + +using Test + +using MathOptInterface +const MOI = MathOptInterface + +import LinearAlgebra + +function runtests() + for name in names(@__MODULE__; all = true) + if startswith("$(name)", "test_") + @testset "$(name)" begin + getfield(@__MODULE__, name)() + end + end + end + return +end + +include("../utilities.jl") + +function test_conic_HermitianPositiveSemidefiniteConeTriangle_1() + T = Float64 + mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{T}()) + bridged_mock = MOI.Bridges.Variable.HermitianToSymmetricPSD{T}(mock) + primal = [ + (T(1) + √T(3)) / T(2), + -T(1) / T(2), + (-T(1) + √T(3)) / T(2), + T(1) / T(2), + ] + soc_primal = [ + primal[1] - one(T), + √T(2) * (primal[2] + one(T)), + primal[3] + one(T), + √T(2) * (primal[4] - one(T)), + ] + t_value = LinearAlgebra.norm(soc_primal) + real_primal = [ + primal[1:3] + zero(T) + -primal[4] + primal[1] + primal[4] + zero(T) + primal[2:3] + t_value + ] + dual = [ + (T(3) - √T(3)) / T(6), + √T(3) / T(6), + (T(3) + √T(3)) / T(6), + -√T(3) / T(6), + ] + soc_dual = [one(T), -dual[1], -dual[2] * √T(2), -dual[3], -dual[4] * √T(2)] + mock.optimize! = + (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!( + mock, + real_primal, + (MOI.VectorAffineFunction{T}, MOI.SecondOrderCone) => + [soc_dual], + (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}) => [ + dual[1] / T(2) + zero(T) + dual[2] + -dual[2] + dual[3] / T(2) + zero(T) + ], + ) + MOI.empty!(bridged_mock) + MOI.Test.test_conic_HermitianPositiveSemidefiniteConeTriangle_1( + bridged_mock, + MOI.Test.Config(), + ) + return +end + +function test_runtests() + MOI.Bridges.runtests( + MOI.Bridges.Variable.HermitianToSymmetricPSDBridge, + """ + variables: a, b, c + constrainedvariable: [r11, r12, r22, c12] in HermitianPositiveSemidefiniteConeTriangle(2) + 1.0 * r11 >= 1.0 + 1.0 * r12 >= 2.0 + 1.0 * r22 >= 3.0 + 1.0 * c12 >= 4.0 + [a, b, c] in PositiveSemidefiniteConeTriangle(2) + """, + """ + variables: a, b, c + constrainedvariable: [v11, v12, v22, v13, v23, v33, v14, v24, v34, v44] in PositiveSemidefiniteConeTriangle(4) + 1.0 * v11 >= 1.0 + 1.0 * v12 >= 2.0 + 1.0 * v22 >= 3.0 + 1.0 * v14 >= 4.0 + 1.0 * v11 + -1.0 * v33 == 0.0 + 1.0 * v13 == 0.0 + 1.0 * v12 + -1.0 * v34 == 0.0 + 1.0 * v23 + 1.0 * v14 == 0.0 + 1.0 * v22 + -1.0 * v44 == 0.0 + 1.0 * v24 == 0.0 + [a, b, c] in PositiveSemidefiniteConeTriangle(2) + """, + ) + return +end + +end # module + +TestVariableHermitianToSymmetricPSD.runtests()