I have written a function that is supposed to output whether the problem is feasible and the Gram matrix.
using DynamicPolynomials
using JuMP
using MosekTools
using SumOfSquares
@polyvar var1 var2
"""
get_max_log_det(p)
Solve the max log det problem for a polynomial.
Returns (isfeasible, Gram matrix)
"""
function get_max_log_det(p)
d = 2 # The degree can be inferred from p
model = SOSModel()
set_optimizer(model, () -> Mosek.Optimizer())
@variable(model, Q[1:d÷2+1, 1:d÷2+1] in PSDCone())
@variable(model, t)
mons = monomials([var1, var2], d÷2)
@constraint(model, mons' * Q * mons == p)
vecQ = [Q[i, j] for j in 1:d÷2+1 for i in 1:j] # triangular portion only
@constraint(model, vcat(t, 1, vecQ) in MOI.LogDetConeTriangle(d÷2+1))
@objective(model, Max, t)
set_silent(model)
optimize!(model)
return primal_status(model) == FEASIBLE_POINT, value.(Q)
end
When I run
get_max_log_det((var1 + var2)^2)
this returns infeasible. If I remove the line
@objective(model, Max, t)
it returns feasible on the same input. What is going on?
I didn’t think about changing the solver until you mentioned it. With SCS I get the expected outcome, so I suppose this is a Mosek specific problem.
I thought one explanation may be that the interior point method of Mosek doesn’t quite work here because the only feasible point for this problem is on the constraint boundary. But that isn’t consistent with how Mosek can detect feasibility when there is no objective.
The problem is that without the objective, Mosek can stop when it finds any primal and dual feasible solution.
With the objective, it has to find the best primal solution, and it stalls when trying to do that. You could try different Mosek settings, but otherwise @blegat may have some suggestions for how to reformulate this into something a bit nicer for Mosek to solve.
The two most common explanations for Mosek running into slow progress / numerical issues are:
Bad scaling of the original problem, e.g., some coefficients are very large and some are very small.
This can lead to numerical/rounding issues and cause the optimizer to stall
Lack of strong duality / constraint qualification. The typical example is a problem with a bounded objective value but the optimum is not reached, e.g., minimize 1/x for x > 0. Mosek has other examples in their modeling cookbook (see Example 8.6, for instance).
It’s not impossible that the objective you provide causes it to loose strong duality and run into numerical issues.