Add complex PSD cone

[araujoms]

Closes #193.

I changed the typedefs of PsdConeTriangle and DensePsdConeTriangle to accept an extra parameter, which determines whether the cone is real or complex. I generalized all relevant code to handle real or complex matrices. Also I added the new cones to the list of supported by MOI.

I did some benchmarking for a simple test problem that computes the minimal eigenvalue of a random complex Hermitian matrix. The orange curve is the times for the native complex version as a function of the dimension, whereas the blue curve refers to the the same problem mapped onto real Hermitian matrices via the usual technique c -> [real(c) imag(c); -imag(c) real(c)]:

As you can see, the difference is brutal. The code I used for benchmark is the one below:

using LinearAlgebra
using COSMO
import Random

function random_herm(::Type{T}, d) where {T}
    a = randn(T, d,d)
    return Hermitian(a + a')
end

function benchmark(::Type{T}, d) where {T}
    Random.seed!(2988)
    c = random_herm(complex(T), d)
    if T <: Real
        real_c = Hermitian([real(c) imag(c); -imag(c) real(c)])
        return mineig_raw(real_c)
    else
        return mineig_raw(c)        
    end
end

tri(d) = div(d*(d+1), 2)
    
function mineig_raw(c::AbstractMatrix{R}) where {R}
    d = size(c, 1)
    T = real(R)
    vec_dim = R <: Complex ? d^2 : tri(d)    

    model = COSMO.Model{T}()
    vec_c = zeros(T, vec_dim)
    COSMO.extract_upper_triangle!(c, vec_c, sqrt(T(2)))

    P = zeros(T, vec_dim, vec_dim)

    id_vec = zeros(T, vec_dim)
    diagonal_indices = tri.(1:d)
    id_vec[diagonal_indices] .= 1

    cs1 = COSMO.Constraint(id_vec', -T(1), COSMO.ZeroSet)
    cs2 = COSMO.Constraint(Matrix(T(1)*I(vec_dim)), zeros(T,vec_dim), COSMO.PsdConeTriangle{T, R}(vec_dim))
    constraints = [cs1; cs2]

    assemble!(model, P, vec_c, constraints, settings = COSMO.Settings{T}(verbose = false, eps_abs = 1e-3, eps_rel = 1e-3, eps_prim_inf = 1e-3, eps_dual_inf = 1e-3))
    result = COSMO.optimize!(model)
#    display(minimum(eigvals(c)))
#    display(result.obj_val)
    return result.times.solver_time
end

I've also checked that it works with JuMP with the following code:

using JuMP
using LinearAlgebra
using COSMO

function mineig(c::AbstractMatrix{R}) where {R}
    d = size(c, 1)
    T = real(R)

    model = GenericModel{T}(COSMO.Optimizer{T})
    if R <: Complex
        @variable(model, rho[1:d,1:d] in HermitianPSDCone())
    else
        @variable(model, rho[1:d,1:d] in PSDCone())
    end
    @constraint(model, tr(rho) == 1)
    @objective(model, Min, real(dot(c,rho)))
    set_attribute(model, "verbose", true)   
    optimize!(model)
    
    display(minimum(eigvals(c)))
    return objective_value(model)
end

I noticed two problems: COSMO often fails to converge, hitting the iteration limit (the peaks of the graph above), and doesn't work for Float32. Both problems were there before my changes, though.

☑️ PR checklist

• I have linked to the issue nr that this PR aims to resolve.
• If this PR is not yet ready to be merged, I have set it to "Draft".
• I have added tests for new features or the bug that this PR aims to fix.
• I have provided docstrings for new functions.
• I updated the repository documentation if these code changes touch the interface.

[CLAassistant]

All committers have signed the CLA.

[migarstka]

Great. Thank you @araujoms for your effort on this.

I am a bit occupied with my day job this week, but I try to review either one evening or on the weekend.

Is there a unit test for the complex cone that you could add?

[araujoms]

Sure, I've added a unit test using the native interface. I could also add one using JuMP, but there doesn't seem to be any.

[migarstka]

I can see one test failure on Julia 1.6 🤔 :

     Testing Running tests...
Complex PSD Cone: Test Failed at /home/runner/work/COSMO.jl/COSMO.jl/test/UnitTests/least_eigenvalue.jl:35
  Expression: ≈(mineig_raw(X), minimum(eigvals(X)), atol = 0.0001, rtol = 0.0001)
   Evaluated: -0.6468784469158599 ≈ -1.9271575272948176 (atol=0.0001, rtol=0.0001)

[araujoms]

Well you could drop support for 1.6. It runs fine on 1.9, and 1.10 will soon be the new LTS anyway.

[araujoms]

As it turns out Julia 1.6 is innocent, I just ran the tests locally, and they went fine. I guess it was just a random instability, I made the test deterministic to avoid that.

[migarstka]

any particular reason why you comment this file out?

[araujoms]

It is already included in runtests.jl. The duplicate inclusion was generating a lot of warnings about function redefinitions.

[migarstka]

It looks good so far. We should also point potential users at this feature in the docs. The section on the PsdConeTriangle is here:

COSMO.jl/docs/src/getting_started.md

Line 35 in c3a04f0

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``

Either by adding a note under the constraints. Or by adding a sentence to the PsdConeTriangle that one can use complex numbers as well (and linking to the test)

[araujoms]

Done. I added another line to the list of cones in Getting Started, and gave more detail in the docstring.

[migarstka]

Thanks a lot

[araujoms]

My pleasure. Thanks for merging.