ClusteredLowRankSolver.jl

ClusteredLowRankSolver.jl provides a primal-dual interior point method for solving clustered low-rank semidefinite programs. This can be used for (semidefinite) programs with polynomial inequality constraints, which can be rewritten in terms of sum-of-squares polynomials.

Clustered Low-Rank Semidefinite Programs

A clustered low-rank semidefinite program is defined as

\[\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_* \rangle + B^T y = c \\ & Y^j \succeq 0, \end{aligned}\]

where $\langle Y^j, A^j_*\rangle$ denotes the vector with entries $\langle Y^j, A^j_p\rangle$ and the matrices $A^j_p$ have the low-rank structure

\[ A_p^j = \bigoplus_{l=1}^{L_j} \sum_{l=1}^{L_j} \sum_{r,s=1}^{R_j(l)} A_p^j(l;r, s) \otimes E_{r,s}^{R_j(l)}.\]

The matrices $A_p^j(l;r, s)$ are of low rank and $A_p^j(l;r, s)^{\sf T} = A_p^j(l;s, r)$. Here $E_{r,s}^n$ is the $n \times n$ matrix with a one at position $(r,s)$ and zeros otherwise.

One example where this structure shows up is when using polynomial constraints which are converted to semidefinite programming constraints by sampling. Such a semidefinite program with low-rank polynomial constraints is defined as

\[\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 $A^j_p(x)$ have the same structure as before but have now polynomials as entries. This can be obtained from polynomial inequality constraints by sum-of-squares characterizations. The interface currently focusses on such semidefinite programs with polynomial constraints.

The implementations contains data types for a clustered low-rank semidefinite programs (ClusteredLowRankSDP) and semidefinite programs with low-rank polynomial constraints (LowRankPolProblem), where a LowRankPolProblem can be converted into a ClusteredLowRankSDP. The implementation also contains data types for representing low-rank (polynomial) matrices as well as functions and data types for working with samples and sampled polynomials.

Installation

The solver is written in Julia, and has been registered as a Julia package. Typing using ClusteredLowRankSolver in the REPL will prompt installation if the package has not been installed yet (from Julia 1.7 onwards).

Documentation

Solver

ClusteredLowRankSolver.solvesdp — Function
solvesdp(sdp; kwargs...)

Solve the clustered SDP with low-rank constraint matrices.

Solve the following sdp:

\[\begin{aligned} \max & ∑_{j=1}^J ⟨ C^j, Y^j⟩ + ⟨ b, y ⟩ && \\ &⟨ A^{j}_*, Y^j⟩ + B^j y = c^j, && j=1,\ldots,J \\ &Y^j ⪰ 0,&& j=1,…,J, \end{aligned}\]

where we optimize over the free variables $y$ and the PSD block matrices $Y^j = diag(Y^{j,1}, ..., Y^{j,L_j})$, and $⟨A_*^j, Y^j⟩$ denotes the vector with entries $⟨A_p^j, Y^j⟩$. The matrices $A^j_p$ have the same block structure as $Y^j$. Every $A^{j,l}$ can have several equal-sized blocks $A^{j,l}_{r,s}$. The smallest blocks have a low rank structure.

Keyword arguments:

  • prec (default: precision(BigFloat)): the precision used
  • gamma (default: 0.9): the step length reduction; a maximum step length of α reduces to a step length of max(gamma*α,1)
  • beta_(in)feasible (default: 0.1 (0.3)): the amount mu is tried to be reduced by in each iteration, for (in)feasible solutions
  • omega_p/d (default: 10^10): the starting matrix variable for the primal/dual is omega_p/d*I
  • maxiterations (default: 500): the maximum number of iterations
  • duality_gap_threshold (default: 10^-15): how near to optimal the solution needs to be
  • primal/dual_error_threshold (default:10^-30): how feasible the primal/dual solution needs to be
  • max_complementary_gap (default: 10^100): the maximum of <X,Y>/#rows(X) allowed
  • need_primal_feasible/need_dual_feasible (default: false): terminate when the solution is primal/dual feasible
  • verbose (default: true): print information after every iteration if true
  • step_length_threshold (default: 10^-7): the minimum step length allowed
  • initial_solutions (default: []): if x,X,y,Y are given, use that instead of omega_p/d * I for the initial solutions

Interface

ClusteredLowRankSolver.ClusteredLowRankSDP — Type
ClusteredLowRankSDP(sos::LowRankPolProblem; as_free::Vector, prec, verbose])

Define a ClusteredLowRankSDP based on the LowRankPolProblem sos.

The PSD variables defined by the keys in the vector as_free will be modelled as extra free variables, with extra constraints to ensure that they are equal to the entries of the PSD variables. Remaining keyword arguments:

  • prec (default: precision(BigFloat)) - the precision of the result
  • verbose (default: false) - print progress to the standard output
ClusteredLowRankSolver.Constraint — Type
Constraint(constant, matrixcoeff, freecoeff, samples[, scalings])

Models a polynomial constaint of the form

\[ c(x) = ∑_l ⟨ Y^l_{r,s}, A^l_{r,s}(x) ⟩ + ∑_i y_i B_i(x)\]

with samples, where r,s are defined by the Block structure with l as first argument

Arguments:

  • constant::MPolyElem
  • matrixcoeff::Dict{Block, LowRankMatPol}
  • freecoeff::Dict{Any, MPolyElem}
  • samples::Vector{Vector}
  • scalings::Vector
ClusteredLowRankSolver.LowRankMatPol — Type
LowRankMat(eigenvalues::Vector, rightevs::Vector{Vector}[, leftevs::Vector{Vector}])

The matrix $∑_i λ_i v_i w_i^T$.

If leftevs is not specified, use leftevs = rightevs.

ClusteredLowRankSolver.Block — Type
Block(l::Any[, r::Int, s::Int])

Specifies a block corresponding to the positive semidefinite variable l.

Specifying r,s makes the Block correspond to the r,s subblock of the variable l.

ClusteredLowRankSolver.approximatefekete — Function
approximatefekete(basis, samples) -> basis, samples

Compute approximate fekete points based on samples and a corresponding orthogonal basis. The basis consists of sampled polynomials, sampled at samples

This preserves a degree ordering of basis if present.