PolynomialOptimization.jl: Global POP Solver

• PolynomialOptimization.jl is a Julia package for global polynomial optimization over real and complex variables, offering automatic relaxation generation and native complex support.
• Its layered architecture combines symbolic modeling, relaxation construction, and direct solver interfaces to efficiently manage large-scale, nonconvex problems.
• Advanced algorithmic techniques like SOS-moment hierarchies and harmonic decompositions enable scalable applications in power systems and quantum information.

PolynomialOptimization.jl is a high-performance Julia package for global polynomial optimization over both real and complex variables. It encompasses several implementations spanning distinct mathematical paradigms—including SDP-based methods, harmonic/polyhedral hierarchies, and direct complex modeling—focused on automatic relaxation generation, computational efficiency, and support for advanced applications such as power systems and quantum information. PolynomialOptimization.jl is sometimes referenced interchangeably with MathProgComplex.jl in certain contexts. Its architecture, algorithms, solver backend integration, and benchmarked scalability have positioned it as a core tool for researchers addressing large, nonconvex, and structure-exploiting polynomial optimization problems (Sliwak et al., 2019, Cristancho et al., 2022, Desef, 17 Dec 2025).

1. Mathematical Frameworks

PolynomialOptimization.jl supports a broad class of polynomial optimization problems (POPs), including both real-variable and complex-variable forms with equality, inequality, and integer constraints. The most general template supported is:

$$\min_{z \in \mathbb{C}^n,\, x \in \{0,1\}^k}\; p(z, \bar{z}) = \sum_{(\alpha, \beta) \in \mathcal{A}} p_{\alpha,\beta} z^\alpha \bar{z}^\beta$$

subject to

$$g_i(z, \bar{z}) = \sum_{(\alpha,\beta) \in \mathcal{B}_i} g_{i,\alpha,\beta} z^\alpha \bar{z}^\beta \geq 0,\quad h_j(z, \bar{z}) = \sum_{(\alpha,\beta) \in \mathcal{C}_j} h_{j,\alpha,\beta} z^\alpha \bar{z}^\beta = 0.$$

This encapsulates real-variable POPs as a special case, and supports modeling standard forms for Lasserre-type SOS-moment hierarchies and their duals (moment relaxations). The package maintains full complex algebra until conversion to real variables is required by a downstream solver. This leads to a more succinct and interpretable representation at the user level and allows cross-experimentation between real and complex formulations at the API level (Sliwak et al., 2019, Desef, 17 Dec 2025).

2. Core Architecture and API

PolynomialOptimization.jl is architected as a layered system consisting of symbolic modeling, automatic relaxation construction, problem-size reduction, and direct solver integration:

• Symbolic Problem Layer: Mathematical objects are implemented via types for variables (over ℝ or ℂ), exponents, polynomials (scalar or matrix-valued), constraints, problems, and evaluation points. Arithmetic overloading enables natural construction of polynomials and the use of complex conjugation, modulus, and other invariant-preserving operations. The API covers both imperative and macro-based model construction (e.g., @polyvar for variable declaration) (Sliwak et al., 2019, Cristancho et al., 2022, Desef, 17 Dec 2025).
• Relaxation Builder Layer: The package provides several builder types (Dense, NewtonPolytope, CorrelativeSparsity, TermSparsity, and user-defined relaxations). Each constructs a relaxation by selecting appropriate monomial bases, groupings, and semidefinite block structure. The system facilitates both exact and inexact reduction strategies depending on problem sparsity, monomial support, and variable interactions (Desef, 17 Dec 2025).
• Solver Interface: Direct, low-level wrappers are provided to leading semidefinite programming solvers—Mosek, Clarabel, Hypatia, SCS, ProxSDP—selecting the optimal matrix and cone format for each solver (e.g., automatic realification of Hermitian blocks). For nonlinear MINLPs, integration is achieved via JuMP, AMPL, and, for structured power network problems, with domain-specific extensions (Sliwak et al., 2019, Desef, 17 Dec 2025).
• Postprocessing and Certificate Extraction: Automated global minimizer recovery is provided via flat-extension tests and extraction from moment matrices. SOS certificates are extracted from dual slack matrices. Direct rational recovery of solutions is implemented for small-scale problems (Desef, 17 Dec 2025).

3. Algorithmic Techniques and Hierarchies

PolynomialOptimization.jl incorporates several advanced algorithmic frameworks:

• SOS-Moment Lasserre Hierarchy: For real and complex POPs, the system generates and solves SDP relaxations, imposing positive semidefiniteness on moment and localizing matrices. For constrained problems, Putinar-style hierarchies are built automatically; moment dual formulations are likewise available (Desef, 17 Dec 2025, Quijorna, 2018).
• Polyhedral Harmonic Hierarchies: An independent framework implements nested polyhedral cones inside and outside the cone of nonnegative forms on the sphere. Optimization-free lower bounds are computable using harmonic decompositions, averaging operators, and explicit cubature integration. These methods avoid large-scale SDP solving and allow systematic tradeoff between accuracy and computational cost (Cristancho et al., 2022).
• Problem Size Reduction: Critical for scalability, the package provides:
  • Newton polytope–based monomial pruning (exact): Monomials outside the half-support convex hull of the data polynomials are discarded prior to relaxation construction.
  • Correlative sparsity and chordal decomposition (exact/inexact): Variables are grouped via variable interaction graphs and block-diagonalized to minimize PSD block sizes.
  • Term sparsity–based hierarchies (inexact): Iterative support union and clique-finding iteratively refine the necessary monomial supports.
  • Inner approximations via DD/SDD cones allow further computational savings for very large problems, with rigorous control over the relaxation gap (Desef, 17 Dec 2025).

4. Complex Variables and Semidefinite Handling

Direct support for complex-variable polynomial optimization is a unique feature:

• Complex variables and their conjugates are maintained as first-class objects in all modeling and arithmetic layers.
• Complex sums-of-squares (SOS) polynomials are supported in full generality, with Hermitian Gram matrix representations and transparent mapping to real SDPs if necessary.
• Moment matrices for complex polynomials are constructed using indices over $(\alpha, \beta)$ pairs, and the |z|² operator is preserved at the symbolic level until realification.
• The system ensures all solver backends can natively process Hermitian PSD cones or automatically perform Schur-Type realifications when a solver is restricted to real semidefinite programming (Desef, 17 Dec 2025, Sliwak et al., 2019).

5. Advanced Applications and Performance Benchmarks

PolynomialOptimization.jl has been applied to a series of benchmark problems, particularly in power systems and quantum information:

• Optimal Power Flow (OPF): The package enables direct modeling of the AC-OPF and preventive-security constrained OPF (PSCOPF) entirely in the complex domain, with automatic conversion to real MINLP or SDP representations as needed. Benchmark results show successful solution of mid-sized IEEE and RTS grid instances, with robust convergence of continuous relaxations and sensitivity of integer convergence to scaling (Sliwak et al., 2019).
• Quantum Information Problems: Application to entanglement distribution protocols over lossy channels involves complex-valued matrix polynomials, semidefinite constraints, and large-scale SDP relaxations (e.g., moment matrices 1000×1000 and larger). The package's exact and approximate dimension-reduction machinery enables efficient solution—reducing dense formulations by orders of magnitude in both construction time and memory usage. Comparisons to MATLAB toolboxes and alternative Julia frameworks show 10–100× gains in relaxation building efficiency (Desef, 17 Dec 2025).

Application Domain	Problem Size (matrix dim.)	Building/Solving Time	Main Features
Power Networks (OPF)	up to 500+ variables	< 5 min	Full complex modeling
Quantum Information	1000–8000 block size	1–5 min/40 GiB RAM	Hermitian PSD, Newton
Harmonic Min. on Sⁿ	few thousand cubature pts	0.1–0.5 s/level	Optimization-free

• Harmonic Sphere Problems: Polyhedral hierarchy approach achieves $\sim$10⁻³-level accuracy in Motzkin and Robinson polynomials at low hierarchy depth, with per-level runtimes under 0.5 s. Sequence selection (pure-power vs. Fang–Fawzi kernels) adjusts tradeoff between computational effort and bound tightness (Cristancho et al., 2022).

6. Extended Features, Limitations, and Future Developments

• Preventive-Security Constraints: Indexed state/auxiliary variables per contingency scenario, PV/PQ switching via binary indicator variables, and complementarity constraints supported natively (KNITRO interface) (Sliwak et al., 2019).
• Generic Element-Based Modeling: Problem elements are extensible Julia types; non-electrical, microgrid, or antenna array applications are possible by composition, requiring no hard-coded network structure (Sliwak et al., 2019).
• Postprocessing and Certificates: Moment-matrix extraction (flat-extension) and SOS certificate extraction available; rational solution refinement supported for problem instances with analytic minima (Desef, 17 Dec 2025).
• Barrier Oracles: Nearly-optimal algorithms for interior-point barrier evaluation and differentiation for SOS/moment cones are implemented using interpolatory/Lagrange bases, avoiding the cubic complexity of Gram-matrix methods (Desef, 17 Dec 2025).
• Limitations:
  • For $n \gtrsim 6$ or $d \gtrsim 20$ complexity grows rapidly; current code is specialized to the sphere (harmonic hierarchy).
  • First-order solvers require accuracy tuning and cannot always guarantee minimizer extraction.
  • Extension to other integration domains (simplex, cube) and to noncommutative POPs is under active development (Cristancho et al., 2022, Desef, 17 Dec 2025).
  • Exploitation of polynomial sparsity and symmetry in cubature routines, GPU acceleration of dimension-reduction, and automatic facial reduction are on the roadmap (Desef, 17 Dec 2025).

7. Context, Usage, and Relationship to Other Tools

PolynomialOptimization.jl fills a gap between abstract, general SOS/moment-based SDP formalisms and specialized, large-scale, efficient implementations tuned for power systems and quantum information. It is distinguished by:

• Complete native modeling for complex variables and Hermitian constraints.
• Automatic and rigorous size reduction for relaxations that, in other frameworks, would be numerically prohibitive or require user intervention.
• Direct interfaces to state-of-the-art solvers, bypassing traditional modeling middleware and providing fine-grained control over block structure, support, and backend features (Desef, 17 Dec 2025).
• A clear, extensible architecture for the addition of new relaxation paradigms (polyhedral, hierarchy, optimization-free, dual moment cones) and post-solution certification.

PolynomialOptimization.jl is referenced in the literature as outperforming SOSTools, YALMIP, SumOfSquares.jl, and TSSOS.jl in terms of preprocessing and relaxation-building speed for structured, large-scale problems. Its design and algorithms draw directly on recent advances in size reduction, harmonic analysis, and quantum convexity (Cristancho et al., 2022, Desef, 17 Dec 2025).