Interface

In this section we will explain the interface, which is focused on semidefinite programs with polynomial constraints. The most general form of such a problem is

\[\begin{aligned} \min \quad & \sum_j \langle Y^j, C^j \rangle + \langle y, b\rangle \\ \text{s.t.} \quad & \langle Y^j, A^j_*(x) \rangle + B^T(x) y = c(x) \\ & Y^j \succeq 0, \end{aligned}\]

where we optimize over free scalar variables $y$ and positive semidefinite variables $Y^j = \mathrm{diag}(Y^{j,1}, \ldots, Y^{j,L_j})$. The polynomial matrices $A^j(x)$ are required to be symmetric and have a block form with low-rank and normal blocks. That is, $A^j(x)$ is block-diagonal with blocks $A^{j,l}(x)$ and

\[ A^{j,l}(x) = \sum_{r,s=1}^{N_{j,l}} A^{j,l}_{r,s}(x) \otimes E_{r,s}^{R_{j}(l)}\]

where $E_{r,s}^N$ is the $N \times N$ matrix with a one in the $(r,s)$ entry and zeros elsewhere. Furthermore, $A^{j,l}_{r,s}(x)$ can be of low rank such that $A^{j,l}_{r,s} = (A^{j,l}_{s,r})^T$. Often, the non-diagonal block structure is not used ($N=1$ for all $j,l$). One example where it is used is in polynomial matrix programs (see, e.g., David de Laat, Nando Leijenhorst (2022)).

This is then converted to a clustered low-rank semidefinite program by Sampling.

We will explain the interface using an example from polynomial optimization. Suppose we have some polynomial, e.g.,

\[f(x,y,z) = x^4 + y^4 + z^4 - 4xyz + x + y + z\]

and we want to find the minimum of this polynomial, or an upper bound on the minimum. We can relax the problem using a sum-of-squares characterization:

\[\begin{aligned} \min \quad & M & \\ \text{s.t.} \quad & f-M & \text{ is a sum-of-squares}. \\ \end{aligned}\]

Given a vector $w(x,y,z)$ of basis polynomials up to degree $d$, we can parametrize a sum-of-squares polynomial $s$ of degree $2d$ by a positive semidefinite matrix $Y$ with

\[ s = ⟨ Y, ww^{\sf T} ⟩.\]

See the examples section for more complicated examples, explaining more intricate parts of the interface.

using ClusteredLowRankSolver, AbstractAlgebra, BasesAndSamples

function min_f(d)
    # Set up the polynomial space and define f
    R, (x,y,z) = PolynomialRing(RealField, ["x", "y", "z"])
    f = x^4 + y^4 + z^4 - 4x*y*z + x + y + z

    # Define the objective
    obj = Objective(0, Dict(), Dict(:M => 1))

    # Define the constraint SOS + M = f
    # free variables
    free_dict = Dict(:M => 1)

    # PSD variables (the sum of squares polynomial) & samples
    psd_dict = Dict()
    w = basis_monomial(2d, x, y, z)
    samples = sample_points_simplex(3,2d)
    basis, samples = approximatefekete(w, samples)

    psd_dict[Block(:Y)] = LowRankMatPol([1], [basis[1:binomial(3+d,d)]])

    # the constraint
    con  = Constraint(f, psd_dict, free_dict, samples)

    # Define the polynomial program (maximize obj s.t. con)
    pol_prob = LowRankPolProblem(true, obj, [con])

    # Convert to a clustered low-rank SDP
    sdp = ClusteredLowRankSDP(pol_prob)

    #solve the SDP
    status, sol, time, errorcode = solvesdp(sdp)
end

Objective

For the Objective, we need the constant, a dictionary with the $C^{j,l}$ matrices, and a dictionary with the entries of $b$. In this case, the constant is 0, the matrix $C$ is the zero matrix, which can be omitted, and the vector $b$ has one entry corresponding to the free variable $M$:

    # Define the objective
    obj = Objective(0, Dict(), Dict(:M => 1))

Constraints

For a constraint, we need the low-rank matrices $A^{j,l}_{r,s}(x)$, the entries of $B^{\sf T}(x)$ and the polynomial $c(x)$. Furthermore, we require a unisolvent set of sample points. In our case, $c = f$, the entry of $B$ corresponding to $M$ is the constant polynomial $1$, and the low-rank matrix corresponding to our PSD matrix variables is the matrix $ww^{\sf T}$. With

    w = monomial_basis(2d, x, y, z)

we use BasesAndSamples to create the monomial basis for w. Then we create a unisolvent set of sample points for polynomials in three variables up to degree $2d$ by

    samples = sample_points_simplex(3,2d)

In general, this is not a very good combination. We can improve the basis with approximatefekete, which orthogonalizes the basis with respect to the sample points. If samples would contain more samples than needed, this would also select a good subset of the samples.

    basis, samples = approximatefekete(w, samples)

returns the basis in evaluated form, which we can only evaluate on samples from samples. Common operations such as multiplications and additions work with these sampled polynomials, but for example extracting the degree is not possible since that requires expressing the polynomials in a graded basis. However, if w is a basis ordered on degree, basis will have the same ordering. To create the LowRankMatPol for the sum-of-squares polynomial, we use

    psd_dict[Block(:Y)] = LowRankMatPol([R(1)], [basis[1:binomial(3+d,d)]])

Using Block(:Y) as key indicates that we use the $(r,s) = (1,1)$ block of the matrix :Y; Block(:Y, r, s) would use the (r,s) subblock. The size of the matrices is implicitely given by the size of the LowRankMatPol, which is defined by the prefactors and the rank one terms. In this case, we want to use the basis up to degree d to get a sum-of-squares polynomial up to degree 2d.

The command

    con  = Constraint(f, psd_dict, free_dict, samples)

creates the constraint.

Low rank polynomial programs

With

    pol_prob = LowRankPolProblem(true, obj, [con])

we create a polynomial problem (LowRankPolProblem). The first argument indicates whether we maximize (true) or minimize (false) the objective (obj) with respect to the constraints (in this case one constraint con).

Clustered low rank semidefinite programs

Finally, we convert the polynomial program to a clustered low-rank semidefinite program by sampling with

    sdp = ClusteredLowRankSDP(pol_prob)

which we can solve with

    status, sol, time, errorcode = solvesdp(sdp)

This function has multiple options; see the section with solver Options.

Retrieving variables from the solution

The solver outputs the solution sol in the form of a CLRSResults struct, which may be used to retrieve variables and the objective from the solver. In the following we want to obtain the (dual) objective of a solution, and the value of a free variable named :a.

    objective = sol.dual_objective
    a_value = sol.freevar[:a]

Similarly, we can retrieve the matrix variable :Y with sol.matrixvar[:Y].