I am implementing a somewhat large SDP with complex-valued variables in JuMP, using ComplexOptInterface. For small test cases, it seems to work, but my real interest is in larger instances (specifically, when
sn > 1 and
n is large in the function below.). For those cases, model generation takes quite a long time, sometimes more than solving time. Is there something that can be improved (i.e., that will make model generation faster or decrease the instance size) in my model description, or is this just too large of a problem?
- The full implementation with examples for some cases is at GitHub - cgois/schmidt_number: Schmidt number certification via symmetric extensions.
- I have chosen JuMP + ComplexOptInterface over Convex.jl since Convex quickly throws an OutOfMemory error for large problems.
- Line 25-29 are definitely not an ideal way of implementing PSD contraints. Is there a better solution?
"""Maximum visibility (w.r.t. random noise) s.t. the DPS criterion is certified. An objective value of 1 means feasibility (unconclusive), and < 1 means entanglement.""" function maximally_mixed_distance(state, local_dim, sn=1, n::Integer=3; ppt::Bool=true) # Constants dim = size(state, 1) noise = eye(dim) / dim^2 aux_dim = local_dim * sn # Dimension with auxiliary spaces A'B' dims = repeat([aux_dim], n + 1) # AA' dim. + `n` copies of BB'. P = kron(eye(aux_dim), symmetric_projection(aux_dim, n)) # Bosonic subspace projector. Qdim = aux_dim * binomial(n + aux_dim - 1, aux_dim - 1) # Dim. extension w/ bosonic symmetries. entangling = kron(eye(local_dim), sn .* ghz(sn), eye(local_dim)) # Entangling between A' and B'. problem = Model(SCS.Optimizer) COI = ComplexOptInterface # Adds complex number support do JuMP. COI.add_all_bridges(problem) # Optimization variables @variable(problem, 0 <= vis <= 1) Q = @variable(problem, [1:Qdim, 1:Qdim] in ComplexOptInterface.HermitianPSDCone()) if ppt # Dummy vars. to enforce PPT (not possible directly in ComplexOptInterface?). fulldim = prod(dims) PSD = @variable(problem, [1:fulldim, 1:fulldim] in ComplexOptInterface.HermitianPSDCone()) end # Constraints noisy_state = vis * state + (1 - vis) * noise @expression(problem, lifted, P * Q * P') @expression(problem, reduced, ptr(lifted, 3:n+1, dims) * entangling) @constraint(problem, noisy_state .== ptr(reduced, [2, 3], [local_dim, sn, sn, local_dim])) if ppt ssys = Int.(1:ceil(n / 2) + 1) @constraint(problem, PSD .== ptransp(lifted, ssys, dims)) end # Solution @objective(problem, Max, vis) @show problem optimize!(problem) @show solution_summary(problem, verbose=true) problem, Q end
Thanks for any suggestions!