S-SOS: Stochastic Sum-Of-Squares for Parametric Polynomial Optimization
Authors: Richard Licheng Zhu
Mathias Oster
Yuehaw Khoo
Venue: NeurIPS 2024 poster (2024)
Keywords: convex optimization
semidefinite programming
sum-of-squares
uncertainty quantification
polynomial optimization
Primary area: optimization
TL;DR: We discuss a stochastic sum-of-squares algorithm for global polynomial optimization and use it for uncertainty quantification on sensor network localization problems.

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.

S-SOS: Stochastic Sun-Of-Squares for Parametric Polynomial Optimization
Richard Zhu
Department of Computational and Applied Mathematics
University of Chicago
Chicago, IL 60637
Mathias Oster
Institute of Geometry and Practical Mathematics
RWTH Aachen
Aachen, Germany
osterOigpm. wwwh--achen.cooo de
Yuehaw Khoo
Department of Statistics University of Chicago Chicago, IL 60637 yuehaw. kahbav.k edu
Abstract
Global polynomial optimization is an important tool across applied mathematics, with many applications in operations research, engineering, and the physical sci- ences. In various settings, the polynomials depend external parameters that may be random. We discuss stochastic sum-of-squares (S-SOS) algorithm based the sum-of-squares hierarchy that constructs series of semidefinite programs to jointly find strict lower bounds on the global minimum and extract 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 em- ploying n-body priors from condensed matter physics to induce sparsity, we can use S-SOS produce solutions and uncertainty intervals for sensor network lo- calization problems containing up to 40 variables and semidefinite matrix sizes surpassing 800 800.
Introduction
Many effective nonlinear and nonconvex optimization techniques local information to identify local minima. But often the case that we want to find global optima. of-uurrsss (SOS)
optimization is powerful and general technique in this setting.
The core idea is follows: suppose we are given polynomials ........mm. f where each function on 1 H and we seek to determine the minimum value of f the closed set S: S {o e 1n (x) 2 V = m. Our optimization problem is then to find infzer [f()) S).
"Corresponding Author
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

An equivalent formulation is to find the largest constant E R (i.e. the tightest lower bound) that can be subtracted from f such that f - 2 0 over the set S. This reduction converts a polynomial optimization problem over semialgebraic set to the problem of checking polynomial non-negativity. This problem NP-hard in general [I] therefore one instead resorts to checking if f - c is a sum-of-squares (SOS) function, e.g. in the unconstrained setting where S = I one seeks to find some polynomials hy : R" + IR such that f -c = >x h If such a decomposition can be found, then we have an easily-checkable certification that f - 2 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 non-negatyv problem.
Notably, if we restrict the fy to have maximum degree d, the search for degree- 2d SOS decomposi- tion of 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. given degree corresponds to particular level of the hierarchy. Solving this SIDP produces a lower bound cæ which has been proven to converge to the true global minimum e* = inf, f(r) as increases, with finite convergence (ca = at finite d) for functions with second-order local optimality conditions [ 6
and asymptotic
convergence with milder assumptions thanks to representation theorems for positive polynomials
from real algebraic geometry 7 8 (]. Further work has elucidated both theoretical implications
7 IO 12 [3 and useful applications of SOS to disparate fields [14 15 16 5 6
18 (see
further discussion in Appendix (A.2).
Motivated by the sum-of-squares certification for lower bound e on a function f(z), we generalize
to the case where the function to be minimized has additional parameters, i.e. [xx, w) where r are variables and w are parameters drawn from some probability distribution w v(w). We
seek a function e(w) that is the tightest lower bound to f(r, w) everywhere: f(a,, w) 2 c(w) with
c(w) - inf, (x, w). This setting was originally presented in [19] as a "Joint and Marginal" approach to parametric polynomial optimization. With the view that w v(w) and seeking to parameterize the minimizers x (w) = argmin,, (xx, we are reminded of some of the prior work in polynomial
chaos, where 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 oq-Squs (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 We then detail how different hierarchies can be constructed (Section 1.2.1]. Finally, in Section 2.3 (complete proof in Appendix A.3.2) we specialize to compact X x 52 and outline the proof for a theorem on quantitative convergence (the gap between the optimal values of the degree-2e 2s S-SOS SDP and the "tightest boundinn" optimization problem goes + 0 as s - oo) of the S-SOS hierarchy for trigonometric polynomials on [0, 1]" x [0, 1] following the kernel formalism of
(Section 2.1
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 x Ω = []1 x [-1, 1j4, 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 simple Monte Carlo-based method.
2 Stochastic Sum-of-squares (S-SOS)
2.0.1
Notation
Let P(S) be the space of polynomials on S. where S E (X, X C R2 and Ω C Rd, respectively, where X and Ω are noot-necessarlly compact) subsets of their respective ambient spaces R and Rd. A polynomial in P(X) can be written as p(x) = E.ezg, cax" E P(X) (substituting n +
d, + w. X + for a polynomial in P(22). Let x == (x)...., kn), == (wi....,wa)) Wa), a be multi-index (size given by context), and ca be the polynomial coefficients. Let (S) for some s E Zzo. E (X denote the subspace of P(S) consisting of polynomials of degree S 8, i.e. polynomials where the multi-indice of the monomial terms satisfy ||a|]1 Psos (X x ②) refers
2

---

to the space of polynomials on X x Ω that can be expressible as a sum-of-squares in r and w jointly, and Poss (-X 52) be the same space restricted to polynomials of degree I s. Additionally. W > for matrix W denotes that W is symmetric positive semidefinite (PSD). Finally, P(12) denotes the set of Lebesgue probability measures on f2 For more details, see Appendix A.1
2.1 Formulation of S-SOS hierarchy
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 find the tightest lower-bounding function and the dual problem seeks to find minimizing probability distribution. Note that the "tightest lower bound" approach is dual the "minimizing distribution" approach, otherwise known as "joint and marginal" moment- based approach originally detailed in [19].
2.1.1 Primal S-SOS: The tightest Nower-boundien function
Consider a polynomial f(z, w) Rn+d + I with I X S R",WE Q C I2 equipped with our optimization variables and as noise parameters, and seek bower-boundngg function e* (w) such that far, w) 2 (w) for all æ, w. particular, we want the tightest lower bound c" (w) = infaex f(r, w) Note that even when f(r, w) is polynomial, the tightest lower bound (w) can be non-polynn non-polynomial. A simple example is the function f(x, w) =
probability measure vww). We interpret
(x: -))² + (wx)) which has c* (w) inf, (zz w) = w /11 + w²) (Appendix A
For to select the "best" lower-bounding function, want maximize the expectation of the swwer-bounding function c(w) under w v(w) while requiring f(z.w) c(w) 2 0, giving us the following optimization problem over L -integrable lowrrbbouddig functions:
/ ((w)du(w)
pz = sup CEL CEL
s.t. f(z,w) c(() 2
Even if we restricted c(w) to be polynomial so that the residual f(x,, ((),,w) c(w) is also polynomial, would still have a challenging nonconvex optimization problem over non-negative polynomials. In SOS optimization, we take relaxation and require the residual to be SOS: f(r,w) - c(w) E Psos(( x Q). Doing the SOS relaxation of the non-negative Equation [1] and restricting c(w), i.e.
f(z, w) - c(w) to polynomials of degree < 2s gives us Equation Q which we call the primal S-SOS degree-2s SDP:
/ (cw))a((w)
P26
sup
(2)
s.t. f(x,w) - c(w) = m, W.r. (B))
Wm, (2,w)
where m,(x,w) w) is a basis function X x Ω Ra(n,ds) containing monomial terms of degree S written as column vector, and W E pa(ndda)xa(nnd..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 lower- bounding function, we require g(2, w) = m m2( WW m, (z. w) to span p2*(X x Q). Selecting all combinations of standard monomial terms of degree I suffices and results in basis function with size a(n, d, s) (n+d+.).
2.1.2
Dual S-SOS: minimizing distribution
The formal dual to Equation IT) (proof of duality in Appendix 1 distribution" ()),w), i.e a probability distribution that places weight the minimizers of f subject to the constraint that the marginal μux (w) matches vw))
seeks to find a "minimizing
f f(a,w)dp(a,
d'
inf (XP(XXX)
(3)
s.l. [ du(z,w) (x(w) = n(w)
3

---

where we have written P(X x Q2) as the space of joint probability distributions on X x Q2 and μax (w) the marginal of μ(x,w) w) with respect to w. obtained via disintegration.
For the primal, we considered polynomials of degree I 28. We do the same here. The formal dual becomes tractable SDP, where the objective turns into moment-minnimization and the con- straints become moment-matching Following B 151, let M Ra(ndd)) xa(n.d.)) be the symmetric PSD moment matrix with entries defined as Ma.5 = Jxxa m (x, (r. ( w)dd((a w) where
m (x,w) w) is the i-th element of the basis function m,. Let E (M(a,,)) be the moment vec- tor of independent moments that completely specifies M, c.g. in the case that use all stan- dard monomials of degree s and have a(n, d, (**d+*), then b(n, d, (**)++2*). We write MVy) the moment matrix that is formed from these independent moments. We have Va(s.)) = -xxxo x m (Ir, w)mm (x,))d/( ()) where the multi-index a(i,j) E Z corresponds to the sum of the multi-indicee corresponding to the i-th entry and the j-th entry of m,(xww).
We write f(x, w) in terms of the monomials f(r, = Elaath fa (z, w|" where [æ, is the concatenation of the n + variables from x,w and a E Zata is multi-index. Note that every monomial [xr, w/ has corresponding moment Yo: [ de, w) Yo. We then observe that the integral in the objective reduces to dot product between the coefficients of ʃ and the moment vector:
[ [le,w))
/ Σ (a
Σ faJk
After converting the distributionn-atchii constraint px (x(w) v(w) in 3 into equality constraints on the moments of w up to degree 2s, we obtain the following dual S-SOS degree-2s SDP:
dz
inf VER
Σ
faya
(4)
||al|:
s.t.
M(y) >
Ya A (a, ma) Mr
We write Mr as the set of (a,, m.) representing the moment- matthnt- constraints on w up to dsgree-2s, i.e. we want to Sxxo Booewaoe dld)(w Sn dddww) m. for all multi-indices
2 with ||l|| S 2s. There are (d42.) multi-indices E Zn0d /|a||| 2s where only the d entries associated with w are non-zero, and therefore the number of momentt-mtchhine constraints is |M, = ( (d=22).. Note that the moment matrix M (yy E [))(n,ds)xa(ndd4, is a symmetric PSD matrix and is the dual variable to the primal W. Observe also that we require the moments of vww) of degree up 2s to be bounded. 4 is often 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
2.2 Variations
In this section, we detail two ways of building 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 Ji..... Zn. trigonometric/Fourier 1 1-periodic monomial basis sin Z1, COS 20 ... sin Fn, COS x,), 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 SOS space of functions: H (EE, (x)
2.2.1 Standard Lasserre hierarchy
In the Lasserre hierarchy, the basis function mx(x)) composed of combinations of monomials up to degree E Z,50 and given level of the hierarchy is set by the maximum degree s. The basis function consists of terms xº with a a multi-index and ||a|]] S 8. The degree-2s SOS function space parameterized by this basis function is that spanned by mx(x)" Wm.(z) for PSD W, i.c. the

---

functions that can result from squaring any linear combination of degree-s polynomials that can be generated from our basis (.(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 T 19]
2.2.2 Cluster basis hierarchy
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 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 triplet interactions between particles.
In the cluster basis hierarchy, a given level of the hierarchy is specified by 2-tuple (b,t), the desired body order b and the maximum degree of single variable t. Body order denotes the maximum number of interacting variables in a given monomial term, e.g. 2 vg would have body order 3 and total degree a + + c. The basis function ma, consists of terms x² with a a multi-index, ||a/lo s (at most b interacting variables can occur in a single term), and ||a|| I (each variable can have up to degree The maximum degree of the basis function m/, is then 8 = bl. If we are to compare mo, from the cluster basis hierarchy with m, from the Lasserre hierarchy, we find that even when bt = we still have strictly fewer terms, c.g. in the case where b = 2, 2, s 4 we have m, containing terms of the form 2 but m,, only has degrec-4 terms of the x².
To expand on this, consider that in the standard Lasserre hierarchy we have mz(x) containing all monomials of degree < in n variables, or (*++)) terms in total. In the proposed cluster basis hierarchy, (x)) has
Arm=
terms. This expression results from the need to sum over body orders k considering that there are (2) ways to choose variables and that each selected variable has l possible degrees so there are t ways to assign degrees I to these & variables. For fixed b, t, s = bt, and n >> b, t we note that (z(x) has Onrht) terms while (z) has O(n") terms the size reduction factor in using the cluster basis asymptotically goes like (t--1) As the number of variables 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 bl + do we might expect asy 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
2.3 Convergence of S-SOS
As we increase the degree s (either in the Lasserre hierarchy or b, f in the cluster basis hierarchy) would expect the SDP objective values p2. (Equation (2) to converge to the optimal value p" and the lower bounding function G2 (w) to converge to the tightest lower bound e" (w) = inf, / (x, In this paper we refer to přz + p² and dž, 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 (w) to approximate c* (w) still allows for asymptotic convergence in L as cc.
2.3.1 Overview of result
We specialize to the particular case of 1 1-periodic trigonometric polynomials f(x w), o(w) on compact X [0, 11mh and compact C d and prove asymptotic convergence of the degree-2s S-SOS hierarchy as + co. Though seemingly restrictive, this choice allows us to leverage Fourier convergence results on a compact domain. For generic [ar, w) where we are only interested in its behavior in some compact set (nearly all practical problems involve a restriction of domain), may
5

---

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 substitution argument justifies the focus on trigonometric polynomials as any convergence result achieved for the trigonometrie polynomial hierarchy lead directly to matching result on regular polynomials (2.2 in (28). Previously, rates of 1 1/8 were attained for the standard Lasserre hierarchy (i.e. without external parameters w) on compact domains with mild assumptions [29 28 2 In our work, a slower convergence rate of In s/s is achieved due to the introduction of external parameters, requiring to approximate the lower-bounding function e² (w) and thereby limiting the obtained convergence rate.
2.3.2 In s/s convergence using polynomial approximation to c" (w)
We would like to bound the gap between the optimal lower bound c (w) in [(x, w) and the lower bound c2, (w) resulting from solving the degree-2s primal S-SOS SDP, i.c.
0 S c² (w) - (w) S e(f, V E
(5)
To that end, we need understand the regularity of e* Without further assumptions, we may assume c" to be Lipschitz continuous, per Proposition 2.1
With Equation (5) we may then integrate
0 $ J inf f(z,w) caa (w)dv(w) s)
where we control E in terms of the degree s. If we can drive e - 0 as s + oo then we are done. Proposition 2.1 (Theorem 2.1 in (3)). Let X X - IR be polynomial. Then + inf g y) Lipschitz continuous.
Theorem 2.1 (Asymptotic convergence of S-SOS). Let f [0, Ω R be trigonometric polynomial of degree 2r, c² ("()) inn f(x the optimal lower bound as a function of w. and v any probability measure on compact c Rd. Let s refer to the degree of the basis in both x, w terms and the degree of the lower-bounding polynomial c(w), Le. ma(1,, w|) Rn+dd -> Ra(n.d,)) is the full basis function of terms w) with ||a|]] I s and c(w) only has terms w with [|al|] I 8.
Let p2, be the solution to the following S-SOS SDP (c.f. Equation (2) with m,(I.w) spanning basis of trigonometric monomials with degree S s:
/ c(w)du(w)
sup CEP
f(r,w) c(w) m,(x,w)) Wm,
Then there constant > depending only on |lf fl1r, 2 1l, n, such that the following holds:
[ (c( du(w) s c.ms
where f denotes the average value of the function f over [0, i.e. f= - f(x)dx and
[If(=)]l== Ee f@)) 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 (I) as we send s + oc.
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
32 Take the lower-bounding minimizer of Equation D e* (w) (not necessarily polynomial) and approximate it with trigonometric polynomial (w) of degree s. We then pass the non-negative component f - o through an integral operator T that is built out of sums of squares of trigonometric polynomials (related to the Christoffel-Darbonn and Jackson kernels) to obtain a srrictly-positioo SOS component:
28 23
Th(r, w) = fr00 11-(x - 1))
gg (w - :) 2h(,, wdrdd

---

(x),, qu (w) can be chosen to be "kernel functions" of bounded degree I so that the output SOS function is of degree 2s. We can show the approximation error ||e* - is small, that the operator T exists and is close to the identity, and the deformation resulting from the SOS projection is small for sufficienlly large degree s (degree of the approximating e (w) and in the SOS kernel). We may conclude that the true strictly-posee positive part may be well-approximated by the SOS hierarchy and find asymptotic convergence along with a convergence rate in the degree 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 p1, converge to p" of the original primal Equation ( 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.
3.1 Simple quadratic SOS function
As simple illustration of S-SOS, we test it on the SOS function
f(r, w) (x w)² (wx)2
(6)
with x E R.w E R. The lower bound c ("(w) = in fxx w) can be computed analytically as (w) /1 w²). Assuming w Uniform (-1, 1), get that the objective value for the
"tightest lower-bounding" primal problem Equation is p" = [ d For further details, see Appendix A.6.
R 0.1187.
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 As gets larger, the basis function m. (x) gets larger and the objective value of the SDP Equation p2s should converge to the theoretical optimal value p"
In Figure we see very good agreement between pz and p2, 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 expect convergence much faster than the worst-case rate. The tapering-ofi of the convergence rate is likely attributed to the numerical tolerance used in solver (CVXPY/MOSEK), we observed that increasing the tolerance shifts the best-achieved gap higher.
much faster convergence. In fact, for most typical functions,
might
3.2 Sensor network localization
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 E &xx positioned in R2 given set of noisy observations of pairwise distances dıj = ||ا - æy between the sensors 15 33 To have a unique global minimum and remove symmetries, censor-anchor distance observations anchored at known locations in the space. This can improve the conditioning of the problem, making it "easier" in some sense.
often added, where several sensors
3.2.1 Definitions
We define SNL problem instance with [|1, as the ground-truth positions for {1, 2, .......... N} sensors, A (-1, 1]xxxt as the ground-truth positions for A 21,2,.... K) anchors, (r) = (dy = ||>; : i E S and dij <r} as the set of observed sensor-sensor distances and Dea (r) {d+R ||ar axll= : E S A and dus S r) as the set of observed sensor-anchor distances, both of which depend on some sensing radius
7

---

P25
p"
P2s
0.12
10a1
0.10
10a2
10a3
10-4
0.06
10-5
0.04
10-6
Pix
10-7
p Jc [o)]d((w)
0.00
14
12
Comparison between the objective value P2 from solving the degree-2s S-SOS SDP and the objective value p" resulting from the best-possible lower bound c"(w) for noise drawn as ~ Uniform(- 1,
Figure 1:
~ 0.1187 is plotted as the line in black and the p2. values are shown blue dots (left) with the gap between the values p" P2. (right).
1). 1) p² = (w)dv(()
Writing 4,, ag E [-1, as the unknown positions of the i-th sensor and the k-th anchor, we can write the potential function to be minimized polynomial:
fxx,w;
Σ -1 2,113 (:-:ddy()) (3)
Σ
(1|=; axl d=((w))3² (7)
dey de)
dsaED.a(r)
censor-sensor interactions
sensor-anchor interactions
The observed sensor-sensor and densor-anchor distances diy(w) (w), (u) can be perturbed arbitrarily, but in this paper we focus linear uniform noise, i.e. for a subset of observed distances have dij. = diy + ewn with Wk Uniform(-11 1), Other noise types may be explored, including those including outliers, which may be 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 global minimum of this function achieved at (XX, 04: A. r) 0 (where the distances have not been perturbed by any noise). In general for non-zero w (measuring distances under noise perturbations) we expect the function minimum to be > 0, as there may not exist 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 known anchors. This corresponds removing all sensor-anchor soft penalty terms from the function and instead selecting NH < N sensors at random to exactly fix known positions via equality constraints in the SDP. The SDP is still large but the effective number of variable sensors has been reduced to N' N NN-
A given SNL problem type is specified by spatial dimension N sensors, K anchors, a sensing radius (0, 2v7), noise type (linear), and anchor type (soft penalty hard equality). Once these are specified, we generate a random problem instance by sampling X ~ Uniform(-1,11) Uniform( 1, 1)d The potential f f(x, w) for a given instance formed (either with sensor-anchor terms not, with terms kept based some sensing radius r, and noise variables appropriately added).
~
The number of anchors is chosen to be few as possible so to still enable exact localization, anchors for a SNL problem in & spatial dimensions. The SDPs are formulated with
iee. K 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

---

3.2.2 Evaluation metrics
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 and the ground-truth sensor positions X, i.e. dist(x, X). Because the S-SOS hierarchy recovers estimates of the sensor positions E[æ.]] along with uncertainty estimates Var[z;], we would like to measure the distance between our ground-truth positions X to our estimated distribution p(x) N (E[x), Var[æ=]). The Mahalanobis distance das (Equation (8)) a modified distance metric that accounts for the uncertainty [38]. We use this as our primary metric for sensor recovery accuracy.
8 (X, 2)) = ((X- -X 14)
(8)
As our baseline method, for each problem instance we apply a basic Monte Carlo method detailed in Algorithm (aarthmmmmppp (Appendix A.7.3) where we sample w ~ vww), use a local optimization solver to find x ("(w) inf, far w), and use this estimate E. Ewww( [z], Var [a]. Note that though this non-SOS method achieves some estimate of the dual SDP objective f fxx w)du(ar w), it is not guaranteed to be a lower bound.
3.2.3 Results
Recovery accuracy. In Table 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 anchors K & 1, and select the sensing radius and the noise scale We then generate I, 20 random instances each problem type, corresponding to a random realization of the ground-truth sensor and anchor configurations X E [-1, xt A E []1,1] producing f(z, w)) that 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 (E[x Cov(x))) which we then compute 8as for (see Equation (8)), treating each dimension as observation of
number of
independent of each other (i.e. X as a flat vector). Each instance solve gives 5y or each method, and we report the median and the +1034% values over the L = 20 instances we generate.
4 Discussion
In this paper, we discuss the stochastic sum-of- (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 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, is relatively straightforward to extend this method to support other noise types (such Gaussian random variates without compact support, which we do in Appendix (.66.4) and support higher-dimensional noise with d >>>>> 1.
Scaling this method to larger problems n > is an open problem for all SOS-type methods. In this paper, we take the approach of sparsification, by making the cluster basis assumption build up block-sparse W. We anticipate that methods that leverage sparsity or other structure in f will be promising avenues of research, 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 fx, w) to exactly decompose into some lower- swwer-bounding e(w) and SOS fsos(a, Relaxing this exact- decomposition assumption and generalizing beyond polynomial f (xx w). c(w) may require novel approaches and would be an exciting area for future work.

---

Table I: Comparison of S-SOS and MCPO solution extraction accuracy. We present the Maha- lanobis distance dAs (Equation ( of the the true sensor positions X to the extracted distribution (E[2,) Var(z]) over solutions recovered from S-SOS for varying SNL problem types. is the spatial dimension, is the sensing radius used to cutoff terms in the potential [x, w). is the noise scale, Nu is the number of hard equality constraints used (sensors fixed known locations). No is the number of clusters used (see Appendix (A.7.4). and is the number of sensors used. Fach SNL problem instance has K L + anchors used in the potential (if Nø 0). The MCPO values are estimated with T = 50 Monte Carlo iterates. Each entry is f + ā where f is the median and robust standard- caniaio-deviation (034%) estimated over 20 runs of the same problem type with varying random initializations of the sensor positions. The entries with the lowest median Đas are bolded. We also compare the number of elements the full basis a,, the cluster basis ags and the reduction multiple when using the cluster basis a When passing to the cluster basis, a a is how much the semidefinite matrix shrinks by.
Parameters
Basis comparison
M-distance (δx)
r
e
NH
Nc
N
a
Gc
af/ag
S-SOS
MCPO
1
0.5
0.3
0
1
10
78
78
1x
0.94 + 0.22
2.61 + 3.86
1
1.0
0.3
0
1
10
78
78
1x
0.29 + 0.16
1.10 0.58
1.5
0.3
0
1
10
78
78
1x
0.11 + 0.11
0.86 - 0.52
1
1.5
0.3
2
1
10
78
78
1x
0.24 + 0.37
1.06 = 1.28
1.5
0.3
4
1
10
78
1x
0.10 0.03
0.61 + 0.41
1
1.5
0.3
6
1
10
78
78
1x
0.06 + 0.04 0.04
0.48 + 0.32
1.5
0.3
8
1
10
78
1x
0.02
0.31 = 0.17
2
1.5
0.1
0
9
9
406
163
2.5x
2.86 + 0.94
1562.39 + 596. 29
2
1.5
0.1
0
9
15
820
317
2.6x
3.25 + 1.19
1848.65 650.45
References
[1] Michael R. Garey and David S. Johnson. Computers and intractability: guide the theory of SP-completeness. A series of books in the mathematical sciences. Freeman, New York [u.a]. 27. print edition, 2009.
[2] Yurii Nesterov. Squared Functional Systems and Optimization Problems. Hans Frenk, Kees Roos, Tamás Terlaky, and Shuzhong Zhang, editors, High Performance Optimization, Applied Optimization, pages 405-440. Springer US, Boston, MA, 2000.
[3] Jean B. Lasserre. Global Optimization with Polynomials and the Problem of Moments. SIAM Journal Optimization, Applied Mathematics.
8an, 2001. Publisher: Society for Industrial and
[4] Monique Laurent. Sums of Squares, Moment Matrices and Optimization Over Polynomials. In Mihai Putinar and Seth Sullivant, editors, Emerging Applications of Algebraic Geometry, The IMA Volumes in Mathematics and its Applications, pages 157-270. Springer, New York, NY, 2009.
[5] Jiawang Nie. Optimality conditions and finite convergence Lasserre's hierarchy. Mathematical Programming, 146(1-2):97-121, August 2014.
[6] Francis Bach and Alessandro Rudi. Exponential convergence of sum-of-squares hierarchies for trigonometric polynomials, April 2023. 1.Xiv:2211.04882 [math].
[7] Mihai Putinar. Positive Polynomials on Compact Semi-algebraic Sets. Indiana University Math- ematics Journal, 42(3):969-999 984, 1993. Publisher: Indiana University Mathematies Department.
(8] Konrad Schmüdgen. The Moment Problem, volume 277 Graduate Texts in Mathematics. Springer International Publishing, Cham, 2017.
10

---

[9] Lucas Slot and Monique Laurent. Sum-of-squares hierarchies for binary polynomial opti- mization. Mathematical Programming, 197(2)6221666 January 2023. Publisher: Springer Verlag.
[10] Jean Lasserre. The Moment-SOS hierarchy, August 2018. arXiv:1808.0344 [math].
[11] Jean-Bernard Lasserre. The Moment-SOS hierarchy: Applications and related topics. To appear in Acta Numerica (2024), September 2023.
[12] Lorenzo Baldi and Bernard Mourrain. On the effective Putinar's Positivstellensatz and moment approximation. Mathematical Programming. 200(1):71-100) June 2023.
[13] Lorenzo Baldi and Lucas Slot. Degree Bounds for Putinar's Positivstellensatz on the Hypercube. SIAM Journal Applied Algebra and Geometry, 8(1):1-25, March 2024.
[14] P.A. Parrilo. Structured Semidefinite Programs and Semialgebraic Geometry Methods in Robustness and Optimization. PhD thesis, California Institute of Technology, Pasadena, CA, 2000.
[15] Jiawang Nie. Sum of squares method for sensor network localization. Computational Optimiza- tion and Applications, 43(2):151-77, June 2009.
[16] Etienne de Klerk. The complexity of optimizing over a simplex, hypercube or sphere: short survey. Central European Journal of Operations Research, 16(2):111-125, June 2008.
[17] Amir Ali Ahmadi and Anirudha Majumdar. DSOS and SDSOS Optimization: More Tractable Alternatives Sum Squares and Semidefinite Optimization. SIAM Journal on Applied Algebra and Geometry, 3(2):193-25) January 2019.
[18] Dávid Papp and Sercan Yildiz. Sum-of-Squares Optimization without Semidefinite Program- ming. SIAM Journal on Optimization, 29()):822-851, January 2019.
[19] Jean B. Lasserre. A "Joint+Marginal"" Approach to Parametric Polynomial Optimization. SIAM Journal on Optimization, 2094):1955-0022 January 2010. Publisher: Society for Industrial and Applied Mathematics.
[20] Bruno Sudret. Global sensitivity analysis using polynomial chaos expansion. Reliability Engineering & System Safery, 93:964 979, July 2008.
[21] Habib N. Najm. Uncertainty Quantification and Polynomial Chaos Techniques in Compu- tational Fluid Dynamics. Annual Review of Fluid Mechanics, 41(1):35-52, 2009. _eprint: https://doi.or//00. 1144///nnure..uui.00000000 165248.
[22] Kun Fang and Hamza Fawzi. The sum-of-squares hierarchy on the sphere, and applications in quantum information theory. Mathematical Programming, 190(1-2):331-366 November 2021. arXiv:1908.00 05155 [quant-ph].
[23] Lucas Slot. Sum-of-Square Hierarchies for Polynomial Optimization and the Christoffel- Darboux Kernel. SIAM Journal on Optimization, 32(4):2612-2633 December 2022.
[24] Ralph Tyrell Rockafellar. Convex Analysis. Princeton Landmarks Mathematics and Physics. Princeton University Press, Princeton, NJ, 2015.
[25] Stephen Boyd and Lieven Vandenberghe. Convex Optimization. Cambridge University Press, I edition, March 2004.
[26] Lieven Vandenberghe. Chordal Graphs and Sparse Semidefinite Optimization, 2017.
[27] Yian Chen, Yuchaw Khoo, and Lek- Lek-Heng Lim. Convex Relaxation for Fokker-Planck June 2023. urXiv:2306.03292 [cs, math].
[28] Francis Bach. Sums-of-squares for dummies: a view from the Fourier domain Machine Learning Research Blog, November 2022.

---

[29] Markus Schweighofer. On the complexity of Schmüdgen's Positivstellensatz. Journal of Complexity, (0(4)):29-544, August 2004.
[30] Frank H. Clarke. Generalized gradients and applications. Transactions of the American Mathematical Society, 205:247-262, 1975.
[31] Corbinian Schlosser, Matteo Tacchi, and Alexey Lazarev. Convergence rates for the moment- SoS hierarchy, February 2024. arXiv:2402.0044 00436 [math].
[32] Milan Korda, Victor Magron, and Rodolfo Rios-Zertuche. Convergence rates for sums-of-
squares hierarchies with correlative sparsity. Mathematical Programming, March 2024.
[33] Anthony Man-Cho So and Yinyu Ye. Theory of semidefinite programming for Sensor Network Localization. Mathematical Programming, 109(2-3):367388, January 2007.
[34] Aaron Meurer, Christopher P. Smith, Mateusz Paprocki, Ondřej Čertík, Sergey B. Kirpichev, Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K. Moore, Sartaj Singh, Thilina Rath- nayake, Sean Vig, Brian E. Granger, Richard P. Muller, Francesco Bonazzi, Harsh Gupta, Shivam Vats, Fredrik Johansson, Fabian Pedregosa, Matthew J. Curry. Andy R. Terrel, Stěpán Roučka, Ashutosh Saboo, Isuru Fernando, Sumith Kulal, Robert Cimrman, and Anthony Sco- patz. Sympy: symbolic computing in python. PeerJ Computer Science, 3:e103, January 2017.
[35] Steven Diamond and Stephen Boyd. CVXPY: Python-embedder modeling language for convex optimization. Journal of Machine Learning Research, 17(83):1-5, 2016.
[36] Akshay Agrawal, Robin Verschueren, Steven Diamond, and Stephen Boyd. A rewriting system for convex optimization problems. Journal of Control and Decision, 5(1):42-66) -0, 2018
[37] MOSEK ApS. The MOSEK optimization toolbox for Python manual. Version 10.0., 2023.
[38] PC Mahalanobis. On the generalized distance in statistics. In Proceedings National Institute of Science of India, volume 49, pages 234-256, 1936. Issue: 2.
[39] Ryan O' O'Donnell. SOS is not obviously automatizable, even approximately. R. 0, 2016.
[40] Matthew B. Hastings, Perturbation Theory and the Sum of Squares, June 2022. arXiv:2205. 12322 [cond-mat, physics:hep-th, physics:quant-phy
[41] M. B. Hastings. Field Theory and The Sum-of- Sum-of-SSuares for Quantum Systems, February 2023. arXiv:2302. 14006 [quant-ph].
[42] Samuel B. Hopkins. Statistical Inference and the Sum of Squares Method. PhD Thesis, Cornell University, August 2018.
[43] Samuel Hopkins and Jerry Li. Mixture models, robustness, and sum of squares proofs. In Proceedings of the 50th Anmual ACM SIGACT Symposium on Theory of Computing, pages 1021-1033 Los Angeles CA USA, June 2018. ACM.
[44] Saeid Sedighi, Kumar Vijay Mishra, M. R. Bhavani Shankar, and Bjorn Ottersten. Localization With Ont Passive Radars in Narrowband Internet-of-Thing Using Multivariate Polynomial Optimization. IEEE Transactions on Signal Processing, 69:2525-2555 2540, January 2021. ADS Bibcode:
[45] Dimitris Bertsimas, Dan Andrei lancu, and Pablo A. Parrilo. A Hierarchy of Near-Optimal Policies for Multistage Adaptive Optimization. IEEE Transactions on Automatic Control, 56 56(12):2809-2224 December 2011.
[46] Luigi Ambrosio, Nicola Gigli, and Savare. Gradient Flows in Metric Spaces and in the Space of Probubility Measures. Birkhäuser, second edition, 2005.
[47] Ivar Ekeland and Roger Temam. Convex analysis and variational problems, 1976.
[48] D. Jackson. The Theory of Approximation. Colloquium Publications. American Mathematical Society, 1930.
12

---

[49] Nathan Krislock and Henry Wolkowicz. Explicit Sensor Network Localization using Semidefi- nite Representations and Facial Reductions, January 2010. arXiv:1002.0013 [math].
[50] Brian D. O. Anderson, Peter N. Belhumeur, Tolga Eren, David K. Goldenberg, A. Stephen Morse, Walter Whiteley, and Y. Richard Yang. Graphical properties of easily localizable sensor networks. Wireless Networks, 15(2)117-111, February 2009.
[51] Brian D. O. Anderson, Iman Shames, Guoqiang Mao, and Bariş Fidan. Formal Theory of Noisy Sensor Network Localization. SIAM Journal Discrete Mathematics, 24(2):684-65) 06444669,444 January 2010.
13

---

A Appendix supplemental material
A.1 Notation
Let P(X) and P(tt) denote the spaces of polynomials on X S R" and Ω S R respectively, where and f2 are (not-necessarily compact) subsets of their respective ambient spaces I and R Specifically, all polynomials of the forms below belong to their respective spaces:
p(x)
Σ ca E P(X), aEZzo
p(w) Σ
Caw E P(22)
oEZ20
where x = (I1,... Fn). "ay....,,W,,, a is a multi-index for the respective spaces, and ca are the polynomial coefficients.
Let (S) for some (X,S2 denote the subspace of P(S) consisting of polynomials of degree < d,, i.. polynomials where the multi-indices of the monomial terms satisfy ||a|| < d. Psos (X x 52) refers to the space polynomials on x Q that can be expressible as a uum-of-squares in x and
jointly, Additionally, W > 0 for matrix W denotes that W is symmetric positive semidefinite (PSD). Finally, (P)) denotes the set of Lebesgue probability measures on f.
A.2 Related work
A.2.1
Sum-of-squares theory and practice
The theoretical justification underlying the SDP relaxations in global optimization here derive from the Positivstellensätzz (positivity certificate) of [7), a representation theorem guaranteeing that strictly positive polynomials certain sets admit sum-of-squares representations, Following this, 10 developed the Moment-SOS hierarchy, describing a hierarchy of primal-dual SDPs (each having fixed degree) of increasing size that provides monotonic non-decreasing sequence of lower bounds.
There is rich theory underlying the SOS hierarchy combining disparate results from algebraic geometry [14] [0] [T], semidefinite programming [15] [8), and complexity theory [I6] 39]. The hierarchy exhibits finite convergence in particular cases where convexity and strict local minimum are guaranteed [5), otherwise converging asymptotically [6]. In practice, the hierarchy often does
better than these guarantees, converging exactly at c for some small
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 high-dimensional sta- tistical 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 PSD matrix W E x=(a.dd)) with ann, d, s) scaling as (*+d+**) for 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 Work attempting to improve the scalability of the core ideas underlying the SOS hierarchy and the SDP method include [7. 18).
A.2.2 Stochastic sum-of-squares and parametric polynomial optimization
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 optimal solutions converges to moment vector of probability measure encoding all information about the globally-optimal solutions (=w = argmin, aegmin_[(xx40 w) and provides proof that the dual problem (our primal) obtains polynomial approximation to the optimal value function that converges almost lifost-unifommly to c" (w).
14

---

A.2.3 Uncertainty quantification and polynomial chaos
Once physical system 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 w) over x where w parameterizes the function and is drawn from some noise distribution v(w).
If only function evaluations (ar, w) are allowed and no other information is known, Monte Carlo is often applied, where one draws w/ ~ v(w) and solves many realizations inf, fx(x) = f(x, we) approximately solve the following stochastic program:
inf Eu (f(r
Standard Monte Carlo methods are ill-suited for integrating high-dimensional functions, so this method is computationally challenging its own right. In addition, 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 (x,
Our approach to quantifying the uncertainty in optimal function value resulting from uncertainty in parameters is to find deterministic lower-bounding c* c"(w) which guarantees (x) 2 (w) no matter the realization of noise. This is reminiscent of the polynomial chaos expansion literature, wherein system of some stochastic variables is expanded into a deterministic function of those stochastic variables, usually in some orthogonal polynomial basis [20 21].
A.3 An example
Example A.1. Let f(x,w) be some polynomial of degree I 2s written in the standard monomial basis, L.e.
[xx
Σ
Jax
||all1520
= Σ
II x
auti
w
||all152e
Let (2, w) E (R((,,d,)) be the basis vector representing the full set of monomials in æ of degree S with a(n, s) = (x+d+*).
For all znd with ||a/|1 I 2s and ag 0 for all le (I. (i.e monomial terms containing only wi,..., ......,,, wa) we must have:
/ Jxxa du((,, - (ddu(w)==
Explicitly, for μ to be valid probability distribution must have:
Sane du(r, w) - =Mopo - 1 = V(1,0....) - =
Suppose ! = [-1, 1], ~ Uniform(-1, so that d = 1, v(w) = 1/2. We require:
1
0
{ Jxx d'dp(,,) [
(ndu(w)
1 a =
3
1-111
0
3
15

---

A.4 Strong duality
To guarantee strong duality theoretically, we need a strictly feasible point in the interior (Slater's condition). For us. this is consequence of Putinar's Positivstellensatz, if [(x w) admits decom- position as [(r, w) = c(w) + g(a, w) where g w) > (i.e. strictly positive), we have strong duality, ie. p- d" and P2. d2. 3 S 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, p2s d2 in all cases encounter where SDP solver returns feasible point.
A.5 Proofs
A.5.1
Primal-dual relationship of S-SOS
Regular SOS Global polynomial optimization can be framed as the following lower-bound maxi- mization problem where need to check global non-negativity:
sup
(9)
cER
s.t. f(z) c 0 Vr
When take the SOS relaxation of the non-negativity constraint in the primal, arrive the SOS primal problem, where we require f(x) c to be SOS which guarantees non-negativity but is a stronger condition than necessary:
sup
(10)
oeR
f(x) - PPoss
The dual to Equation
9 is the following moment- commenttmiiimizatioon problem:
/ ((x)du(a))
inf mEP(X)
(11)
with fants)= =1
Taking some spanning basis m,(x)) R Re(n,s) of monomials up degree 8, have the moment matrix M (a(n)))xa().));
Mr,
m. (1) = Yo
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 moment-miniminn minimization problem, which is now in solvable numerical form:
Elaya inf
(12)
y
with M(y) >
M(y)11 =1
where we write M (y) as the matrix formed by placing the moments from into their appropriate places and the first element of (2) be 1, hence M1,1 dd(()) is simply the normalization constraint. For further reading, see 15 3).
Now let us lift this problem into the stochastic setting with parameters w sampled from given distribution L.e. replacing (x, w). We need to make some choice for the objective. The expectation of the lower bound under v(w) is a reasonable choice, i.e.
Stochastic SOS
/ c(w)du(w)
16

---

but we could also make other choices, such as ones that encourage more robust lower bounds. In this paper however, formulate the primal S-SOS below (same as Equation (T):
/ ((w)ddl())
sup cEL)(
(13)
s.t. [(x,w) c(w) 2 0
Note that if the ansatz space for the function c(w) is general enough, the maximization the curve e is equivalent to pointwise maximization, i.e. we recover the best approximation for almost all w.
Then the dual problem has a very similar form the soo-stochastic case.
Theorem A.2. The dual Equation 13 is the following moment minimization where pla w) is probability measure on x Qz
inf MEP(XxS) xX)) with
fxx w'dp(o,w)
I ddd((w) for all N
Remark A.3. Notice, that the condition Sxxa dulll, w) = Sa duw() implies that the first marginal of μ is the noise distribution и. Let (A. denote the disintegration of with respect to v, [46]. Then the moment matching condition is equivalent to M (X) 1 for almost all w and j being Young measure w.r.t. v. The idea is that (x) a minimizing density for every single configuration of w.
Proof. We P2o(X Ω) denote the space of non-negative polynomials on X Ω. Given measure v on Ω and polynomial function p X x + R consider
/ ?(w)du(w)
sup TELL (0.v) gEP2a(XXX))
12
s.t (e(w)->((c)) = g(x,w)
This is equivalent to
inf (R.)) (Xx2)
f(r,g) + 9(7,4)
with
((..g = 1/ 7(w)du(w)
and
f
g(7,9) -X(---0=-0)0 -00 else
ie.. g is the characteristic function enforcing non-negativity. Denote by h' the Legendre dual, i.e.
h (y) sup(r, u) - 1(x)).
x
Then by Rockafellar duality, 47 24], and noting that signed Borel measures are the dual to continuous functions, the dual problem reads
(
sup TELM(R,))
and we would have
sup TEL EB
f(%1) + g(r.)). Ter
inf
qEP
17

---

The Legendre duals of f and g can be explicitly calculated as
(T) = [0ifr= = -1 and S 0
(o0 else
and
9 (T) 1 ((x,w)))(w,z) f P20(X ②) and [(w) (X)
else
since
((W))(W))a) + SOEX ((a, ()d,p( - s(7,2)
(T( = sup
(Tw) + ()dr(w)+ + faxx ((a,w)d
sup
7.0
nn
0 ifT -1 and 1 00
-00 else
and
g (T, sup / (w))(() + I glo,w)ddp(e, + X(S-1-9=0) 7.
sup / [(W)[(w)d nxx
(f(x,w) - (wW))d((wxx)
- =E P20(X x 22)
else
[[))([(w) (x)) + / (f(a,))da(,xz)
sup n
E x 52)
nxx
else
(d((ww) if f P2o(X x 22) and I'()) (X)
5xx
else
Altogether, we get
[ f(x,w)du(wxx)) (x) =
2
else.
A.5.2
Convergence of S-SOS hierarchy
(w) at best Lipschitz continuous By Proposition 2.1 we argue that c² (w) inf, (x, w) Lipschitz continuous. One cannot expect much more the tollowing example shows:
Example A.4. Consider g R R R defined by
g(a,p,4) = (u² + px + 2)
Then we have for every (p. q) IR? that
0
2
inf
g p)
*ER
[ 9) else.
Therefore, (p.q) + infeeR g(z, (. is once differentiable but not twice.
18

---

Lemma approximating polynomials
Lemma A.S. Let Q be compact subset of Rn and Ω - IR be Lipschitz continuous. Then there
exists a trigonometric polynomial 9s of degree s and a constant C >> depending only on S and such that
g 2 9s
and
+ In(s)
llg -0 -
Cl
Proof. By Jackson's inequality for Lipschitz functions polynomial g of degree s with
48
we have the existence of trigonometric
C¹
lgg - gf ('nn S
| m*(2))
as well as
lgg 9 (1)000) S T.g In(s)
Then we define ys = 11g - g g'llz= and hence 2 gs- Furthermore,
llg - gallea(a) S (C" + +
Lg
Writing C(Q) ma | have the desired form where ||2] the volume of Z.
Convergence at s/s/s rate using integral operator methodology
TR. be trigonometric polynomial of degree 2r, c² (w) = inf fxx w) the optimal lower bound as a function of w. and any probability measure on compact C IR.. Let s = (82, Sc). 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(w).
Theorem A.1 (Asymptotic convergence of S-SOS). Let [0, x S
Let oz (w) be the lower bounding function obtained from the primal S-SOS SDP with my(ar, w) spanning basis of trigonometric monomials with degree S 8x terms and of degree 8w in w terms:
/ c(w)dv(w)
P2s
sup (1),
f(z,w) c(w) " Wm,(zw
Then there is a constant C > O depending only on S2 d, and n such that for all ows 8x 2 max(3r 38.} the following holds:
/
(c((w) cta dulw) & s)
12
2r
dd
6r =
((.8) S|1 - 5||
1
22
s?
6r
-d
+ cc In((2s))
+ 1|ee" -c'lle
1
28c
where f denotes the average value of the function f over [0, , t.e. f = S0 f(x) and [)(())lr=E> f(() denotes the norm of the Fourier coefficients.
e(f,s) bounds the expected error, giving us asymptotic convergence
min(s Swi sc) OO. Note the first two terms give O ( convergence rate. However, the overall error will be dominated by the degree of clw) (from the third term) hence convergence rate is O mm.)
19

---

Proof of Theorem A.] Let ? C R be compact and f R" x 1 + R be I 1-periodic trigonometric polynomial (t.p.) of degree Σ 2r. We then make isomorphic [0, 1 and hereafter consider
Q 0, 1 and f [0, 1 x [0, 1]d R. Let and 5. Let the best lower bound be
c¹ (w) = inf (x.
Proof outline. We split the error into two parts. First, we use the fact that there a lower-bounding Lp. ca of degree such that
||e - call S in
Se
and
c
This will provide us with a degree-8, t.p. approximation to the lower bounding function, which in general is only known to Lipschitz continuous.
Next, we show, that for any b > 0 there is degree-2s SOS tp. fsos(I, w) such that
fsos f (c -).
We write 8 = (82, Sw) where 8z, Sw denotes the respective max degrees in the variables w. Once we have constructed this, we can compute fsos = cg - and since we know that Ssos 0 everywhere and eg is some degree-s, tp. we have found a degree-s. lower-bounding t.p. The construction of this SOS t.p. adds another error term. If we can drive - 0 as 5 min (Sps Swi se) oo then are done.
Proof continued. To that end, let c Ω - ER be the best degree-s, trigonometric approximation of e" with respect to I such that
By 30 we know that c* is locally Lipschitz continuous with Lipschitz constant Le and hence, by get that there is C(1)) > 0 such that
Lemma A.S
1 In Sc
|c cal S C(22)
Le
Sc
Next we introduce cz (w) which some degree-2s t.p. After an application of the triangle inequality and Cauchy- -echwa- the integrated error term Sn |e c2 ldw we have
FEX La inf f(x,w) 6 (el S / (c](w) c2. (w)) [allic - callz>
S. inf (x, - czolw)) - #EX
[ (c((w)- (d)
1 In Se 8c
cz
C(S2)
Le
gap between some SDP solution ei, (w) and tp. (w)) approx. of L-contin.
Now we want to show that for any & > 0 we can construct a degree-2s SOS trigonometric polynomial fsos (x. w) such that
fsos = f - c² b.
with 6 = 2/2 and = (sn, sw) > r. We can then set f - fsos = c b c₁ as the degree-2s lower-bounding function. If we can drive b e/2 + 0 as 8, lc
are done, as by construction
c5. =
Observe that by assumption f c b a t.p. in (x,w) where f is degree-2r and c is degree & 2r. Denote by (f - f + b)w its coefficients w.r.t the w basis. Note that the coefficients are functions in z. Following the integral operator proof methodology in [6), define the integral operator
to be
Th(x, = f
Jxxn (w (1 2) 2h( @ dzddo
where y is a trigonometric polynomial in of degree S and ga is trigonometric polynomial in x of degree S 8z- The intuition is that this integral operator explicitly builds SOS function of degrees (8z, sw) out of any non-negative function h by hitting it against the kernels qu, qw-
20

---

We want to find a positive function h X x Ω - IR such that
Th = f - c2 + b.
In frequency space, the Fourier transform turns a convolution into pointwise multiplication so we have:
Th () = سw * qw (w)) 4- * (z) h(a, w)). In the Fourier domain it is easy to write down the coefficients of h: if || لل > (2r, 28c)
h(t,w) { ((, (0)) b1==0la=0 i (w0) Ou * (3)
otherwise.
Computing Th - h gives:
f(2, w) (w(w) + b- h(z,w)
Σ (2)
(1- Aw (1- (00) - 2u
exp exp (2imm" w) exp( (im x)
(x)
Σ (bl==0 -(a) (1-
exp(2imm w)
fm * 9 (3))
and thus after requiring &. * (0) = d * ân = we have:
max (xx ca( (w)) + b h(2,,w) w)l
IIlf fllF max max 370 870
Fw qw (w)
Ox * (i)
+ max |,o² - E11 2x4
i qw (w)
As reminder, because cº 2 c everywhere we have f - 2 forrf - b > 0, since
1 = ==/2: > O. Since Th = f - cg + b> and is SOS, we need to guarantee 0.
maxx, ((x - (w) + Mxx w & then
max Th h| b. 2,4
Since Th 2 and b > 0 have
= Th + h - Th 2 Th || - Th 2b-b20 0
and hence h > 0 if ensure max
[Th h < b.
Now let us show that
max f(z,w) c% (w) + b - h(x 5
can be ensured s (s., sw) is large enough.
Using the same kernel and bounds as in [6]. we choose for E (x, wh the triangular kernel such that
(z(3) = (1 -1--
6r
d
d
|2
2
$2,0
Note that (x)+ = max (2r 0). Then we have
max [f(I, (w) + b h(x,
1
1
<IIʃ fllF max
max
1
+ ||0 - Eall max
ii (w) qx (+)
qu
q (w))
(1
2r -d
2
+ ||c; - -|1r - (1
62
<|f flr 1
1
$2
22
2
21

---

Therefore, by choosing and 8x large enough such that
Il - fll 1- -
6r
-
1
6r
+ l
-Elll
1
62 -
S b
8 2
2
s
we have
120 0
and thus Th is SOS. By design we have
ca - 62 <b
and thus
J lca Ggldd dww S
Recalling
La inf [xx - ( 5 =ex
1 (c((w) -
1 In sc 3c
C(S2)
Lee
gap between some SDP solution cz (w) and (w(w
approx. ervor of L-contio.
we can additionally choose & large enough to guarantee
In 8c
C(S2)
8c
2
and then we are done.
Setting 8a, 8w, 8C and sending 8 - oo we have asymptotic behavior final error expression:
J inf zex [(x,w) (2) dw S si 82
+ +2- + C3 Ins
In
with the constants C1, C22, C33 depending on r, n, |f flle. llco E. Ω and Le
Convergence at 1/s rate using picsewise-constant approximation to c"(w) Prior work [6] achieves 1/s² convergence for the regular SOS hierarchy without further assumptions. our work, we could only achieve sns/s for S-SOS due to the need to first approximate the tightest lower-bounding function c"(u (w) with a polynomial approximation, which converges slower rate.
To accelerate the convergence rate, we want to control the regularity of c* (w). We can achieve 1/s by approximating the c² (w) pointwise instead of using smooth parameterized polynomial. By constructing domain decomposition of Ω and finding a SOS approximation & for each domain, we stitch these together to build piecewise-constant approximation to the lower-bounding function c" This lets us leverage the s/8² rate for the standard SOS hierarchy, however note that we lose lower bound everywhere. Another
the guarantee that the lower-bounding function eg (w) is disadvantage of this approach is that the computational complexity scales exponentially in as we required to solve for the SOS lower bound with degree pointwise at (ss)) points.
For Ω [0, 1]2t Idd we achieve the following:
Proposition A.6. Let Ω [0, 1jd and f [0, 1] - R be trigonometric polynomial of degree 2r. Let {w,) be grid points where 8p is the number of equally-spaced points along each dimension for total of (s,)) points in Q. Denote by c% (wi) the best SOS approximation of degree of x + f (x, wi)
each grid point we and define
of
(w(w Σ (w(l))leBBBll
where Boolwir e) denotes the Loe ball of radius centered at wi with (2s,) so that the balls are non-overlapping. Then we have for some constant C' depending only on maxu; [[f(wis f(wi, ))], n, spi
/ (c((w) (w)) & max (If(w.) F(wi, Alle 1 2
2r
-n
Cvd
< C
wi
$2
Sp
22

---

Proof of Proposition A.6 Let (w) 10, + R be the best piecewise-constanl approximation of
the true lower-bounding function c* (w) inf, f(x,w) on equidistant grid-points, i.c. Sp evenly-
spaced points in each dimension for total of o points.
Now also define a piecewise-constant approximation using the best pointwise SOS lower bounds:
c (w)
C (w
where (w,, is the best lower bound (resulting from regular SOS) of degree of x: f w.) and 3o (Wi, e) denotes the Lee ball of radius e centered at Wi with e = (2sp) -1 so that the balls are non-overlapping.
By we have that (w:) - cg (wss can be bounded by
max (1xx))) f(we, 1 11-(1-- 6r )
2
Then have:
J [c'(w) (w]]] &
Σ (w)(w)) c, + +Ic" - -c
Using the same bound we get for the first term from the proof of Theorem A.I term a 22) dependence and approximation to 1-periodic trigonometric polynomials from [48] for the second:
we can reduce the first the theorem the L convergence of piecewise-constant
S [lo(w)ss - - (ww)) d()]]d & max [lf(;, - f(,ww)lle 1 (1-
6r
Cvd
|
+
wi
s²
Sp
straightforward to extend this to arbitrary compact Ω C Rd, as one can construct a smooth homeomorphism between Q and [0, 1] and therefore can achieve the same convergence rates general compact Ω as we can for [0, 1]d.
A.6 S-SOS for a simple quadratic potential
We provide simple application of S-SOS to simple quadratic potential that admits closed-form solution so as to demonstrate its usage and limitations.
A.6.1 Analytic solution for the lower bounding function c* (w) with (iifform-1,11
Let E IR and w Uniform(11,1 Suppose that we have
((x,w) (x w) (wx)²
In this case we may explicitly evaluate the exact minimum function c"(w) = inf, f(x; Note that f(r:w) x² 2wx + w² + 2²
Explicilly evaluating the zeros of the first derivative have
Oz (xf(x;w)
2x 2w + 2w? 0
x² (1 + w² =
x² + wz
and, thus,
c² ('(w) = inf f(; = 1+
Note that despite f(x w) being simple d SOS polynomial, the tightest lower-bound (w) = inf, [(x, w) is explicitly not polynomial. However, it is algebraic, as it is defined implicitly as the root of the polynomial equation
c (w) w²) wt =0
23

---

(س(WW)
c (w) - -c,(w)
0.5
(w)
0.12
4
2s
0.4
0.10
0.3
0.08
0.2
0.06
0.04
0.02
0.00
-1.00
-0.75
-0.50
-0.25
0.00
0.25
0.50
0.75
1.00
-1.00
-0.75
-0.50
-0.25
0.00
0.25
0.50
0.75
1.00
Figure 2: Lower bound functions for basis function degree d = 2, (left) and the optimality gap to the true lower bound c" c"(w) c2. (w) (right)
A.6.2 Degree-2s S-SOS to find polynomial lower-bounding function (w)
Observe that the tightest lower-bounding function ("(w) is not polynomial even in this simple setting. However, we can relax the problem to trying to find C2. E (2)(()) to obtain a weaker bound with inf, f(x, w) c* (w) 2 (w))
now proceed with formulating and solving the degree-2s primal S-SOS SDP (Equation (2). We assume that 62
(2s(w) (w) is parameterized by polynomial of degree < 28 in w. Observe that this class of functions is not large enough to contain the true function (w).
We choose E {2. 4} and use the standard monomial basis in æ, x,w, we have the feature maps m2(I, w) : 2 + g and m4(x,w) w) : _ + 1s since there are (m+ts unique monomials of up to degree-s variables. These assumptions together enable us to explicitly write SOS SDP in terms of coefficient matching. Note that we must assume some noise distribution vww). For this section, we present results assuming w Legendre quadrature with k zeroes on [-1, to evaluate the objective f c(w)) In fact, sample points suffice to exactly integrate polynomials of degree 2k 1.
(iiffrmm-1,111 We solve the resulting SDP in CVXPY using
We solve the SDP for two different levels of the hierarchy, = 2 and 8 = 4 (producing lower-bound polynomials of degree and 8 respectively). and plot the lower bound functions 02. (w) vs the true lower bound c* (w) = w /(1 + w²) as well as the optimality gap to the true lower bound in Fig
A.6.3 Convergence of lower bound as degree & increases
To solve the S-SOS SDP in practice, must choose maximum degree 2s for the SOS function mz(x, w)' WW)T (z, w) and the lower-bounding function c(w), which are both restricted to be poly- nomials. Indeed, larger not only increases the dimension of our basis function m, (x, w) but also the complexity of the resulting SDP. We would expect that džs + d" as 8 + co, i.e. the optimal value of the degree-2s S-SOS SDP (Equation (4h) converges to that of the "minimizing distribution" optimization problem (Equation (3)
In particular, note that in the standard SOS hicrarchy typically find finite convergence (exact agreement at some degree 2s* < oo). However, in S-SOS, we thus far have only a guarantee of asymptotic convergence, as each finite-degree S-SOS SDP solves for polynomial approximation to the optimal lower bound c* (w) = infzex f(x, w). In Figure [ objective values
we illustrate the primal S-SOS SDP
/ c(w)dr(w) with f(x,w) c(w) E Pios (X x Q!
P/ = sup
cEP
24

---

8
clw)). 0= 0.001
c((w), 0.01
c4(w), 0= 0.1
6
c((w), o= 1.0
cl(w), 0= 10.0
(w) 1+
4
2
0
-3
2
Figure 3: Different lower-bounding functions for degree-4 S-SOS done on the simple quadratic
potential [(x, w) (x - w)² + (wx)²² The true lower-bounding function c* (w) is plotted in black.
for given level of the hierarchy (a chosen degree for the basis m, (z, w)) and their convergence towards the optimal objective value
f (w()du(w))
2 0.1187
3
for the simple quadratic potential, assuming v(w) with w Uniform(-1, 1). We that in the log- linelinear plot (right) we have a "hinge" -type curve, with linear decay (in logspace) and then flattening completely. This suggests perhaps that realistic scenarios the degree needed achieve close approximation is very low, lower than suggested by our bounds. The flattening that occurs here is likely due the numerical tolerance used in our solver (CVXPY/MOSEK), as increasing the tolerance also increases the asymptotic gap and decreases the degree which the gap flattens out.
A.6.4 Effect of different noise distributions
In the previous sections. we assumed that w Uniform(-1, This enabled solve the primal exactly using Legendre quadrature of polynomials. Note that in Figure I we see that the lower-bounding ez (w), eq (w) for Uniform(-1. 1) is 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 g present the lower-buuding functions (w) achieved by degree-4 S-SOS by solving the dual for w Normal(O, o²) for varying widths a. We can see that for small a < 1, the primal solution only cares about the lower-bound accuracy within small region of w 0, and the lower-bounding curve fails to "generalize" effectively outside the region of consideration.
A.7 S-SOS for sensor network localization
A.7.1 SDP formulation
Recall the form of ((x, w):
fx, X, A, =
Σ
(1|=; - 1,113 22
(112 2x (w) du (r.())
Σ
(ryED.r(r
sensor-sensor interactions
sensor-anchor interactions
25

---

Note that the function f(z,w) is exactly degree-4 SOS polynomial, so it suffices to choose the degree-2 monomial basis containing
( elements as m2, w) RNeed - IR" That is, we have N sensor positions in spatial dimensions and d parameters for total of N + d variables. Let the moment matrix be M E Raxe with elements defined as
Mij = [ mi (I,, w))m 2 (x, w)dp(rww)
for J E { a). which fully specifies the minimizing distribution p(z, w) in Equation Our SDP is then the form
di = inf Σ faya V
s.t. M(y) x 0 Yo ma A (a, (a, M2
yo A (a, E 7
where Yo Ma corresponds to the moment- noment-matching constraints of Equation (4) and yo y correspond to any possible hard equality constraints required to set the exact position (and uncertainty) of sensor [[xx] x E[= [[a;]² O for all w. MM represents the (d+2*) moment-matching constraints necessary for all moments w.r.t. w and H represents the 28n constraints needed to set exact positions of known sensor positions in R (i.e. constraint per sensor and dimension, 2 each
for mean and variance).
A.7.2 Noise types
In this paper focus on the linear uniform noise case, as it is more accurate reflection of measurement noise true SNL problems. Special robust estimation approaches may needed to properly handle the outlier noise case.
Linear uniform noise: for subset of edges we write dij.a(w) = -d + owk, We ~ Uniform(-1, 1), and c 2 some noise scale we set. The same random variate we may perturb any number of edges. Otherwise the observed distances are the true distances.
Outlier uniform noise: for subset of edges we ignore any information in the actual measurement dij.k Wk, wk Uniform(D. 2v() where & is the physical dimension of the problem, i.e. 24 E R
A.7.3 Algorithms: S-SOS and MCPO
Here we explicitly formulate MCPO and S-SOS as algorithms. Let X = R Rd and use the
standard monomial basis. We write [x1,..... ., Fn W1,..., w.l] Our objective is to approximate (w) = inf, [(x w) for all w, with view towards maximizing ( (w)dv(w) for w sampled from some probability density v(w).
MCPO (Algorithm) simply samples w; and finds a set tuples (-² (wy), wi)) where the optimal minimizer (x" (wp, W is computed using a local optimization scheme (we use BFGS).
S-SOS (Algorithm via solving the dual (Equation (4) is also detailed below.
A.7.4 Cluster basis hierarchy
Recall from Section 2.2.2 that we defined the cluster basis hierarchy using body order b and maximum degree per variable In this section, review the additional modifications needed scale S-SOS for SNL.
In SNL, far w) is by design degree polynomial in = [z, w). with interactions of body order b = 2 (due to the (x;,2;) interactions) and maximum individual variable degree t = 4. Written this way, we want to only consider monomial terms [zr w/ with |||]]] s. ||a||oo and ||a|]o 2. To sparsify our problem, we start with some k-clustering (k clusters, mutually-exclusive) of the sensor set C (CI...., C.). This clustering can be considered leveraging some kind of "coarse"
26

---

Algorithm I Monte Carlo Point Optimization (MCPO)
1: Input: Function (); sampler for distribution v(w), number of samples T
2: Output: Approximate integral I, empirical distribution pr(œ), empirical mean empirical covariance Σ
3: for to 7 do
4: Sample w, v(w)
5: Find minimizer T, min, f(x;ws) using BFGS
6: end for
7: Estimate integral R / (Tzww)
8: Construct empirical distribution
T pD( = T 1=1
Σ 8(x xz)
9: Calculate empirical mean j = + Σ 3y and covariance &
T-I Στ
(xz - [)) - -2)²
Algorithm 2 Stochastic Sum-of-squares (S-SOS), Dual formulation
1: Input: Maximum basis function degree s E Z>o, complete basis function m, (,(x,w) Rn+d R( ( function (x;w) 12+ R represented as a dictionary mapping multi-index a E Zafd + coefficient fas probability density function for v(w) with known moments I du())( < oo v ||a/11 & 2s, any hard equality constraints where we want to set 2% = zx for some K.
2: Let i1, 12, i4 be the lexicographically-orreee arrays
2 ((++d+))((xtd)
2//0 ("+4+*))x ((++d) (
(n+d+2)) (n++d) -20
which correspond to the arrays of multi-indices for all degree- 1, degree-s, and degree-2s mono- mials in the variables 2,
("*d"")) ("+d+*) as a matrix of variables to be estimated.
3: Create M IR
(++++))) 20
4: Create y E R
as vector of variables to be estimated, corresponding to the vector of independent moments that fully specifies M.
5: Add M > 0 constraint.
6: for in length(i2) do
Require M to be formed from the elements of 8: Compute 0ij = :2[i] + (2]]] as the multi-index corresponding to the sum of the multi-
7: for j in length(i2) do
indices i2 i2]...
9: Add constraint M. = ya.y
10: end for
11: end for
12: for each row in 14 do
Require yo moments to equal the known moments of w
13: if EP= O; = 0 then
14: Add constraint ya = f 2 dv(w) (all)d)du(w).
end
16: end for
for do
b Handle any hard equality constraints in our variables Form multi-index an Zo+d where the entry for æy is set to I and everything else is zero.
19: Form multi-index 022 E Z where the entry for x2 is set to and everything else is zero.
20: Add constraint Ya, xT.
E[xx]] 2/-
21: Add constraint You (xt)².
Var[ur]] = B[r
E B(ra]* P 0.
22: end for
23: Form the objective to be minimized: F F=
Eaest Sall..
24: Solve SDP where we compute inf F subject to above constraints.
25: Output: If the problem is feasible (i.c. there exists degree 2s decomposition of f into fsos and
(w)), return moment matrix M E IR( (""**=")xx****= ("+d+*)
dual objective value dz Otherwise.
terminate and return failed/infeasible SDP solve.
27

---

information about which sensors are close to each other. For example, just looking at the polynomial (z,w) enables us see which sensors (i. j) must be interacting.
Assume that there is some a priori clustering given us. We denote 2:(1) as the subset of the variables restricted to the cluster C₁, i.e. _²() {x, j E Ci). Moreover, let G = (V, E) be graph where the vertices V = {1....,k] k} correspond to the k clusters and the edges correspond to known cluster-cluster interactions.
{(i.j) : i,j E V}
The SOS part of the function f() 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²(x)" m,(x)"Wmo(an (m)(a))
Σ (m (ellitTTooom w + VEV
Σ (KJ)E
(ol0)/TWoea (
where (kk) are symmetric PSD matrices and ((.) are rectangular matrices where we require W(()) = (WG.+)TT mx(x) for x E R" here behaves as before and denotes the basis function generated by all ("++) combinations of monomials with degree I 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, (z) on the full variable set will contain terms that mix variables from two different clusters that may not have edge connecting them.
Efficiency gains in the SDP solve occur when we constrain certain of the off-diagonal W(+)) blocks be zero. i.c. 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 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 [xr, w) that correspond to cross-terms that are "ignored" by the particular cluster basis we chose. The resulting polynomial f(x, w) 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 that have Nc ground-truth sensor positions and use those sensor labels to create our partitioning of the sensors. We connect every cluster using plus-one c; C++1 (including the wrap-around one) connections, so that the cluster-cluster connectivity graph has No: edges. We then use this information to throw out observed distances from the set D.s and from the full basis function m2( w).. See our code for complete details.
we do Nc-means clustering of the
A.7.5 Hard equality constraints
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 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 anchor if it is within the sensing radius).
To deal with the equality constraints where we set the exact position of sensor x; = xT, we solve the dual Equation (4) and implement them as equality constraints on the moment matrix, i.e. for the basis element m2 [z., w); = Ii we may set E[r;]] x. = Mo.; x; = 0. Note that we also need to set Var( = 0 so for m2(x, w)) of we add the equality constraint Var(z;) [[= E[r,
Mo.j MG. O.
A.7.6 Solution extraction
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 (ar by inspecting the appropriate entries Mo, corresponding to m2(2, w)); x; and Mo.; corresponding to m²(x,w), x².
28

---

M-distance
MCPO 5-50S
100
10-1
10
50
100 150
200
250
300
Figure 4: Comparison of the performance of MCPO and S-SOS (degree-4) for sensor recovery accuracy in 1D SNL with varying number of samples T used the estimate of empirical f. L. M-distance is dM. our metric for sensor recovery accuracy per Equation (8). The problem type here is NN sensor,
spatial dimension, |S2] 2 noise variables, 0.1 noise scale, sensing radius problem. The full basis is used here for the S-SOS SDP.
A.7.7 Impact of using MCPO with varying numbers of samples T'
In Figure 4 we can see how das 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 30 seconds for the T 300 MCPO point. Despite taking 100x longer, the MCPO solution recovery still dramatically underperforms S-SOS in 6 This reflects the poor performance of local optimization methods vs. global optimization method (when it is available).
A.7.8 Scalability
The largest 2D SNL experiment we could run had N 15 sensors, No clusters, and noise parameters. This generated Nł d = 39 variables and 820 basis elements in the naive m2(#, w) construction, which was reduced to 317 after our application of the cluster basis, giving us W, M E R317x317 A single solve in CVXPY (MOSEK) took 30 minutes on our workstation (2x Intel Xeon 6130 Gold and 256GB RAM). We attempted run with = 20 sensors and Ne = clusters and d = 9 noise parameters, but the process failed due to OOM constraints. Thus, report the largest experiment that succeeded.
A.8 Limitations
As discussed in the main text, our work provides both a theoretical convergence rate well as numerical understanding of the stochastic sum-of-squares (S-SOS). Notably, our convergence rate bottlenecked on the approximation quality of ct (w), only achieving In s/s where regular SOS results achieve 1/s2. In practice, we see much faster convergence, which may be a clue that other more restrictive assumptions more powerful methods may enable one to find better convergence results.
In addition, our numerical work focuses 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 N 15. This seems much smaller than the state-of-the-arr i.e. 44 where the authors propose an algebraic reduction of the noiseless SNL problem so that they can analytically obtain the range of PSD matrix. This simplifies the SDP dramatically, however sensor-sensor distances are perturbed with observation noise) this approach is unusable, requiring significant modification. his reduction used to solve SNL problems of 10k-100k sensors, but they develop a highly specialized algorithm that does not any SDP solvers.
noisy setting (where observed
29

---

Note that noiseless SNL simplifies the problem. In the noiseless setting, localizing even small number sensors near an anchor will propagate the correctly localized positions 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] [1].
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 solution for the sensor positions for every possible configuration of noise via the probability distribution μuz,, w). This additional complexity places further limits on the size of the problems that can be solved with S-SOS.
30

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the
paper's contributions and scope?
Answer: [Yes]
Justification: The abstract and introduction (Section 1 summarize the contributions which are principally an elucidation of a stochastic SOS hierarchy (an optimization procedure), proof (with listed assumptions) giving quantitative rate of convergence, and numerical experiments. The full proof, numerical plots and tables, as well as code necessary to replicate results are all provided.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize other settings.
It is fine to include aspirational goals as motivation as long as is clear that these goals attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: See the Discussion section (Section where we highlight the non-scalability of all algorithms (including ours) that require optimization over the SDP cone, and dis- cuss potential solutions, some of which we rectify in this paper by using sparsification assumptions such as a "cluster basis". Despite a quantitative convergence proof with clear rate, we also highlight the assumptions made, including the fact that we restrict to 1 --periodic trigonometric polynomials over [0, 1|" x Ω with & compact (Section 2.3.2 and Appendix A.5.2..
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested few datasets with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they seale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimess.
31

---

While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, 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 impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: full proof of our core result provided in the supplementary material (Appendix A.5.2 along with 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 cross- referenced.
All assumptions should be clearly stated or referenced 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 provide 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.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: All information regarding experiments 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 the figures, tables, and other results we report.
Guidelines:
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 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 specific model and empirical evaluation, it may be necessary to either make 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., 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 submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm.
32

---

(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is new model (e.g.. a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome 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., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions 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.
Please see the NeurIPS code and data submission guidelines (https://nips... cc/
gublic/guidess htte///oooSSmmssoooooooooooo
for more details.
While we encourage the release of code and data, we understand that this might be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed run reproduce the results. See the NeurlPS code and data submission guidelines (https /nips.cc/paut httpppplwlpppppppppppppppppaaaep Stw.:iddee//oeSSdmbmeeeeeeeeee
ssioioolicyyy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, 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 the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can provided either with the code, appendix, or as supplemental material.
7. Experiment Statistical Significance
33

---

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, report confidence intervals solution recovery metrics for SNL which are computed over multiple runs. In particular, we report robust mean and I-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, confi- dence 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, 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.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2- 2-sigma error bar than state that they have 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 range (e.g. negative eтror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Specific details regarding the compute resources used detailed in Section as well 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 compute required for each of the individual experimental runs well 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 into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics http::/://auurp..ccoooo/////ttta
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.
34

---

If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. there a special consid- eration due to laws or regulations their jurisdiction).
10. 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 no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact 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., deploymennt 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 tied to particular applications, let alone deployments. However, if there is 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.
there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
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)?
Question:
Answer: [NA]
Justification: Our contributions are primarily algorithmic and do not have high risk for misuse.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have 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 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 best faith effort.
12. Licenses for existing assets
35

---

Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: All open-source source software packages such as scientific computing libraries and optimization libraries were cited.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset used and, possible, include a URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g.. website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperssi thcode. com datasets has curated licenses for some datasets. Their licensing guide can help determine the license dataset.
For existing datasets that are re- re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. If this information is available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: No from scratch.
used, all data necessary for experiments generated
Guidelines:
The answer NA means that the paper does release new assets.
Researchers should communicate the details the dataset/codeemoodll as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether how consent was obtained from people whose asset used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: Human subjects were not involved.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
36

---

15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/teview based the requirements of your country institution) were obtained?
Answer: [NA]
Justification: Human subjects were not involved.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
37


Reviews:
[{"review_id": "ASLP6F7J62", "reviewer": "Reviewer_MyRu", "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).", "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": "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.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- One strength lies in the convergence results on S-SOS, which is a vital generalization of the standard SOS. The proof and arguments in the Appendix seem non-trivial and thorough.\n- The cluster-basis hierarchy is novel and practically relevant. Although SOS has good theoretical performance, the SDPs could be expensive to solve in real-world engineering problems. The proposed hierarchy reduces the sizes of the SDPs that need to be solved. \n- Applications to the SNL problem complement the theoretical results in convergence and demonstrate the advantage of leveraging the S-SOS technique over standard Monte Carlo methods.", "weaknesses": "- 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. \n\n- 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.\n\n- 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.\n\n- 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer MyRu for their detailed comments.\n\nThe reviewer brings up a good point about presenting the assumptions in a subsection. We will set aside a section in the appendix for this purpose.\n\nAs for novelty in our proof techniques or lack thereof, it was not obvious to us that one could take the existing results and just apply them directly to a new framework. The way a proof is presented also matters and we hope that our detailed presentation in the appendix is instructive for similar paths taken in the future.\n\nWe believe our global response covers the remaining points brought up in the strengths/weaknesses sections, which center on the limitations of our theory and numerics:\n- Limitations of trigonometric polynomial assumption and the possibility of extending to more general polynomial families\n- Expanded cluster basis hierarchy discussion\n- Small size of experiments\n\nThank you again for your consideration.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "k3tvtQYBVG", "reviewer": "Reviewer_pRMv", "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.\n\nThe 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.\n\nRegarding 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.\n\nIn 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.", "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.\n\nCould 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": "NA", "rating": 3, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "-- Modification of the SOS hierarchy to incorporate noisy parametric problems \n-- Error bounds under assumptions", "weaknesses": "-- Relevance \n-- Novelty", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer pMRv for their detailed comments.\n\nConcerning our work’s novelty and relevance, we believe that our proposal is one of specific interest to the ML/AI community here.\n\nApplied mathematicians have studied SDPs and convex optimization for some time, while ML/AI scientists are more interested in non-convex/scalable methods as they can be applied to deep learning and other topics of recent interest. We believe fruitful work is to be done at the intersection of these two. To give just two examples of papers that fuse ideas across polynomial optimization, SOS methods, and neural networks, consider:\n- Agrawal et al 2019 [[link](http://arxiv.org/abs/1910.12430)]: Introduces differentiable optimization layers in neural networks, which combines ideas from practical convex optimization and produces a new tool to be used in deep learning.\n- Jaini et al 2019 [[link](http://arxiv.org/abs/1905.02325)]: A framework for high-dimensional density estimation inspired by sum-of-squares polynomials.\n\nOur work assumes that both $f(x, \\omega)$ and the lower-bounding function $c(\\omega)$ are polynomials. Once we characterize how the SOS hierarchy works in this simple case, we can generalize to the cases where both functions are more complex, perhaps even parameterized by neural networks. Insights from understanding the rate of convergence and the speed of solving a given finite-degree SDP can help us design algorithms in the general case where the functions are more complex.\n\nConsider that we propose a general primal/dual program whose solution provides a tight lower bound. To solve it in practice, one may make the polynomial assumption and then truncate at finite degree, which gives rise to the SDPs we analyze here. Another direction may be to solve the entire problem in the neural network setting: parameterize the basis function $m(x, \\omega)$ with a neural network, use a smooth interpolating function class for the lower bound $c(\\omega)$, and solve the program with an iterative gradient method. Finally, inspired by Agrawal et al 2019, imagine using the lower bound (or extracted solutions) obtained as output from our solved SDP as input to another layer in a neural network.\n\nThe reviewer may not be impressed by the size of the problems that can be solved here, but note our discussion on scalability and size in the global response. We will call the reviewer’s attention to the fact that noiseless v. noisy SNL are two very different problems - in the former we can solve SNL instances of ~100k sensors, in the latter it becomes much more difficult. The resulting output of our program actually solves for the global optimum, which gives sensor positions for every possible configuration of parameter $\\omega$. Finally, we only use general-purpose SDP solvers in this work. Many additional methods can be explored, such as exploiting more specialized solvers, low-rank structure, sparsity, and subspaces of the PSD cone. We also expect that leveraging GPUs and other methods of acceleration can immensely increase the size and practicality of the framework we propose.\n\nOur work is also relevant to many problems commonly seen in the AI + Science area. Recent works often seek to learn some kind of molecular potential with deep learning, to do something like “deep molecular dynamics”. The goal here is to have some energy function that depends on some external parameters (e.g. noise). We believe that present approaches lack a principled understanding, and that grounding it in a polynomial optimization setting and sum-of-squares methods can lead to generally useful results, if not a new way of thinking about such problems. To such an end, our paper also proposes a cluster basis hierarchy, introduced and used to scale the size of problems we can solve.\n\nThank you again for your consideration.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "33S4vORgAG", "reviewer": "Reviewer_Cxce", "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.", "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?\n\n2. 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. \n\n3. Is there any intuition as to how the result stated in Prop. 2.2. can be extended to the multivariate case?", "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.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. This proposes S-SOS, a novel heirarchy of relaxations for solving stochastic polynomial optimization problems.\n2. The paper provides asymptotic convergence results for the sequence of relaxations - that is, as the degree of the SOS program increases, the error between the lower bound obtained via the relaxation and the true optimal value reduces to 0. In particular, the finding that by using piecewise constant approximation of $c(\\omega)$ achieves the $1/s^2$ convergence rate, outperforming the polynomial parameterization of $c(\\omega)$, is quite interesting and surprising to me. \n3. The paper is written quite clearly, with no major typos that I could find. \n4. The proofs appear to be correct, with a couple of small typos such as in equation 6.", "weaknesses": "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)$.\n2. 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer Cxce for their detailed comments.\n\nWith regards to the typos identified, we appreciate the flagging of these and will correct them in the next version. The reviewer is indeed correct in that Equation (6) should have the plus sign.\n\nWe thank the reviewer for flagging additional approaches to polynomial optimization and will include those citations. In particular, we find the DSOS and SDSOS papers of particular interest. It seems that the core idea is to consider further relaxing the SDP to a SOCP or LP problem, which constrains the PSD cone further. We didn’t consider further relaxations as our goal was to propose a framework using the tightest possible yet still solvable bound (hence the SDP relaxation). We expect that LP and SOCP relaxations would lead to more efficient and scalable algorithms, though at the expense of accuracy.\n\nWe 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$. This is true in all cases in that setting, but once we pass to the parameterized/stochastic setting of S-SOS, this no longer holds true generically (or even in most practical cases). \n\nJust as an example, consider the case discussed in Section 3.1, where $f(x, \\omega) = (x - \\omega)^2 + \\omega^2 x^2$. The cost is a degree-4 polynomial in two variables but we show no finite-degree polynomial will achieve the exact cost $c^*(\\omega) = \\omega^4 / (1 + \\omega^2)$. Instead, we see that it can be well-approximated by a polynomial of degree $s$ with $s$ small. And indeed, in many practical scenarios, we anticipate that this is the case. Why this is true (that simple polynomials in general have non-polynomial best lower bounds but may be well-approximated by low-degree polynomials nonetheless) is out of scope in our analysis but it is a very interesting practical fact.\n\nThe reviewer has a good point that it is possible to extend the idea behind Prop 2.2 to the multivariate case. We neglected this in our work, but it appears that as in the 1D case of Prop 2.2, one can create a piecewise grid approximation to the lower-bounding function and do a SOS decomposition at each point. However, the constants that prefix the rates will get worse and have a dependence on the spatial dimension, i.e. we will have a rate of the form $\\frac{C \\sqrt{d}}{s_p}$ where $d$ is the dimension of the noise space $\\Omega$ and $s_p$ is the number of grid points per dimension. One should also observe that the number of degrees of freedom for this approximation scales exponentially in the space dimension. We shall adapt the paper accordingly.\n\nThank you for your consideration.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "g5m3PxUkQG", "reviewer": "Reviewer_Dyr8", "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.", "questions": "Could you explain the application of the joint+marginal method as a starting point? \n\nCould you comment on the rates of convergence of the  joint+marginal method?  \n\nCould 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": "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.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "-- The convergence rate results are neat, although limited to 1-periodic trigonometric polynomials over compact sets. \n\n-- The manuscript cites a number of interesting recent preprints (such as https://theses.hal.science/LAAS-POP/hal-04201167v1).", "weaknesses": "-- 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. \n\n-- It is not clear how the objective of the heuristically extracted solution (A.7.6) differs from the objective function of the relaxation.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank reviewer Dyr8 for their detailed comments. In particular, we appreciate all the external citations they flagged.\n\nThe reviewer highlights that the sensor network localization (SNL) problem can be reduced to trivial sizes citing work solving instances of 20k-100k sensors, while our numeric experiments only have up to 15 sensors. A few things need to be clarified on this front.\n\nIn the cited paper (Krislock and Wolkowicz 2018, https://arxiv.org/abs/1002.0013), 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. Krislock and Wolkowicz use this reduction to solve SNL problems of 10k-100k sensors, but they develop a highly specialized algorithm that does not use any SDP solvers, as per the abstract.\n\nNoiseless SNL dramatically simplifies the problem. The intuition here is that in the noiseless setting, localizing even a small number of sensors near an anchor will propagate the correctly localized positions to sensors nearby. If one can find small groups of sensors that can be well-oriented w.r.t. each other, one can pursue localization of these groups in parallel and 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 (c.f. https://link.springer.com/article/10.1007/s11276-007-0034-9, https://epubs.siam.org/doi/10.1137/100792366).\n\nThe solution extraction in A.7.6 uses the moment matrix to find the mean/variance of the sensor positions. We did not clarify this in sufficient detail in the paper but Lasserre 2010 [[link](https://epubs.siam.org/doi/10.1137/090759240)] discusses and proves the convergence of a sequence of semidefinite relaxations (here, our dual). Our dual (their primal) has moment vector solutions converging to that of a probability distribution encoding all globally optimal solutions. Thus, as long as we solve the SDP with large enough degree $s$ to high accuracy, our extracted solutions will be “close” to the true solution of the infinite-degree relaxation. In our cluster basis hierarchy, we do not prove the convergence of the same moment vector. Numerics seem to suggest that it behaves similarly as the full basis hierarchy, but this is presently a limitation of our theory, and we hope to find a similar convergence proof.\n\nDue to the limited space available, we could not add much detail on the joint+marginal method. The joint+marginal method is essentially our dual semifinite program (our Eq 3 and 4). In that paper, no quantitative convergence rates are provided, although later works (such as the ones we reference) have them.\n\nThank you for bringing up the Bellon et al paper. In Bellon et al, a parametric SDP has a trajectory of solutions (a sequence of PSD matrices $X$ parameterized by $t$) and seeks to characterize the geometry of the trajectory $(X, t)$. In our work, the primal SDP has as its objective $\\int c(\\omega) d\\nu(\\omega)$ where $c(\\omega)$ is guaranteed to be a lower-bound to the cost function $f(x, \\omega)$ and $\\nu(\\omega)$ is the density of $\\omega$. When we solve the finite-degree primal SDP (Eq 2), we obtain a PSD matrix $W$ that describes the “non-negative part” of the function $f(x, \\omega)$ for all $\\omega$. As formulated we cannot directly compare our results with that of Bellon et al.\n\nWe instead consider the related problems\n$$ \\max c s.t. f(x, \\omega) - c = m(x)^T W(\\omega) m(x), W(\\omega) \\succcurlyeq 0 $$\nwith dual\n$$ \\min \\langle M(\\omega), K \\rangle s.t. \\langle K_a, M(\\omega) \\rangle = 0, M(\\omega) \\succcurlyeq 0 $$\nwhich describes a series of $(c, W, M)$ values that depend on the parameter $\\omega \\in \\Omega$. We must further restrict to an open interval $\\Omega = (\\omega_i, \\omega_f) \\subseteq \\mathbb{R}$. Now we have a dual SDP that matches the parametric SDP of Bellon et al.\n\nLet’s first analyze the simple function $f(x, \\omega) = (x - \\omega)^2 + (\\omega x)^2$ with $\\omega \\in [-1, 1]$. The parametric SDP amounts to finding a SOS approximation to the function at every $\\omega$. We can analytically solve for the lower-bounding path $c(\\omega)$ and the SOS residual $m(x)^T W(\\omega) m(x)$. We find:\n$$ m(x)^T W(\\omega) m(x) = f(x, \\omega) - c(\\omega) = \\frac{1}{\\omega^2 + 1} (\\omega^2 x - \\omega + x)^2 $$\nThe path $(W, c)$ taken as $\\omega$ varies in $[-1, 1]$ sweeps out a smooth curve and we conclude that all points $(W, c)$ are regular for $\\omega \\in [-1, 1]$.\n\nIt is more difficult for SNL, although we have a few observations to start with:\n- In the noiseless setting (parameter $\\omega = 0$) there exists a unique configuration of the sensor positions.\n- In the noisy setting, for small $\\omega$, the recovered sensor positions will be “close” to the correct unique sensor positions.\n- Consider 3 anchors (A, B, C) and 1 sensor (X) in 2D SNL. If the 3 anchors are not degenerate, then when $\\omega=0$ we have the sensor positioned at the intersection of three circles, each having radius equal to the noiseless distance. If we perturb one of the observed distances (X-C) along a line, we note that there are two perturbing values that would lead to loss = 0, leading to consistent distances. This corresponds to the circle intersecting the other two at two points, with mirror symmetry about the segment intersecting A and B.\n\nThis last observation helps us see that we should have a discontinuity when the perturbation places the sensor at the segment intersecting A and B, and thus we should have multi-valued $W, M$ on either side of that point.\n\nIt’s possible that we misunderstood the work somehow. Regardless, it’s quite interesting and we shall incorporate this into our paper.\n\nThank you for your consideration.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "FEQTlqdIEq", "reviewer": "Reviewer_iKyk", "summary": "This paper is concerned with the following natural generalization of polynomial optizimation: \nGiven 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.\n\nGiven the computational hardness of the problem, the authors focus on the sum-of-squares relaxation (of a given degree $s$).\n\nFor 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.\nFinally, they complement these results with experimental evaluations (in which they apply certain SDP sparsification techniques to improve the scalability).", "questions": "What can the authors say for other classes of polynomials?", "limitations": "I encourage the authors to further discuss the limitations of their theoretical results.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "I find the question  investigated in this paper and the two main theorems (2.1 and 2.2) valuable and non-obvious.\n\nThe \"Cluster basis hierarchy\" is also a non-trivial original contribution.\n\nFinally, I believe the experiments included in the paper will be of interest to the wider NeurIPS' audience.", "weaknesses": "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.\n\nThe 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.\n\nThe 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer iKyk for their detailed comments.\n\nWe believe that the reviewer’s primary concerns are addressed by our global response:\n- Limitations of trigonometric polynomial assumption and the possibility of extending to more general polynomial families\n- Expanded cluster basis hierarchy discussion\n- Small size of experiments\n\nRegarding the core ideas behind our results and a more careful discussion of the theorems, limitations, and insight - we hope to improve and expand our discussion to make the journey easier and more fruitful.\n\nThank you again for your consideration.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We saw the following points come up in multiple reviews so we thought it would make sense to address it as a global response.

The reviewers chiefly seemed to have concerns about the following:
- The limitations of our theoretical assumptions (particularly the trigonometric polynomial one) and the possibility of extending to more general polynomial families
- Clarity, motivation, and further discussion (e.g. more on limitations, cluster basis hierarchy)
- The limitations of our numerical experiments, such as the relevance of sensor network localization and the small size of our experiments

**Limitations of theory.**
With respect to our theory, we agree wholeheartedly that the specialization to 1-periodic trigonometric polynomials is a severe limitation. To review, our work proves quantitative convergence of the S-SOS hierarchy in the setting of (i) a compact domain $X, \Omega$ and (ii) 1-periodic trigonometric polynomials. Periodic trigonometric polynomials are a natural choice to take advantage of Fourier convergence results on a compact domain. If the function $f(x, \omega)$ can be assumed to be only of interest for $(x, \omega)$ in some compact set, then we can rescale the compact set to be a 1-periodic compact domain and apply the 1-periodic trigonometric polynomial results. To use other polynomial families, one can apply a substitution argument, i.e. any result for the trigonometric polynomial hierarchy leads to a matching result for regular polynomials (2.2 in Bach Rudi 2023, link). We regret that we didn’t mention this and will add this to our paper.

**Improved clarity and discussion, esp. on the cluster basis hierarchy.**
The reviewers also mentioned that our paper could have been written more clearly and with much more discussion and elaboration, including on the motivations and intuitions behind the theory as well as the cluster basis hierarchy. Due to the limited space of the Neurips venue, we squeezed what we could into the main text and put everything else in the appendix. We hope to greatly expand the supplement as we continue our revisions.

As for the cluster basis hierarchy, we agree that more discussion is needed — as it is a core technique behind our goal of scaling up the S-SOS approach, particularly when the cost function is well-structured. We will expand our discussion of it in the supplement and hope to cover it in more detail in future work. Unfortunately, it is difficult to prove similar theoretical results for the cluster basis. But here we can also provide numerical support for the convergence of the cluster basis.

**Limitations of numerics: relevance and small scale.**
Finally, with respect to our numerics, the reviewers commented that the core problem (sensor network localization, SNL) is interesting but not of relevance to the broader Neurips community and we only demonstrated it on $N=15$ sensors, a small number compared to the state-of-the-art. We commented in response to another reviewer but will reproduce our response here.

*In the cited paper (Krislock and Wolkowicz 2018, https://arxiv.org/abs/1002.0013), 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. Krislock and Wolkowicz use this reduction to solve SNL problems of 10k-100k sensors, but they develop a highly specialized algorithm that does not use any SDP solvers, as per the abstract.*

*Noiseless SNL dramatically simplifies the problem. The intuition here being that 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 (c.f. https://link.springer.com/article/10.1007/s11276-007-0034-9, https://epubs.siam.org/doi/10.1137/100792366).*

Note also that the $N=15$ number is deceptive. **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. This means that we get a solution for the sensor positions for every possible configuration of noise, via the probability distribution $\mu(x, \omega)$.

As for the relevance to Neurips: sensor network localization is a problem that this conference is quite unfamiliar with. It is an old problem in polynomial optimization but generally challenging in the noisy setting. Once lifted into the stochastic setting, we find that our S-SOS framework is a natural fit for this but also many other possible problems in "AI + Science". We believe that taking this approach can be quite fruitful and we hope to see this line of work through.

We want to thank all reviewers for their consideration.

Meta-review: 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.
Decision: Accept (poster)