What software can I use to express a polynomial as a sum of squares?

Question

I have a 169 term degree 8 homogeneous polynomial in 14 variables that I believe can be expressed as a sum of 8 squares. I discovered the other day that methods exist for expressing a polynomial as a sum of squares, but I'm wondering if there is any software available that can do this calculation. Does MATLAB, Mathematica, or Sage have this capability?

I saw this MATLAB library called SOSTOOLS that solves optimization problems having to do with polynomial sums of squares, but I can't tell if it has methods that do what I'm trying to do.

Here's the math behind the method. Let $x$ be a vector of $s$ polynomials that form a basis of the space of polynomials you're working in. Then any polynomial $P$ can be expressed in the form $P=x^tMx$ for some matrix $M$. Let $L$ be the set of all matrices satisfying $x^tLx=0$. Then $\{x^t(M+l)x|l\in L\}$ is the set of all representations of $P$ in this form. If $M+l$ is positive semidefinite with rank $r$, then we can find an $r\times s$ matrix $Q$ satisfying $Q^tQ=M+l$. Then we get $P=(Qx)^t(Qx)$, so we have expressed $P$ as a sum of $r$ squares. $L$ can be parameterized linearly, so this is some sort of linear optimization problem I guess.

In my case, we can let $x$ be the vector of all monomials of degree 4 in my 14 variables. Then we are trying to find positive semidefinite $M+l$ with rank 8.

Answer 1

I don't know about performance for polynomials of that size, but first thing I would try is to solve the SDP directly.

Since the question was about software, I will show what I wrote using Sage, but I'll omit the math behind it, as that can be found in many other places. I wrote the SDP into SDPA format to run CSDP on it.

First, some utility functions. sos.py:

from typing import Union, List

from sage.matrix.all import Matrix
from sage.matrix.all import matrix
from sage.rings.monomials import monomials
from sage.rings.polynomial.multi_polynomial_element import MPolynomial
from sage.rings.polynomial.polynomial_element_generic import Polynomial_generic_field as Polynomial

AnyPolynomial = Union[MPolynomial, Polynomial]


def get_all_monomials(p: AnyPolynomial) -> List[AnyPolynomial]:
    """ Get all monomials that could be terms of p. """
    n_variables = p.nvariables() if (isinstance(p, MPolynomial)) else 1
    return list(filter(
        lambda m: m.degree() <= p.degree(),
        monomials(list(p.variables()), [p.degree() + 1] * n_variables)
    ))


def get_monomial_matrices(p: AnyPolynomial, monomial_factors: List[AnyPolynomial], all_monomials: List[AnyPolynomial]) -> List[
    Matrix]:
    """
    Construct matrices F_d, d=0,...,len(all_monomials), for a vector of monomials v, such that
    entry (i,j) is one if v_i v_j is a polynomial of degree d - 1.

    :param p: Polynomial.
    :param monomial_factors: Possible monomials of the squares of p.
    :param all_monomials: All possible monomials that could appear in p.
    """
    # All monomials get a matrix, plus the objective
    number_of_matrices = len(all_monomials) + 1
    f_i = [matrix(p.parent(), len(monomial_factors), len(monomial_factors), sparse=True) for _ in
           range(number_of_matrices)]

    for i, monomial_1 in enumerate(monomial_factors):
        for j, monomial_2 in enumerate(monomial_factors):
            product = monomial_1 * monomial_2
            if product not in all_monomials:
                continue
            index = all_monomials.index(product)
            f_i[index + 1][i, j] = 1

    return f_i


def get_sdpa_format(matrices: List[Matrix], block_number=1) -> List[str]:
    """ Convert the matrices to a list of lines in sdpa format. """
    matrix_entries = []

    for mat_index, mat in enumerate(matrices):
        for (i, j), entry in mat.dict().items():
            # We only need to write entries above the diagonal
            if i > j:
                continue

            # Matrix number, block number, i, j, entry
            matrix_entries.append(f'{mat_index} {block_number} {i + 1} {j + 1} {entry}')

    return matrix_entries

Now, the construction of the SDP:

from math import sqrt
from pathlib import Path
from typing import Union

from sage.modules.free_module_element import vector
from sage.rings.all import RDF
from sage.rings.polynomial.multi_polynomial_element import MPolynomial
from sage.rings.polynomial.polynomial_element_generic import Polynomial_generic_field as Polynomial
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing

from csdp import run_csdp, read_csdp_solution
from sos import get_all_monomials, get_monomial_matrices, get_sdpa_format


def is_sos(p: Union[MPolynomial, Polynomial]) -> bool:
    """
    Is p SOS? If yes, print the SOS decomposition.

     :param p: Univariate or multivariate polynomial.
     """
    if p.degree() % 2 == 1:
        return False

    all_monomials = get_all_monomials(p)

    # All possible factors which could contribute to terms of the squares of p
    monomial_factors = list(filter(lambda m: m.degree() <= p.degree() / 2, all_monomials))

    # F_0 is the objective
    f_i = get_monomial_matrices(p, monomial_factors, all_monomials)

    # Now we need to make sure the coefficients of p match with the monomials list
    p_coeff = [str(p.monomial_coefficient(m)) for m in all_monomials]

    assert len(f_i) - 1 == len(p_coeff)

    matrix_entries = '\n'.join(get_sdpa_format(f_i))
    filename = 'sos'

    # Only one block here
    block_sizes = [len(monomial_factors)]

    Path(f'{filename}.sdpa').write_text(
        f'{len(all_monomials)} Number of constraint matrices\n'
        f'1 Number of blocks\n'
        f'{" ".join(str(s) for s in block_sizes)} Block sizes\n'
        f'{" ".join(p_coeff)} Right-hand side\n'
        f'{matrix_entries}\n'
    )

    csdp_result = run_csdp(f'{filename}.sdpa', f'{filename}.sol')
    print(csdp_result.msg)

    if csdp_result.ret_code != 0:
        return False

    solution = read_csdp_solution(f'{filename}.sol', block_sizes=block_sizes)

    # We have only one block here, get it
    block = solution[0]

    # Compute the squares
    f = [sqrt(a) * v[0] for a, v, _ in block.eigenvectors_right()]
    v = vector(p.parent(), monomial_factors)

    # Rounding the coefficients for readability, and do not print zero terms
    precision = 2
    sos_terms = [u.inner_product(v) for u in f if
                 not u.inner_product(v).map_coefficients(lambda c: round(c, precision)).is_zero()]
    sos_terms_str = ")^2 + (".join(str(s.map_coefficients(lambda c: round(c, precision))) for s in sos_terms)

    # This should be the original p
    q = sum(t ^ 2 for t in sos_terms).map_coefficients(lambda c: round(c, precision))
    print(f'p = ({sos_terms_str})^2')
    print(f'  = {q}')

    return True


if __name__ == '__main__':
    PR = PolynomialRing(RDF, 'x, y')
    x, y = PR.gens()
    motzkin: MPolynomial = 1 + x^2 * y^4 + x^4 * y^2 - 3 * x^2 * y^2

    poly: MPolynomial = (1.3 * x ^ 2 + y) ^ 2
    result = is_sos(poly)

Example:

    poly: MPolynomial = 1.69*x^4 + 2.6*x^2*y + y^2
    is_sos(poly)  # this will print p = (-1.3*x^2 - y)^2 = 1.69*x^4 + 2.6*x^2*y + y^2

As you can see I used some methods to run csdp, but all they do is call csdp on the command line and read back the result, so I've omitted the definitions.