Connecting optimization with spectral analysis of tri-diagonal matrices

TL;DR

This paper establishes a novel connection between optimization and spectral analysis by approximating function extrema through eigenvalues of tri-diagonal matrices, linking optimization, orthogonal polynomials, and linear algebra.

Contribution

It introduces a method to approximate global extrema of functions using eigenvalues of tri-diagonal matrices, bridging multiple mathematical fields.

Findings

Extends the spectral analysis approach to optimization problems.

Provides a hierarchy of matrices for increasingly accurate approximations.

Links polynomial roots to optimization extrema.

Abstract

We show that the global minimum (resp. maximum) of a continuous function on a compact set can be approximated from above (resp. from below) by computing the smallest (rest. largest) eigenvalue of a hierarchy of (r x r) tri-diagonal univariate moment matrices of increasing size. Equivalently it reduces to computing the smallest (resp. largest) root of a certain univariate degree-r orthonormal polynomial. This provides a strong connection between the fields of optimization, orthogonal polynomials, numerical analysis and linear algebra, via asymptotic spectral analysis of tri-diagonal symmetric matrices.

Tables

Table 1. Table 1. Comparing τ r ℓ subscript superscript 𝜏 ℓ 𝑟 \tau^{\ell}_{r} and θ ¯ r subscript ¯ 𝜃 𝑟 \underline{\theta}_{r} on 4 toy examples

pb	Moztkin	Matyas	Booth	Three-hump camel
${\underset{¯}{\theta }}_5$	0.801	3.69	69.81	9.58
${\tau }_5^{\ell }$	0.873	2.06	56.64	15.07
${\underset{¯}{\theta }}_6$	0.801	2.99	63.54	4.439
${\tau }_6^{\ell }$	0.808	1.68	45.49	12.68

Equations

$$\underline{f}=\underset{\mathbf{x}}{\inf }\{f(\mathbf{x}):\mathbf{x}\in \mathbf{\Omega }\};\overline{f}=\underset{\mathbf{x}}{\sup }\{f(\mathbf{x}):\mathbf{x}\in \mathbf{\Omega }\}.$$

$${\underline{\theta }}_r=\underset{\sigma }{\inf }\{\int f\sigma d\lambda :\int \sigma d\lambda =1;\sigma \in \Sigma [\mathbf{x}{]}_r\}$$

$${\overline{\theta }}_r=\underset{\sigma }{\sup }\{\int f\sigma d\lambda :\int \sigma d\lambda =1;\sigma \in \Sigma [\mathbf{x}{]}_r\},$$

$${\tau }_r^{\ell }$$

$${\tau }_r^u$$

$${\tau }_r^{\ell }=\underset{\sigma }{\inf }\{\int z\sigma d\#\lambda :\int \sigma d\#\lambda =1;\sigma \in \Sigma [z{]}_r\},$$

$${\lambda }_{\min }(\mathbf{A},\mathbf{C})=\underset{\tau }{\sup }\{\tau :\mathbf{A}⪰\tau \mathbf{C}\}.$$

$$\mathbf{H}_{2}(\boldsymbol{\phi})\,=\,\left[\begin{array}[]{ccc}\phi_{0}&\phi_{1}&\phi_{2}\\ \phi_{1}&\phi_{2}&\phi_{3}\\ \phi_{2}&\phi_{3}&\phi_{4}\end{array}\right].$$

$${\mathbf{H}}_r(g;\mathbf{\varphi })(\alpha ,\beta )=\sum _{\gamma }{g}_{\gamma }{\varphi }_{\gamma +\alpha +\beta },\forall \alpha ,\beta \in {\mathbb{N}}_{.}^n$$

$$\mathbf{H}_{2}(g\,;\boldsymbol{\phi})\,=\,\left[\begin{array}[]{ccc}\phi_{1}&\phi_{2}&\phi_{3}\\ \phi_{2}&\phi_{3}&\phi_{4}\\ \phi_{3}&\phi_{4}&\phi_{5}\end{array}\right].$$

$$\#\lambda (C)=\lambda (f^{-1}(C)),\forall C\in \mathcal{B}(\mathbb{R}).$$

$$\#{\lambda }_k={\int }_{[0+\infty )}{z}^{k}d\#\lambda (z)={\int }_{\mathbf{\Omega }}f(\mathbf{x}{)}^{k}d\lambda (\mathbf{x}),k=0,1,\dots$$

$$f(\mathbf{x}{)}^k=\sum _{\alpha \in {\mathbb{N}}_{2kd}^n}{f}_{k\alpha }{\mathbf{x}}^{\alpha },k=0,1,\dots ,$$

$$\#{\lambda }_k=\sum _{\alpha \in {\mathbb{N}}_{2kd}^n}{f}_{k\alpha }\prod _{i=1}^n\left({\int }_{-1}^1{x}^{{\alpha }_i}dx\right),k\in \mathbb{N}.$$

$$\#{\lambda }_k={\int }_{\mathbf{\Omega }}g(\mathbf{x}{)}^{k}d\lambda$$

$$\#{\lambda }_k=\frac{1}{2\pi }{\int }_{{\mathbb{R}}^2}((xy-1{)}^2+y^2{)}^k\exp (-x^2/2)\exp (-y^2/2)dxdy,k\in \mathbb{N},$$

$$\underline{f}=\min \{x:x\in \mathrm{supp}(\#\lambda )\};\overline{f}=\max \{x:x\in \mathrm{supp}(\#\lambda )\}.$$

$${\tau }_r^{\ell }$$

$${\tau }_r^u$$

$$\underline{f}\le {\tau }_r^{\ell }\le {\tau }_r^u\le \overline{f}\forall r\in \mathbb{N}.$$

$$\underset{r\to \infty }{\lim }{\tau }_r^{\ell }=\underline{f};\underset{r\to \infty }{\lim }{\tau }_r^u=\overline{f}.$$

$${\tau }_r^{\ell }$$

$${\tau }_r^u$$

$$\underset{a}{\sup }\{a:{\mathbf{H}}_r(x;\#\lambda )⪰a{\mathbf{H}}_r(\#\lambda )\}={\lambda }_{\min }({\mathbf{H}}_r(x;\#\lambda ),{\mathbf{H}}_r(\#\lambda )).$$

$$0<a_z=\lambda (f^{-1}(\{z\}))$$

$$x\mapsto D_{1}(x)={\rm det}\left[\begin{array}[]{cc}\#\lambda_{0}&\#\lambda_{1}\\ 1&x\end{array}\right]\,=\,\#\lambda_{0}\,x-\#\lambda_{1}$$

$$x\mapsto D_{2}(x)={\rm det}\left[\begin{array}[]{ccc}\#\lambda_{0}&\#\lambda_{1}&\#\lambda_{2}\\ \#\lambda_{1}&\#\lambda_{2}&\#\lambda_{3}\\ 1&x&x^{2}\end{array}\right]$$

$$x{T}_j(x)=a_j{T}_{j+1}(x)+b_j{T}_j(x)+a_{j-1}{T}_{j-1}(x),\forall x\in \mathbb{R},j\in \mathbb{N},$$

$$d_{j}={\rm det}\left[\begin{array}[]{ccccc}\#\lambda_{0}&\#\lambda_{1}&\cdots&\cdots&\#\lambda_{j-1}\\ \#\lambda_{1}&\#\lambda_{2}&\cdots&\cdots&\#\lambda_{j}\\ \cdots&\cdots&\cdots&\cdots&\cdots\\ \#\lambda_{j-1}&\#\lambda_{j}&\dots&\cdots&\#\lambda_{2j-2}\end{array}\right],\quad j\in\mathbb{N}.$$

$$J=\left[\begin{array}[]{cccccc}b_{0}&a_{0}&0&\cdots&\cdots&0\\ a_{0}&b_{1}&a_{1}&0&\cdots&0\\ 0&a_{1}&b_{2}&a_{2}&\cdots&0\\ 0&0&\cdots&\cdots&\cdots&0\end{array}\right]$$

$$J_{r}=\left[\begin{array}[]{cccccc}b_{0}&a_{0}&0&\cdots&\cdots&0\\ a_{0}&b_{1}&a_{1}&0&\cdots&0\\ 0&a_{1}&b_{2}&a_{2}&\cdots&0\\ 0&\cdots&\cdots&\cdots&\cdots&0\\ 0&0&\cdots&a_{r-2}&b_{r-1}&a_{r-1}\\ 0&0&\cdots&\cdots&a_{r-1}&b_{r}\end{array}\right]$$

$${\hat{\mathbf{H}}}_r(\#\lambda )(i,j)=\int T_i(x)T_j(x)d\#\lambda ={\delta }_{i=j},\forall i,j=0,1,\dots ,r$$

Full text

This work was supported by the European Research Council (ERC) under the ERC-Advanced Grant for the TAMING project, and also by the AI Interdisciplinary Institute ANITI through the French “Investing for the Future PIA3” program under the Grant agreement n*∘*ANR-19-PI3A-0004.

1. Introduction

The goal of this paper is show that the global minimum (resp. maximum) of a continuous function $f$ on a compact set $\mathbf{\Omega}\subset\mathbb{R}^{n}$ can be approximated as closely as desired from above (resp. from below) by the smallest eigenvalues $(\tau^{\ell}_{r})_{r\in\mathbb{N}}$ (resp. largest eigenvalues $(\tau^{u}_{r})_{r\in\mathbb{N}}$) of a sequence of tri-diagonal univariate “moment” matrices of increasing size $r$. Equivalently, $\tau^{\ell}_{r}$ (resp. $\tau^{u}_{r}$) is the smallest (resp. largest) root of a univariate orthogonal polynomial of increasing degree $r$. Thus it reveals a (perhaps suprising) connection between the fields of optimization and the asymptotic spectral analysis of tri-diagonal symmetric matrices (also related to the asymptotic analysis of univariate orthogonal polynomials).

Notice that there is a large body of literature on numerical analysis of tri-diagonal symmetric matrices for which efficient specialized algorithms exist (for instance the characteristic polynomial can be computed efficiently and roots of univariate polynomials can also be computed efficiently); see e.g. [9, 19].

Our result is valid in a quite general context, namely if $f$ is a continuous function and $\mathbf{\Omega}$ is compact. However to turn it into a practical algorithm requires effective computation of such matrices. This in turn requires computing integrals $\int_{\mathbf{\Omega}}f(\mathbf{x})^{k}d\lambda$ for arbitrary $k$, where $\lambda$ is a measure whose support is exactly $\mathbf{\Omega}$ (and whose choice is left to the user). If $f$ is a polynomial, $\mathbf{\Omega}$ is a “simple” compact set like a box, an ellipsoid, a simplex, a sphere, the discrete hypercube $\{0,1\}^{n}$ (or their image by an affine mapping) and $\lambda$ is a distinguished measure (e.g., Lebesgue measure (or other families of polynomial densities), rotation invariant (for the sphere), counting measure for $\{0,1\}^{n}$) then such integrals are obtained in closed-form.

In addition, after the initial version of this paper was released, Laurent and Slot [15] proved a $O(\log^{2}r/r^{2})$ rate of convergence for the monotone sequence of upper bounds $(\tau^{\ell}_{r})_{r\in\mathbb{N}}$ on the minimum (resp. lower bounds $(\tau^{u}_{r})_{r\in\mathbb{N}}$ on the maximum) under a weak geometric condition on $\mathbf{\Omega}$ (satisfied by the above simple sets).

Background

Let $f:\mathbf{\Omega}\to\mathbb{R}$ be a continuous function on a compact set $\mathbf{\Omega}\subset\mathbb{R}^{n}$, and consider the optimization problem:

[TABLE]

In [14] the author showed that one may approximate $\underline{f}$ from above (resp. $\overline{f}$ from below) by solving the following hierarchy111The hierarchy (1.2) should not be confused with the Moment-SOS hierarchy described in [12] to solve the same problem (1.1), but which yields a converging sequence of lower bounds $\rho_{r}\uparrow\underline{f}$ as $r\to\infty$. of optimization problems indexed by $r\in\mathbb{N}$:

[TABLE]

where $\lambda$ is a fixed measure whose support is exactly $\mathbf{\Omega}$ and $\Sigma[\mathbf{x}]_{r}$ is the convex cone of SOS polynomials of degree at most $2r$. Indeed $\underline{\theta}_{r}\downarrow\underline{f}$ (resp. $\overline{\theta}_{r}\uparrow\overline{f}$) as $t\to\infty$. If $f$ is a polynomial and one knows all moments of the measure $\lambda$ on $\mathbf{\Omega}$, then each problem (1.2) is a very specific semidefinite program. As a matter of fact it reduces to a generalized eigenvalue problem involving two moment-like matrices whose size ${n+r\choose r}$ increases with $r$. For instance this is the case whenever $\mathbf{\Omega}$ is a simple set (e.g. the box $[-1,1]^{n}$, the Euclidean unit ball, the sphere, the simplex, the discrete hypercube, and their affine transformation) and $\lambda$ is an appropriate measure (e.g., Lebesgue measure on unit box, Euclidean unit ball, simplex, rotation invariant measure on the sphere, counting measure on the discrete hypercube, etc.).

In a recent series of papers [3, 4, 5, 6, 20], de Klerk, Laurent and their co-workers have been able to analyze the convergence $\underline{\theta}_{r}\downarrow\underline{f}$ of such a hierarchy by appropriate clever choices of the reference measure (as indeed $\lambda$ it can be any measure whose support is exactly $\mathbf{\Omega}$). Ultimately they could provide rates of convergence. In particular and remarkably, they show that for certain important sets and reference measures (e.g. the box $\mathbf{\Omega}=[-1,1]^{n}$ and the sphere $\mathbb{S}^{n-1}$) a convergence rate $O(1/r^{2})$ is achieved.

However all variants of (1.2) and (1.3) consider density polynomials $\sigma\in\mathbb{R}[\mathbf{x}]$ in $n$ variables which results in eigenvalue problems with multivariate Hankel-type matrices whose size ${n+r\choose n}$ grows very quickly with $r$. Therefore and so far, this hierarchy of upper (resp. lower) bounds has not been proved to be efficient in practice (even for small size problems) and its main interest is rather theoretical as it provides an algorithm with proven rate of convergence $O(1/r^{2})$ to the global optimum on some simple sets.

Contribution

Our contribution is twofold:

I. We first provide an alternative converging hierarchy of upper bounds in the same spirit as (1.2) and (1.3) but following a different path. The main distinguishing feature is to reduce the $n$-dimensional initial problem to a one-dimensional equivalent problem by using the pushforward measure $\#\lambda$ of the measure $\lambda$ on $\mathbf{\Omega}$, by the mapping $f:\mathbf{\Omega}\to\mathbb{R}$. It results in solving again a hierarchy of eigenvalue problems but with the major advantage of considering Hankel moment matrices in just ONE variable (hence of size $r$ in contrast to ${n+r\choose n}$).

To achieve this result we exploit the fact that $\underline{f}$ and $\overline{f}$ are the left and right endpoints of the support of $\#\lambda$ on the real line, and therefore by invoking Lasserre [13, Theorem 3.3], one may approximate $\underline{f}$ from above by solving:

[TABLE]

where the real symmetric $(r+1)\times(r+1)$ matrix $\mathbf{H}_{r}(\#\lambda)$ (resp. $\mathbf{H}_{r}(x;\#\lambda)$) is the moment matrix of the pushforward $\#\lambda$ (resp. the localizing matrix associated with $\#\lambda$ and the univariate linear polynomial $x\mapsto x$). Similarly, one may approximate $\overline{f}$ from below by solving:

[TABLE]

and indeed $\tau^{\ell}_{r}\downarrow\underline{f}$ (resp. $\tau^{u}_{r}\uparrow\overline{f}$) as $r$ increases; see [13].

Remark 1.1**.**

Equivalently, by duality in convex optimization:

[TABLE]

where $\Sigma[z]_{r}$ is the convex cone of univariate SOS polynomials of degree at most $2r$; similarly for $\tau^{u}_{r}$ just replace “$\inf$” by “$\sup$”. The formulation (1.5) resembles (1.2) but with the important difference that in (1.2) one searches over SOS polynomials of degree at most $r$ in $n$-variables whereas in (1.5) one searches over SOS polynomials of degree at most $r$ in ONE variable.

Also notice that in contrast to (1.2), in (1.5) (or in (1.4)) the function $f$ to minimize does not appear explicitly; it is encoded in the pushforward measure $\#\lambda$. So if one is able to compute (or approximate) the moments of $\#\lambda$, then both matrices $\mathbf{H}_{r}(x;\,\#\lambda)$ and $\mathbf{H}_{r}(\#\lambda)$ are known and (1.4) can be solved in practice. Also to solve (1.2) in practice one needs $f$ to be a polynomial.

II. We next further simplify the analysis by considering an orthonormal basis $(T_{j})_{j\in\mathbb{N}}$ of polynomials w.r.t. the pushforward measure $\#\lambda$. Recall that a family of $T_{j}$’s, $j\leq r$, can be obtained from the moment matrix $\mathbf{H}_{r}(\#\lambda)$ by simple determinant computations. In this new basis $(T_{j})_{j\in\mathbb{N}}$, the moment matrix $\hat{\mathbf{H}}_{r}(\#\lambda)$ becomes the identity and the localizing matrix $\hat{\mathbf{H}}_{r}(x;\,\#\lambda)$ now becomes a tri-diagonal (banded) symmetric matrix whose coefficients have a direct expression in terms of the parameters defining the classical three-term recurrence relation satisfied by the $T_{j}$’s.

Therefore the convergence $\tau^{\ell}_{r}\downarrow\underline{f}$ (resp. $\tau^{u}_{r}\uparrow\overline{f}$) is simply the asymptotic behavior of the smallest (resp. largest) eigenvalue of tri-diagonal $r\times r$ “moment” matrices, as $r$ increases; equivalently the asymptotic behavior of the smallest (rest. largest) root of a certain univariate polynomial orthogonal w.r.t. $\#\lambda$.

This reveals a strong and perhaps surprising connection between the fields of (global) optimization and the spectral analysis of tri-diagonal Hankel matrices (or extremal roots of a family of orthogonal polynomials). Actually such a link already appeared in de Klerk and Laurent to analyze convergence of upper bounds (1.2) for the specific univariate (trivial) optimization problem $\min\,\{x:x\in[-1,1]\}$ and specific reference measure $\lambda$. Then they used this univariate problem as a building block to prove rates of convergence of the bounds in (1.2) in case where $\mathbf{\Omega}=[-1,1]^{n}$ or $\mathbb{S}^{n-1}$ (and with specific measures $\lambda$).

There is a large body of literature on various aspects of tri-diagonal symmetric matrices, including practical efficient algorithms; see for instance Aurentz et al. [1], Businger [2], Ford [9], Kiliç [11], Mallik [16], Osipov [18], and Routh [19]. Therefore this may also open the door to practical algorithms for good approximations of $\underline{f}$ and $\overline{f}$ in non trivial cases, as soon as one can obtain moments of the measure $\#\lambda$ for reasonably large degrees.

If the sequence $(\tau^{\ell}_{r})_{r\in\mathbb{N}}$ (or $(\tau^{u}_{r})_{r\in\mathbb{N}}$) has obvious numerical advantages when compared with $(\underline{\theta}_{r})_{r\in\mathbb{N}}$ (or $(\overline{\theta}_{r})_{r\in\mathbb{N}}$) for a same fixed $r$, the main drawback of $(\tau^{\ell}_{r})_{r\in\mathbb{N}}$ is the computation of moments of $\#\lambda$ which may become tedious for non modest dimension $n$. However, sparsity of $f$ (i.e. when $f$ has a few monomials only) can be exploited. Also the rate of convergence for $\tau^{\ell}_{r}\downarrow\underline{f}$ as $r\to\infty$, is more difficult to analyze because we do not know the density of the pushforward measure $\#\lambda$ with respect to Lebesgue measure on $\mathbf{\Omega}$.

Actually, recently Laurent and Slot [15] have provided a further detailed analysis of some relative merits of the sequences $(\tau^{\ell}_{r})_{r\in\mathbb{N}}$ and $(\underline{\theta}_{r})_{r\in\mathbb{N}}$ beyond the scope the present paper. As already mentioned, they prove that remarkably, $\tau^{\ell}_{r}\downarrow\underline{f}$ at a $O(\log^{2}r/r^{2})$ rate under a weak geometric condition on $\mathbf{\Omega}$ (satisfied for the simple sets $\mathbf{\Omega}$ mentioned previously).

2. Main result

2.1. Notation, definitions and preliminary results

Let $\mathbb{R}[\mathbf{x}]$ denote the ring of polynomials in the $n$ variables $\mathbf{x}=(x_{1},\ldots,x_{n})$ and $\mathbb{R}[\mathbf{x}]_{t}\subset\mathbb{R}[\mathbf{x}]$ denote the vector space of polynomials of degree at most $t$, hence of dimension $s(t)={n+t\choose n}$. Let $\Sigma[\mathbf{x}]\subset\mathbb{R}[\mathbf{x}]$ denote the space of polynomials that are sums-of-squares (in short SOS polynomials) and let $\Sigma[\mathbf{x}]_{t}\subset\mathbb{R}[\mathbf{x}]_{2t}$ denote the space of SOS polynomials of degree at most $2t$. For univariate polynomials in the variable $x$, we use the notation $\mathbb{R}[x]$, $\Sigma[x]$, $\mathbb{R}[x]_{t}$ and $\Sigma[x]_{t}$.

With $\alpha\in\mathbb{N}^{n}$ and $\mathbf{x}\in\mathbb{R}^{n}$, the notation $\mathbf{x}^{\alpha}$ stands for $x_{1}^{\alpha_{1}}\cdots x_{n}^{\alpha_{n}}$. Also for every $\alpha\in\mathbb{N}^{n}$, let $|\alpha|:=\sum_{i}\alpha_{i}$ and $\mathbb{N}^{n}_{t}:=\{\alpha\in\mathbb{N}^{n}:|\alpha|\leq t\}$, where $\mathbb{N}=\{0,1,2,\ldots\}$.

The support of a Borel measure $\mu$ on $\mathbb{R}^{n}$ is the smallest closed set $\mathbf{\Omega}$ such that $\mu(\mathbb{R}^{n}\setminus\mathbf{\Omega})=0$. Denote by $\mathcal{B}(\mathbf{X})$ the Borel $\sigma$-field associated with a topological space $\mathbf{X}$, and $\mathscr{M}(\mathbf{X})_{+}$ the space of finite nonnegative Borel measures on $\mathbf{X}$.

Generalized eigenvalue

Given two real symmetric matrices $\mathbf{A},\mathbf{C}\in\mathbb{R}^{n\times n}$ denote by $\lambda_{\min}(\mathbf{A},\mathbf{C})$ the smallest generalized eigenvalue with respect to the pair $(\mathbf{A},\mathbf{C})$, that is, the smallest scalar $\tau$ such that $\mathbf{A}\mathbf{x}=\tau\,\mathbf{C}\mathbf{x}$ for some nonzero vector $\mathbf{x}\in\mathbb{R}^{n}$. When $\mathbf{C}$ is the identity matrix then $\lambda_{\min}(\mathbf{A},\mathbf{C})$ is just the smallest eigenvalue of $\mathbf{A}$. Computing $\lambda_{\min}(\mathbf{A},\mathbf{C})$ can be done via a pure and efficient linear algebra routine. The notation $\mathbf{A}\succeq 0$ (resp. $\mathbf{A}\succ 0$) stands for $\mathbf{A}$ is positive semidefinite (resp. positive definite). If $\mathbf{A},\mathbf{C}\succ 0$ then:

[TABLE]

Moment matrix

Given a real sequence $\boldsymbol{\phi}=(\phi_{\alpha})_{\alpha\in\mathbb{N}^{n}}$, let $\mathbf{H}_{r}(\boldsymbol{\phi})$ denote the multivariate (Hankel-type) moment matrix defined by $\mathbf{H}_{r}(\boldsymbol{\phi})(\alpha,\beta)=\phi_{\alpha+\beta}$ for all $\alpha,\beta\in\mathbb{N}^{n}_{r}$. For instance, in the univariate case $n=1$, with $r=2$, $\mathbf{H}_{2}(\boldsymbol{\phi})$ is the Hankel matrix

[TABLE]

If $\boldsymbol{\phi}=(\phi_{j})_{j\in\mathbb{N}}$ is the moment sequence of a Borel measure $\phi$ on $\mathbb{R}$ then $\mathbf{H}_{r}(\boldsymbol{\phi})\succeq 0$ for all $r=0,1,\ldots$. Conversely, if $\mathbf{H}_{r}(\boldsymbol{\phi})\succeq 0$ for all $r\in\mathbb{N}$, then $\boldsymbol{\phi}$ is the moment sequence of some finite (nonnegative) Borel measure $\phi$ on $\mathbb{R}$. The converse result is not true anymore in the multivariate case.

Localizing matrix

Given a real sequence $\boldsymbol{\phi}=(\phi_{\alpha})_{\alpha\in\mathbb{N}^{n}}$ and $g\in\mathbb{R}[\mathbf{x}]$, $\mathbf{x}\mapsto g(\mathbf{x})=\sum_{\alpha}g_{\alpha}\,\mathbf{x}^{\alpha}$, $\mathbf{H}_{r}(g\,;\boldsymbol{\phi})$ denote the multivariate (Hankel-type) localizing matrix defined by

[TABLE]

For instance, in the univariate case $n=1$, with $r=2$ and $g(x)=x$, $\mathbf{H}_{2}(g\,;\boldsymbol{\phi})$ is the Hankel matrix

[TABLE]

Equivalently, if $\boldsymbol{\phi}$ is the moment sequence of a Borel measure $\phi$ on $\mathbb{R}^{n}$ then $\mathbf{H}_{r}(g\,;\boldsymbol{\phi})$ is the moment matrix $\mathbf{H}_{r}(\boldsymbol{\nu})$ associated with the measure $d\nu=g\,d\phi$.

Pushforward measure

Let $\lambda\in\mathscr{M}(\mathbf{\Omega})_{+}$ be a finite Borel measure on $\mathbf{\Omega}\subset\mathbb{R}^{n}$ whose support is exactly $\mathbf{\Omega}$, that is, $\mathbf{\Omega}$ is the smallest closed set such that $\lambda(\mathbb{R}^{n}\setminus\mathbf{\Omega})=0$. Let $\#\lambda$ be the pushforward measure on $\mathbb{R}$ of $\lambda$ with respect to (w.r.t.) the mapping $f:\mathbf{\Omega}\to\mathbb{R}$. That is:

[TABLE]

In particular, its moments $\boldsymbol{\#\lambda}=(\#\lambda_{k})_{k\in\mathbb{N}}$ read:

[TABLE]

It is straightforward to see that the support of $\#\lambda$ is contained in the interval $[\underline{f},\overline{f}]$ with $\underline{f}$ and $\overline{f}$ as its left and right endpoints.

Remark 2.1** (Computing the moments $\boldsymbol{\#\lambda}$).**

While (2.2) is quite general, exact numerical computation of $\#\lambda_{k}$ is not possible in full generality. However, for specific combinations of ($f$, $\mathbf{\Omega}$, $\lambda$) one may obtain $\boldsymbol{\#\lambda}$ in closed-form (let alone combinations where numerical approximations schemes, e.g. cubature formula or Monte-Carlo, can be used). This is the case when $f$ is a polynomial and $\mathbf{\Omega}$ is a “simple” set like a box, an ellipsoid, a simplex, a sphere, the discrete hypercube $\{0,1\}^{n}$, or their image by an affine map. Then several choices of $\lambda$ are possible: Lebesgue measure (or specific families of polynomial densities) if $\mathbf{\Omega}$ is a box, simplex or ellipsoid, the rotation invariant measure if $\mathbf{\Omega}$ is the sphere, the discrete counting measure if $\mathbf{\Omega}$ is the hyper cube, etc. For instance, if $\mathbf{\Omega}=[-1,1]^{n}$ and $\lambda$ is Lebesgue measure then writing

[TABLE]

one obtains

[TABLE]

Similarly if $\mathbf{\Omega}=\{\mathbf{x}:\|\mathbf{x}\|_{2}^{2}\leq 1\,\}$ then for every $k\in\mathbb{N}$:

[TABLE]

which is also obtained in closed form. Notice that for large $k$, computing $\#\lambda_{k}$ can be time and space consuming.

Remark 2.2** (On compactness assumption).**

In fact non compact sets, e.g., $\mathbf{\Omega}=\mathbb{R}^{n}$ with $\lambda$ being the normal distribution $\mathcal{N}(0,I)$, or $\mathbf{\Omega}=\mathbb{R}^{n}_{+}$ with $\lambda$ being the exponential distribution $d\lambda=\prod_{i}\exp(-x_{i})\,d\mathbf{x}$, can be treated as well; see [14].

In particular, and interestingly, it permits to provide upper bounds that converge to the global minimum $f^{*}$ of $f$ on $\mathbb{R}^{n}$, even if $f^{*}$ is not attained (a difficult case in practice because then $\|\mathbf{x}\|\to\infty$ if $f(\mathbf{x})\to f^{*}$). For instance with $f(x,y)=(xy-1)^{2}+y^{2}$ on $\mathbb{R}^{2}$, $f^{*}=0$ is not attained and $f(n,1/n)\to 0$ as $n\to\infty$. Then letting $d\lambda=(2\pi)^{-1}\exp(-(x^{2}-y^{2})/2)dxdy$, the pushforward $\#\lambda$ is supported on $[0,+\infty)$, and

[TABLE]

is obtained in closed-form. However, for ease of exposition and to avoid technicalities, we restrict to the compact case.**

2.2. Main result

We have just seen that the support of the pushforward $\#\lambda$ is precisely contained in the interval $[\underline{f},\overline{f}]$, with $\underline{f},\overline{f}$ as its left and right endpoints. Therefore the problem of approximating $\underline{f}$ and $\overline{f}$ reduces to approximate the endpoints of the support of $\#\lambda$ from the sole knowledge of its moments. That is,

[TABLE]

In [13] the author has already considered the more general problem of bounding the support of a measure $\mu$ on $\mathbb{R}^{n}$ from knowledge of its marginal moments. In our case $\mu$ is the push forward measure $\#\lambda$ on $\mathbb{R}$ and therefore we can invoke Theorem 3.3 in [13]. More precisely:

Let $\mathbf{\Omega}\subset\mathbb{R}^{n}$ be compact with nonempty interior, $\lambda\in\mathscr{M}(\mathbf{\Omega})_{+}$ and consider the hierarchy of optimization problems indexed by $r\in\mathbb{N}$:

[TABLE]

where $\mathbf{H}_{r}(x;\,\#\lambda)$ is the (univariate) localizing matrix associated with the polynomial $x\mapsto x$ and the measure $\#\lambda$ on $\mathbb{R}$, and $\mathbf{H}_{r}(\#\lambda)$ is the (univariate) moment matrix associated with $\#\lambda$.

Theorem 2.3**.**

Let $\lambda\in\mathscr{M}(\mathbf{\Omega})_{+}$ be such that ${\rm supp}(\lambda)=\mathbf{\Omega}$. Then:

[TABLE]

In addition the sequence $(\tau_{r}^{\ell})_{r\in\mathbb{N}}$ (resp. $(\tau_{r}^{u})_{r\in\mathbb{N}}$) is monotone non-increasing (resp. non-decreasing), and:

[TABLE]

Finally, for all $r\in\mathbb{N}$:

[TABLE]

Proof.

That the sequence $(\tau^{\ell}_{r})_{r\in\mathbb{N}}$ is monotone non-increasing is straightforward as the feasible set in (2.4) shrinks with $r$. The same argument shows that $(\tau^{u}_{r})_{r\in\mathbb{N}}$ is monotone non-decreasing. Next, the support of $\#\lambda$ is contained in the interval $[\underline{f},\overline{f}]$ with $\underline{f}$ and $\overline{f}$ as its left and right endpoints, and we know all moments of $\#\lambda$. Therefore (2.6)-(2.8) follow from [13, Theorem 3.3, p. 3379]. ∎

Even though $\mathbf{H}_{r}(x;\,\#\lambda)$ may not be positive definite we still have

[TABLE]

This is because $\mathbf{H}_{r}(x-a;\,\#\lambda)=\mathbf{H}_{r}(x;\,\#\lambda)-a\,\mathbf{H}_{r}(\#\lambda)$; see Lemma 3.1.

So after one has reduced the $n$-dimensional problem (1.2) to the one-dimensional problem (1.5), Theorem 2.3 shows that one thus has to handle Hankel moment matrices of size $r$ only whereas $\underline{\theta}_{r}$ requires to handle Hankel-like moment matrices of size ${n+r\choose r}$. However the moment information needed to build up the moment and localizing matrices $\mathbf{H}_{r}(\#\lambda)$, $\mathbf{H}_{r}(x;\,\#\lambda)$, still requires computing $\int_{\mathbf{\Omega}}f(\mathbf{x})^{2k}d\lambda$ with $k\leq 2r$ (hence handling $n$-variate moments up to degree $2d$). But once this has been done, the eigenvalue problem is much easier.

2.3. On finite convergence

For same reasons as for the convergence $\underline{\theta}_{r}\downarrow\underline{f}$, in (2.7) the convergence $\tau^{\ell}_{r}\downarrow\underline{f}$ as $r\to\infty$, is only asymptotic and not finite in general (and similarly for $\tau^{u}_{r}\uparrow\underline{f}$ in (2.9)). Indeed the rationale behind (1.5) is to approximate the Dirac measure $\delta_{z^{*}}$ at the minimizer $z^{*}=\underline{f}\in\mathbb{R}$ by a measure of the form $\sigma d\#\lambda$, and such an approximation cannot be exact because the pushforward measure $\#\lambda$ has no atom if $\mathbf{\Omega}$ is not finite.

An exception is when $\mathbf{\Omega}$ is a finite set (e.g., the discrete hypercube $\{0,1\}^{n}$ or $\{-1,1\}^{n}$) and $\lambda$ is the counting measure $\sum_{\mathbf{v}\in\mathbf{\Omega}}\delta_{\mathbf{v}}$, with $\delta_{\mathbf{v}}$ being the Dirac measure at $\mathbf{v}\in\mathbf{\Omega}$. Introduce the finite set $\mathcal{Z}:=\{f(\mathbf{v}):\mathbf{v}\in\mathbf{\Omega}\}\subset\mathbb{R}$.

Lemma 2.4**.**

Let $\mathbf{\Omega}\subset\mathbb{R}^{n}$ be a finite set and let $\lambda$ be the counting measure on $\mathbf{\Omega}$. Then finite convergence takes place at most at $r=\bar{r}=\mathrm{card}(\mathcal{Z})-1$. That is, $\tau^{\ell}_{\bar{r}}=\underline{f}$ and $\tau^{u}_{\bar{r}}=\overline{f}$.

Proof.

The support of the pushforward measure $\#\lambda$ is the finite set $\mathcal{Z}$ with smallest element $z^{*}:=\underline{f}\in\mathcal{Z}$ and largest element $\overline{f}\in\mathcal{Z}$. In addition, the pushforward measure $\#\lambda$ satisfies $\#\lambda=\sum_{z\in\mathcal{Z}}a_{z}\,\delta_{z}$, where:

[TABLE]

Let $p\in\mathbb{R}[z]_{\bar{r}}$ be the degree-$\bar{r}$ polynomial $z\mapsto\prod_{a\in\mathcal{Z};a\neq z^{*}}\frac{z-a}{z^{*}-a}$, so that $p(z)=0$ for all $z^{*}\neq z\in\mathcal{Z}$, and $p(z^{*})=1$. Then letting $\sigma^{*}:=a_{z^{*}}^{-1}p^{2}\in\Sigma[z]_{\bar{r}}$, one obtains $\int\sigma^{*}d\#\lambda=\sum_{z\in\mathcal{Z}}a_{z}\,\sigma^{*}(z)=p(z^{*})^{2}=1$, and $\int z\sigma^{*}\,d\#\lambda\,=\,z^{*}=\underline{f}$. This proves that $\sigma^{*}$ is an optimal solution of (1.5) when $r=\bar{r}$ (recall that $\tau^{\ell}_{r}\geq\underline{f}$ for all $r$) and so $\tau^{\ell}_{\bar{r}}=\underline{f}$ (and with a similar argument, $\tau^{u}_{\bar{r}}=\overline{f}$). ∎

So Lemma 2.4 states that $\bar{r}=\mathrm{card}(\mathcal{Z})-1$ determines when finite convergence is guaranteed, i.e., $\tau^{\ell}_{\bar{r}}=\underline{f}$. In other words, the larger is the set $\mathcal{Z}$ the more difficult is to compute $\underline{f}$ (or $\overline{f}$).

2.4. Link with tri-diagonal matrices

Let $(T_{j})_{j\in\mathbb{N}}$ be a family of orthonormal polynomials with respect to $\#\lambda$. For instance, such a family can be computed from the moments $(\#\lambda_{k})_{k\in\mathbb{N}}$ as follows. $T_{0}(x)=1=D_{0}(x)/\#\lambda_{0}$ for all $x\in\mathbb{R}$. Then compute the degree-one polynomial:

[TABLE]

and normalize $T_{1}(x)=a\,D_{1}(x)$ to obtain $\int T_{1}(x)^{2}d\#\lambda=1$, i.e., $T_{1}(x)=a\,D_{1}(x)$ with $a=\#\lambda_{0}^{-1/2}(\#\lambda_{0}\#\lambda_{2}-\#\lambda_{1}^{2})^{-1/2}$. Next, to obtain $T_{2}$ compute

[TABLE]

and again normalize $T_{2}(x)=a\,D_{2}(x)$ to obtain $\int T_{2}^{2}d\#\lambda=1$, etc. Next, the orthonormal polynomials satisfy the so-called three-term recurrence relation:

[TABLE]

where $a_{j}=(d_{j}\,d_{j+2}/d_{j+1}^{2})^{1/2}$, $b_{j}=\int xT_{j}(x)^{2}\,d\#\lambda$, and

[TABLE]

The tri-diagonal infinite matrix:

[TABLE]

is called the Jacobi matrix associated with the orthonormal polynomials $(T_{j})_{j\in\mathbb{N}}$; see e.g. Dunkl and Xu [8, pp. 10–11].

Proposition 2.5**.**

Expressed in the basis of orthonormal polynomials $(T_{j})_{j\in\mathbb{N}}$, the moment matrix $\hat{\mathbf{H}}_{r}(\#\lambda)$ is the identity matrix while the $(r+1)\times(r+1)$ localizing matrix $\hat{\mathbf{H}}_{r}(x;\,\#\lambda)$ is the $r$-truncation

[TABLE]

of the Jacobi matrix (2.11).

Proof.

That the moment matrix $\hat{\mathbf{H}}_{r}(\#\lambda)$ expressed in the basis $(T_{j})_{j\in\mathbb{N}}$ is the identity matrix follows from

[TABLE]

Next, in this basis the localizing matrix $\hat{\mathbf{H}}_{r}(x;\,\#\lambda)$ associated with $\#\lambda$ and the polynomial $x\mapsto x$, reads:

[TABLE]

for all $i,j=0,1,\ldots,r$, where we have used (2.10). Hence $\hat{\mathbf{H}}_{r}(x;\,\#\lambda)$ is a tri-diagonal matrix where at row $i$ the three elements are $(a_{i-1},b_{i},a_{i+1})$. Therefore, $\hat{\mathbf{H}}_{r}(x;\#\lambda)$ is the $r$-truncation of the Jacobi matrix (2.11). ∎

As a consequence we thus obtain:

Corollary 2.6**.**

Let $\tau^{\ell}_{r}$ and $\tau^{u}_{r}$ be as in Theorem 2.3 and let $J_{r}$ be the tri-diagonal matrix in Proposition 2.5. Then $\tau^{\ell}_{r}=\lambda_{\min}(J_{r})$ and $\tau^{u}_{r}=\lambda_{\max}(J_{r})$. Therefore:

[TABLE]

Also for every $r\in\mathbb{N}$, $\tau^{\ell}_{r}$ (resp. $\tau^{u}_{f}$) is the smallest (resp. largest) root of the univariate polynomial $T_{r+1}$.

Proof.

Follows from Theorem 2.3 and the definition of $\tau^{\ell}_{r}$ and $\tau^{u}_{r}$. The last statement can be found in Dunkl and Xu [8, Theorem 1.3.12]. ∎

Remark 2.7**.**

The use of an othonormal polynomial basis to reduce the initial $n$-dimensional problem (1.2) to a standard ($n$-dimensional) eigenvalue problem was already proposed in [14] and in de Klerk et al. [5] and de Klerk and Laurent [6] but for the original $n$-dimensional problem and the reference measure $\lambda$ (and not on $\mathbb{R}$ for the pushforward measure $\#\lambda$).

In addition, in de Klerk and Laurent have used the univariate problem $\min\{x:x\in[-1,1]\,\}$ as a building block to prove the $O(1/r^{2})$ rate of convergence for the bounds (1.2) and (1.3) in the multivariate case of the Sphere, the unit box and for some different reference measures $\lambda$. They observed that if $f$ is the univariate polynomial $x$ then solving the resulting eigenvalue problem is computing the smallest eigenvalue of the Jacobi matrix associated with $\lambda$ (or equivalently, the smallest root of a certain orthogonal polynomial) as in Corollary 2.6. For specific reference measures $\lambda$, the associated orthogonal polynomials have been well-studied (e.g. Chebyshev or Legendre polynomials), in particular the asymptotic behavior of their smallest (or largest) root used by the authors in . For more details the interested reader is referred to [3, 4, 5, 6]. However, in Corollary 2.6 the underlying univariate problem $\min\{\,x:x\in{\rm supp}(\#\lambda)\,\}$ in (2.3) is equivalent to the original multivariate problem (1.1). The price to pay is that the density of $\#\lambda$ is not known explicitly and makes the analysis of the rate of convergence more intricate.

Discussion

As already mentioned, computing the entries of $\mathbf{H}_{r}(\#\lambda)$ (and hence of $\mathbf{H}(x;\,\#\lambda)$ as well) is easy but tedious for large $n$. For a fixed $r\in\mathbb{N}$, and once the moment matrix $\mathbf{H}_{r}(\#\lambda)$ has been computed, computing the scalar $\tau^{u}_{r}$ or $\tau^{\ell}_{r}$ in (2.4)-(2.5) is definitely easier than computing $\underline{\theta}_{r}$ or $\overline{\theta}_{r}$ as in the former one handles univariate moment matrices of size $r$ instead of $n$-variate moment matrices of size ${n+d\choose n}$ in the latter. However we have not proved any rate of convergence for $\tau^{\ell}_{r}\downarrow\underline{f}$ whereas $\underline{\theta}_{r}\downarrow\underline{f}$ at a rate $O(1/r^{2})$ for some simple sets $\mathbf{\Omega}$ and appropriate measures $\lambda$. As mentioned earlier, convergence analysis of the sequence $(\tau^{\ell}_{r})_{r\in\mathbb{N}}$ is difficult because we do not have an explicit expression of the density of $\#\lambda$ w.r.t. Lebesgue measure on $[\underline{f},\overline{f}]$.

For illustration purpose, for $r=5,6$, we have considered four toy problems in $n=2$ variables to compare the upper bounds $\tau^{\ell}_{r}$ on $\underline{f}$ obtained in (2.4) with the upper bounds $\underline{\theta}_{r}$ obtained in (1.2) as described in de Klerk and Laurent [5]. Hence for the same $r$, the former are obtained by solving eigenvalue problems with matrices of size $6$ for $r=5$ (resp. size $7$ for $r=6$) as opposed to matrices of size ${2+r\choose 2}=21$ for $r=5$ (resp. size $28$ for $r=6$) for the latter.

Motzkin polynomial: $f(\mathbf{x})=64\,(x_{1}^{4}x_{2}^{2}+x_{1}^{2}x_{2}^{4})-48\,x_{1}^{2}x^{2}_{2}+1$

Matyas function: $f(\mathbf{x})=26\,(x_{1}^{2}+x_{2}^{2})-48\,x_{1}x_{2}$

Booth function: $f(\mathbf{x})=(10\,x_{1}+20\,x_{2}-7)^{2}+(20\,x_{1}+10\,x_{2}-5)^{2}$

Three-hum-camel function: $f(\mathbf{x})=5^{6}x_{1}^{6}/6-1.05*5^{4}x_{1}^{4}+50\,x_{1}^{2}+245\,x_{1}x_{2}+25\,x_{2}^{2}$

In Table 1 are displayed the results. Except for the Motzkin polynomial the bounds $\tau^{\ell}_{r}$ are comparable and even better than the bounds $\underline{\theta}_{r}$.

Convergence rate

After the initial version of this paper was released on a public repository, Laurent and Slot [15] have analyzed the convergence $\tau^{\ell}_{r}\downarrow\underline{f}$ as $r\to\infty$, of Theorem 2.3. Quite remarkably they prove a $O(\log^{2}r/r^{2})$ rate of convergence under a weak geometric condition on $\mathbf{\Omega}$ (satisfied for the “simple” sets $\mathbf{\Omega}$). Their proof is in the same vein and spirit as the one in Slot and Laurent [20] where the authors prove convergence $\underline{\theta}_{r}\downarrow\underline{f}$ at a $O(\log^{2}r/r^{2})$ rate for general convex bodies. For more details the interested reader is referred to [15].

3. Appendix

Lemma 3.1**.**

Let $\mathbf{H}_{r}(x;\,\#\lambda)$ be the Hankel matrix associated with $\#\lambda$ and the polynomial $x\mapsto x$. Then

[TABLE]

Proof.

Let $c<\underline{f}$ be arbitrary, fixed. As $x\geq\underline{f}>c$ for all $x$ in the support of $\#\lambda$, it follows that $\mathbf{H}_{r}(x-c;\,\#\lambda)\succ 0$, and since $\mathbf{H}_{r}(\#\lambda)\succ 0$,

[TABLE]

Notice also that since $\mathbf{H}_{r}(x-c;\,\#\lambda)=\mathbf{H}_{r}(x;\,\#\lambda)-c\,\mathbf{H}_{r}(\#\lambda)$:

[TABLE]

[TABLE]

Next,

[TABLE]

and the proof is complete. ∎

4. Conclusion

We have exhibited a strong (and perhaps surprising) connection between global optimization and spectral analysis of tri-diagonal univariate moment matrices (equivalently, roots of some sequence of univariate orthogonal polynomials). Essentially computing the global minimum (resp. maximum) of a function $f$ on a compact set $\mathbf{\Omega}\subset\mathbb{R}^{n}$ reduces to a one-dimensional problem, namely computing the limit of the smallest (resp. largest) eigenvalue of a sequence of tri-diagonal moment matrices whose size $r$ is independent of the dimension $n$. Of course the entries of these matrices require computing integrals $\int_{\mathbf{\Omega}}f(x)^{k}d\lambda$, $k\in\mathbb{N}$, for some choice of a measure $\lambda$ whose support is exactly $\mathbf{\Omega}$. When $\mathbf{\Omega}$ is a simple set then this is theoretically easy but becomes tedious for large $n$.

On the theoretical side, the question of how fast such bounds converge to the global minimum was not addressed in this paper but was addressed in the recent work of Laurent and Slot [15]. In [15] the authors prove that the bounds converge to the global minimum (resp. global maximum) at a $O(\log^{2}r/r^{2})$ rate, under a weak geometric condition on $\mathbf{\Omega}$ (satisfied by all of simple sets except for the discrete hypercube).

On a more practical side, whether such an approach may provide good bounds for not too large $r$ and significant dimension $n$, is a topic of future research. Interestingly, it is definitely related to a problem in numerical analysis, namely how efficiently can be computed integrals $\int_{\mathbf{\Omega}}f(\mathbf{x})^{k}\,d\lambda$ when $f$ is a polynomial, $\mathbf{\Omega}$ is a simple set and $\lambda$ an appropriate measure whose support is $\mathbf{\Omega}$.

At last but not least, observe that if the global optimum can be approximated at a $O(\log^{2}r/r^{2})$ rate, so far one cannot approximate global minimizers as we are essentially addressing a related one-dimensional problem. Therefore our result in this paper suggests the following question: For the particular class of NP-hard polynomial optimization problems studied in this paper, is the problem of computing (only) the optimal value of same computational complexity as the problem of computing the optimal value and global minimizers ?