How to Minimize Operator/Spectral Norm of a Density Matrix with JuMP + Hypatia

Specific Domains Optimization (Mathematical)
,
1

I am coming from MATLAB/Python but Julia is becoming more attractive, especially for higher precision. Currently I am just trying out some toy problems to familiarize myself with Julia, JuMP, and associated solvers.

I am trying to minimize the operator norm (spectral/inf-schatten norm) of a density matrix \rho, subject to the constraint that some (fixed) test POVM element F_i is valid: tr(F_i\rho)\in[0, 1].

Quantum info jargon aside, the important part here is that \rho has the properties of being a Hermitian and positive semi-definite 2x2 matrix with trace 1. It has complex entries. F_i can be some fixed 2x2 matrix that constrains \rho. Here is my first attempt at doing so with with JuMP.

using LinearAlgebra
using JuMP
using Hypatia

model = Model(() -> Hypatia.Optimizer(verbose = false))

#Introduces rho (H) as a 2d PSD Hermitian
@variable(model, H[1:2, 1:2] in HermitianPSDCone())

#Unit trace constraint
@constraint(model, LinearAlgebra.tr(H) == 1)

#op norm
@variable(model, t)
@constraint(model, [t; vec(H)] in MOI.NormSpectralCone(2, 2))

#A test POVM element
F = [1+0*im 0; 0 0]

#Probability must be a probability
@constraint(model, real(LinearAlgebra.tr(F * H)) in MOI.Interval(0, 1))

#Minimize the op norm
@objective(model, Min, t)

println("Starting solver.")

optimize!(model)
println(termination_status(model))
println(objective_value(model))
println(value.(H))

However, I obtain the error:

ERROR: LoadError: MethodError: no method matching promote_operation(::typeof(vcat), ::Type{Float64}, ::Type{MathOptInterface.ScalarAffineFunction{ComplexF64}}, ::Type{Float64})
The function `promote_operation` exists, but no method is defined for this combination of argument types.
...

Can JuMP handle complex numbers in this case?

1 Like
2

Hi @Menixm, welcome to the forum :smile:

That’s an unfortunate error. I’ll see if we can improve it.

The underlying issue is that Hypatia hasn’t implemented support for MOI.ScalarAffineFunction{Complex{Float64}} in the MOI.NormSpectralCone set. The error happens because we try to reformulate the problem into something that Hypatia does support, but this fails.

I think you can skip MOI.NormSpectralCone and use the Hypatia cones directly. Something like:

julia> using LinearAlgebra

julia> using JuMP

julia> using Hypatia

julia> model = Model(Hypatia.Optimizer)
A JuMP Model
├ solver: Hypatia
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none

julia> @variable(model, H[1:2, 1:2] in HermitianPSDCone())
2×2 Hermitian{GenericAffExpr{ComplexF64, VariableRef}, Matrix{GenericAffExpr{ComplexF64, VariableRef}}}:
 real(H[1,1])                    real(H[1,2]) + imag(H[1,2]) im
 real(H[1,2]) - imag(H[1,2]) im  real(H[2,2])

julia> @constraint(model, LinearAlgebra.tr(H) == 1)
real(H[1,1]) + real(H[2,2]) = 1

julia> @variable(model, t)
t

julia> @constraint(
           model,
           [t; vec(real.(H)); vec(imag.(H))] in Hypatia.EpiNormSpectralCone{Float64,ComplexF64}(2, 2),
       )
[t, real(H[1,1]), real(H[1,2]), real(H[1,2]), real(H[2,2]), 0, -imag(H[1,2]), imag(H[1,2]), 0] ∈ EpiNormSpectralCone{Float64, ComplexF64}(2, 2, false)

julia> F = [1+0*im 0; 0 0]
2×2 Matrix{Complex{Int64}}:
 1+0im  0+0im
 0+0im  0+0im

julia> @constraint(model, 0 <= real(LinearAlgebra.tr(F * H)) <= 1)
real(H[1,1]) ∈ [0, 1]

julia> @objective(model, Min, t)
t

julia> optimize!(model)

 iter        p_obj        d_obj |  abs_gap    x_feas    z_feas |      tau       kap        mu | dir_res     prox  step     alpha
    0   1.7321e+00  -1.8409e+00 | 7.00e+00  3.66e-01  3.86e-01 | 1.00e+00  1.00e+00  1.00e+00 |
    1   6.0358e-01  -2.1827e-01 | 2.02e+00  8.36e-02  8.83e-02 | 1.31e+00  2.93e-01  3.00e-01 | 2.2e-16  7.8e-01  co-a  7.00e-01
    2   7.2189e-01   6.3339e-01 | 2.37e-01  7.66e-03  8.08e-03 | 1.43e+00  1.03e-02  3.14e-02 | 1.8e-15  9.1e-01  co-a  9.00e-01
    3   7.0848e-01   6.9513e-01 | 3.05e-02  1.25e-03  1.32e-03 | 1.32e+00  2.98e-03  4.31e-03 | 4.7e-16  5.5e-01  co-a  8.50e-01
    4   7.0728e-01   7.0587e-01 | 2.89e-03  1.32e-04  1.39e-04 | 1.25e+00  2.96e-04  4.07e-04 | 3.0e-15  6.5e-01  co-a  9.00e-01
    5   7.0713e-01   7.0698e-01 | 2.72e-04  1.41e-05  1.49e-05 | 1.17e+00  2.83e-05  3.81e-05 | 4.6e-14  9.4e-01  co-a  9.00e-01
    6   7.0711e-01   7.0708e-01 | 3.67e-05  2.32e-06  2.45e-06 | 1.06e+00  4.58e-06  5.20e-06 | 3.5e-13  4.6e-01  co-a  8.50e-01
    7   7.0711e-01   7.0710e-01 | 3.50e-06  2.44e-07  2.58e-07 | 1.01e+00  4.59e-07  4.95e-07 | 2.2e-12  4.4e-01  co-a  9.00e-01
    8   7.0711e-01   7.0711e-01 | 3.33e-07  2.56e-08  2.70e-08 | 9.65e-01  4.62e-08  4.72e-08 | 1.6e-11  3.9e-01  co-a  9.00e-01
    9   7.0711e-01   7.0711e-01 | 3.19e-08  2.67e-09  2.82e-09 | 9.25e-01  4.70e-09  4.53e-09 | 5.5e-11  2.9e-01  co-a  9.00e-01
   10   7.0711e-01   7.0711e-01 | 1.23e-08  1.08e-09  1.15e-09 | 9.04e-01  2.00e-09  1.77e-09 | 5.6e-11  1.7e-01  co-a  6.00e-01
optimal solution found; terminating

status is Optimal after 10 iterations and 0.007 seconds


julia> termination_status(model)
OPTIMAL::TerminationStatusCode = 1

julia> objective_value(model)
0.7071067808809698

julia> value.(H)
2×2 Matrix{ComplexF64}:
 0.5+0.0im  0.0+0.0im
 0.0-0.0im  0.5+0.0im

I don’t know if I got the ordering of the real and imaginary variables correct in [t; vec(real.(H)); vec(imag.(H))]. Perhaps @chriscoey or @lkapelevich can confirm.

1 Like
3

It might actually be this:

julia> using LinearAlgebra

julia> using JuMP

julia> using Hypatia

julia> mat_to_vec(h::AbstractMatrix) = mat_to_vec(real.(h), imag.(h))
mat_to_vec (generic function with 3 methods)

julia> function mat_to_vec(HR::AbstractMatrix{T}, HI::AbstractMatrix{T}) where {T}
           y = Vector{T}(undef, 2 * length(HR))
           for (k, h_k) in enumerate(zip(HR, HI))
               y[2k-1], y[2k] = h_k
           end
           return y
       end
mat_to_vec (generic function with 3 methods)

julia> begin
           model = Model(Hypatia.Optimizer)
           @variable(model, H[1:2, 1:2] in HermitianPSDCone())
           @constraint(model, LinearAlgebra.tr(H) == 1)
           @variable(model, t)
           @constraint(
               model,
               [t; mat_to_vec(H)] in Hypatia.EpiNormSpectralCone{Float64,ComplexF64}(2, 2),
           )
           F = [1+0*im 0; 0 0]
           @constraint(model, 0 <= real(LinearAlgebra.tr(F * H)) <= 1)
           @objective(model, Min, t)
           optimize!(model)
       end

 iter        p_obj        d_obj |  abs_gap    x_feas    z_feas |      tau       kap        mu | dir_res     prox  step     alpha
    0   1.7321e+00  -1.8409e+00 | 7.00e+00  3.66e-01  3.86e-01 | 1.00e+00  1.00e+00  1.00e+00 |
    1   4.9215e-01  -4.0034e-03 | 1.27e+00  5.39e-02  5.69e-02 | 1.36e+00  2.40e-01  2.00e-01 | 2.2e-16  9.5e-01  co-a  8.00e-01
    2   5.2733e-01   4.5012e-01 | 2.31e-01  7.04e-03  7.43e-03 | 1.56e+00  1.67e-02  3.22e-02 | 1.4e-15  6.1e-01  co-a  8.50e-01
    3   5.0125e-01   4.9710e-01 | 1.04e-02  3.82e-04  4.03e-04 | 1.44e+00  8.94e-04  1.47e-03 | 5.3e-16  4.3e-01  co-a  9.50e-01
    4   5.0004e-01   4.9991e-01 | 2.95e-04  1.21e-05  1.28e-05 | 1.36e+00  2.74e-05  4.15e-05 | 7.3e-15  4.4e-01  co-a  9.70e-01
    5   5.0000e-01   5.0000e-01 | 8.29e-06  3.91e-07  4.13e-07 | 1.26e+00  7.86e-07  1.16e-06 | 2.1e-13  6.3e-01  co-a  9.70e-01
    6   5.0000e-01   5.0000e-01 | 3.72e-07  2.15e-08  2.27e-08 | 1.15e+00  4.24e-08  5.26e-08 | 4.5e-12  3.3e-01  co-a  9.50e-01
    7   5.0000e-01   5.0000e-01 | 1.06e-08  6.84e-10  7.20e-10 | 1.09e+00  1.26e-09  1.50e-09 | 1.9e-11  3.7e-01  co-a  9.70e-01
optimal solution found; terminating

status is Optimal after 7 iterations and 0.007 seconds


julia> termination_status(model)
OPTIMAL::TerminationStatusCode = 1

julia> objective_value(model)
0.5000000022401856

julia> value.(H)
2×2 Matrix{ComplexF64}:
 0.5+0.0im  0.0+0.0im
 0.0-0.0im  0.5+0.0im
4

Hi @odow, thanks for the overview!

I was indeed looking at the Hypatia cones and did not realize you could use them directly in the @constraint macros in this way. That is very cool.

I will try some other toy problems in this vein. Your example produces the \rho I expected for this example. But out of an abundance of caution for now, I guess I’ll hold off on marking it as a solution until I can test how the Hypatia cone expects the point to be passed in (i.e., is the real/imag vector stacking there right). Still working out Julia syntax and reading the Hypatia/conref docs…

Thanks so much!

1 Like
5

I’ll give this a go. I will have to take a look at that mat_to_vec function.

I am aware of some vector stacking function that may be related by @araujoms in his QKD cone. Will puzzle over these things.

6

I think I’ve convinced myself it should be the interleaved mat_to_vec approach:

julia> using LinearAlgebra

julia> using JuMP

julia> using Hypatia

julia> mat_to_vec(h::AbstractMatrix) = mat_to_vec(real.(h), imag.(h))
mat_to_vec (generic function with 3 methods)

julia> function mat_to_vec(HR::AbstractMatrix{T}, HI::AbstractMatrix{T}) where {T}
           y = Vector{T}(undef, 2 * length(HR))
           for (k, h_k) in enumerate(zip(HR, HI))
               y[2k-1], y[2k] = h_k
           end
           return y
       end
mat_to_vec (generic function with 3 methods)

julia> begin
           model = Model(Hypatia.Optimizer)
           @variable(model, H[1:2, 1:2] in HermitianPSDCone())
           @constraint(model, H .== [4 -2+2im; -2-2im 5])
           @variable(model, t)
           @constraint(
               model,
               [t; mat_to_vec(H)] in Hypatia.EpiNormSpectralCone{Float64,ComplexF64}(2, 2),
           )
           @objective(model, Min, t)
           optimize!(model)
       end
2 of 6 primal equality constraints are dependent

 iter        p_obj        d_obj |  abs_gap    x_feas    z_feas |      tau       kap        mu | dir_res     prox  step     alpha
    0   1.7321e+00  -7.5681e+00 | 5.00e+00  3.66e-01  8.33e-01 | 1.00e+00  1.00e+00  1.00e+00 |
    1   6.1752e+00   3.5155e+00 | 6.76e-01  1.92e-01  4.36e-01 | 2.87e-01  7.83e-01  1.50e-01 | 7.1e-15  9.7e-01  co-a  8.50e-01
    2   6.7379e+00   6.1896e+00 | 8.84e-02  3.58e-02  8.15e-02 | 2.30e-01  1.06e-01  1.88e-02 | 2.8e-16  8.5e-01  co-a  8.50e-01
    3   7.3623e+00   7.2929e+00 | 1.00e-02  3.63e-03  8.27e-03 | 2.27e-01  7.46e-03  1.96e-03 | 2.3e-15  6.8e-01  co-a  9.00e-01
    4   7.3701e+00   7.3627e+00 | 9.01e-04  4.00e-04  9.11e-04 | 2.06e-01  8.10e-04  1.78e-04 | 5.6e-16  4.1e-01  co-a  9.00e-01
    5   7.3723e+00   7.3718e+00 | 4.33e-05  2.12e-05  4.83e-05 | 1.94e-01  3.65e-05  8.39e-06 | 3.9e-15  8.2e-01  co-a  9.50e-01
    6   7.3723e+00   7.3722e+00 | 3.82e-06  2.37e-06  5.40e-06 | 1.73e-01  3.92e-06  7.50e-07 | 6.1e-14  4.9e-01  co-a  9.00e-01
    7   7.3723e+00   7.3723e+00 | 3.53e-07  2.54e-07  5.79e-07 | 1.62e-01  4.14e-07  7.01e-08 | 2.2e-13  2.2e-01  co-a  9.00e-01
    8   7.3723e+00   7.3723e+00 | 1.04e-08  7.88e-09  1.79e-08 | 1.57e-01  1.17e-08  2.03e-09 | 1.9e-12  5.1e-01  co-a  9.70e-01
    9   7.3723e+00   7.3723e+00 | 2.83e-09  2.52e-09  5.74e-09 | 1.47e-01  4.12e-09  5.72e-10 | 4.0e-11  3.0e-01  co-a  7.00e-01
optimal solution found; terminating

status is Optimal after 9 iterations and 0.127 seconds


julia> value(t)
7.372281291214951

julia> maximum(LinearAlgebra.svdvals(value.(H)))
7.372281323269014
7

yes

1 Like