Polynomials

SAGE relaxations are well-suited to sparse polynomials, or polynomials of high degree. Sageopt provides a symbolic representation for such functions with the sageopt.Polynomial class. The Polynomial class thinks in terms of the following expression:

\[x \mapsto \sum_{i=1}^m c_i \prod_{j=1}^n {x_j}^{\alpha_{ij}}\]

i.e. with parameters \(\alpha \in \mathbb{N}^{m \times n}\), and \(c \in \mathbb{R}^m\). This page contains (1) the documentation for this class, (2) discussion on the concept of conditioning, (3) documentation for pre-built functions which assist in polynomial optimization, and (4) some advanced topics.

Polynomial objects

class sageopt.symbolic.polynomials.Polynomial(alpha, c)

A concrete representation for a multivariate polynomial \(x \mapsto \sum_{i=1}^m c_i \prod_{j=1}^n x_j^{\alpha_{ij}}\).

Polynomial objects provide the same operator overloading as the Signomial class. In particular, Polynomials implement +, -, *, /, and **.

Parameters

  • alpha (ndarray) – The returned Polynomial will have one monomial for each row in alpha. The number alpha[i, j] represents the power of variable j on the i-th monomial.

  • c (ndarray) – The number c[i] os the coefficient on the i-th monomial, where the i-th monomial is defined by alpha[i, :].

Examples

There are three ways to make Polynomial objects. The first way is to call the constructor:

alpha = np.array([[1, 0], [0, 1], [1, 1]])
c = np.array([1, 2, 3])
f = Polynomial(alpha, c)
x = np.array([-4, 7])
val = f(x)  # val = 1 * (-4) + 2 * (7) + 3 * (-4 * 7)
print(val)  # -74

You can also use Polynomial.from_dict, which maps exponent vectors (represented as tuples) to scalars:

alpha_and_c = {(1,): 2}
f = Polynomial.from_dict(alpha_and_c)
print(f(1))  # equal to 2.
print(f(-3))  # equal to -6.

The final way to construct a Polynomial is with algebraic syntax, like:

x = standard_poly_monomials(3)
f = (x[0] ** 2) * x[2] + 4 * x[1] * x[2] * x[0]
z = np.array([1, 2, 3])
val = f(z)  # val = (1**2) * 3 + 4 * (2 * 3 * 1)
print(val)  # 27.

Polynomial objects are not limited to numeric problem data for alpha and c. In fact, it’s very common to have c contain a coniclifts Expression. For example, if we started with a Polynomial f and then updated

gamma = sageopt.coniclifts.Variable()
f =  f - gamma

then f.c would be a coniclifts Expression depending on the variable gamma.

Notes

The Polynomial class subclasses the Signomial class. This is done because most algebraic operations between polynomials and signomials are identical. However it is important to remember that polynomials and signomials evaluate in very different ways.

Polynomial objects have a dictionary attribute called metadata. You can store any information you’d like in this dictionary. However, the information in this dictionary will not automatically be propogated when creating new Polynomial objects (as happens when performing arithmetic on Polynomials).

property n

The dimension of the space over which this Polynomial is defined; this Polynomial accepts inputs in \(\mathbb{R}^{n}\).

property m

The number of monomial basis functions \(x \mapsto \prod_{j=1}^n x_j^{a_j}\) used by this Polynomial.

property alpha

Has shape (m, n). Each row specifies the exponents of a monomial which appear in this Polynomial. The rows are ordered for consistency with the property c.

property c

Has shape (m,). The scalar c[i] is this Polynomial’s coefficient on the basis function lambda x: np.prod(np.power(x, alpha[i, :])).

property alpha_c

The keys of alpha_c are tuples of length n, containing real numeric types (e.g int, float). These tuples define monomial exponents. This Polynomial could be evaluated by the code snippet lambda x: np.sum([ alpha_c[a] * np.prod(np.power(x,a)) for a in alpha_c]).

property sig_rep

Return the signomial representative of the current Polynomial, as well as a list of constraints needed to enforce the relationship between the current Polynomial and the signomial representative.

Returns

  • sr (Signomial) – If this Signomial is globally nonnegative, then the current Polynomial is also globally nonnegative.

  • sr_cons (list of coniclifts Constraints) – If the current Polynomial has nonconstant coefficients (i.e. some entries of c are coniclifts Variables), then sr will also have nonconstant coefficients. In order to enforce the relationship between sr and the current Polynomial, we may require constraints between c and sr.c. Any such constraints are in this list.

property grad

A numpy ndarray of shape (n,) whose entries are Polynomials. For a numpy ndarray x, grad[i](x) is the partial derivative of this Polynomial with respect to coordinate i, evaluated at x. This array is constructed only when necessary, and is cached upon construction.

property hess

A numpy ndarray of shape (n, n), whose entries are Polynomials. For a numpy ndarray x, hess[i,j](x) is the (i,j)-th partial derivative of this Polynomial, evaluated at x. This array is constructed only when necessary, and is cached upon construction.

without_zeros()

Return a Polynomial which is symbolically equivalent to self, but which doesn’t track basis functions alpha[i,:] for which c[i] == 0.

query_coeff(a)

Returns the coefficient of the monomial lambda x: np.prod(np.power(x, a)) for this Polynomial.

constant_location()

Return the index i so that alpha[i, :] is the zero vector.

even_locations()

Return the largest ndarray evens, so that np.all(alpha[evens,:] % 2 == 0).

standard_multiplier()

Returns a Polynomial which is globally nonnegative by construction, for use as a modulator in SAGE hierarchies. The particular polynomial has exponents alpha = alpha[even_locations(), :], and a coefficient vector of all ones.

as_signomial()
Returns

f – For every elementwise positive vector x, we have self(x) == f(np.log(x)).

Return type

Signomial

static from_dict(d)

Construct a Polynomial object which represents the function:

lambda x: np.sum([ d[a] * np.prod(np.power(x, a)) for a in d]).
Parameters

d (Dict[Tuple[Float], Float]) –

Returns

s

Return type

Polynomial

polynomials.standard_poly_monomials()

Return x where x[i](z) = z[i] for all numeric z of length n.

Parameters

n (int) – The polynomials will be defined on \(R^n\).

Returns

x – An array of length n, containing Polynomial objects.

Return type

ndarray

Example

This function is useful for constructing Polynomials using algebraic syntax.

x = standard_poly_monomials(3)
f = (x[0] ** 2) * x[2] - 4 * x[1] * x[2] * x[0]

Conditioning

SAGE can naturally handle certain structured constraints in optimization problems, or nonnegativity problems. The process by which these constraints are handled is known as partial dualization. You can think of partial dualization as a type of “conditioning”, in the sense of conditional probability.

In the polynomial case, the “nice” sets \(X\) are those satisfying three properties

  1. invariance under reflection about the \(n\) hyperplanes \(H_i = \{ (x_1,\ldots,x_n) : x_i = 0 \}\).

  2. the set \(X \cap \mathbb{R}^n_{++}\) is log-convex, and

  3. the closure of \(X \cap \mathbb{R}^n_{++}\) equals \(X \cap \mathbb{R}^n_+\)

Take a look at Section 4 of MCW2019 to see why this is the case.

Sageopt is designed so users can take advantage of partial dualization without being experts on the subject. Towards this end, sageopt includes a function which can infer a suitably structured \(X\) from a given collection of polynomial equations and inequalities. That function is described below.

sageopt.relaxations.sage_polys.infer_domain(f, gts, eqs, check_feas=True)

Identify a subset of the constraints in gts and eqs which can be incorporated into conditional SAGE relaxations for polynomials. Construct a PolyDomain object from the inferred constraints.

Parameters

  • f (Polynomial) – The objective in a desired optimization problem. This parameter is only used to determine the dimension of the set defined by constraints in gts and eqs.

  • gts (list of Polynomials) – For every g in gts, there is a desired constraint that variables x satisfy g(x) >= 0.

  • eqs (list of Polynomials) – For every g in eqs, there is a desired constraint that variables x satisfy g(x) == 0.

  • check_feas (bool) – Indicates whether or not to verify that the returned PolyDomain is nonempty.

Returns

X

Return type

PolyDomain or None

The function above captures a small portion of what is possible with conditional SAGE certificates for polynomials. In order to take more full advantage of the possibilities, you will need to describe the set yourself. Refer to the Advanced topics section for more information.

Optimization

Here are sageopt’s convenience functions for polynomial optimization:

We assume the user has already read the section on polynomial conditioning. Newcomers to sageopt might benefit from reading this page in one browser window, and keeping our page of Examples open in an adjacent window. It might also be useful to have a copy of MCW2019 at hand, since that article is referenced throughout this section.

A remark: The functions described here are reference implementations. Significant speed improvements are possible if you build variants of these functions directly with sageopt’s backend: coniclifts is sageopt’s backend.

Optimization with structured constraints

poly_relaxation requires that any and all constraints are incorporated into a set \(X\) , which satisfies the properties for polynomial conditioning.

sageopt.poly_relaxation(f, X=None, form='dual', **kwargs)

Construct a coniclifts Problem instance for producing a lower bound on

\[\min\{ f(x) \,:\, x \in X \}\]

where \(X=R^{\texttt{f.n}}\) by default.

Parameters

  • f (Polynomial) – The objective function to be minimized.

  • X (PolyDomain or None) – If X is None, then we produce a bound on f over \(R^{\texttt{f.n}}\).

  • form (str) – Either form='primal' or form='dual'.

Returns

prob

Return type

sageopt.coniclifts.Problem

Notes

This function also accepts the following keyword arguments:

poly_ellint

Controls the complexity of a polynomial modulating function. Must be nonnegative.

sigrep_ellint

Controls the complexity of the SAGE certificate applied to the Lagrangian’s signomial representative.

The dual formulation does not allow that both poly_ell and sigrep_ell are greater than zero (such functionality is not implemented at this time).

Sageopt does not currently implement solution recovery for dual-form relaxations generated by this function. If you want to recover solutions, call poly_constrained_relaxation with empty lists gts=[] and eqs=[].

Optimization with arbitrary polynomial constraints

In more general settings, you can use poly_constrained_relaxation. In addition to allowing sets \(X\) described above, this function allows explicit polynomial inequality constraints (\(g(x) \geq 0\)) and equality constraints (\(g(x) = 0\)).

sageopt.poly_constrained_relaxation(f, gts, eqs, X=None, form='dual', **kwargs)

Construct a coniclifts Problem representing a SAGE relaxation for the polynomial optimization problem

\[\begin{split}\begin{align*} \min\{ f(x) :~& g(x) \geq 0 \text{ for } g \in \text{gts}, \\ & g(x) = 0 \text{ for } g \in \text{eqs}, \\ & \text{and } x \in X \} \end{align*}\end{split}\]

where \(X = R^{\texttt{f.n}}\) by default. The optimal value of this relaxation will produce a lower bound on the minimization problem described above. When form='dual', a solution to this relaxation can be used to help recover optimal solutions to the problem described above.

Parameters

  • f (Polynomial) – The objective to be minimized.

  • gts (list of Polynomials) – For every g in gts, there is a desired constraint that variables x satisfy g(x) >= 0.

  • eqs (list of Polynomials) – For every g in eqs, there is a desired constraint that variables x satisfy g(x) == 0.

  • X (PolyDomain) – If X is None, then we produce a bound on f subject only to the constraints in gts and eqs.

  • form (str) – Either form='primal' or form='dual'.

Returns

prob

Return type

sageopt.coniclifts.Problem

Notes

This function also accepts the following keyword arguments:

pint

Controls the complexity of Lagrange multipliers on explicit polynomial constraints gts and eqs. Defaults to p=0, which corresponds to scalar Lagrange multipliers.

qint

The lists gts and eqs are replaced by lists of polynomials formed by all products of <= q elements from gts and eqs respectively. Defaults to q = 1.

ellint

Controls the strength of the SAGE proof system, as applied to the Lagrangian. Defaults to ell=0, which means the primal Lagrangian must be an X-SAGE polynomial.

slacksbool

For dual relaxations, determines if constraints “mat @ vec in dual SAGE cone” is represented by “mat @ vec == temp, temp in dual SAGE cone”. Defaults to False.

For further explanation of the parameters p, q, and ell in the function above, we refer the user to the advanced topics section.

Solution recovery for polynomial optimization

Section 4.2 of MCW2019 introduces two solution recovery algorithms for dual SAGE relaxations. The main algorithm (“Algorithm 2”) is implemented by sageopt’s function poly_solrec, and the second algorithm (“Algorithm 2L”) is simply to use a local solver to refine the solution produced by the main algorithm. The exact choice of local solver is not terribly important. For completeness, sageopt includes a generic sageopt.local_refine() function which relies on the COBYLA solver.

sageopt.poly_solrec(prob, ineq_tol=1e-08, eq_tol=1e-06, skip_ls=False, **kwargs)

Recover a list of candidate solutions from a dual SAGE relaxation. Solutions are guaranteed to be feasible up to specified tolerances, but not necessarily optimal.

Parameters

  • prob (coniclifts.Problem) – A dual-form SAGE relaxation, from poly_constrained_relaxation.

  • ineq_tol (float) – The amount by which recovered solutions can violate inequality constraints.

  • eq_tol (float) – The amount by which recovered solutions can violate equality constraints.

  • skip_ls (bool) – Whether or not to skip least-squares solution recovery.

Returns

sols – A list of feasible solutions, sorted in increasing order of objective function value. It is possible that this list is empty, in which case no feasible solutions were recovered.

Return type

list of ndarrays

Notes

This function accepts the following keyword arguments:

zero_tolfloat

Used in magnitude recovery. If a component of the Lagrangian’s moment vector is smaller than this (in absolute value), pretend it’s zero in the least-squares step. Defaults to 1e-20.

heuristic_signsbool

Used in sign recovery. If True, then attempts to infer variable signs from the Lagrangian’s moment vector even when a completely consistent set of signs does not exist. Defaults to True.

all_signsbool

Used in sign recovery. If True, then consider returning solutions which differ only by sign. Defaults to True.

This function is implemented only for poly_constrained_relaxation (not poly_relaxation).

When faced with a polynomial optimization problem over nonnegative variables, one should formulate the problem in terms of signomials. This reformulation is without loss of generality from the perspective of solving a SAGE relaxation, but the local-refinement stage of solution recovery is somewhat different. The following function may be important if a polynomial optimization problem has some variable equal to zero in an optimal solution.

sageopt.local_refine_polys_from_sigs(f, gts, eqs, x0, **kwargs)
This is a helper function which …

  1. accepts signomial problem data (representative of a desired polynomial optimization problem),

  2. transforms the signomial data into equivalent polynomial data, and

  3. performs local refinement on the polynomial data, via the COBYLA solver.

Parameters

  • f (Signomial) – Defines the objective function to be minimized. From “f” we will construct a polynomial “p” where p(y) = f(np.log(y)) for all positive vectors y.

  • gts (list of Signomial) – Each defining an inequality constraint g(x) >= 0. From this list, we will construct a list of polynomials gts_poly, so that every g0 in gts has a polynomial representative g1 in gts_poly, satisfying g1(y) = g0(np.log(y)) for all positive vectors y.

  • eqs (list of Signomial) – Each defining an equality constraint g(x) == 0. From this list, we will construct a list of polynomials eqs_poly, so that every g0 in gts has a polynomial representative g1 in eqs_poly, satisfying g1(y) = g0(np.log(y)) for all positive vectors y.

  • x0 (ndarray) – An initial condition for the signomial optimization problem min{ f(x) |  g(x) >= 0 for g in gts, g(x) == 0 for g in eqs }.

Other Parameters

  • rhobeg (float) – Controls the size of COBYLA’s initial search space around y0 = exp(x0).

  • rhoend (float) – Termination criteria, controlling the size of COBYLA’s smallest search space.

  • maxfun (int) – Termination criteria, bounding the number of COBYLA’s iterations.

Returns

y – The output of COBYLA for the polynomial optimization problem min{ p(y) | g(y) >= 0 for g in gts_poly, g(y) == 0 for g in eqs_poly, y >= 0 } with initial condition y0 = exp(x0).

Return type

ndarray

Advanced topics

PolyDomain objects

class sageopt.symbolic.polynomials.PolyDomain(n, **kwargs)

Represent a set \(X \subset R^n\) satisfying the following properties:

  1. invariance under reflection about the \(n\) hyperplanes \(H_i = \{ (x_1,\ldots,x_n) : x_i = 0 \}\).

  2. the set \(X \cap R^n_{++}\) is is log convex, and

  3. the closure of \(X \cap R^n_{++}\) equals \(X \cap R^n_+\)

Such sets are used in polynomial conditional SAGE relaxations.

Parameters

n (int) – The dimension of the space in which this set lives.

Other Parameters

  • logspace_cons (list of coniclifts.constraints.Constraint) – Constraints over the variable y := log(|x|), which define this PolyDomain.

  • gts (list of callable) – Inequality constraint functions (g(x) >= 0) which can be used to represent X.

  • eqs (list of callable) – Equality constraint functions (g(x) == 0) which can be used to represent X.

  • check_feas (bool) – Whether or not to check that X is nonempty. Defaults to True.

  • log_AbK (tuple) – Specify a convex set in the coniclifts standard. log_AbK[0] is a SciPy sparse matrix. The first n columns of this matrix correspond to the variables over which this set is supposed to be defined. Any remaining columns are for auxiliary variables.

Notes

The constraint functions in gts and eqs should allow arguments where some components equal to zero. These functions can be Polynomial objects, but are not required to be.

Only one of log_AbK and logspace_cons can be provided upon construction. If more than one of these value is provided, the constructor will raise an error.

check_membership(x_val, tol)

Evaluate self.gts and self.eqs at x_val, to check if x_val belongs to this PolyDomain.

Parameters

  • x_val (ndarray) – Check if x_val belongs in this domain.

  • tol (float) – Infeasibility tolerance.

Returns

res – True iff x_val belongs to the domain represented by self, up to infeasibility tolerance tol.

Return type

bool

parse_coniclifts_constraints(logspace_cons)

Modify this PolyDomain object, so that it the log of its intersection with the positive orthant is the set of values satisfying constraints in logspace_cons.

Parameters

logspace_cons (list of coniclifts.Constraint) – The provided constraints must be defined over a single coniclifts Variable.

Reference hierarchy parameters

Here we describe the precise meanings of parameters p and ell in poly_constrained_relaxation. In primal form, poly_constrained_relaxation operates by moving explicit polynomial constraints into a Lagrangian, and attempting to certify the Lagrangian as nonnegative over X; this is a standard combination of the concepts reviewed in Section 2 of MCW2019. Parameter ell is essentially the same as in poly_relaxation: to improve the strength of the SAGE proof system, modulate the Lagrangian L - gamma by powers of the polynomial t = Polynomial(2 * L.alpha, np.ones(L.m)). Parameters p and q affect the unmodulated Lagrangian seen by poly_constrained_relaxation; this unmodulated Lagrangian is constructed with the following function.

sageopt.relaxations.sage_polys.make_poly_lagrangian(f, gts, eqs, p, q)

Given a problem

\[\begin{split}\begin{align*} \min\{ f(x) :~& g(x) \geq 0 \text{ for } g \in \text{gts}, \\ & g(x) = 0 \text{ for } g \in \text{eqs}, \\ & \text{and } x \in X \} \end{align*}\end{split}\]

construct the q-fold constraints q-gts and q-eqs, by taking all products of <= q elements from gts and eqs respectively. Then form the Lagrangian

\[L = f - \gamma - \sum_{g \, \in \, \text{q-gts}} s_g \cdot g - \sum_{g \, \in \, \text{q-eqs}} z_g \cdot g\]

where \(\gamma\) is a coniclifts Variable of dimension 1, and the coefficients on Polynomials \(s_g\) and \(z_g\) are coniclifts Variables of a dimension determined by p.

Parameters

  • f (Polynomial) – The objective in a desired minimization problem.

  • gts (list of Polynomials) – For every g in gts, there is a desired constraint that variables x satisfy g(x) >= 0.

  • eqs (list of Polynomials) – For every g in eqs, there is a desired constraint that variables x satisfy g(x) == 0.

  • p (int) – Controls the complexity of s_g and z_g.

  • q (int) – The number of folds of constraints gts and eqs.

Returns

  • L (Polynomial) – L.c is an affine expression of coniclifts Variables.

  • ineq_dual_polys (a list of pairs of Polynomials.) – If the pair (s_g, g) is in this list, then s_g is a generalized Lagrange multiplier to the constraint g(x) >= 0.

  • eq_dual_polys (a list of pairs of Polynomials.) – If the pair (z_g, g) is in this list, then z_g is a generalized Lagrange multiplier to the constraint g(x) == 0.

  • gamma (coniclifts.Variable.) – In primal-form SAGE relaxations, we want to maximize gamma. In dual form SAGE relaxations, gamma induces a normalizing equality constraint.

Notes

The Lagrange multipliers s_g and z_g share a common matrix of exponent vectors, which we call alpha_hat.

When p = 0, alpha_hat consists of a single row, of all zeros. In this case, s_g and z_g are constant Polynomials, and the coefficient vectors s_g.c and z_g.c are effectively scalars. When p > 0, the rows of alpha_hat are initially set set to all p-wise sums of exponent vectors appearing in either f, or some g in gts, or some g in eqs. Then we replace

alpha_hat = np.vstack([2 * alpha_hat, alpha_hat])
alpha_multiplier = np.unique(alpha_hat, axis=0)

This has the effect of improving performance for problems where alpha_hat would otherwise contain very few rows in the even integer lattice.