README.md

TSSOS

TSSOS aims to provide a user-friendly and efficient tool for solving optimization problems with polynomials, which is based on the structured moment-SOS hierarchy. To use TSSOS in Julia, run

pkg> add https://github.com/wangjie212/TSSOS

Documentation

Dependencies

• JuMP
• Mosek
• CliqueTrees

Usage

Unconstrained polynomial optimization

An unconstrained polynomial optimization problem could be formulized as

$$\mathrm{inf}_{\mathbf{x}\in\mathbb{R}^n}\ f(\mathbf{x}),$$

where $f\in\mathbb{R}[\mathbf{x}]$ is a polynomial.

Taking $f=1+x_1^4+x_2^4+x_3^4+x_1x_2x_3+x_2$ as an example, to compute the first TS step of the TSSOS hierarchy, run

using TSSOS
using DynamicPolynomials
@polyvar x[1:3]
f = 1 + x[1]^4 + x[2]^4 + x[3]^4 + x[1]*x[2]*x[3] + x[2]
opt,sol,data = tssos(f, x, TS="MD")

By default, the monomial basis computed by the Newton polytope method is used. If one sets newton=false:

opt,sol,data = tssos(f, x, newton=false, TS="MD")

then the standard monomial basis will be used.

To compute higher TS steps of the TSSOS hierarchy, repeatedly run

opt,sol,data = tssos(data, TS="MD")

Options
nb: specify the first nb variables to be $\pm1$ binary variables
TS: "block" by default (maximal chordal extension), "signsymmetry" (sign symmetries), "MD" (approximately smallest chordal extension), false (invalidating term sparsity iterations)
solution: true (extract optimal solutions), false

Output
basis: monomial basis
blocksize: sizes of blocks
blocks: block structrue
GramMat: Gram matrices (set Gram=true)
flag: 0 if global optimality is certified; 1 otherwise

Constrained polynomial optimization

A constrained polynomial optimization problem could be formulized as

$$\mathrm{inf}_{\mathbf{x}\in\mathbf{K}}\ f(\mathbf{x}),$$

where $f\in\mathbb{R}[\mathbf{x}]$ is a polynomial and $\mathbf{K}$ is the basic semialgebraic set

$$\mathbf{K}\coloneqq\lbrace \mathbf{x}\in\mathbb{R}^n \mid g_i(\mathbf{x})\ge0, i=1,\ldots,m,\ h_j(\mathbf{x})=0, j=1,\ldots,\ell\rbrace,$$

for some polynomials $g_i,h_j\in\mathbb{R}[\mathbf{x}]$.

Taking $f=1+x_1^4+x_2^4+x_3^4+x_1x_2x_3+x_2$ and $\mathbf{K}\coloneqq\lbrace \mathbf{x}\in\mathbb{R}^2 \mid g(\mathbf{x})=1-x_1^2-2x_2^2\ge0, h(\mathbf{x})=x_2^2+x_3^2-1=0\rbrace$ as an example, to compute the first TS step of the TSSOS hierarchy, run

@polyvar x[1:3]
f = 1 + x[1]^4 + x[2]^4 + x[3]^4 + x[1]*x[2]*x[3] + x[2]
g = 1 - x[1]^2 - 2*x[2]^2
h = x[2]^2 + x[3]^2 - 1
pop = [f, g, h]
d = 2 # set the relaxation order
opt,sol,data = tssos(pop, x, d, numeq=1, TS="block", solution=true)

To compute higher TS steps of the TSSOS hierarchy, repeatedly run

opt,sol,data = tssos(data, TS="MD", solution=true)

Options
nb: specify the first nb variables to be $\pm1$ binary variables
TS: "block" by default (maximal chordal extension), "signsymmetry" (sign symmetries), "MD" (approximately smallest chordal extension), false (invalidating term sparsity iterations)
eqTS: the same as TS by default (type of term sparsity for equality constraints), false
normality: true (impose normality condtions), false
GroebnerBasis: true (work in the quotient ring by computing a Gröbner basis), false
solution: true (extract optimal solutions), false

One could also exploit correlative sparsity and term sparsity simultaneously.

using DynamicPolynomials
n = 6
@polyvar x[1:n]
f = 1 + sum(x.^4) + x[1]*x[2]*x[3] + x[3]*x[4]*x[5] + x[3]*x[4]*x[6] + x[3]*x[5]*x[6] + x[4]*x[5]*x[6]
pop = [f, 1 - sum(x[1:3].^2), 1 - sum(x[3:6].^2)]
order = 2 # set the relaxation order
opt,sol,data = cs_tssos(pop, x, order, numeq=0, TS="MD", solution=true) # compute the first TS step of the CS-TSSOS hierarchy
opt,sol,data = cs_tssos(data, TS="block", solution=true) # compute higher TS steps of the CS-TSSOS hierarchy

Options
nb: specify the first nb variables to be $\pm1$ binary variables
CS: "MF" by default (approximately smallest chordal extension), "NC" (not performing chordal extension), false (invalidating correlative sparsity exploitation)
TS: "block" by default (maximal chordal extension), "signsymmetry" (sign symmetries), "MD" (approximately smallest chordal extension), false (invalidating term sparsity iterations)
eqTS: the same as TS by default (type of term sparsity for equality constraints), false
MomentOne: true (add a first-order moment PSD constraint for each variable clique), false
solution: true (extract an approximately optimal solution), false

Output
basis: monomial basis
blocksize: sizes of blocks
blocks: block structrue
GramMat: Gram matrices (set Gram=true)
moment: moment matrices
flag: 0 if global optimality is certified; 1 otherwise

Exploiting symmetries

TSSOS also supports exploiting symmetries for polynomial optimization problems (thanks to SymbolicWedderburn.jl).

using DynamicPolynomials
using PermutationGroups
@polyvar x[1:3]
f = sum(x) + sum(x.^4)
pop = [f, 1 - sum(x.^2)]
order = 2 # set the relaxation order
G = PermGroup([perm"(1,2,3)", perm"(1,2)"]) # define the permutation group acting on variables
opt,basis,Gram = tssos_symmetry(pop, x, order, G)

Options
TS: "block" by default (maximal chordal extension), "MD" (approximately smallest chordal extension), false (invalidating term sparsity iterations)
eqTS: the same as TS by default (type of term sparsity for equality constraints), false
numeq: number of equality constraints

Output
basis: symmetry adapted basis Gram: Gram matrices

The AC-OPF problem

Check out example/runopf.jl and example/modelopf.jl.

Sum-of-squares optimization

TSSOS supports more general sum-of-squares optimization (including polynomial optimization as a special case):

$$\mathrm{inf}_{\mathbf{y}\in\mathbb{R}^n}\ \mathbf{c}^{\intercal}\mathbf{y}$$

$$\mathrm{s.t.}\ a_{k0}+y_1a_{k1}+\cdots+y_na_{kn}\in\mathrm{SOS},\ k=1,\ldots,m.$$

where $\mathbf{c}\in\mathbb{R}^n$ and $a_{ki}\in\mathbb{R}[\mathbf{x}]$ are polynomials. SOS constraints could be handled with the routine add_psatz!:

info = add_psatz!(model, nonneg, vars, ineq_cons, eq_cons, order, TS="block", SO=1, GroebnerBasis=false)

where nonneg is a nonnegative polynomial constrained to admit a Putinar's style SOS representation on the semialgebraic set defined by ineq_cons and eq_cons, and SO is the sparse order.

The following is a simple exmaple.

$$\mathrm{sup}\ \lambda$$

$$\mathrm{s.t.}\ x_1^2 + x_1x_2 + x_2^2 + x_2x_3 + x_3^2 - \lambda(x_1^2+x_2^2+x_3^2)=\sigma+\tau_1(x_1^2+x_2^2+y_1^2-1)+\tau_2(x_2^2+x_3^2+y_2^2-1),$$

$$\sigma\in\mathrm{SOS},\deg(\sigma)\le2d,\ \tau_1,\tau_2\in\mathbb{R}[\mathbf{x}],\deg(\tau_1),\deg(\tau_2)\le2d-2.$$

using JuMP
using MosekTools
using DynamicPolynomials
using MultivariatePolynomials
using TSSOS

@polyvar x[1:3]
f = x[1]^2 + x[1]*x[2] + x[2]^2 + x[2]*x[3] + x[3]^2
d = 2 # set the relaxation order
@polyvar y[1:2]
h = [x[1]^2 + x[2]^2 + y[1]^2-1, x[2]^2 + x[3]^2 + y[2]^2 - 1]
model = Model(optimizer_with_attributes(Mosek.Optimizer))
@variable(model, lower)
nonneg = f - lower*sum(x.^2)
info = add_psatz!(model, nonneg, [x; y], [], h, d, TS="block", GroebnerBasis=true)
@objective(model, Max, lower)
optimize!(model)

Check out example/sosprogram.jl for a more complicated example.

Compute a locally optimal solution

Moment-SOS relaxations provide lower bounds on the optimum of the polynomial optimization problem. As the complementary side, one could compute a locally optimal solution which provides an upper bound on the optimum of the polynomial optimization problem. The upper bound is useful in evaluating the quality (tightness) of those lower bounds provided by moment-SOS relaxations. In TSSOS, for a given polynomial optimization problem, a locally optimal solution could be obtained via the nonlinear programming solver Ipopt:

obj,sol,status = local_solution(data.n, data.m, data.supp, data.coe, numeq=data.numeq, startpoint=rand(data.n))

Complex polynomial optimization

TSSOS also supports solving complex polynomial optimization. See Exploiting Sparsity in Complex Polynomial Optimization for more details.

A general complex polynomial optimization problem could be formulized as

$$\mathrm{inf}_{\mathbf{z}\in\mathbf{K}}\ f(\mathbf{z},\bar{\mathbf{z}}),$$

with

$$\mathbf{K}\coloneqq\lbrace \mathbf{z}\in\mathbb{C}^n \mid g_i(\mathbf{z},\bar{\mathbf{z}})\ge0, i=1,\ldots,m,\ h_j(\mathbf{z},\bar{\mathbf{z}})=0, j=1,\ldots,\ell\rbrace,$$

where $\bar{\mathbf{z}}$ stands for the conjugate of $\mathbf{z}:=(z_1,\ldots,z_n)$, and $f, g_i, i=1,\ldots,m, h_j, j=1,\ldots,\ell$ are real-valued complex polynomials satisfying $\bar{f}=f$ and $\bar{g}_j=g_j$.

Consider the following example:

$$\mathrm{inf}\ 3-|z_1|^2-0.5\mathbf{i}z_1\bar{z}_2^2+0.5\mathbf{i}z_2^2\bar{z}_1$$

$$\mathrm{s.t.}\ z_2+\bar{z}_2\ge0,|z_1|^2-0.25z_1^2-0.25\bar{z}_1^2=1,|z_1|^2+|z_2|^2=3,\mathbf{i}z_2-\mathbf{i}\bar{z}_2=0.$$

using DynamicPolynomials
n = 2 # set the number of complex variables
@complex_polyvar z[1:n]
f = 3 - z[1]*conj(z[1]) - 0.5im*z[1]*conj(z[2])^2 + 0.5im*z[2]^2*conj(z[1])
g1 = z[2] + conj(z[2])
g2 = z[1]*conj(z[1]) - 0.25*z[1]^2 - 0.25*conj(z[1])^2 - 1
g3 = z[1]*conj(z[1]) + z[2]*conj(z[2]) - 3
g4 = im*z[2] - im*conj(z[2])
pop = [f, g1, g2, g3, g4]
order = 2 # set the relaxation order
opt,sol,data = complex_tssos(pop, z, order, numeq=3, TS="block", solution=true) # no correlative sparsity
opt,sol,data = complex_cs_tssos(pop, z, order, numeq=3, TS="block", solution=true)

Options
nb: specify the first nb complex variables to be of unit norm (satisfying $|z_i|=1$)
CS: "MF" by default (approximately smallest chordal extension), "NC" (not performing chordal extension), false (invalidating correlative sparsity exploitation)
TS: "block" by default (maximal chordal extension), "MD" (approximately smallest chordal extension), false (invalidating term sparsity iterations)
eqTS: the same as TS by default (type of term sparsity for equality constraints), false
ConjugateBasis: include conjugate variables in monomial bases (false by default)
normality: specify the normal order
MomentOne: true (add a first-order moment PSD constraint for each variable clique), false

Sums of rational functions optimization

The sum-of-rational-functions optimization problem could be formulized as

$$\mathrm{inf}_{\mathbf{x}\in\mathbf{K}}\ \sum_{i=1}^N\frac{p_i(\mathbf{x})}{q_i(\mathbf{x})},$$

where $p_i,q_i\in\mathbb{R}[\mathbf{x}]$ are polynomials and $\mathbf{K}$ is the basic semialgebraic set

$$\mathbf{K}\coloneqq\lbrace \mathbf{x}\in\mathbb{R}^n \mid g_i(\mathbf{x})\ge0, i=1,\ldots,m,\ h_j(\mathbf{x})=0, j=1,\ldots,\ell\rbrace,$$

for some polynomials $g_i,h_j\in\mathbb{R}[\mathbf{x}]$.

Taking $\frac{p_1}{q_1}=\frac{x^2+y^2-yz}{1+2x^2+y^2+z^2}$, $\frac{p_2}{q_2}=\frac{y^2+x^2z}{1+x^2+2y^2+z^2}$, $\frac{p_3}{q_3}=\frac{z^2-x+y}{1+x^2+y^2+2z^2}$, and $\mathbf{K}\coloneqq\lbrace \mathbf{x}\in\mathbb{R}^2 \mid g=1-x^2-y^2-z^2\ge0\rbrace$ as an example, run

@polyvar x y z
p = [x^2 + y^2 - y*z, y^2 + x^2*z, z^2 - x + y] # define the vector of denominators
q = [1 + 2x^2 + y^2 + z^2, 1 + x^2 + 2y^2 + z^2, 1 + x^2 + y^2 + 2z^2] # define the vector of numerator
g = [1 - x^2 - y^2 - z^2]
d = 2 # set the relaxation order
opt = SumOfRatios(p, q, g, [], [x;y;z], d, QUIET=true, SignSymmetry=true) # Without correlative sparsity
opt = SparseSumOfRatios(p, q, g, [], [x;y;z], d, QUIET=true, SignSymmetry=true) # With correlative sparsity

Options
SignSymmetry: true (exploit sign symmetries), false GroebnerBasis: true (work in the quotient ring by computing a Gröbner basis), false

Polynomial matrix optimization

The polynomial matrix optimization problem aims to minimize the smallest eigenvalue of a polynomial matrix subject to a tuple of polynomial matrix inequalties (PMIs), which could be formulized as

$$\mathrm{inf}_{\mathbf{x}\in\mathbf{K}}\ \lambda_{\mathrm{min}}(F(\mathbf{x})),$$

where $F\in\mathbb{S}[\mathbf{x}]^p$ is a $p\times p$ symmetric polynomial matrix and $\mathbf{K}$ is the basic semialgebraic set

$$\mathbf{K}\coloneqq\lbrace \mathbf{x}\in\mathbb{R}^n \mid G_j(\mathbf{x})\succeq0, j=1,\ldots,m\rbrace,$$

for some symmetric polynomial matrices $G_j\in\mathbb{S}[\mathbf{x}]^{q_j}, j=1,\ldots,m$. Note that when $p=1$, $\lambda_{\min}(F(\mathbf{x}))=F(\mathbf{x})$. More generally, one may consider

$$\mathrm{inf}_{\mathbf{y}\in\mathbb{R}^t}\ \mathbf{c}^{\intercal}\mathbf{y}$$

$$\mathrm{s.t.}\ F_{0}(\mathbf{x})+y_1F_{1}(\mathbf{x})+\cdots+y_tF_{t}(\mathbf{x})\succeq0 \textrm{ on } K,$$

where $F_i\in\mathbb{S}[\mathbf{x}]^{p}, i=0,1,\ldots,t$ are a tuple of symmetric polynomial matrices.

The following is a simple exmaple.

using DynamicPolynomials
using TSSOS

@polyvar x[1:5]
F = [x[1]^4 x[1]^2 - x[2]*x[3] x[3]^2 - x[4]*x[5] x[1]*x[4] x[1]*x[5];
x[1]^2 - x[2]*x[3] x[2]^4 x[2]^2 - x[3]*x[4] x[2]*x[4] x[2]*x[5];
x[3]^2 - x[4]*x[5] x[2]^2 - x[3]*x[4] x[3]^4 x[4]^2 - x[1]*x[2] x[5]^2 - x[3]*x[5];
x[1]*x[4] x[2]*x[4] x[4]^2 - x[1]*x[2] x[4]^4 x[4]^2 - x[1]*x[3];
x[1]*x[5] x[2]*x[5] x[5]^2 - x[3]*x[5] x[4]^2 - x[1]*x[3] x[5]^4]
G = Vector{Matrix{Poly}}(undef, 2)
G[1] = [1 - x[1]^2 - x[2]^2 x[2]*x[3]; x[2]*x[3] 1 - x[3]^2]
G[2] = [1 - x[4]^2 x[4]*x[5]; x[4]*x[5] 1 - x[5]^2]
@time opt,sol,data = tssos(F, G, x, 3, TS="MD") # compute the first TS step of the TSSOS hierarchy
@time opt,sol,data = tssos(data, TS="MD") # compute higher TS steps of the TSSOS hierarchy

Options
CS: "MF" by default (approximately smallest chordal extension), "NC" (not performing chordal extension), false (invalidating correlative sparsity exploitation)
TS: "block" by default (maximal chordal extension), "MD" (approximately smallest chordal extension), false (invalidating term sparsity iterations)

For more examples, please check out example/pmi.jl.

Choose an SDP solver

By default, Mosek is chosen as the backend SDP solver. To tune parameters, create a JuMP model with the desired attributes:

tuned_optimizer = optimizer_with_attributes(Mosek.Optimizer, 
    "MSK_DPAR_INTPNT_CO_TOL_PFEAS" => 1e-8,  # primal feasibility tolerance
    "MSK_DPAR_INTPNT_CO_TOL_DFEAS" => 1e-8,  # dual feasibility tolerance
    "MSK_DPAR_INTPNT_CO_TOL_REL_GAP" => 1e-8,  # relative primal-dual gap tolerance
    "MSK_DPAR_OPTIMIZER_MAX_TIME" => 1e4,  # limit of running time
    "MSK_IPAR_NUM_THREADS" => 0 # number of threads available for Mosek
)

and run for instance tssos(..., optimizer=tuned_optimizer)

You may also choose other SDP solvers supported by JuMP:

using DynamicPolynomials
using TSSOS
using JuMP
using COSMO
using SCS

@polyvar x[1:6]
f = 1+sum(x.^4)+x[1]*x[2]*x[3]+x[3]*x[4]*x[5]+x[3]*x[4]*x[6]+x[3]*x[5]*x[6]+x[4]*x[5]*x[6]
pop = [f, 1-sum(x[1:3].^2), 1-sum(x[3:6].^2)]
d = 2
optimizer = optimizer_with_attributes(COSMO.Optimizer)
opt,sol,data = cs_tssos(pop, x, d, numeq=1, TS=false, Gram=true, solution=true, QUIET=true, optimizer=optimizer)
optimizer = optimizer_with_attributes(SCS.Optimizer)
opt,sol,data = cs_tssos(pop, x, d, numeq=1, TS=false, Gram=true, solution=true, QUIET=true, optimizer=optimizer)

High-precision polynomial optimization

For tssos, cs_tssos, complex_tssos, complex_cs_tssos, tssos_symmetry, and complex_tssos_symmetry, pass an empty JuMP GenericModel whose numeric type matches its solver. For example:

using DynamicPolynomials, JuMP, TSSOS, ClusteredLowRankSolver

setprecision(BigFloat, 256) do
    @polyvar x[1:1]
    model = GenericModel{BigFloat}(ClusteredLowRankSolver.Optimizer)
    set_optimizer_attribute(model, "duality_gap_threshold", big"1e-30")
    pop = [(x[1] - big"0.25")^2 + big"0.1", 1 - x[1]^2]
    opt, sol, data = cs_tssos(pop, x, 1; model, Gram=true, solution=true)
end

The same interface accepts GenericModel{BigFloat}(Clarabel.Optimizer{BigFloat}). Attainable accuracy depends on the solver and its settings; changing the model type alone does not tighten solver tolerances. Construct input coefficients with BigFloat (for example, big"0.1") or exact integers/rationals; converting an already rounded Float64 cannot recover precision. Keep model construction and solving inside the desired setprecision scope.

opt and the available Gram matrices, equality multipliers, and moment matrices retain the model's numeric type (Complex{BigFloat} for complex data). The symmetry routines also construct their adapted bases at that precision, using exact character projections followed by generic dense QR. This can use more memory than the Float64 sparse factorization. Begin the symmetry hierarchy at the desired precision: reusing a basis built at lower precision cannot recover its missing digits.

add_psatz_symmetry! takes its numeric type from its positional model argument and accepts polynomial coefficients containing that model's affine expressions, such as f - lower. Unlike the optimization entry points, this certificate builder can add constraints to an existing model.

The optional solution=true extraction and Ipopt refinement remain Float64/ComplexF64 heuristics: they use separate working copies of the moments and do not reduce the precision of the stored SDP results.

Pass either model or optimizer, and supply a fresh empty model for each call, including higher term-sparsity steps. With a supplied model, configure any dualization in its optimizer instead of passing dualize=true. Existing optimizer=... calls and default solver selection continue to use Float64.

Tips for modelling polynomial optimization problems

• When possible, explictly include a sphere/ball constraint (or multi-sphere/multi-ball constraints).
• When the feasible set is unbounded, try the homogenization technique introduced in Homogenization for polynomial optimization with unbounded sets.
• Scale the coefficients of the polynomial optimization problem to $[-1, 1]$.
• Scale the variables so that they take values in $[-1, 1]$ or $[0, 1]$.
• Try to include more (redundant) inequality constraints.

Christoffel-Darboux Kernels

Given a measure $\mu$ supported on $\Omega \in \mathbb{R}^n$, and a degree $d \in \mathbb{N}^{*}$, one can define, whenever the moment matrix $M_d^{\mu}$ is invertible, the Christoffel-Darboux kernel of order $d$:

$$ K_d^\mu: \mathbb{R}^n\times\mathbb{R}^n \to \mathbb{R}, \quad (x,y)\mapsto v_d(x)^\top (M_d^{\mu})^{-1} v_d(y), $$

where $v_d$ is the standard monomial basis of order $d$.

Then, for any $\mathbf{x} \in \mathbf{R}^n$, the Christoffel polynomial $\Lambda_d^{\mu}$ is an SOS polynomial of degree $2d$ defined as $\Lambda_d^{\mu}(\mathbf{x}):=K_d^\mu(x,x)$.

In practice, when solving a particular POP instance via moment-SOS relaxations, we have access to a sequence of pseudo-moments $y$, to which we can associate the Christoffel polynomial $\Lambda_d^{y}$. The sublevel sets of $\Lambda_d^{y}$ provide an approximation of the support of a measure that is concentrated on the global minimizers of the original POP.

For more information on how to construct Christoffel polynomials and how to use them to strengthen the bounds from (correlatively-sparse) Moment-SOS relaxations, visit CDK_Bound_Strengthening or docs/src/technique.md.

Non-commutative polynomial optimization

Visit NCTSSOS

Joint spetral radii

Visit SparseJSR

MATLAB implementation

Visit SPOT

References

[1] TSSOS: A Moment-SOS hierarchy that exploits term sparsity
[2] Chordal-TSSOS: a moment-SOS hierarchy that exploits term sparsity with chordal extension
[3] CS-TSSOS: Correlative and term sparsity for large-scale polynomial optimization
[4] TSSOS: a Julia library to exploit sparsity for large-scale polynomial optimization
[5] Sparse polynomial optimization: theory and practice
[6] Strengthening Lasserre's Hierarchy in Real and Complex Polynomial Optimization
[7] Exploiting Sign Symmetries in Minimizing Sums of Rational Functions
[8] Leveraging Christoffel-Darboux Kernels to Strenghten Moment-SOS Relaxations

Contact

Jie Wang: wangjie212@amss.ac.cn

Victor Magron: vmagron@laas.fr

Srećko Ðurašinović: srecko001@e.ntu.edu.sg