Show HN: Sostactic – polynomial inequalities using sums-of-squares in Lean

This is a collection of Lean4 tactics for proving polynomial inequalities via sum-of-squares (SOS) decompositions, powered by a Python backend. You can use it via Python or Lean. These tactics are significantly more powerful than nlinarith and positivity -- i.e., they can prove inequalities they cannot. In theory, they can be used to prove any of the following types of statements - prove that a polynomial is nonnegative globally - prove that a polynomial is nonnegative over a semialgebraic set (i.e., defined by a set of polynomial inequalities) - prove that a semialgebraic set is empty, i.e., that a system of polynomial inequalities is infeasible To see more about how it works, read the blog post. Full documentation is provided towards the end of this README ; feel free to point your coding agent at it. import Sostactic -- AM-GM inequality example (x y : ℝ) : 2 * x * y ≤ x^2 + y^2 := by sos_decomp -- Motzkin polynomial: ratio of SOS example (x y z : ℝ) : 0 ≤ x^4 * y^2 + x^2 * y^4 + z^6 - 3 * x^2 * y^2 * z^2 := by sos_decomp (degree := 2) -- 4x³ - 3x + 1 ≥ 0 on [0, 1] (touches 0 at x = 1/2) -- linarith, positivity, and nlinarith all fail example (x : ℝ) (h1 : 0 ≤ x) (h2 : 0 ≤ 1 - x) : 0 ≤ 4 * x^3 - 3 * x + 1 := by putinar_decomp (order := 2) -- two disjoint disks centered at (0,0) and (3,3) -- linarith and nlinarith fail example : ¬∃ x y : ℝ, 0 ≤ 1 - x^2 - y^2 ∧ 0 ≤ 1 - (x - 3)^2 - (y - 3)^2 := by rintro ⟨x, y, h1, h2⟩ schmudgen_empty (order := 1) Look at Sostactic/Examples.lean for more interesting examples! - elan (Lean 4 toolchain manager) - Python 3.10+ If your project uses lakefile.toml : [[require]] name = "sostactic" git = "https://github.com/mmaaz-git/sostactic.git" rev = "main" Or if your project uses lakefile.lean : require sostactic from git "https://github.com/mmaaz-git/sostactic.git" @ "main" Then fetch and build: lake update sostactic lake build Note: Sostactic tracks Mathlib v4.29.0-rc8 . If your project uses a different Mathlib version, you may need to align them. When creating a new project, the default Mathlib version should match. Create a virtual environment in your project root and install the dependencies: python3 -m venv .venv .venv/bin/python3 -m ensurepip .venv/bin/python3 -m pip install -r .lake/packages/Sostactic/python/requirements.txt That's it. Sostactic auto-discovers .venv/bin/python3 in your project root. Tip: Add .venv/ to your.gitignore . If your Python lives elsewhere, set the SOSTACTIC_PYTHON environment variable: export SOSTACTIC_PYTHON=/path/to/your/python3 This package has 3 possible ways to use it: as a Python API, as a Python CLI, and as a Lean package. The Lean tactics call the Python CLI to generate the certificate, and then the certificate is checked inside Lean. The Python program uses cvxpy , a convex optimization solver, due to the well-established correspondence between SOS polynomials and semidefinite programming (SDP). One issue is that the SDP solution is numerical: we use a rationalization and projecting step, with some additional numerical tricks, to get an exact answer, which is then checked in Lean. A caveat is that, due to numerical issues, exactification may fail, for a variety of reasons. In the spirit of interactive theorem provers, the solver provides suggestions about how to fix the solution by passing additional arguments. For more details on the theory and algorithms, see the blog post. The rest of this README introduces the Python API, the Python CLI, and then the Lean tactics. The entire Python program is contained in python/sos.py , with tests in python/test_sos.py . There are 5 main functions: Prove that a polynomial p is globally nonnegative by finding an exact SOS certificate. from sympy import symbols, Poly from python.sos import sos_decomp x, y = symbols('x y') result = sos_decomp(Poly(x**2 + y**2, x, y)) print(result["success"]) # True print(result["exact_sos"]) # weighted SOS terms: [(weight, poly), ...] With denom_degree_bound=0 (default), finds a pure SOS decomposition: p = w₁s₁² + w₂s₂² + ... With denom_degree_bound > 0 , finds SOS polynomials d and n such that d·p = n , handling polynomials like Motzkin's that are nonneg but not themselves SOS: x, y, z = symbols('x y z') motzkin = Poly(x**4*y**2 + x**2*y**4 + z**6 - 3*x**2*y**2*z**2, x, y, z) result = sos_decomp(motzkin, denom_degree_bound=2) print(result["exact_num_sos"]) # numerator SOS decomposition print(result["exact_den_sos"]) # denominator SOS decomposition You can also pass an explicit denom_template as a sympy expression with free parameters (e.g., a*x**2 + b*y**2 + c*z**2 ). On failure, result["diagnostics"] contains rank info and result["suggestion"] has a hint for what to try next. Prove that target >= 0 on the semialgebraic set {x : g₁(x) >= 0, ..., gₘ(x) >= 0} using Putinar's Positivstellensatz. from python.sos import putinar x = symbols('x') result = putinar( target=4*x**3 - 3*x + 1, constraints=[x, 1 - x], order=2 ) print(result["success"]) # True print(result["exact_certificate_blocks"]) # certificate blocks with multipliers and SOS The order is the global relaxation order r . Each block with multiplier h is truncated by deg(h * sigma) = 0, ..., gₘ(x) >= 0} is empty. Equivalent to putinar(-1, constraints, order=order) . from python.sos import putinar_empty x = symbols('x') result = putinar_empty(constraints=[x, -x - 1], order=1) print(result["success"]) # True — the set {x >= 0, x = 0 on [0, 1] python python/sos.py putinar --poly "4*x**3 - 3*x + 1" --constraints "x, 1-x" \ --vars "x" --order 2 # prove emptiness (--poly defaults to -1) python python/sos.py putinar --constraints "x, -x-1" --vars "x" --order 1 # prove x(1-x) >= 0 on [0, 1] python python/sos.py schmudgen --poly "x*(1-x)" --constraints "x, 1-x" \ --vars "x" --order 1 # prove emptiness python python/sos.py schmudgen --constraints "1-x**2-y**2, x**2+y**2-2" \ --vars "x y" --order 1 | Argument | Description | |---|---| --poly | Target polynomial (default "-1" for putinar /schmudgen ) | --constraints | Comma-separated list of constraint polynomials | --vars | Space-separated list of variable names | --denom-degree-bound | Denominator degree bound (for sos_decomp ) | --denom-template | Explicit denominator template (for sos_decomp ) | --order | Global relaxation order for putinar /schmudgen | --basis-overrides | Per-block SOS degree caps, e.g. "0:2,3:0" | --out | Path to write JSON certificate | The CLI outputs a JSON certificate. For a successful sos_decomp : { "success": true, "kind": "sos", "status": "optimal", "exact_sos": [ {"weight": "1", "poly": {"expr": "y", "gens": ["x", "y"]}}, {"weight": "1", "poly": {"expr": "x", "gens": ["x", "y"]}} ], "exact_identity": {"expr": "0", "gens": ["x", "y"]}, "exact_residual": {"expr": "0", "gens": ["x", "y"]}, ... } This says x² + y² = 1·y² + 1·x² . For Positivstellensatz certificates, the result includes exact_certificate_blocks with multiplier/SOS pairs for each constraint block. On failure, the output includes diagnostics and a suggestion: { "success": false, "status": "infeasible", "suggestion": "not SOS — try adding a denominator with denom_degree_bound=2" } Proves that a polynomial f is nonnegative globally by finding an SOS decomposition, optionally with a denominator to decompose it into a ratio of SOS polynomials. By Artin's theorem, any globally nonnegative polynomial can be decomposed into a ratio of SOS. Goals can be in any of these forms: 0 ≤ f , f ≥ 0 , a ≤ b , or a ≥ b . -- plain SOS example (x y : ℝ) : 0 ≤ x^2 + y^2 := by sos_decomp -- with denominator search (degree bound for the denominator) example (x y z : ℝ) : 0 ≤ x^4 * y^2 + x^2 * y^4 + z^6 - 3 * x^2 * y^2 * z^2 := by sos_decomp (degree := 2) -- with explicit denominator template example (x y z : ℝ) : 0 ≤ x^4 * y^2 + x^2 * y^4 + z^6 - 3 * x^2 * y^2 * z^2 := by sos_decomp (degree := 2) (denom_template := "a*x^2 + b*y^2 + c*z^2") -- also handles a ≤ b and a ≥ b goals example (x y : ℝ) : 2 * x * y ≤ x^2 + y^2 := by sos_decomp Proves 0 ≤ f given polynomial inequality hypotheses, using Putinar's or Schmüdgen's Positivstellensatz. Hypotheses can be in any of these forms: 0 ≤ g , g ≥ 0 , a ≤ b , or a ≥ b . -- Putinar: requires Archimedean (compact) constraints example (x : ℝ) (h1 : 0 ≤ x) (h2 : x ≤ 1) : 0 ≤ 4 * x^3 - 3 * x + 1 := by putinar_decomp (order := 2) -- Schmüdgen: requires compact constraints, uses products of gᵢ example (x : ℝ) (h1 : 0 ≤ x) (h2 : x ≤ 1) : 0 ≤ x * (1 - x) := by schmudgen_decomp (order := 1) Proves False from contradictory polynomial inequality hypotheses, showing the feasible set is empty. -- {x ≥ 0, x ≤ -1} is empty example (x : ℝ) (h1 : 0 ≤ x) (h2 : x ≤ -1) : False := by putinar_empty (order := 1) -- disjoint disks example (x y : ℝ) (h1 : x^2 + y^2 ≤ 1) (h2 : x^2 + y^2 ≥ 2) : False := by schmudgen_empty (order := 1) | Parameter | Description | |---|---| order := n | Relaxation order for putinar_decomp / schmudgen_decomp / emptiness tactics | degree := n | Degree bound for sos_decomp | denom_template := "..." | Explicit denominator template (for sos_decomp ) | basis_overrides := "0:2,3:0" | Per-block SOS degree caps used to shrink exactification search space | cert := "path.json" | Load a pre-generated certificate from file | The SDP solver can be slow for large problems. You can pre-generate certificates and load them from JSON files, so the solver doesn't re-run on every file reload: example (x y z : ℝ) : 0 ≤ x^4 * y^2 + x^2 * y^4 + z^6 - 3 * x^2 * y^2 * z^2 := by sos_decomp (cert := "certs/motzkin.json") Lean still independently verifies the certificate — the JSON file is untrusted. If one of the sub-goals does not close automatically, the tactic leaves it as a named goal for you to close interactively. For sos_decomp and the Positivstellensatz decomp tactics, the goals are: sos_identity /pstv_identity — the polynomial identity (normally closed byring_nf )sos_nonneg /pstv_nonneg — the nonnegativity claim (normally closed bypositivity ) For sos_decomp with a denominator, the goals are hD_ne , hmul , hD_nonneg , and hN_nonneg . You can close them with case : example (x y : ℝ) : 0 ≤ x^2 + y^2 := by sos_decomp case sos_identity => ring case sos_nonneg => positivity The SDP solver can be slow for large polynomials. You can use the Python CLI to generate a certificate once (e.g., overnight), save it to a JSON file, and then load it in Lean so the solver doesn't re-run on every file reload: # generate the certificate python python/sos.py sos_decomp \ --poly "x**4*y**2 + x**2*y**4 + z**6 - 3*x**2*y**2*z**2" \ --vars "x y z" --denom-degree-bound 2 --out certs/motzkin.json # or for Positivstellensatz certificates python python/sos.py putinar --poly "4*x**3 - 3*x + 1" \ --constraints "x, 1-x" --vars "x" --order 2 --out certs/chebyshev.json Then in Lean: example (x y z : ℝ) : 0 ≤ x^4 * y^2 + x^2 * y^4 + z^6 - 3 * x^2 * y^2 * z^2 := by sos_decomp (cert := "certs/motzkin.json") example (x : ℝ) (h1 : 0 ≤ x) (h2 : x ≤ 1) : 0 ≤ 4 * x^3 - 3 * x + 1 := by putinar_decomp (cert := "certs/chebyshev.json") Lean still independently verifies the certificate — the JSON file is untrusted. The backend can fail either because the relaxation is too small (the SDP is infeasible at the chosen parameters) or because a numerical SDP solution does not exactify cleanly. In both cases, it returns diagnostics and suggestions. Increase the relaxation order. If the solver returns infeasible, the truncation may be too small. Try increasing the order: -- order 1 might not be enough example ... := by putinar_decomp (order := 2) -- try order 3 example ... := by putinar_decomp (order := 3) Use a denominator. For sos_decomp , if the polynomial is nonneg but not SOS (like Motzkin's), you need a denominator: -- fails without denominator -- example ... := by sos_decomp -- works with denominator degree 2 example ... := by sos_decomp (degree := 2) Use basis overrides. When exactification fails after the SDP is already feasible, the solver suggests basis_overrides to shrink selected SOS blocks. Run once without them, then rerun explicitly with the suggested override: example ... := by putinar_decomp (order := 3) (basis_overrides := "1:2,2:3") The format is "block:degree,block:degree,..." where each pair restricts a specific block to a smaller SOS degree cap. For example, 0:2 means block 0 may only use monomials of degree at most 1. Use a denominator template. For sos_decomp with a denominator, you can guide the solver by specifying the form of the denominator with free parameters: example (x y z : ℝ) : 0 ≤ x^4 * y^2 + x^2 * y^4 + z^6 - 3 * x^2 * y^2 * z^2 := by sos_decomp (degree := 2) (denom_template := "a*x^2 + b*y^2 + c*z^2") Try Schmüdgen instead of Putinar. For constrained problems, Schmüdgen's theorem is more powerful (it uses products of all subsets of constraints). If Putinar fails, Schmüdgen may succeed: -- putinar might fail here -- example ... := by putinar_decomp (order := 2) -- schmudgen uses product multipliers example ... := by schmudgen_decomp (order := 2) Pre-generate the certificate. If the solver is too slow to run interactively, generate the certificate with the Python CLI and load it in Lean (see Pre-generating certificates above). For large polynomials, Lean's proof reconstruction can hit the default heartbeat limit. You can increase it locally: set_option maxHeartbeats 800000 in example ... := by sos_decomp (degree := 2) If one of the sub-goals does not close automatically, it will be left as a named goal (e.g., sos_identity or sos_nonneg ) for you to close manually. See Interactive goals. Lean: Mathlib (multivariate polynomials, positivity tactic) Python: sympy, numpy, cvxpy, clarabel MIT