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

**URL:** https://discourse.julialang.org/t/how-to-minimize-operator-spectral-norm-of-a-density-matrix-with-jump-hypatia/129006
**Category:** Optimization (Mathematical)
**Tags:** question, jump
**Created:** [May 14, 2025, 6:10pm UTC](https://discourse.julialang.org/t/how-to-minimize-operator-spectral-norm-of-a-density-matrix-with-jump-hypatia/129006 "2025-05-14T18:10:22Z")
**Posts on this page:** 8
**Page:** 1

<div class="post-metadata">

### Author: ![Menixm](https://avatars.discourse-cdn.com/v4/letter/m/ed655f/32.png) [@Menixm](https://discourse.julialang.org/u/Menixm)
#### Post date: [May 14, 2025, 6:10pm UTC](https://discourse.julialang.org/t/how-to-minimize-operator-spectral-norm-of-a-density-matrix-with-jump-hypatia/129006/1 "2025-05-14T18:10:22Z")

</div>

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.

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

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

---

<div class="post-metadata">

### Author: ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [May 14, 2025, 9:43pm UTC](https://discourse.julialang.org/t/how-to-minimize-operator-spectral-norm-of-a-density-matrix-with-jump-hypatia/129006/2 "2025-05-14T21:43:01Z")

</div>

Hi @Menixm, welcome to the forum 😄

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

---

<div class="post-metadata">

### Author: ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [May 14, 2025, 9:57pm UTC](https://discourse.julialang.org/t/how-to-minimize-operator-spectral-norm-of-a-density-matrix-with-jump-hypatia/129006/3 "2025-05-14T21:57:46Z")

</div>

It might actually be this:

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

```

---

<div class="post-metadata">

### Author: ![Menixm](https://avatars.discourse-cdn.com/v4/letter/m/ed655f/32.png) [@Menixm](https://discourse.julialang.org/u/Menixm)
#### Post date: [May 14, 2025, 9:59pm UTC](https://discourse.julialang.org/t/how-to-minimize-operator-spectral-norm-of-a-density-matrix-with-jump-hypatia/129006/4 "2025-05-14T21:59:29Z")

</div>

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!

---

<div class="post-metadata">

### Author: ![Menixm](https://avatars.discourse-cdn.com/v4/letter/m/ed655f/32.png) [@Menixm](https://discourse.julialang.org/u/Menixm)
#### Post date: [May 14, 2025, 10:04pm UTC](https://discourse.julialang.org/t/how-to-minimize-operator-spectral-norm-of-a-density-matrix-with-jump-hypatia/129006/5 "2025-05-14T22:04:15Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [May 14, 2025, 10:12pm UTC](https://discourse.julialang.org/t/how-to-minimize-operator-spectral-norm-of-a-density-matrix-with-jump-hypatia/129006/6 "2025-05-14T22:12:19Z")

</div>

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

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

```

---

<div class="post-metadata">

### Author: ![lkapelevich](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lkapelevich/32/5672_2.png) [@lkapelevich](https://discourse.julialang.org/u/lkapelevich)
#### Post date: [May 15, 2025, 2:16am UTC](https://discourse.julialang.org/t/how-to-minimize-operator-spectral-norm-of-a-density-matrix-with-jump-hypatia/129006/7 "2025-05-15T02:16:21Z")

</div>

> [@odow](#):
>
> I think I’ve convinced myself it should be the interleaved

yes

---

<div class="post-metadata">

### Author: ![araujoms](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/araujoms/32/217734_2.png) [@araujoms](https://discourse.julialang.org/u/araujoms)
#### Post date: [May 16, 2025, 11:57am UTC](https://discourse.julialang.org/t/how-to-minimize-operator-spectral-norm-of-a-density-matrix-with-jump-hypatia/129006/8 "2025-05-16T11:57:05Z")

</div>

Note that Hypatia has two spectral norm cones, `EpiNormSpectral` for full matrices, and `EpiNormSpectralTri` for Hermitian matrices. Since here `H` is Hermitian you should use the latter, as it’s much more efficient. Although `EpiNormSpectral` uses the interleaved vectorization that @odow wrote, `EpiNormSpectralTri` uses a scaled vectorization of the upper triangular.

I do provide a function `svec` to compute it in my [QKD cone](https://github.com/araujoms/ConicQKD.jl), but it’s just a wrapper around Hypatia’s `smat_to_svec`, made to work conveniently with JuMP variables.