Polynomials in SDPSL

Before discussing polynomial constraints, let us briefly discuss how polynomials can be created and manipulated.

Recall that we actually use Laurent polynomials, that is, polynomials in which variables can assume positive as well as negative powers. In SDPSL, Laurent polynomials are instances of class Polynomial. The coefficients are instances of Number, and therefore are either double-precision floating-points, or istances of mpf_class, depending on whether the double-precision version or the GMP version of the library is used.

Polynomials can be either univariate or multivariate. They are basically maps associating monomials (instances of class Monomial) with coefficients. Many operations are defined for polynomials, such as addition and multiplication.

To create polynomials, first variables have to be defined. Variables have names, but these need not be unique, since they are uniquely identified by indices.

Since polynomials can be constructed from variables, variables can be multiplied and added together in order to construct polynomials. Monomials, on the other hand, have to be explicitly constructed from variables. This is made so as to avoid ambiguity when variables are multiplied together, for instance.

Iteration over polynomials is carried out over pairs of monomials and coefficients, that is, over structures std::pair<Monomial, Number>.

The following code illustrates the basics of polynomial creation and iteration.

 #include <sdpsl.h>
 
 using namespace std;
 using namespace sdp;
 
 int main()
 {
   PolynomialVariable x("x"), y("y");
   
   /* The () operator is used to take powers: x(2) = x^2 */
   Monomial mon = Monomial(x(2)) * Monomial(y(3));
   Polynomial p = (x + 1) * (x + 2);
   Polynomial q = x(3) * (y(-2) + 3);
   Polynomial r = p * q;
 
   /* Finds out the maximum exponent of variable x in r, and the
      maximum coefficient of any monomial
   */
   int max_exp = 0;
   double max_coeff = 0;
 
   for (auto& s : r) {
     if (max_exp < s.first[x]) max_exp = s.first[x];
     if (max_coeff < s.second) max_coeff = s.second;
   }
 
   cout << "For the polynomial r = " << r << ":" << endl;
   cout << "  Maximum degree of x = " << max_exp << endl;
   cout << "  Maximum coefficient = " << max_coeff << endl;
 
   return 0;
 }

See Also: For more information on how polynomials work and are implemented, see PolynomialVariable, Monomial, and Polynomial.

Adding polynomial constraints

Adding polynomial constraints is an operation very similar to adding linear constraints, except that we might have to specify a basis used to represent the polynomial constraint as linear constraints (cf. Polynomial constraints).

Specifying constraints

The left-hand side of a polynomial constraint is an expression like

\[ \sum_{j=1}^r \langle A_j, X_j \rangle \]

where the $A_j$ are matrices whose entries are Laurent polynomials and the $X_j$ are the variables of the problem.

Such expressions are represented in SDPSL by instances of SDPPolynomialCombination, which is basically a map associating matrices of polynomials (instances of Matrix<Polynomial>) to variables of the problem.

Polynomial combinations work very similarly to linear combinations and can be specified using a similar syntax. Masks can also be specified for each matrix of polynomials using SDPMaskedPolyMatrix.

A polynomial constraint is represented by an instance of SDPPolynomialConstraint. Such instances contain the left-hand side of the constraint, which is an instance of SDPPolynomialCombination, and the right-hand side, which is a polynomial. Instances of SDPPolynomialConstraint can be created from instances of SDPPolynomialCombination and Polynomial using the operator ==.

using namespace sdp;  // We use the double-precision version
SDPVariable X(5), Y(7);
Matrix<Polynomial> A(5), B(7);
MatrixMask mask(7);
SDPPolynomialCombination C = A * X + SDPMaskedPolyMatrix(B, mask) * Y;
PolynomialVariable x("x");
Polynomial p = x + 1;
SDPPolynomialConstraint constraint = C == p;

Adding polynomial constraints to a container

We can add polynomial constraints to a constraint container without having to specify a basis for the constraints representation. In this case, the standard monomial basis is used, that is, the basis consisting of all monomials occurring on the left-hand side and on the right-hand side is used.

As an example of the use of polynomial constraints, consider the problem of fiding the global minimum of a univariate polynomial. The global minimum of a univariate polynomial $p$ equals

\[ \max \{\, \lambda : \text{$p - \lambda$ is nonnegative}\,\} = \max \{\, \lambda : \text{$p - \lambda$ is a sum-of-squares}\,\}. \]

The problem of finding the maximum above can be formulated as an SDPSL semidefinite program in a standard way. Say $p$ has degree $2d$ and consider the vector $v(x) = (1, x, x^2, \ldots, x^d)$ of monomials. Then $v(x) v(x)^{\sf T}$ is a moment matrix and the maximum above is equal to the optimal value of the problem

\[ \begin{array}{rl} {\rm maximize}&s - t\\ &\langle v(x) v(x)^{\sf T}, X\rangle + s - t = p,\\ &\text{$X$ is positive semidefinite, $s$, $t \geq 0$.} \end{array} \]

Here, we have three variables, $X$, $s$, and $t$. We see $s$ and $t$ as 1-by-1 matrix variables (we use this trick to represent the free variable $\lambda$ as the difference $s-t$ of nonnegative numbers).

The program below uses SDPSL to output the above SDP in sparse SDPA format.

 /*
   polymin.cpp
 */
 
 #include <cstdlib>
 #include <iostream>
 #include <sdpsl.h>
 
 using namespace std;
 using namespace sdp;
 
 //
 // Prototypes
 //
 static void poly_min_sdp(const PolynomialVariable& x, const Polynomial& p);
 
 /*
   main
 */
 int main()
 {
   PolynomialVariable x("x");
   Polynomial p = (x + 1) * (x + 2) * (x + 3) * (x - 2) * (x - 4) * (x - 5);
   
   poly_min_sdp(x, p);
   
   return 0;
 }
 
 
 /*
   poly_min_sdp
 
   Given a univariate polynomial p on variable x, outputs an SDP whose
   optimal value is the global minimum of p. Here we use the standard
   monomial basis.
 */
 void poly_min_sdp(const PolynomialVariable& x, const Polynomial& p)
 {
   /* We first find out the degree of the polynomial */
   int deg = 0;
   
   for (auto& c : p)
     if (c.first[x] > deg)
       deg = c.first[x];
     else if (c.first[x] < 0) {
       cerr << "Negative exponents are not supported." << endl;
       exit(-1);
     }
   
   if (deg % 2 != 0) {
     cerr << "Degree must be even!" << endl;
     exit(-1);
   }
   
   /* Then we build the problem */
   SDPVariable X(1 + deg / 2);
   SDPVariable s(1), t(1);
 
   /* First we construct a moment matrix to express the polynomial p as
    * a sum of squares. 
    */
   Matrix<Polynomial> R(1 + deg / 2);
   
   for (int i = 0; i <= deg / 2; i++)
     for (int j = 0; j <= deg / 2; j++)
       R[i][j] = x(i + j);
   
   /* Then we create the constraints */
   SDPConstraintContainer cc;
   
   cc.add_constraint(R * X + Matrix<Polynomial>(1, 1, 1) * s
                     + Matrix<Polynomial>(1, 1, -1) * t == p);
   
   /* And the problem with its objective function */
   SDPProblem problem;
 
   problem.set_objective(Matrix<double>(1, 1, 1) * s + Matrix<double>(1, 1, -1) * t);
   problem << cc;
     
   /* And we write out the SDPA file */
   cout << SDPAWriter(problem);
 }

Using other bases

Other bases can be specified when a polynomial constraint is added to a constraint container. This can have a great effect on the numerical stability of the resulting semidefinite programming problem.

Different bases are specified using classes that implement the interface of SparseBasisCoefficients. An instance of such a class must implement two functions:

basis_size: a function that returns the size of the basis in question;
basis_coefficients: a function that takes a polynomial p, a reference to a vector of numbers (either doubles or instances of mpf_class), and a reference to a VectorMask, and puts in the vector of numbers the expansion of the polynomial on the basis, putting in the mask the nonzero positions of the expansion. Function basis_coefficients returns true if it could successfully compute the expansion, and false otherwise. When it returns false, an exception is thrown by the constraint container.

In the previous example, we did not specify the monomial basis explicitly; it was deduced from the polynomials given. What really happens behind the scenes is that an instance of class SparseMonomialBasis is then used. This class is derived from SparseBasisCoefficients and implements a basis of monomials.

Such bases can also be explicitly specified. For instance, lines 72-73 of the code above could be equivalently changed to the following:

SparseMonomialBasis sc;
for (int i = 0; i <= deg; i++)
    sc.add_monomial(Monomial(x(i)));
cc.add_constraint(R * X + Matrix<Polynomial>(1, 1, 1) * s
                  + Matrix<Polynomial>(1, 1, -1) * t == p, sc);

In some applications, for efficiency reasons, it is useful to implement a specific class derived from SparseBasisCoefficients to deal with polynomial constraints. In general, however, one may use SparseBasisChanger.

This class takes a monomial basis and a basis of polynomials spanning the same subspace and computes a basis change matrix between the two bases. Later, it expands polynomials given to it in the polynomial basis specified.

As an example, suppose that instead of using the standard monomial basis in the previous example, we want to use a basis of Laguerre polynomials. This can be done by substituting lines 69-73 by the following code:

   /* We create a basis of Laguerre polynomials */
   vector<Polynomial> laguerre_basis;
   compute_laguerre_polys(deg, x, 0, laguerre_basis);
 
   /* And we use it to add the polynomial constraint */
   SDPPolynomialConstraint C = R * X + Matrix<Polynomial>(1, 1, 1) * s
     + Matrix<Polynomial>(1, 1, -1) * t == p;
 
   SparseBasisChanger sc(SparseMonomialBasis(C), laguerre_basis);
 
   SDPConstraintContainer cc;
   cc.add_constraint(C, sc);

Function compute_laguerre_polys computes all Laguerre polynomials of degree 0 upto deg and parameter 0, and puts them in the vector laguerre_basis. The call SparseMonomialBasis(C) creates a monomial basis taking all monomials occurring in the polynomials that appear in the constraint C. Finally, when we add the constraint C, we specify the basis changer sc.

Table of Contents

Polynomials in SDPSL

Adding polynomial constraints

Specifying constraints

Adding polynomial constraints to a container

Using other bases