Add complex PSD cone (#194)

* add complex psd cone

* silence test warnings

* add unit test for complex psd cone

* disable chordal decomposition for complex psd cones

* make test deterministic

* mention complex psd cone in getting started

* better docstring

docs/src/getting_started.md

@@ -31,8 +31,9 @@ ZeroSet | The set ``\{ 0 \}^{dim}`` that contains the origin
 Nonnegatives | The nonnegative orthant ``\{ x \in \mathbb{R}^{dim} : x_i \ge 0, \forall i=1,\dots,\mathrm{dim} \}``
 Box(l, u) | The hyperbox ``\{ x \in \mathbb{R}^{dim} : l \leq x \leq u\}`` with vectors ``l \in \mathbb{R}^{dim} \cup \{-\infty\}`` and ``u \in \mathbb{R}^{dim} \cup \{+\infty\}``
 SecondOrderCone | The second-order (Lorenz) cone ``\{ (t,x) \in \mathbb{R}^{dim}  :  \|x\|_2   \leq t \}``
-PsdCone | The vectorized positive semidefinite cone ``\mathcal{S}_+^{dim}``. ``x`` is the vector obtained by stacking the columns of the positive semidefinite matrix ``X``, i.e. ``X \in \mathbb{S}^{\sqrt{dim}}_+ \rarr \text{vec}(X) = x \in \mathcal{S}_+^{dim}``
-PsdConeTriangle | The vectorized positive semidefinite cone ``\mathcal{S}_+^{dim}``. ``x`` is the vector obtained by stacking the columns of the upper triangular part of the positive semidefinite matrix ``X``, i.e. ``X \in \mathbb{S}^{d}_+ \rarr \text{svec}(X) = x \in \mathcal{S}_+^{dim}`` where ``d=\sqrt{1/4 + 2 \text{dim}} - 1/2``
+PsdCone | The vectorized positive real semidefinite cone ``\mathcal{S}_+^{dim}``. ``x`` is the vector obtained by stacking the columns of the positive semidefinite matrix ``X``, i.e. ``X \in \mathbb{S}^{\sqrt{dim}}_+ \rarr \text{vec}(X) = x \in \mathcal{S}_+^{dim}``
+PsdConeTriangle (real) | The vectorized real positive semidefinite cone ``\mathcal{S}_+^{dim}``. ``x`` is the vector obtained by stacking the columns of the upper triangular part of the positive semidefinite matrix ``X`` and scaling the off-diagonals by ``\sqrt{2}``, i.e. ``X \in \mathbb{S}^{d}_+ \rarr \text{svec}(X) = x \in \mathcal{S}_+^{dim}`` where ``d=\sqrt{1/4 + 2 \text{dim}} - 1/2``
+PsdConeTriangle (complex) | The vectorized complex positive semidefinite cone ``\mathcal{S}_+^{dim}``. ``x`` is the vector obtained by stacking the real and imaginary parts of the columns of the upper triangular part of the positive semidefinite matrix ``X`` and scaling the off-diagonals by ``\sqrt{2}``. The ordering follows [the one from MathOptInterface](https://jump.dev/MathOptInterface.jl/stable/reference/standard_form/#MathOptInterface.HermitianPositiveSemidefiniteConeTriangle). I.e. ``X \in \mathbb{H}^{\sqrt{dim}}_+ \rarr \text{svec}(X) = x \in \mathcal{S}_+^{dim}``.
 ExponentialCone | The exponential cone ``\mathcal{K}_{exp} = \{(x, y, z) \mid y \geq 0,  ye^{x/y} ≤ z\} \cup \{ (x,y,z) \mid   x \leq 0, y = 0, z \geq 0 \}``
 DualExponentialCone | The dual exponential cone ``\mathcal{K}^*_{exp} = \{(x, y, z) \mid x < 0,  -xe^{y/x} \leq e^1 z \} \cup \{ (0,y,z) \mid   y \geq 0, z \geq 0 \}``
 PowerCone(alpha) | The 3d power cone ``\mathcal{K}_{pow} = \{(x, y, z) \mid x^\alpha y^{(1-\alpha)} \geq  \|z\|, x \geq 0, y \geq 0 \}`` with ``0 < \alpha < 1``

src/COSMO.jl

@@ -10,6 +10,7 @@ using Reexport
 
 export assemble!, warm_start!, empty_model!, update!
 
+const RealOrComplex{T <: Real} = Union{T, Complex{T}}
 const DefaultFloat = Float64
 const DefaultInt   = LinearAlgebra.BlasInt
 

src/MOI_wrapper.jl

@@ -20,7 +20,8 @@ const IntervalConvertible = Union{GreaterThan, Interval}
 const Zeros = MOI.Zeros
 const SOC = MOI.SecondOrderCone
 const PSDT = MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle}
-const SupportedVectorSets = Union{Zeros, MOI.Nonnegatives, SOC, PSDT, MOI.ExponentialCone, MOI.DualExponentialCone, MOI.PowerCone, MOI.DualPowerCone}
+const ComplexPSDT = MOI.Scaled{MOI.HermitianPositiveSemidefiniteConeTriangle}
+const SupportedVectorSets = Union{Zeros, MOI.Nonnegatives, SOC, PSDT, ComplexPSDT, MOI.ExponentialCone, MOI.DualExponentialCone, MOI.PowerCone, MOI.DualPowerCone}
 const AggregatedSets = Union{Zeros, MOI.Nonnegatives, MOI.GreaterThan}
 aggregate_equivalent(::Type{<: MOI.Zeros}) = COSMO.ZeroSet
 aggregate_equivalent(::Type{<: Union{MOI.GreaterThan, MOI.Nonnegatives}}) = COSMO.Nonnegatives
@@ -449,7 +450,12 @@ end
 
 
 function processSet!(b::Vector{T}, rows::UnitRange{Int}, cs, s::PSDT) where {T <: AbstractFloat}
-    push!(cs, COSMO.PsdConeTriangle{T}(length(rows)))
+    push!(cs, COSMO.PsdConeTriangle{T, T}(length(rows)))
+    nothing
+end
+
+function processSet!(b::Vector{T}, rows::UnitRange{Int}, cs, s::ComplexPSDT) where {T <: AbstractFloat}
+    push!(cs, COSMO.PsdConeTriangle{T, Complex{T}}(length(rows)))
     nothing
 end
 

src/algebra.jl

@@ -222,23 +222,23 @@ function symmetrize_full!(A::AbstractMatrix{T}) where {T <: AbstractFloat}
 	nothing
 end
 
-# this function assumes real symmetric X and only considers the upper triangular part
-function is_pos_def!(X::AbstractMatrix{T}, tol::T=zero(T)) where T
+# this function assumes real or complex Hermitian X and only considers the upper triangular part
+function is_pos_def!(X::AbstractMatrix{<:RealOrComplex{T}}, tol::T=zero(T)) where {T}
 	# See https://math.stackexchange.com/a/13311
 	@inbounds for i = 1:size(X, 1)
 		X[i, i] += tol
 	end
-	F = cholesky!(Symmetric(X), check = false)
+	F = cholesky!(Hermitian(X), check = false)
 	return issuccess(F)
 end
 
-function is_neg_def!(X::AbstractMatrix{T}, tol::T=zero(T)) where T
+function is_neg_def!(X::AbstractMatrix{<:RealOrComplex{T}}, tol::T=zero(T)) where {T}
 	@. X *= -one(T)
 	return is_pos_def!(X, tol)
 end
 
-is_pos_def(X::AbstractMatrix{T}, tol::T=zero(T)) where T = is_pos_def!(copy(X), tol)
-is_neg_def(X::AbstractMatrix{T}, tol::T=zero(T)) where T = is_pos_def!(-X, tol)
+is_pos_def(X::AbstractMatrix{<:RealOrComplex{T}}, tol::T=zero(T)) where {T} = is_pos_def!(copy(X), tol)
+is_neg_def(X::AbstractMatrix{<:RealOrComplex{T}}, tol::T=zero(T)) where {T} = is_pos_def!(-X, tol)
 
 
 "Round x to the closest multiple of N."

src/chordal_decomposition/chordal_decomposition.jl

@@ -75,7 +75,8 @@ function _analyse_sparsity_pattern(ci::ChordalInfo{T}, csp::Array{Int, 1}, sets:
 end
 
 DenseEquivalent(C::COSMO.PsdCone{T}, dim::Int) where {T <: AbstractFloat} = COSMO.DensePsdCone{T}(dim)
-DenseEquivalent(C::COSMO.PsdConeTriangle{T}, dim::Int) where {T <: AbstractFloat} = COSMO.DensePsdConeTriangle{T}(dim)
+DenseEquivalent(C::COSMO.PsdConeTriangle{T, T}, dim::Int) where {T <: AbstractFloat} = COSMO.DensePsdConeTriangle{T, T}(dim)
+DenseEquivalent(C::COSMO.PsdConeTriangle{T, Complex{T}}, dim::Int) where {T <: AbstractFloat} = COSMO.DensePsdConeTriangle{T, Complex{T}}(dim)
 
 function nz_rows(a::SparseMatrixCSC{T}, ind::UnitRange{Int}, DROP_ZEROS_FLAG::Bool) where {T <: AbstractFloat}
   DROP_ZEROS_FLAG && dropzeros!(a)

src/constraint.jl

@@ -93,7 +93,17 @@ function Constraint{T}(
 	# this constructor doesnt work with cones that need special arguments like the power cone
 	set_type <: ArgumentCones && error("You can't create a constraint by passing the convex set as a type, if your convex set is a $(set_type). Please pass an object.")
 	# call the appropriate set constructor
-	convex_set = set_type{T}(size(A, 1))
+    n = size(A, 1)
+    # we can deduce whether the PsdCone must be real or complex from the dimension
+	if set_type <: PsdConeTriangles
+	    if n == 1 || isqrt(n)^2 != n
+            convex_set = set_type{T, T}(n)
+        else
+            convex_set = set_type{T, Complex{T}}(n)
+        end
+    else
+        convex_set = set_type{T}(n)
+    end
 	Constraint{T}(A, b, convex_set, args...)
 end