From d37c005e914ad89d29b442b74e704eb736837a13 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mateus=20Ara=C3=BAjo?= Date: Sat, 17 Feb 2024 00:03:50 +0100 Subject: [PATCH] Support complex/hermitian relative entropy cone (#831) Co-authored-by: Lea Kapelevich Co-authored-by: Chris Coey --- examples/entanglementassisted/JuMP.jl | 59 ++-- examples/nearestcorrelation/JuMP.jl | 10 +- examples/relentrentanglement/JuMP.jl | 29 +- src/Cones/arrayutilities.jl | 13 +- src/Cones/epipersepspectral/matrixcsqr.jl | 1 - src/Cones/epirelentropy.jl | 6 +- src/Cones/epitrrelentropytri.jl | 375 ++++++++++++++-------- src/MathOptInterface/cones.jl | 18 +- src/linearalgebra/dense.jl | 22 ++ test/cone.jl | 20 +- test/moicones.jl | 11 +- test/nativeinstances.jl | 96 +++--- test/runconetests.jl | 4 +- 13 files changed, 409 insertions(+), 255 deletions(-) diff --git a/examples/entanglementassisted/JuMP.jl b/examples/entanglementassisted/JuMP.jl index c153dd839..9e7ff484a 100644 --- a/examples/entanglementassisted/JuMP.jl +++ b/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) diff --git a/examples/nearestcorrelation/JuMP.jl b/examples/nearestcorrelation/JuMP.jl index 826bbac7b..04df3f3f5 100644 --- a/examples/nearestcorrelation/JuMP.jl +++ b/examples/nearestcorrelation/JuMP.jl @@ -14,17 +14,17 @@ 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) @@ -32,7 +32,7 @@ function build(inst::NearestCorrelationJuMP{T}) where {T <: Float64} 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 diff --git a/examples/relentrentanglement/JuMP.jl b/examples/relentrentanglement/JuMP.jl index b22d55475..60d11eea0 100644 --- a/examples/relentrentanglement/JuMP.jl +++ b/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) diff --git a/src/Cones/arrayutilities.jl b/src/Cones/arrayutilities.jl index e529c6af1..036951ddf 100644 --- a/src/Cones/arrayutilities.jl +++ b/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 diff --git a/src/Cones/epipersepspectral/matrixcsqr.jl b/src/Cones/epipersepspectral/matrixcsqr.jl index fac73e6e5..23b4b9780 100644 --- a/src/Cones/epipersepspectral/matrixcsqr.jl +++ b/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 diff --git a/src/Cones/epirelentropy.jl b/src/Cones/epirelentropy.jl index abcdcac4e..47f064a90 100644 --- a/src/Cones/epirelentropy.jl +++ b/src/Cones/epirelentropy.jl @@ -397,10 +397,14 @@ function get_central_ray_epirelentropy(w_dim::Int) 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 diff --git a/src/Cones/epitrrelentropytri.jl b/src/Cones/epitrrelentropytri.jl index 0284c9e48..94162e522 100644 --- a/src/Cones/epitrrelentropytri.jl +++ b/src/Cones/epitrrelentropytri.jl @@ -5,8 +5,6 @@ This Julia package Hypatia.jl is released under the MIT license; see LICENSE file in the root directory or at https://github.com/jump-dev/Hypatia.jl =# -# TODO add hess_prod, generalize for complex - """ $(TYPEDEF) @@ -14,10 +12,11 @@ Epigraph of matrix relative entropy cone of dimension `dim` in svec format. $(FUNCTIONNAME){T}(dim::Int, use_dual::Bool = false) """ -mutable struct EpiTrRelEntropyTri{T <: Real} <: Cone{T} +mutable struct EpiTrRelEntropyTri{T <: Real, R <: RealOrComplex{T}} <: Cone{T} use_dual_barrier::Bool dim::Int d::Int + is_complex::Bool point::Vector{T} dual_point::Vector{T} @@ -28,6 +27,7 @@ mutable struct EpiTrRelEntropyTri{T <: Real} <: Cone{T} feas_updated::Bool grad_updated::Bool hess_updated::Bool + hess_aux_updated::Bool inv_hess_updated::Bool hess_fact_updated::Bool dder3_aux_updated::Bool @@ -41,18 +41,18 @@ mutable struct EpiTrRelEntropyTri{T <: Real} <: Cone{T} vw_dim::Int V_idxs::UnitRange{Int} W_idxs::UnitRange{Int} - V::Matrix{T} - W::Matrix{T} - V_fact::Eigen{T} - W_fact::Eigen{T} - Vi::Matrix{T} - Wi::Matrix{T} - W_sim::Matrix{T} + V::Matrix{R} + W::Matrix{R} + V_fact::Eigen{R} + W_fact::Eigen{R} + Vi::Matrix{R} + Wi::Matrix{R} + W_sim::Matrix{R} V_λ_log::Vector{T} W_λ_log::Vector{T} - V_log::Matrix{T} - W_log::Matrix{T} - WV_log::Matrix{T} + V_log::Matrix{R} + W_log::Matrix{R} + WV_log::Matrix{R} z::T Δ2_V::Matrix{T} Δ2_W::Matrix{T} @@ -63,20 +63,25 @@ mutable struct EpiTrRelEntropyTri{T <: Real} <: Cone{T} d2zdV2::Matrix{T} d2zdW2::Matrix{T} d2zdVW::Matrix{T} - mat::Matrix{T} - mat2::Matrix{T} - mat3::Matrix{T} - ten3d::Array{T, 3} - matd2::Matrix{T} - - function EpiTrRelEntropyTri{T}(dim::Int; use_dual::Bool = false) where {T <: Real} + tempvec::Vector{T} + mat::Matrix{R} + mat2::Matrix{R} + mat3::Matrix{R} + mat4::Matrix{R} + Wsim_Δ3::Array{R, 3} + + function EpiTrRelEntropyTri{T, R}( + dim::Int; + use_dual::Bool = false, + ) where {T <: Real, R <: RealOrComplex{T}} @assert dim > 2 @assert isodd(dim) - cone = new{T}() + cone = new{T, R}() cone.use_dual_barrier = use_dual cone.dim = dim + cone.is_complex = (R <: Complex) cone.vw_dim = div(dim - 1, 2) - cone.d = svec_side(cone.vw_dim) + cone.d = svec_side(R, cone.vw_dim) return cone end end @@ -86,27 +91,30 @@ function reset_data(cone::EpiTrRelEntropyTri) cone.feas_updated = cone.grad_updated = cone.hess_updated = - cone.inv_hess_updated = - cone.hess_fact_updated = cone.dder3_aux_updated = false + cone.hess_aux_updated = + cone.inv_hess_updated = + cone.hess_fact_updated = cone.dder3_aux_updated = false ) end -function setup_extra_data!(cone::EpiTrRelEntropyTri{T}) where {T <: Real} +function setup_extra_data!( + cone::EpiTrRelEntropyTri{T, R}, +) where {T <: Real, R <: RealOrComplex{T}} vw_dim = cone.vw_dim d = cone.d cone.rt2 = sqrt(T(2)) cone.V_idxs = 2:(vw_dim + 1) cone.W_idxs = (vw_dim + 2):(cone.dim) - cone.V = zeros(T, d, d) - cone.W = zeros(T, d, d) - cone.Vi = zeros(T, d, d) - cone.Wi = zeros(T, d, d) - cone.W_sim = zeros(T, d, d) + cone.V = zeros(R, d, d) + cone.W = zeros(R, d, d) + cone.Vi = zeros(R, d, d) + cone.Wi = zeros(R, d, d) + cone.W_sim = zeros(R, d, d) cone.V_λ_log = zeros(T, d) cone.W_λ_log = zeros(T, d) - cone.V_log = zeros(T, d, d) - cone.W_log = zeros(T, d, d) - cone.WV_log = zeros(T, d, d) + cone.V_log = zeros(R, d, d) + cone.W_log = zeros(R, d, d) + cone.WV_log = zeros(R, d, d) cone.Δ2_V = zeros(T, d, d) cone.Δ2_W = zeros(T, d, d) cone.Δ3_V = zeros(T, d, d, d) @@ -116,11 +124,12 @@ function setup_extra_data!(cone::EpiTrRelEntropyTri{T}) where {T <: Real} cone.d2zdV2 = zeros(T, vw_dim, vw_dim) cone.d2zdW2 = zeros(T, vw_dim, vw_dim) cone.d2zdVW = zeros(T, vw_dim, vw_dim) - cone.mat = zeros(T, d, d) - cone.mat2 = zeros(T, d, d) - cone.mat3 = zeros(T, d, d) - cone.ten3d = zeros(T, d, d, d) - cone.matd2 = zeros(T, d^2, d^2) + cone.tempvec = zeros(T, vw_dim) + cone.mat = zeros(R, d, d) + cone.mat2 = zeros(R, d, d) + cone.mat3 = zeros(R, d, d) + cone.mat4 = zeros(R, d, d) + cone.Wsim_Δ3 = zeros(R, d, d, d) return end @@ -128,8 +137,9 @@ get_nu(cone::EpiTrRelEntropyTri) = 2 * cone.d + 1 function set_initial_point!( arr::AbstractVector{T}, - cone::EpiTrRelEntropyTri{T}, -) where {T <: Real} + cone::EpiTrRelEntropyTri{T, R}, +) where {T <: Real, R <: RealOrComplex{T}} + incr = (cone.is_complex ? 2 : 1) arr .= 0 # at the initial point V and W are diagonal, equivalent to epirelentropy (arr[1], v, w) = get_central_ray_epirelentropy(cone.d) @@ -137,12 +147,14 @@ function set_initial_point!( for i in 1:(cone.d) arr[1 + k] = v arr[cone.vw_dim + 1 + k] = w - k += i + 1 + k += incr * i + 1 end return arr end -function update_feas(cone::EpiTrRelEntropyTri{T}) where {T <: Real} +function update_feas( + cone::EpiTrRelEntropyTri{T, R}, +) where {T <: Real, R <: RealOrComplex{T}} @assert !cone.feas_updated point = cone.point @@ -150,6 +162,7 @@ function update_feas(cone::EpiTrRelEntropyTri{T}) where {T <: Real} for (X, idxs) in zip((cone.V, cone.W), (cone.V_idxs, cone.W_idxs)) @views svec_to_smat!(X, point[idxs], cone.rt2) end + VH = Hermitian(cone.V, :U) WH = Hermitian(cone.W, :U) if isposdef(VH) && isposdef(WH) @@ -176,7 +189,9 @@ function update_feas(cone::EpiTrRelEntropyTri{T}) where {T <: Real} return cone.is_feas end -function update_grad(cone::EpiTrRelEntropyTri{T}) where {T <: Real} +function update_grad( + cone::EpiTrRelEntropyTri{T, R}, +) where {T <: Real, R <: RealOrComplex{T}} @assert cone.is_feas rt2 = cone.rt2 zi = inv(cone.z) @@ -206,30 +221,43 @@ function update_grad(cone::EpiTrRelEntropyTri{T}) where {T <: Real} Δ2!(Δ2_V, V_λ, cone.V_λ_log) - spectral_outer!(W_sim, V_vecs', Symmetric(cone.W, :U), mat) + spectral_outer!(W_sim, V_vecs', Hermitian(cone.W, :U), mat) @. mat = W_sim * Δ2_V - spectral_outer!(mat2, V_vecs, Symmetric(mat, :U), mat3) + spectral_outer!(mat2, V_vecs, Hermitian(mat, :U), mat3) + @. mat2 *= zi @views smat_to_svec!(cone.dzdV, mat2, rt2) - axpby!(-1, Vi, -zi, mat2) + axpby!(-1, Vi, -1, mat2) @views smat_to_svec!(g[cone.V_idxs], mat2, rt2) cone.grad_updated = true return cone.grad end -function update_hess(cone::EpiTrRelEntropyTri{T}) where {T <: Real} +function update_hess_aux(cone::EpiTrRelEntropyTri) @assert cone.grad_updated + Δ2!(cone.Δ2_W, cone.W_fact.values, cone.W_λ_log) + Δ3!(cone.Δ3_V, cone.Δ2_V, cone.V_fact.values) + @inbounds for k in 1:(cone.d) + @. @views cone.Wsim_Δ3[:, :, k] = cone.W_sim * cone.Δ3_V[:, :, k] + end + + cone.hess_aux_updated = true + return cone.hess_aux_updated +end + +function update_hess(cone::EpiTrRelEntropyTri) + cone.hess_aux_updated || update_hess_aux(cone) isdefined(cone, :hess) || alloc_hess!(cone) rt2 = cone.rt2 V_idxs = cone.V_idxs W_idxs = cone.W_idxs zi = inv(cone.z) - (V_λ, V_vecs) = cone.V_fact - (W_λ, W_vecs) = cone.W_fact + V_vecs = cone.V_fact.vectors + W_vecs = cone.W_fact.vectors Δ2_V = cone.Δ2_V Δ2_W = cone.Δ2_W - Δ3_V = cone.Δ3_V + dzdV = cone.dzdV dzdW = cone.dzdW d2zdV2 = cone.d2zdV2 d2zdW2 = cone.d2zdW2 @@ -240,31 +268,26 @@ function update_hess(cone::EpiTrRelEntropyTri{T}) where {T <: Real} H = cone.hess.data # u - @views dzdVzi = H[V_idxs, 1] - @. dzdVzi = zi * cone.dzdV - H[1, 1] = abs2(zi) - @. H[1, V_idxs] = zi * dzdVzi + @. H[1, V_idxs] = zi * dzdV @. H[1, W_idxs] = zi * dzdW # vv - Δ3!(Δ3_V, Δ2_V, V_λ) - d2zdV2!(d2zdV2, V_vecs, cone.W_sim, Δ3_V, cone.ten3d, cone.matd2, mat, mat2, rt2) + d2zdV2!(d2zdV2, V_vecs, cone.Wsim_Δ3, mat, mat2, mat3, rt2) @. d2zdV2 *= -1 @views Hvv = H[V_idxs, V_idxs] symm_kron!(Hvv, cone.Vi, rt2) - mul!(Hvv, dzdVzi, dzdVzi', true, true) + mul!(Hvv, dzdV, dzdV', true, true) @. Hvv += zi * d2zdV2 # vw eig_dot_kron!(d2zdVW, Δ2_V, V_vecs, mat, mat2, mat3, rt2) @views Hvw = H[V_idxs, W_idxs] @. Hvw = -zi * d2zdVW - mul!(Hvw, dzdVzi, dzdW', true, true) + mul!(Hvw, dzdV, dzdW', true, true) # ww - Δ2!(Δ2_W, W_λ, cone.W_λ_log) eig_dot_kron!(d2zdW2, Δ2_W, W_vecs, mat, mat2, mat3, rt2) @views Hww = H[W_idxs, W_idxs] symm_kron!(Hww, cone.Wi, rt2) @@ -275,15 +298,79 @@ function update_hess(cone::EpiTrRelEntropyTri{T}) where {T <: Real} return cone.hess end +function hess_prod!(prod::AbstractVecOrMat, arr::AbstractVecOrMat, cone::EpiTrRelEntropyTri) + cone.hess_aux_updated || update_hess_aux(cone) + rt2 = cone.rt2 + V_idxs = cone.V_idxs + W_idxs = cone.W_idxs + dzdW = cone.dzdW + dzdV = cone.dzdV + z = cone.z + tempvec = cone.tempvec + temp = cone.mat + temp2 = cone.mat2 + Varr_simV = cone.mat3 + arr_W_mat = cone.mat4 + (V_λ, V_vecs) = cone.V_fact + (W_λ, W_vecs) = cone.W_fact + + @inbounds for i in 1:size(arr, 2) + @views V_arr = arr[V_idxs, i] + @views W_arr = arr[W_idxs, i] + + const1 = arr[1, i] / z + dot(dzdV, V_arr) + dot(dzdW, W_arr) + prod[1, i] = const1 / z + + @views @. prod[W_idxs, i] = dzdW * const1 + # Hwv * a_v + arr_V_mat = svec_to_smat!(temp, V_arr, rt2) + spectral_outer!(Varr_simV, V_vecs', Hermitian(arr_V_mat, :U), temp2) + @. temp = Varr_simV * cone.Δ2_V + spectral_outer!(temp, V_vecs, Hermitian(temp, :U), temp2) + @. temp /= -z + @views prod[W_idxs, i] .+= smat_to_svec!(tempvec, temp, rt2) + # Hww * a_w + svec_to_smat!(arr_W_mat, W_arr, rt2) + Warr_simW = spectral_outer!(temp2, W_vecs', Hermitian(arr_W_mat, :U), temp) + @. temp = Warr_simW * cone.Δ2_W / z + @. Warr_simW /= W_λ' + ldiv!(Diagonal(W_λ), Warr_simW) + @. temp += Warr_simW + spectral_outer!(temp, W_vecs, Hermitian(temp, :U), temp2) + @views prod[W_idxs, i] .+= smat_to_svec!(tempvec, temp, rt2) + + @views @. prod[V_idxs, i] = dzdV * const1 + # Hvv * a_v + for k in 1:(cone.d) + @views mul!(temp2[:, k], cone.Wsim_Δ3[:, :, k], Varr_simV[:, k]) + end + @. temp = temp2 + temp2' + # destroys arr_W_mat + Warr_simV = spectral_outer!(arr_W_mat, V_vecs', Hermitian(arr_W_mat, :U), temp2) + @. temp += Warr_simV * cone.Δ2_V + @. temp /= -z + @. Varr_simV /= V_λ' + ldiv!(Diagonal(V_λ), Varr_simV) + @. temp += Varr_simV + spectral_outer!(temp, V_vecs, Hermitian(temp, :U), temp2) + @views prod[V_idxs, i] .+= smat_to_svec!(tempvec, temp, rt2) + end + + return prod +end + function update_dder3_aux(cone::EpiTrRelEntropyTri) @assert !cone.dder3_aux_updated - @assert cone.hess_updated + cone.hess_aux_updated || update_hess_aux(cone) Δ3!(cone.Δ3_W, cone.Δ2_W, cone.W_fact.values) cone.dder3_aux_updated = true return end -function dder3(cone::EpiTrRelEntropyTri{T}, dir::AbstractVector{T}) where {T <: Real} +function dder3( + cone::EpiTrRelEntropyTri{T, R}, + dir::AbstractVector{T}, +) where {T <: Real, R <: RealOrComplex{T}} cone.dder3_aux_updated || update_dder3_aux(cone) rt2 = cone.rt2 V_idxs = cone.V_idxs @@ -303,26 +390,10 @@ function dder3(cone::EpiTrRelEntropyTri{T}, dir::AbstractVector{T}) where {T <: @views v_dir = dir[V_idxs] @views w_dir = dir[W_idxs] - # TODO in-place - Vvd = Symmetric(cone.d2zdV2, :U) * v_dir - Wwd = Symmetric(cone.d2zdW2, :U) * w_dir - d2zdVW = Symmetric(cone.d2zdVW, :U) - VWwd = d2zdVW * w_dir - - const0 = zi * (u_dir + dot(v_dir, dzdV)) + dot(w_dir, dzdW) - const1 = - abs2(const0) + zi * (-dot(v_dir, VWwd) + (dot(v_dir, Vvd) + dot(w_dir, Wwd)) / 2) - - V_part_1a = const0 * (Vvd - VWwd) - W_part_1a = Wwd - d2zdVW * v_dir - - # u - ziconst1 = dder3[1] = zi * const1 - # v, w # TODO prealloc - V_dir = Symmetric(zero(mat), :U) - W_dir = Symmetric(zero(mat), :U) + V_dir = Hermitian(zero(mat), :U) + W_dir = Hermitian(zero(mat), :U) V_dir_sim = zero(mat) W_dir_sim = zero(mat) VW_dir_sim = zero(mat) @@ -339,6 +410,35 @@ function dder3(cone::EpiTrRelEntropyTri{T}, dir::AbstractVector{T}) where {T <: spectral_outer!(W_dir_sim, W_vecs', W_dir, mat) spectral_outer!(VW_dir_sim, V_vecs', W_dir, mat) + for k in 1:(cone.d) + @views mul!(mat2[:, k], cone.Wsim_Δ3[:, :, k], V_dir_sim[:, k]) + end + @. mat = mat2 + mat2' + spectral_outer!(mat, V_vecs, Hermitian(mat, :U), mat2) + Vvd = smat_to_svec!(zero(w_dir), mat, rt2) + + @. mat = W_dir_sim * cone.Δ2_W + spectral_outer!(mat, W_vecs, Hermitian(mat, :U), mat2) + Wwd = smat_to_svec!(zero(w_dir), mat, rt2) + + @. mat = VW_dir_sim * cone.Δ2_V + spectral_outer!(mat, V_vecs, Hermitian(mat, :U), mat2) + VWwd = smat_to_svec!(zero(w_dir), mat, rt2) + + @. mat = V_dir_sim * cone.Δ2_V + spectral_outer!(mat, V_vecs, Hermitian(mat, :U), mat2) + VWvd = smat_to_svec!(zero(w_dir), mat, rt2) + + const0 = zi * u_dir + dot(v_dir, dzdV) + dot(w_dir, dzdW) + const1 = + abs2(const0) + zi * (-dot(v_dir, VWwd) + (-dot(v_dir, Vvd) + dot(w_dir, Wwd)) / 2) + + V_part_1a = -const0 * (Vvd + VWwd) + W_part_1a = Wwd - VWvd + + # u + dder3[1] = zi * const1 + @inbounds @views for j in 1:(cone.d) Vds_j = V_dir_sim[:, j] Wds_j = W_dir_sim[:, j] @@ -347,39 +447,39 @@ function dder3(cone::EpiTrRelEntropyTri{T}, dir::AbstractVector{T}) where {T <: Wds_i = W_dir_sim[:, i] DΔ3_V_ij = Diagonal(Δ3_V[:, i, j]) DΔ3_W_ij = Diagonal(Δ3_W[:, i, j]) - diff_dot_V_VV[i, j] = dot(Vds_j, DΔ3_V_ij, Vds_i) - diff_dot_W_WW[i, j] = dot(Wds_j, DΔ3_W_ij, Wds_i) + diff_dot_V_VV[i, j] = dot(Vds_i, DΔ3_V_ij, Vds_j) + diff_dot_W_WW[i, j] = dot(Wds_i, DΔ3_W_ij, Wds_j) end for i in 1:(cone.d) VWds_i = VW_dir_sim[:, i] DΔ3_V_ij = Diagonal(Δ3_V[:, i, j]) - diff_dot_V_VW[i, j] = dot(Vds_j, DΔ3_V_ij, VWds_i) + diff_dot_V_VW[i, j] = dot(VWds_i, DΔ3_V_ij, Vds_j) end end # v d3WlogVdV!(d3WlogVdV, Δ3_V, V_λ, V_dir_sim, cone.W_sim, mat) svec_to_smat!(V_part_1, V_part_1a, rt2) - rdiv!(V_dir_sim, Diagonal(sqrt.(V_λ))) + @. V_dir_sim /= sqrt(V_λ)' ldiv!(Diagonal(V_λ), V_dir_sim) V_part_2 = d3WlogVdV @. V_part_2 += diff_dot_V_VW + diff_dot_V_VW' mul!(V_part_2, V_dir_sim, V_dir_sim', true, zi) - mul!(mat, Symmetric(V_part_2, :U), V_vecs') + mul!(mat, Hermitian(V_part_2, :U), V_vecs') mul!(V_part_1, V_vecs, mat, true, zi) @views dder3_V = dder3[V_idxs] smat_to_svec!(dder3_V, V_part_1, rt2) - @. dder3_V += ziconst1 * dzdV + @. dder3_V += const1 * dzdV # w svec_to_smat!(W_part_1, W_part_1a, rt2) - spectral_outer!(mat2, V_vecs, Symmetric(diff_dot_V_VV, :U), mat) + spectral_outer!(mat2, V_vecs, Hermitian(diff_dot_V_VV, :U), mat) axpby!(true, mat2, const0, W_part_1) - rdiv!(W_dir_sim, Diagonal(sqrt.(W_λ))) + @. W_dir_sim /= sqrt(W_λ)' ldiv!(Diagonal(W_λ), W_dir_sim) W_part_2 = diff_dot_W_WW mul!(W_part_2, W_dir_sim, W_dir_sim', true, -zi) - mul!(mat, Symmetric(W_part_2, :U), W_vecs') + mul!(mat, Hermitian(W_part_2, :U), W_vecs') mul!(W_part_1, W_vecs, mat, true, zi) @views dder3_W = dder3[W_idxs] smat_to_svec!(dder3_W, W_part_1, rt2) @@ -440,53 +540,48 @@ function Δ3!(Δ3::Array{T, 3}, Δ2::Matrix{T}, λ::Vector{T}) where {T <: Real} return Δ3 end +# TODO consider refactor with eig_dot_kron! function d2zdV2!( d2zdV2::Matrix{T}, - vecs::Matrix{T}, - inner::Matrix{T}, - Δ3::Array{T, 3}, - ten3d::Array{T, 3}, # temp - matd2::Matrix{T}, # temp - mat::Matrix{T}, # temp - mat2::Matrix{T}, # temp + vecs::Matrix{R}, + Wsim_Δ3::Array{R, 3}, + mat::Matrix{R}, # temp + mat2::Matrix{R}, # temp + mat3::Matrix{R}, # temp rt2::T, -) where {T <: Real} +) where {T <: Real, R <: RealOrComplex{T}} d = size(vecs, 1) - - @inbounds for i in 1:d - @views ten3d_i = ten3d[i, :, :] - @views Δ3_i = Δ3[i, :, :] - @. mat2 = inner * Δ3_i - spectral_outer!(ten3d_i, vecs, Symmetric(mat2, :U), mat) - end - @inbounds for j in 1:d, i in 1:j - @views temp_ij = matd2[block_idxs(d, i), block_idxs(d, j)] - @views D_ij = ten3d[:, i, j] - spectral_outer!(temp_ij, vecs, D_ij, mat) - end - copytri!(matd2, 'U') - + V = copyto!(mat, vecs') + V_views = [view(V, :, i) for i in 1:d] rt2i = inv(rt2) - row_idx = 1 - @inbounds for j in 1:d, i in 1:j - col_idx = 1 - ijeq = (i == j) - di = d * (i - 1) - dj = d * (j - 1) - for l in 1:d, k in 1:l - rho = (k == l ? (ijeq ? T(0.5) : rt2i) : (ijeq ? rt2i : one(T))) - dlk = d * (l - 1) + k - dkl = d * (k - 1) + l - dji = dj + i - dij = di + j - d2zdV2[row_idx, col_idx] = - rho * - (matd2[dji, dlk] + matd2[dij, dlk] + matd2[dji, dkl] + matd2[dij, dkl]) + scals = (R <: Complex{T} ? [rt2i, rt2i * im] : [rt2i]) + + col_idx = 1 + @inbounds for j in 1:d + for i in 1:(j - 1), scal in scals + mul!(mat3, V_views[j], V_views[i]', scal, false) + @. mat2 = mat3 + mat3' + for k in 1:d + @views mul!(mat3[:, k], Wsim_Δ3[:, :, k], mat2[:, k]) + end + # mat2 = vecs * (mat3 + mat3) * vecs' + @. mat2 = mat3 + mat3' + mul!(mat3, Hermitian(mat2, :U), V) + mul!(mat2, V', mat3) + @views smat_to_svec!(d2zdV2[:, col_idx], mat2, rt2) col_idx += 1 end - row_idx += 1 - end + mul!(mat2, V_views[j], V_views[j]') + for k in 1:d + @views mul!(mat3[:, k], Wsim_Δ3[:, :, k], mat2[:, k]) + end + @. mat2 = mat3 + mat3' + mul!(mat3, Hermitian(mat2, :U), V) + mul!(mat2, V', mat3) + @views smat_to_svec!(d2zdV2[:, col_idx], mat2, rt2) + col_idx += 1 + end return d2zdV2 end @@ -495,30 +590,30 @@ function VdWs_element( j::Int, k::Int, l::Int, - Vds::Matrix{T}, - Ws::Matrix{T}, -) where {T <: Real} + Vds::Matrix{R}, + Ws::Matrix{R}, +) where {T <: Real, R <: RealOrComplex{T}} @inbounds begin - a = Vds[l, j] * Ws[i, k] + Vds[l, i] * Ws[j, k] - b = Vds[k, i] * Vds[l, j] * Ws[k, l] - c = Vds[k, l] * a + b + a = Ws[i, k] * Vds[k, l] + Vds[i, k] * Ws[k, l] + b = Vds[i, k] * Vds[k, l] * Ws[l, j] + c = Vds[l, j] * a + b end return c end function d3WlogVdV!( - d3WlogVdV::Matrix{T}, + d3WlogVdV::Matrix{R}, Δ3::Array{T, 3}, λ::Vector{T}, - Vds::Matrix{T}, - Ws::Matrix{T}, - Δ4_ij::Matrix{T}, # temp -) where {T <: Real} + Vds::Matrix{R}, + Ws::Matrix{R}, + Δ4_ij::Matrix{R}, +) where {T <: Real, R <: RealOrComplex{T}} d = length(λ) @inbounds for j in 1:d, i in 1:j Δ4_ij!(Δ4_ij, i, j, Δ3, λ) - t = zero(T) + t = zero(R) for l in 1:d for k in 1:(l - 1) t += @@ -530,7 +625,7 @@ function d3WlogVdV!( d3WlogVdV[i, j] = t end - copytri!(d3WlogVdV, 'U') + copytri!(d3WlogVdV, 'U', true) return d3WlogVdV end diff --git a/src/MathOptInterface/cones.jl b/src/MathOptInterface/cones.jl index 8220cfc7b..d4c67a35a 100644 --- a/src/MathOptInterface/cones.jl +++ b/src/MathOptInterface/cones.jl @@ -586,20 +586,23 @@ See [`Cones.EpiTrRelEntropyTri`](@ref). $(TYPEDFIELDS) """ -struct EpiTrRelEntropyTriCone{T <: Real} <: MOI.AbstractVectorSet +struct EpiTrRelEntropyTriCone{T <: Real, R <: RealOrComplex{T}} <: MOI.AbstractVectorSet dim::Int use_dual::Bool end export EpiTrRelEntropyTriCone -function EpiTrRelEntropyTriCone{T}(dim::Int) where {T <: Real} - return EpiTrRelEntropyTriCone{T}(dim, false) +function EpiTrRelEntropyTriCone{T, R}(dim::Int) where {T <: Real, R <: RealOrComplex{T}} + return EpiTrRelEntropyTriCone{T, R}(dim, false) end -MOI.dimension(cone::EpiTrRelEntropyTriCone where {T <: Real}) = cone.dim +MOI.dimension(cone::EpiTrRelEntropyTriCone) = cone.dim -function cone_from_moi(::Type{T}, cone::EpiTrRelEntropyTriCone{T}) where {T <: Real} - return Cones.EpiTrRelEntropyTri{T}(cone.dim, use_dual = cone.use_dual) +function cone_from_moi( + ::Type{T}, + cone::EpiTrRelEntropyTriCone{T, R}, +) where {T <: Real, R <: RealOrComplex{T}} + return Cones.EpiTrRelEntropyTri{T, R}(cone.dim, use_dual = cone.use_dual) end """ @@ -754,7 +757,8 @@ const HypatiaCones{T <: Real} = Union{ HypoPerLogdetTriCone{T, Complex{T}}, EpiPerSepSpectralCone{T}, EpiRelEntropyCone{T}, - EpiTrRelEntropyTriCone{T}, + EpiTrRelEntropyTriCone{T, T}, + EpiTrRelEntropyTriCone{T, Complex{T}}, WSOSInterpNonnegativeCone{T, T}, WSOSInterpNonnegativeCone{T, Complex{T}}, WSOSInterpPosSemidefTriCone{T}, diff --git a/src/linearalgebra/dense.jl b/src/linearalgebra/dense.jl index 494ad3e76..f4fc10ade 100644 --- a/src/linearalgebra/dense.jl +++ b/src/linearalgebra/dense.jl @@ -119,6 +119,17 @@ function spectral_outer!( return mat end +function spectral_outer!( + mat::AbstractMatrix{R}, + vecs::Union{Matrix{R}, Adjoint{T, Matrix{R}}}, + diag::AbstractVector{T}, + temp::Matrix{R}, +) where {T <: Real, R <: RealOrComplex{T}} + mul!(temp, vecs, Diagonal(diag)) + mul!(mat, temp, vecs') + return mat +end + function spectral_outer!( mat::AbstractMatrix{T}, vecs::Union{Matrix{T}, Adjoint{T, Matrix{T}}}, @@ -130,6 +141,17 @@ function spectral_outer!( return mat end +function spectral_outer!( + mat::AbstractMatrix{R}, + vecs::Union{Matrix{R}, Adjoint{R, Matrix{R}}}, + symm::Hermitian{R}, + temp::Matrix{R}, +) where {T <: Real, R <: RealOrComplex{T}} + mul!(temp, vecs, symm) + mul!(mat, temp, vecs') + return mat +end + #= nonsymmetric square: LU =# diff --git a/test/cone.jl b/test/cone.jl index 433f1ac67..285934f06 100644 --- a/test/cone.jl +++ b/test/cone.jl @@ -766,28 +766,30 @@ end show_time_alloc(C::Type{<:Cones.EpiRelEntropy}) = show_time_alloc(C(9)) # EpiTrRelEntropyTri -function test_oracles(C::Type{<:Cones.EpiTrRelEntropyTri}) +function test_oracles(C::Type{<:Cones.EpiTrRelEntropyTri{T, R}}) where {T, R} for dW in [1, 2, 4] - test_oracles(C(1 + 2 * Cones.svec_length(dW)), init_tol = 1e-4) + test_oracles(C(1 + 2 * Cones.svec_length(R, dW)), init_tol = 1e-4) end for dW in [6, 10] - test_oracles(C(1 + 2 * Cones.svec_length(dW)), init_tol = 1e-1, init_only = true) + test_oracles(C(1 + 2 * Cones.svec_length(R, dW)), init_tol = 1e-1, init_only = true) end end -function test_barrier(C::Type{Cones.EpiTrRelEntropyTri{T}}) where {T} +function test_barrier(C::Type{Cones.EpiTrRelEntropyTri{T, R}}) where {T, R} dW = 3 - dw = Cones.svec_length(dW) + dw = Cones.svec_length(R, dW) function barrier(s) (u, v, w) = (s[1], s[1 .+ (1:dw)], s[(2 + dw):end]) - V = new_herm(v, dW, T) - W = new_herm(w, dW, T) - return -log(u - dot(W, log(W) - log(V))) - logdet_pd(V) - logdet_pd(W) + V = new_herm(v, dW, R) + W = new_herm(w, dW, R) + return real(-log(u - dot(W, log(W) - log(V)))) - logdet_pd(V) - logdet_pd(W) end return test_barrier(C(1 + 2 * dw), barrier, TFD = BigFloat) end -show_time_alloc(C::Type{<:Cones.EpiTrRelEntropyTri}) = show_time_alloc(C(13)) +function show_time_alloc(C::Type{Cones.EpiTrRelEntropyTri{T, R}}) where {T, R} + return show_time_alloc(C(Cones.svec_length(R, 4) * 2 + 1)) +end # WSOSInterpNonnegative function test_oracles(C::Type{Cones.WSOSInterpNonnegative{T, R}}) where {T, R} diff --git a/test/moicones.jl b/test/moicones.jl index c7696d450..48d6ddc8b 100644 --- a/test/moicones.jl +++ b/test/moicones.jl @@ -360,10 +360,15 @@ function test_moi_cones(T::Type{<:Real}) end @testset "EpiTrRelEntropyTriCone" begin - moi_cone = Hypatia.EpiTrRelEntropyTriCone{T}(3) + moi_cone = Hypatia.EpiTrRelEntropyTriCone{T, T}(7) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.EpiTrRelEntropyTri{T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 + @test hyp_cone isa Cones.EpiTrRelEntropyTri{T, T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 7 + + moi_cone = Hypatia.EpiTrRelEntropyTriCone{T, Complex{T}}(9) + hyp_cone = Hypatia.cone_from_moi(T, moi_cone) + @test hyp_cone isa Cones.EpiTrRelEntropyTri{T, Complex{T}} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 9 end @testset "WSOSInterpNonnegative" begin diff --git a/test/nativeinstances.jl b/test/nativeinstances.jl index 249b275b8..f04a15baf 100644 --- a/test/nativeinstances.jl +++ b/test/nativeinstances.jl @@ -562,7 +562,7 @@ function doublynonnegativetri1(T; options...) A = T[1 0 0; 0 0 1] b = ones(T, 2) G = SparseMatrixCSC(-one(T) * I, q, n) - for use_dual in [false, true] + for use_dual in (false, true) use_dual = false h = (use_dual ? T[1, 0, 1] : zeros(T, q)) cones = Cone{T}[Cones.DoublyNonnegativeTri{T}(q, use_dual = use_dual)] @@ -2194,7 +2194,7 @@ function epipersepspectral_matrix1(T; options...) tol = test_tol(T) Random.seed!(1) rt2 = sqrt(T(2)) - for d in [1, 3], is_complex in [false, true] + for d in [1, 3], is_complex in (false, true) R = (is_complex ? Complex{T} : T) dim = 2 + Cones.svec_length(R, d) W = rand(R, d, d) @@ -2223,7 +2223,7 @@ function epipersepspectral_matrix2(T; options...) tol = test_tol(T) Random.seed!(1) rt2 = sqrt(T(2)) - for d in [1, 3], is_complex in [false, true] + for d in [1, 3], is_complex in (false, true) R = (is_complex ? Complex{T} : T) dim = 2 + Cones.svec_length(R, d) c = T[1] @@ -2257,8 +2257,7 @@ end function epipersepspectral_matrix3(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) - for d in [2, 4], is_complex in [false, true] + for d in [2, 4], is_complex in (false, true) R = (is_complex ? Complex{T} : T) dim = 2 + Cones.svec_length(R, d) c = zeros(T, dim) @@ -2386,56 +2385,60 @@ function epitrrelentropytri1(T; options...) Random.seed!(1) rt2 = sqrt(T(2)) side = 4 - svec_dim = Cones.svec_length(side) - dim = 2 * svec_dim + 1 - c = T[1] - A = zeros(T, 0, 1) - b = T[] - W = rand(T, side, side) - W = W * W' - V = rand(T, side, side) - V = V * V' - h = vcat( - zero(T), - Cones.smat_to_svec!(zeros(T, svec_dim), V, rt2), - Cones.smat_to_svec!(zeros(T, svec_dim), W, rt2), - ) - G = zeros(T, dim, 1) - G[1, 1] = -1 - cones = Cone{T}[Cones.EpiTrRelEntropyTri{T}(dim)] + for is_complex in (false, true) + R = (is_complex ? Complex{T} : T) + svec_dim = Cones.svec_length(R, side) + dim = 2 * svec_dim + 1 + c = T[1] + A = zeros(T, 0, 1) + b = T[] + W = randn(R, side, side) + W = W * W' + V = randn(R, side, side) + V = V * V' + h = vcat( + zero(T), + Cones.smat_to_svec!(zeros(T, svec_dim), V, rt2), + Cones.smat_to_svec!(zeros(T, svec_dim), W, rt2), + ) + G = zeros(T, dim, 1) + G[1, 1] = -1 + cones = Cone{T}[Cones.EpiTrRelEntropyTri{T, R}(dim)] - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ dot(W, log(W) - log(V)) atol = tol rtol = tol + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ dot(W, log(W) - log(V)) atol = tol rtol = tol + end end function epitrrelentropytri2(T; options...) tol = test_tol(T) rt2 = sqrt(T(2)) side = 3 - svec_dim = Cones.svec_length(side) - dim = 2 * svec_dim + 1 - c = vcat(zeros(T, svec_dim), -ones(T, svec_dim)) - A = Matrix{T}(I, svec_dim, 2 * svec_dim) - b = zeros(T, svec_dim) - b[[sum(1:i) for i in 1:side]] .= 1 - h = vcat(T(5), zeros(T, 2 * svec_dim)) - g_svec = Cones.scale_svec!(-ones(T, svec_dim), sqrt(T(2))) - G = vcat(zeros(T, 2 * svec_dim)', Diagonal(vcat(g_svec, g_svec))) - cones = Cone{T}[Cones.EpiTrRelEntropyTri{T}(dim)] + for is_complex in (false, true) + R = (is_complex ? Complex{T} : T) + svec_dim = Cones.svec_length(R, side) + dim = 2 * svec_dim + 1 + c = vcat(zeros(T, svec_dim), -ones(T, svec_dim)) + A = Matrix{T}(I, svec_dim, 2 * svec_dim) + b = zeros(T, svec_dim) + Cones.smat_to_svec!(b, Matrix{R}(I, side, side), rt2) + h = vcat(T(5), zeros(T, 2 * svec_dim)) + G = vcat(zeros(T, 2 * svec_dim)', Diagonal(-ones(T, 2 * svec_dim))) + cones = Cone{T}[Cones.EpiTrRelEntropyTri{T, R}(dim)] - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - W = Hermitian( - Cones.svec_to_smat!(zeros(T, side, side), r.s[(svec_dim + 2):end], rt2), - :U, - ) - @test dot(W, log(W)) ≈ T(5) atol = tol rtol = tol + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + W = Hermitian( + Cones.svec_to_smat!(zeros(R, side, side), r.s[(svec_dim + 2):end], rt2), + :U, + ) + @test dot(W, log(W)) ≈ T(5) atol = tol rtol = tol + end end function epitrrelentropytri3(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) side = 3 svec_dim = Cones.svec_length(side) dim = 2 * svec_dim + 1 @@ -2444,7 +2447,7 @@ function epitrrelentropytri3(T; options...) b = [zero(T)] h = zeros(T, dim) G = Diagonal(-one(T) * I, dim) - cones = Cone{T}[Cones.EpiTrRelEntropyTri{T}(dim)] + cones = Cone{T}[Cones.EpiTrRelEntropyTri{T, T}(dim)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal @@ -2455,8 +2458,7 @@ end function epitrrelentropytri4(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) - side = 3 + side = 4 svec_dim = Cones.svec_length(side) dim = 2 * svec_dim + 1 c = vcat(zero(T), ones(T, svec_dim), zeros(T, svec_dim)) @@ -2464,7 +2466,7 @@ function epitrrelentropytri4(T; options...) b = [zero(T)] h = zeros(T, dim) G = Diagonal(-one(T) * I, dim) - cones = Cone{T}[Cones.EpiTrRelEntropyTri{T}(dim)] + cones = Cone{T}[Cones.EpiTrRelEntropyTri{T, T}(dim)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal diff --git a/test/runconetests.jl b/test/runconetests.jl index 702bdcbe7..491579d04 100644 --- a/test/runconetests.jl +++ b/test/runconetests.jl @@ -45,7 +45,8 @@ function cone_types(T::Type{<:Real}) Cones.EpiPerSepSpectral{Cones.MatrixCSqr{T, T}, T}, Cones.EpiPerSepSpectral{Cones.MatrixCSqr{T, Complex{T}}, T}, Cones.EpiRelEntropy{T}, - Cones.EpiTrRelEntropyTri{T}, + Cones.EpiTrRelEntropyTri{T, T}, + Cones.EpiTrRelEntropyTri{T, Complex{T}}, Cones.WSOSInterpNonnegative{T, T}, Cones.WSOSInterpNonnegative{T, Complex{T}}, Cones.WSOSInterpPosSemidefTri{T}, @@ -117,5 +118,4 @@ sep_spectral_funs = [ # end # println() # end - end;