omnibias.sos.rounding
Sum-of-squares constructions, polynomial positivity checks, and rigorous validation of proposed Gram matrices.
import omnibias.sos.roundingModule contract & mathematical conventions
Peyrl-Parrilo rational projection of a float Gram onto exact coefficients.
The SDP proposes a floating-point Gram Q_float with z^T Q_float z approx p.
To turn that into a *proof* the Gram must match p's coefficients **exactly**
(over the rationals) while staying positive definite. This module does the exact
part: round Q_float to a rational matrix at a chosen denominator, then project
that rational matrix orthogonally (in the trace inner product) onto the affine set
{Q : z^T Q z = p}.
The projection is unusually clean here because each Gram entry (i, j) feeds
**exactly one** product monomial alpha = basis[i] + basis[j], so the
coefficient-matching constraints act on **disjoint** groups of entries. The
orthogonal projection therefore decouples per monomial alpha:
delta_alpha = (p_alpha - <A_alpha, Q>) / <A_alpha, A_alpha>
is added to every entry of group alpha (A_alpha is the 0/1 symmetric
indicator of the group), a single exact rational correction per monomial. All of
this is fractions.Fraction arithmetic -- no floating point, so the
resulting matrix matches the coefficients to the last bit. Whether it is still
positive definite is decided downstream by the rigorous interval LDL^T check.
Shared formulas do not imply universal floating-point identity. Verify dtype, shape, derivative-order support, and the stated numerical assumptions for the path you use.
RationalGram
omnibias.sos.rounding.RationalGramSource ↗RationalGram = list[list[Fraction]]This exported constant is defined by the expression above. Its exact value and annotation come from the pinned source; use it with the module contract.
Read the implementation
RationalGram = list[list[Fraction]]exact_coefficient_residual
omnibias.sos.rounding.exact_coefficient_residualSource ↗exact_coefficient_residual(gram: Sequence[Sequence[Fraction]], polynomial: CoefficientSource, basis: Sequence[Exponent]) -> FractionMax |coeff(z^T Q z) - coeff(p)| as an exact rational (0 iff matched).
A soundness guard: certification only proceeds when this is exactly 0, so a
bug in the projection can never let an unmatched Gram masquerade as a proof.
Parameters
| Name | Type & default | Contract |
|---|---|---|
grampositional-or-keyword | Sequence[Sequence[Fraction]]Required | See the operation contract and implementation below. |
polynomialpositional-or-keyword | CoefficientSourceRequired | See the operation contract and implementation below. |
basispositional-or-keyword | Sequence[Exponent]Required | See the operation contract and implementation below. |
FractionExecuted regression · is sos helper agrees
Recorded coverage confirms this test executed the API body. This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.
def test_is_sos_helper_agrees() -> None:
assert is_sos(X * X - X * Y + Y * Y + ONE)
assert not is_sos(Polynomial(2, {(4, 2): 1.0, (2, 4): 1.0, (2, 2): -3.0, (0, 0): 1.0}))python -m pytest "packages/omnibias-sos/tests/test_certify.py::test_is_sos_helper_agrees" -qRead the implementation
def exact_coefficient_residual(
gram: Sequence[Sequence[Fraction]],
polynomial: CoefficientSource,
basis: Sequence[Exponent],
) -> Fraction:
r"""Max ``|coeff(z^T Q z) - coeff(p)|`` as an exact rational (``0`` iff matched).
A soundness guard: certification only proceeds when this is exactly ``0``, so a
bug in the projection can never let an unmatched Gram masquerade as a proof.
"""
residual = Fraction(0)
products = gram_products(basis)
monomials = set(products) | set(polynomial.support)
for alpha in monomials:
pairs = products.get(alpha, [])
value = sum((mult * gram[i][j] for i, j, mult in pairs), Fraction(0))
diff = abs(value - Fraction(polynomial.coefficient(alpha)))
residual = max(residual, diff)
return residualmin_norm_correction
omnibias.sos.rounding.min_norm_correctionSource ↗min_norm_correction(matrix: Sequence[Sequence[Fraction]], residual: Sequence[Fraction]) -> list[Fraction] | NoneExact delta with matrix @ delta = residual in the row space of matrix.
Solves the normal equations (A A^T) z = residual exactly and returns
delta = A^T z (the correction lying in A's row space, hence the small
perturbation used to project rounded Gram entries back onto the exact
coefficient-matching set). None if the normal equations are inconsistent.
Parameters
| Name | Type & default | Contract |
|---|---|---|
matrixpositional-or-keyword | Sequence[Sequence[Fraction]]Required | See the operation contract and implementation below. |
residualpositional-or-keyword | Sequence[Fraction]Required | See the operation contract and implementation below. |
list[Fraction] | NoneExecuted regression · interval constraint is proved and exact
Recorded coverage confirms this test executed the API body. This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.
def test_interval_constraint_is_proved_and_exact() -> None:
p = 2.0 * ONE1 - X1 * X1
gs = [ONE1 - X1 * X1] # [-1, 1]
cert = certify_nonneg_on_set(p, gs)
assert cert.certified
assert _identity_holds(cert, p, gs)
assert _min_in_set(p, gs) >= -1e-9
assert is_nonneg_on_set(p, gs)python -m pytest "packages/omnibias-sos/tests/test_positivstellensatz.py::test_interval_constraint_is_proved_and_exact" -qRead the implementation
def min_norm_correction(
matrix: Sequence[Sequence[Fraction]], residual: Sequence[Fraction]
) -> list[Fraction] | None:
r"""Exact ``delta`` with ``matrix @ delta = residual`` in the row space of ``matrix``.
Solves the normal equations ``(A A^T) z = residual`` exactly and returns
``delta = A^T z`` (the correction lying in ``A``'s row space, hence the small
perturbation used to project rounded Gram entries back onto the exact
coefficient-matching set). ``None`` if the normal equations are inconsistent.
"""
n_eq = len(matrix)
if n_eq == 0:
return []
n_var = len(matrix[0])
gram = [
[sum((matrix[i][k] * matrix[j][k] for k in range(n_var)), Fraction(0)) for j in range(n_eq)]
for i in range(n_eq)
]
z = solve_rational_system(gram, list(residual))
if z is None:
return None
return [sum((matrix[i][c] * z[i] for i in range(n_eq)), Fraction(0)) for c in range(n_var)]project_coefficients_exact
omnibias.sos.rounding.project_coefficients_exactSource ↗project_coefficients_exact(rounded: Sequence[Fraction], matrix: Sequence[Sequence[Fraction]], target: Sequence[Fraction]) -> list[Fraction] | NoneProject rounded onto {u : matrix @ u = target} exactly (min-norm step).
Parameters
| Name | Type & default | Contract |
|---|---|---|
roundedpositional-or-keyword | Sequence[Fraction]Required | See the operation contract and implementation below. |
matrixpositional-or-keyword | Sequence[Sequence[Fraction]]Required | See the operation contract and implementation below. |
targetpositional-or-keyword | Sequence[Fraction]Required | See the operation contract and implementation below. |
list[Fraction] | NoneExecuted regression · interval constraint is proved and exact
Recorded coverage confirms this test executed the API body. This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.
def test_interval_constraint_is_proved_and_exact() -> None:
p = 2.0 * ONE1 - X1 * X1
gs = [ONE1 - X1 * X1] # [-1, 1]
cert = certify_nonneg_on_set(p, gs)
assert cert.certified
assert _identity_holds(cert, p, gs)
assert _min_in_set(p, gs) >= -1e-9
assert is_nonneg_on_set(p, gs)python -m pytest "packages/omnibias-sos/tests/test_positivstellensatz.py::test_interval_constraint_is_proved_and_exact" -qRead the implementation
def project_coefficients_exact(
rounded: Sequence[Fraction],
matrix: Sequence[Sequence[Fraction]],
target: Sequence[Fraction],
) -> list[Fraction] | None:
r"""Project ``rounded`` onto ``{u : matrix @ u = target}`` exactly (min-norm step)."""
residual = [
Fraction(target[i]) - sum((matrix[i][c] * rounded[c] for c in range(len(rounded))), Fraction(0))
for i in range(len(matrix))
]
delta = min_norm_correction(matrix, residual)
if delta is None:
return None
return [rounded[c] + delta[c] for c in range(len(rounded))]project_to_exact_gram
omnibias.sos.rounding.project_to_exact_gramSource ↗project_to_exact_gram(gram: Sequence[Sequence[float]], polynomial: CoefficientSource, basis: Sequence[Exponent], *, denominator: int) -> RationalGramRound gram at denominator and exactly project onto z^T Q z = p.
Returns a symmetric rational Gram matrix whose expansion z(x)^T Q z(x)
equals polynomial **exactly** (verify with exact_coefficient_residual,
which is 0 by construction). Positive-definiteness is *not* checked here.
Parameters
| Name | Type & default | Contract |
|---|---|---|
grampositional-or-keyword | Sequence[Sequence[float]]Required | See the operation contract and implementation below. |
polynomialpositional-or-keyword | CoefficientSourceRequired | See the operation contract and implementation below. |
basispositional-or-keyword | Sequence[Exponent]Required | See the operation contract and implementation below. |
denominatorkeyword-only | intRequired | See the operation contract and implementation below. |
RationalGramExecuted regression · is sos helper agrees
Recorded coverage confirms this test executed the API body. This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.
def test_is_sos_helper_agrees() -> None:
assert is_sos(X * X - X * Y + Y * Y + ONE)
assert not is_sos(Polynomial(2, {(4, 2): 1.0, (2, 4): 1.0, (2, 2): -3.0, (0, 0): 1.0}))python -m pytest "packages/omnibias-sos/tests/test_certify.py::test_is_sos_helper_agrees" -qRead the implementation
def project_to_exact_gram(
gram: Sequence[Sequence[float]],
polynomial: CoefficientSource,
basis: Sequence[Exponent],
*,
denominator: int,
) -> RationalGram:
r"""Round ``gram`` at ``denominator`` and exactly project onto ``z^T Q z = p``.
Returns a symmetric rational Gram matrix whose expansion ``z(x)^T Q z(x)``
equals ``polynomial`` **exactly** (verify with :func:`exact_coefficient_residual`,
which is ``0`` by construction). Positive-definiteness is *not* checked here.
"""
size = len(basis)
rational: RationalGram = [
[_round_fraction(gram[i][j], denominator) for j in range(size)] for i in range(size)
]
# Symmetrise exactly (rounding may perturb the two triangles differently).
for i in range(size):
for j in range(i + 1, size):
avg = (rational[i][j] + rational[j][i]) / 2
rational[i][j] = rational[j][i] = avg
for alpha, pairs in gram_products(basis).items():
target = Fraction(polynomial.coefficient(alpha))
current = sum((mult * rational[i][j] for i, j, mult in pairs), Fraction(0))
norm_sq = sum(1 if i == j else 2 for i, j, _mult in pairs)
delta = (target - current) / norm_sq
for i, j, _mult in pairs:
rational[i][j] = rational[i][j] + delta
if i != j:
rational[j][i] = rational[i][j]
return rationalsolve_rational_system
omnibias.sos.rounding.solve_rational_systemSource ↗solve_rational_system(matrix: Sequence[Sequence[Fraction]], rhs: Sequence[Fraction]) -> list[Fraction] | NoneExact solution of matrix @ x = rhs over the rationals (free vars = 0).
Gauss-Jordan elimination in ~fractions.Fraction arithmetic; returns a
particular solution (non-pivot variables set to zero) or None when the
system is inconsistent. Used for the exact coefficient-matching projections in
the Positivstellensatz and auxiliary-functional methods, where a Gram entry can
feed several monomials (so the disjoint-group shortcut does not apply).
Parameters
| Name | Type & default | Contract |
|---|---|---|
matrixpositional-or-keyword | Sequence[Sequence[Fraction]]Required | See the operation contract and implementation below. |
rhspositional-or-keyword | Sequence[Fraction]Required | See the operation contract and implementation below. |
list[Fraction] | NoneExecuted regression · interval constraint is proved and exact
Recorded coverage confirms this test executed the API body. This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.
def test_interval_constraint_is_proved_and_exact() -> None:
p = 2.0 * ONE1 - X1 * X1
gs = [ONE1 - X1 * X1] # [-1, 1]
cert = certify_nonneg_on_set(p, gs)
assert cert.certified
assert _identity_holds(cert, p, gs)
assert _min_in_set(p, gs) >= -1e-9
assert is_nonneg_on_set(p, gs)python -m pytest "packages/omnibias-sos/tests/test_positivstellensatz.py::test_interval_constraint_is_proved_and_exact" -qRead the implementation
def solve_rational_system(
matrix: Sequence[Sequence[Fraction]], rhs: Sequence[Fraction]
) -> list[Fraction] | None:
r"""Exact solution of ``matrix @ x = rhs`` over the rationals (free vars = 0).
Gauss-Jordan elimination in :class:`~fractions.Fraction` arithmetic; returns a
particular solution (non-pivot variables set to zero) or ``None`` when the
system is inconsistent. Used for the exact coefficient-matching projections in
the Positivstellensatz and auxiliary-functional methods, where a Gram entry can
feed several monomials (so the disjoint-group shortcut does not apply).
"""
n_eq = len(matrix)
n_var = len(matrix[0]) if n_eq else 0
rows = [[Fraction(v) for v in row] + [Fraction(rhs[i])] for i, row in enumerate(matrix)]
pivot_cols: list[int] = []
pivot_row = 0
for col in range(n_var):
pivot = next((i for i in range(pivot_row, n_eq) if rows[i][col] != 0), None)
if pivot is None:
continue
rows[pivot_row], rows[pivot] = rows[pivot], rows[pivot_row]
inv = rows[pivot_row][col]
rows[pivot_row] = [v / inv for v in rows[pivot_row]]
for i in range(n_eq):
if i != pivot_row and rows[i][col] != 0:
factor = rows[i][col]
rows[i] = [rows[i][k] - factor * rows[pivot_row][k] for k in range(n_var + 1)]
pivot_cols.append(col)
pivot_row += 1
if pivot_row == n_eq:
break
for i in range(n_eq):
if all(rows[i][c] == 0 for c in range(n_var)) and rows[i][n_var] != 0:
return None
solution = [Fraction(0)] * n_var
for idx, col in enumerate(pivot_cols):
solution[col] = rows[idx][n_var]
return solution