Thanks for checking out the issue.
I’ve since tried simplifying the optimisation to the following minimisation to isolate where the issue stems from:
\begin{aligned}
\min_{t,\rho} \;&\; t,\\
\text{for} \;&\; J := \sum_{j,k} E_{jk} \otimes \left\{ U^\dagger \left[
K\left(E_{jk} \otimes I\right)K^\dagger - E_{jk} \otimes \rho \right]
U\right\} \\
\text{subject to}\;&\;\operatorname{tr}\rho = 1,\\
&\; \rho \succeq 0, \\
&\; tI - J \succeq 0, \\
&\; tI + J \succeq 0.
\end{aligned}
Here E_{jk} is the matrix with 1 on the j-th row, k-th column, and zero everywhere else; I is the identity matrix. This explicitly minimises the spectral norm of J. The spectral norm is unitarily invariant, so it should give the same result no matter what U is.
My findings are
- JuMP gives consistent results using the same optimiser and tolerances, so it’s a problem with Convex.jl
- The issue only occurs when both U and K are complex. It doesn’t happen when either of them is purely real or purely imaginary
Here’s my code:
using Convex, SCS, LinearAlgebra, JuMP
function opt_Convex(U,K;d₁,d₂)
ρ = Semidefinite(d₂)
t = Variable()
J = sum(kron(
((1:d₁).==j)*((1:d₁).==k)',
U'*(
K*kron(((1:d₁).==j)*((1:d₁).==k)', I(d₂))*K' -
kron(((1:d₁).==j)*((1:d₁).==k)', ρ)
)*U
) for j ∈ 1:d₁, k ∈ 1:d₁)
constraints = [
tr(ρ) == 1,
t*I(d₁^2*d₂)-J ⪰ 0,
t*I(d₁^2*d₂)+J ⪰ 0,
]
optimizer = SCS.Optimizer()
Convex.MOI.set(optimizer, Convex.MOI.RawOptimizerAttribute("eps_abs"), 1E-9)
Convex.MOI.set(optimizer, Convex.MOI.RawOptimizerAttribute("eps_rel"), 1E-12)
problem = minimize(t, constraints)
solve!(problem, optimizer)
return problem.optval, opnorm(evaluate(J), 2)
end
function opt_JuMP(U,K;d₁,d₂)
model = Model(SCS.Optimizer)
@variable(model, ρ[1:d₂,1:d₂] ∈ HermitianPSDCone())
@variable(model, t)
J = sum(kron(
((1:d₁).==j)*((1:d₁).==k)',
U'*(
K*kron(((1:d₁).==j)*((1:d₁).==k)', I(d₂))*K' -
kron(((1:d₁).==j)*((1:d₁).==k)', ρ)
)*U
) for j ∈ 1:d₁, k ∈ 1:d₁)
@constraint(model, tr(ρ)==1)
@constraint(model, Hermitian(t*I(d₁^2*d₂)-J) ∈ HermitianPSDCone())
@constraint(model, Hermitian(t*I(d₁^2*d₂)+J) ∈ HermitianPSDCone())
MOI.set(model, MOI.RawOptimizerAttribute("eps_abs"), 1E-9)
MOI.set(model, MOI.RawOptimizerAttribute("eps_rel"), 1E-12)
@objective(model, Min, t)
optimize!(model)
return objective_value(model), opnorm(value.(J), 2)
end
d₁,d₂ = 2,2;
# K[1] is purely real, K[2] is purely imaginary, K[3] is complex
K = Any[reshape(1:(d₁*d₂)^2,d₁*d₂,d₁*d₂) + im*reshape((d₁*d₂)^2+1:2*(d₁*d₂)^2,d₁*d₂,d₁*d₂) for _ ∈ 1:3];
K[3] *= sqrt(d₁/tr(K[3]'K[3]));
K[1] = real.(K[3]);
K[1] *= sqrt(d₁/tr(K[1]'K[1]));
K[2] = K[1]*im;
# The unitaries are identities with complex phases
# U[1] is purely real, U[2] is purely imaginary, U[3] is complex
U = Any[I(d₁*d₂) for _ ∈ 1:3];
U[2] = U[1]*im;
U[3] = U[1]*exp(im*π/4);
# Here you’ll find that Δarr[:,:,:,2] (JuMP) is consistent
# between the optimal and evaluated values, while there is
# a discrepancy in Δarr[:,:,:,1] (Convex.jl), but only at
# Δarr[3,3,:,1] when both U and K are complex
Δarr = zeros(3,3,2,2);
for j ∈ 1:3, k ∈ 1:3
Δarr[j,k,:,1] .= opt_Convex(U[j],K[k];d₁,d₂);
Δarr[j,k,:,2] .= opt_JuMP(U[j],K[k];d₁,d₂);
end
With the output
julia> Δarr
3×3×2×2 Array{Float64, 4}:
[:, :, 1, 1] =
1.52194 1.52194 1.44642
1.52194 1.52194 1.44642
1.52194 1.52194 1.44032
[:, :, 2, 1] =
1.52194 1.52194 1.44642
1.52194 1.52194 1.44642
1.52194 1.52194 1.44759
[:, :, 1, 2] =
1.52194 1.52194 1.44633
1.52194 1.52194 1.44633
1.52194 1.52194 1.44633
[:, :, 2, 2] =
1.52194 1.52194 1.44633
1.52194 1.52194 1.44633
1.52194 1.52194 1.44633
Each column should be consistent as they are optimised over the same K. The only discrepancy happens for Convex.jl ([:,:,1,1] and [:,:,2,1]) with both U and K complex.
Maybe it has to do with some routine that separates the real and imaginary parts when Convex.jl preprocesses the variables before sending it to the optimiser?
)
. Thanks anyway!