S-SOS: Stochastic Sum-Of-Squares for Parametric Polynomial Optimization

Authors: richard l zhu
mathias  oster
yuehaw  khoo
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-09 17:24:29

Abstract: Global polynomial optimization is an important tool across applied mathematics, with many applications in operations research, engineering, and the physical sciences. In various settings, the polynomials depend on external parameters that may be random. We discuss a stochastic sum-of-squares (S-SOS) algorithm based on the sum-of-squares hierarchy that constructs a series of semidefinite programs to jointly find strict lower bounds on the global minimum and extracts candidates for parameterized global minimizers. We prove quantitative convergence of the hierarchy as the degree increases and use it to solve unconstrained and constrained polynomial optimization problems parameterized by random variables. By employing n-body priors from condensed matter physics to induce sparsity, we can use S-SOS to produce solutions and uncertainty intervals for sensor network localization problems containing up to 40 variables and semidefinite matrix sizes surpassing 800 x 800.

Main paper: 
 introduction: 
 Many effective nonlinear and nonconvex optimization techniques use local information to identify local minima. But it is often the case that we want to find global optima. Sum-of-squares (SOS) optimization is a powerful and general technique in this setting.The core idea is as follows: suppose we are given polynomials g 1 , . . . , g m , f where each function is on R n → R and we seek to determine the minimum value of f on the closed set S: S = {x ∈ R n | g i (x) ≥ 0 ∀ i = 1, . . . , m}. Our optimization problem is then to find inf x∈R n {f (x)|x ∈ S}.An equivalent formulation is to find the largest constant c ∈ R (i.e. the tightest lower bound) that can be subtracted from f such that f -c ≥ 0 over the set S. This reduction converts a polynomial optimization problem over a semialgebraic set to the problem of checking polynomial non-negativity. This problem is NP-hard in general [1] , therefore one instead resorts to checking if f -c is a sum-of-squares (SOS) function, e.g. in the unconstrained setting where S = R n one seeks to find some polynomials h k : R n → R such that f -c = k h 2 k . If such a decomposition can be found, then we have an easily-checkable certification that f -c ≥ 0, as all sum-of-squares are non-negative. Since the converse is not true (not all non-negative functions are sum-of-squares), this is a relaxation of the original non-negativity problem.Notably, if we restrict the h k to have maximum degree d, the search for a degree-2d SOS decomposition of a function can be automated as a semidefinite program (SDP) [2, 3, 4] . Solving this SDP for varying degrees d generates the well-known Lasserre (SOS) hierarchy. A given degree d corresponds to a particular level of the hierarchy. Solving this SDP produces a lower bound c d which has been proven to converge to the true global minimum c * = inf x f (x) as d increases, with finite convergence (c d = c * at finite d) for functions with second-order local optimality conditions [5, 6] and asymptotic convergence with milder assumptions thanks to representation theorems for positive polynomials from real algebraic geometry [7, 8, 9] . Further work has elucidated both theoretical implications [7, 3, 10, 11, 12, 13] and useful applications of SOS to disparate fields [14, 15, 16, 5, 6, 17, 18] (see further discussion in Appendix A.2).Motivated by the sum-of-squares certification for a lower bound c on a function f (x), we generalize to the case where the function to be minimized has additional parameters, i.e. f (x, ω) where x are variables and ω are parameters drawn from some probability distribution ω ∼ ν(ω). We seek a function c(ω) that is the tightest lower bound to f (x, ω) everywhere: f (x, ω) ≥ c(ω) with c(ω) → inf x f (x, ω). This setting was originally presented in [19] as a "Joint and Marginal" approach to parametric polynomial optimization. With the view that ω ∼ ν(ω) and seeking to parameterize the minimizers x * (ω) = argmin x f (x, ω), we are reminded of some of the prior work in polynomial chaos, where a system of stochastic variables is expanded into a deterministic function of those stochastic variables [20, 21] .Contributions and outline. Our primary contributions are a quantitative convergence proof for the Stochastic Sum-of-Squares (S-SOS) hierarchy of semidefinite programs (SDPs), a formulation of a new hierarchy (the cluster basis hierarchy) that uses the structure of a problem to sparsify the SDP, and numerical results on its application to the sensor network localization problem.In Section 2, we review the S-SOS hierarchy of SDPs [19] and its primal and dual formulations (Section 2.1). We then detail how different hierarchies can be constructed (Section 2.2.1). Finally, in Section 2.3 (complete proof in Appendix A. 5 .2) we specialize to compact X × Ω and outline the proof for a theorem on quantitative convergence (the gap between the optimal values of the degree-2s S-SOS SDP and the "tightest lower-bounding" optimization problem goes → 0 as s → ∞) of the S-SOS hierarchy for trigonometric polynomials on [0, 1] n × [0, 1] d following the kernel formalism of [22, 6, 23] .In Section 3 we review the hierarchy's applications in parametric polynomial minimization and uncertainty quantification, focusing on several variants of sensor network localization on X × Ω = [-1, 1] n × [-1, 1] d . We present numerical results for the accuracy of the extracted solutions that result from S-SOS, comparing to other approaches to parametric polynomial optimization, including a simple Monte Carlo-based method.2 Stochastic Sum-of-squares (S-SOS) 
notation: 
 Let P(S) be the space of polynomials on S, where S ∈ {X, Ω}. X ⊆ R n and Ω ⊆ R d , respectively, where X and Ω are (not-necessarily compact) subsets of their respective ambient spaces R n and R d . A polynomial in P(X) can be written as p(x) = α∈Z n ≥0 c α x α ∈ P(X) (substituting n → d, x → ω, X → Ω for a polynomial in P(Ω)). Let x := (x 1 , . . . , x n ), ω := (ω 1 , . . . , ω d ), α be a multi-index (size given by context), and c α be the polynomial coefficients. Let P s (S) for some s ∈ Z ≥0 , S ∈ {X, Ω} denote the subspace of P(S) consisting of polynomials of degree ≤ s, i.e. polynomials where the multi-indices of the monomial terms satisfy ||α|| 1 ≤ s. P SOS (X × Ω) refers to the space of polynomials on X × Ω that can be expressible as a sum-of-squares in x and ω jointly, and P s SOS (X × Ω) be the same space restricted to polynomials of degree ≤ s. Additionally, W ≽ 0 for a matrix W denotes that W is symmetric positive semidefinite (PSD). Finally, P(Ω) denotes the set of Lebesgue probability measures on Ω. For more details, see Appendix A.1.We present two formulations of the S-SOS hierarchy that are dual to each other in the sense of Fenchel duality [24, 25] . The primal problem seeks to find the tightest lower-bounding function and the dual problem seeks to find a minimizing probability distribution. Note that the "tightest lower bound" approach is dual to the "minimizing distribution" approach, otherwise known as a "joint and marginal" moment-based approach originally detailed in [19] .Consider a polynomial f (x, ω) : R n+d → R with x ∈ X ⊆ R n , ω ∈ Ω ⊆ R d equipped with a probability measure ν(ω). We interpret x as our optimization variables and ω as noise parameters, and seek a lower-bounding function c * (ω) such that f (x, ω) ≥ c * (ω) for all x, ω. In particular, we want the tightest lower bound c * (ω) = inf x∈X f (x, ω). Note that even when f (x, ω) is polynomial, the tightest lower bound c * (ω) can be non-polynomial. A simple example is the function f (x, ω) = (x -ω) 2 + (ωx) 2 , which has c * (ω) = inf x f (x, ω) = ω 4 /(1 + ω 2 ) (Appendix A.6.1).For us to select the "best" lower-bounding function, we want to maximize the expectation of the lower-bounding function c(ω) under ω ∼ ν(ω) while requiring f (x, ω) -c(ω) ≥ 0, giving us the following optimization problem over L 1 -integrable lower-bounding functions:EQUATIONs.t. f (x, ω) -c(ω) ≥ 0Even if we restricted c(ω) to be polynomial so that the residual f (x, ω) -c(ω) is also polynomial, we would still have a challenging nonconvex optimization problem over non-negative polynomials.In SOS optimization, we take a relaxation and require the residual to be SOS: f (x, ω) -c(ω) ∈ P SOS (X × Ω). Doing the SOS relaxation of the non-negative Equation (1) and restricting c(ω), i.e. f (x, ω) -c(ω) to polynomials of degree ≤ 2s gives us Equation ( 2), which we call the primal S-SOS degree-2s SDP:EQUATIONs.t. f (x, ω) -c(ω) = m s (x, ω) T W m s (x, ω)where m s (x, ω) is a basis function X × Ω → R a(n,d,s) containing monomial terms of degree ≤ s written as a column vector, and W ∈ R a(n,d,s)×a(n,d,s) a symmetric PSD matrix. Here, a(n, d, s) represents the dimension of the basis function, which depends on the degree s and on the dimensions n, d. For this formulation to find the best degree-2s approximation to the lowerbounding function, we require g(x, ω) = m s (x, ω) T W m s (x, ω) to span P 2s (X × Ω). Selecting all combinations of standard monomial terms of degree ≤ s suffices and results in a basis function with size a(n, d, s) = n+d+s s .The formal dual to Equation (1) (proof of duality in Appendix A.5.1) seeks to find a "minimizing distribution" µ(x, ω), i.e. a probability distribution that places weight on the minimizers of f (x, ω) subject to the constraint that the marginal µ X (ω) matches ν(ω):EQUATIONs.t. X dµ(x, ω) = µ X (ω) = ν(ω)where we have written P(X × Ω) as the space of joint probability distributions on X × Ω and µ X (ω) is the marginal of µ(x, ω) with respect to ω, obtained via disintegration.For the primal, we considered polynomials of degree ≤ 2s. We do the same here. The formal dual becomes a tractable SDP, where the objective turns into moment-minimization and the constraints become moment-matching. Following [3, 15] , let M ∈ R a(n,d,s)×a(n,d,s) be the symmetric PSD moment matrix with entries defined asM i,j = X×Ω m (i) s (x, ω)m (j)s (x, ω)dµ(x, ω) where m (i) s (x, ω) is the i-th element of the basis function m s . Let y ∈ R b(n,d,s) be the moment vector of independent moments that completely specifies M , e.g. in the case that we use all standard monomials of degree ≤ s and have a(n, d, s) = n+d+s s , then b(n, d, s) = n+d+2s. We write M (y) as the moment matrix that is formed from these independent moments. We havey α(i,j) = X×Ω m (i) s (x, ω)m (j)s (x, ω)dµ(x, ω) where the multi-index α(i, j) ∈ Z n+d ≥0 corresponds to the sum of the multi-indices corresponding to the i-th entry and the j-th entry of m s (x, ω).We write f (x, ω) in terms of the monomials f (x, ω) = ||α||1≤2s f α [x, ω] α , where [x, ω] is the concatenation of the n + d variables from x, ω and α ∈ Z n+d ≥0 is a multi-index. Note that every monomial [x, ω] α has a corresponding moment y α : [x, ω] α dµ(x, ω) = y α . We then observe that the integral in the objective reduces to a dot product between the coefficients of f and the moment vector:f (x, ω)dµ(x, ω) = α f α [x, ω] α dµ(x, ω) = α f α y αAfter converting the distribution-matching constraint µ X (ω) = ν(ω) in (3) into equality constraints on the moments of ω up to degree 2s, we obtain the following dual S-SOS degree-2s SDP:d * 2s = inf y∈R b(n,d,s) ||α||1≤2s f α y α (4) s.t. M (y) ≽ 0 y α = m α ∀ (α, m α ) ∈ M νWe write M ν as the set of (α, m α ) representing the moment-matching constraints on ω α up to degree-2s, i.e. we want to set X×Ω ω α dµ(x, ω) = Ω ω α dν(ω) = m α for all multi-indices α ∈ Z d ≥0 with ||α|| 1 ≤ 2s. There are d+2s 2s multi-indices α ∈ Z n+d ≥0 , ||α|| 1 ≤ 2s where only the d entries associated with ω are non-zero, and therefore the number of moment-matching constraints is|M ν | = d+2s2s . Note that the moment matrix M (y) ∈ R a(n,d,s)×a(n,d,s) is a symmetric PSD matrix and is the dual variable to the primal W . Observe also that we require the moments of ν(ω) of degree up to 2s to be bounded. ( 4) is often a more convenient form than (2), especially when working with additional equality or inequality constraints, as we will see in Section 3. For concrete examples of the primal and dual SDPs with explicit constraints, see Appendix A.3.In this section, we detail two ways of building a hierarchy, one based on the maximum degree of monomial terms in the basis function (Lasserre) and a novel one based on the maximum number of interactions occurring in the terms of the basis function (cluster basis). To define any SOS hierarchy, we first select a monomial basis. Some examples include the standard monomial basis x 1 , . . . , x n , trigonometric/Fourier 1-periodic monomial basis sin x 1 , cos x 1 , . . . , sin x n , cos x n ), or others. Using this basis, we write down a basis function m(x) which comprises some combinations of monomials. Squared linear combinations of the basis functions then span a SOS space of functions:H : {( i h i m i (x)) 2 }.In the Lasserre hierarchy, the basis function m s (x) is composed of all combinations of monomials up to degree s ∈ Z >0 and a given level of the hierarchy is set by the maximum degree s. The basis function consists of terms x α with α a multi-index and ||α|| 1 ≤ s. The degree-2s SOS function space parameterized by this basis function is that spanned by m s (x) T W m s (x) for PSD W , i.e. the functions that can result from squaring any linear combination of degree-s polynomials that can be generated from our basis m s (x). As we increase the degree s, our basis function gets larger and our S-SOS SDP objective values converge to the optimal value of the "tightest lower-bounding" problem Equation (1) [19] .In this section, we propose a cluster basis hierarchy, wherein we utilize possible spatial organization of the problem to sparsify the problem and reduce the size of the SDP that must be solved [26, 27] . The cluster basis is a physically motivated prior often used in statistical and condensed matter physics, where we assume that our degrees of freedom can be arrayed in space, with locally close variables interacting strongly (kept in the model) and globally separated variables interacting weakly (ignored). Moreover, one may also keep only the terms with interactions between a small number of degrees of freedom, such as considering only pairwise or triplet interactions between particles.In the cluster basis hierarchy, a given level of the hierarchy is specified by a 2-tuple (b, t), the desired body order b and the maximum degree of a single variable t. Body order denotes the maximum number of interacting variables in a given monomial term, e.g. x a i x b j x c k would have body order 3 and total degree a + b + c. The basis function m b,t consists of terms x α with α a multi-index, ||α|| 0 ≤ b (at most b interacting variables can occur in a single term), and ||α|| ∞ ≤ t (each variable can have up to degree t. The maximum degree of the basis function m b,t is then s = bt. If we are to compare m b,t from the cluster basis hierarchy with m s from the Lasserre hierarchy, we find that even when bt = s we still have strictly fewer terms, e.g. in the case where b = 2, t = 2, s = 4 we have m s containing terms of the form x 4 i but m b,t only has degree-4 terms of the x 2 i x 2 j . To expand on this, consider that in the standard Lasserre hierarchy we have m s (x) containing all monomials of degree ≤ s in n variables, or n+s s terms in total. In the proposed cluster basis hierarchy, m b,t (x)has b k=0 n k t kterms. This expression results from the need to sum over body orders k, considering that there are n k ways to choose k variables and that each selected variable has k possible degrees so there are t k ways to assign degrees ≤ t to these k variables. For fixed b, t, s = bt, and n ≫ b, t we note that m s (x) has O(n bt ) terms while m b,t (x) has O(n b ) terms so the size reduction factor in using the cluster basis asymptotically goes like n b(t-1) . As the number of variables n in the problem grow, we have asymptotic dominance in using the cluster basis to reduce the size of the SDP that must be solved. As bt → ∞ we might expect asymptotic convergence of the SDP hierarchy just like the standard Lasserre hierarchy, however, a full proof of convergence is out of scope for this paper. For further details, see discussion in Appendix A.7.4.As we increase the degree s (either s in the Lasserre hierarchy or b, t in the cluster basis hierarchy) we would expect the SDP objective values p * 2s (Equation ( 2)) to converge to the optimal value p * and the lower bounding function c * 2s (ω) to converge to the tightest lower bound c * (ω) = inf x f (x, ω). In this paper we refer to p * 2s → p * and d * 2s → d * interchangeably as strong duality occurs in practice despite being difficult to formally verify (Appendix A.4). This convergence is a common feature of SOS hierarchies.In this section we show that using polynomial c * 2s (ω) to approximate c * (ω) still allows for asymptotic convergence in L 1 as s → ∞.We specialize to the particular case of 1-periodic trigonometric polynomials f (x, ω), c(ω) on compact X = [0, 1] n and compact Ω ⊂ R d and prove asymptotic convergence of the degree-2s S-SOS hierarchy as s → ∞. Though seemingly restrictive, this choice allows us to leverage Fourier convergence results on a compact domain. For generic f (x, ω) where we are only interested in its behavior in some compact set (nearly all practical problems involve a restriction of domain), we may rescale the compact set to this 1-periodic compact domain and apply the results found here. To use other families of polynomials, we note that a substitution argument justifies the focus on trigonometric polynomials as any convergence result achieved for the trigonometric polynomial hierarchy lead directly to a matching result on regular polynomials (2.2 in [28] ). Previously, rates of 1/s 2 were attained for the standard Lasserre hierarchy (i.e. without external parameters ω) on compact domains with mild assumptions [29, 28, 9] . In our work, a slower convergence rate of ln s/s is achieved due to the introduction of external parameters, requiring us to approximate the lower-bounding function c * (ω) and thereby limiting the obtained convergence rate.2.3.2 ln s/s convergence using a polynomial approximation to c * (ω)We would like to bound the gap between the optimal lower bound c * (ω) = inf x∈X f (x, ω) and the lower bound c * 2s (ω) resulting from solving the degree-2s primal S-SOS SDP, i.e.EQUATIONTo that end, we need to understand the regularity of c * . Without further assumptions, we may assume c * to be Lipschitz continuous, per Proposition 2.1.With Equation ( 5) we may then integrate0 ≤ Ω inf x f (x, ω) -c * 2s (ω)dν(ω) ≤ |Ω|ε(f, s)where we control ε in terms of the degree s. If we can drive ϵ → 0 as s → ∞ then we are done. Proposition 2.1 (Theorem 2.1 in [30] ). Let g : X ×Y → R be polynomial. Then y → inf x∈X g(x, y) is Lipschitz continuous. Theorem 2.1 (Asymptotic convergence of S-SOS). Let f : [0, 1] n × Ω → R be a trigonometric polynomial of degree 2r, c * (ω) = inf x f (x, ω) the optimal lower bound as a function of ω, and ν any probability measure on compact Ω ⊂ R d . Let s refer to the degree of the basis in both x, ω terms and the degree of the lower-bounding polynomial c(ω), i.e. m s ([x, ω]) : R n+d → R a(n,d,s) is the full basis function of terms [x, ω] α with ||α|| 1 ≤ s and c(ω) only has terms ω α with ||α|| 1 ≤ s.Let p * 2s be the solution to the following S-SOS SDP (c.f. Equation (2)) with m s (x, ω) a spanning basis of trigonometric monomials with degree ≤ s:p * 2s = sup c∈P 2s (Ω),W ≽0 c(ω)dν(ω) s.t. f (x, ω) -c(ω) = m s (x, ω) T W m s (x, ω)Then there is a constant C > 0 depending only on ||f -f || F , ||c * -c * ||, r, Ω, n, d such that the following holds:Ω [c * (ω) -c * 2s (ω)] dν(ω) ≤ C ln s swhere f denotes the average value of the functionf over [0, 1] n , i.e. f = [0,1] n f (x)dx and ||f (x)|| F = x | f (x)| denotes the norm of the Fourier coefficients. Thus we have asymptotic convergence of the S-SOS SDP hierarchy to the optimal value p * of Equation (1) as we send s → ∞.Proof. The following is an outline of the proof. For complete details, including the full theorem and proof, please see Appendix A.5.2.The key argument used follows that of [22, 28, 23, 31, 32] . Take the lower-bounding minimizer of Equation (1) c * (ω) (not necessarily polynomial) and approximate it with a trigonometric polynomial c * a (ω) of degree s. We then pass the non-negative component f -c * a through an integral operator T that is built out of sums of squares of trigonometric polynomials (related to the Christoffel-Darboux and Jackson kernels) to obtain a strictly-positive SOS component:T h(x, ω) = X×Ω |q x (x -x)| 2 |q ω (ω -ω)| 2 h(x, ω)dxdω.q x (x), q ω (ω) can be chosen to be "kernel functions" of bounded degree ≤ s so that the output SOS function is of degree 2s. We can show the approximation error ||c * -c * a || is small, that the operator T exists and is close to the identity, and the deformation resulting from the SOS projection is small for sufficiently large degree s (degree of the approximating c * a (ω) and in the SOS kernel). We may conclude that the true strictly-positive part may be well-approximated by the SOS hierarchy and find asymptotic convergence along with a convergence rate in the degree s of the hierarchy. 
numerical experiments: 
 We present two numerical studies of S-SOS demonstrating its use in applications. The first study (Section 3.1) numerically tests how the optimal values of the SDP Equation (2) p * 2s converge to p * of the original primal Equation (1) as we increase the degree. The second study (Section 3.2) evaluates the performance of S-SOS for solution extraction and uncertainty quantification in various sensor network localization problems.As a simple illustration of S-SOS, we test it on the SOS functionf (x, ω) = (x -ω) 2 + (ωx) 2 (6) with x ∈ R, ω ∈ R. The lower bound c * (ω) = inf x f (x, ω) can be computed analytically as c * (ω) = ω 4 /(1 + ω 2 ).Assuming ω ∼ Uniform(-1, 1), we get that the objective value for the "tightest lower-bounding" primal problem Equation ( 1)is p * = 1 -1 ω 42(1+ω 2 ) dω = π 4 -2 3 ≈ 0.1187. For further details, see Appendix A.6.We are interested in studying the quantitative convergence of the S-SOS hierarchy numerically. The idea is to solve the primal (dual) degree-2s SDP to find the tightest polynomial lower bound (the minimizing probability distribution) for varying degrees s. As s gets larger, the basis function m s (x) gets larger and the objective value of the SDP Equation (2) p * 2s should converge to the theoretical optimal value p * .In Figure 1 we see very good agreement between p * and p * 2s with exponential convergence as s increases. This is much faster than the rate we found in Section 2.3.2, but agrees with the exponential convergence results from [6] achieved with local optimality assumptions. Due to the simplicity of (6), it is not surprising that we see much faster convergence. In fact, for most typical functions, we might expect convergence much faster than the worst-case rate. The tapering-off of the convergence rate is likely attributed to the numerical tolerance used in our solver (CVXPY/MOSEK), as we observed that increasing the tolerance shifts the best-achieved gap higher.Sensor network localization (SNL) is a common testbed for global optimization and SDP solvers due to the high sensitivity and ill-conditioning of the problem. In SNL, one seeks to recover the positions of N sensors X ∈ R N ×ℓ positioned in R ℓ given a set of noisy observations of pairwise distances d ij = ||x i -x j || between the sensors [15, 33] . To have a unique global minimum and remove symmetries, sensor-anchor distance observations are often added, where several sensors are anchored at known locations in the space. This can improve the conditioning of the problem, making it "easier" in some sense.We define a SNL problem instance with X ∈ [-1, 1] N ×ℓ as the ground-truth positions for S = {1, 2, . . . , N } sensors, A ∈ [-1, 1] K×ℓ as the ground-truth positions for A = {1, 2, . . . , K} anchors, D ss (r) = {d ij = ||x i -x j || : i, j ∈ S and d ij ≤ r} as the set of observed sensor-sensor distances and D sa (r) = {d ik = ||x i -a k || : i ∈ S, k ∈ A and d ik ≤ r} as the set of observed sensor-anchor distances, both of which depend on some sensing radius r. Writing x i , a k ∈ [-1, 1] ℓ as the unknown positions of the i-th sensor and the k-th anchor, we can write the potential function to be minimized as a polynomial:f (x, ω; X, A, r) = dij ∈Dss(r) (||x i -x j || 2 2 -d ij (ω) 2 ) 2sensor-sensor interactions+ d ik ∈Dsa(r) (||x i -a k || 2 2 -d ik (ω) 2 ) 2sensor-anchor interactions (7) The observed sensor-sensor and sensor-anchor distances d ij (ω), d ik (ω) can be perturbed arbitrarily, but in this paper we focus on linear uniform noise, i.e. for a subset of observed distances we haveEQUATION). Other noise types may be explored, including those including outliers, which may be a better fit for robust methods (Appendix A.7.2). Equation ( 7) contains soft penalty terms for sensor-sensor terms and sensor-anchor terms. We can see that this is a degree-4 polynomial in the standard monomial basis elements, and a global minimum of this function is achieved at f (X, 0 d ; X, A, r) = 0 (where the distances have not been perturbed by any noise). In general for non-zero ω (measuring distances under noise perturbations) we expect the function minimum to be > 0, as there may not exist a configuration of sensors X that is consistent with the observed noisy distances.We can also support equality constraints in our solution, in particular hard equality constraints on the positions of certain sensors relative to known anchors. This corresponds to removing all sensor-anchor soft penalty terms from the function and instead selecting N H < N sensors at random to exactly fix in known positions via equality constraints in the SDP. The SDP is still large but the effective number of variable sensors has been reduced toN ′ = N -N H .A given SNL problem type is specified by a spatial dimension ℓ, N sensors, K anchors, a sensing radius r ∈ (0, 2 √ ℓ), a noise type (linear), and anchor type (soft penalty or hard equality). Once these are specified, we generate a random problem instance by sampling X ∼ Uniform(-1, 1) n , A ∼ Uniform(-1, 1) d . The potential f (x, ω) for a given instance is formed (either with sensor-anchor terms or not, with terms kept based on some sensing radius r, and noise variables appropriately added).The number of anchors is chosen to be as few as possible so as to still enable exact localization, i.e. K = ℓ + 1 anchors for a SNL problem in ℓ spatial dimensions. The SDPs are formulated with the help of SymPy [34] and solved using CVXPY [35, 36] and Mosek [37] on a server with two Intel Xeon 6130 Gold processors (32 physical cores total) and 256GB of RAM. For an expanded discussion and further details, see Appendix A.7.The accuracy of the recovered solution is of primary interest, i.e. our primary evaluation metric should be the distance between our extracted sensor positions x and the ground-truth sensor positions X, i.e. dist(x, X). Because the S-SOS hierarchy recovers estimates of the sensor positions E[x i ] along with uncertainty estimates Var[x i ], we would like to measure the distance between our ground-truth positions X to our estimated distribution p(x) = N (E[x], Var[x]). The Mahalanobis distance δ M (Equation ( 8)) is a modified distance metric that accounts for the uncertainty [38] . We use this as our primary metric for sensor recovery accuracy.EQUATIONAs our baseline method, for each problem instance we apply a basic Monte Carlo method detailed in Algorithm 1 (Appendix A.7.3) where we sample ω ∼ ν(ω), use a local optimization solver to find x * (ω) = inf x f (x, ω), and use this to estimateE ω∼ν [x], Var ω∼ν [x].Note that though this non-SOS method achieves some estimate of the dual SDP objective f (x, ω)dµ(x, ω), it is not guaranteed to be a lower bound.Recovery accuracy. In Table 1 we see a comparison of the S-SOS method and the MCPO baseline.Each row corresponds to one SNL problem type, i.e. we fix the physical dimension ℓ, the number of anchors K = ℓ + 1, and select the sensing radius r and the noise scale ϵ. We then generate L = 20 random instances of each problem type, corresponding to a random realization of the ground-truth sensor and anchor configurationsX ∈ [-1, 1] N ×ℓ , A ∈ [-1, 1]K×ℓ , producing a f (x, ω) that we then solve the SDP for (in the case of S-SOS) or do pointwise optimizations for (in the case of MCPO). Each method outputs estimates for the sensor positions and uncertainty around it as aN (E[x], Cov[x]), which we then compute δ M for (see Equation ( 8)), treating each dimension as independent of each other (i.e. X as a flat vector). Each instance solve gives us one observation of δ M or each method, and we report the median and the ±1σ 34% values over the L = 20 instances we generate.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes]Justification: A full proof of our core result is provided in the supplementary material (Appendix A.5.2), along with an outline of the core ideas in the main body (Section 2.3.2).Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
discussion: 
 In this paper, we discuss the stochastic sum-of-squares (S-SOS) method to solve global polynomial optimization in the presence of noise, prove two asymptotic convergence results for polynomial f and compact Ω, and demonstrate its application to parametric polynomial minimization and uncertainty quantification along with a new cluster basis hierarchy that enables S-SOS to scale to larger problems.In our experiments, we specialized to sensor network localization and low-dimensional uniform random noise with small n, d. However, it is relatively straightforward to extend this method to support other noise types (such as Gaussian random variates without compact support, which we do in Appendix A.6.4) and support higher-dimensional noise with d ≫ 1.Scaling this method to larger problems n ≫ 1 is an open problem for all SOS-type methods. In this paper, we take the approach of sparsification, by making the cluster basis assumption to build up a block-sparse W . We anticipate that methods that leverage sparsity or other structure in f will be promising avenues of research, as well as approximate solving methods that avoid the explicit materialization of the matrices W, M . For example, we assume that the ground-truth polynomial possesses the block-sparse structure because our SDP explicitly requires the polynomial f (x, ω) to exactly decompose into some lower-bounding c(ω) and SOS f SOS (x, ω). Relaxing this exactdecomposition assumption and generalizing beyond polynomial f (x, ω), c(ω) may require novel approaches and would be an exciting area for future work. We also compare the number of elements in the full basis a f , the cluster basis a c , and the reduction multiple when using the cluster basis a f /a c . When passing to the cluster basis, a f /a c is how much the semidefinite matrix shrinks by. A Appendix / supplemental material A.1 NotationParameters Basis comparison M-distance (δ M ) ℓ r ϵ N H N C N a f a c a f /a c S-Let P(X) and P(Ω) denote the spaces of polynomials on X ⊆ R n and Ω ⊆ R d , respectively, where X and Ω are (not-necessarily compact) subsets of their respective ambient spaces R n and R d . Specifically, all polynomials of the forms below belong to their respective spaces:p(x) = α∈Z ≥0 c α x α ∈ P(X), p(ω) = α∈Z ≥0 c α ω α ∈ P(Ω)where x = (x 1 , . . . , x n ), ω = (ω 1 , . . . , ω d ), α is a multi-index for the respective spaces, and c α are the polynomial coefficients.Let P d (S) for some S ∈ {X, Ω} denote the subspace of P(S) consisting of polynomials of degree ≤ d, i.e. polynomials where the multi-indices of the monomial terms satisfy ||α|| 1 ≤ d. P SOS (X × Ω) refers to the space of polynomials on X × Ω that can be expressible as a sum-of-squares in x and ω jointly. Additionally, W ≽ 0 for a matrix W denotes that W is symmetric positive semidefinite (PSD). Finally, P(Ω) denotes the set of Lebesgue probability measures on Ω.The theoretical justification underlying the SDP relaxations in global optimization we use here derive from the Positivstellensätz (positivity certificate) of [7] , a representation theorem guaranteeing that strictly positive polynomials on certain sets admit sum-of-squares representations. Following this, [3, 10, 11] developed the Moment-SOS hierarchy, describing a hierarchy of primal-dual SDPs (each having fixed degree) of increasing size that provides a monotonic non-decreasing sequence of lower bounds.There is rich theory underlying the SOS hierarchy combining disparate results from algebraic geometry [14, 10, 11] , semidefinite programming [15, 18] , and complexity theory [16, 39] . The hierarchy exhibits finite convergence in particular cases where convexity and a strict local minimum are guaranteed [5] , otherwise converging asymptotically [6] . In practice, the hierarchy often does even better than these guarantees, converging exactly at c * s for some small s. The SOS hierarchy has found numerous applications in wide-ranging fields, including: reproducing certain results of perturbation theory and providing useful lower-bound certifications in quantum field theory and quantum chemistry [40, 41] , providing better provable guarantees in high-dimensional statistical problems [42, 43] , useful applications in the theory and practice of sensor network localization [15, 44] and in robust and stochastic optimization [45] .Due to the SDP relaxation, the SOS hierarchy is quite powerful. This flexibility comes at a cost, primarily in the form of computational complexity. The SDP prominently features a PSD matrix W ∈ R a(n,d,s)×a(n,d,s) with a(n, d, s) scaling as n+d+s s for n dimensions and maximum degree s. Without exploiting the structure of the polynomial, such as locality (coupled terms) or sparsity, solving the SDP using a standard interior point method becomes prohibitively expensive for moderate values of s or n. Work attempting to improve the scalability of the core ideas underlying the SOS hierarchy and the SDP method include [17, 18] .The S-SOS hierarchy we present in this paper as a solution to parametric polynomial optimization was presented originally by [19] as a "Joint + Marginal" approach. That work provides the same hierarchy of semidefinite relaxations where the sequence of optimal solutions converges to the moment vector of a probability measure encoding all information about the globally-optimal solutions x * (ω) = argmin x f (x, ω) and provides a proof that the dual problem (our primal) obtains a polynomial approximation to the optimal value function that converges almost-uniformly to c * (ω).Once a physical system or optimization problem is characterized, sensitivity analysis and uncertainty quantification seek to quantify how randomness or uncertainty in the inputs can affect the response. In our work, we have the parametric problem of minimizing a function f (x, ω) over x where ω parameterizes the function and is drawn from some noise distribution ν(ω).If only function evaluations f (x, ω) are allowed and no other information is known, Monte Carlo is often applied, where one draws ω k ∼ ν(ω) and solves many realizations of inf x f k (x) = f (x, ω k ) to approximately solve the following stochastic program:f * = inf x E ω∼ν [f (x, ω)]Standard Monte Carlo methods are ill-suited for integrating high-dimensional functions, so this method is computationally challenging in its own right. In addition, we have no guarantees on our result except that as we take the number of Monte Carlo iterates T → ∞ we converge to some unbiased estimate of E ω∼ν [f (x, ω)].Our approach to quantifying the uncertainty in optimal function value resulting from uncertainty in parameters ω is to find a deterministic lower-bounding c * (ω) which guarantees f (x, ω) ≥ c * (ω) no matter the realization of noise. This is reminiscent of the polynomial chaos expansion literature, wherein a system of some stochastic variables is expanded into a deterministic function of those stochastic variables, usually in some orthogonal polynomial basis [20, 21] .Example A.1. Let f (x, ω) be some polynomial of degree ≤ 2s written in the standard monomial basis, i.e.f (x, ω) = ||α||1≤2s f α x α = ||α||1≤2s f (α1,...,α n+d ) n i=1 x αi i d i=1 ω αn+i iLet m s (x, ω) ∈ R a(n,d,s) be the basis vector representing the full set of monomials in x, ω of degree ≤ s with a(n, d, s) = n+d+s s . For all α ∈ Z n+d ≥0 with ||α|| 1 ≤ 2s and α k = 0 for all k ∈ {1, . . . , n} (i.e. monomial terms containing only ω 1 , . . . , ω d ) we must have:X×Ω ω α dµ(x, ω) - Ω ω α dν(ω) = 0Explicitly, for µ to be a valid probability distribution we must have:X×Ω dµ(x, ω) -1 = M 0,0 -1 = y (1,0,...) -1 = 0 Suppose Ω = [-1, 1], ω ∼ Uniform(-1, 1) so that d = 1, ν(ω) = 1/2. We require: X×Ω ω α dµ(x, ω) = [-1,1] ω α dν(ω) =              1 α = 0 0 α = 1 1 3 α = 2 0 α = 3 1 5 α = 4To guarantee strong duality theoretically, we need a strictly feasible point in the interior (Slater's condition). For us, this is a consequence of Putinar's Positivstellensatz, if f (x, ω) admits a decomposition as f (x, ω) = c(ω) + g(x, ω) where g(x, ω) > 0 (i.e. is strictly positive), we have strong duality, i.e. p * = d * and p * 2s = d * 2s [3, 8] . However, it is difficult to verify the conditions analytically. In practice, strong duality is observed in most cases, so in this paper we refer to solving the primal and dual interchangeably, as p * 2s = d * 2s in all cases we encounter where a SDP solver returns a feasible point.Regular SOS Global polynomial optimization can be framed as the following lower-bound maximization problem where we need to check global non-negativity:EQUATIONs.t. f (x) -c ≥ 0 ∀xWhen we take the SOS relaxation of the non-negativity constraint in the primal, we now arrive at the SOS primal problem, where we require f (x) -c to be SOS which guarantees non-negativity but is a stronger condition than necessary:EQUATIONs.t. f (x) -c ∈ P SOS (X).The dual to Equation ( 9) is the following moment-minimization problem:EQUATIONwith dµ(x) = 1.Taking some spanning basis m s (x) : R n → R a(n,s) of monomials up to degree s, we have the moment matrix M ∈ R a(n,s)×a(n,s) :M i,j = m i (x)m j (x)dµ(x) = y αwhere we introduce a moment vector y whose elements correspond to the unique moments of the matrix M . Then we may write the degree-2s moment-minimization problem, which is now in a solvable numerical form:EQUATION)with M (y)1,1 = 1 M (y) ≽ 0where we write M (y) as the matrix formed by placing the moments from y into their appropriate places and we set the first element of m s (x) to be 1, hence M 1,1 = dν(x) = 1 is simply the normalization constraint. For further reading, see [15, 3] .Stochastic SOS Now let us lift this problem into the stochastic setting with parameters ω sampled from a given distribution ν, i.e. replacing x → (x, ω). We need to make some choice for the objective. The expectation of the lower bound under ν(ω) is a reasonable choice, i.e.but we could also make other choices, such as ones that encourage more robust lower bounds. In this paper however, we formulate the primal S-SOS as below (same as Equation ( 1)):EQUATIONs.t. f (x, ω) -c(ω) ≥ 0Note that if the ansatz space for the function c(ω) is general enough, the maximization of the curve c is equivalent to a pointwise maximization, i.e. we recover the best approximation for almost all ω.Then the dual problem has a very similar form to the non-stochastic case.Theorem A.2. The dual to Equation (13) is the following moment minimization where µ(x, ω) is a probability measure on X × Ω:inf µ∈P(X×Ω) f (x, ω)dµ(x, ω) with X×Ω ω α dµ(x, ω) = Ω ω α dν(ω) for all α ∈ N d .Remark A.3. Notice, that the condition X×Ω ω α dµ(x, ω) = Ω ω α dν(ω) implies that the first marginal of µ is the noise distribution ν. Let µ ω denote the disintegration of µ with respect to ν, [46] . Then the moment matching condition is equivalent to µ ω (X) = 1 for almost all ω and µ being a Young measure w.r.t. ν. The idea is that µ ω (x) is a minimizing density for every single configuration of ω.Proof. We use P ≥0 (X × Ω) to denote the space of non-negative polynomials on X × Ω. Given measure ν on Ω and polynomial function p :X × Ω → R consider sup γ∈L 1 (Ω,ν) q∈P ≥0 (X×Ω). Ω γ(ω)dν(ω) s.t p(x, ω) -γ(ω) = q(x, ω) This is equivalent to - inf γ∈L 1 (Ω,µ) q∈P ≥0 (X×Ω) f (γ, q) + g(γ, q) with f (γ, q) = - Ω γ(ω)dν(ω) and g(γ, q) = -χ {f -γ-q=0} = 0 if f -γ -q = 0 -∞ else ,i.e. g is the characteristic function enforcing non-negativity.Denote by h * the Legendre dual, i.e.h * (y) = supx ⟨x, y⟩ -h(x).Then by Rockafellar duality, [47, 24] , and noting that signed Borel measures B are the dual to continuous functions, the dual problem readssup Γ∈L ∞ (Ω,µ),µ∈B -f * (Γ, µ) -g * (-(Γ, ν))and we would havesup Γ∈L ∞ (Ω,µ),µ∈B -f * (Γ, µ) -g * (-(Γ, µ)) = - inf γ∈L 1 (Ω,µ) q∈P ≥0 (X×Ω) f (γ, q) + g(γ, q).The Legendre duals of f and g can be explicitly calculated asf * (Γ, µ) = 0 if Γ = -1 and µ ≤ 0 ∞ else and g * (Γ, µ) =    Ω×X f (x, ω)dµ(ω, x) if f -γ ∈ P ≥0 (X × Ω) and Γ(ω) = µ ω (X) ∞ else since f * (Γ, µ) = sup γ,q Ω γ(ω)Γ(ω)dν(ω) + Ω×X q(x, ω)dµ(x, ω) -f (γ, q) = sup γ,q Ω γ(ω)(Γ(ω) + 1)dν(ω) + Ω×X q(x, ω)dµ(x, ω) = 0 if Γ = -1 and µ ≤ 0 ∞ else and g * (Γ, µ) = sup γ,q Ω γ(ω)Γ(ω)dν(ω) + Ω×X q(x, ω)dµ(ω, x) + χ {f -γ-q=0} =    sup γ Ω γ(ω)Γ(ω)dν(ω) + Ω×X (f (x, ω) -γ(ω))dµ(ω, x) if f -γ =∈ P ≥0 (X × Ω) ∞ else =    sup γ Ω γ(ω)(Γ(ω) -µ ω (X))dν(ω) + Ω×X (f (x, ω)dµ(ω, x) if f -γ ∈ P ≥0 (X × Ω) ∞ else =    Ω×X (f (x, ω)dµ(ω, x) if f -γ ∈ P ≥0 (X × Ω) and Γ(ω) = µ ω (X) ∞ elseAltogether, we get-f * (Γ, µ) -g * (-Γ, -µ) =    Ω×X f (x, ω)dµ(ω, x) if µ ω (X) = 1 ∞ else.A.5.2 Convergence of S-SOS hierarchy c * (ω) is at best Lipschitz continuous By Proposition 2.1, we argue that c * (ω) = inf x f (x, ω) is Lipschitz continuous. One cannot expect much more as the following example shows:Example A.4. Consider g : R × R 2 → R defined by g(x, p, q) = (x 2 + px + q) 2 .Then we have for every(p, q) ∈ R 2 that inf x∈R g(x, p, q) = 0 if p 2 4 ≥ q ( p 2 4 -q) 2 else.Therefore, (p, q) → inf x∈R g(x, p, q) is once differentiable but not twice.Lemma A.5. Let Ω be a compact subset of R n and g : Ω → R be Lipschitz continuous. Then there exists a trigonometric polynomial g s of degree s and a constant C > 0 depending only on Ω and n such that g ≥ g s and ∥g -g s ∥ L 1 (Ω) ≤ 1 + ln(s) s C∥g∥ H 1 (Ω) .Proof. By Jackson's inequality for Lipschitz functions [48] we have the existence of a trigonometric polynomial g ′ of degree s with∥g -g ′ ∥ L 1 (Ω) ≤ C ′ s ∥g∥ H 1 (Ω)as well as∥g -g ′ ∥ L ∞ (Ω) ≤ L g ln(s) s .Then we define g s = g ′ -∥g -g ′ ∥ L ∞ (Ω) and hence g ≥ g s . Furthermore,∥g -g s ∥ L 1 (Ω) ≤ (C ′ + |Ω| ln s) s L g .Writing C(Ω) = max{C ′ , |Ω|} we have the desired form where |Ω| is the volume of Ω.Convergence at ln s/s rate using integral operator methodology Theorem A.1 (Asymptotic convergence of S-SOS). Let f : [0, 1] n × Ω → R be a trigonometric polynomial of degree 2r, c * (ω) = inf x f (x, ω) the optimal lower bound as a function of ω, and ν any probability measure on compact Ω ⊂ R d . Let s = (s x , s ω , s c ), referring separately to the degree of the basis in x terms, the degree of the basis in ω terms, and the degree of the lower-bounding polynomial c(ω).Let c * 2s (ω) be the lower bounding function obtained from the primal S-SOS SDP with m s (x, ω) a spanning basis of trigonometric monomials with degree ≤ s x in x terms and of degree ≤ s ω in ω terms:p * 2s = sup c∈P 2sc (Ω),W ≽0 c(ω)dν(ω) s.t. f (x, ω) -c(ω) = m s (x, ω) T W m s (x, ω)Then there is a constant C > 0 depending only on Ω, d, and n such that for all s ω , s x ≥ max{3r, 3s c } the following holds:Ω [c * (ω) -c * 2s (ω)] dν(ω) ≤ |Ω|ϵ(f, s) ε(f, s) ≤∥f -f ∥ F 1 -1 - 6r 2 s 2 ω -d 1 - 6r 2 s 2 x -n + ∥c * -c * ∥ F 1 -1 - 6r 2 s 2 ω -d + C (1 + ln(2s c )) 2s c .where f denotes the average value of the function f over ϵ(f, s) bounds the expected error, giving us asymptotic convergence as s = min(s x , s ω , s c ) → ∞.Note the first two terms give a O( 1 s 2 ) convergence rate. However, the overall error will be dominated by the degree of c(ω) (from the third term) hence our convergence rate is O( ln s s ).Proof of Theorem A.1. Let Ω ⊂ R d be compact and f : R n × Ω → R be a 1-periodic trigonometric polynomial (t.p.) of degree ≤ 2r. We then make Ω isomorphic to [0, 1] d and hereafter considerΩ = [0, 1] d and f : [0, 1] n × [0, 1] d → R. Let ε > 0 and b = ε 2 . Let the best lower bound be c * (ω) = inf x∈X f (x, ω).Proof outline. We split the error into two parts. First, we use the fact that there is a lower-bounding t.p. c * a of degree s c such that ∥c * -c * a ∥ ≤ C 1 + ln s c s c and c * ≥ c * a . This will provide us with a degree-s c t.p. approximation to the lower bounding function, which in general is only known to be Lipschitz continuous.Next, we show, that for any b > 0 there is a degree-2s SOS t.p. f SOS (x, ω) such thatf SOS = f -(c * a -b).We write s = (s x , s ω ) where s x , s ω denotes the respective max degrees in the variables x, ω. Once we have constructed this, we can compute f -f SOS = c * a -ε and since we know that f SOS ≥ 0 everywhere and c * a -ε is some degree-s c t.p. we have found a degree-s c lower-bounding t.p. The construction of this SOS t.p. adds another error term. If we can drive ε → 0 as s = min(s x , s ω , s c ) → ∞ then we are done.Proof continued. To that end, let c * a : Ω → R be the best degree-s c trigonometric approximation of c * with respect to L 1 such that c * ≥ c * a .By [30] , we know that c * is locally Lipschitz continuous with Lipschitz constant L c * and hence, by Lemma A.5 we get that there is C(Ω) > 0 such that∥c * -c * a ∥ L 1 Ω) ≤ C(Ω) 1 + ln s c s c L c * .Next we introduce c * 2s (ω) which is some degree-2s t.p. After an application of the triangle inequality and Cauchy-Schwarz on the integrated error term Ω |c * -c * 2s |dω we haveΩ inf x∈X f (x, ω) -c * 2s (ω) dω ≤ Ω |c * a (ω) -c * 2s (ω)|dω + |Ω|∥c * -c * a ∥ L 2 (Ω) Ω inf x∈X f (x, ω) -c * 2s (ω) dω ≤ Ω |c * a (ω) -c * 2s (ω)|dωgap between some SDP solution c * 2s (ω) and t.p. c * a (ω)+ C(Ω) 1 + ln s c s c L c * approx. error of L-contin. fn.Now we want to show that for any ε > 0 we can construct a degree-2s SOS trigonometric polynomial Observe that by assumption f -c * a + b is a t.p. in (x, ω) where f is degree-2r and c * a is degree s c ≥ 2r. Denote by (f -f a * + b) ω its coefficients w.r.t the ω basis. Note that the coefficients are functions in x. Following the integral operator proof methodology in [6] , define the integral operator T to bef SOS (x, ω) such that f SOS = f -c * a + b. with b = ε/2 and s = (s x , s ω ) > r. We can then set f -f SOS = c * a -b = c *T h(x, ω) = X×Ω |q ω (ω -ω)| 2 |q x (x -x)| 2 h(x, ω)dxdω,where q ω is a trigonometric polynomial in ω of degree ≤ s ω and q x is a trigonometric polynomial in x of degree ≤ s x . The intuition is that this integral operator explicitly builds a SOS function of degrees (s x , s ω ) out of any non-negative function h by hitting it against the kernels q x , q ω .We want to find a positive function h : X × Ω → R such thatT h = f -c * a + b.In frequency space, the Fourier transform turns a convolution into pointwise multiplication so we have:T h(x, ω) = qω * qω (ω) • qx * qx (x) • ĥ(x, ω).In the Fourier domain it is easy to write down the coefficients of ĥ:ĥ(x, ω) =    0 if ∥x, ω∥ ∞ > max{2r, 2s c } f (x, ω) -ĉ * a (ω)1 x=0 + b1 x=0 1 ω=0 qω * qω (ω) • qx * qx (x)otherwise.Computing T h -h gives:f (x, ω) -c * a (ω) + b -h(x, ω) = ω,x f (x, ω) 1 - 1 qω * qω (ω) • qx * qx (x) exp(2iπ ωT ω) exp(2iπx T x) + ω (b1 ω=0 -c * a ) 1 - 1 qω * qω (ω) exp(2iπ ωT ω)and thus after requiring qω * qω (0) = qx * qx (0) = 1 we have:max x,ω |f (x, ω) -c * a (ω) + b -h(x, ω)| ≤∥f -f ∥ F max ω̸ =0 max x̸ =0 1 - 1 qω * qω (ω) • qx * qx (x) + max ω̸ =0 ∥c * a -c * a ∥ F 1 - 1 qω * qω (ω).As a reminder, because c * ≥ c * a everywhere we have f -c a ≥ f -c * ≥ 0 or f -c * a + b > 0, since b = ε/2 > 0. Since T h = f -c * a + b > 0 and it is a SOS, we need to guarantee h > 0. If max x,ω |f (x, ω) -f a * (ω) + b -h(x, ω)| ≤ b then max x,ω |T h -h| < b.Since T h ≥ b and b > 0 we haveh = T h + h -T h ≥ T h -∥h -T h∥ ∞ ≥ b -b ≥ 0 and hence h > 0 if we ensure max x,ω |T h -h| ≤ b. Now let us show that max x,ω |f (x, ω) -c * a (ω) + b -h(x, ω)| ≤ b can be ensured if s = (s x , s ω ) is large enough.Using the same kernel and bounds as in [6] , we choose for z ∈ {x, ω} the triangular kernel such thatqz (ẑ) = 1 - 6r 2 z 2 d + d i=1 1 - |ẑ i | s x,ω + .Note that (x) + = max(x, 0). Then we havemax x |f (x, ω) -c * a (ω) + b -h(x, ω)| ≤∥f -f ∥ F max ω,x 1 - 1 qω * qω (ω) • qx * qx (x) + ∥c * a -c * a ∥ F max ω 1 - 1 qω * qω (ω) ≤∥f -f ∥ F 1 -1 - 6r 2 s 2 ω -d 1 - 6r 2 s 2 x -n + ∥c * a -c * a ∥ F 1 -1 - 6 2 s 2 ω -dTherefore, by choosing s ω and s x large enough such that∥f -f ∥ F 1 -1 - 6r 2 s 2 ω -d 1 - 6r 2 s 2 x -n + ∥c * a -c * a ∥ F 1 -1 - 6 2 s 2 ω -d ≤ b = ε 2we have h ≥ 0 and thus T h is SOS. By design we havec * a -c * 2s ≤ b and thus Ω |c * a -c * 2s |dω ≤ ε 2 .RecallingΩ inf x∈X f (x, ω) -c * 2s (ω) dω ≤ Ω |c * a (ω) -c * 2s (ω)|dωgap between some SDP solution c * 2s (ω) and t.p. c * a (ω)+ C(Ω) 1 + ln s c s c L c * approx. error of L-contin. fn.we can additionally choose s c large enough to guaranteeEQUATIONand then we are done.Setting s x , s ω , s c = s and sending s → ∞ we have asymptotic behavior of the final error expression:Ω inf x∈X f (x, ω) -c * 2s (ω) dω ≤ C 1 1 s 2 + C 2 1 s + C 3 ln s s = O ln s s with the constants C 1 , C 2 , C 3 depending on r, n, d, ∥f -f ∥ F , ∥c a -c * a ∥ F , Ω and L c * .Convergence at 1/s rate using piecewise-constant approximation to c * (ω) Prior work [6] achieves 1/s 2 convergence for the regular SOS hierarchy without further assumptions. In our work, we could only achieve ln s/s for S-SOS due to the need to first approximate the tightest lower-bounding function c * (ω) with a polynomial approximation, which converges at a slower rate.To accelerate the convergence rate, we want to control the regularity of c * (ω). We can achieve 1/s by approximating the c * (ω) pointwise instead of using a smooth parameterized polynomial. By constructing a domain decomposition of Ω and finding a SOS approximation in x for each domain, we can stitch these together to build a piecewise-constant approximation to the lower-bounding function c * . This lets us leverage the 1/s 2 rate for the standard SOS hierarchy, however note that we lose the guarantee that the lower-bounding function c * s (ω) is a true lower bound everywhere. Another disadvantage of this approach is that the computational complexity scales exponentially in d as we are required to solve for the SOS lower bound with degree s pointwise at O(s d p ) points. For Ω = [0, 1] d ⊂ R d we achieve the following: Proposition A.6. Let Ω = [0, 1] d and f : [0, 1] n → R be a trigonometric polynomial of degree 2r. Let {ω i } be grid points where s p is the number of equally-spaced points along each dimension for a total of (s p ) d points in Ω. Denote by c * s (ω i ) the best SOS approximation of degree s of x → f (x, ω i ) at each grid point ω i and definec * s (ω) = s d p i=1 c * s (ω i )1 ω∈B∞(ωi,ϵ)where B ∞ (ω i , ϵ) denotes the L ∞ ball of radius ϵ centered at ω i with ϵ = (2s p ) -1 so that the balls are non-overlapping. Then we have for some constant C ′ depending only onmax ωi ||f (ω i , •) - f (ω i , •)|| F , r, n, d, s p : Ω [c * (ω) -c * s (ω)]dω ≤ max ωi ∥f (ω i , •) -f (ω i , •)∥ F 1 -1 - 6r 2 s 2 -n + C √ d s p ≤ C ′ s .Proof of Proposition A.6. Let c * a (ω) : [0, 1] d → R be the best piecewise-constant approximation of the true lower-bounding function c * (ω) = inf x f (x, ω) on equidistant grid-points, i.e. s p evenlyspaced points in each dimension for a total of s d p points. Now also define a piecewise-constant approximation using the best pointwise SOS lower bounds:c * s (ω) = s d p i=1 c * s (ω i )1 ω∈B∞(ωi,ϵ)where c * s (ω i ) is the best lower bound (resulting from regular SOS) of degree s of x → f (x, ω i ) and B ∞ (ω i , ϵ) denotes the L ∞ ball of radius ϵ centered at ω i with ϵ = (2s p ) -1 so that the balls are non-overlapping.By [6] we have that c * a (ω i ) -c * s (ω i ) can be bounded bymax ωi ∥f (ω i , •) -f (ω i , •)∥ F 1 -1 - 6r 2 s 2 -n.Then we have:Ω [c * (ω) -c * s (ω)]dω ≤ s d p i=1 |c * a (ω i ) -c * s (ω i )|(2ϵ) d + ∥c * -c * a ∥ L 1 (Ω) .Using the same bound we get for the first term from the proof of Theorem A.1, we can reduce the first term to a O(1/s 2 ) dependence and we use the theorem on the L 1 convergence of piecewise-constant approximation to 1-periodic trigonometric polynomials from [48] for the second:Ω [c * (ω)s -c * s (ω)]dω ≤ max ωi ∥f (•, ω i ) -f (•, ω i )∥ F 1 -1 - 6r 2 s 2 -n |Ω| + C √ d s pIt is straightforward to extend this to arbitrary compact Ω ⊂ R d , as one can construct a smooth homeomorphism between Ω and [0, 1] d and therefore we can achieve the same convergence rates on general compact Ω as we can for [0, 1] d .We provide a simple application of S-SOS to a simple quadratic potential that admits a closed-form solution so as to demonstrate its usage and limitations.A.6.1 Analytic solution for the lower bounding function c * (ω) with ω ∼ Uniform(-1, 1) Let x ∈ R and ω ∼ Uniform(-1, 1). Suppose that we havef (x, ω) = (x -ω) 2 + (ωx) 2In this case we may explicitly evaluate the exact minimum function c * (ω) = inf x f (x; ω). Note thatf (x; ω) = x 2 -2ωx + ω 2 + ω 2 x 2Explicitly evaluating the zeros of the first derivative we have∂ x f (x; ω) = 2x * -2ω + 2ω 2 x * = 0 x * (1 + ω 2 ) = ω x * = ω 1 + ω 2 and, thus, c * (ω) = inf x f (x; ω) = ω 4 1 + ω 2 .Note that despite f (x, ω) being a simple degree-2 SOS polynomial, the tightest lower-bound c * (ω) = inf x f (x, ω) is explicitly not polynomial. However, it is algebraic, as it is defined implicitly as the root of the polynomial equation Observe that the tightest lower-bounding function c * (ω) is not polynomial even in this simple setting. However, we can relax the problem to trying to find c 2s ∈ P 2s (Ω) to obtain a weaker bound withc * (ω)(1 + ω 2 ) -ω 4 = 0inf x f (x, ω) = c * (ω) ≥ c 2s (ω).We now proceed with formulating and solving the degree-2s primal S-SOS SDP (Equation ( 2)). We assume that c 2s (ω) is parameterized by a polynomial of degree ≤ 2s in ω. Observe that this class of functions is not large enough to contain the true function c * (ω).We choose s ∈ {2, 4} and use the standard monomial basis in x, ω, we have the feature maps m 2 (x, ω) : R 2 → R 6 and m 4 (x, ω) : R 2 → R 15 , since there are n+s s unique monomials of up to degree-s in n variables. These assumptions together enable us to explicitly write a SOS SDP in terms of coefficient matching. Note that we must assume some noise distribution ν(ω). For this section, we present results assuming ω ∼ Uniform(-1, 1). We solve the resulting SDP in CVXPY using Legendre quadrature with k = 5 zeroes on [-1, 1] to evaluate the objective c(ω)dν(ω). In fact, k sample points suffice to exactly integrate polynomials of degree ≤ 2k -1.We solve the SDP for two different levels of the hierarchy, s = 2 and s = 4 (producing lower-bound polynomials of degree 4 and 8 respectively), and plot the lower bound functions c 2s (ω) vs the true lower bound c * (ω) = ω 4 /(1 + ω 2 ) as well as the optimality gap to the true lower bound in Fig. 2 .To solve the S-SOS SDP in practice, we must choose a maximum degree 2s for the SOS function m 2 (x, ω) T W m 2 (x, ω) and the lower-bounding function c(ω), which are both restricted to be polynomials. Indeed, a larger s not only increases the dimension of our basis function m s (x, ω) but also the complexity of the resulting SDP. We would expect that d * 2s → d * as s → ∞, i.e. the optimal value of the degree-2s S-SOS SDP (Equation ( 4)) converges to that of the "minimizing distribution" optimization problem (Equation ( 3)).In particular, note that in the standard SOS hierarchy we typically find finite convergence (exact agreement at some degree 2s * < ∞). However, in S-SOS, we thus far have only a guarantee of asymptotic convergence, as each finite-degree S-SOS SDP solves for a polynomial approximation to the optimal lower bound c * (ω) = inf x∈X f (x, ω). In Figure 1 , we illustrate the primal S-SOS SDP objective values for a given level of the hierarchy (a chosen degree s for the basis m s (x, ω)) and their convergence towards the optimal objective valuep * 2s = sup c∈P 2s (Ω) c(ω)dν(ω) with f (x, ω) -c(ω) ∈ P 2s SOS (X × Ω)c * (ω)dν(ω) = π 4 - 2 3 ≈ 0.1187for the simple quadratic potential, assuming ν(ω) = 1 2 with ω ∼ Uniform(-1, 1). We note that in the log-linear plot (right) we have a "hinge"-type curve, with a linear decay (in logspace) and then flattening completely. This suggests perhaps that in realistic scenarios the degree needed to achieve a close approximation is very low, lower than suggested by our bounds. The flattening that occurs here is likely due to the numerical tolerance used in our solver (CVXPY/MOSEK), as increasing the tolerance also increases the asymptotic gap and decreases the degree at which the gap flattens out.In the previous two sections, we assumed that ω ∼ Uniform(-1, 1). This enabled us to solve the primal exactly using Legendre quadrature of polynomials. Note that in Figure 1 we see that the lower-bounding c * 2 (ω), c * 4 (ω) for ω ∼ Uniform(-1, 1) is a smooth polynomial that has curvature (i.e. sign matching that of the true minimum). This is actually not guaranteed, as we will see shortly.In Figure 3 , we present the lower-bounding functions c * 4 (ω) achieved by degree-4 S-SOS by solving the dual for ω ∼ Normal(0, σ 2 ) for varying widths σ. We can see that for small σ ≪ 1, the primal solution only cares about the lower-bound accuracy within a small region of ω = 0, and the lower-bounding curve fails to "generalize" effectively outside the region of consideration.Recall the form of f (x, ω):f (x, ω; X, A, r) = dij ∈Dss(r) (||x i -x j || 2 2 -d ij (ω) 2 ) 2sensor-sensor interactions+ d ik ∈Dsa(r) (||x i -a k || 2 2 -d ik (ω) 2 ) 2sensor-anchor interactions Note that the function f (x, ω) is exactly a degree-4 SOS polynomial, so it suffices to choose the degree-2 monomial basis containing a = N ℓ+d+2 2 elements as m 2 (x, ω) : R N ℓ+d → R a . That is, we have N sensor positions in ℓ spatial dimensions and d parameters for a total of N ℓ + d variables.Let the moment matrix be M ∈ R a×a with elements defined asM i,j := m (i) 2 (x, ω)m (j) 2 (x, ω)dµ(x, ω)for i, j ∈ {1, . . . , a}, which fully specifies the minimizing distribution µ(x, ω) as in Equation (4).Our SDP is then of the formd * 4 = inf y α f α y α s.t. M (y) ≽ 0 y α = m α ∀ (α, m α ) ∈ M ν y α = y * α ∀ (α, y * α )∈ H where y α = m α corresponds to the moment-matching constraints of Equation ( 4) and y α = y * α correspond to any possible hard equality constraints required to set the exact position (and uncertainty) of a sensorE[x i ] = x * i , E[x 2 i ] -E[x i ] 2 =0 for all ω. M ν represents the d+2s 2s moment-matching constraints necessary for all moments w.r.t. ω and H represents the 2ℓn constraints needed to set the exact positions of n known sensor positions in R ℓ (i.e. 1 constraint per sensor and dimension, 2 each for mean and variance).In this paper we focus on the linear uniform noise case, as it is a more accurate reflection of measurement noise in true SNL problems. Special robust estimation approaches may be needed to properly handle the outlier noise case.• Linear uniform noise: for a subset of edges we writed ij,k (ω) = d * ij + ϵω k , ω k ∼ Uniform(-1, 1), and ϵ ≥ 0 some noise scale we set. The same random variate ω k may perturb any number of edges. Otherwise the observed distances are the true distances.• Outlier uniform noise: for a subset of edges we ignore any information in the actual measurementd ij,k = ω k , ω k ∼ Uniform(0, 2 √ ℓ)where ℓ is the physical dimension of the problem, i.e. x i ∈ R ℓ .Here we explicitly formulate MCPO and S-SOS as algorithms. Let X = R n , Ω = R d and use the standard monomial basis. We write z = [x 1 , . . . , x n , ω 1 , . . . , ω d ]. Our objective is to approximate c * (ω) = inf x f (x, ω) for all ω, with a view towards maximizing c * (ω)dν(ω) for ω sampled from some probability density ν(ω).MCPO (Algorithm 1) simply samples ω t and finds a set of tuples (x * (ω t ), ω t ) where the optimal minimizer (x * (ω t , ω t ) is computed using a local optimization scheme (we use BFGS). S-SOS (Algorithm 2) via solving the dual (Equation ( 4)) is also detailed below.Recall from Section 2.2.2 that we defined the cluster basis hierarchy using body order b and maximum degree per variable t. In this section, we review the additional modifications needed to scale S-SOS for SNL.In SNL, f (x, ω) is by design a degree s = 4 polynomial in z = [x, ω], with interactions of body order b = 2 (due to the (x i , x j ) interactions) and maximum individual variable degree t = 4. Written this way, we want to only consider monomial terms [x, ω] α with ||α|| 1 ≤ s, ||α|| ∞ ≤ 4, and ||α|| 0 ≤ 2.To sparsify our problem, we start with some k-clustering (k clusters, mutually-exclusive) of the sensor set C = {C 1 , . . . , C k }. This clustering can be considered as leveraging some kind of "coarse" Algorithm 1 Monte Carlo Point Optimization (MCPO) 1: Input: Function f (x; ω), sampler for distribution ν(ω), number of samples T 2: Output: Approximate integral Î, empirical distribution p D (x), empirical mean µ, empirical covariance Σ 3: for t = 1 to T do 4:Sample ω t ∼ ν(ω) 5:Find minimizer x t = min x f (x; ω t ) using BFGS 6: end for7: Estimate integral Î ≈ 1 T T t=1 f (x t ; ω t ) 8: Construct empirical distribution p D (x) = 1 T T t=1 δ(x -x t ) 9: Calculate empirical mean μ = 1 T T t=1 x t and covariance Σ = 1 T -1 T t=1 (x t -μ)(x t -μ) T .Algorithm 2 Stochastic Sum-of-squares (S-SOS), Dual formulation ) , function f (x; ω) : R n+d → R represented as a dictionary mapping multi-indexα ∈ Z n+d ≥0→ coefficient f α , probability density function for ν(ω) with known moments ω α dν(ω) < ∞ ∀ ||α|| 1 ≤ 2s, any hard equality constraints where we want to set x k = x * k for some k ∈ K. 2: Let i1, i2, i4 be the lexicographically-ordered arraysEQUATION≥0 which correspond to the arrays of multi-indices for all degree-1, degree-s, and degree-2s monomials in the variables z. ) as a vector of variables to be estimated, corresponding to the vector of independent moments that fully specifies M . 5: Add M ≽ 0 constraint. 6: for i in length(i2) do 7:for j in length(i2) do ▷ Require M to be formed from the elements of y.Computeα ij = i2[i] + i2[j]as the multi-index corresponding to the sum of the multiindices i2[i], i2[j].Add constraint M i,j = y αij . 10: end for 11: end for 12: for each row α in i4 do ▷ Require y α moments to equal the known moments of ω α .13:if n i=1 α i = 0 then 14: Add constraint y α = z α dν(ω) = ω α[-d:] dν(ω).15:end if 16: end for 17: for k in K do ▷ Handle any hard equality constraints in our variables x.Form multi-index α 1 ∈ Z n+d ≥0 where the entry for x k is set to 1 and everything else is zero.19:Form multi-index α 2 ∈ Z n+d ≥0 where the entry for x 2 k is set to 1 and everything else is zero.Add constraint information about which sensors are close to each other. For example, just looking at the polynomial f (x, ω) enables us to see which sensors (i, j) must be interacting.y α1 = x * k . ▷ E[x k ] = x * k . 21: Add constraint y α2 = (x * k ) 2 . ▷ Var[x k ] = E[x 2 k ] -E[x k ] 2 = 0.Assume that there is some a priori clustering given to us. We denote x (i) as the subset of the variables restricted to the cluster C i , i.e. x (i) = {x j : j ∈ C i }. Moreover, let G = (V, E) be a graph where the vertices V = {1, . . . , k} correspond to the k clusters and the edges E = {(i, j) : i, j ∈ V } correspond to known cluster-cluster interactions.The SOS part of the function f (x) may then be approximated as the sum of dense intra-cluster interactions and sparse inter-cluster interactions, where the cluster-cluster interactions are given exactly by edges in the graph G:m s (x) T W m s (x) ≈ i∈V m s (x (i) ) T W (i) m s (x (i) ) + (i,j)∈E m s (x (i) ) T W (i,j) m s (x (j) )where W (k) are symmetric PSD matrices and W (i,j) are rectangular matrices where we require W (i,j) = (W (j,i) ) T . m s (x) for x ∈ R n here behaves as before and denotes the basis function generated by all n+s s combinations of monomials with degree ≤ s. Notice that this is a strict reduction from the standard Lasserre hierarchy at the same degree s, since in general the standard basis m s (x) on the full variable set will contain terms that mix variables from two different clusters that may not have an edge connecting them.Efficiency gains in the SDP solve occur when we constrain certain of the off-diagonal W (i,j) blocks to be zero, i.e. the graph G is sparse in cluster-cluster interactions. As we can see from the block decomposition written above, this resembles block sparsity on the matrix W . We may interpret the above scheme as having a hierarchical structure out to depth 2, where we have dense interactions at the lowest level and sparse interactions aggregating them. In full generality, the resulting hierarchical sparsity in W may be interpreted as generating a chordal W , which is known to admit certain speed-ups in SDP solvers [26] .When attempting to solve an SNL problem in the cluster basis instead of the full basis, we need to throw away terms in the potential f (x, ω) that correspond to cross-terms that are "ignored" by the particular cluster basis we chose. The resulting polynomial f (x, ω) has fewer terms and produces a cluster basis SDP that is easier to solve, but generally less accurate due to the sparser connectivity.In particular, for the rows in Table 1 that have N C > 1, we do a N C -means clustering of the ground-truth sensor positions and use those sensor labels to create our partitioning of the sensors. We connect every cluster using plus-one c i , c i+1 (including the wrap-around one) connections, so that the cluster-cluster connectivity graph has N C edges. We then use this information to throw out observed distances from the set D ss and from the full basis function m 2 (x, ω). See our code for complete details.The sensor-anchor terms in Equation ( 7) are added to make the problem easier, because by adding them now each sensor no longer needs to rely only on a local neighborhood of sensors to localize itself, but can also use its position relative to some known anchor. When we remove them entirely, we need to incorporate hard equality constraints between certain sensors and known "anchor" positions. This fixes certain known sensors but lets every other sensor be unrooted, defined only relative to other sensors (and potentially an anchor if it is within the sensing radius).To deal with the equality constraints where we set the exact position of a sensor x i = x * i , we solve the dual Equation (4) and implement them as equality constraints on the moment matrix, i.e. for the basis element m 2 (x, ω) i = x i we may set E[x i ] -x * i = M 0,i -x * i = 0. Note that we also need to set Var(x i ) = 0 so for m 2 (x, ω) j = x 2 i we add the equality constraint Var(x i ) = E[x 2 i ] -E[x i ] 2 = M 0,j -M 2 0,i = 0.Once the dual SDP has been solved, we extract the moment matrix M and can easily recover the point and uncertainty estimates for the sensor positions E[x], Var[x] by inspecting the appropriate entries M 0,i corresponding to m 2 (x, ω) i = x i and M 0,j corresponding to m 2 (x, ω) j = x 2 i . A.7.7 Impact of using MCPO with varying numbers of samples TIn Figure 4 we can see how δ M varies as we scale the number of samples T used in the MCPO estimate of the empirical mean/covariance of the recovered solutions. In this particular example, the runtime of the S-SOS estimate was 0.3 seconds, comparing to 30 seconds for the T = 300 MCPO point. Despite taking 100x longer, the MCPO solution recovery still dramatically underperforms S-SOS in δ M . This reflects the poor performance of local optimization methods vs. a global optimization method (when it is available).The largest 2D SNL experiment we could run had N = 15 sensors, N C = 9 clusters, and d = 9 noise parameters. This generated N ℓ + d = 39 variables and 820 basis elements in the naive m 2 (x, ω) construction, which was reduced to 317 after our application of the cluster basis, giving us W, M ∈ R 317×317 . A single solve in CVXPY (MOSEK) took 30 minutes on our workstation (2x Intel Xeon 6130 Gold and 256GB of RAM). We attempted a run with N = 20 sensors and N C = 9 clusters and d = 9 noise parameters, but the process failed due to OOM constraints. Thus, we report the largest experiment that succeeded.As discussed in the main text, our work provides both a theoretical convergence rate as well as numerical understanding of the stochastic sum-of-squares (S-SOS). Notably, our convergence rate is bottlenecked on the approximation quality of c * (ω), only achieving ln s/s where regular SOS results achieve 1/s 2 . In practice, we see much faster convergence, which may be a clue that other more restrictive assumptions or more powerful methods may enable one to find better convergence results.In addition, our numerical work focuses on the sensor network localization problem. Due to the time complexity of SDPs, the maximum number of sensors we could solve in any one instance is limited to N = 15. This seems much smaller than the state-of-the-art, i.e. [49] where the authors propose an algebraic reduction of the noiseless SNL problem so that they can analytically obtain the range of the PSD matrix. This simplifies the SDP dramatically, however in the noisy setting (where observed sensor-sensor distances are perturbed with observation noise) this approach is unusable, requiring significant modification. his reduction is used to solve SNL problems of 10k-100k sensors, but they develop a highly specialized algorithm that does not use any SDP solvers.Note that noiseless SNL simplifies the problem. In the noiseless setting, localizing even a small number of sensors near an anchor will propagate the correctly localized positions to sensors nearby.As such, if one can find small groups of sensors that can be well-oriented with respect to each other, one can pursue localization of these groups in parallel and then orient them globally at the very end.In the noisy setting, any error in localization can easily propagate to the whole instance, dramatically reducing the size of problems that can be effectively solved [50, 51] .Thus, we observe that uncertainty makes the SNL problem considerably more difficult than the noiseless setting. S-SOS handles this naturally while also solving the problem for its global optimum. We get a solution for the sensor positions for every possible configuration of noise via the probability distribution µ(x, ω). This additional complexity places further limits on the size of the problems that can be solved with S-SOS.• While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes] Justification: All information regarding experiments is disclosed, with core details in the main draft and full experimental details in the appendix. In addition, we provide code in our submission (and will publicly release it upon paper decision) to reproduce all of the figures, tables, and other results we report.• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example (a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.(b) If the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully. (c) If the contribution is a new model (e.g., a large language model), then there should either be a way to access this model for reproducing the results or a way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 
open access to data and code: 
 Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes] Justification: Our code has been provided with our submission and a cleaned-up version will be released publicly once the preprint is public.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: All necessary details to replicate results are summarized in the core body and explicated in full in supplementary material (Appendix A.7).Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: Where necessary, we report confidence intervals on our solution recovery metrics for SNL which are computed over multiple runs. In particular, we report robust mean and 1-sigma robust std estimates and details can be found at Section 3.2.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: Specific details regarding the compute resources used are detailed in Section 3, as well a section on scalability in Appendix A.7.8.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes]Justification: All authors have reviewed the NeurIPS Code of Ethics and certify that the research presented conforms with all aspects of it.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA] Justification: As foundational research on algorithms and optimization, we do not see substantial societal impact positive or negative tied to possible applications of this research.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: Our contributions are primarily algorithmic and do not have a high risk for misuse.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets 
 
 Appendix:  


Decision: This paper considers an important extension of sum-of-squares for parametric polynomial optimization. In particular, the authors consider its stochastic counterpart and provide quantitative convergence results of the hierarchy as the degree of the polynomial increases. An application to sensor network localization problems is shown.

My own assessment is that the paper is a valuable and unique contribution the NeurIPS community, focusing on core stochastic optimization theory within the important SOS framework. There were some concerns about the size of the numerical experiment, but to the AC, the theoretical findings moved the quality needle across the bar for acceptance.

Reviews: review: paper_summary: paper_summary: The Sum of Squares (SOS) technique is well-known in polynomial optimization studies. This paper studies its variant, in which the function to be lower bounded has random parameters drawn from some probability distribution. The main contribution is a cluster-based SDP hierarchy for the Stochastic-SOS (S-SOS) method and a convergence proof for this and the standard Lasserre. Experiments demonstrate the effectiveness of the new hierarchy in the sensor network localization problem (SNL).
strength_weakness: strength_weakness: - It would be better if the authors presented the assumptions in a subsection and explained their restrictiveness. For example, Theorem 2.1 assumes the objective polynomial to be trigonometric. I wonder how challenging it is to extend the results to more general polynomial families. 

- The paper uses only two paragraphs to introduce and illustrate the cluster-basis hierarchy, a vital part of its contributions. A more comprehensive discussion of the hierarchy could enhance the paper's content and readability. Besides, I would suggest the authors polish the main paper more.

- Novelties in the proof techniques would often attract independent interests. I wonder if the convergence proof is established by applying only the standard techniques and existing results.

- The experiments seem to contain, at most, $N=15$ sensors, which are suitable for synthetic validation but may not be sufficient to demonstrate the technique's practical advantage.
questions: questions: Please check the comments on the paper's strengths and weaknesses. The authors should also feel free to point out any inaccuracy or misunderstanding. Thanks.
limitations: limitations: It may be beneficial if the authors could explain more about the paper's limitations in theory (restrictiveness of the assumptions) and practice (practicality of the experimental settings). But I do not find the current presentation concerning. Similarly, I do not think the paper would have any negative societal impact.
paper_summary: paper_summary: In their study, the authors investigate parametric polynomial optimization where the function to be minimized is \( f(x, \omega) \). Here, \( x \) represents the decision variable, and \( \omega \) signifies a noise parameter. The primary goal is to approximate the best lower bound, \( c^*(\omega) = \inf_x f(x, \omega) \), for each value of \( \omega \). To manage this setup, the Sum-Of-Squares (SOS) hierarchy is adapted, an approach first introduced by Lasserre in his seminal ``Joint and Marginal'' work in the late 2000s.

The paper elaborates on the derivation of this hierarchy and its dual, as well as discussing the rates of convergence under suitable assumptions about the optimal solution \( c^* \), and includes some applications towards the end of the paper. From a technical standpoint, the work seems correct and presents a rigorous approach to handling stochastic variables in polynomial optimization. However, it does not significantly deviate from established methodologies, which may limit its appeal in terms of novelty.

Regarding its relevance to the NeurIPS audience, while the modified hierarchy certainly adds practical value, the paper does not address applications that align closely with the core interests of the community. The experiments, focusing on sensor network localization, seem peripheral to the main areas of interest at NeurIPS, which typically centers around more direct applications to machine learning and artificial intelligence technologies.

In conclusion, although the paper is technically proficient and might captivate a niche audience interested in theoretical optimization, it appears to fall short of the high innovation standards and relevance to ML typically expected for NeurIPS publications.
strength_weakness: strength_weakness: -- Relevance 
-- Novelty
questions: questions: In your paper, you discuss the adaptation of the Sum-Of-Squares (SOS) hierarchy to handle stochastic parameters in polynomial optimization, which is a significant theoretical advancement. However, the practical applications presented, such as sensor network localization, seem somewhat tangential to the core interests of the NeurIPS community, which often focuses on direct applications in machine learning and artificial intelligence.

Could you elaborate on how the methodologies developed in your study could be applied to more central problems in machine learning? Additionally, are there potential modifications or extensions to the SOS hierarchy that could make it more directly applicable to common challenges in neural network training or optimization under uncertainty?
limitations: limitations: NA
paper_summary: paper_summary: This work proposes a new sum of squares (SOS) based relaxation for stochastic polynomial optimization,called Stochastic SOS (S-SOS). For a class of problems where the cost function $f$ can depend on variables $x$ as well as randomly drawn parameters $\omega$, the authors propose a hierarchy of relaxations whose solutions yield strict lower bounds on the optimal value of $f$. For a degree $s$ relaxation, it is shown that the obtained lower bound is $\log(s)/s$ and $1/s^2$ for different parameterizations of the parametric lower bound function $c(\omega)$. Finally, the proposed methods are applied to sensor network localization problems, where S-SOS outperforms the stated baseline.
strength_weakness: strength_weakness: 1. In equation (6), I believe the polynomial in question should be $f(x,\omega) = (x-\omega)^2 + (\omega x)^2$. Otherwise, the correct solution should be $c^*(\omega) = -\omega^4/(1-\omega^2)$.
2. Perhaps in the related work section, a little space could be devoted to other approaches to polynomial optimization, such as the LP/SDP relaxations constructed using the Polya and Handelman theorems, and the LP/SOCP relaxations obtained by using DSOS and SDSOS polynomials.
questions: questions: 1. Are there problems for which the exact lower bound can be achieved in a "single step"? That is, suppose the polynomial cost and constraint functions are of degree $s$, are there problems for which a degree $s$ relaxation suffices to obtain the true minimum? Do such problems arise in practical circumstances?

2. While the DSOS/SDSOS hierarchies do not, in general, converge to global optima, have the authors considered using the associated LP relaxations, particularly in the context of the sensor localization problem? This is salient since the LP/SOCP relaxations can be used to solve much larger problems than the SDP solutions. 

3. Is there any intuition as to how the result stated in Prop. 2.2. can be extended to the multivariate case?
limitations: limitations: The limitations of the work are discussed, though they aren't stated in a single section. It would be helpful to add such a section, perhaps in the appendix.
paper_summary: paper_summary: The manuscript considers variants of stochastic-SOS hierarchy for parametric polynomial optimization (POP), which had been previously considered using similar (joint+marginal) methods (https://arxiv.org/abs/0905.2497). Unfortunately the (joint+marginal) methods did not work very well, perhaps because parametric POP is known to have a very complicated structure (https://doi.org/10.1287/moor.2021.0097). While the numerical experiments are not exactly conclusive, the present manuscript provides guarantees on the rate of convergence, which improve upon the joint+marginal method somewhat. The "cluster basis" variant of the hierarchy could be of independent interest.
strength_weakness: strength_weakness: -- The numerical experiments are poorly designed. The sensor network localization problem is known to be reducible (https://arxiv.org/abs/1002.0013) to trivial sizes. Even early papers (https://link.springer.com/article/10.1007/s10107-009-0338-x) considered n=10000 senosrs. The 2010 paper tests on instances with n = 20000 to 100000 sensors, while the present authors consider instances on less than n=15 sensors. This is 4-5 orders of magnitude less than the state of the art. 

-- It is not clear how the objective of the heuristically extracted solution (A.7.6) differs from the objective function of the relaxation.
questions: questions: Could you explain the application of the joint+marginal method as a starting point? 

Could you comment on the rates of convergence of the  joint+marginal method?  

Could you demonstrate the behaviour of your S-SOS on the more complicated behaviours from https://doi.org/10.1287/moor.2021.0097? Cf. Definition 3.9 (Irregular accumulation point), but also Definition 3.7 (Discontinuous non-isolated multiple point) and Definition 3.6 (Discontinuous isolated multiple point).
limitations: limitations: The discussion of the limitations is misleading. It claims that all SDP-based algorithms fail to scale, while there clearly are (https://arxiv.org/abs/1002.0013) SDP-based algorithms that scale well.
paper_summary: paper_summary: This paper is concerned with the following natural generalization of polynomial optizimation: 
Given a distribution $\nu(\omega)$ and a polynomial $f(x,\omega)$, find the tightest (in expectation) function $c(\omega)$ that lower bounds $f(x,\omega)$  everywhere.

Given the computational hardness of the problem, the authors focus on the sum-of-squares relaxation (of a given degree $s$).

For trigonometric polynomials, they show a $O(\frac{\log s}{s})$ convergence of the error of the relaxation w.r.t. the optimum (here we hide factors depending on the input in the $O$ notation). Under stronger assumption on the optimal solution they further achieve fast error convergence.
Finally, they complement these results with experimental evaluations (in which they apply certain SDP sparsification techniques to improve the scalability).
strength_weakness: strength_weakness: In my opinion, the main weakness of the paper lies in the fact that it  fails to explain the  *core* ideas behind its results. There is not a careful discussion of the theorems, their limitations and the insight required to prove them. This makes the paper easy to follow but also quite frustrating to read.  I believe significant improvements can be made in this sense. Such improvements  could greatly improve the quality of the manuscript.

The theoretical analysis result is limited to trigonometric polynomials. This is an important restriction, but a discussion addressing this topic, possible extensions and related obstacles is not present.

The cluster basis hierarchy (henceforth CBH) used to scale up the SDP  is non-trivial and interesting. However it is never discussed in-depth and no theoretical result is provided. I understand that convergence results for CBH may be hard to prove, but a more comprehensive discussion seems needed.
questions: questions: What can the authors say for other classes of polynomials?
limitations: limitations: I encourage the authors to further discuss the limitations of their theoretical results.

Comments: title: Reply to the "Huh?" by Reviewer Dyr8
Thank you!
Huh?
Response to Author Rebuttal
comment: We thank reviewer Dyr8 for their remark. In particular, we appreciate the external citation they flagged.

Indeed, our comment was not precise and fails to be true in general. To be clearer, we were referring to the unconstrained case $$\min_{x\in\mathbb R^n} p(x)$$ for $p$ being already in sum of squares form. Hence, the observations of https://arxiv.org/abs/2403.08329 do not quite fit to the setting we wanted to stress.

We hoped for exact convergence with the same degree in the unconstrained setting. This will be true for example in the cases where the objective assumes $0$ as minimum. We would like to point out that the latter is the case for many physical objectives that are in the focus of our research and, thus, we made the claim. One can also observe that if $p(x)-p_{min}$ allows for an SOS representation it will be of same degree as the original objective in the unconstrained case since no cancellation effects can be observed.

So the remaining question is if we can shift SOS polynomials by constants and stay SOS. This, indeed, seems not to be true in general as the example $$p(x,y,z) = x^2y^2z^2+(x^2-1)^2+(y^2-1)^2+(z^2-1)^2$$ potentially can show (numerically verified). We would like to follow up on this question.

As a general remark why assuming $\min p =0$ is of interest, we would like to point out that eventually we do not know the exact form of the objective but are given only function values at samples. However, one can assume a physical model that in many cases would allow an exact SOS recovery in one step as the minimum will be zero.

Furthermore, observe that in the parametric case the situation is more complex.  Here we have to approximate the lower bounding function by polynomials, which in general is a severe restriction and, thus, we cannot expect to obtain finite convergence in any relevant case. Additionally, the degree we use for the SOS hierarchies is lower bounded by the degree we use to approximate the lower bounding function and the question of finite convergence seem less relevant.
Many thanks for having worked out the example. I see the distinction between the separated and non-separated parametric problems, where Bellon considered the former and you are consider the latter. Thank you!
The authors write:

**We observe that in the regular SOS hierarchy, if one starts with a degree-2s SOS polynomial then you will find exact convergence when using the hierarchy at degree-2s.**

I would suggest that the authors review https://arxiv.org/abs/2403.08329, which analyzes the following: 
\[
\begin{array}{rl}
\min_{x \in R} & x \\
\mathrm{s.t.} &1-x^2 \geq 0 \textrm{ and }
 & x+(1-\varepsilon)x^2 \geq 0 ,\\
\end{array}
\]
parametrized by a scalar $\varepsilon \in [0,1]$. There, the convergence of the moment-SOS hierarchy is finite, but arbitrarily slow as $\varepsilon$ goes to zero.
I thank the authors for their thoughtful rebuttal. At this point, I'm happy to maintain my positive score for this work.

**Edit:**

I would like to elaborate on my justifications for keeping the score as is (weak accept): the work, while interesting, (generally) technically correct, and novel, has a few drawbacks that prevent me from granting a higher score. First, the scale of experiments is far lower than previous work, as pointed out by Reviewer Dyr8. Second, the lack of a multivariate bound (which I believe nearly all problems, including the experiments stated here are), makes the analysis a little handwavy (the dependence on the square root of the dimension in a possible multivariate version of Prop. 2.2., as pointed out by the authors, weakens the result slightly from an applicability standpoint).