Solving an SDP repeatedly where the matrix data changes in each repetition

Specific Domains Optimization (Mathematical)
1

Dear All,

I am trying to solve the following semidefinite optimization problem in JuMP that I want to solve repeatedly:

\begin{align*} p^{\star}(\alpha)=\text{minimize}_{X,y} &\; \mathbf{tr}CX+\sum_{i=1}^{n}c_{i}y_{i}\\ \text{subject to } & X\succeq0,\\ & X+\sum_{i=1}^{n}y_{i}A_{i}(\alpha) + B(\alpha) \succeq0 ,\\ & X\in\mathbb{S}^{m\times m},\;y\in\mathbb{R}^{n}, \end{align*}

where in each repetition, the model stays the same except the symmetric matrices A_{i}(\alpha), B_i(\alpha) change, which in turn depends on a vector \alpha\in\mathbb{R}^{p}. The vector \alpha changes for each repetition and the rules that generate the matrices A_{i}(\alpha), B_{i}(\alpha) from \alpha is fairly complicated (solution to another SDP). In my current implementation, I am rebuilding the model from scratch every time I solve it for a different \alpha, which seems to be expensive. What would be the best way of solving such a problem in JuMP without recreating the model every time?

3

What solver? Have you benchmarked the overhead of creating the JuMP model?

5

I am using Mosek. No, I have not benchmarked, but this SDP solution is being passed to Gurobi as part of a heuristic solution, so I want the solution time to be as small as possible. Building the model every time I call the heuristic solution will be a bottleneck for my use case.

6

I’d benchmark it first. You can also improve things by formulating the model close to the form that Mosek expects.

Use print_active_bridges(model) to see how we’re reformulating. Here’s a tutorial: Example: ellipsoid approximation · JuMP

You probably need something like this:

model = Model(Mosek.Optimizer)
@variable(model, X[1:m, 1:m], PSD)
@variable(model, y[1:n])
@variable(model, Z[1:m, 1:m], PSD)
@constraint(model, c_con, X + sum(y[i] * A[i] for i in 1:m) + B .== Z)

You could also benchmark doing something like:

for i in 1:m
    set_normalized_coefficient.(c_con, y[i], A_new[i])
end
set_normalized_rhs.(c_con, -B_new)

Note that you need -B because JuMP will move it to the right-hand side.

1 Like
7

Okay sounds great, thanks very much @odow !

1 Like