An Efficient Framework for Global Non-Convex Polynomial Optimization over the Hypercube]An Efficient Framework for Global Non-Convex Polynomial Optimization over the Hypercube

[1]Pierre-David<EMAIL_ADDRESS>1]Dalton<EMAIL_ADDRESS>1]Matthew<EMAIL_ADDRESS>1]M. Harper<EMAIL_ADDRESS>
[1]Qualcomm AI Research[Qualcomm AI Research is an initiative of Qualcomm Technologies, Inc.], 5775 Morehouse Dr, San Diego, 92121, CA, USA

We present a novel efficient theoretical and numerical framework for solving global non-convex polynomial optimization problems. We analytically demonstrate that such problems can be efficiently reformulated using a non-linear objective over a convex set; further, these reformulated problems possess no spurious local minima (i.e., every local minimum is a global minimum). We introduce an algorithm for solving these resulting problems using the augmented Lagrangian and the method of Burer and Monteiro [7]. We show through numerical experiments that polynomial scaling in dimension and degree is achievable for computing the optimal value and location of previously intractable global polynomial optimization problems in high dimension.

§ INTRODUCTION

For polynomials $p(x)$ in $D \in \mathbb{N}$ dimensions of degree at most $d \in \mathbb{N}$, consider  global non-convex polynomial optimization problems of the form,

\begin{align}
\label{eq:problemcontinuous}
\min_{[-1,1]^D} & \;\; p(x).
\end{align}

Based on an evolution of the ideas proposed by Lasserre in [16], we have developed a novel theoretical and numerical framework for efficiently solving problems in the form of (<ref>); specifically, such problems can be efficiently reformulated as problems involving a non-linear objective over a convex cone of semi-definite matrices in $O(D\cdot d^2)$ dimensions.  Further, our reformulation possesses no spurious local minima (i.e., every local minimum is a global minimum).  Our numerical framework for solving the non-convex polynomial optimization problems from (<ref>) is based on the augmented Lagrangian and the “Burer-Monteiro Method” [7] for expressing semi-definite constraints efficiently.

Our contributions can be summarized as:

* Efficient reformulations of global polynomial optimization problems: a reformulation of (<ref>) is introduced as a problem involving a canonical non-linear objective over a cone of semi-definite matrices. Both formulations share the same optimal value (shown in Theorem <ref>) and location (discussed in Section <ref>).

* Efficient characterization of product measures: a complete and efficient characterization is constructed of the set of product measures supported over $[-1,1]^D$ through their moments using semi-definite constraints (shown in Proposition <ref>).

* Efficient numerical algorithms: an efficient algorithmic framework and numerical implementation prototype is provided for solving our efficient reformulation of  (<ref>) with low computational complexity (discussed in Section <ref>).  Through examples, we demonstrate polynomial scaling on previously intractable global optimization problems (shown for dimensions up to $250$ in Section <ref>).

The remainder of this paper is structured as follows: Previous work is discussed in Section <ref>.  Section <ref> introduces notation, theory, and our novel reformulation of (<ref>). Section <ref> discusses algorithmic and numerical aspects associated with our specific computational treatment of the reformulation. Section <ref> presents numerical results demonstrating correctness and performance of our framework. Section <ref> presents conclusions and a discussion of future work. Proofs can be found in Appendix <ref>.

§.§ Previous Work

Our reformulation and efficient framework for solving problems of the form (<ref>) is inspired by the polynomial optimization framework introduced by
Lasserre [16], who proposed reformulating (<ref>) as a Semi-Definite Program (SDP). Specifically, instead of evaluating polynomials directly, the optimization can be done by integrating over specific multidimensional measures, the moments of which then replace the monomials of the original polynomial in the objective. The primary advantage of this approach is that the resulting convex problem can be solved more efficiently than the original formulation. However, rather than provide an exact reformulation, [16] results in a convex relaxation, where the solution provides a lower bound on the optimal minimum.  The quality of this bound depends on the degree of the Sums-Of-Squares (SOS) [17, 26] polynomials used to generate the relaxation such that better relaxations increase the size of the resulting SDP, adversely affecting scaling. In [16, 24] it was demonstrated that any sequence of problems of increasing SOS degree (also known as the “Lasserre Hierarchy”) will eventually converge to an exact solution in finite time. Although this hierarchy has been found to converge relatively rapidly on some specific problems (see, e.g., [10]), theoretical estimates indicate that super-exponential sizes are required in the worst case [9].

Proofs of convergence and degree estimates for the approach of [16] rely on results by Pitunar [27])(“Putinar Positivstellensatz”) and Nie and Schweighofer [23]. The results of [27, 23] address the problem of expressing positive polynomials over compact semi-algebraic subsets of Euclidean spaces through SOS, or SOS-like, representations of finite degree. They demonstrate that this can be done in any dimension; however, their estimates demand a super-exponential SOS degree of the original problem in both dimension and degree, resulting in intractable estimates for the cost of solving polynomial optimization problems using [16]'s framework in the general case.

Our framework and results overcome the drawbacks of [16] by restricting the space of feasible positive polynomials (and measures) to combinations of products of positive one-dimensional polynomials (and product measures). In fact, Proposition <ref> demonstrates that such polynomials can be represented efficiently through an SOS representation (specifically, at a cost linear in the dimension and degree of the original problem).  In this restricted case, the moment problem[i.e., the problem of extending elements of the dual space of polynomials to positive product measures.] can be characterized and solved efficiently, requiring the knowledge of only linearly many moments in the dimension and degree of the original problem.  Such conclusions are much stronger than the general case considered in [16] and highlight the efficiency of our novel reformulation.

From an algorithmic standpoint, it is well-known that SDPs are tractable [22, 1]. Indeed, being convex problems, SDPs may be solved in polynomial time [6, 3]. Despite such characteristics, however, the treatment of semi-definite constraints by existing numerical solvers often carry a prohibitively high cost in practice [15, 9]. For example, using an eigenvalue decomposition for enforcing positive-definiteness may carry a cost of $O(N^3)$ where $N$ is the number of primal variables [13]; this computational cost is prohibitively high for many interesting problems arising from the Lasserre framework. In practice, semi-definite relaxations are often solved using interior-point methods (e.g., [2, 14, 22]). Modern implementations, however, are limited to problems of relatively small sizes due to their associated computational costs [21], especially for increasing interest in problems that are large-scale and potentially ill-posed. [12]

One recent approach that reduces this cost is that or Burer and Monteiro (BM, [7]) who proposed to replace semi-definite constraints of the form $A \succeq 0$ by equality constraints of the form: $A = R \, R^T$. Boumal [4] showed that when the rank of $R$ is large enough (i.e., roughly the square root of the number of constraints in dual form), the resulting (non-linear) BM problem possess the property that every local minimum is a global minimum (second-order stationarity is, in fact, necessary and sufficient). When the rank of $R$ is appropriately small, this leads to linear scaling estimates for pragmatically solving SDPs.  Another recent approach for mitigating the restrictions of large SDPs is presented by Dressler et. al in  [11].  The authors present a Positivstellensatz based on sums of nonnegative circuit polynomials (SONC), which coupled with geometric programming (GP) allows for computing faster lower bounds than Lassere's SOS/SDP approach on certain classes of problems.  Further, for a polynomial that is sparse and SONC, there is a guarantee of a certificate with the same sparsity.

While the approach of Dressler et. al [11] is a promising new approach, leveraging BM [7] allows for potentially greater generality for larger classes of polynomial optimization. Polynomials optimization problems on compact sets can be tackled exactly using the BM approach within our framework to enforce semi-definite constraints.  By reformulating problems of the form in Equation (<ref>), we can solve the resulting formulations using BM and the augmented Lagrangian.  More details are discussed in Sections <ref>- <ref>.

§ THEORY

In Section <ref>, we introduce the notation used throughout the remainder of the paper.  In Section <ref>, we present our reformulation and theoretical results demonstrating the correctness and efficiency of our proposed approach.

§.§ Notation

We represent the dimension of Problem <ref> as $D \in \mathbb{N}$, and use $d$ to represent the problem degree (i.e., the degree of the objective polynomial). A $D$-dimensional multi-index $n=(n_1, n_2, ..., n_D)$ is an element of $\mathbb{N}^D$, where $\lvert n \rvert = \sum_{i=1}^D n_i$ refers to the multi-index degree. Unless otherwise stated, a basis of monomials of the form,

\begin{equation}
x^n = \prod_{i=1}^D x_i^{n_i}
\end{equation}

is used to represent elements of $\mathbb{P}_{D,d}$, the vector space of polynomials of degree $d$ in dimension $D$, i.e.,

\begin{equation}
\label{eq:P_Dd}
\mathbb{P}_{D,d} := \left \{ \sum_{\lvert n \rvert \leq d} \, p_n \, x^n : \, p_n \in \mathbb{R} \; \forall \; n \in \mathbb{N}^D \, \mathrm{s.t.} \, \lvert n \rvert \leq d \right \},
\end{equation}

where the short-hand notation

\begin{equation*}
\sum_{\lvert n \rvert \leq d} a_n= \sum_{ n \in \{ n : \lvert n \rvert \leq d \} } a_{(n_1, n_2, ..., n_D)}
\end{equation*}

represents summations over all multi-indices of degree less than or equal to $d$.  In particular, note that $\mathbb{P}_{D,d}$ has dimension

\begin{equation}
\# \{ n \in \mathbb{N}^D \, : \, \lvert n \rvert \leq d \} = \binom{D+d}{d} \leq \min \left ( D^d, d^D \right ).
\end{equation}

Hence, for a fixed degree $d$, the dimension of the vector space of polynomials space $\mathbb{P}_{D,d}$ grows polynomially in the problem dimension $D$. The notation $\mathbb{P}_{D,\infty}$ is used to represent polynomials of any finite degree. The support of a polynomial $p \in \mathbb{P}_{D,\infty}$ corresponds to the locations of its nonzero coefficients; i.e.,

\begin{equation}
\mathrm{supp}(p) := \{ n \in \mathbb{N}^D \, : \, p_n \not = 0 \},
\end{equation}

and the cardinality (number of elements) of the support is denoted by

$$N_p := \# \, \mathrm{supp}(p).$$

The set of regular Borel measures over $\mathbb{R}^D$ ([8]) is denoted as $\mathbb{B}_{D}$, and
the support of a measure $\mu(\cdot) \in \mathbb{B}_D$ is designated by $\mathrm{supp}(\mu(\cdot))$, representing the smallest set such that $\mu \left ( \mathbb{R}^D \backslash \mathrm{supp}(\mu(\cdot)) \right )  = 0$.
For a fixed multi-index $n$, a $D$-dimensional Borel measure  $\mathbb{B}_{D}$ has an $n^{th}$ moment defined to be the multi-dimensional Lebesgue integral,

\begin{equation}
\mu_n := \int_{\mathbb{R}^D} x^n \, \mathrm{d} \mu(x).
\end{equation}

A regular product measure is a regular Borel measures of the form,

\begin{equation}
\mu(\cdot) = \prod_{i=1}^D  \mu_i(\cdot) ,
\end{equation}

where each factor $\mu_i(\cdot)$ belongs to $\mathbb{B}_1$. In this context, the moments of a product measure up to degree $d$ correspond to a $D$-tuple of real vectors $( \mu_1, \mu_2, ..., \mu_D ) \in \mathbb{R}^{(2d+1) \times D}$, where each element is defined as

\begin{equation}
\mu_{i,n_i} := \int_{\mathbb{R} } x_i^{n_i} \, \mathrm{d} \mu_i (x_i).
\end{equation}

The vectors $\mu_i$ refer to the vector of moments of the 1D measure $\mu_i$.  With a single sub-index, $\mu_i \in \mathbb{R}^{d+1}$ represents the vector of moments, and with two sub-indices, $\mu_{i,n_i} \in \mathbb{R}$  represents scalar moments (vector elements).  A set of $L$ moments of product measures is represented by $\mu = \{\mu_i^{(l)}\}_{l,i} \in \mathbb{R}^{(d+1) \times D \times L}$. When discussing a measure, parentheses are employed (e.g., $\mu_i(\cdot)$), whereas the absence of parentheses indicates that we are discussing vectors.

Given a sequence of moments $\mu_i$ of a 1D measure, the degree-$d$ 1D moment matrix (or just moment matrix when the context is clear) associated with the one-dimensional polynomial $g(x_i) \in \mathbb{P}_{1,d^\prime}$, $0 \leq d^\prime < \infty$ is defined to be the $(d +1) \times (d +1) $ matrix $\mathcal{M}_d(\mu_i, g(x_i))$ with elements,

\begin{equation}
[ \mathcal{M}_d(\mu_i, g(x)) ]_{m_i,n_i} = \int g(x_i) \, x_i^{m_i + n_i} \, \mathrm{d} \mu_i(x_i) = \sum_{k_i=0}^{d^\prime} g_{k_i} \, \mu_{k_i+m_i+n_i},
\end{equation}

for $0 \leq m_i, n_i \leq d $. When $g(\cdot ) \equiv 1$, we use the short-hand notation,

\begin{equation}
\mathcal{M}_d(\mu_i; 1) = \mathcal{M}_d(\mu_i).
\end{equation}

The construction of such matrices requires the knowledge of only the first $2 d+ d^\prime + 1$ moments of $\mu_i (\cdot) $, i.e., $\{ \mu_{i,n_i}\}_{n_i=0}^{2d + d^\prime} $. As shown in Proposition <ref>, positive definiteness of the matrix $\mathcal{M}_d(\mu_i)$ is necessary and sufficient to ensure the existence of positive measures with such moments. This, however, is true only in the 1D case, and generally false in higher dimensions.

§.§ Efficient Reformulation

Given some integer $L$, our reformulation of (<ref>) takes the form,

\begin{align}
\label{eq:reformulationcontinuous}
\min_{  \mu \in \mathbb{R}^{(d+1)\times D \times L} }&\;\; \sum_{n \in \mathrm{supp}(p)} \, p_n \, \phi_n (\mu) \nonumber \\
\mathrm{subject\, to} &  \;\; \mathcal{M}_d (\mu^{(l)}_i) \succeq 0,  \\
&\;\; \mathcal{M}_{d-1} (\mu^{(l)}_i; \, 1-x_i^2) \succeq 0,  \;&  i=1,..., D, l=1,..., L \nonumber \\
& \;\; \mu^{(l)}_{i,0} \geq 0, \nonumber \\
&\;\; \phi_{(0,...,0)}  (\mu)  = 1, \nonumber  \\
& \;\;  \phi_n (\mu)  = \sum_{l=1}^L  \prod_{i=1}^D \mu^{(l)}_{i,n_i}. \nonumber
\end{align}

This reformulation corresponds to a nonlinear objective subject to (convex) semi-definite constraints. Further, because the feasible set lies in a space of size $L \cdot D \cdot (d + 1) $, which is linear in both degree and dimension (see Table <ref> for more details on computational costs), our reformulation is efficient.

Following Proposition <ref>, the problem may also be interpreted as the minimization of a linear objective over the set consisting in the normalized (with mass equal to $1$) convex combinations of $L$ positive product measures. This interpretation is key in demonstrating how the reformulated problem <ref> is in fact equivalent to Problem <ref> as shown in Theorem <ref>-<ref>.

Type     Expression      Type      Size     Eval. Cost
O     $\sum_{n \in \mathrm{supp}(p)} \, p_n \left ( \sum_{l=1}^L  \prod_{i=1}^D \mu^{(l)}_{i,n_i} \right )$     Scalar     $1$     $ O\left ( N_p \, L \, D \right )$
C     $ \mathcal{M}_d \left ( \mu^{(l)}_i \right ) \succeq 0$     PSD     $(d+1) \times (d+1)$     $  O\left ( L \, D \, d^2 \right )$
C     $ \mathcal{M}_{d-1} \left (\mu^{(l)}_i; \, 1-x_i^2 \right )   \succeq 0 $     PSD     $ d \times d$     $  O\left ( L \, D \, d^2 \right )$
C     $  \sum_{l=1}^L \prod_{i=1}^D \mu^{(l)}_{i,0} = 1 $     Scalar     $1$     $  O\left ( L \, d \right )$

Size and evaluation costs of each component (Objective (O) or Constraints (C) Category Type) of the proposed reformulation.

The key advantages of this reformulation can be found below in Theorems <ref> and <ref>; Theorem <ref> demonstrates the equivalence between the original and reformulated problem, and Theorem <ref> demonstrates that the reformulated problem may be solved efficiently using local descent techniques only.

In order to prove Theorem <ref>, we rely on Proposition <ref> below, which demonstrates the possibility of efficiently (i.e., at a polynomial cost) extending linear functionals belonging to the dual space $\mathbb{P}_{D,d}^*$ when the latter are in product form. This is in contrast to the general case which rules out such extension, a consequence related to the fact that positive polynomials cannot generally be written as combination of sums-of-squares having a sufficiently small degree (see, e.g., Motzkin polynomial [29]). This underlies key differences between our framework and that introduced by Lasserre (e.g., [18]). Indeed, the lack of existence of an efficient representation in Lasserre's formulation leads to reformulations (i.e., semi-definite relaxations) that provide a lower bound. To obtain a solution to the original problem under Lasserre's framework, it is necessary to consider a hierarchy of problems (the commonly-called Lasserre hierarchy) of increasing (and potentially exponential) size. As shown below, our proposed framework does not suffer from such drawbacks; in fact, (1) the dimension of our reformulated problem is fixed and of size at most polynomial in dimension and degree, (2) the optimal value obtained matches that of the original problem exactly, and (3) no hierachy is required. The proof for Proposition <ref> is provided in the Appendix <ref>.

Let $D,d \in \mathbb{N}$ and $( \mu_1, \mu_2, ..., \mu_D ) \in \mathbb{R}^{(d+1) \times D}$ be such that for each $i = 1, ..., D$, $\mu_{i,0} = 1$, and

\begin{align*}
\mathcal{M}_d(\mu_i) &\succeq 0,\\
\mathcal{M}_{d-1}(\mu_i;  1-x_i^2) &\succeq 0 .
\end{align*}

Then, there exists a regular Borel product measure,

\begin{equation*}
\mu( \cdot ) := \prod_{i=1}^D  \mu_i(\cdot),
\end{equation*}

supported over $[-1,1]^D$ such that $\mu \left ( [-1,1]^D \right )  = 1$, and

\begin{equation*}
\int_{[-1,1]^D} \, x^n \, \left (  \prod_{i=1}^D \mu_i(x_i) \right ) \, \mathrm{d} x =  \prod_{i=1}^D \mu_{i,n_i}
\end{equation*}

for all multi-index $n \in \mathbb{N}^D$ such that $0 \leq n_i  \leq d$ for all $i$.

Proposition <ref> specifically states that the feasible set of Problem <ref> corresponds to the set of normalized convex combination of $L$ product measures supported over the hypercube $[-1,1]^D$. This, in turn, implies the following results in Theorem <ref> and Theorem <ref> (proofs provided in Appendix <ref>).

Problem <ref> and Problem <ref> share the same global minimum value.

Let $ \left \{  \mu_i^{(0)},\mu_i^{(1)}  \right \}_{i=1}^{D} \in \mathbb{R}^{(d+1) \times D \times 2 }  $ be a feasible point of Problem <ref> corresponding to a local minimum, and let $\alpha^{(j)} = \prod_{i=1}^D \mu^{(j)}_{i,0}$

for $j \in \{0,1\}$. Then, such a point corresponds to a global minimum if and only if $\alpha^{(0)} \not = 0$ and the point,

$$\left \{  \frac{ \prod_{i=1}^D\mu^{(0)}_{i,n_i} }{\alpha^{(0)}}  \right \}_{\lvert n \rvert \leq d},$$

is a local minimum, or $\alpha^{(1)} \not = 0$ and the point,

$$\left \{  \frac{ \prod_{i=1}^D\mu^{(1)}_{i,n_i} }{\alpha^{(1)}}  \right \}_{\lvert n \rvert \leq d},$$

is a local minimum.

Theorem <ref> implies it is possible to solve Problem <ref> by solving Problem <ref>, while Theorem <ref> implies that every local minimum of Problem <ref> is a global minimum. Indeed, using Theorem <ref>, a descent algorithm can be devised for computing a local minimum of Problem <ref> (see Section <ref>):

* Following the computation of a local minimum, at least one of $\alpha^{(0)}$ or $\alpha^{(1)}$ must be nonzero by construction.  Depending on the nonzero status of $\alpha^{(i)}$, $i\in\{0,1\}$, one may substitute the corresponding point specified in Theorem <ref> and relaunch the descent algorithm:

* If the algorithm fails to descend, it implies that the new point is a local minimum and one concludes that the original point is a global minimum following Theorem Theorem <ref>.

* If the algorithm does successfully descend, another local minimum with lower value will be computed, and the above process is repeated until converging to a global minimum.

Finally, even though Theorem <ref> is stated using a combination of only two ($L=2$) product measures, the conclusion still holds for any $L \geq 2$; the choice of $L=2$ is made for the purpose of simplifying the analysis; however, the analysis holds for larger $L > 2$. In practice, it is common to set $2 \leq L \leq 20$ (see Section <ref>). Adding more product measures increases the number of degrees of freedom and often improves numerical convergence.

§ ALGORITHM

Theorem <ref> demonstrates the feasibility of recovering the optimal value of Problem <ref> by solving Problem <ref>. In addition, Theorem <ref> dictates that such optimal values can be reached using local descent methods only. Based on these theoretical results and techniques, we have implemented a numerical algorithm for solving the reformulated problem. The remainder of this section provides details about our algorithm.[There exists many alternative algorithms for solving Problem <ref>. Our algorithm is a single instance and represents our attempt at devising an efficient way to solve the reformulated problem. It is by no means unique.] Further, technical details specific to our implementation can be found in Section <ref>, and numerical results are discussed in Section <ref>.

Our algorithm is an extension of the method of Burer and Monteiro (BM), presented in [7]. The main idea behind the BM method is to replace positive-definite constraints of the form $A \succeq 0$ with constraints of the form: $A = R R^T$ for some matrix $R$ of appropriate size. Problem <ref> can thus be written as [In order to simplify our implementation, we fix $\mu^{(l)}_{i,0} = 1$ for $i\geq 2$, while allowing $\mu^{(l)}_{0,0}$ to remain variable.]

\begin{align}
\label{eq:reformulationcontinuous_BM}
\min_{  \mu \in \mathbb{R}^{  (d+1) \times D \times L}   }&\;\; \sum_{n \in \mathrm{supp}(p)} \, p_n \, \phi_n (\mu),\nonumber \\
\mathrm{subject\, to} & \;\; \mathcal{M}_d(\mu^{(l)}_i) = R^{(l)}_i R^{(l)^T}_i ,  \\
&\;\; \mathcal{M}_{d-1}(\mu^{(l)}_i; \, 1-x_i^2) = S^{(l)}_i S^{(l)^T}_i , \; & i=1,..., D, l=1,..., L \nonumber\\
&\;\; R_{i}^{(l)} \in \mathbb{R}^{(d+1) \times (d+1)} , \, S_{i}^{(l)}  \in \mathbb{R}^{ d \times d }, \nonumber\\
& \; \; \mu_{i,0} \geq 0, \nonumber \\
& \;\;  \phi_{(0,...,0)} (\mu) = 1, \nonumber \\
& \;\;   \phi_n (\mu) =  \sum_{l=1}^L  \prod_{i=1}^D \mu^{(l)}_{i,n_i} .\nonumber
\end{align}

Upon re-expressing the (element-wise) matrix equality constraints as a sequence of $K = L\cdot D ( (d+1)^2 + d^2 + 1) + 1$ scalar equality constraints, we obtain a generalized BM problem reformulation of the form,

\begin{align}
\label{eq:reformulationcontinuous_BM_general}
\min_{ \mu \in \mathbb{R}^{  (d+1) \times D \times L}   }&\;\; \sum_{n \in \mathrm{supp}(p)} \, p_n \, \phi_n(\mu) \\
%\left ( \sum_{l=1}^L  \prod_{i=1}^D \mu^{(l)}_{i,n_i} \right )\\
\mathrm{subject\, to} &\;\; c_k\left ( \mu, \{R_i^{(l)}, S_i^{(l)}\} \right ) = 0 , &\;\;  \, k=1,..., K \nonumber
\end{align}

We consider this form in Equation (<ref>) for numerical optimization. Our numerical implementation constructs the augmented Lagrangian as proposed in [7]:

\begin{align}
\label{eq:augmented_lagrangian}
\mathcal{L}_\gamma (\mu, \{R_i^{(l)}, S_i^{(l)}\}; \lambda) &:=  \sum_{n \in \mathrm{supp}(p)} \, p_n \, \phi_n(\mu) +  \sum_{k=1}^K \lambda_k \cdot c_k\left ( \mu, \{R_i^{(l)}, S_i^{(l)}\} \right ) \\
&+ \frac{\gamma}{2} \, \sum_{k=1}^K  c^2_k\left ( \mu, \{R^{(l)}, S^{(l)}\} \right ), \nonumber
\end{align}

where $\gamma > 0$ is a penalty parameter. This is similar to the Lagrangian introduced in [7], with the difference that the objective exhibits non-linearities.

We solve this problem numerically using a technique similar to one proposed by Boyd et al. [5], which involves finding a saddle-point (i.e., a simultaneous minimum and maximum in the primal and dual variables respectively) to the augmented Lagrangian.  That is, our approach focuses on finding $\mu^*$ and $\lambda^*$, such that

\begin{equation}
\mathcal{L}_\gamma (\mu^*, \{R_i^{(l)^*}, S_i^{(l)^*}\}; \lambda^*) = \max_\lambda \,\min_{\mu, \{R_i^{(l)}, S_i^{(l)}\}}   \mathcal{L}_\gamma (\mu, \{R_i^{(l)}, S_i^{(l)}\} ; \lambda).
\end{equation}

Specifically, our implementation alternatively solves the primal and dual problems until convergence is achieved. The primal problem takes the form

\begin{equation}
\min_{\mu, \{R_i^{(l)}, S_i^{(l)}\}} \, \mathcal{L}_\gamma (\mu, \{R^{(l)}, S^{(l)}\}; \lambda),
\end{equation}

where the Lagrange multipliers (dual variables) $\lambda$ are fixed. This is an unconstrained optimization problem with an objective that is bounded below (thanks to the penalty term).  This can be solved using a descent technique; for this purpose, our algorithm leverages L-BFGS [19].  Meanwhile, the dual problem takes the form

\begin{equation}
\max_{\lambda} \, \mathcal{L}_\gamma (\mu, \{R^{(l)}, S^{(l)}\}; \lambda),
\end{equation}

where the primal variables $\mu$,  $R^{(l)}$  and $S^{(l)}$ are fixed. Because this problem is generally unbounded at non-optimal points, we use an update rule proposed by [5]:

\begin{eqnarray}
\lambda_{\mathrm{next}} &=& \lambda_{\mathrm{current}} + \gamma \, \nabla_{\lambda} \,\mathcal{L}_\gamma (\mu, \{R_i^{(l)}, S_i^{(l)}\}; \lambda)\\
&=& \lambda_{\mathrm{current}} + \gamma \, \sum_{k=1}^K  c^2_k(\mu, \{R_i^{(l)}, S_i^{(l)}\})
\nonumber
\end{eqnarray}

This process is summarized in Algorithm <ref>. We have implemented this algorithm in the C/C++ language. Section <ref> presents numerical results.

Alternating descent/ascent algorithm for computing a saddle-point of the augmented Lagrangian

[1]
AlternatingSaddlePoint$\mu_0, \lambda_0; \gamma$
$\mu \gets \mu_0$
$\lambda \gets \lambda_0$
Not Converged
$\mu \gets \mathrm{argmin}_{\mu} \, \mathcal{L}_\gamma (\mu; \lambda)$Using L-BFGS
$\lambda \gets \lambda + \gamma \, \nabla_{\lambda} \,\mathcal{L}_\gamma (\mu; \lambda) $

return $(\mu^*,\lambda^*)$

§.§ Optimal Location

When the global minimum $x^* \in [-1,1]^D$ of Problem <ref> is unique, it can be shown that the solution of the reformulated problem (Problem <ref>) must consists in the (redundant) convex combinations of delta measures centered at the point $x^*$. Under these conditions, it is found that,

\begin{equation}
\frac{\mu^{(l)}_{i,n_i} }{\prod_{i=1}^D \mu_{i,0}^{(l)} } = \int_{-1}^1 x_i^{n_i} \,  \delta(x_i - x_i^*) \, \mathrm{d} x_i = (x_i^*)^{n_i},
\end{equation}

for every fixed index $l$ such that $ \forall i, \mu_{i,0}^{(l)} \not = 0 $ for every $i$.  In this special case, the above scheme does not merely provide a framework for finding the global optimal value, but it also provides a simple way of computing a globally-optimal location. Indeed, it suffices to pick any $l$ such the above condition holds and then compute

\begin{equation}
x_i^* = \frac{\mu^{(l)}_{i,1} }{ \prod_{i=1}^D  \mu_{i,0}^{(l)} },\quad{   }\forall i = 1, ..., D.
\end{equation}

This approach is summarized in Algorithm <ref>.

Optimal location extraction (assuming unique global minimum)

[1]
OptimalLocation$\mu$Input is output of Algorithm <ref>
$l^* \gets \mathrm{argmax}_{l=1,...,L} \, \prod_{i=1}^D \mu_{i,0}^{(l)} $
$i = 1, ..., D$
$ x_i^* \gets \frac{\mu^{(l^*)}_{i,1} }{\prod_{i=1}^D \mu_{i,0}^{(l^*)} } $
return $x^*$Global minimum location

When the global minimum is not unique (or not known to be so), a small amount of noise is added to the objective, making it unique with high probability.

§.§ Further Implementation details

This section provides greater details and implementation technicalities, based on the high-level algorithmic information provided in Section <ref>.

§.§.§ Constraint simplification

As implied from Equation (<ref>), our reformulation dictates that the semi-definite constraints should be expressed as

\begin{align*}
&\mathcal{M}_d(\mu^{(l)}_i) \succeq 0,  \\
& \mathcal{M}_{d-1}(\mu^{(l)}_i; \, 1-x_i^2) \succeq 0
\end{align*}

for $l=1, ..., L$ and $i=1,...,D$. Within our implementation, however, we use the equivalent constraints,

\begin{align*}
&\mathcal{M}_d(\mu^{(l)}_i) \succeq 0, \; -1 \leq \mu^{(l)}_{i,n_i} \leq 1,
\end{align*}

for $i=1,...,D $ and $0 \leq n_i \leq d$. This reduces computational complexity and memory costs by a factor of roughly $2$. The scalar inequality constraints are enforced using the BM method as described in Section <ref>. In this case, a BM factor $R_k R_k^T$ is degenerate and corresponds to a squared scalar. For instance, the constraint $\mu^{(l)}_{i,n_i} \leq 1$ becomes

$$1-\mu^{(l)}_{i,n_i}  = (t^{(l)}_{i,n_i})^2, \;\; t^{(l)}_{i,n_i} \in \mathbb{R}.$$

§.§.§ Numerical stability

The aforementioned constraints are sufficient to ensure the moments of the product measures are bounded within $[-1,1]$; i.e.,

\begin{equation}
\label{eq:redundant_constraint}
-1 \leq   \prod_{i=1}^D \mu^{(l)}_{i,n_i} \leq 1
\end{equation}

for $l=1,...,L$ and $n \in \mathrm{supp}(p)$. Throughout the computations, however, it may occur that such constraints are violated, resulting in large values for the product moments and potential numerical instabilities. This is particularly dire in higher dimensions where small violations may compound to create large values through multiplication.

In order to mitigate this potential numerical instability issue and improve the numerical behavior, our implementation  explicitly enforces the (redundant) constraints found in Equation (<ref>) by adding them to the formulation found in Equation (<ref>).

§.§.§ Polynomial basis

Our choice of monomials as a basis for the polynomials in Section <ref> is made to facilitate the exposition; however, any polynomial basis may have been used. In our implementation, we leverage a Chebyshev basis $\{T_{n_i}(x_i) \}$ (e.g.,[20]).  This choice offers better numerical stability without significantly increasing computational costs. The only differences are that the 1D moments are defined as,

\begin{equation}
\mu_{i,n_i} := \int_{-1}^1 \, T_{n_i}(x_i) \, \mathrm{d} \mu_i (x_i),
\end{equation}

and the moment matrices $ \mathcal{M}_d(\mu_i)$ have the form (Hankel + Toeplitz),

\begin{equation}
[ \mathcal{M}_d(\mu_i) ]_{m_i,n_i} = \int_{\mathbb{R}} \, T_{m_i}(x) \,  T_{n_i} (x_i) \, \mathrm{d} \mu_i(x_i) =  \frac{1}{2} \, \mu_{m_i+n_i} + \frac{1}{2} \, \mu_{\lvert m_i - n_i \rvert}.
\end{equation}

This is a consequence of the fact that Chebyshev polynomials satisfy the identity,

\begin{equation}
T_{m_i}(x_i) \,  T_{n_i}(x_i) = T_{m_i + n_i}(x_i) + T_{\lvert m_i - n_i \rvert} (x_i).
\end{equation}

§.§.§ Initialization

Our current implementation uses a simple initialization scheme (no warm start) for the primal and dual variables of the augmented Lagrangian (Equation (<ref>)) which can be described as follows:

* Primal Initialization: each scalar element of the primal variables, $\mu , \, \{R_i^{(l)}, S_i^{(l)}\}$,
is initialized randomly using a distribution $U\left ( [-1,1] \right )$ uniform over $[-1,1]$.

* Dual Initialization: the dual variables ($\lambda$) are initialized to be zero ($0$).

§.§.§ L-BFGS and Line Search

As previously discussed, we use L-BFGS [19] to compute descent directions (${\bf g}$) at each iteration and solve the primal problem. Iterates are then computed using a line search with backtracking and a Wolfe condition; i.e., the scheme computes the next point so it corresponds to the smallest $k$, such that

\begin{equation}
f({\bf x}_{\mathrm{current}} + 0.5^k \, {\bf d}) \leq  f({\bf x}_{\mathrm{current}}) + 0.5^k \, \beta \, ( \nabla f({\bf x}_{\mathrm{current}}) \cdot {\bf g })
\end{equation}

for some parameter $0 < \beta < 1$ (referred to as the Wolfe factor). The order of approximation of the Hessian typically ranges between $20$ and $200$. We further can leverage a mechanism by which the algorithm reverts to a gradient descent direction whenever the line search fails using the direction ${\bf g}$ computed using L-BFGS.

§.§.§ Stopping criteria

We use two (2) types of stopping criteria both for the inner (based on L-BFGS)  and outer (saddle-point) iterations:

* A criterion based on the size of the relative gradient;

* A criterion based on the absolute change in value between the previous and current iterates.

We also enforce feasibility; i.e.,

$$ \lvert c_k (x) \rvert < \epsilon_{f},$$

for some stopping parameter $0 < \epsilon_f < 1$ (absolute feasibility) as a part of the stopping criterion for the outer iteration.

§ NUMERICAL RESULTS

This section presents results demonstrating correctness and performance of our framework and algorithm for optimizing polynomials over the hypercube $[-1,1]^D$. For this purpose, we construct classes of highly non-convex functions with known global optimal value and location (see Figure <ref> and <ref>). The resulting functions possess characteristics that make global optimization intractable using traditional techniques, namely, non-convexity, exponentially many local minima and exponentially small (in relative volume) basin of attraction as a function of dimension $D$.

The numerical experiments were carried out on a Dell 1U server with 2x AMD EPYC 7452 32-core (128 system threads) and 256 GB DRAM and were limited to the use of a single (1) core. Timing represents wall-time for the whole computations, excluding I/O and the time needed to setup the problem (i.e. pre-computations). Relative error between a computed value ${\bf x}_{comp}$ and a known value ${\bf x}_{exact}$ refers to the quantity,

\begin{equation}
\frac{\lvert \lvert {\bf x}_{comp} - {\bf x}_{exact} \rvert \rvert_2 }{ 10^{-15} + \lvert \lvert {\bf x}_{exact} \rvert \rvert_2 }
\end{equation}

where $\lvert \lvert \cdot \rvert \rvert_2$ is the $\ell^2$-norm.

§.§ Numerical Example 1

The first example is derived from a class of problems that consists of objective polynomials of the form,

\begin{align*}
f_D (x) &= \frac{1}{D} \, \sum_{i=1}^D T_2(x_i) - \prod_{i=1}^D T_8(x_i)
=  \frac{1}{D} \, \sum_{i=1}^D (2 x_i^2-1)  - \prod_{i=1}^D \cos(8 \, \cos^{-1}(x_i)).
\end{align*}

Figure <ref> depicts $f_2(\cdot)$ in two dimensions (i.e., $D=2$).

Top view

3-D View

Two-dimensional (2-D) rendering of the  non-convex function $f_2(x)$ from Example 1 used for numerical experiments.  Note the exponential number of local minima. Global optimization of this function is intractable using traditional techniques.

For every fixed dimension $D$, the function $f_D (x)$ possess a unique global minimum at the origin (i.e., ${\bf 0} = (0,0,...,0)$) with value equal to

$$f_D^* = f_D({\bf 0}) = -2.0.$$

Each function exhibits more than $3^D$ local minima, and the basin of attraction of the global minimum has volume smaller than $\frac{1}{5^D}$. Together, these characteristics imply that the use of any descent technique is bound to fail. Indeed, any uniformly initialized descent technique will succeed in finding the global minimum with probability
at most $ \frac{1}{5^D}$. In dimension $D=10$, this corresponds to a probability of approximately one in $10^{7}$; furthermore, in dimension $D=100$, this probability corresponds to approximately one in $10^{70}$, making the problem intractable for descent techniques in general.

Numerical experiments for this example were performed using the algorithm described in Section <ref> together with the parameters shown in Table <ref>. With such a choice of parameters, we expect a relative error smaller than $10^{-2}$. The relative error for the optimal value and the optimal location, as well as the time required before reaching the stopping criterion are shown in Figure <ref> and Figure <ref> respectively for dimension $D$ from $1$ to $250$.

Parameters used for performing numerical experiments related to Example 1.
4*Formulation     Size of moment matrices ($d \times d$)     $8$

Penalty parameter $(\gamma)$       $8.0$

Burer-Monteiro rank     $9$ (full)

Number product measures ($L$)     $4$
3*Solver (KKT)     Stopping criterion (relative gradient)     $10^{-2}$

Stopping criterion (absolute value)     $10^{-4}$

Stopping criterion (absolute feasibility)     $10^{-2}$
4*Solver (Primal/L-BFGS)     Approximation order       $100$

Wolfe Beta factor ($\beta$)     $0.3$

Stopping criterion (relative gradient)     $ 10^{-3} $

Stopping criterion (absolute value)     $10^{-4}$

Figure <ref> shows that our algorithm succeeds in identifying both the global minimum value and the global optimal location within the expected level of accuracy in every problem instance up to $250$ dimensions.

Example 1 Relative error vs problem dimension ($D$) for dimension $1$ to $250$. Top: relative error associated with computed optimal value, Bottom: relative error associated with computed optimal location. The error lies below the expected upper bound (black dotted line) of $10^{-2}$ in every case.

As shown in Figure <ref>. our algorithm further achieves the desired accuracies in a computationally-efficient manner, exhibiting a highly-advantageous quadratic scaling of the time-to-solution as a function of dimension (i.e. $O(D^2)$ scaling). This is in line with expectations since the support of $f_D$ has $N_{f_D} = D+1$ elements, implying that the most expensive component of the algorithm (the evaluation of the objective) carries a cost on the order of $O(LD^2)$ (see Table <ref>).

Example 1 computation (wall) time (s) vs problem dimension ($D$). Our implementation exhibit a polynomial (quadratic; (black line)) scaling on problems previously considered intractable.

§.§ Numerical Example 2

The second numerical example comes from a class of functions of the form,

\begin{align*}
g_D (x) &= \frac{1}{D} \, \sum_{i=1}^D T_4(x_i) + \left ( \frac{1}{D} \, \sum_{i=1}^D  T_1(x_i)\right )^3 \\
&= \frac{1}{D} \, \sum_{i=1}^D \left ( 8x_i^4 - 8x_i^2 + 1 \right ) + \left ( \frac{1}{D} \, \sum_{i=1}^D  x_i \right )^3.
\end{align*}

Top view

3-D View

Two-dimensional (2-D) rendering of the  non-convex function $g_2(x)$ from Example 2 used for numerical experiments. Note the exponential number of local minima as well as the asymmetry. Global optimization of this function is intractable using traditional techniques.

As with Example 1, the functions of Example 2 are also highly non-convex and have exponentially many local mimina ($2^D$ in dimension $D$). Every function possesses a unique minimum with value approximately equal to $-1.3911$ located near $(-0.75553, ..., -0.75553)$; i.e.,

\begin{equation}
g^*_D := g_D(-0.75553,..., -0.75553) \approx -1.3911.
\end{equation}

Figure <ref> depicts $g_2(x)$ in $D=2$ dimensions. By contrast with Example 1, Example 2 is asymmetric about the origin. Further, the cost of evaluating the objective in this case is significantly higher; the functions $g_D(\cdot)$ each have $N_{g_D} = O(D^3)$ non-zero coefficients compared with $O(D)$ for Example 1.  As discussed below, this results in greater computational complexity.  As in Example 1, experiments were performed using the Algorithm described in Section <ref> with parameters detailed in Table <ref>. Given greater computational costs, we consider dimensions from $1$ to $45$. As before, a relative error of $10^{-2}$ is expected using such parameters.

Parameters used for performing numerical experiments related to Numerical Example 2
4*Formulation     Size of moment matrices ($d \times d$)     $4$

Penalty parameter $(\gamma)$       $10.0$

Burer-Monteiro rank     $5$ (full)

Number product measures ($L$)     $6$
3*Solver (KKT)     Stopping criterion (relative gradient)     $10^{-2}$

Stopping criterion (absolute value)     $10^{-3}$

Stopping criterion (absolute feasibility)     $10^{-2}$
4*Solver (Primal/L-BFGS)     Approximation order       $40$

Wolfe Beta factor ($\beta$)     $0.4$

Stopping criterion (relative gradient)     $ 10^{-2} $

Stopping criterion (absolute value)     $10^{-5}$

Figure <ref> shows how our method  accurately (i.e., within $10^{-2}$) computes the relative error in the optimal value and location for numerical example 2.  Figure <ref> depicts the time taken to compute the latter and exhibits a scaling of $O(D^4)$. As the driving computational component (the objective evaluation) carries a cost of $O(D^3)$ (resulting in a total cost of $O(D^4)$ according to Table <ref>), the results are in line with expectations.

Example 2 relative error vs problem dimension ($D$) for dimension $1$ to $45$. Top: optimal value; Bottom: optimal location. The error lies below the expected upper bound (dotted black line) of $10^{-2}$ in every case.

Computation (wall) time (s) vs problem dimension ($D$). Our implementation exhibits a polynomial (quartic, degree-$4$; black line) scaling on problems previously considered intractable.

§ CONCLUSION

In this paper, we have introduced a new framework for efficiently reformulating general polynomial optimization problems over the hypercube $[-1,1]^D$ as non-linear problems with no spurious local minima. We have rigorously demonstrated the correctness of our approach, and have introduced an efficient algorithm for the computational treatment of the reformulation.  We have further demonstrated the correctness, performance, and practicality of our implementation using difficult, previously-intractable polynomial optimization problems.

From a theoretical standpoint, future work will target the treatment of more complex polynomial (i.e., semi-algebraic) constraints. Performance improvements, both algorithmic and from a software optimization perspective (including parallelism), will also take a prominent role in the near future.

§ DATA AVAILABILITY STATEMENT

The data produced for the redaction of this paper will be made available under reasonable requests satisfying the policy of Qualcomm Technologies Inc. for the public release of internal research material.

§ PROOFS

Let $p \in \mathbb{P}_{1,d}$ be a 1D polynomial of degree $d$ positive on $[-1,1]$. Then,

\begin{equation}
p(x) = (f(x))^2 + (1-x^2) \, (g(x))^2
\end{equation}

for some polynomials $f(x)$ and $g(x)$ of degree at most $d$ and $d-1$ respectively.

Let $E$ be a real vector space, $F\subset E$ be a subspace and $K \subset E$ a convex cone. Assume further that: $E = K+F$. Let $\phi : F \rightarrow \mathbb{R}$ be a linear functional over $F$ that is positive over $K \cap F$; i.e., $\phi(y) \geq 0 \;\; \forall \;\; y \in K \cap F$  Then, there exists a linear functional,

$$\psi : E \rightarrow \mathbb{R},$$

such that

\begin{align}
\psi(x) &= \phi(x) \;\; \forall \;\; x \in F \\
\psi(y) &\geq 0 \;\; \forall \;\; y \in K.
\end{align}

[The statement is adapted from [28], Theorem IV.14, and the preceding remarks on the one-to-one correspondence of Baire measures and regular Borel measures.]

Let $X$ be a compact Hausdorff space and assume that $\psi : C (X) \rightarrow \mathbb{R}$  is a continuous linear functional over $C(X)$, the space of functions continuous on $X$. Assume further that $\psi(\cdot)$ is positive over $C(X)$; i.e.,

\begin{equation}
\psi(f) \geq 0 \;\;  \forall \;\; f \in C(X) \, : \, f(x) \geq 0 \;\;\forall \;\; x \in X.
\end{equation}

Then, there exists a unique regular Borel measure $\mu (\cdot)$, such that

\begin{equation}
\psi(f) = \int_X f(x) \, \mathrm{d} \mu(x).
\end{equation}

Suppose $T:V_1 \rightarrow V_2$ is a bounded linear transformation from a normed linear space $\langle V_1, \lvert \lvert \cdot \rvert \rvert_1\rangle$ to a complete normed linear space $\langle V_2, \lvert \lvert \cdot \rvert \rvert_2\rangle$. Then, $T$ can be uniquely extended to a bounded linear transformation (with the same bound), $\tilde{T}$, from the completion of $V_1$ to $\langle V_2, ||\cdot||_2\rangle$.

The polynomials are dense in $C[a,b]$, the space of real-valued functions continuous over the interval $[a,b]$, for any finite real number $a < b$.

Let $F := \mathbb{P}_{1,d} $ be the vector space of 1D polynomials of degree at most $d$, $E = \mathbb{P}_{1,\infty}:= \mathbb{R}[x]$, the space of 1D polynomials with real coefficients of any finite degree, and

\begin{equation}
K := \left \{ p \in  \mathbb{P}_{1,\infty} \, : \, p(x) \geq 0 \;\; \forall \;\; x \in [-1,1]  \right \}
\end{equation}

be the cone of 1D polynomials of any finite degree positive on $[-1,1]$. We first demonstrate that

$$ E = K+F.$$

Given $p \in E$, let $h(x) \equiv \min_{x \in [-1,1]} p(x) > -\infty$. This constant polynomial is well-defined since continuous functions reach their extrema over compact sets. Further, $h \in F$ as $F$ contains the constant polynomials by definition. Finally, the function

\begin{equation}
p(x) - h(x) = p(x) - \min_{x \in [-1,1]} p(x) \geq 0 \;\; \forall \;\; x \in [-1,1]
\end{equation}

implies that $p-h \in K$ by construction. Since this holds for every element $p\in E$, this demonstrates the claim.

Now, consider the following family of linear functionals over $F$: $\forall i=1,..., D$ and $q = \sum_{k=0}^d q_k \,x^k \in F$, define

\begin{equation}
\phi_i (q(x_i)) :=  \sum_{n_i=0}^d   q_{n_i}  \, \mu_{i,n_i}.
\end{equation}

Note in particular that: $ \phi_i (x_i^{n_i}) = \mu_{i,n_i}$. We claim that each functional is positive on $F \cap K$. Indeed, assume $q$ is positive; i.e., $q \in F \cap K$. Then, by Theorem <ref> there exists polynomials $f \in \mathbb{P}_{D,d}$ and $g \in \mathbb{P}_{D,d-1}$, such that

\begin{align}
\phi_i (q(x_i)) &= \phi_i \left (f^2(x_i) + (1-x_i^2) \, g^2(x_i) \right )\\
&= \sum_{m_i=0, n_i=0}^d f_{m_i} f_{n_i} \, \phi_i(x^{m_i+n_i}) + \sum_{m_i=0, n_i=0}^{d-1} g_{m_i} g_{n_i} \, \phi_i( (1-x_i^2) \, x^{m_i+n_i}) \\
&= \sum_{m_i=0, n_i=0}^d f_{m_i} f_{n_i} \, \mu_{i,m_i+n_i} + \sum_{m_i=0, n_i=0}^{d-1} g_{m_i} g_{n_i} \, \left  (\mu_{i,m_i+n_i} - \mu_{i,m_i+n_i+2} \right ) \\
&= f^T \, \mathcal{M}_d ( \mu_i ) \, f + g^T \,  \mathcal{M}_{d-1} ( \mu_i;(1-x_i^2) )   \, g \\
& \geq 0,
\end{align}

where the last inequality follows by assumption since $\mathcal{M}_d ( \mu_i ) \succeq 0$ and $\mathcal{M}_{d-1} ( \mu_i; \, 1-x_i^2 ) \succeq 0 $ for every $i$. Therefore, $\phi_i()$ are positive functionals on $F \cap K$ and thus satisfies the hypothesis of the Riesz extension theorem (Theorem <ref>). This implies that for each $\phi_i(\cdot)$ there exists an extension $\psi_i : E \rightarrow \mathbb{R} $ that is positive on $K$.

Recall that $E$ corresponds to $\mathbb{P}_{1,\infty}$, which is a dense subspace of the space $C[-1,1]$ of functions continuous on $[-1,1]$ by Theorem <ref> (Stone-Weierstrass). For each $i$, we can therefore invoke the B.L.T. theorem (Theorem <ref>) with $\psi_i(\cdot)$ a linear transformation from the normed space $\langle \mathbb{P}_{1,\infty}, \lvert \lvert \cdot \rvert \rvert_\infty \rangle$ to the complete normed space $\langle \mathbb{R}, \lvert \cdot \rvert \rangle$. In this way, we can assume w.l.o.g. that

$$\psi_i : C([-1,1]) \rightarrow \mathbb{R}$$

is a continuous linear functional defined over the entire space $C([-1,1])$ and positive over the cone of continuous functions positive on $[-1,1]$.   Under these conditions, the Riesz-Markov theorem (Theorem <ref>) implies the existence of unique regular Borel measures $\{\mu_i (\cdot)\}$ supported on $[-1,1]$, such that

\begin{equation}
\mu_{i,n_i} = \phi_i (x_i^{n_i}) = \int_{[-1,1]} x_i^{n_i} \, \mathrm{d} \mu_i (x_i),
\end{equation}

where $\mu_{i,0}  =1$. Finally, it follows as a consequence of Caratheodory theorem (see, e.g., [8], Theorem 5.1.3) that there exists a unique regular product measure

$$\mu(\cdot) = \prod_{i=1}^D \mu_i(\cdot)$$

with moments

\begin{equation}
\mu_n  =  \int_{[-1,1]^D} x^n \, \mathrm{d} \mu(x)   =  \prod_{i=1}^D \, \int_{[-1,1]} x_i^{n_i} \, \mathrm{d} \mu_i (x_i) = \prod_{i=1}^D \mu_{i,n_i}, \; 0 \leq n_i \leq d,
\end{equation}

and the constraints $\mu_{i,0} = 1 \, \forall \, i$ further imply that

\begin{equation}
\mu_n([-1,1]^D) = \prod_{i=1}^D \mu_{i}([-1,1]) = \prod_{i=1}^D \mu_{i,0} = 1,
\end{equation}

which proves the desired result.

We now introduce some quantities that will prove important for the remainder of the proof. For $0 < L \in \mathbb{N}$, define

and $\mathcal{F}_L$ and $\mathcal{S}_L$ to be

\begin{align}
\begin{split}
\mathcal{F}_L :=
\bigg\{   \bigg\{\mu_n:= {}& \sum_{l=1}^L \prod_{i=1}^D \mu_{i, n_i}^{(l)}  \bigg\}_{\lvert n \rvert \leq d}  \, : \, \{\mu_i^{(l)} \}_{i,l=1}^{D,L} \in \mathbb{R}^{D \times (d+1) \times L},\, \\
& \, \mathcal{M}_d(\mu^{(l)}_i), \, \mathcal{M}_{d-1}(\mu^{(l)}_i ; (1-x_i^2)) \succeq 0 \; \forall \, i,l, \, \mu_{0} = 1 \bigg\},
\end{split}
\end{align}

\begin{align}
\mathcal{S}_L &:= \left \{ \mu(\cdot) := \sum_{l=1}^L \prod_{i=1}^D  \mu_i^{(l)} (\cdot)  \, :  \,  \mu^{(l)}_i (\cdot) \in \mathbb{B}_1  \; \forall \; i,l, \, \mu\left( [-1,1]^D \right ) = 1\right \},
\end{align}

where as before, the symbol $\mathbb{B}_D$ represents the set of finite regular Borel measures supported over $[-1,1]^D$. With this notation, we may now proceed.

There exists a continuous surjective map,

$$\zeta :\mathcal{S}_L \rightarrow \mathcal{F}_L,$$

between $ \mathcal{S}_L$ and $ \mathcal{F}_L$ that satisfies

\begin{equation}
\left [  \zeta \left ( \mu(\cdot) \right ) \right ]_n =  \sum_{l=1}^L \prod_{i=1}^D \mu_{i, n_i}^{(l)}
\end{equation}

for every $  \lvert n \rvert  \leq d$.

Given an element

\begin{equation}
\mu(\cdot) := \sum_{l=1}^L \prod_{i=1}^D  \mu_i^{(l)} (\cdot)  \in \mathcal{S}_L,
\end{equation}

we define the map $\zeta(\cdot)$ through

\begin{equation}
\left [  \zeta \left ( \mu(\cdot) \right ) \right ]_n := \int_{[-1,1]^D}  \, x^{n} \, \mathrm{d} \mu (x) = \sum_{l=1}^L \prod_{i=1}^D  \int_{[-1,1]^D}  \, x_i^{n_i} \, \mathrm{d} \mu_i^{(l)} (x_i).
\end{equation}

First, we demonstrate that $\zeta \left ( \mu(\cdot) \right )$ belongs to $\mathcal{F}_L$. To do so, let

\begin{equation}
\mu_{i,n_i}^{(l)} := \int_{[-1,1]^D}  \, x_i^{n_i} \, \mathrm{d} \mu_i^{(l)} (x_i) = [\zeta \left ( \mu(\cdot) \right )]_n,
\end{equation}

and note that,

\begin{equation}
\sum_{l=1}^L \prod_{i=1}^D  \mu^{(l)}_{i,0} = \left [  \zeta \left ( \mu(\cdot) \right ) \right ]_0 =  \int \, \mathrm{d} \left ( \sum_{l=1}^L \prod_{i=1}^D     \mu_i^{(l)} (x_i) \right ) = \mu\left ( [-1,1]^D \right )  = 1.
\end{equation}

Thus, it suffices to shows that $\mathcal{M}_d(\mu^{(l)}_i)$ and $\mathcal{M}_{d-1}(\mu^{(l)}_i ;\, (1-x_i^2))$ are positive semi-definite for each $i=1,...,D$ and $l=1,...,L$. To see how this is so, note that following our definition,

\begin{align}
f^T \, \mathcal{M}_d(\mu^{(l)}_i) \, f &= \int_{-1}^1 \, f^2(x_i) \, \mathrm{d}\mu^{(l)}_i(x_i) \geq 0,\\
g^T \, \mathcal{M}_{d-1}(\mu^{(l)}_i; (1-x_i^2)) \, g &= \int_{-1}^1 \, (1-x_i^2) \, g^2(x_i) \, \mathrm{d}\mu^{(l)}_i(x_i) \geq 0,
\end{align}

for every $f\in\mathbb{P}_{1,d} $ and $g\in\mathbb{P}_{1,d-1} $, where positivity follows from the fact the integrands are non-negative over $[-1,1]$ and the measures are positive and supported over $[-1,1]$ by definition. Since this holds for every vector of coefficients $f$ and $g$, we conclude that all matrices are indeed positive semi-definite. This proves the claim.

Regarding surjectivity, consider an element $\mu \in \mathcal{F}_L$.
Following Proposition <ref>, for each vector of moments $\mu^{(l)}_{i}$ there exists of a regular Borel measure $\mu^{(l)}_i(\cdot)$ supported over $[-1,1]$, such that

\begin{equation}
\int_{[-1,1]^D} \, x_i^{n_i} \, \mathrm{d} \mu_i^{(l)} (x_i) = \mu^{(l)}_{i,n_i}
\end{equation}

for every $\lvert n \rvert  \leq d$ that also satisfies: $\sum_{l=1}^L \prod_{i=1}^D \mu^{(l)}_{i,0} = 1 $. The element,

\begin{equation}
\mu(\cdot) := \sum_{l=1}^L \prod_{i=1}^D  \mu_i^{(l)} (\cdot)
\end{equation}

therefore belongs to $\mathcal{S}_L$ and maps to $\mu$ through $\zeta(\cdot)$ as may be seen through substitution. Since this is true for every element of $\mathcal{F}_L$, this demonstrates surjectivity. Finally, continuity follows from the fact that given any two element $\mu(\cdot), \nu(\cdot) \in \mathcal{S}_L$, one finds that,

\begin{align}
&\left \lvert   \left [  \zeta \left ( \sum_{l=1}^L \prod_{i=1}^D  \mu_i^{(l)} (\cdot)   \right ) \right ]_n  -    \left [  \zeta \left ( \sum_{l=1}^L \prod_{i=1}^D  \nu_i^{(l)} (\cdot)   \right ) \right ]_n  \right \rvert  \\
&= \left \lvert \int_{[-1,1]^D} \, x^n  \, \mathrm{d}  \left ( \sum_{l=1}^L \prod_{i=1}^D  \mu_i^{(l)} (\cdot) - \sum_{l=1}^L \prod_{i=1}^D  \nu_i^{(l)} (\cdot) \right )  \right \rvert \\
&\leq \left  \lvert  \sum_{l=1}^L \prod_{i=1}^D  \mu_i^{(l)} (\cdot) - \sum_{l=1}^L \prod_{i=1}^D  \nu_i^{(l)} (\cdot) \right \rvert
\end{align}

by the triangle inequality and since $\lvert x^n \rvert \leq 1$ for all $x \in [-1,1]^D$, where $\lvert  \cdot \rvert$ indicates total variation (see, e.g., [8]).

Consider the problems,

\begin{equation}
\label{eq:problemmoments}
\min_{ \mu \in \mathcal{F}_L }   \sum_{ n \in \mathrm{supp}(p)} \, p_n \left ( \sum_{l=1}^L  \prod_{i=1}^D \mu^{(l)}_{i,n_i} \right )
\end{equation}

\begin{equation}
\label{eq:problemmeasure}
\min_{ \mu(\cdot) \in \mathcal{S}_L } \; \int p(x) \, \mathrm{d} \left (  \sum_{l=1}^L  \prod_{i=1}^D \mu^{(l)}_i (x_i) \right ).
\end{equation}

Then, the following holds:

* Both problem share the same global minimum value.

* $ \mu(\cdot) \in \mathcal{S}_L$ corresponds to a global minimum of Problem <ref> if and only if $ \zeta(\mu(\cdot)) \in \mathcal{F}_L$ corresponds to a global minimum of Problem <ref>.

* If $ \mu \in \mathcal{F}_L$ corresponds to a local minimum of Problem <ref>, then any element of the non-empty set $\zeta^{-1}(\mu) \in \mathcal{S}_L$ consists of local minima of Problem <ref>.

First, note that the objective functional satisfies,

\begin{align}
\label{eq:invfunc}
\int p(x) \, \mathrm{d} \left(\sum_{l=1}^L  \prod_{i=1}^D \mu^{(l)}_i(x_i) \right ) &=  \sum_{n \in \mathrm{supp}(p)} \, p_n \left (  \left [  \zeta \left(\sum_{l=1}^L  \prod_{i=1}^D \mu^{(l)}_i(\cdot) \right )  \right ]_{n} \right )
% \sum_{n \in \mathrm{supp}(p)} \, p_n \left ( \sum_{l=1}^L  \prod_{i=1}^D \mu^{(l)}_{i,n_i} \right ) \\
\end{align}

by definition. This implies that $ \mu^*(\cdot) \in \mathcal{S}_L$ is a global minimum of Problem <ref> if and only if $\mu^* := \zeta(\mu^*(\cdot))$ is be a global minimum of Problem <ref>. Indeed, if a point in $\mathcal{F}_L$ achieves a strictly lower value, then any point belonging to the pre-image of that point under $\zeta(\cdot)$  (which is non-empty by surjectivity. See Proposition <ref>) achieves a lower value and vice-versa. This is a contradiction and demonstrates the first two points.

For the third point, consider a local minimum $\mu \in \mathcal{F}_L $. Let $\mu(\cdot)\in \mathcal{S}_L $ be any element of the pre-image of this local minimum under $\zeta(\cdot)$; i.e., $\mu(\cdot) \in    \zeta^{-1} \left (  \mu \right )  $ which is non-empty thanks to surjectivity. Since $\mu$ is a local minimum, there exists an open neighborhood $\mathcal{N}$ of $\mu$ in $\mathcal{F}_L$ where the functional may only achieve values that are either greater than or equal to the current value. Consider the set,

\begin{equation}
\zeta^{-1} (\mathcal{N}) \subset \mathcal{S}_L
\end{equation}

which is open and contains $\zeta^{-1} \left (  \mu \right )$ since $\zeta(\cdot)$ is continuous and surjective. From the invariance of the functional value (Equation (<ref>)), it follows that every element of this set must lead to a value in Problem <ref> that is either greater than or equal to the current value, i.e., the set $\zeta^{-1} (\mu ) $  must consist of local mimina of Problem <ref>. This proves the third claim.

First, note that $\mathcal{F}_L$ is in a trivial bijective correspondence with the feasible set of Problem <ref>. Then, following Corollary <ref>, it suffices to consider the problem,

\begin{align}
\min_{\mu(\cdot) \in \mathcal{S}_L}\; \int_{[-1,1]^D} p(x) \, \mathrm{d} \left(\sum_{l=1}^L \, \prod_{i=1}^D \mu_i(x_i) \right )
\end{align}

which shares the same global optimal value. We claim that the optimal value $\pi^*$ of this problem (and thus that of Problem Problem <ref>) equals,

\begin{equation}
\min_{x \in [-1,1]^D} \, p(x)
\end{equation}

To see why, let $x^* \in [-1,1]^D$ be an optimal location (which exists since $p(x)$ is continuous and $[-1,1]^D$ is compact) and note that,

\begin{equation}
\prod_{i=1}^D \delta(\cdot - x_i^*)
\end{equation}

where $\delta(\cdot)$ is a point (delta) measure centered at the origin, belongs to $\mathcal{S}_L$. Therefore, we find that,

\begin{equation}
\pi^* \leq \int_{[-1,1]^D} p(x) \, \prod_{i=1}^D \delta(x_i - x_i^*)  \, \mathrm{d} x = p(x^*) =  \min_{x \in [-1,1]^D} \, p(x)
\end{equation}

Then, consider the fact that,

\begin{align}
\label{eq:equivmin_upperbound}
\int_{[-1,1]^D} &p(x) \, \mathrm{d} \left(\sum_{l=1}^L \prod_{i=1}^D \mu_i(x_i) \right ) \\
&\geq \left ( \min_{x \in [-1,1]^D} \, p(x)  \right ) \, \int_{[-1,1]^D} \mathrm{d} \left ( \sum_{l=1}^L \prod_{i=1}^D \mu^{(l)}_{i}(x_i) \right )  = \min_{x \in [-1,1]^D} \, p(x)
\end{align}

for every element of $\mathcal{S}_L$. Taking the infimum over $\mathcal{S}_L$ on both sides imply that,

\begin{equation}
\pi^* \geq \min_{x \in [-1,1]^D} \, p(x)
\end{equation}

which together with Equation (<ref>) imply that,

\begin{equation}
\pi^* = \min_{x \in [-1,1]^D} \, p(x)
\end{equation}

as claimed.

Consider an element of $\mathcal{S}_2$ of the form,

\begin{equation}
\mu^{(0)}(\cdot)  +  \mu^{(1)}(\cdot ) :=  \prod_{i=1}^D \mu_i^{(0)}(x_i)  +  \prod_{i=1}^D \mu_i^{(1)}(x_i)
\end{equation}

Assume such element is not a global minimum of the problem,

\begin{equation}
\min_{ \mu(\cdot) \in \mathcal{S}_2}\; \int p(x) \, \mathrm{d} \left( \prod_{i=1}^D \mu_i^{(0)}(x_i)  + \prod_{i=1}^D \mu_i^{(1)}(x_i) \right ).
\end{equation}

Then, if either 1) $\mu^{(0)}\left (  [-1,1]^D \right )  \in \{0,1\}$, or 2) $\mu^{(0)}\left (  [-1,1]^D \right ) \in (0,1)$ and,

\begin{equation}
\frac{1}{\mu^{(0)}([-1,1]^D) } \, \int p(x) \, \mathrm{d}  \mu^{(0)}(x) \not =  \frac{1}{\mu^{(1)}([-1,1]^D) } \, \int p(x) \, \mathrm{d}  \mu^{(1)}(x),
\end{equation}

there exists a path $\pi(t)$, $t \in [0,1]$ lying entirely in $\mathcal{S}_2$ and originating from this point along which $\int p(x) \, \mathrm{d} \pi(t)(x)$ monotonically decreases.

Consider a sub-optimal element $\mu(\cdot)$ of $\mathcal{S}_2$ of the form,

\begin{equation}
\mu(\cdot) =  \mu^{(0)}(\cdot) +  \mu^{(1)}(\cdot) = \prod_{i=1}^D \mu^{(0)}_i(\cdot) +   \prod_{i=1}^D \mu^{(1)}_i(\cdot).
\end{equation}

From Theorem <ref>, we know that the optimal value is given by,

\begin{equation}
\min_{x \in [-1,1]^D} p(x) = p(x^*)
\end{equation}

for some $x^* \in [-1,1]^D$. Therefore, sub-optimality implies that,

\begin{equation}
\int p(x) \, \mathrm{d} \mu(x) > p(x^*)
\end{equation}

We consider two cases: 1) the case where $\mu^{(0)}\left (  [-1,1]^D \right ) \in \{0,1\}$ and 2) the case where $\mu^{(0)}\left (  [-1,1]^D \right ) \in (0,1)$.

Case 1: $\mu^{(0)}\left (  [-1,1]^D \right ) \in \{0,1\}$

Assume w.l.o.g. that $\mu^{(0)}\left (  [-1,1]^D \right ) = 1$, so that $\mu(\cdot) = \prod_{i=1}^D \mu^{(0)}_i(\cdot)$.  Consider the measure,

\begin{equation}
\delta_{x^*}(\cdot) = \prod_{i=1}^D \delta(\cdot - x_i^*) .
\end{equation}

where $\delta(\cdot)$ is a point (delta) measure supported at the origin. Now, consider the parametrized path,

\begin{equation}
\pi(t) := \prod_{i=1}^D \pi^{(0)}_{i} (t) + \prod_{i=1}^D \pi^{(1)}_{i} (t)
\end{equation}

\begin{align}
\pi^{(0)}_{i} (t) &:= (1-t)^{1/D}  \,  \mu^{(0)}_{i}(\cdot) \\
\pi^{(1)}_{i} (t) &:= t^{1/D}\, \delta (\cdot - x_i^*) \\
\end{align}

for $t \in [0,1]$. In particular, note that this path starts at $ \pi(0) = \mu^{(0)}(\cdot) = \mu(\cdot)$. We further claim that this path lies entirely in $\mathcal{S}_2$ and that the objective monotonically decreases along it. To see why this is so, first note that $\delta(\cdot - x_i^*) \in \mathbb{B}_1$ for all $i$ by construction, and $\mu_i^{(0)}(\cdot) \in \mathbb{B}_1$ for all $i$ by assumption. Then, substitute to get,

\begin{align}
\pi(t) =   (1-t) \, \prod_{i=1}^D \mu_i^{(0)}(\cdot) + t \prod_{i=1}^D \delta(\cdot - x_i^*)   = (1-t) \, \mu^{(0)}_i (\cdot) + t \, \delta_{x^*} (\cdot)
\end{align}

and note that,

\begin{align}
& \pi(t)([-1,1]^D) = (1-t) \, \mu^{(0)} \left ( [-1,1]^D\right ) + t \, \delta_{x^*} \left ( [-1,1]^D \right ) =  (1-t) + t = 1,
\end{align}

for all $t \in [0,1]$, which implies that $\pi(t)$ belongs to $\mathcal{S}_2$ as claimed. For the second part of the claim, consider the functional value along the path, namely,

\begin{align}
\int p(x) \; \mathrm{d}\pi(t)(x) &=  (1-t)\,  \int p(x) \, \mathrm{d} \mu^{(0)}(x) + t\,  \int p(x)\, \delta_{x^*}(x) \, \mathrm{d}x \\
&= \int p(x) \, \mathrm{d} \mu^{(0)}(x) + t\, \left ( \int p(x) \, \delta_{x^*}(x) \, \mathrm{d}x  - \int p(x) \, \mathrm{d} \mu^{(0)}(x) \right )
\end{align}

This implies that,

\begin{equation}
\frac{d}{d\, t} \, \int p(x) \; \mathrm{d}\pi(t)(x) =  p(x^*) -  \int p(x) \, \mathrm{d} \mu^{(0)}(x) < 0
\end{equation}

since $\int p(x) \, \mathrm{d} \mu^{(0)}(x) > p(x^*)  $ by assumption of sub-optimality. This demonstrates monotonic decrease of the functional along the path as claimed.

Case 2: $\mu^{(0)}\left (  [-1,1]^D \right ) \in (0,1)$.

\begin{equation}
\alpha := \mu^{(0)}\left (  [-1,1]^D \right ).
\end{equation}

Under our hypothesis, one can assume w.l.o.g. that,

\begin{equation}
\label{eq:descent_case2}
\frac{1}{1-\alpha} \, \int p(x) \, \mathrm{d} \mu^{(1)}(x) < \frac{1}{\alpha} \, \int p(x) \, \mathrm{d} \mu^{(0)}(x) .
\end{equation}

Again, consider a parametrized path $\pi(t)$, $t\in[0,1]$ such that,

\begin{align}
\pi(t) :=   \prod_{i=1}^D \pi^{(0)}_{i} (t) +  \prod_{i=1}^D \pi^{(1)}_{i} (t)
\end{align}

\begin{align}
\pi^{(0)}_{i} (t) & =  \left ( 1 - t \right )^{1/D} \, \mu_{i}^{(0)}(\cdot) \\
\pi^{(1)}_{i} (t) &= \left( 1 + t \, \frac{\alpha}{1-\alpha} \right )^{1/D} \, \mu_{i}^{(1)} (\cdot)
\end{align}

for $t \in [0,1]$, which begins at $ \pi(0) = \mu^{(0)}(\cdot) +  \mu^{(1)}(\cdot) =  \mu(\cdot)$. We claim that this path lies entirely in $\mathcal{S}_2$ and that the objective monotically decreases along it. Indeed, $\mu_{i}^{(0)}(\cdot), \, \mu_{i}^{(1)}(\cdot) \in \mathbb{B}_1$ for every $i$ by assumption. In addition,

\begin{align}
\pi(t) &=   (1-t) \, \prod_{i=1}^D \mu^{(0)}_i (\cdot) +   \left (1 + t \, \frac{\alpha}{1-\alpha} \right )  \, \prod_{i=1}^D \mu_i^{(1)} (\cdot)\\
&= (1-t) \,  \mu^{(0)} (\cdot) +   \left (1 + t \, \frac{\alpha}{1-\alpha} \right )  \,  \mu^{(1)} (\cdot)
\end{align}

which implies that

\begin{align}
\pi(t)([-1,1]^D) &=  (1-t) \, \alpha  +  \left (1 + t\, \frac{\alpha}{1-\alpha} \right ) \, (1-\alpha)\\
&= \alpha +  (1-\alpha)  - t\, \alpha + t \alpha \\
&= 1
\end{align}

for all $t \in [0,1]$. This shows that $\pi(t)$ belongs to $\mathcal{S}_2$ as claimed. Finally, we have that the functional value along the path is,

\begin{align}
\int p(x) \, \mathrm{d} \pi (x)  &=  \int p(x) \, \mathrm{d} \left ( (1-t) \, \mu^{(0)}(x) +  \left(1+    t \,  \frac{\alpha}{1-\alpha} \right)\,  \mu^{(1)}(x) \right )\\
&=   \int p(x) \, \mathrm{d}  \mu^{(0)} (x) +   \int p(x) \, \mathrm{d} \mu^{(1)}  (x)  \\
&+ t\, \alpha \, \left ( \frac{1}{1-\alpha} \, \int_{\mathbb{R}^D } p(x) \, \mathrm{d}  \mu^{(1)}  (x) - \frac{1}{\alpha} \, \int_{\mathbb{R}^D } p(x) \, \mathrm{d}\mu^{(0)}(x)  \right )
\end{align}

which has derivative,

\begin{equation}
\frac{d}{dt} \int p (x)\, \mathrm{d} \pi(t)(x)  = \alpha \, \left (  \frac{1}{1-\alpha} \, \int_{\mathbb{R}^D } p(x) \, \mathrm{d} \mu^{(1)}(x) -  \frac{1}{\alpha} \, \int_{\mathbb{R}^D } p(x) \, \mathrm{d} \mu^{(0)}(x) \right ) < 0
\end{equation}

since $0 < \alpha < 1 $ and by assumption (Equation (<ref>)). This proves the statement in the second case.

Under the hypotheses of Proposition <ref>, an element of $\mathcal{S}_2$ of the form,

\begin{equation}
\mu^{(0)}(\cdot) +   \mu^{(1)}(\cdot)  :=  \prod_{i=1}^D \mu^{(0)}_i(\cdot) +  \prod_{i=1}^D \mu^{(1)}_i(\cdot)
\end{equation}

is a global minimum of,

\begin{equation}
\min_{ \mu(\cdot) \in \mathcal{S}_2 }\, \int p(x) \, \mathrm{d} \left (    \prod_{i=1}^D \mu_i^{(0)}(\cdot) +  \prod_{i=1}^D \mu^{(1)}_i(\cdot)   \right )
\end{equation}

if and only if either of $\frac{1}{\alpha^{(0)}} \, \prod_{i=1}^D \mu^{(0)}_i(\cdot)$ or $\frac{1}{\alpha^{(1)}} \, \prod_{i=1}^D \mu^{(1)}_i(\cdot)$ is itself a local minimum, where $\alpha^{(j)} := \prod_{i=1}^D \mu^{(j)}_i([-1,1])$.

First, note that if either of $\frac{1}{\alpha^{(0)}} \, \mu^{(0)} (\cdot) $ or $\frac{1}{\alpha^{(1)}} \, \mu^{(1)} (\cdot) $ is feasible and is a local minimum, then we can invoke Case 1) of Proposition <ref> (i.e., with $\alpha^{(0)} \in \{0,1\}$) with the combination $\frac{1}{\alpha^{(0)}} \, \mu^{(0)} (\cdot) + 0  \in \mathcal{S}_2$ or $0 + \frac{1}{\alpha^{(1)}} \, \mu^{(1)} (\cdot) \in \mathcal{S}_2 $ respectively. This implies that there exists a descent path from each point or that the measure is a global minimum. Since it is a local minimum by assumption, there cannot be a feasible descent path originating from that point and it is therefore a global minimum.

Next, assume the measure is a global minimum with global minimum value $p^*$. If $\mu^{(0)} \left( [-1,1]^D \right ) \in \{0,1\}$, the result follows trivially as the point must also be a local minimum. Otherwise, $\mu^{(0)} \left( [-1,1]^D \right ) \in (0,1)$ and Proposition <ref> implies that,

\begin{equation}
\label{eq:measure_equal_val}
\frac{1}{\alpha^{(0)}} \, \int p(x) \, \mathrm{d} \prod_{i=1}^D \mu^{(0)}_i(x) =  \frac{1}{\alpha^{(1)}} \,\int p(x) \, \mathrm{d} \prod_{i=1}^D \mu^{(1)}_i(x)
\end{equation}

for without an equality, this corresponds to Case 2) of Proposition <ref> and there must exists a descent path, which is a contradiction to local optimality. Equation (<ref>) then implies that,

\begin{align}
p^* &= \int p(x) \, \mathrm{d} \left (  \prod_{i=1}^D \mu^{(0)}_i(x)  +   \prod_{i=1}^D \mu^{(1)}_i(x)\right )  \, \mathrm{d} x\\
&= \frac{1}{\alpha^{(0)}} \, \int p(x) \, \mathrm{d} \prod_{i=1}^D \mu^{(0)}_i(x)
\end{align}

and similarly for $\prod_{i=1}^D \mu^{(1)}_i(x) $. Thus, we find that both normalized product measures achieve the optimal value and must be themselves global (and therefore local) minima, which proves the claim.

Consider a global minimum,

\begin{equation}
\left \{ \sum_{l=1}^2 \prod_{i=1}^D \mu^{(l)}_{i, n_i} \right \}_{\lvert n \rvert \leq d } \in \mathcal{F}_2
\end{equation}

of Problem <ref>.  Let,

$$\mu^{(l)}(\cdot) \in \zeta^{-1} \left ( \left \{ \mu^{(l)}_n \right \}_{\lvert n \rvert \leq d } \right )$$

for $l=1,2$ be a regular Borel product measure supported over $[-1,1]^D$ and associated with the respective moments, which existence is guaranteed by Proposition <ref>. Then, invoking Corollary <ref> (Case 2), we find that,

\begin{equation}
\mu^{(0)}(\cdot) + \mu^{(1)} (\cdot)
\end{equation}

must be a global minimum of Problem <ref>. Being such, Corollary <ref> implies that either,

\begin{equation}
\frac{\mu_i^{(0)}(
\cdot)}{\alpha^{(0)}}    \; \mathrm{, \,or}     \; \frac{\mu_i^{(1)}(
\cdot)}{\alpha^{(1)}}
\end{equation}

is a global minimum of Problem <ref>.
Applying the forward map $\zeta(\cdot)$ of Proposition <ref> and invoking Case 2 of Problem <ref> again implies that either,

\begin{equation}
\left \{ \frac{ \prod_{i=1}^D\mu_{i,n_i}^{(0)} }{\alpha^{(0)}}  \right \}_{\lvert n \rvert \leq d}  \; \mathrm{, \,or}     \; \left \{ \frac{ \prod_{i=1}^D\mu_{i,n_i}^{(1)} }{\alpha^{(1)}}  \right \}_{\lvert n \rvert \leq d}
\end{equation}

are global minima or Problem <ref> as claimed. This proves the first direction.

For the converse, assume w.l.o.g. that $    \left \{ \frac{ \prod_{i=1}^D\mu_{i,n_i}^{(0)} }{\alpha^{(0)}}  \right \}_{\lvert n \rvert \leq d}  $ is a local minimum. In particular, note that Proposition <ref> guarantees the existence of a measure,

$$ \frac{\mu^{(0)}(\cdot) }{\alpha^{(0)}}  := \frac{\prod_{i=1}^D \mu_i^{(0)}(\cdot)}{\alpha^{(0)}}  \in \mathcal{S}_2 $$
associated with these moments. Further, note that Corollary <ref> (Case 3) implies that $ \frac{\mu^{(0)}(\cdot) }{\alpha^{(0)}} $ must be a
local minima. In turn, Corollary <ref> implies that it must actually be a global minimum. Finally Corollary <ref> (Case 2) implies that the moments,

\begin{equation}
\left \{ \left [ \zeta \left ( \frac{\mu^{(0)}(\cdot) }{\alpha^{(0)}} \right ) \right ]_n \right \}_{\lvert n \rvert \leq  d} = \left \{  \frac{ \prod_{i=1}^D \mu_{i,n_i}^{(0)} }{\alpha^{(0)}}  \, \right \}_{\lvert n \rvert \leq  d}
\end{equation}

must be associated with a global minimum of Problem <ref>. This proves the desired result.

[1] Aardal, K. and  G.L. Nemhauser and R. Weismantel.
“Semidefinite Programming and Integer Programming.” Handbooks in Operations Research and Management Science 12 (2005): 393-514.
[2] Andersen, Martin, et al. “Interior-point methods for large-scale cone programming.” Optimization for machine learning 5583 (2011).
[3] Ben-Tal, Aharon, and Arkadi Nemirovski. Lectures on modern convex optimization: analysis, algorithms, and engineering applications. Society for industrial and applied mathematics, 2001.
[4] Boumal, Nicolas, Vlad Voroninski, and Afonso Bandeira. “The non-convex Burer-Monteiro approach works on smooth semidefinite programs.” Advances in Neural Information Processing Systems 29 (2016).
[5] Boyd, Stephen, et al. “Distributed optimization and statistical learning via the alternating direction method of multipliers.” Foundations and Trends in Machine learning 3.1 (2011): 1-122.
[6] Bubeck, Sebastien. “Convex Optimization: Algorithms and Complexity.” arXiv:1405.4980v2. (2014)
[7] Burer, Samuel, and Renato DC Monteiro. “A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization.” Mathematical Programming 95.2 (2003): 329-357.
[8] Cohn, Donald L. Measure theory. Vol. 1. New York: Birkhäuser, 2013.
[9]
de Klerk, Etienne and Monique Laurent. “On the Lasserre Hierarchy of Semidefinite Programming Relaxations of Convex Polynomial Optimization Problems.” SIAM Journal on Optimization 21 (2011):824-832.
[10] de Klerk, Etienne, and Monique Laurent. “Convergence analysis of a Lasserre hierarchy of upper bounds for polynomial minimization on the sphere.” Mathematical Programming 193.2 (2022): 665-685.
[11]
Dressler, Mareike, Sadik Iliman, and Timo de Wolff. “A Positivstellensatz for Sums of Nonnegative Circuit Polynomials.” SIAM Journal of Applied Algebraic Geometry vol 1 (2017):536-555.
[12]
Fujisawa, K., Fukuda, M., Kojima, M., Nakata, K. “Numerical Evaluation of SDPA (Semidefinite Programming Algorithm).” In: Frenk, H., Roos, K., Terlaky, T., Zhang, S. (eds) High Performance Optimization. Applied Optimization vol 33. Springer, Boston, MA.
[13] Golub, Gene H., and Charles F. Van Loan. Matrix computations. JHU press, 2013.
[14] Grant, Michael, Stephen Boyd, and Yinyu Ye. “CVX: Matlab software for disciplined convex programming.” (2011).
[15] Kannan, Hariprasad and Nikos Komodakis and Nikos Paragios. “Chapter 9 - Tighter continuous relaxations for MAP inference in discrete MRFs: A survey.” Handbook of Numerical Analysis: Processing, Analyzing and Learning of Images, Shapes, and Forms: Part 2 (editors: Ron Kimmel and Xue-Cheng Tai)
Elsevier, Vol. 20 (2019): 351-400.
[16] Lasserre, Jean B. “Global optimization with polynomials and the problem of moments.” SIAM Journal on optimization 11.3 (2001): 796-817.
[17]
Lasserre, Jean B. “Sum of Squares Approximation of Polynomials, Nonnegative on a Real Algebraic Set.” SIAM Journal on Optimization 16.2 (2005): 610-628.
[18] Lasserre, Jean B. “Moments and sums of squares for polynomial optimization and related problems.” Journal of Global Optimization 45.1 (2009): 39-61.
[19] Liu, Dong C., and Jorge Nocedal. “On the limited memory BFGS method for large scale optimization.” Mathematical programming 45.1-3 (1989): 503-528.
[20] Mason, John C., and David C. Handscomb. Chebyshev polynomials. CRC press, 2002.
[21] Mittelmann, H. “An independent benchmarking of SDP and SOCP solvers.” Mathematical Programming Ser. B 95, (2003):4017-430.
[22]
Nesterov, Yurii E. and Arkadii Nemirovskii. “Interior-point polynomial algorithms in convex programming.” SIAM studies in applied mathematics 13 (1994).
[23] Nie, Jiawang, and Markus Schweighofer. “On the complexity of Putinar's Positivstellensatz.” Journal of Complexity 23.1 (2007): 135-150.
[24] Parrilo, Pablo A. “Semidefinite programming relaxations for semialgebraic problems.” Mathematical Programming 96:2 (2003):293-320.
[25] Powers, Victoria, and Bruce Reznick. "Polynomials that are positive on an interval." Transactions of the American Mathematical Society 352.10 (2000): 4677-4692.
[26] Prajna, Stephen and Papachristodoulou, Antonis and Parrilo, Pablo A. “Introducing SOSTOOLS: A general purpose sum of squares programming solver.”, Proceedings of the 41st IEEE Conference on Decision and Control 1 (2002): 741-746
[27] Putinar, Mihai. “Positive polynomials on compact semi-algebraic sets.” Indiana University Mathematics Journal 42.3 (1993): 969-984.
[28] Reed, Michael, and Barry Simon. I: Functional analysis. Vol. 1. Gulf Professional Publishing, 1980.
[29] Theodore, S. “Motzkin. The arithmetic-geometric inequality.” Inequalities (Proc. Sympos. Wright-Patterson Air Force Base, Ohio, 1965): 205-224.
[30] Rudin, Walter. Real and complex analysis: 3rd. Ed.. McGraw-Hill, 1987.