tests and examples for the complex relative entropy cone

examples/entanglementassisted/JuMP.jl

@@ -19,8 +19,8 @@ struct EntanglementAssisted{T <: Real} <: ExampleInstanceJuMP{T}
     ne::Int
 end
 
-function build(inst::EntanglementAssisted{T}) where {T <: Float64}
-    gamma = 0.2
+function build(inst::EntanglementAssisted{T}) where {T <: Real}
+    gamma = T(1) / 5
     ampl_damp = [
         1 0
         0 sqrt(gamma)
@@ -33,46 +33,57 @@ function build(inst::EntanglementAssisted{T}) where {T <: Float64}
     ne = inst.ne
     @assert nb * ne == ampl_dim
     rt2 = sqrt(T(2))
-    sa = Cones.svec_length(ampl_dim)
-    sb = Cones.svec_length(nb)
+    sbe = Cones.svec_length(Complex, ampl_dim)
+    sb = Cones.svec_length(Complex, nb)
 
-    model = JuMP.Model()
+    model = JuMP.GenericModel{T}()
     JuMP.@variables(model, begin
-        ρ[1:na, 1:na], PSD
-        cond_epi
-        qe_epi
+        ρA[1:na, 1:na], Hermitian
+        conditional
+        von_neumann
     end)
 
-    Q1 = Symmetric(ampl_damp * ρ * ampl_damp')
-    Q2 = zeros(JuMP.AffExpr, nb * ne, nb * ne)
-    kron!(Q2, I(nb), partial_trace(Q1, 1, [nb, ne]))
-    Q3 = partial_trace(Q1, 2, [nb, ne])
-    Q1_vec = Cones.smat_to_svec!(zeros(JuMP.AffExpr, sa), Q1, rt2)
-    Q2_vec = Cones.smat_to_svec!(zeros(JuMP.AffExpr, sa), Q2, rt2)
-    Q3_vec = Cones.smat_to_svec!(zeros(JuMP.AffExpr, sb), Q3, rt2)
-    RE_cone = Hypatia.EpiTrRelEntropyTriCone{T}(1 + 2 * sa)
-    E_cone =
-        Hypatia.EpiPerSepSpectralCone{T}(Cones.NegEntropySSF(), Cones.MatrixCSqr{T, T}, nb)
+    ρBE = Hermitian(ampl_damp * ρA * ampl_damp')
+    IρE = Matrix{JuMP.GenericAffExpr{Complex{T}, JuMP.GenericVariableRef{T}}}(
+        undef,
+        nb * ne,
+        nb * ne,
+    )
+    kron!(IρE, I(nb), partial_trace(ρBE, 1, [nb, ne]))
+    ρB = partial_trace(ρBE, 2, [nb, ne])
 
+    ρBE_vec = Vector{JuMP.GenericAffExpr{T, JuMP.GenericVariableRef{T}}}(undef, sbe)
+    IρE_vec = Vector{JuMP.GenericAffExpr{T, JuMP.GenericVariableRef{T}}}(undef, sbe)
+    ρB_vec = Vector{JuMP.GenericAffExpr{T, JuMP.GenericVariableRef{T}}}(undef, sb)
+
+    Cones._smat_to_svec_complex!(ρBE_vec, ρBE, rt2)
+    Cones._smat_to_svec_complex!(IρE_vec, IρE, rt2)
+    Cones._smat_to_svec_complex!(ρB_vec, ρB, rt2)
+    RE_cone = Hypatia.EpiTrRelEntropyTriCone{T, Complex{T}}(1 + 2 * sbe)
+    E_cone = Hypatia.EpiPerSepSpectralCone{T}(
+        Cones.NegEntropySSF(),
+        Cones.MatrixCSqr{T, Complex{T}},
+        nb,
+    )
     JuMP.@constraints(model, begin
-        tr(ρ) == 1
-        vcat(cond_epi, Q1_vec, Q2_vec) in RE_cone
-        vcat(qe_epi, 1, Q3_vec) in E_cone
+        tr(ρA) == 1
+        vcat(conditional, IρE_vec, ρBE_vec) in RE_cone
+        vcat(von_neumann, 1, ρB_vec) in E_cone
     end)
 
-    JuMP.@objective(model, Max, (cond_epi + qe_epi) / -log(T(2)))
+    JuMP.@objective(model, Max, (conditional + von_neumann) / -log(T(2)))
 
     return model
 end
 
 # partial trace of Q over system i given subsystem dimensions subs
-function partial_trace(Q::Symmetric, i::Int, subs::Vector{Int})
+function partial_trace(Q::AbstractMatrix, i::Int, subs::Vector{Int})
     @assert 1 <= i <= length(subs)
     @assert size(Q, 1) == prod(subs)
     return sum(partial_trace_j(j, Q, i, subs) for j in 1:subs[i])
 end
 
-function partial_trace_j(j::Int, Q::Symmetric, i::Int, subs::Vector{Int})
+function partial_trace_j(j::Int, Q::AbstractMatrix, i::Int, subs::Vector{Int})
     X1 = sparse(I, 1, 1)
     X2 = copy(X1)
     for (k, dim) in enumerate(subs)

examples/nearestcorrelation/JuMP.jl

@@ -14,25 +14,25 @@ struct NearestCorrelationJuMP{T <: Real} <: ExampleInstanceJuMP{T}
     side::Int
 end
 
-function build(inst::NearestCorrelationJuMP{T}) where {T <: Float64}
+function build(inst::NearestCorrelationJuMP{T}) where {T <: Real}
     side = inst.side
     M = randn(T, side, side)
     M = M * M'
     vec_dim = Cones.svec_length(side)
-    m_vec = zeros(T, vec_dim)
+    m_vec = Vector{T}(undef, vec_dim)
     Cones.smat_to_svec!(m_vec, M, sqrt(T(2)))
 
-    model = JuMP.Model()
+    model = JuMP.GenericModel{T}()
     JuMP.@variable(model, x_vec[1:vec_dim])
-    X = zeros(JuMP.AffExpr, side, side)
+    X = Matrix{JuMP.GenericAffExpr{T, JuMP.GenericVariableRef{T}}}(undef, side, side)
     Cones.svec_to_smat!(X, one(T) * x_vec, sqrt(T(2)))
     JuMP.@constraint(model, diag(X) .== 1)
 
     JuMP.@variable(model, y)
     JuMP.@objective(model, Min, y)
     JuMP.@constraint(
         model,
-        vcat(y, x_vec, m_vec) in Hypatia.EpiTrRelEntropyTriCone{T}(1 + 2 * vec_dim)
+        vcat(y, x_vec, m_vec) in Hypatia.EpiTrRelEntropyTriCone{T, T}(1 + 2 * vec_dim)
     )
 
     return model

examples/relentrentanglement/JuMP.jl

@@ -16,36 +16,41 @@ struct RelEntrEntanglementJuMP{T <: Real} <: ExampleInstanceJuMP{T}
     nb::Int
 end
 
-function build(inst::RelEntrEntanglementJuMP{T}) where {T <: Float64}
+function build(inst::RelEntrEntanglementJuMP{T}) where {T <: Real}
     (na, nb) = (inst.na, inst.nb)
     side = na * nb
-    Rho = randn(T, side, side)
+    Rho = randn(Complex{T}, side, side)
     Rho = Rho * Rho'
-    Rho = Symmetric(Rho / tr(Rho))
-    vec_dim = Cones.svec_length(side)
-    rho_vec = zeros(T, vec_dim)
+    Rho = Hermitian(Rho / tr(Rho))
+    vec_dim = Cones.svec_length(Complex, side)
+    rho_vec = Vector{T}(undef, vec_dim)
     Cones.smat_to_svec!(rho_vec, Rho, sqrt(T(2)))
 
-    model = JuMP.Model()
+    model = JuMP.GenericModel{T}()
     JuMP.@variable(model, tau_vec[1:vec_dim])
-    Tau = zeros(JuMP.AffExpr, side, side)
-    Cones.svec_to_smat!(Tau, one(T) * tau_vec, sqrt(T(2)))
+    Tau = Matrix{JuMP.GenericAffExpr{Complex{T}, JuMP.GenericVariableRef{T}}}(
+        undef,
+        side,
+        side,
+    )
+    Cones._svec_to_smat_complex!(Tau, one(T) * tau_vec, sqrt(T(2)))
     JuMP.@constraint(model, tr(Tau) == 1)
 
     JuMP.@variable(model, y)
     JuMP.@objective(model, Min, y / log(T(2)))
     JuMP.@constraint(
         model,
-        vcat(y, tau_vec, rho_vec) in Hypatia.EpiTrRelEntropyTriCone{T}(1 + 2 * vec_dim)
+        vcat(y, tau_vec, rho_vec) in
+        Hypatia.EpiTrRelEntropyTriCone{T, Complex{T}}(1 + 2 * vec_dim)
     )
-    pt = partial_transpose(Symmetric(Tau), 2, [na, nb])
-    JuMP.@constraint(model, Symmetric(pt) in JuMP.PSDCone())
+    pt = partial_transpose(Hermitian(Tau), 2, [na, nb])
+    JuMP.@constraint(model, Hermitian(pt) in JuMP.HermitianPSDCone())
 
     return model
 end
 
 # partial transpose of Q over system i given subsystem dimensions subs
-function partial_transpose(Q::Symmetric, i::Int, subs::Vector{Int})
+function partial_transpose(Q::AbstractMatrix, i::Int, subs::Vector{Int})
     @assert 1 <= i <= length(subs)
     @assert size(Q, 1) == prod(subs)
     n = length(subs)

src/Cones/arrayutilities.jl

@@ -202,11 +202,13 @@ $(SIGNATURES)
 Copy a complex Hermitian matrix upper triangle to a svec-scaled real vector
 in-place.
 """
-function smat_to_svec!(
+smat_to_svec!(
     vec::AbstractVector{T},
     mat::AbstractMatrix{Complex{T}},
     rt2::Real,
-) where {T}
+) where {T} = _smat_to_svec_complex!(vec, mat, rt2)
+
+function _smat_to_svec_complex!(vec, mat, rt2)
     k = 1
     m = size(mat, 1)
     @assert m == size(mat, 2)
@@ -258,6 +260,10 @@ function svec_to_smat!(
     vec::AbstractVector{T},
     rt2::Real,
 ) where {T}
+    return _svec_to_smat_complex!(mat, vec, rt2)
+end
+
+function _svec_to_smat_complex!(mat, vec, rt2)
     k = 1
     m = size(mat, 1)
     @assert m == size(mat, 2)
@@ -266,7 +272,7 @@ function svec_to_smat!(
             mat[i, j] = vec[k]
             k += 1
         else
-            mat[i, j] = Complex(vec[k], -vec[k + 1]) / rt2
+            mat[i, j] = (vec[k] - im * vec[k + 1]) / rt2
             k += 2
         end
     end
@@ -375,7 +381,6 @@ function eig_dot_kron!(
 ) where {T <: Real, R <: RealOrComplex{T}}
     @assert issymmetric(inner) # must be symmetric (wrapper is less efficient)
     rt2i = inv(rt2)
-    d = size(inner, 1)
     copyto!(V, vecs') # allows fast column slices
     V_views = [view(V, :, i) for i in 1:size(inner, 1)]
     scals = (R <: Complex{T} ? [rt2i, rt2i * im] : [rt2i]) # real and imag parts

src/Cones/epipersepspectral/matrixcsqr.jl

@@ -492,7 +492,6 @@ function dder3(cone::EpiPerSepSpectral{<:MatrixCSqr{T}}, dir::AbstractVector{T})
     w_λi = cache.w_λi
     σ = cache.σ
     ∇h = cache.∇h
-    ∇2h = cache.∇2h
     Δ2h = cache.Δ2h
     r_X = cache.w1
     ξ_X = cache.w2

src/Cones/epirelentropy.jl

@@ -397,10 +397,14 @@
         u = 1.2023 / rtw_dim - 0.015
         v = 0.432 / rtw_dim + 1.0125
         w = -0.3057 / rtw_dim + 0.972
-    else
+    elseif w_dim <= 300
         u = 1.1513 / rtw_dim - 0.0069
         v = 0.4873 / rtw_dim + 1.0008
         w = -0.4247 / rtw_dim + 0.9961
+    else # use asymptotic expansion for the highest dimensions
+        u = 1 / rtw_dim + 0.75 / w_dim
+        v = 1 + 0.5 / rtw_dim
+        w = 1 - 0.5 / rtw_dim + 0.25 / w_dim
     end
     return [u, v, w]
 end