A new Algorithm for Overcomplete Tensor Decomposition based on Sums-of-Squares Optimisation

TL;DR

This paper introduces a novel algorithm leveraging Sums-of-Squares Optimization for overcomplete tensor decomposition, enabling the extraction of individual components from high-degree polynomials through quadratic form reduction.

Contribution

It develops a new class of algorithms based on Sums of Squares Programming to efficiently decompose overcomplete tensors by reducing polynomial degrees.

Findings

Effective component extraction via eigenvalue decomposition.

Reduction of high-degree polynomials to quadratic forms.

Potential for improved tensor decomposition accuracy.

Abstract

In this thesis, a new class of algorithms based on Sums of Squares Programming is developed. These allow to reduce a degree-$d$ homogeneous polynomial $T = \sum_{i = 1}^m \langle a_i, X \rangle^d $ to a quadratic form being close to a rank-$1$ form via a low-degree reduction polynomial $W\in\sum \mathbb{R}[X]^2$. $W$ can be thought of as a `weight function' attaining high values on merely one of the components $a_i$. The component can then be extracted by running an eigenvalue decomposition on the quadratic form $\sum_{i=1}^m W(a_i) \langle a_i, X \rangle^2$.

Equations

$$\{a_1,\dots ,a_m\}\mapsto \sum _{i=1}^m⟨a_i,X{⟩}^d$$

$$\sum _{|\alpha |=k}{\mathbb{E}}_{\mu }[Y^{\alpha }]X^{\alpha }$$

$$\sum _{|\alpha |=k}{\mathbb{E}}_{\mu }[Y^{\alpha }]X^{\alpha }=\frac{1}{m}\sum _{i=1}^m⟨a_i,X{⟩}^k$$

$$M=\sum _{i=1}^m{a}_i{a}_i^T=A{A}^T=\left(AR\right){\left(AR\right)}^T$$

$${\mathbb{E}}_{\mu }[f]$$

$${\mathbb{E}}_{\mu }[p]$$

$$\mu (\{x\})={\mathbb{E}}_{\mu }[f_x]={\mathbb{E}}_{\nu }[f_x]\stackrel{f_x\ge 1_x}{\ge }{\mathbb{E}}_{\nu }[1_x]=\nu (\{x\})\text{ for }x\in {\mathbb{R}}^n$$

$${\mathbb{E}}_{\mu }[X^{\alpha }]$$

$${\mathcal{M}}_{\mu }=\sum _{|\alpha |\le d}{\mathbb{E}}_{\mu }[X^{\alpha }]X^{\alpha }$$

$${\mathcal{M}}_{\mu }=\sum _{k=0}^d\sum _{i=1}^m⟨a_i,X{⟩}^k$$

$$\tilde{T}⟨X,Y,Z⟩=\sum _{i=1}^m⟨a_i,X⟩⟨a_i,Y⟩⟨a_i,Z⟩$$

$$M\text{mathchar}58=\tilde{T}⟨X,Y,v⟩=\sum _{i=1}^m⟨a_i,v⟩⟨a_i,X⟩⟨a_i,Y⟩=X^T\left(\sum _{i=1}^m⟨a_i,v⟩a_i{a}_i^T\right)Y$$

$$c_i=\sqrt[3]{\frac{{\mu }_i}{⟨u_i,v⟩}}{u}_i=\sqrt[3]{\frac{⟨a_i,v⟩\Vert a_i{\Vert }^2}{⟨\frac{\pm a_i}{\Vert a_i\Vert },v⟩}}{u}_i=\sqrt[3]{\pm \Vert a_i{\Vert }^3}{u}_i=\pm \Vert a_i\Vert u_i=a_i$$

$$⟨x,y{⟩}_{{\mathcal{M}}_2}\text{mathchar}58=x^T{\mathcal{M}}_2^{-1}y$$

$$⟨A{e}_i,A{e}_j{⟩}_{{\mathcal{M}}_2}={\left(A{e}_i\right)}^T(A^{-1}{)}^T{A}^{-1}\left(A{e}_j\right)=⟨e_i,e_j⟩$$

$$P_{\mu }\left[|\overline{Y^{\alpha }}-{\mathbb{E}}_{\mu }[Y^{\alpha }]|\ge \tau \right]\le \frac{{\sigma }^{2d}}{N{\tau }^2}$$

$$N\ge \frac{100{\sigma }^{2d}}{{\tau }^2}\in \mathcal{O}(\frac{{\sigma }^{2d}}{{\tau }^2})$$

$$(1)\mathbb{E}[1]=1$$

$$(2)\mathbb{E}[P^2]\ge 0\text{ for all square polynomials }{P}^2\in \mathbb{R}[X{]}_{\le d}$$

$$\mathbb{E}[X^{\otimes k}]\text{mathchar}58=\sum _{|\alpha |=k}\mathbb{E}[X^{\alpha }]X^{\alpha }\approx \sum _{i=1}^m⟨a_i,X{⟩}^k$$

$$T=\sum _{i=1}^m⟨a_i,X{⟩}^3$$

$$T^2=\sum _{i,j=1}^m⟨a_i,X{⟩}^3⟨a_j,X{⟩}^3$$

$$W_0\text{mathchar}58=\prod _{l=1}^d\frac{1}{M}⟨G^{(l)},X{⟩}^2$$

$$\mathbb{E}[W_0⟨a_j,X{⟩}^2]\ge (1-\mathcal{O}(\epsilon ))\mathbb{E}[W_0]$$

$${\mathbb{E}}_{\xi \sim \mathcal{N}(0,1)}[{\xi }^2{1}_{\{\xi \ge {\tau }_M\}}]=M$$

$$\left(\genfrac{}{}{0ex}{}{d}{\alpha }\right)\text{mathchar}58=\frac{d!}{{\alpha }_1!\cdots {\alpha }_n!}$$

$$⟨a,X{⟩}^d=\sum _{|\alpha |=d}\left(\genfrac{}{}{0ex}{}{d}{\alpha }\right)a^{\alpha }{X}^{\alpha }$$

$$P=\sum _{|\alpha |\le d}\left(\genfrac{}{}{0ex}{}{|\alpha |}{\alpha }\right)P_{\alpha }{X}^{\alpha }.$$

$$T=\sum _{i=1}^m⟨a_i,X{⟩}^d$$

$${\delta }_x\text{mathchar}58\mathcal{P}({\mathbb{R}}^n)\to {\mathbb{R}}_{\ge 0},A\mapsto \left\{\begin{array}{cc}1,& \text{x∈A}\\ 0,& \text{otherwise}\end{array}\right\}$$

$${\delta }_{ij}=\left\{\begin{array}{cc}1,& \text{if i=j}\\ 0,& \text{otherwise}\end{array}\right\}$$

Full text

1. 1 Introduction

2. 2 Overview

3. 2.1 Outline of this Thesis

4. 2.2 Results and Related Work

5. 3 Preliminaries

6. 3.1 Moment and Tensor Decompositions

7. 3.2 Tensor Notation and Operations

8. 3.3 SOS Optimisation

9. 3.4 Multilinear Algebra on the Space of Polynomials

10. 4 The v-Algorithm for exact Moment and Tensor Decomposition

11. 4.1 v-Algorithm for Moment Decomposition

12. 4.2 v-Algorithm for Tensor Decomposition

13. 4.3 Qualitative Aspects

14. 5 Efficient Tensor Decomposition via SOS

15. 5.1 Finding a Testimony

16. 5.2 Guarantees for the Optimiser

17. 5.3 Approximate Recovery of One Component

18. 5.4 Choice of $v$

19. 5.5 Recovering all Components

20. 5.6 Noise Stability and Tensor Decomposition

21. 6 Conclusions and Open Problems

22. Supplement I - Existence of Optimal Solutions

1 Introduction

What is tensor decomposition?

Tensor decomposition is the problem of finding a preimage111The $a_{i}$ can only be recovered up to reordering and double occurences  –  this is suggested by the use of set notation. It is further reasonable to assume generally that no two $a_{i}$ are multiples of each other. Note that if $d$ is even, then we must further accept that we will not of course be able to distinguish between $a_{i}$ and $-a_{i}$. Even taking these effects aside, there will be many preimages in general. in the assignment

[TABLE]

where $a_{i}\in\mathbb{R}^{n}$ are (distinct) vectors not being multiples of each other and the image $T=\sum_{i=1}^{m}\langle{a}_{i},X\rangle^{d}\in\mathbb{R}[X]_{=d}=\mathbb{R}[{X}_{1},\ldots,{X}_{n}]_{=d}$ is a homogeneous polynomial in the variable vector $X=\left({X}_{1},\ldots,{X}_{n}\right)$. The polynomial $T$ is called a tensor222The usual approach from category theory defines tensors as members of an abstract space $\underbrace{\mathbb{R}^{n}\otimes\ldots\otimes\mathbb{R}^{n}}_{d\text{ times}}$ fulfilling a universal property. It can be shown that the space of degree-$d$ homogeneous polynomials in noncommuting variables is indeed a model of this property. However, in this whole thesis we will always work with symmetric $d$-tensors. It can be shown that these correspond precisely to the homogeneous polynomials in commuting variables, i.e. the space $\mathbb{R}[X]_{=d}$ and the smallest such $m$ is called the rank of $T$. $\langle\cdot,\cdot\rangle$ denotes the scalar product of two real vectors.

To get a better intuition for the problem, think of a probability measure $\mu=\frac{1}{m}\sum_{i=1}^{m}\delta_{a_{i}}$ which is finitely supported on the (distinct) nodes $a_{i}\in\mathbb{R}^{n}$. Suppose we do not know the $a_{i}$ and the only thing we are given are samples from a vector $Y\sim\mu$ being a $\mu$-distributed random variable. Can we find out where $\mu$ is supported  –  i.e.  –  can we compute the $a_{i}$ from that information?

Note that the degree $k$-moments $\mathbb{E}_{\mu}[Y^{\alpha}]$ of $Y$ (where $\alpha\in\mathbb{N}_{0}^{n}$ is a multi-index of length $|\alpha|=k$) essentially form a tensor for all $k\in\{0,\ldots,d\}$.

[TABLE]

It can be shown that the following identity holds:

[TABLE]

by expanding the right hand side and explicitly computing the moments on the left hand side.

Now the key point is the following: If we can draw sufficiently many samples from $Y$, then we may assume that we know the left hand side (by the central limit theorems of probability theory)  –  up to a small noise term. If we can then run a tensor decomposition algorithm on the left hand side, we will indeed get the support nodes approximately out of $\mu$.

At first sight, tensor decomposition might seem like an impossible task, especially to people used to matrix decompositions. This is due to the fact that if a quadratic form $Q$ (corresponding to a symmetric matrix $M\in\mathbb{R}^{n\times n}\cong\mathbb{R}^{n}\otimes\mathbb{R}^{n}$) admits a decomposition $Q=\sum_{i=1}^{m}\langle{a_{i}},X\rangle^{2}$ (corresponding to $M=\sum_{i=1}^{m}a_{i}a_{i}^{T}$), then for most of the time this decomposition will not be unique. The reason for this is that we are able to perform an orthogonal transformation: Write $A=({a}_{1},\ldots,{a}_{m})\in\mathbb{R}^{n\times m}$. Note that $M=AA^{T}$. Then for any orthogonal matrix $R$,

[TABLE]

Hence the columns of $AR$ will, too, form a matrix decomposition of $M$. If $m,n>1$ and the $a_{i}$ are generic vectors, then the columns of $AR$ will always form another decomposition, unless $R$ was a permutation matrix (in that case only the order of the $a_{i}$ would change). Thus in general the operation $A\mapsto AR$ will yield a decomposition distinct from $\{{a}_{1},\ldots,{a}_{m}\}$.

However, there is an obvious reason why unicity fails: The second-order moments are just not enough data. If only these are known, then $\mu$ could as well be, for instance, the Gaussian measure $\mu\sim\mathcal{N}(\frac{1}{m}\sum_{i=1}^{m}a_{i},\frac{1}{m}M)$ with mean $\frac{1}{m}\sum_{i=1}^{m}a_{i}$ and psd covariance matrix $\frac{1}{m}M$.333Usually, covariance matrices are required to be positive definite. Yet there is a reasonable extension of the notion of Gaussians for the psd case.

This points us to the question: Can more data help restoring unicity? There is a natural argument suggesting this: $\mu$ is uniquely determined by the values of

[TABLE]

for all measurable and integrable functions $f$ on $\mathbb{R}^{n}$ (since $\mu$ is finitely supported, we can regard all functions $\mathbb{R}^{n}\to\mathbb{R}$ as measurable and integrable). Actually, we can replace “all functions $f$” by

[TABLE]

for all polynomials $p\in\mathbb{R}[X]$. Indeed this is true, since $\mu$ is supported on a compact set, we may approximate any function on $K\mathrel{\mathop{\mathchar 58\relax}}=\{{a}_{1},\ldots,{a}_{m}\}$ by polynomials. Due to polynomial interpolation, we may even restrict to polynomials $p\in\mathbb{R}[X]_{\leq d}$ of degree at most $d=2m$: For any $x\in\mathbb{R}^{n}$, the indicator function $\mathbb{1}_{x}$ coincides with some nonnegative polynomial $f_{x}$ on $K\cup\{x\}$, where $f_{x}$ can be chosen of degree less or equal $d$.444$2m$ is the degree of the nonnegative multivariate interpolation polynomial $\mathcal{I}_{x}\mathrel{\mathop{\mathchar 58\relax}}=\prod_{i=1}^{m}\frac{\|X-a_{i}\|^{2}}{\|x-a_{i}\|^{2}}$ satisfying $\mathcal{I}_{x}(y)=\mathbb{1}_{x}(y)$ for $y\in K\cup\{x\}$. We may thus choose $f_{x}\mathrel{\mathop{\mathchar 58\relax}}=\mathcal{I}_{x}$. At this point, $\mu$ is uniquely determined among all finitely supported measures on $\mathbb{R}^{n}$ by the values $\mathbb{E}_{\mu}[f]$ for all polynomials $f$ with $\deg(f)\leq d$:

Indeed, if $\nu$ was another finitely supported measure such that $\mathbb{E}_{\mu}[f]=\mathbb{E}_{\nu}[f]$ for all polynomials $f$ with $\deg(f)\leq d$, then

[TABLE]

Hence $\nu$ is supported on $K$, too (which implies $\mathbb{E}_{\nu}[\sum_{i=1}^{m}\mathbb{1}_{a_{i}}]=\mathbb{E}_{\nu}[1]$), and the weight $\nu$ puts on $a_{i}$ is less than $\mu(\{a_{i}\})$. Since $\sum_{i=1}^{m}\nu(\{a_{i}\})=\mathbb{E}_{\nu}[1]=\mathbb{E}_{\mu}[1]=\sum_{i=1}^{m}\mu(\{a_{i}\})$, this already implies $\mu(\{a_{i}\})=\nu(\{a_{i}\})$. Thus $\mu=\nu$.

Now $\mathbb{R}[X]_{\leq d}$ is finite-dimensional, whence we can further reduce to knowing the moments

[TABLE]

for any multi-index $\alpha$ of length $|\alpha|\leq d$. Grouping all these moments together, we see that at this point, $\mu$ is uniquely determined among all measures on $\mathbb{R}^{n}$ by the degree-$d$ polynomial

[TABLE]

Our problem then translates to finding the vectors ${a}_{1},\ldots,{a}_{m}$ such that

[TABLE]

which can again be shown by expanding the right hand side using the multinomial theorem. To be fair, this is not a tensor decomposition problem as stated in (1), but it may be seen as a “dehomogenised variant” of (1). From the above considerations, we know that we can simultaneously decompose all of those moment tensors $\mathcal{M}_{k}=\sum_{i=1}^{m}\langle{a}_{i},X\rangle^{k}$ by the same $a_{i}$ and that this simultaneous decomposition is unique. We call such a problem a moment decomposition problem. Moment decomposition and tensor decomposition are closely related problems.

One might argue that we cheated to get uniqueness by requiring multiple moments, and that the issue would look differently if we started with just one tensor. However, as we will see in §4.2, this is essentially not the case: One such tensor of degree $d\geq 2m$ suffices to recover $\pm a_{i}$. The reason for this is that, given a sole tensor of high degree, we can generate lower degree tensors $\sum_{i=1}^{m}\langle a_{i},v\rangle^{d-k}\langle{a_{i}},X\rangle^{k}$ (where $v\in\mathbb{S}^{n-1}$ is some random vector) that we can use to cast a tensor decomposition problem into a moment decomposition problem. These lower degree tensors can thus be seen as some sort of “fake moments”.555The terminology follows [BS16], where the authors considered “fake moments” of higher degree.

Undercomplete tensors

There is another case where it is very easy to see that uniqueness holds:

Proposition 1.1**.**

Suppose $T=\sum_{i=1}^{m}\langle{a_{i}},X\rangle^{3}$ is a $3$-tensor with $m\leq n$ orthogonal666This assumption can be relaxed to linear independence of the $a_{i}$ components $a_{i}\in\mathbb{R}^{n}$. Then ${a}_{1},\ldots,{a}_{m}$ is the only tensor decomposition of $T$ with $m$ or less components and we can compute it by Jennrich’s algorithm ([Har70], see algorithm 1), which is a classical result of tensor decomposition.777The algorithm was attributed to R. Jennrich in R. Harshman’s publication cited above. It seems that R. Jennrich contributed the proof of uniqueness and the basic form of the algorithm, but did not make an own publication as of our knowledge.

Proof.

Let us quickly review why this algorithm will produce correct results. Note that with probability $1$ the matrix $M$ in line 3 will have $m$ distinct eigenvalues. Indeed, due to orthogonality the $a_{i}$ are eigenvectors of $M$ corresponding to eigenvalues $\mu_{i}\mathrel{\mathop{\mathchar 58\relax}}=\langle a_{i},v\rangle\|a_{i}\|^{2}$. Since $v$ was uniformly random, hitting the lower-dimensional set $\{v\in\mathbb{S}^{n-1}\mid\langle a_{i},v\rangle\|a_{i}\|^{2}=\langle a_{j},v\rangle\|a_{j}\|^{2}\}$ has probability zero over the choice of $v$ for every two $i\neq j$ (note that the $a_{i}$ are nonzero). Hence every eigenspace of $M$ is one-dimensional and contains thus precisely two eigenvectors $\pm u$ of unit length. Since $a_{i}$ was an eigenvector of $M$, we get $u=\pm\frac{a_{i}}{\|a_{i}\|}$.

This shows that up to reordering, the eigendecomposition computed in line 4 satisfies $u_{i}=\pm\frac{a_{i}}{\|a_{i}\|}$ and $\mu_{i}=\langle a_{i},v\rangle\|a_{i}\|^{2}$. Hence the vector $c_{i}$ computed in line 5 satisfies

[TABLE]

Note that the unknown factor $\pm 1$ of $u_{i}$ cancels with the one present in $\langle u_{i},v\rangle$ below the root. Of course the eigenvalue decomposition could yield the eigenvectors in a different order. But this does not matter, since the claim $\>T=\sum_{i=1}^{m}\langle{c_{i}},X\rangle^{3}\>$ we had on the output ${c}_{1},\ldots,{c}_{m}$ does not depend at all on the order of the $c_{i}$. Uniqueness follows from the fact that the resulting $c_{i}$ are identical to the $a_{i}$ no matter what representation we chose. ∎

Note that such a procedure wouldn’t be possible if we only had access to the second order input, i.e. $\sum_{i=1}^{m}\langle{a}_{i},X\rangle^{2}$. We need that extra degree of freedom a $3$-tensor admits by allowing us to scale the eigenvalues of its matrix reductions.888The eigenvalues of a random symmetric matrix are distinct with probability $1$, but $m\leq n$ randomly chosen $a_{i}$ will almost-surely not be orthogonal. Jennrich’s Algorithm can though be extended to work for at most $n$ linearly independent components: If we get both the $2$nd order moment and the degree-$3$ moment of the $a_{i}$ with, the requirement of orthogonality for uniqueness is, as we will see, without loss of generality, contrary to the matrix case. Another issue here is the requirement of orthogonality, which allows $T$ to be just of very low rank, i.e. $m\leq n$. In the applications, though, the case $m>n$ is common. Due to the qualitative impact Jennrich’s Algorithm has, it is common practice to distinguish the following cases of tensor decomposition:

The undercomplete case: Here $m\leq n$ and thus we can hope for the $a_{i}$ to be linearly independent (if, for instance, the $a_{i}$ are generic vectors, then we’d expect them to be). If we get access to the second-order moment $\mathcal{M}_{2}=\sum_{i=1}^{m}a_{i}a_{i}^{T}$ as well (this is (approximately) the case e.g. if our tensor stems from a probability distribution from which we can take samples) we may even reduce to the orthogonal case: For simplicity, assume $m=n$.999The case $m<n$ can be dealt with by learning a projection to a lower-dimensional space (for example out of the eigendecomposition of $\mathcal{M}_{2}$) before employing this argument. Then the matrix $\mathcal{M}_{2}$ will be positive definite and the $a_{i}$ will be orthogonal w.r.t. the scalar product defined by

[TABLE]

Indeed, note that $A=({a}_{1},\ldots,{a}_{m})$ factorises $\mathcal{M}_{2}$, that is, $AA^{T}=\mathcal{M}_{2}$. Thus $\mathcal{M}_{2}^{-1}=(A^{-1})^{T}A^{-1}$ and $a_{i}=Ae_{i}$ satisfy

[TABLE] 2. 2.

The overcomplete case: Here $n<m\leq\binom{n+d-1}{d}$ holds101010In the applications, the components typically correspond to features of a dataset. To ensure that these features are actually meaningful, it is therefore recommendable to search for a decomposition with a small number $m$ of components. From that point of view, it becomes clear that the case $m\geq\binom{n+d-1}{d}=\dim(\mathbb{R}[X]_{=d})$ isn’t of much interest: Then the $\langle{a_{i}},X\rangle^{d}$ will be linearly dependent in the space $\mathbb{R}[X]_{=d}$. By Caratheodory’s theorem from convex geometry, the cone spanned by the forms $\langle{a_{i}},X\rangle^{d}$ would then also be spanned by $\dim(\mathbb{R}[X]_{=d})=\binom{n+d-1}{d}$-many forms (Cor. 7.4.21 in [RAG]). Thus it could be possible to change the size of some of the components or even make some of them vanish or reappear, which wouldn’t fit well with the intuition of them being “relevant features”. and therefore the $a_{i}$ will always be linearly dependent, even if the $\langle{a_{i}},X\rangle^{d}$ are not. In particular, we can’t obtain the $a_{i}$ as eigenvectors of the same matrix. Overcomplete tensor decomposition is usually the more interesting case, though also way more difficult.

Overcomplete tensors and the variance problem

As of now, people are trying to find time-efficient, noise stable and easily implementable algorithms for the overcomplete case. Since we have seen that the solution will get unique once the input degree $d$ is high in relation to the number of components (e.g. $d\geq 2m$ will do), it is natural to conjecture that the overcomplete case will become easier if we have the possibility to “increase” the value of $d$ (while fixing $m$). This is possible e.g. in the empirical case where we can estimate any moment from given samples. However, there is at least one problem with that:

Suppose $Y\sim\mu$ is a $\mu$-distributed random vector whose variances are bounded by some $\sigma^{2}\in\mathbb{R}_{>0}$. To generate the moment $\mathbb{E}_{\mu}[Y^{\alpha}]$ for some multi-index $\alpha$ of length $d$ by averaging, we need a number of samples which is, in general, exponential in $d$. Indeed, by Chebyshev’s law of large numbers, we have

[TABLE]

Here we denote by $\overline{Y^{\alpha}}$ the average over $N$ i.i.d. copies of the random variable $Y^{\alpha}$. If we want to get the moment $\mathbb{E}_{\mu}[Y^{\alpha}]$ up to noise of the magnitude of $\tau$ with at least $0.99$ certainty, then we need (by setting the right hand side equal to $0.01$)

[TABLE]

samples.111111The claim of exponentiality follows from the fact that Chebyshev’s bound is tight for some random variables with variance $\sigma>1$. If the data is obtained e.g. by pricey physical measurements, generating higher-order moments may literally get expensive.

History of moment generation procedures

The research of B. Barak, D. Steurer and J. Kelner [BKS15] together with R. Ge and T. Ma’s [GM15] suggests a way around that. The authors of [BKS15] proposed that, given a degree-$d$ tensor $T=\sum_{i=1}^{m}\langle{a}_{i},X\rangle^{d}$, it might be possible to generate higher-order “fake moments” of the solution vectors $a_{i}$ by optimising over the cone of pseudo-expectations.

A pseudo-expectation of degree $d$ is a linear functional $\mathbb{E}$ on the space $\mathbb{R}[X]_{\leq d}$ satisfying

[TABLE]

This optimisation can be done by a powerful tool known as sums of squares (SOS) programming, which we will discuss in §3.3. Now, if such a pseudo-expectation would have anything to do with $\mu=\sum_{i=1}^{m}\delta_{a_{i}}$, then we could hope that

[TABLE]

and thus we could try to run a noise-stable tensor decomposition algorithm on $\mathbb{E}[X^{\otimes k}]$.

In [BKS15], the authors demonstrated this technique: To get a $\theta(\varepsilon)$-approximate decomposition of a measure supported on $m$ points, they take a tensor $T$ of degree $d\in\Omega(\frac{\log(\sigma)+\tau}{\varepsilon})$ (where $\tau$ is some noise parameter. Note that this requirement still assumes that we have a sufficiently high-order tensor $T$ to begin with  –  for arbitrary accuracy $\varepsilon$ we’d need to be able to generate actual moments of arbitrary high degree). From that, they compute a degree-$k$ pseudo-expectation $\mathbb{E}$ where $k\approx\max(4d,\frac{12\log(m)}{\varepsilon})$. Then the authors applied a brute-force yet noise stable decomposition algorithm to get one of the $a_{i}$ (approximately) out of their pseudo-distribution $\mathbb{E}$. However, we encountered an issue with their decomposition algorithm (“sampling from pseudo-distributions”  –  Lemma 5.1ff in [BKS15]), which we will address in §2.2. An error in one of the key lemmas broke their decomposition algorithm for the general case. Fortunately, they gave a second algorithm working particularly for the Dictionary Learning problem (“refined sampling from pseudo-distributions” – see §7 in [BKS15]. It assumes that samples from a $\mu$-distributed vector $Y$ are given where the distribution should satisfy certain assumptions).

Such an approach can, as the authors pointed out, only work if we assume that we can generate the higher order moments by a reasonable amount of samples. This is why the authors assumed “niceness conditions” on the distribution which imply that moments of high order $d$ are known up to a noise constant $\tau$ independent of $d$. Recall that in the worst case $\tau$ would grow exponentially with $d$. Still, this shows how much can be done with high order moments and therefore the question remained whether it would still be possible to efficiently generate such higher degree fake moments even if only a 3-tensor is given. In particular, consider the case where we are left with a slightly overcomplete degree-$3$ tensor $T$ of rank $m$ in between $n$ and $n^{1.5}$

[TABLE]

Ge and Ma then showed that a quasi-polynomial time procedure can be realised for “average” degree-$3$ input tensors $T$, restricting to the case where the components are chosen randomly from an $n$-dimensional hypercube, precisely $a_{i}\in\{\pm\frac{1}{\sqrt{n}}\}^{n}$. Alas, to this end, they used the sampling procedure of [BKS15] as a key part of their algorithm without giving a proof on their own. Therefore, the only publication we know of where higher degree fake moments have successfully be used to help with the decomposition of 3-tensors is [HSS16], where the tensor

[TABLE]

is reshaped and reweighed in a sophisticated manner to generate a proxy for the moment of order 4. But this method does not even use pseudo-expectations anymore.

Despite these issues in the current meta of research, the reader should have got an idea how valuable high order moments are and that sums of squares programming can be a valuable tool in generating them. In this thesis, we will show that sums of squares programming can also be used to decompose high order tensors directly.

This is interesting because it indicates that in situations very similar to the setting of [BKS15], we can work directly on the given actual moments instead of generating higher degree fake moments. This does not bridge the gap in [BKS15] (since to this end, one would need to verify that all of the arguments used work (at least qualitatively) for the fake moments as well, which is likely not the case in the current formulation), but it essentially allows to solve moment decomposition problems with the components lying on the unit sphere when sufficient data is given.

What’s even more interesting is that if it would be possible to design two compatible sums of squares based procedures, one for the generation part and another one for the decomposition part, then this could open a whole lot of possibilities. But we are not quite there yet (and we do not even know if or for which tensors this can possibly work) and this is an interesting subject for future research.

Acknowledgements

This thesis would not have been possible without the many valuable discussions and the support of my supervisor Markus Schweighofer. Also, I wish to thank Adam Kurpisz for telling me a lot about sums of squares programming and my parents for “funding my research”.

Disclaimer

This thesis was written just by myself. When the term “we” occurs, it is either supposed to include the reader or a matter of habit.

2 Overview

2.1 Outline of this Thesis

We will start by introducing the basic notions and notations of tensor decomposition in §3.1. Section §3.2 will cover the connection between homogeneous polynomials and symmetric tensors. In §3.3 we give a very brief and rudimentary introduction to sums of squares programming.

A common technique in tensor decomposition is to use linear shrinking maps that take a high order tensor and reduce it to a lower order tensor. In our case, we work a lot with “matrix reductions”, where the initial tensor $T=\sum_{i=1}^{m}\langle{a}_{i},X\rangle^{d}$ is reduced to a weighted quadratic form $\sum_{i=1}^{m}W(a_{i})\langle{a_{i}},X\rangle^{2}$. In §3.4 we introduce several such linear “shrinking maps” connected to polynomial evaluation that will allow us to write down the linear constraints of our sums of squares programmes.

Done with the preliminaries, we will start §4 by giving an algorithm that can compute the exact solution of a moment decomposition problem in exponential time from the first $d\geq 2m$ moments. In §4.2, we will then show how this algorithm can be adapted to work for tensor decomposition by feeding the original algorithm with “fake moments” of lower degree.

The algorithms of Section §4 aren’t very efficient neither with respect to computation time nor with respect to the amount of data needed, but they illustrate some of the main ideas and techniques that we are going to use in the approximate setting of Section §5: Here we investigate what approximation results the same kind of algorithms can achieve when we restrict to polynomially sized SOS programmes.

2.2 Results and Related Work

Based on ideas and similar concepts present in [BKS15], we develop a new class of algorithms for tensor decomposition that can be seen as a generalisation of Jennrich’s Algorithm to the case of overcomplete tensor decomposition. We show exact recovery guarantees in the case that the input tensor $T$ is of sufficiently high order ($d\geq 2m$ with $m$ being the number of components).

These algorithms are in some sense “matrix reduction algorithms”, since they reduce high order tensors to quadratic forms which correspond to symmetric matrices and then recover the components by running eigenvalue decompositions on the reduced matrices. The broad concept of such matrix reductions is very old – note that Jennrich’s classical result Alg. 1 from the introductory section can actually be seen an example of such a matrix reduction algorithm – at least in the broad sense: In Jennrich’s case we reduced the input tensor via the polynomial $W\mathrel{\mathop{\mathchar 58\relax}}=\langle{v},X\rangle$ and then performed an eigenvalue decomposition.

The main conceptual novelty is in showing that it’s possible to use sums of squares programming in order to find such $W$ which attains high values only on one of the components $a_{i}$ (and acts thus as some sort of a weight function on the $a_{i}$), which yields a quadratic form being approximately of rank 1. The sums of squares condition is needed to ensure that the weights are nonnegative. It will turn out that this nonnegativity condition for the values $W(a_{i})$ is essentially what enables us to find such $W$ via optimisation, since it allows to cap the maximum weight by a simple linear constraint, e.g. $\sum_{i=1}^{m}W(a_{i})=1$. This may sound quite different from Jennrich’s Algorithm, since there we needed no such thing as an SOS constraint on the reduction polynomial. However, we will see in §4.3 that Jennrich’s classical Algorithm admits an equivalent SOS based formulation.

The authors of [BKS15] were already using techniques which can, in our terminology, be described as matrix reduction via SOS polynomials: After computing some degree $k\geq d$ pseudo-distribution $\mathbb{E}$ satisfying certain constraints (in particular that $\mathbb{E}[\langle{a_{j}},X\rangle^{k}]$ is not too small for some $j\in[m]$), they wanted to reduce $\mathbb{E}$ to a degree 2 pseudo-distribution (which corresponds to a psd matrix) via conditioning with a sum of squares polynomial $W$. The authors hoped to show that with sufficiently high probability the polynomial

[TABLE]

which is a product of $d$ squares of independent Gaussian linear forms $\langle G^{(l)},X\rangle$ given by Gaussian vectors $G^{(l)}\sim\mathcal{N}(0,I_{n})$, would satisfy

[TABLE]

where $\varepsilon$ is an (unknown) approximation constant depending on $d,m$ and the condition of the problem. This was Lemma 5.2 in [BKS15].

Unfortunately though, there was an error in the proof of Lemma 5.2: In lines 2 - 3 on page 19, they choose pairs $(\tau_{M},M)\in\mathbb{R}^{2}$ of real numbers such that a standard (expected value 0 and standard deviation 1) Gaussian scalar variable $\xi$ satisfies

[TABLE]

Of course this is feasible precisely for any $M\in(0,1)$.

On the same page, in line 21, they chose $M=(1/\varepsilon)\cdot\log(1/\varepsilon)\gg 1$ for some very small $\varepsilon>0$. Therefore $1>M\gg 1$, rendering the choice of $M$ infeasible. This essentially breaks the proof. The statement also has a surprising qualitative aspect, since we would expect that the choice of the scaling factor $M$ would not make any difference on the quality of the estimation (4).

Even more unfortunate is that Ge and Ma cited particularly this Lemma without repeating the proof in [GM15]. Using Lemma 5.2 as an integral part of their decomposition algorithm, they designed the procedure which we described in the introductory chapter (1) and which achieves the currently best known recovery guarantees for average case tensors (in the sense that they achieve quasipolynomial time while allowing $m$ to be quite large, almost $m=n^{1.5}$).

3 Preliminaries

Notation 3.1**.**

Vectors of polynomial unknowns and random vectors will always be denoted by capital letters such as $X,Y$, whereas variables with numerical values such as $x,y\in\mathbb{R}^{n}$ will be denoted by lower case letters. In the case of polynomial unknowns, we implicitly understand that $X=\left({X}_{1},\ldots,{X}_{n}\right)$, where the ${X}_{1},\ldots,{X}_{n}$ are algebraically independent (scalar) unknowns. A single scalar unknown (for univariate polynomials) will be denoted by $\Lambda$. We write $\mathbb{N}=\{1,2,3,\ldots\}$ and $[m]\mathrel{\mathop{\mathchar 58\relax}}=\{1,\ldots,m\}$ for $m\in\mathbb{N}$. $e_{i}\in\mathbb{R}^{n}$ will always denote the $i$-th unit vector, such that $(e_{i})_{j}=\begin{Bmatrix}1,&\text{if i = j}\\ 0,&\text{otherwise}\end{Bmatrix}=\delta_{ij}$. $\|A\|_{\operatorname{spec}}\mathrel{\mathop{\mathchar 58\relax}}=\sup_{x\in\mathbb{S}^{n-1}}\langle Ax,Ax\rangle^{1/2}$ will denote the spectral norm of the matrix $A$. For $d\in\mathbb{N}_{0}$, let $\mathbb{R}[X]_{\leq d}$ denote the space of all polynomials in $X=({X}_{1},\ldots,{X}_{n})$ up to degree $d$ and $\mathbb{R}[X]_{=d}$ the subspace of all homogeneous polynomials of degree precisely $d$. For $\alpha\in\mathbb{N}_{0}^{n}$, we denote the multinomial coefficients

[TABLE]

occurring in the the important multinomial theorem: For $a\in\mathbb{R}^{n}$ we have:

[TABLE]

Notation 3.2**.**

Each polynomial may be represented by its coefficients w.r.t. the canonical basis $(X^{\alpha})_{|\alpha|\leq d}$, where $\alpha\in\mathbb{N}_{0}^{n}$ denotes a multi-index. We, however, use the convention to write a polynomial $P$ as

[TABLE]

Hence, we represent a polynomial w.r.t. the basis consisting of all scaled monomials $\binom{|\alpha|}{\alpha}X^{\alpha}$. This convention is technical convenience  –  mainly for compatibility with the multinomial theorem and polynomial evaluation. Details will become clear from §3.1, §3.2 and §3.4.

Remark 3.3**.**

During this whole thesis, we will use the term “tensor” both to describe a multilinear map $({v}_{1},\ldots,{v}_{d})\mapsto T\langle{v}_{1},\ldots,{v}_{d}\rangle\in\text{Hom}(\mathbb{R}^{n},\ldots,\mathbb{R}^{n};\mathbb{R})$ and a homogeneous polynomial $T=\sum_{|\alpha|=d}\binom{d}{\alpha}T_{\alpha}X^{\alpha}\in\mathbb{R}[X]_{=d}$. For most of the time, we will work with polynomials instead of multilinear maps. We will justify this and explain the connection in §3.2.

3.1 Moment and Tensor Decompositions

Definition 3.4**.**

(Tensor decomposition)  Let $T\in\mathbb{R}[X]_{=d}$ for some $d\in\mathbb{N}$. ${a}_{1},\ldots,{a}_{m}$ is called an (unweighted symmetric) tensor decomposition of $T$, if

[TABLE]

The smallest $m\in\mathbb{N}\cup\{\infty\}$ for which such a tensor decomposition exists, is called the rank of $T$ (where the rank is $\infty$ iff no tensor decomposition exists) and the $a_{i}$ are called components of the decomposition. $\|a_{i}\|$ is called the magnitude of $a_{i}$ in the decomposition.

Remark 3.5**.**

When talking about tensor decompositions, we will implicitly assume that the $a_{i}$ are all nonzero and no nonnegative multiples of each other. This makes sense for otherwise we could group $a_{i}$ and $\lambda a_{i}$ together to one summand $\sqrt[d]{(1+\lambda)}a_{i}$. If $d$ is odd, then we can extend this assumption to the $a_{i}$ being no real multiples of each other. Slightly imprecise, we will often call the $a_{i}$ components of $T$. This practice is justified if there exists only one tensor decomposition (which is the usually the case for low-rank decompositions of high degree tensors). However, we will also use this practice when we merely hope that there could be uniqueness. Note that for even $d$ it’s not possible to distinguish between components $a_{i}$ and $-a_{i}$.

Notation 3.6**.**

Let $x\in\mathbb{R}^{n}$. Then by $\delta_{x}$ we denote the Dirac measure of $x$, that is

[TABLE]

This is not to be confused with the Kronecker delta: For $i,j\in\mathbb{N}$ we denote

[TABLE]

Definition 3.7**.**

Let $\mu$ be a measure on some sub-sigma algebra $\mathcal{A}\subseteq\mathcal{P}(\mathbb{R}^{n})$ such that $\mathbb{E}_{\mu}[P]$ exists for every homogeneous polynomial of degree $k$. The polynomial

[TABLE]

is called the $k$-th moment tensor of $\mu$ (or, for brevity, the $k$-th moment).121212For technical reasons that will become apparent from the following proposition, our definition of moments is slightly different than in the introduction: The summands are now rescaled by some multinomial coefficients.

Proposition 3.8**.**

Let $\mu=\sum_{i=1}^{m}\lambda_{i}\delta_{a_{i}}$ be a measure finitely supported on the set $\{{a}_{1},\ldots,{a}_{m}\}$ with weights $\lambda_{i}>0$. Then for any $k\in\mathbb{N}$ the $k$-th moment exists and satisfies

[TABLE]

Proof.

Note that

[TABLE]

Furthermore, according to the multinomial theorem we know for each $a\in\{{a}_{1},\ldots,{a}_{m}\}$ that

[TABLE]

By expanding the right hand side of (5), we get

[TABLE]

where the last step used (6). The right hand side of equation (7) is precisely $\mathcal{M}_{k}$. ∎

Hence $\sqrt[k]{\lambda_{1}}a_{1},\ldots,\sqrt[k]{\lambda_{m}}a_{m}$ is a tensor decomposition of $\mathcal{M}_{k}$ for each $k\in\{0,\ldots,d\}$. This motivates the following definition.

Definition 3.9**.**

(Moment decomposition)  Let $T_{0},\ldots,T_{d}$ be given Tensors of orders $k_{0}<\ldots<k_{d}$. We call ${\lambda}_{1},\ldots,{\lambda}_{m},{a}_{1},\ldots,{a}_{m}$ a moment decomposition of $T_{0},\ldots,T_{d}$, if for $j\in\{0,\ldots,d\}$

[TABLE]

For simplicity, we will usually write $\lambda_{1}a_{1},\ldots,\lambda_{m}a_{m}$ instead of ${\lambda}_{1},\ldots,{\lambda}_{m},{a}_{1},\ldots,{a}_{m}$. Note that this is slight abuse of notation, since the $a_{i}$ and $\lambda_{i}$ are interpreted as distinct variables.

Moment decompositions may be seen as simultaneous (weighted) decompositions of several tensors at once. As the name suggests, the most important special case is when the first $d$ moment tensors of a measure are given:

Proposition 3.10**.**

Let $\mu=\sum_{i=1}^{m}\lambda_{i}\delta_{a_{i}}$ be a measure finitely supported on the set $\{{a}_{1},\ldots,{a}_{m}\}$ with weights $\lambda_{i}>0$. Then $\lambda_{1}a_{1},\ldots,\lambda_{m}a_{m}$ is a moment decomposition of $\mathcal{M}_{0},\ldots,\mathcal{M}_{d}$ for any $d$.

Proof.

This was shown already in Prop. 3.8. ∎

Definition 3.11**.**

Let $\mu$ a measure on some sub-sigma algebra $\mathcal{A}\subseteq\mathcal{P}(\mathbb{R}^{n})$ such that its first $d$ moments $\mathcal{M}_{0},\ldots,\mathcal{M}_{d}$ exist. Using terminology from the theory of numerical integration, a moment decomposition of $\mathcal{M}_{0},\ldots,\mathcal{M}_{d}$ is called a quadrature formula for $\mu$ up to degree $d$. This is due to the fact that for any such decomposition $\lambda_{1}a_{1},\ldots,\lambda_{m}a_{m}$ we have

[TABLE]

for all polynomials of degree up to $d$. The largest $d\in\mathbb{N}\cup\{\infty\}$ for which the actual integral with respect to $\mu$ can be replaced by this quadrature rule is called the exactness degree of the quadrature formula $\lambda_{1}a_{1},\ldots,\lambda_{m}a_{m}$ for $\mu$.

If $\mu$ is a finitely supported measure131313In that case we can always assume $\mathcal{A}=\mathcal{P}(\mathbb{R}^{n})$ and suppress the sigma-algebra in the notation., then there exists a quadrature formula of infinite exactness degree and our objective is to find it. This case will be the focus of this thesis.

Remark 3.12**.**

An obvious difference between tensor decomposition and moment decomposition is that the first problem does not ask for weights. This is due to the fact that for $\lambda_{i}\in\mathbb{R}_{\geq 0}$, $a_{i}\in\mathbb{R}^{n}$ we have that

[TABLE]

Thus, there is an ambiguity in the input data: We will never know if the component we search for is $a_{i}$ and the corresponding weight is $\lambda_{i}$ or if the component is $\sqrt[d]{\lambda_{i}}a_{i}$ with weight $1$. From that it becomes clear that we need more information if we want to distinguish between these cases. One way is the above method to require multiple moments of our solution vectors. Another way would be to impose constraints on the $a_{i}$ that restrict where the $a_{i}$ can be located (such as $\|a_{i}\|=1$).141414For even $d$, we will see in §4.3 that this particular constraint is equivalent to knowing all even-degree moments of degree lower than $d$

Let us define the corresponding notion of magnitude for quadrature formulae.

Definition 3.13**.**

Let $\lambda_{1}a_{1},\ldots,\lambda_{m}a_{m}$ a quadrature formula for the measure $\mu$. Then we call $\lambda_{i}a_{i}$ the $i$-th weighted component and $\lambda_{i}\|a_{i}\|$ its magnitude.

In some sense, moment decomposition seems to be the easier task, since the additional information we have breaks for instance the input ambiguity described above. On the other hand, tensor decomposition seems to be the more general task, since it relies on less information. Indeed, once we solved the tensor decomposition problem for the moment tensor $\mathcal{M}_{d}$ of even degree $d$, then we can search for each $k$ for some weights $\lambda_{k,i,\pm}$ minimising

[TABLE]

by solving a linear least squares problem.151515Here, $\|P\|_{F}$ denotes the Frobenius norm of $P$, see Def. 3.16. Note there is a problem when starting with some odd degree $d$: If the components $a_{i}$ and $-a_{i}$ occur both with the same weight in $\mu$, they cancel in the odd-degree tensors, whence they are lost. It is not too hard to see that you can get the true weights and the true component lengths by comparison of the $\lambda_{k,i,\pm}$.

Surprisingly, it works also the other way around: Once we have an algorithm that can solve the moment decomposition problem, but are given only one tensor of degree $d$, then we can generate some fake moments of lower degree and feed them to the moment decomposition algorithm. We’ll discuss methods to generate artificial lower degree moments, given a sole tensor $T$, in Section §4.2.

We will go with this approach and thus first design an algorithm for moment decomposition, that we will later adapt to work for tensor decomposition as well. But first, here are some more basic facts about tensor decompositions.

Remark 3.14**.**

1. (a)

For even $d$, if $T\in\mathbb{R}[X]_{=d}$ admits a tensor decomposition, then this implies that $T\geq 0$ on all of $\mathbb{R}^{n}$. 2. (b)

There is no efficient way (i.e. no procedure running in polynomial time) to decide the rank of a general $3$-tensor (assuming $P\neq NP$). 3. (c)

In fact, many problems that admit efficient algorithms in the matrix case turn out to be NP-hard when generalised to the tensor case. In **[HL13]**, C. Hillar and L.-H. Lim pointed out that “Most tensor problems are NP-hard” (which was also the title of their paper)

The results of [HL13] do not tell us whether or not there could be efficient approximation algorithms for (special cases of) the tensor decomposition problem. To be able to analyse this, let us first define how we measure the quality of an approximate solution as well as all of the related notions.

Definition 3.15**.**

(Hausdorff distance) The finite sets $M,N$ (w.l.o.g. with $\#M\geq\#N$) are called $\varepsilon$-close, if there exists a surjective map $\sigma\>\mathrel{\mathop{\mathchar 58\relax}}M\to N$ such that

[TABLE]

for all $a\in M$. The smallest such $\varepsilon\in\mathbb{R}_{\geq 0}$ is called the Hausdorff distance of $M$ and $N$. We denote it as

[TABLE]

For the special case that the sets $M=\{{a}_{1},\ldots,{a}_{m}\}$, $N=\{{b}_{1},\ldots,{b}_{m}\}$ have the same cardinality $m=\#M=\#N$, $\varepsilon$-closeness is equivalent to

[TABLE]

That is, such sets are $\varepsilon$-close if we can match them to pairs $(a_{i},b_{\sigma(i)})$ of “neighbours” such that $\|a_{i}-b_{\sigma(i)}\|<\varepsilon$.

Definition 3.16**.**

(Frobenius Norm) Let $P=\sum_{|\alpha|\leq d}\binom{|\alpha|}{\alpha}P_{\alpha}X^{\alpha}\in\mathbb{R}[X]_{\leq d}$. Then

[TABLE]

is called the Frobenius norm of $P$. It corresponds to the 2-norm of $P$ when $P$ is seen as a vector of its coefficients with respect to the basis $(\sqrt{\binom{|\alpha|}{\alpha}}X^{\alpha})_{|\alpha|\leq d}$.

Definition 3.17**.**

(Forward and Backward Error) Let $T\in\mathbb{R}[X]_{=d}$ a tensor and ${b}_{1},\ldots,{b}_{k}\in\mathbb{R}^{n}$ an approximate solution to the tensor decomposition problem for $T$. Then we call

[TABLE]

the backward error of the approximate solution ${b}_{1},\ldots,{b}_{k}$. Furthermore, with respect to an exact solution ${a}_{1},\ldots,{a}_{m}$ (such that $T=\sum_{i=1}^{m}\langle{a}_{i},X\rangle^{d}$), we define

[TABLE]

as the forward error of ${b}_{1},\ldots,{b}_{k}$ at ${a}_{1},\ldots,{a}_{m}$. The nomenclature is historical convention. Note that the latter notion will usually depend on the choice of ${a}_{1},\ldots,{a}_{m}$. However, if the true decomposition happens to be unique, then the forward error is a property of ${b}_{1},\ldots,{b}_{k}$ and $T$, too.

Our approximation algorithms in §5 will focus on minimising the forward error of the approximate solution with respect to all exact solutions which works of course best if the solution is unique.

Condition

A widespread concept in numerics is condition. This notion captures the fact that the quality of the results that an approximation algorithm produces can depend gravely on some properties of the input data / the exact solution. In fact, condition can be seen as a notion of “the exact solution’s quality”. Specifically in tensor decomposition, problematic instances could be e.g. such where some components sit very closely together. To deal with this, we will briefly introduce some parameters that will naturally appear in our following computations. Notice that all of these notions are defined with respect to a fixed set of solution vectors. In the case where the decomposition is not unique, we do not know of a canonical definition of these parameters with respect to just $T$.

Notation 3.18**.**

Let $({a}_{1},\ldots,{a}_{m})$ a list of $m$ distinct vectors in $\mathbb{R}^{n}$. We define:

[TABLE]

We’ll hide their dependency of $({a}_{1},\ldots,{a}_{m})$ in the notation, but they will always be defined with respect to $({a}_{1},\ldots,{a}_{m})$.

In the case that all the components have unit length,

[TABLE]

is one minus the maximum squared correlation between two components, which can be seen by explicitly computing the characteristic polynomial of a rank 2 matrix (the characteristic polynomial of a rank 2 matrix can be computed by looking at its trace and at the sum of its $2\times 2$ principal minors). These parameters all have in common that they measure how well-separated the components are: This can either be done by looking at the length of their differences $\|a_{i}-a_{j}\|^{2}$ or by looking at their correlation $1-\frac{\langle a_{i},a_{j}\rangle}{\|a_{i}\|\|a_{j}\|}$ together with their length differences $\|a_{j}\|^{2}-\|a_{i}\|^{2}$. Both approaches are connected via the identity

[TABLE]

For convenience, we defined all of these parameters since they will naturally appear in our estimations in §5. The correlation metric can be defined slightly more generally:

Definition and Proposition 3.19**.**

(Correlation Metric) Let $x,y\in\mathbb{R}^{n}$. We define

[TABLE]

as the correlation metric which is a metric on the unit sphere, and

[TABLE]

as the square correlation metric which is a metric on the unit sphere modulo $\pm 1$, that is, on the real projective space $\mathbb{P}_{\mathbb{R}}^{n-1}$.

Proof.

For $x,y\in\mathbb{S}^{n-1}$, note that $\|x-y\|^{2}=\|x\|^{2}+\|y\|^{2}-2\langle x,y\rangle=2(1-\langle x,y\rangle)=2\cdot d_{\mathbb{S}^{n-1}}(x,y)^{2}$. This shows that $d_{\mathbb{S}^{n-1}}$ is a metric and hence the first claim. For the second claim, recall that for unit vectors we have

[TABLE]

showing the triangular inequality and whenever $1-\langle x,y\rangle^{2}=0$ then by Cauchy-Schwarz $x=\pm y$. ∎

3.2 Tensor Notation and Operations

In the introductory chapter, we defined tensors in a way that the reader might not be familiar with. In fact, the much more common approach to tensors is to see them as multilinear maps. In this section, we will state this second, more standard definition of tensors and show that both viewpoints are equivalent for symmetric tensors.

Definition 3.20**.**

An (order-$d$) tensor is a multilinear map

[TABLE]

*We denote $T\langle{v}_{1},\ldots,{v}_{d}\rangle$ instead of $T({v}_{1},\ldots,{v}_{d})$ to stress the linearity in each of the arguments surrounded by angular brackets.

$T$ is called (totally) symmetric, if $T\langle{v}_{1},\ldots,{v}_{d}\rangle=T\langle v_{\sigma(1)},\ldots,v_{\sigma(d)}\rangle$ for all permutations $\sigma\in S_{n}$.*

Definition 3.21**.**

An order-$d$ tensor $T$ is called simple, if it is the product of $d$ linear forms, that is, if there exist linear maps $l_{i}\>\mathrel{\mathop{\mathchar 58\relax}}\>\mathbb{R}^{n}\to\mathbb{R}$ such that

[TABLE]

for all ${v}_{1},\ldots,{v}_{d}\in\mathbb{R}^{n}$. We denote such a simple tensor by $l_{1}\otimes\cdots\otimes l_{d}$

Every Tensor can be written as a linear combination of simple tensors, since the $n^{d}$ tensors given by

[TABLE]

form a basis of the space of all linear maps $\mathbb{R}^{n}\times\ldots\times\mathbb{R}^{n}\to\mathbb{R}$, which we denote by $\text{Hom}(\mathbb{R}^{n},\ldots,\mathbb{R}^{n};\mathbb{R})$ (while the subspace of all symmetric tensors is denoted by $\text{Sym}(\mathbb{R}^{n},\ldots,\mathbb{R}^{n};\mathbb{R})$). We write $w=(i_{1},\ldots,i_{d})\in[n]^{d}$ and call $w$ a word made of the “letters” $i_{1},\ldots,i_{d}$. With respect to this basis, we can write each tensor as a $d$-fold indexed array, by looking at the coordinate representation

[TABLE]

Also, we can define $X^{w}\mathrel{\mathop{\mathchar 58\relax}}=X_{i_{1}}\cdot\ldots\cdot X_{i_{d}}$ for each word $w=(i_{1},\ldots,i_{d})\in[n]^{d}$.

Proposition 3.22**.**

Let $T$ be a simple symmetric tensor. Then $T$ or $-T$ is the $d$-th power of some linear form, i.e.

[TABLE]

Proof.

Let $T=l_{1}\otimes\cdots\otimes l_{d}$. Fix some $v$ and consider the linear form

[TABLE]

which is a multiple of the form $l_{d}$. We can clearly choose $v$ in such a way that $l_{i}\langle v\rangle\neq 0$ for $i\in\{1,\ldots,d\}$, since these linear forms will be nonzero and thus vanish on lower-dimensional subspaces.161616It is well known that $\mathbb{R}^{n}$ is not a union of finitely many lower-dimensional subspaces. Then the above form is even a nonzero multiple of $l_{d}$. Now by symmetry we may exchange the arguments of $T$ and get that $\varphi$ is exactly the same as

[TABLE]

whence $\varphi$ is a nonzero multiple of $l_{d-1}$, too. Therefore, $l_{d-1}$ must be a nonzero multiple of $l_{d}$. By repeating this argument, we may see that actually all of the $l_{i}$ must be nonzero multiples of each other. We can therefore choose scalar multiples ${\lambda}_{1},\ldots,{\lambda}_{d-1}$ such that $l_{i}=\lambda_{i}l_{d}$ for $i\neq d$. Set $a\mathrel{\mathop{\mathchar 58\relax}}=\frac{l_{d}}{\sqrt[d]{|\prod_{i=1}^{d-1}\lambda_{i}|}}$. Then

[TABLE]

just as we claimed. ∎

Definition 3.23**.**

The smallest number $m\in\mathbb{N}\cup\{\infty\}$ for which a symmetric tensor $T$ admits a representation

[TABLE]

as a sum of $m$ $d$-th (tensor) powers of linear forms is called the (symmetric) rank of $T$.

Please be aware that the argument of Proposition 3.22 works only if the tensor really is simple. In a general symmetric tensor, we can usually not exchange a non-symmetric simple summand by $\pm 1$ times a power of a linear form. Requiring the factors of $\pm 1$ to be $1$ is a serious restriction, too, in particular since for even $d$ it implies that the tensor is positive semidefinite, that is, $T\langle{v}_{1},\ldots,{v}_{d}\rangle\geq 0$ for all $({v}_{1},\ldots,{v}_{d})\in\mathbb{R}^{n}\times\ldots\times\mathbb{R}^{n}$. Therefore the concept of rank we use in this thesis might or might not be different from what the reader could have read before in other works on tensor decomposition: Some people define the rank of a tensor $T$ as the smallest number of simple tensors such that $T$ can be written as a linear combination of them.

Due to these restrictions, a tensor decomposition does not need to exist in general. If no tensor decomposition exists, we understand the rank to be infinite. Still, these particular decompositions of symmetric tensors are an interesting object of study: As we saw in the introductory chapter, moment tensors of finitely supported measures will always admit a tensor decomposition.

Now, let us just briefly point out the connection between symmetric multilinear forms and homogeneous polynomials that we both called tensors.

A multi-index $\alpha$ with $|\alpha|=d$ corresponds to an equivalence class of $d$ indices $({i}_{1},\ldots,{i}_{d})$ modulo permutation. In particular, when writing down a polynomial in its coefficient representation, the coefficients of the monomials $X^{\alpha}$ correspond to the “upper triagonal” entries $T_{({i}_{1},\ldots,{i}_{d})}$ of a traditionally indexed tensor (where by upper triagonal we mean the entries corresponding to words $({i}_{1},\ldots,{i}_{d})$ such that $i_{1}\leq\ldots\leq i_{d}$  –  note that a symmetric tensor is determined by its upper triagonal). Thus there is a linear map

[TABLE]

which gets bijective when restricted to $\text{Sym}(\mathbb{R}^{n},\ldots,\mathbb{R}^{n};\mathbb{R})$. Indeed, if $T$ is symmetric, the coefficients corresponding to $w$ and $w^{\prime}$ are the same if the words $w$ and $w^{\prime}$ are permutations of each other, which is the case if and only if $X^{w}=X^{w^{\prime}}$. From that we see that

[TABLE]

for any word $w^{\prime}$ satisfying $X^{w^{\prime}}=X^{\alpha}$. Here we denote by $\binom{d}{\alpha}$ the number of words $w$ satisfying $X^{w}=X^{\alpha}$. It can be shown that

[TABLE]

(thus this is consistent with our definition in 3.1). For that reason, we call $\binom{d}{\alpha}$ the multinomial coefficients. Since $\binom{d}{\alpha}>0$, the identity (15) shows injectivity as well as surjectivity, whence the two spaces are isomorphic.

Under the map $\varphi$, simple tensors are mapped in the following way

[TABLE]

This implies that the notions of this section and §3.1 correspond.

Viewing tensors as polynomials has several advantages over the traditional representation as multilinear forms or arrays with $d$ indices, respectively. One is that we can add two tensors of different orders, since they’re both members of the polynomial ring. Another one is that the tensor product of two simple tensors $\langle a,X\rangle^{k}$ and $\langle{b},X\rangle^{l}$ becomes the regular product $\langle a,X\rangle^{k}\langle{b},X\rangle^{l}$ of two polynomials. However, the most important advantage is for sure that polynomials are the native environment for SOS optimisation (see next section).

3.3 SOS Optimisation

Sums of squares optimisation is a powerful tool that allows us to check whether or not a polynomial $W$ can be written as a sum of squares of polynomials. Furthermore, it allows us to optimise over all such sums of squares (SOS) polynomials $W$ satisfying certain linear constraints.

As such, it is a generalisation of linear programming (LP) and it serves, in some sense, as an approximation to the condition that a polynomial is nonnegative. The main novelty that SOS programming brings is that it allows to solve optimisation problems of the form:

[TABLE]

for some function $f$ being linear in the coefficients of the polynomial $W$ and some linear subspace $\mathcal{L}\subseteq\mathbb{R}[X]_{\leq d}$. Here, $\sum\mathbb{R}[X]_{\leq\frac{d}{2}}^{2}$ denotes the convex cone of all SOS polynomials of degree $\leq d$. Here, $d$ is called the degree of the programme.

The motivation behind this is that we want to optimise over polynomials which satisfy a global nonnegativity constraint $W\geq 0$. Alas, nonnegative polynomials are hard to optimise over, so we need to replace the nonnegativity condition $\geq 0$ by something which is easier to verify:

Definition 3.24**.**

Let $P,Q\in\mathbb{R}[X]_{\leq d}$. We write

[TABLE]

Of course each sums of squares polynomial will be globally nonnegative. The converse direction is in general not true, however, there are some special cases where we know by classical results dating back to D. Hilbert that equivalence holds between “$\geq 0$” and “$\succeq 0$” (e.g. 2.1.1, 2.3.2 and 2.3.5 in [RAG]).

The first case are univariate polynomials, which is really easy to see by computing the polynomial factorisation over $\mathbb{C}$ and then grouping the non-real factors by pairs of complex conjugates. The remaining real roots then have to be two-fold due to nonnegativity. The second case is the one of degree two polynomials in an arbitrary number of variables. Here the equivalence can be proven by diagonalisation of symmetric matrices (after homogenising). The reason why the sums of squares condition is so much easier to check is that it can be reduced to the task of finding some psd matrix:

Notation 3.25**.**

Let $x\in\mathbb{R}^{n}$. We denote

[TABLE]

Furthermore, let $\mathfrak{m}_{\leq d}\mathrel{\mathop{\mathchar 58\relax}}=\{X^{\alpha}\mid\>{|\alpha|\leq d}\}$. With respect to any fixed ordering of this set of cardinality $\binom{n+d}{d}$, we may identify $\mathfrak{m}_{\leq d}\equiv[\binom{n+d}{d}]$ and treat elements of $\mathbb{R}^{\mathfrak{m}_{\leq d}\times\mathfrak{m}_{\leq d}}$ as matrices.

Definition 3.26**.**

Let $W\in\mathbb{R}[X]_{\leq d}$. We call $G=(G_{\alpha,\beta})\in\mathbb{R}^{\mathfrak{m}_{\leq d/2}\times\mathfrak{m}_{\leq d/2}}$ a Gram matrix of $W$, if

[TABLE]

Since any monomial of degree less or equal to $d$ is a product of two monomials of degree less or equal to $\frac{d}{2}$, each polynomial has a Gram matrix representation. Gram matrices are in general highly non-unique since for instance any monomial of degree greater or equal than 2 can be written in at least two different ways. Indeed, we can write the monomial as $X_{i}\cdot X^{\alpha}$ for some $i\in[n]$ and some $\alpha\neq(0,\ldots,0)$ but also as $1\cdot\left(X_{i}X^{\alpha}\right)=X^{(0,\ldots,0)}\cdot X^{\beta}$ for $\beta\mathrel{\mathop{\mathchar 58\relax}}=\alpha+e_{i}$. Hence the coefficient of this monomial in $W$ may be distributed over several entries of $G$.

It can be shown that a polynomial is a sum of squares if and only if it has a positive semidefinite Gram matrix (see for instance §2.6 in [RAG]). This condition can be encoded within a semidefinite programme (SDP) and then be checked by an SDP solver. The condition that $G$ is a Gram matrix of $W$ can likewise be encoded by adding linear constraints between the coefficients of $G$ and the coefficients of $W$.

Sums of squares programmes can be, in some sense, efficiently solved numerically by SDP solvers. However, there are some caveats related to feasibility and the coefficients which one would need to take care of in order to get precise formulations of the guarantees SOS solvers give. For simplicity, we will not do this but instead formulate what we call the “magical black box rule” of SOS programming:

Magical Black Box Rule of SOS Programming:

An $n$ variate SOS programme of degree $d$ can be solved in time $n^{\theta(d)}$. 171717By “solved” we mean that if $f$ is the linear objective function and $W^{\ast}$ is some optimal solution then to any fixed numerical accuracy $\gamma>0$ we can find in time $n^{\theta(d)}$ a feasible polynomial $W$ which achieves $|f(W^{\ast})-f(W)|<\gamma$, where the hidden constants in the Landau notation may depend on $\gamma$. For simplicity, we will pretend that this would work for $\gamma=0$ as well.

Now, this rule is wrong: There exist at least the following caveats:

1. (I)

There is an issue with the exact complexity of general feasibility problems, that we will not get into. 2. (II)

We ignored the encoding length of the coefficients we used to write down the linear constraints. Also, it is recommendable to have a bound on the diameter of the feasible space. 3. (III)

The number of linear constraints has an impact on the running time, but following P. Parillo [Par00], it is usually ignored since the cost of the linear constraints is dominated by the cost of the SOS constraints.

However, for most cases these will not be an issue and therefore people usually pretend that you can solve an SOS programme of fixed degree $d$ exactly in polynomial time.

In order to not go beyond the scope of this thesis, we will leave it with this very coarse introduction to SOS programming. It’s needless to say that we merely scratched the surface of a huge topic. There is a long and detailed survey of M. Laurent [Lau09] that we recommend to anyone interested in the details. Since we did not give a theorem with rigid runtime guarantees, we will refrain from claiming any such runtime guarantees (in the sense of the classical complexity classes of theoretical computer science) on the algorithms we develop. Instead, we will use a simpler notion of efficiency:

We say that a problem can be solved efficiently, if there exists a polynomial time SOS algorithm for it:

Definition 3.27**.**

Suppose $\mathcal{B}$ is a black box that can solve all $n$-variate SOS optimisation programmes of degree $d$ in time $n^{d}$. An SOS algorithm is an algorithm which may make arbitrarily many calls to the black box $\mathcal{B}$ in addition to the usual algorithmic operations (for each call, the time it takes the black box to process the call is of course counted towards the running time of the SOS algorithm).

In the next section, we will develop an SOS algorithm for the moment decomposition problem which calls the black box $m$ times with $d=2m$. This algorithm is not efficient with respect to the notion we just defined. This is due to the fact that the black box calls take time exponential in $m$.

In Section §5, we will then show that basically the same algorithm can be used with a lower value of $d$ (e.g. $d$ so small that we get a polynomial time SOS algorithm), as long as we accept that the algorithm will return only an approximate solution. In practice, the value of $d$ will usually be determined by the amount of data that is given (at least if we are in the empirical case). An increase in the value of $d$ will therefore not just mean an increase in computation time, but also require more data, as outlined in the introductory chapter.181818Note that increasing the accuracy at will is thus usually not possible with this algorithm. However, if we get sufficiently close to the actual solution, then we could try to find the latter by running a local searching procedure afterwards. Therefore it makes sense to analyse what guarantees we get for a fixed value of $d$ instead of going the other way around and say what value of $d$ is needed to get a guarantee that the approximation error is smaller than a fixed constant $\varepsilon$. Our approach is also the one which is technically more convenient.

As a last note on SOS programming, recall that in the introductory chapter we defined pseudo-expectation operators $\mathbb{E}$ of degree $d$ as linear functionals on the space $\mathbb{R}[X]_{\leq d}$ satisfying

[TABLE]

and claimed that we could efficiently optimise over them. This is due to the fact that the cone of pseudo expectation operators of degree $d$ is the dual cone of $\sum\mathbb{R}[X]_{\leq\frac{d}{2}}^{2}$. The literature on SOS programming seems to have a bias towards the dual point of view, but we will prefer to work directly with optimisation over SOS polynomials.

3.4 Multilinear Algebra on the Space of Polynomials

In this chapter, we will introduce some basic notions and notations that will facilitate the operations we have to employ when dealing with (symmetric) tensors.

We will have to redefine the usual methods for tensor manipulation in terms of polynomials. For those already used to tensor decomposition, but unused to see tensors as members of a polynomial ring: We will, in particular, define the multilinear multiplication in such a way that it can be computed in terms of the monomial coefficients (i.e. the upper triagonal) only. It turns out that this was the reason why we had to scale the standard basis of monomials by a factor $\binom{d}{\alpha}$, which corresponds to the size of one of the equivalence classes described in §3.2.

Notation 3.28**.**

For $P\in\mathbb{R}[X]_{\leq d}$ and $k\leq d$, we denote by $P_{=k}$ the $k$-th homogeneous part of $P$, that is

[TABLE]

We have of course $P_{=k}\in\mathbb{R}[X]_{=k}$ and

[TABLE]

The reason why this representation is more suitable for our purposes is that it fits well with the following inner product on $\mathbb{R}[X]_{\leq d}$.

Definition 3.29**.**

(Reznick’s Scalar Product, [Rez92]) Let $P,Q\in\mathbb{R}[X]_{\leq d}$. Consider their representations

[TABLE]

We define

[TABLE]

which is bilinear in $P$ and $Q$. Fixing the left hand side, we get a linear form

[TABLE]

It is immediate that $\langle\cdot\mid\cdot\rangle_{F}$ defines an inner product on $\mathbb{R}[X]_{\leq d}$. In fact, it is actually just a rescaled version of the standard scalar product on $\mathbb{R}[X]_{\leq d}$ w.r.t. the monomial basis. The following proposition shows that the Reznick inner product allows us to express polynomial evaluation.

Proposition 3.30**.**

Let $P\in\mathbb{R}[X]_{=k}$ homogeneous and $a\in\mathbb{R}^{n}$ (where $k\in\mathbb{N}_{0}$). Then

[TABLE]

More general, if $P\in\mathbb{R}[X]_{\leq d}$ is not necessarily homogeneous, then for $k\leq d$

[TABLE]

Proof.

Let first $P\in\mathbb{R}[X]_{=k}$. By the multinomial theorem we have

[TABLE]

Using this representation, we can write down both sides of the equation as:

[TABLE]

The last step in the second line uses that all non-degree-$k$ coefficients of $P$ vanish due to homogeneity. The claim for general $P$ is an immediate consequence. ∎

We’ll need a slightly more general version of this fact:

Proposition 3.31**.**

Let $P\in\mathbb{R}[X]_{\leq d}$ and $a\in\mathbb{R}^{n}$. Then

[TABLE]

Proof.

Applying the previous proposition on all homogeneous parts $P_{=k}$, we get

[TABLE]

∎

Sometimes we want to apply the Reznick scalar product only partially to a given Polynomial $P$. This is achieved by the following “tensoring”.191919For homogeneous $Q$, this corresponds to the notion of the tensor product of linear maps. We will use (or, more accurately, abuse) this notation even for the case when $Q$ is not homogeneous.

Definition 3.32**.**

Fix an integer $k\in\{0,\ldots,d\}$, which we will suppress in the notation along with $d$. Let $Q\in\mathbb{R}[X]_{\leq k}$. Denote by $V=({V}_{1},\ldots,{V}_{n})$ a new vector of algebraic unknowns. Define

[TABLE]

for each $P\in\mathbb{R}[X]_{\leq d}$. Here, the algebraic vector $V$ is plugged in formally as if it was an ordinary real vector $v\in\mathbb{R}^{n}$. This operation yields a homogeneous polynomial in $\mathbb{R}[V]_{=d-k}$. Note that the map

[TABLE]

is linear in both $P$ and $Q$.

The following proposition shows that $\phi_{Q\otimes id}$ can be seen as a linear shrinking map that replaces powers of $\langle{a},X\rangle$ by “evaluation terms” $Q(a)$.

Proposition 3.33**.**

For any $a\in\mathbb{R}^{n}$ and any polynomial $Q$ of degree at most $k\leq d$ we have

[TABLE]

Proof.

From Prop. 3.31 it follows that $\langle Q\langle X,v\rangle^{d-k}\mid\sum_{l=0}^{d}\langle a,X\rangle^{l}\rangle_{F}=Q(a)\langle a,v\rangle^{d-k}$ for any concrete value $v\in\mathbb{R}^{n}$ of $V$. Since $\mathbb{R}^{n}$ is Zariski-open, this identity must hold on the level of variables, too. ∎

Note that all terms $\langle a,X\rangle^{l}$ of degree less than $d-k$ are redundant: Their contribution will vanish, since $Q\langle X,v\rangle^{d-k}$ has no monomials of degree less than $d-k$. After the reduction by the polynomial $Q$ is done, there is no need to keep the additional variable $V$. Therefore, we use the convention to silently replace $V$ again by $X$.

Corollary 3.34**.**

Let $T_{0},\ldots,T_{d}$ be given Tensors of orders $0,\ldots,d$ having a moment decomposition $\lambda_{1}a_{1},\ldots,\lambda_{m}a_{m}$, that is,

[TABLE]

Let $Q\in\mathbb{R}[X]_{\leq d}$ be of degree $k$. Then we have

[TABLE]

Proof.

[TABLE]

∎

4 The v-Algorithm for exact Moment and Tensor Decomposition

4.1 v-Algorithm for Moment Decomposition

In this section, we consider a finitely supported measure $\mu=\sum_{i=1}^{m}\lambda_{i}\delta_{a_{i}}$ and we want to recover the $\lambda_{i}$ and the distinct vectors $a_{i}$. Our main result in this section is that if we are given sufficiently many moments $\mathcal{M}_{0},\ldots,\mathcal{M}_{d}$ of $\mu$ (where $d\geq 2m$ will suffice), we can compute them exactly by SOS programming as a the unique moment decomposition of $\mathcal{M}_{0},\ldots,\mathcal{M}_{d}$. Precisely:

Theorem 4.1**.**

*(Exact decomposition from moments up to $2m$)

Let $m>1$, ${a}_{1},\ldots,{a}_{m}\in\mathbb{R}^{n}$ be distinct vectors, ${\lambda}_{1},\ldots,{\lambda}_{m}\in\mathbb{R}_{>0}$. Suppose we are given the $k-$th moment*

[TABLE]

of the quadrature formula $\lambda_{1}a_{1},\ldots,\lambda_{m}a_{m}$ for each $k\in\{1,\ldots,2m\}$. Then $\lambda_{1}a_{1},\ldots,\lambda_{m}a_{m}$ is the only moment decomposition of ${M}_{1},\ldots,{M}_{2m}$ and there exists an SOS algorithm to compute $\lambda_{1}a_{1},\ldots,\lambda_{m}a_{m}$ from the input ${M}_{1},\ldots,{M}_{2m}$ in time $n^{\theta(m)}$.

We’ve seen in section §3.4, in particular from 3.34, that with any polynomial $W$ of degree $d-2$ we can perform some sort of a “reduction to a matrix”

[TABLE]

The basic step in the following algorithm is to optimise for a nonnegative “weight function” $W$ which will allow us to compute the matrix202020Technically, this is a quadratic form. As we did with tensors before, we aren’t going to distinguish between symmetric matrices and quadratic forms for simplicity.

[TABLE]

We want to find a weight function which concentrates on one of the $a_{i}$: If the weight function (approximately) satisfies W(a_{i})=\begin{Bmatrix}\frac{1}{\lambda_{i}},&\text{if i=j}\\ 0,&\text{otherwise}\end{Bmatrix} for some $j\in[m]$, then $M$ will be (close to) the rank-$1$-matrix $\langle{a_{j}},X\rangle^{2}$ and we will be able to get $a_{j}$ out of $M$. In order to decide which component $a_{j}$ we are going to recover in a single step, we choose $v\in\mathbb{S}^{n-1}$ at random and then solve for the $a_{j}$ which maximises $\langle a_{i},v\rangle$. Let’s put this all together:

Now, this is in some sense an “illicit” formulation of an algorithm, since we wrote down the operations, functions and constraints involved in terms of the solution $\lambda_{1}a_{1},\ldots,\lambda_{m}a_{m}$. To get an actually implementable algorithm (and thus get a proof of Theorem 4.1), we need to show that all of these can be rewritten in terms of the input data ${M}_{1},\ldots,{M}_{2m}$.

Let us first, though, review why this algorithm produces the correct result. To this end, we formulate correctness in the following theorem and prove it.

Theorem 4.2**.**

Let $W^{\ast}$ be an optimiser of

[TABLE]

over all SOS polynomials $W\in\mathbb{R}[X]_{\leq 2m}$ satisfying

[TABLE]

Then for $j=\operatorname{argmax}_{i\in[m]}\langle a_{i},v\rangle$:

[TABLE]

and therefore $\sum_{i=1}^{m}\lambda_{i}W^{\ast}(a_{i})a_{i}a_{i}^{T}=a_{j}a_{j}^{T}$. From that we can compute $\pm a_{j}$ and $\lambda_{j}$.

Proof.

Note that since we are optimizing over SOS polynomials of degree $d-2\geq 2m-2$, the feasible space contains the interpolation polynomials

[TABLE]

which satisfy

[TABLE]

Our claim is that any optimal solution $W^{\ast}$ must satisfy the property in (23) as well for $j=\operatorname{argmax}_{i\in[m]}\langle a_{i},v\rangle$. (Note that with probability one over the choice of $v$, the values $\langle a_{i},v\rangle$ will all be distinct. Hence there will be only one maximiser $j$ which justifies the usage of the $\operatorname{argmax}$ function.) Indeed, we see that for any feasible $W$, since $W$ is a sum of squares, we have $\lambda_{i}W(a_{i})\geq 0$ and thus

[TABLE]

where strict inequality holds if the term $\lambda_{j}W(a_{j})$ is strictly smaller than $1$.212121Here we use that the values $\langle a_{i},v\rangle$ are all distinct, whence there is precisely one maximum value. Since $I_{j}$ attains the upper bound (24), $W^{\ast}$ must, too, and by the above we get $\lambda_{i}W^{\ast}(a_{i})=\delta_{ij}$ as claimed. Hence the matrix $M$ is precisely

[TABLE]

Any eigenvector of $M$ with nonzero eigenvalue is contained in $\operatorname{im}(M)$ and therefore a multiple of $a_{j}$. Since $u$ has unit length, $u=\pm\frac{a_{i}}{\|a_{i}\|}$. It remains to recover $\|a_{j}\|$ and the sign. Note that for the corresponding eigenvalue $\mu$ we know

[TABLE]

Hence $\mu=\|a_{j}\|^{2}$ and thus $\sqrt{\mu}u$ will be either $a_{j}$ or $-a_{j}$. It remains to show that the sign is correct. Note that

[TABLE]

will also be equal to $a_{j}$. Hence $\operatorname{sgn}(\langle u,w\rangle)=1$ if and only if $\sqrt{\mu}u=a_{j}$ and $\operatorname{sgn}(\langle u,w\rangle)=\operatorname{sgn}(-\|a_{j}\|^{2})=-1$ otherwise. Now since we recovered $a_{j}=\sqrt{\mu}u$, we might plug it into the identity $\lambda_{j}W^{\ast}(a_{j})=1$ to recover $\lambda_{j}=\frac{1}{W^{\ast}(\sqrt{\mu}u)}$ as well. ∎

Remark 4.3**.**

At this point the observant reader is probably questioning why we didn’t just output $w$, which would have been much simpler. There is no actual reason for this except for the fact that reducing to a matrix will be slightly more natural for the case of even degree tensor decomposition in the next section: When we can not distinguish between $\pm a_{i}$ anyway, computing the matrix $a_{i}a_{i}^{T}$ will save us from having to make a sign-decision. This decision is instead transferred to the computation of the top eigenvectors.222222By that we mean the eigenvectors $u\in\mathbb{S}^{n-1}$ of unit length that correspond to the largest eigenvalue. For some matrices, e.g. for (real) symmetric matrices with $n$ distinct eigenvalues, there are precisely two top eigenvectors $\pm u$.

Remark 4.4**.**

*On invalid input (that is, input which does not fulfill the requirements), weird things might happen: If for instance the degree of $W$ is too low, then $M$ will likely not be a rank 1 matrix (but it might still be very close to one, which is essentially what we show in §5). And if for instance $\mu=\sum_{i=1}^{m}\lambda_{i}\delta_{a_{i}}$ is a signed measure (where the weights $\lambda_{i}\in\mathbb{R}\setminus\{0\}$ may be negative) then we can show by the same interpolation argument that the decomposition is unique for $d\geq 2m$ (assuming the $m$ components are distinct). However, we can in general not compute this decomposition by the $v$-Algorithm: If one of the components has a negative weight $\lambda_{i}<0$, then the optimisation problem will immediately become unbounded and thus there will be no optimiser. This result will then also certify that $\mu$ is not a measure of finite support size (for otherwise we would have a bound on the optimal value).

Still, this property of the algorithm reveals an issue problematic for noise stability: If the moments of $\mu$ are disturbed by some small noise tensors, this might result in an unbounded problem, since e.g. an arbitrarily small but negative $0\approx\lambda_{i}<0$ might have an arbitrarily high impact on the optimal value of the SOS programme. This can be dealt with by introducing a “complexity bound” on $W$ implying e.g. an upper bound on $\|W\|_{\textrm{F}}$. Bounding the complexity of $W$ comes with the usual issues known from Machine Learning problems: Having the bound too small results in a sub-optimal solution (“underfitting”), while leaving it too high makes the algorithm more vulnerable to noise (“overfitting”).*

Let us now summarise what is left to do. We need to show that…

1. (a)

The algorithm can be written down in a way that relies only on the input data. 2. (b)

We can repeat the algorithm $m$ times in order to recover all components. In particular, we need a constraint which ensures that we will not recover the same component twice.

We’ll start by writing down the algorithm in rigid manner, showing (a) and (b) simultaneously (see Algorithm 3). Now we’re ready to prove Theorem 4.1, in particular that Algorithm 3 produces the correct results.

Proof.

(of Theorem 4.1) Let us show first that the functions and constraints are indeed the same as before. Suppose the $T_{k}$ admit a simultaneous decomposition

[TABLE]

Let us fix this decomposition. By the results of §3.4, in particular by Prop. 3.31 and Cor. 3.34, we have that

[TABLE]

Note that all of these (the objective function, the left side of the constraint and $M$) are linear functions in the variable $W$ and $W$ is the only variable, since $v$ is fixed and the $T_{k}$ are all given constant tensors. In particular, $(\ast)$ is indeed an SOS programme.

The equations (26) show that the first loop of the algorithm will indeed yield a component of the $T_{k}$. In fact, it’ll yield the component $a_{i}$ which maximises $\langle a_{i},v\rangle$. Note that the constraint $W(c_{i})=\sum_{|\alpha|\leq d}W_{\alpha}c_{i}^{\alpha}=0$ is linear in $W$, too (while $c_{i}$ is known at that time), and can thus be expressed by computing the powers $c_{i}^{\alpha}$ of $c_{i}$.

To see that we can recover all components, we claim that:

1. (i)

in each step we get a component that we didn’t have before. 2. (ii)

the algorithm terminates after precisely $m$ steps.

(i) is actually quite clear: In the first step, since there are no restrictions of the kind $W(c_{i})=0$ yet, $W^{\ast}$ will concentrate on the component maximizing $\langle a_{i},v\rangle$ over all $i\in[m]$ as we’ve seen before. Suppose that in the first $k$ steps we’ve got each time a distinct vector with corresponding weight, say $\rho_{1}c_{1},\ldots,\rho_{k}c_{k}$. By induction suppose that we have a matching:

[TABLE]

for some permutation $\sigma$ of $[m]$ (and likewise for the weights). Let $c_{k+1}$ denote the output of the k+1-st iteration which was subject to the constraints $W(c_{1})=\ldots=W(c_{k})=0$. Now some of the terms in the estimation (24) will vanish, allowing us to refine our bound to:

[TABLE]

which is again attained by the interpolation polynomial corresponding to the maximising $a_{i}\notin\{{c}_{1},\ldots,{c}_{k}\}$. By the same arguments that were used in the proof of Theorem 4.2 we can then again reduce our tensor to a rank $1$ matrix and see that we get a new component out of it (which is precisely the maximising $a_{i}\notin\{{c}_{1},\ldots,{c}_{k}\}$ and equal to $c_{k+1}$) together with the corresponding weight.

(ii): After $m$ iterations, by our previous considerations, we will have added the constraints $W(a_{1})=\ldots=W(a_{m})=0$. But then the constraint $\sum_{i=1}^{m}\lambda_{i}W(a_{i})=1$ of the SOS problem becomes infeasible. ∎

Remark 4.5**.**

Note that when choosing the decomposition ${a}_{1},\ldots,{a}_{m}$, it might seem like we had a freedom of choice. However, we get that the resulting output $\{c_{1},\ldots,c_{m}\}$ is equal to $\{{a}_{1},\ldots,{a}_{m}\}$ independently of the choice of the $a_{i}$. This shows that there is in fact only one possible choice for the $a_{i}$ (up to renumbering). Thus the algorithm provides us with a proof of uniqueness in this setting.

Remark 4.6**.**

Note that the algorithm would also work without resampling the random vector $v$ in every iteration (as long as $v$ is not orthogonal to any of the components which is the case with probability 1). However, then the optimal value might get very low on the last iterations which is bad for several reasons that will become slightly more obvious from section §5. Also, if one of the components is the zero vector, it will be recovered if all other remaining components attain a negative value of $\langle a_{i},v\rangle$. In this case, $M$ will be the zero matrix.

4.2 v-Algorithm for Tensor Decomposition

The algorithm developed in §4.1 can be adapted to work for tensor decomposition as well, that is, when only one tensor $T$ of degree $d\in 2\mathbb{N}_{>2}$ (and rank $m\leq d/2$) is given. This is possible even though the $v$-algorithm makes explicit use of all the lower-degree-moments to write down the constraints and the objective function.

To achieve this, we would need to generate something alike lower degree moments of the $a_{i}$. But generating the actual lower degree moments from the input $T$ has to be impossible, particularly since the $a_{i}$ will never be unique in the setting of even degree tensor decomposition. But we have a chance in generating some “fake moments” by abusing the input ambiguity which we described in 3.12: Suppose again for this section that

[TABLE]

is a tensor decomposition of $T$. Pick some random $w\in\mathbb{S}^{n-1}$ and compute the tensors

[TABLE]

by applying the linear maps $\phi_{\langle{w},X\rangle^{d-k}\otimes\operatorname{id}}$ to $T$. Note that $T=T_{d}$.

Now, the $T_{k}$ aren’t necessarily proper moments of the $a_{i}$. But in some sense, they are fake moments: Indeed, write $\lambda_{i}\mathrel{\mathop{\mathchar 58\relax}}=\langle a_{i},w\rangle^{d}$ and $b_{i}\mathrel{\mathop{\mathchar 58\relax}}=\frac{a_{i}}{\langle a_{i},w\rangle}$. Then we have

[TABLE]

From that we see that the tensors $T_{k}$ admit a simultaneous weighted decomposition with weights $\lambda_{i}$ and components $b_{i}$. Once we have both the weights $\lambda_{i}$ and the components $b_{i}$, we might get the true components $a_{i}$ out of them. Let’s see what will happen if we plug those fake moments into the $v$–algorithm.

Algorithm 4 shows a version of the $v$–algorithm that has been “hacked” in such a way that whenever it gets fed with the fake moments described above it outputs the correct components up to a factor of $\pm 1$. We do not have to prove much here since most of the work has already been done in the previous section. In particular, we know that $(\lambda_{i},b_{i})$ in line 8 will be one component-weight pair of the unique solution to the moment decomposition problem

[TABLE]

that we described above. Since we know uniqueness, there is no other possibility than

[TABLE]

and $b_{i}=\frac{a_{j}}{\langle a_{j},v\rangle}$ for some $a_{j}$. Furthermore, this $a_{j}$ will be exactly the one maximising the ratio $\frac{\langle a_{j},v\rangle}{\langle a_{j},w\rangle}$, since the objective function satisfies

[TABLE]

for $j=\operatorname{argmax}_{i\in[m]}\frac{\langle a_{j},v\rangle}{\langle a_{j},w\rangle}$. Note that we used again that the values $\frac{\langle a_{j},v\rangle}{\langle a_{j},w\rangle}$ are distinct and well-defined due to randomness of $w$. Actually, $v$ does not necessarily need to be random any more, since $w$ already is: The only thing we have to ensure is that $v\neq w$. (However, low correlation between $v$ and $w$ is of course preferable, e.g. $v\perp w$).

Remark 4.7**.**

(Success probability in floating point systems) Note that the $v$–algorithm has success probability $1$ even though it contains a randomised step which is the choice of $v$. Therefore we can say that it is just a “mildly randomised” algorithm contrary to the procedures in [BKS15] and [GM15], where even a single component recovery step can have a very low success probability $\ll 1$. However, the choice of $v$ will start to matter in the approximate setting of §5. Though even here there is an issue when performing this algorithm in a finite-accuracy floating point system: We have used the argument that the values of the scalar products $\langle a_{i},v\rangle$ are all pairwise distinct. Now with a certain positive probability it could happen that these values are indistinguishable up to the numerical accuracy. In this case, our recovery step could fail and we would have to repeat it. Now what is the probability for this event? Let us consider the set

[TABLE]

which corresponds to the good event that two values of the scalar products are distinguishable up to $\epsilon$. Geometrically spoken, this set can be imagined as the intersection of the two half spaces $\{v\in\mathbb{R}^{n}\mid\pm\langle a_{i}-a_{j},v\rangle>\epsilon\}$ with the $n$– dimensional unit sphere. It can be shown that the probability to hit $\mathcal{E}_{ij}$ with the choice of $v$ is $\frac{2}{\pi}\arccos(\frac{\epsilon}{\|a_{i}-a_{j}\|})$ and that $\arccos(x)\to\leavevmode\nobreak\ \frac{\pi}{2}\quad(x\to\leavevmode\nobreak\ 0)$ with convergence order $\theta(x)$. The probability that not all the values are pairwise distinct is then by the union bound \leq m\big{(}1-\frac{2}{\pi}\arccos(\frac{\epsilon}{\min_{i\neq j}\|a_{i}-a_{j}\|})\big{)}. Thus recovery works if $\epsilon$ is several orders of magnitude smaller than $\|a_{i}-a_{j}\|$ and $m$ is reasonable.

4.3 Qualitative Aspects

This short section is dedicated to the discussion of special cases as well as some ideas, optimisations and alternative approaches that we will not elaborate in full detail. Let us start with some special situations where exact recovery is possible from less data than $d=2m$. In particular, we want to point out the connection between the $v$–algorithm and Jennrich’s algorithm.

Orthogonal Components and Connection to Jennrich’s Algorithm

The above algorithm implies that there is only one quadrature formula fulfilling the requirements and thus gives a proof of uniqueness in this setting. However, note that knowing $d\geq 2m$ moments is quite a heavy requirement.

In special instances, the $v$-algorithm might achieve exactness for way smaller values of $d$. As a friendly example, consider the tensor

[TABLE]

where $e_{i}\in\mathbb{R}^{n}$ denotes the $i-$th standard basis vector. On the set $\{e_{1},\ldots,e_{n}\}$, the polynomial $X_{i}$ evaluates the same way as the $0/1$ indicator function $\mathbb{1}_{e_{i}}$. Well, $X_{i}$ is no SOS polynomial, but $X_{i}^{2}$ is and from that we get that the $v-$algorithm will recover the exact solution even for $d=4$. A more detailed analysis shows that the case $d=3$ can also be made to work: The objective function

  

[TABLE]

can of course be written down with a degree-2+1 polynomial $W\langle{v},X\rangle$ on the left hand side of Reznick’s inner product. However, since $W$ consumes two degrees, we will not have enough data to compute the matrix $M$. But we can still do the “vector reduction” outlined in Equation (25). In the case $a_{i}=e_{i}$, we may now use $W=X_{i}^{2}$ as a testimony to show that this adapted $v-$algorithm can recover the components of $T$ exactly even for $d=3$. More generally, if the $a_{i}$ are linearly independent, then there exists an invertible matrix $S$ such that $Sa_{i}=e_{i}$. Hence $\tilde{W}\mathrel{\mathop{\mathchar 58\relax}}=W(SX)$ testimonies that the $v$-algorithm succeeds for the linearly independent case.

Its theoretical properties in this setting correspond precisely to those of Jennrich’s algorithm, where we needed the moments of degree $2$ and $3$ (the second moment was used to transform the components to orthogonal ones). These moments are precisely the ones that the $v-$algorithm needs in this setting: Note that our testimony $W$ is homogeneous of degree $2$, whence we do not need the moments of degree [math] and $1$. Also, it seems that the lower-order moment generation procedure isn’t able to generate the second order moment: This is due to the fact that the weights $\lambda_{i}=\langle a_{i},v\rangle^{d}$ have to be positive for which we get a guarantee only if $d$ is even. This analysis demonstrates that the $v$-algorithm can indeed be seen as a generalisation of Jennrich’s algorithm to the overcomplete setting.

Using Symmetry

Another friendly example are the vertices of a regular simplex with edge length $r$ in $n$ dimensions. Here $\|a_{j}-a_{i}\|$ attains only two values: [math] and $r$. Therefore we can choose a univariate SOS interpolation polynomial232323E.g. $p=(1-\frac{1}{r}\Lambda)^{2}$ $p\in\mathbb{R}[\Lambda]$ of degree $2$ with $p(0)=1,p(r)=0$ and we see that $p(\|X\leavevmode\nobreak\ -\leavevmode\nobreak\ a_{j}\|^{2})$ testimonies that exact recovery is possible already from the moments of degree less or equal $4$, independent of $n$. It’s one of the beautiful properties of this algorithm that it can exploit such symmetry with seemingly no additional effort.242424Note that tensor decomposition over the simplex is relatively easy: Since there are only $n+1$ points ${a}_{1},\ldots,{a}_{n+1}$ in question, we can make the ansatz $T_{k}=\sum_{i=1}^{m}\lambda_{i}\langle{a_{i}},X\rangle^{k}$ and just search for the weights via linear least squares. This shall thus just be seen as an illustrative example. Note further that the latter algorithm needs a priori knowledge on where the components are situated whereas the $v$-algorithm does not.

Components with Differing Lengths and Discriminatory Polynomials

An obvious downside that comes with the randomness is that we will not get the components in any particular order. Depending on the choice of $v$, they could get recovered in any order, with the odds slightly favouring components of high norm $\|a_{i}\|$ to come first.252525Note that the actual chance for a component to become $\operatorname{argmax}_{i}\langle a_{i},v\rangle$ depends not only on the norm, but also on the geometry of the $a_{i}$: Isolated components are more likely to be picked than components which have a “neighbour” being highly correlated with them.

However, there are deterministic variants of the algorithm working for certain special cases, of which we want to present one. Suppose the $a_{i}$ all have different lengths. Now, what we could try to do in the case of moment decomposition is to replace the objective function $(\ast)$ by

[TABLE]

Then we can again characterise the optimizing weight functions like we did before: Indeed, they will concentrate on the highest-norm component. Be careful: The requirement of distinct norms $\|a_{i}\|$ is inherently necessary. Otherwise in general we will not get a rank $1$ matrix! (Note further that this will recover the highest norm components first and not necessarily the highest magnitude components).

In the general setting of tensor decomposition, it’s possible (and it may or may not be useful) to replace $\langle b_{i},v\rangle$ by $\|b_{i}\|^{2}=\frac{\|a_{i}\|^{2}}{\langle a_{i},w\rangle^{2}}$, since sufficient randomness is already given by the denominator of $\frac{\|a_{i}\|^{2}}{\langle a_{i},w\rangle^{2}}$. More generally, we can actually plug in any non-constant polynomial $f$ of degree at most $2$ in place of $\langle v,X\rangle^{2}$, since the differences $f(a_{j})-f(a_{i})$ will be non-zero with probability one (for moment decompositions this is also true as long as the polynomial $f$ is sampled from a continuous distribution. If we have some prior knowledge about the domain of the $a_{i}$, it might be possible to choose $f$ deterministically in a way that guarantees that the differences are always non-zero).

The analysis suggests that it makes sense to introduce the notion of discriminatory polynomials, which are polynomials achieving high values of $\min_{i\neq j}f(a_{j})-f(a_{i})$ for the component $a_{j}$ maximising $f$. Once we move to the setting of approximate recovery, the choice of $f$ becomes crucial. We will see that we need to make a choice which provides us with a sufficiently big gap $\min_{i\neq j}f(a_{j})-f(a_{i})$. While it can be shown that random choices achieve this gap with substantial probability, finding good choices of $f$ could be an interesting subject of further research.

Canonical Procedures for Lower-Degree Moment Generation

In the setting of tensor decomposition, where we do not have the lower degree moments, a common constraint is to require that all components have unit length, i.e. $\|a_{i}\|=1$. The lengths are then either put into some weights $\lambda_{i}$, which yields an equivalent formulation, or the weights are still assumed to be one, which results in a proper special case and a slightly simplified problem. (For some problems related to tensor decomposition this is justified, see for example [BKS15]). This implies that for each even $k$ we can generate

[TABLE]

which is then the true moment of the $a_{i}$. In the general case, we can still compute the even degree fake moments

[TABLE]

and feed them as input to the algorithm. This will force a decision on the input ambiguity 3.12: It will automatically ensure that the algorithm interprets the components as points of unit length $\frac{a_{i}}{\|a_{i}\|}$ and the weights as $\lambda_{i}=\|a_{i}\|^{d}$.

Unfortunately, it’s not clear what we can do in order to obtain the corresponding $\|a_{i}\|^{k}$-scaled moments when $k$ is odd. It’s very plausible though that there is no canonical choice here to make. Note that having these fake moments for odd degree would automatically enable us to distinguish between $a_{i}$ and $-a_{i}$. But these are indistinguishable in an even degree tensor! Here, the moment generation procedure essentially makes a random, noncanonical decision. One could ask though whether it’s possible to generate

[TABLE]

for some $s_{i}\in\{\pm 1\}$, which is a question we have to leave open. This corresponds to the question whether we could potentially plug in $\sqrt{X_{1}^{2}+\ldots+X_{n}^{2}}^{k}$ into (33) in a way that makes sense.

Of course we could still try to run the algorithm with only the even degree inputs. As we’ve seen in section 3.4, this means that we’d essentially optimise over weight polynomials with vanishing odd degree homogeneous parts. If in that way we get a solution vector $c$, it might be an inferior approximation compared to the vector we’d have gotten out of the $v$–algorithm.

However, note that there exist reasonable testimonies for this case, e.g. we could concatenate some univariate interpolation polynomial $p$ with $\|a_{j}\|^{2}\|X\|^{2}-\langle{a_{j}},X\rangle^{2}$: All monomials occurring in $p(\|a_{j}\|^{2}\|X\|^{2}-\langle{a_{j}},X\rangle^{2})$ have even degree. The values of $\|a_{j}\|^{2}\|a_{i}\|^{2}-\langle a_{j},a_{i}\rangle^{2}$ for $i\neq j$ will be distinct from [math] by Cauchy-Schwarz (as long as no other component $a_{i}$ is a multiple of $a_{j}$) and therefore we can choose an interpolation polynomial $p$ such that $p(\|a_{j}\|^{2}\|a_{i}\|^{2}-\langle a_{j},a_{i}\rangle^{2})=\delta_{ij}$. This shows that the approach of this paragraph is a valid alternative to the generation of randomised lower degree fake moments in §4.2.

5 Efficient Tensor Decomposition via SOS

Throughout this section, let again

[TABLE]

a weighted decomposition of the tensors $T_{0},{T}_{1},\ldots,{T}_{d}$. Let also $v\in\mathbb{S}^{n-1}$ be random and reconsider the SOS optimisation problem

[TABLE]

We’ve seen that in the setting $d\geq 2m$ the algorithm allows for a very clean analysis. However, having access to moments of that high order is a very heavy requirement and results in very bad running time: Solving the degree $d$ SOS programmes involved will take our black box $n^{2m}$ time steps, which is exponential in $m$. The obvious next question to ask is what guarantees we will get if we restrict to SOS programmes which are polynomial–sized in both $m$ and $n$. (This would be the case if for instance $d$ depended logarithmically on $m$). In this regime, we can’t of course hope for $W$ to vanish on all but one $a_{i}$. However, we can try to search for a weight polynomial which attains very small values on most of the components while being of relatively high magnitude on one particular component.

The matrix $M$ will then be only very close to a rank $1$ matrix. But as long as this approximation is good enough, we will be able to show that the top eigenvectors of $M$ are still very close to one of the components. To this end, we will make use of the following lemma:

Lemma 5.1**.**

*(Eigenvectors of gapped matrices - Adapted from Lemma A.3 in [HSS16])

Let $M\in\mathbb{R}^{n\times n}$ be a symmetric matrix and let $a\in\mathbb{R}^{n}$ , $\gamma>0$ such that*

[TABLE]

Then for each top eigenvector $u$ of $M$ we have $\langle u,\frac{a}{\|a\|}\rangle^{2}\geq\gamma$.

Proof.

Denote by $\mu_{\max}$ the eigenvalue of $M$ of largest absolute value. Recall the characterisation

[TABLE]

for any top eigenvector $u$ of $M$. Thus $\|M\|_{\operatorname{spec}}=\sup_{x\in\mathbb{S}^{n-1}}|\langle x,Mx\rangle|$. Let $u$ denote a top eigenvector of $M$. In particular, $\|M\|_{\operatorname{spec}}=|\langle u,Mu\rangle|$ and therefore

[TABLE]

Hence $\|M\|_{\operatorname{spec}}-\langle u,a\rangle^{2}\leq\|M\|_{\operatorname{spec}}-\gamma\|a\|^{2}$, implying $\gamma\|a\|^{2}\leq\langle u,a\rangle^{2}$. ∎

This will give us a guarantee that the top eigenvector is close to the true component $a_{j}$ as long as the $W$-weighted Matrix $M$ is sufficiently close to a rank $1$ matrix. Now we are interested in the question how close the matrix $M$ can actually become to a rank 1 matrix for a given value of $d$. From the preceding discussion it should be very clear that there are two essential issues to answer here:

How close can a degree $d\leavevmode\nobreak\ <\leavevmode\nobreak\ 2m$ polynomial come to one of the 0/1 Kronecker delta functions $\delta_{ij}$ on the set $K\leavevmode\nobreak\ \mathrel{\mathop{\mathchar 58\relax}}=\leavevmode\nobreak\ \{a_{1},\ldots,a_{m}\}$? 2. 2.

How close will the optimiser $W^{\ast}$ come to one of the $\delta_{ij}$?

The first question asks for the best testimony in the feasible space, while the second question asks how reliably the optimisation solves for a polynomial which has similar properties as the testimony. In the exact case, recall that the optimiser and the testimony (which was an interpolation polynomial) were equally good for our purposes. In the approximate setup, this will not be the case anymore: Our results will depend on the choice of $v$, so we have to make a good choice. Let us first address though the testimony problem:

5.1 Finding a Testimony

Univariate Chebyshev Polynomials

During this paragraph, let us switch to the univariate case, where $\Lambda$ is a scalar algebraic unknown. In this case, there is already a lot of common knowledge about similar optimisation problems that will help us.

Definition 5.2**.**

The degree-$d$ polynomials $\mathcal{T}_{d}$ uniquely determined by

[TABLE]

are well-known as the Chebyshev polynomials.

These polynomials have several characterising extremal properties one of which is shown in the following lemma.

Lemma 5.3**.**

Let $x\in\mathbb{R}$ not in the interval $[-1,1]$. Consider the optimisation problem

[TABLE]

The unique optimiser to this problem is the Chebyshev polynomial $\mathcal{T}_{d}$.

Proof.

(taken from [FP10]) For any $d+1$ values ${t}_{1},\ldots,{t}_{d+1}$ we can express any $p\in\mathbb{R}[X]_{\leq d}$ by the interpolation formula

[TABLE]

Evaluating this formula in $x$ and using that $|p(t_{i})|\leq 1$ yields the bound

[TABLE]

This bound is attained if and only if $p(t_{i})=\operatorname{sgn}(u_{i})$ for all $i\in[d+1]$ where $u_{i}\mathrel{\mathop{\mathchar 58\relax}}=\prod_{i\neq j}\frac{|x-t_{i}|}{|t_{i}-t_{j}|}$. Since $x\notin[-1,1]$, the terms $|x-t_{i}|$ have the same sign for all $i$. Together with the fact that $t_{1}<\ldots<t_{d+1}$ were sorted in ascending order, we get $\operatorname{sgn}(u_{i+1})=(-1)\operatorname{sgn}(u_{i})$, that is, the signs of $p$ have to alternate between $1$ and $-1$ on $d+1$ values. For a polynomial of degree $d$ bounded by $1$, this is only possible if $t_{2},\ldots,t_{d}$ are the $d-1$ local maximisers of $p$ and $t_{1},t_{d+1}\in\{\pm 1\}$ are boundary maximisers of $p$ on $[-1,1]$.

All $y\in[-1,1]$ may be represented as $y=\cos(\vartheta)$ for some $\vartheta\in[0,\pi)$. Using the identity

[TABLE]

and the properties of the cosine function, we see that over the course of one semi-period of $\cos(\vartheta)$ (that is, as $\vartheta$ increases from [math] to $\pi$) the values of $\mathcal{T}_{d}$ perform $d$ semi-periods, attaining $d+1$ times a value of $\pm 1$. Hence the $d$–th Chebyshev polynomial has indeed this oscillating property and is thus the solution to this optimisation problem. ∎

Definition 5.4**.**

Let $\mathcal{I}=[a,b]$ an interval (where $a<b$). Consider the linear polynomial

[TABLE]

mapping $[a,b]$ to $[-1,1]$ monotonically. Then

[TABLE]

is called the $d$-th Chebyshev polynomial on $\mathcal{I}$.

Corollary 5.5**.**

On the interval $\mathcal{I}=[a,b]$, where $b>a>0$, the Chebyshev polynomial $\mathcal{T}_{\mathcal{I},d}$ on $\mathcal{I}$ is the optimiser of

[TABLE]

The optimum value is

[TABLE]

where we used the abbreviation $\kappa\mathrel{\mathop{\mathchar 58\relax}}=\frac{b+a}{b-a}$. The optimum value grows exponentially with $d$.

Proof.

The optimality property follows trivially from Lemma 5.3. For the explicit value of the optimum we used the well-known fact that

[TABLE]

for $|x|>1$ and that $\kappa=\frac{b+a}{b-a}>1$ since $0<a<b$. ∎

Multivariate Chebyshev Polynomials

Now, from univariate Chebyshev polynomials it’s just a small step to multivariate but rotation symmetric polynomials:

Proposition 5.6**.**

Let $d\equiv 2\mod 4$, $x\in\mathbb{R}^{n}$ and $q\mathrel{\mathop{\mathchar 58\relax}}=\|X-x\|^{2}$. We have that

[TABLE]

is a rotation symmetric262626that is, w.r.t rotations around $x$ SOS polynomial.

Proof.

Since $d\equiv 2\mod 4$ we know that $\mathcal{T}_{\mathcal{I},(d-2)/2}$ is an even degree polynomial. At any of its local minima it attains the value $-1$. Therefore, $\mathcal{T}_{\mathcal{I},(d-2)/2}+1$ is a nonnegative polynomial in one variable – hence a sum of squares. The concatenation with another SOS polynomial $q$ yields of course again an SOS polynomial. Symmetry under rotations around $x$ is clear. ∎

Note that this polynomial has a local maximum at $X=x$. It attains high values on $y$ when $\|y-x\|$ is either very small or very big. For moderately sized $\|y-x\|^{2}\in[a,b]=\mathcal{I}$ it will attain low values. Therefore, such a polynomial does just what we want if $x=a_{j}$ and $\|a_{j}-a_{i}\|\in[a,b]$ for all $i\neq j$. We could thus take

[TABLE]

(with $\kappa_{\min},\kappa_{\max}$ as defined in 3.18), but it might be wise to be cautious and leave some extra space, e.g. to take

[TABLE]

for some small $\delta\geq 0$.272727We will need this $\delta$-wide safety margin in section §5.5. The factor of 2 is needed due to an unfortunate suboptimality of the SOS proof we will give there. Let us therefor write

[TABLE]

These SOS polynomials aren’t all members of the feasible space, since in general they will not fulfill

[TABLE]

Of course this can be fixed by a simple rescaling, but we do not know the exact scaling factor. Let $W_{j,\delta,\text{test}}$ denote the correctly rescaled version of $W_{j,\delta}$. What we’re actually interested in is

[TABLE]

The following proposition deals with the unknown scaling factor by estimating $\mathfrak{w}$. Let us in the following suppress $j$ and $\delta$ (and sometimes also $\mathcal{I}$) in the notation.

Proposition 5.7**.**

Let $d\equiv 2\mod 4$. Fix some $j\in[m]$. For the polynomials $q\mathrel{\mathop{\mathchar 58\relax}}=\|X-a_{j}\|^{2}$ and $\mathcal{C}_{d}\mathrel{\mathop{\mathchar 58\relax}}=\mathcal{T}_{\mathcal{I},(d-2)/2}+1$ we have that

[TABLE]

is a rotation symmetric282828that is, w.r.t rotations around $a_{j}$ member of the feasible space of $\eqref{eq:sosprog}$ after rescaling. Let us denote by $W_{\text{test}}$ the unique multiple of $\mathcal{C}_{d}(q)$ such that

[TABLE]

and write

[TABLE]

We claim that

[TABLE]

where $\kappa=\frac{2\kappa_{\max}+\frac{1}{2}\kappa_{\min}}{2\kappa_{\max}-\frac{1}{2}\kappa_{\min}+2\delta}$

Proof.

Note that the ratio $\mathfrak{r}=\frac{\mathcal{C}_{d}(0)}{\|\mathcal{C}_{d}\|_{\infty,\mathcal{I}}}$ is invariant under rescaling, whence we clearly have $\mathfrak{r}=\frac{W_{\text{test}}(a_{j})}{\mathfrak{w}}$. Unlike $\mathfrak{w}$ though, we can explicitly compute $\mathfrak{r}$. For the estimate, note that

[TABLE]

This shows the first inequality in (48). The second inequality follows by plugging in the known explicit value of $\mathfrak{r}$ or $\mathcal{C}_{d}(0)$, respectively. ∎

We will usually not work with the monstrosity that is Equation (48), but use instead the qualitative fact that $\mathfrak{w}\leq\frac{1}{\mathfrak{r}\lambda_{j}}$ drops exponentially with $d$.

5.2 Guarantees for the Optimiser

Now that we know that there is a testimony in the feasible space which has the properties we want, we still have to show that the optimiser $W^{\ast}$ will concentrate on precisely one of the $a_{i}$ (like the testimony did). Alas, this problem becomes much more difficult in the approximate setting: When we had the perfect testimony satisfying (23) in the feasible space, we could easily deduce (22). This is due to the fact that if $\mathcal{I}_{j}$ does not waste any concentration on a sub-optimal component, then neither can $W^{\ast}$. Now however, our testimony will likely waste some concentration on components with a very low value of $\langle a_{i},v\rangle$ and the optimiser could theoretically beat the testimony by putting e.g. all of that concentration on the component achieving the second largest value of $\langle a_{i},v\rangle$, provided that $\Delta_{j,i,v}\mathrel{\mathop{\mathchar 58\relax}}=\langle a_{j},v\rangle-\langle a_{i},v\rangle$ is very small.

The latter event is unlikely to happen, if $\|a_{j}\|\geq\|a_{i}\|$ and the components are not too correlated, but the choice of $v$, whether it’s done randomly or deterministically, is still a non-trivial issue which we will address in its own section. For the moment, let us just see what guarantees we get for a fixed value of $\Delta_{j,i,v}$. To be slightly more general, let us further replace the discriminatory polynomial $\langle{v},X\rangle$ by an arbitrary polynomial $f$ of degree at most two. Likewise, $\Delta_{j,i,v}$ is replaced by $\Delta_{j,i,f}$.

Lemma 5.8**.**

Let $W^{\ast}$ denote the optimiser of

[TABLE]

where $f$ is some polynomial in $X$ of degree less or equal $2$ (e.g. $f=\langle{v},X\rangle$ or $f=\langle{v},X\rangle^{2}$ or $f=(X_{1}^{2}+\ldots+X_{n}^{2})$). Then

[TABLE]

Proof.

Let us write again for brevity

[TABLE]

where we suppose again that $a_{j}$ is the maximiser of $f$ among all $i$. Since $W^{\ast}$ is the optimiser of $(\ast)$, we have

[TABLE]

If $W_{\text{test}}(a_{i})\geq W^{\ast}(a_{i})$ for some $i\neq j$, then $W_{\text{test}}(a_{i})-W^{\ast}(a_{i})\leq\mathfrak{w}$. Otherwise we use $W_{\text{test}}(a_{i})\leavevmode\nobreak\ -\leavevmode\nobreak\ W^{\ast}(a_{i})\leavevmode\nobreak\ \leq\leavevmode\nobreak\ \mathfrak{w}\leavevmode\nobreak\ -\leavevmode\nobreak\ W^{\ast}(a_{i})$. Combining these, we can upper-bound the left hand side and get

[TABLE]

Hence

[TABLE]

or, equivalently

[TABLE]

∎

Note that for $f=\langle{v},X\rangle^{2}$ or $f=\langle{v},X\rangle$, the maximum delta value can easily be bounded by $\max_{i\neq j}\|a_{i}-a_{j}\|\|a_{i}+a_{j}\|$ or $\max_{i\neq j}\|a_{i}-a_{j}\|$, respectively. The problem is thus to find a polynomial $f$ or a choice of $v$ such that the minimum delta value can be lower bounded. This seems to be a very hard task to make deterministically, but the chance of getting a good choice of $v$ at random is not too bad, as we will see. In certain special cases, there might be systematic choices, especially when prior knowledge to the domain of the $a_{i}$ is given (e.g. if they lie on the $\{\pm 1\}^{n}$-hypercube).

5.3 Approximate Recovery of One Component

Now we’ve got everything we need in order to formulate the recovery step for a single component, supposed we are given the correct choice of $v$ as a parameter. To be precise, the parameter is actually some discriminatory polynomial $f$ which can of course be chosen as $\langle{v},X\rangle$. This kernel step is formulated in Algorithm 5. It will return some component up to an error term that can be made infinitely small as $d$ grows.

The procedure itself did not alter that much. We will start again by proving that the first round of the algorithm will indeed yield an approximation to some component of the $T_{k}$. To this end, we will prove the following theorem.

Lemma 5.9**.**

Let $j\in[m]$ and $W$ an SOS polynomial such that $\sum_{i=1}^{m}W(a_{i})=1$ and $\sum_{i\neq j}W(a_{i})\leq C\in(0,1)$. Then each top eigenvector $u$ of $M\mathrel{\mathop{\mathchar 58\relax}}=\sum_{i=1}^{m}\lambda_{i}W(a_{i})a_{i}a_{i}^{T}$ satisfies

[TABLE]

for $\gamma=1-\frac{2C\rho_{\text{spec}}}{\|a_{j}\|^{2}}$. Here, $\rho_{\text{spec}}$ is one of the condition parameters from 3.18. In particular, if $\langle a_{i},u\rangle\geq 0$ then

[TABLE]

Proof.

Denote by $\mu_{\max}$ the eigenvalue of $M$ of largest absolute value. Recall that $\|M\|_{\operatorname{spec}}$ is the absolute value of the largest eigenvalue of $M$. We have thus292929$M$ is psd. Hence the $|\cdot|$ is actually unnecessary.

[TABLE]

We want to use Lemma 5.1. To this end, we need a lower bound for $\|M\|_{\operatorname{spec}}$ and an upper bound for $\|M-a_{j}a_{j}^{T}\|_{\operatorname{spec}}$:

[TABLE]

Dividing by $\|a_{j}\|^{2}$, we get the bound

[TABLE]

Likewise we get, using $\lambda_{j}W(a_{j})-1=-\sum_{i\neq j}\lambda_{i}W(a_{i})$

[TABLE]

Therefore

[TABLE]

and thus we may apply lemma 5.1 with

[TABLE]

∎

It remains to put all of these bounds together to prove our guarantees for the output.

Theorem 5.10**.**

Let $j\in[m]$ such that $f(a_{j})$ is maximal among all $a_{i}$. Then each top eigenvector $u$ of $M\mathrel{\mathop{\mathchar 58\relax}}=\sum_{i=1}^{m}\lambda_{i}W^{\ast}(a_{i})a_{i}a_{i}^{T}$ satisfies

[TABLE]

Proof.

By Lemma 5.8, we have $\sum_{i\neq j}W^{\ast}(a_{i})\leq\mathfrak{w}\cdot\|\lambda\|_{1}\frac{\max_{i\neq j}\Delta_{j,i,f}}{\min_{i\neq j}\Delta_{j,i,f}}$. Therefore we may choose $C\mathrel{\mathop{\mathchar 58\relax}}=\mathfrak{w}\cdot\|\lambda\|_{1}\frac{\max_{i\neq j}\Delta_{j,i,f}}{\min_{i\neq j}\Delta_{j,i,f}}$ in Lemma 5.9 and get thus a vector $u$ satisfying

[TABLE]

out of $M$. ∎

Let us quickly discuss the meaning of the constants occurring in our bound of the error term: $\mathfrak{w}$ is the friendly approximation constant that we can make as small as we want as we let $d\to\infty$ – or perhaps it’s better to say $d\to 2m$. Note that $\mathfrak{w}$ will decrease exponentially with $d$ and that it depends only on $d$ and the condition parameters $\kappa_{\min}$ and $\kappa_{\max}$.303030As long as we are recovering only one component, we do not yet need the safety margin $\varepsilon$, but this will come in soon. As we saw already in the preliminaries, $\rho_{\text{spec}}$ is part of the condition parameters as well and we would have expected something related to the weights like $\|\lambda\|_{1}$ to occur in there anyway. The term $\|a_{j}\|^{2}$ in the denominator is slightly weirder: It tells us that we can’t recover components of small norm. First, we’d actually expect the weight $\lambda_{j}$ to occur there, too, since it should be ok to recover a vector of small length as long as the weight is sufficiently big. We can indeed pull $\lambda_{j}$ off the hat by using the estimate $\mathfrak{w}\leq(\lambda_{j}\mathfrak{r})^{-1}$.

But there is still a problem with vectors close to the zero vector. This seems to be an issue with our recovery technique: These eigenvector methods all rely on vectors of unit length and it’s probably due to the side effects of rescaling that everything gets more unstable for small length components. I’m not entirely sure with this, though. It would thus be preferable if we could without loss of generality assume that the components were of unit length and rescale the weights accordingly. However, to be able to do so while maintaining full generality, we would need to be able to generate the fake moments

[TABLE]

for some $s_{i}\in\{\pm 1\}$, as outlined in §4.3.

Of course the most problematic term is $\frac{1}{\min_{i\neq j}\Delta_{j,i,f}}$. If we chose the wrong discriminatory polynomial, then this term might actually cast a well-conditioned problem into an ill-conditioned one. Therefore, we have to take care that at least we do not make one of the worst possible choices. We dedicate the next section to this problem.

But first we have to show that the length and the factor of $\pm 1$ are correct, too. Note that from Eq. (52) in the proof of Lemma 5.8 we already know that the top eigenvalue $\mu_{1}$ satisfies

[TABLE]

By a very similar calculation, we also get the upper bound

[TABLE]

Hence

[TABLE]

Therefore $|\|c\|^{2}-\|a_{j}\|^{2}|=|\mu_{1}-\|a_{j}\|^{2}|=|\|M\|_{\operatorname{spec}}-\|a_{j}\|^{2}|\leq\mathfrak{w}\cdot\rho_{\text{spec}}\|\lambda\|_{1}\frac{\max_{i\neq j}\Delta_{j,i,f}}{\min_{i\neq j}\Delta_{j,i,f}}$ as claimed in Alg. 5, which shows approximate correctness of the length. For the sign, let again

[TABLE]

and note that

[TABLE]

It’s clear that the scalar product $\langle L,u\rangle$ will be positive if and only if $L$ is closer to $u$ than to $-u$. Equivalently, we can of course replace $L$ by $L\cdot\|a_{j}\|$. Set $\eta\mathrel{\mathop{\mathchar 58\relax}}=\frac{a_{j}}{\|a_{j}\|}$. The correct sign choice for $u$ achieves $\|u-\eta\|\leq\|u-(-\eta)\|$. By Lemma 5.9 we know that there exists a sign choice $s\in\{\pm 1\}$ such that

[TABLE]

Now by the triangle inequality, (55) and (56), we have

[TABLE]

whereas for the wrong sign choice

[TABLE]

where the last step holds due to $\frac{\|L\|}{\|a_{j}\|}\geq 1-\frac{\|L-a_{j}\|}{\|a_{j}\|}\geq 1-C\frac{\sqrt{\kappa_{\max}}}{\|a_{j}\|}$. If $C$ is small enough such that $\tilde{C}\leq 1$ clearly we can distinguish the correct and the wrong $s$ since for the correct one $0\approx\|su-L/\|a_{j}\|\|$ and $\|su+L/\|a_{j}\|\|\gg 1$. We apologise for the sloppiness in this one, but since this is about making a binary decision, we do not have to care about how the exact error looks like as long as it will decrease fast with $d$.

Note that we didn’t give any guarantee regarding correctness of the weight. This is very much intentional: $W$ will usually be a (moderately) high degree polynomial and thus if $c$ is only slightly off from the true component $a_{j}$, $W(c)$ might already be a terrible approximation for $W(a_{j})$. Hence, trusting the continuity of $W$ would be madness! A much better idea is to recover all components $c_{i}$ first and then get the weights by solving a linear least squares programme of the kind

[TABLE]

It’s of course also possible to boost the accuracy of the components by running some local searching procedure afterwards. As Hopkins et al. pointed out in [HSS16], there is a lot of potential for synergy here, since these SOS based procedures output exactly what a local searching procedure needs: Some approximate solution near the global minimum of

[TABLE]

5.4 Choice of $v$

One thing we had to leave open up to now is the choice of $v$, or, more generally, the choice of a suitable discriminatory polynomial $f$. Up to now, we’ve been always trying to keep full generality. However, we currently do not know of any choice of discriminatory polynomials that is universally good in the fully general case. Let us switch therefore to the case $\|a_{i}\|=1$ for all $i\in[m]$. Also, it turns out that the choice $f=\langle{v},X\rangle^{2}$ instead of $\langle{v},X\rangle$ will make our life easier. Let us therefore redefine

[TABLE]

Now, we should first estimate the chances to find a good $v\sim\mathcal{U}$ at random. Precisely, we’re interested in the conditional probability

[TABLE]

that $\Delta_{j,i,v}$ attains a value greater or equal to some $r\in\mathbb{R}_{\geq 0}$ subject to the event $\langle a_{j},v\rangle^{2}\geq\langle a_{i},v\rangle^{2}\geq\langle a_{k},v\rangle^{2}$ for all $k\neq i,j$. We will brute-force-estimate this probability by just considering the event

[TABLE]

that $a_{j}$ is $(1-\gamma)$-square-correlated with $v$ for some small $\gamma\in(0,1)$. This event has probability

[TABLE]

which is constant in $n$ and $m$. From that, we will indeed be able to show that whenever $\mathcal{E}$ happens, then $\langle a_{j},v\rangle^{2}-\langle a_{i},v\rangle^{2}\geq r$ for some $r$ (here we use $\|a_{i}\|^{2}=\|a_{j}\|^{2}$). This is precisely formulated in the following lemma. Note that this estimation is really brute, since we didn’t even condition on the event $\forall k\neq i,j\mathrel{\mathop{\mathchar 58\relax}}\langle a_{j},v\rangle^{2}\geq\langle a_{i},v\rangle^{2}\geq\langle a_{k},v\rangle^{2}$.

Lemma 5.11**.**

Suppose all $\|a_{i}\|=1$ and $v\in\mathcal{E}\mathrel{\mathop{\mathchar 58\relax}}=\{v\in\mathbb{S}^{n-1}\mid\langle a_{j},v\rangle^{2}\geq(1-\gamma)\}$ for some fixed $j\in[m]$ and $\gamma\in(0,1)$. Let $\rho_{\text{minspec}}\mathrel{\mathop{\mathchar 58\relax}}=\min_{i,j\in[m],i\neq j}d_{\mathbb{P}^{n-1}}(a_{i},a_{j})^{2}$. Then

[TABLE]

Proof.

By the “square” triangular inequality for $d_{\mathbb{P}^{n-1}}(\cdot,\cdot)^{2}$:

[TABLE]

we have that

[TABLE]

where for the fist estimate we choose $x=a_{i},y=a_{j},z=v$. ∎

This shows $\langle a_{j},v\rangle^{2}-\langle a_{i},v\rangle^{2}>0$ for sufficiently small $\gamma$. We get a bound just in terms of the correlation.

It’s an interesting question to ask to what extent we need prior knowledge on the domain of the components in order to construct a good discriminatory polynomial. We’ll see that we can use the bound above to search for moment decompositions where the components lie on the unit sphere. This suggests that prior knowledge on the domain can do a lot, but recall that we didn’t need it in the exact case.

Conditioning

Now, of course we could try to run the SOS programme with a random choice of $v$ and see afterwards whether or not it worked. However, SOS programmes are expensive to solve and therefore we will precondition the choice of $v$ such that we get a success probability of 1 again, provided certain conditions are met. To this end, let us observe that if $k$ is even, any candidate for $v$ that is highly correlated with some $a_{i}$ (i.e. $\langle a_{i},v\rangle^{2}\geq(1-\delta)\|a_{i}\|^{2}$ for some small $\delta\in(0,1)$) will attain a high value of $T_{k}(v)$. Indeed, this is due to the estimate that for any component $a_{i}$ and any even $k$, we have

[TABLE]

Hence an obvious precondition would be to filter out all choices that attain a low value of $T_{k}(v)$. Once we’re left with a unit vector $v$ on which $T_{k}$ attains a high value, we can do the following, quantitatively weaker, converse implication:

By a simple averaging / convexity argument that we adapted from [BKS15], we know that

[TABLE]

This is true since $T_{0}^{-1}T_{k}(v)=\sum_{i=1}^{m}\frac{\lambda_{i}}{\|\lambda\|_{1}}\langle a_{i},v\rangle^{k}$ is a convex combination of the $\langle a_{i},v\rangle^{k}$. It implies that $|\langle a_{j},v\rangle|\geq\sqrt[k]{\frac{T_{k}(v)}{\|\lambda\|_{1}}}$. We have thus seen that it is feasible to condition on

[TABLE]

for any fixed $\delta\in(0,1)$ and any $R\in[\min_{i\in[m]}\sqrt[k]{\lambda_{i}}\|a_{i}\|,\>\max_{i\in[m]}\sqrt[k]{\lambda_{i}}\|a_{i}\|]$ with constant success probability at least $\frac{2}{\pi}\arccos(\sqrt{1-\delta})$. Recall that we had the unit sphere assumption. Thus we can choose $R$ s.t. $R^{k}=\min_{i\in[m]}\lambda_{i}=\mathrel{\mathop{\mathchar 58\relax}}\lambda_{\min}$.313131It might be surprising that we do not take the maximum. The reason for this is that in the final algorithm, we will need that any $i\in[m]$ can beat the threshold.

Of course the process of sampling and discarding takes time, but first of all it’s constant time and secondly it is likely not near as much time as solving an SOS programme. For that reason, we will ignore the effect this conditioning has on the running time. We will choose $k$ as large as possible, that is, $k=d$. We get that

[TABLE]

is feasible and happens with probability at least $\frac{2}{\pi}\arccos(\sqrt{1-\delta})$. Note that this is an asymptotically optimal bound, since for $d\to\infty$ we get $r_{\delta}\to 1-\delta$. Choosing $\gamma$ such that $1-\gamma=r_{\delta}$ and plugging this into Lemma 5.11, we get

[TABLE]

Let us now forget about the parameter $\delta$ by setting $r\mathrel{\mathop{\mathchar 58\relax}}=r_{10^{-3}}$.323232The $\delta$ of the next section will have a different meaning.

5.5 Recovering all Components

Note that for any vectors of variables $X,Y,Z$ the SOS triangle inequality

[TABLE]

holds. Indeed, by taking the triangle inequality

[TABLE]

to the power 2 it follows by some easy algebraic operations that (61) holds thus true for all $x,y,z\in\mathbb{R}^{n}$. Since the identity (61) is globally valid and of degree 2, it can be written as a sum of squares.

Lemma 5.12**.**

Let $W_{\text{test}}$ be the testimony polynomial we chose in Prop. 5.7 w.r.t. some $a_{i}\neq a_{j}$ and w.r.t. the interval $[a,b]$ (where $0<a<\frac{1}{2}\kappa_{\min}$ and $b>2\kappa_{\max}$). It satisfies the constraint

[TABLE]

for some $g\geq 0$ with $\deg(g)=d-4$ provided

[TABLE]

Proof.

Let us first search for a univariate identity of the kind

[TABLE]

for some $g\succeq 0$ Note that such an identity would certify

[TABLE]

for all $x\in[a,b]$. Let us strengthen (62) to

[TABLE]

This is satisfied by the rational function

[TABLE]

which is actually a polynomial: Indeed, recall that $a$ and $b$ are the boundary maximisers of the Chebyshev polynomial with respect to $\mathcal{I}$. Hence, they are roots of $\mathfrak{w}-\mathcal{C}_{d}$. It remains to show that $g$ is a sum of squares, which is the equivalent to being non-negative in the case of one variable. For example by plotting $\mathfrak{w}-\mathcal{C}_{d}$ (or by an easy argumentation) it is clear that this polynomial attains non-negative values precisely on $[a,b]$ (recall $d\equiv 2\mod 4$). The same is true for $(\Lambda-a)(b-\Lambda)$. Hence, their quotient is globally non-negative, i.e. $g\succeq 0$. Concatenated with the squared norm, we get

[TABLE]

where $\tilde{g}\mathrel{\mathop{\mathchar 58\relax}}=g(\|X-a_{i}\|^{2})$ Now, let us lower bound the right hand side. By the SOS triangle inequality we have

[TABLE]

Therefore,

[TABLE]

Let us consider any $\delta>0$ such that $\|c-a_{i}\|^{2}-2a\geq 2\delta$ and $b-2\|c-a_{i}\|^{2}\geq 2\delta$.333333The factor of 2 is ugly and much likely due to the fact that the norm is not a polynomial. It infers that we can’t take a Chebyshev polynomial with respect to the interval $[\kappa_{\min}-\delta,\kappa_{\max}+\delta]$, but that we have to leave a bigger safety margin $[\frac{1}{2}\kappa_{\min}-\delta,2\kappa_{\max}+\delta]$. Then we can further estimate:

[TABLE]

Thus by multiplying two SOS polynomials and by expanding:

[TABLE]

The last step uses $6\delta\leq b-4a$, which we get by adding up the two constraints on $\delta$ we imposed in such a way that $\|c-a_{i}\|^{2}$ cancels. We get

[TABLE]

and therefore

[TABLE]

[TABLE]

∎

Thus, we have to show that we can choose $\delta$ greater or equal to the distance $\|a_{j}-c\|^{2}$ such that any feasible polynomial of the second round has to attain a small value at the already recovered $a_{j}$. Now we are left with the problem that the approximation accuracy depends on the interval size $b-a$ which depends on $\delta$ which again depends on the approximation accuracy. So we need to get rid of all the cross-dependencies! We’ll do so by bounding everything with respect to a constant $\rho_{\text{minspec}}$ chosen in such a way that343434The factor of $2$ is technical convenience.

[TABLE]

Also, we will sometimes need to convert distances $x$ w.r.t. $d_{\mathbb{P}^{n-1}}^{2}$ to distances with respect to $d_{\mathbb{S}^{n-1}}^{2}$.353535Note that $\|\cdot\|^{2}=2d_{\mathbb{S}^{n-1}}^{2}$ on $\mathbb{S}^{n-1}$, which we will also use. This is done via the function $F(x)\mathrel{\mathop{\mathchar 58\relax}}=1-\sqrt{1-x}$. It might be difficult to read, but there is no way around it, so please do not be confused.

This turns all out to be very technical, but it’s done in Algorithm 6. The parameters shall be seen as follows: For each choice of parameters $(\rho_{\text{minspec}},\lambda_{\min})$ we get a different algorithm where lenient parameter choices can decompose a bigger class of tensors but require a higher value of $d$ to achieve the same approximation accuracy. E.g. for the maximally restrictive choice $2\rho_{\text{minspec}}=1$ we’d get an algorithm for orthogonal tensor decomposition.

The parameter $\rho_{\text{minspec}}$ is needed just for the “distinct recovery constraint” which ensures that we get each time a different component. Precisely, we need $\rho_{\text{minspec}}$ and $\lambda_{\min}$ to compute the constant $\mathfrak{w}=(\mathfrak{r}\lambda_{\min})^{-1}$, which is in this algorithm the same for all components, and an estimate $\tilde{\varepsilon}$ for the approximation error ($\rho_{\text{minspec}}$ is needed for nothing else!). Thus, it would be interesting to know if we could work completely without $\rho_{\text{minspec}}$. On the other hand, $\lambda_{\min}$ is also used to have an “accepting threshold” $T_{d}(v)\geq(1-1/1000)^{d/2}\lambda_{\min}$ such that we can condition on

[TABLE]

Now, let us deal with the cross-dependencies: For fixed $\rho_{\text{minspec}}$ define for brevity $\rho\mathrel{\mathop{\mathchar 58\relax}}=\frac{1}{2}F(\rho_{\text{minspec}})=\frac{1}{2}(1-\sqrt{1-\rho_{\text{minspec}}})\in\mathcal{O}(\rho_{\text{minspec}})$.363636The factor of $\frac{1}{2}$ is again technical convenience.

For fixed $m,\rho_{\text{minspec}},\lambda_{\min},T_{0}\>(=\|\lambda\|_{1})$ we say that $d$ satisfies the three useful constraints, if:

1. (I)

$0.25\rho_{\text{minspec}}\geq 2(1-r)$, where $r=(1-1/1000)\sqrt[d/2]{\frac{\lambda_{i}}{\|\lambda\|_{1}}}$ 2. (II)

$\rho\geq 4F(\tilde{\varepsilon})$, where $\tilde{\varepsilon}\mathrel{\mathop{\mathchar 58\relax}}=\frac{16T_{0}}{\rho_{\text{minspec}}\cdot\mathfrak{r}\cdot\lambda_{\min}}$  and

$\mathfrak{r}\mathrel{\mathop{\mathchar 58\relax}}=\frac{1}{2}\mathcal{C}_{d}(\frac{8+2\rho}{8-\rho})\in\mathcal{O}((1+\frac{3\rho}{8-\rho})^{d})=\mathcal{O}((1+\rho_{\text{minspec}})^{d})$ 3. (III)

$0.25\rho_{\text{minspec}}\geq 4m\sqrt{2-2(1-\tilde{\varepsilon})^{d/2}}$

The constraints are needed to ensure that the approximation error is small enough such that the algorithm stays stable over $m$ rounds. All of these constraints are feasible for sufficiently high $d$ (supposed that $0.25\rho_{\text{minspec}}>2/1000$ for the first one. Note that $(1-(1-\tilde{\varepsilon})^{d})\in\mathcal{O}(\tilde{\varepsilon})$ for the last one and that $\tilde{\varepsilon}$ decreases exponentially with $d$). But $d$ does not even need to be very high (up to constants): We get a logarithmic dependency of $d$ on the parameters.

Note that the ratio $\mathfrak{r}$ is correctly estimated: With the testimony from Section §5.1 we get that the interval $\mathcal{I}=[a,b]$ can actually be chosen as

[TABLE]

as long as the approximation accuracy w.r.t. $d_{\mathbb{P}^{n-1}}^{2}$ is less or equal $\tilde{\varepsilon}$. Using the estimates $\kappa_{\max}=\max_{i\neq j}\|a_{i}-a_{j}\|^{2}\leq 2^{2}=4$ and $\frac{1}{2}\kappa_{\min}=\frac{1}{2}\min_{i\neq j}\|a_{i}-a_{j}\|^{2}=\min_{i\neq j}d_{\mathbb{S}^{n-1}}(a_{i},a_{j})^{2}\geq(1-\sqrt{1-\rho_{\text{minspec}}})=2\rho\in\mathcal{O}(\rho_{\text{minspec}})$, we see that

[TABLE]

would actually be valid373737This shall mean: It is attained for a testimony in the feasible space. if the true approximation error $\varepsilon$ (w.r.t. $d_{\mathbb{P}^{n-1}}^{2}$) was less or equal $\tilde{\varepsilon}$. We simplified this to

[TABLE]

which will be valid if and only if we can show that it is possible to choose $d$ so high that the “converted” approximation error $4F(\varepsilon)$ is less or equal to $\rho$. The second useful constraint states $\rho\geq 4F(\tilde{\varepsilon})$, so it remains to show that the approximation error gets smaller than $\tilde{\varepsilon}$ for sufficiently high $d$.

Now, let’s show this: Suppose the parameters $\rho_{\text{minspec}}$ and $\lambda_{\min}$ were correctly specified for given input tensors $T_{0},\ldots T_{d}$ with exact decomposition $T_{k}=\sum_{i=1}^{m}\lambda_{i}\langle{a}_{i},X\rangle^{k}$. Furthermore, suppose that $d$ satisfies the 3 useful constraints with respect to the given input and parameter choices.

Choose the first round’s $v\in\mathbb{S}^{n-1}$ such that $T_{d}(v)\geq(1-\frac{1}{1000})^{d/2}\lambda_{\min}$. By the results of §5.3, we can recover a proxy $c_{1}=u$ to the component $a_{j}$ maximising $\langle a_{j},v\rangle^{2}$ such that

[TABLE]

Note that $\max_{i\neq j}\Delta_{j,i,v}\leq\max_{i\neq l}\langle a_{l},v\rangle^{2}-\langle a_{i},v\rangle^{2}\leq\rho_{\text{spec}}$ by definition and since we are on the unit sphere even $\rho_{\text{spec}}\leq\max_{i\neq l}(\|a_{i}a_{i}^{T}\|_{\operatorname{spec}}+\|a_{l}a_{l}^{T}\|_{\operatorname{spec}})=2$. Since $\min_{i\neq j}\Delta_{j,i,v}\geq\rho_{\text{minspec}}-2(1-r)$, we get

[TABLE]

and hence $d_{\mathbb{P}^{n-1}}(c_{1},a_{j})^{2}\leq\varepsilon$ with accuracy at least $\varepsilon=\frac{8\cdot\|\lambda\|_{1}}{\mathfrak{r}\cdot\lambda_{\min}\cdot(\rho_{\text{minspec}}-2(1-r))}$. From the first useful constraint we get that

[TABLE]

and therefore383838Note that in this estimation we’ve left some extra space “$0.75\rho_{\text{minspec}}\geq 0.5\rho_{\text{minspec}}$” which we will need for the subsequent rounds to work.

[TABLE]

will indeed be sufficiently small (at least in the first round). We have already argued before that the decision between $\pm u$ is not the most critical issue here, so let us just assume that it is correct for sufficiently small $\tilde{\varepsilon}$.393939It would be a problem though if the moment decomposition could have distinct components satisfying $a_{i}=-a_{l}$ for $i\neq l$. (Note that this case is excluded due to $\min_{i\neq j}d_{\mathbb{P}^{n-1}}(a_{i},a_{j})^{2}\geq\rho_{\text{minspec}}$). Then $W^{\ast}$ could concentrate on both of them equally, since the discriminatory polynomial $\langle{v},X\rangle^{2}$ would satisfy $\langle a_{j},v\rangle^{2}=\langle a_{l},v\rangle^{2}$. This problem is a suboptimality that stems from the fact that we weren’t able to prove in §5.4 that we could require $\exists j\mathrel{\mathop{\mathchar 58\relax}}\langle a_{j},v\rangle\geq\sqrt{r}$ instead of $\exists j\mathrel{\mathop{\mathchar 58\relax}}\langle a_{j},v\rangle^{2}\geq r$. But then we have

[TABLE]

Knowing $\mathfrak{w}$ and $\tilde{\varepsilon}$, we can write down the constraint

[TABLE]

by adding a new “slack variable” $g$ to the SOS programme. This constraint ensures indeed that the optimiser of the 2nd round satisfies $W^{\ast}(a_{j})\leq\mathfrak{w}$ (such that we will not recover the same component twice), since the distance of $c_{1}$ and $a_{j}$ is less or equal $2-\leavevmode\nobreak\ 2\sqrt{1-\tilde{\varepsilon}}$. As we saw in §5.5, we needed a safety margin of at least $2F(\varepsilon)\leavevmode\nobreak\ =\leavevmode\nobreak\ 2-\leavevmode\nobreak\ 2\sqrt{1-\varepsilon}$ in the interval of the Chebyshev polynomial to ensure that the testimony satisfies the distinct recovery constraint. This is satisfied since we saw $\varepsilon\leq\tilde{\varepsilon}$.

It is left to show that this argumentation does not break until the $m$-th round. The only thing that changes from round 2 and onwards is that we can’t just take some $v$ which achieves just a large value of $T_{d}(v)$: We need to ensure that $v$ is highly correlated with some component $a_{i}\neq a_{j}$. This would be easy if we had access to the tensor $S_{\text{correct}}\mathrel{\mathop{\mathchar 58\relax}}=\sum_{i\neq j}\lambda_{i}\langle{a_{i}},X\rangle^{d}$. To this end, we need the 3rd useful constraint which essentially guarantees us that $S-\langle{u},X\rangle^{d}=S-\langle{c_{1}},X\rangle^{d}$ is a sufficiently good approximation to $\sum_{i\neq j}\lambda_{i}\langle{a_{i}},X\rangle^{d}$. We can show that they have to attain similar values on $v$ by using the evaluation property of Reznick’s scalar product (Prop. 3.31):

[TABLE]

In particular, since

[TABLE]

we get that in the 2nd round we can sample with the same success probability as in round 1 some $v\in\mathbb{S}^{n-1}$ such that $S(v)\geq r-\mathcal{O}(\sqrt{\tilde{\varepsilon}})$ and any such $v$ satisfies $S_{\text{correct}}(v)\geq r-2\mathcal{O}(\sqrt{\tilde{\varepsilon}})$.

This shows that in the end, what has changed is that we have to start our estimations for $\tilde{\varepsilon}$ with a slightly lower value of $r$. We’ll show that our estimations for $\tilde{\varepsilon}$ were conservative enough such that they do not have to be changed. For the sake of simplicity, we demonstrated the above for round 2. Generally, in round $i$ we are left with some $v\in\mathbb{S}^{n-1}$ such that

[TABLE]

Set $r^{\prime}\mathrel{\mathop{\mathchar 58\relax}}=r-2m\sqrt{2-2(1-\tilde{\varepsilon})^{d/2}}$. This is where the 3rd useful constraint comes in: It guarantees that $0.25\rho_{\text{minspec}}\geq 4m\sqrt{2-2(1-\tilde{\varepsilon})^{d/2}}=2(r-r^{\prime})$, whence in each round we can still deduce the estimate

[TABLE]

from Equation (63). Hence, the approximation guarantee $\tilde{\varepsilon}$ will not get lower over the course of the $m$ rounds of the algorithm.

5.6 Noise Stability and Tensor Decomposition

In the current state, the algorithm is not very stable to noise: For instance, a single component with negative weight could make the optimisation problem unbounded. This can be changed, though: We just need to cap the Frobenius norm $\|W\|_{F}\leq D$ for some constant $D$.404040Then we also have $\|{W\langle V,X\rangle^{2}}\|_{F}\leq\tilde{D}$ for some polynomial $\tilde{D}\in\mathbb{R}[V]_{\leq 2}$. The cap needs of course to be chosen such that the testimonies are still members of the feasible space. Capping can be realised for instance by adding linear constraints on the coefficients of $W$ or by constraining the trace of the Gram matrices in the SDP. Indeed, if the correct input polynomial $\sum_{k=0}^{d}T_{k}$ is disturbed by some noise polynomial $E$ where $\|E\|_{F}$ is several orders of magnitude smaller than the Frobenius norm of the input polynomial, then we have by Cauchy-Schwarz, that

[TABLE]

Of course the Lemmas giving guarantees on the optimiser would have to be reproven while carrying the additional noise term with us. A full-fledged noise analysis would go beyond the scope of this thesis. Philosophically spoken, capping the Frobenius norm limits the damage that a wrong summand can do: In the special case that the noise term is, for instance, a single component $a_{i}$ of unit length with small negative weight $-1\ll\lambda_{i}<0$, we get that

[TABLE]

Hence $\lambda_{i}W(a_{i})$ will be very small.

Regarding tensor decomposition, it should be said that the fake moment generation procedure we gave in §4.2 is not the most ideal choice for the setting of approximate recovery. The problem is again that the random choice of $w$ can cast a well-conditioned problem into an ill-conditioned one. Therefore it’s better to use the canonical norm-scaled fake moments of even degree that we presented in §4.3.

This also allows us to assume without loss of generality that the components have unit length, which we have seen to be very practical. However, this means that we’d have to do some more work for the case of tensor decomposition: First of all, we need a new testimony having only even degree coefficients. We can employ essentially the same idea as in §4.3: Take a univariate Chebyshev polynomial $\mathcal{T}_{\mathcal{I},d}+1$ with respect to the interval $\mathcal{I}\mathrel{\mathop{\mathchar 58\relax}}=[0,\max_{i\neq j}\leavevmode\nobreak\ \|a_{i}\|^{2}\leavevmode\nobreak\ \|a_{j}\|^{2}\leavevmode\nobreak\ -\leavevmode\nobreak\ \langle a_{i},a_{j}\rangle^{2}]$ and concatenate it with the “Cauchy-Schwarz gap” $\|X\|^{2}\|a_{j}\|^{2}-\langle X,a_{j}\rangle^{2}$. This will separate the components up to their length and a factor of $\pm 1$. Of course we’d still have to prove that a similar sums of squares constraint as in §5.5 exists such that we can avoid to recover the same component twice.

These topics may be addressed in detail in some future research.

6 Conclusions and Open Problems

Sums of squares optimisation has shown to be a valuable tool for moment decomposition when sufficient data is given. The algorithmic approximation scheme we developed in Section §5 can compute the decomposition to arbitrary high accuracy, as long as we have access to all the moments of arbitrary high degree. If only a finite number of moments is accessible, then the approximation accuracy is fixed, but can eventually be increased by running some local searching procedure afterwards. The error depends on the condition of the problem, whence decompositions with well-separated components are easier to recover. This is a very much expected phenomenon, even though it does not occur in more elementary procedures such as Jennrich’s, where the working requirement is just a relation between $m$ and $n$.

On the other hand, on many instances the algorithm might produce results way better than the worst case bounds suggest. In particular, we’ve seen that for certain instances (e.g. highly symmetric ones such as the components lying on a simplex) exact recovery is possible already for very low degree.

If this thesis got them interested, the reader is encouraged to dedicate further research to this topic. The question which is probably the most interesting is whether some variant of the $v$-algorithm could be combined with a procedure for the generation of higher order moments to further improve on the approximation accuracy – at the expense of larger computation time. In the following, we’ve collected some opportunities related directly to the $v$-algorithm.

Related to the $v$-Algorithm

Although the algorithm works, there is still a lot of potential for optimisation. We would like to collect the most notable suboptimalities of the algorithm to motivate further research:

The issue likely to be the most important is to find strategies for the choice of $v$ other than just sampling it at random. To this end, it can be helpful to require some prior knowledge on the domain of the $a_{i}$. We would particularly be interested in deterministic procedures for the choice of discriminatory polynomials. 2. 2.

In hindsight, some of our techniques seem to be way more suitable for even degree tensor decomposition than for moment decomposition: This is particularly true for the conditioning we performed. It seems also that the radially symmetric testimonies we chose do not work too well with most of the other correlation-based estimates. While the arguments could be made to work, many technical estimates were needed.

Everything could potentially become simpler when considering tensor decomposition. We’ve seen that the unit sphere assumption is not that much of a restriction in this case due to the possibility to generate the even degree norm-scaled fake moments. Furthermore, then we can concatenate our univariate Chebyshev polynomials with the Cauchy-Schwarz gap, which has the advantage that it depends directly on the square correlation. This has the potential to make the analysis much easier. 3. 3.

The choice of radially symmetric Chebyshev polynomials seems weirdly heuristic. On the one hand, restricting to radially symmetric polynomials is a potential waste of the valuable and computationally expensive degree: After the substitution with the squared norm we can only plug-in polynomials of degree less or equal $d-2$. On the other hand, it’s likely that multivariate polynomials will achieve better separation results. 4. 4.

It also comes with the seemingly unnecessary restriction $d\equiv 2\mod 4$: If we would have designed the algorithm from the very beginning for moment decomposition on the unit sphere, then we could have used that on the unit sphere the metric $d_{\mathbb{S}^{n-1}}(x,a_{j})^{2}$ can be represented by the linear polynomial $1-\langle X,a_{j}\rangle$. If we had concatenated the Chebyshev polynomials with these linear polynomials instead of the norm, our life would have been much easier in the setting of our main result.

Declaration of Authorship

I declare that the submitted thesis

*A new Algorithm for Overcomplete Tensor Decomposition based on Sums-of-Squares Optimisation

is my own unaided work. Direct or indirect sources are acknowledged as references. This thesis was not previously presented to another examination board and hasn’t been published before.

Konstanz, 5 October 2018

.

Alexander Taveira Blomenhofer

Supplement I - Existence of Optimal Solutions

We did not yet give an argument why we can always assume the occuring SOS optimisation problems to have an optimal solution. While in the exact case, we were able to give explicit examples (the interpolation polynomials), for the general case it is actually not clear whether or not optimal solutions exist. However, note that in all of our analysis we did merely use one property of the optimal solution $W^{\ast}$, namely that the objective function achieves on $W^{\ast}$ some value at least as high as the values on all of the testimony polynomials

[TABLE]

Let us reconsider the basic SOS optimisation problem from above

[TABLE]

With the expectation operator $\mathbb{E}_{\mu}$ of $\mu=\sum_{i=1}^{m}\lambda_{i}\delta_{a_{i}}$, this can be rewritten as

[TABLE]

One problem is that the feasible space is in general not compact: Indeed, assume there is some sum of squares polynomial $W_{0}$ vanishing on every $a_{i}$. Then

[TABLE]

Thus it could be possible to traverse the feasible space in a direction where the value of the objective function does not change. Doing so would of course yield an unnecessary increase of the optimisers “complexity” (complexity can be measured by an appropriate norm on $\mathbb{R}[X]_{\leq d-2}$). This can easily be prevented by adding a complexity bound.414141It is not clear whether such a complexity bound will affect the quality of the optimal solution. However, we are fine with an optimal solution that is “at least as good” as the testimonies. The most basic approach in doing so would be to cap the absolute values of the coefficients of $W$ using linear constraints. We will use a slightly different approach: With respect to any continuous measure $\nu$, $\mathbb{E}_{\nu}[W]$ will be strictly positive for any non-zero SOS polynomial $W$. Therefore we can take for instance $\nu=\mathcal{U}$ as the uniform probability measure on the set $B\mathrel{\mathop{\mathchar 58\relax}}=\operatorname{conv}\mathbb{S}^{n-1}$ and cap the size of $\mathbb{E}_{\mathcal{U}}[W]$. We just have to verify that every testimony polynomial still lies in the feasible space. Hence, we should upper bound

[TABLE]

for an arbitrary unit vector $a\in\mathbb{S}^{n-1}$ and any reasonable choice of the interval $I$. (Recall that for $\mu=\sum_{i=1}^{m}\lambda_{i}\delta_{a_{i}}$ we used $I=[\min_{i\neq j}\|a_{i}-a_{j}\|^{2},\>\max_{i\neq j}\|a_{i}-a_{j}\|^{2}]\subseteq(0,4]$ if all $a_{i}\in\mathbb{S}^{n-1}$). W.l.o.g. we may assume $a=e_{1}$ since $\mathcal{U}$ is rotation-symmetric. $\mathcal{T}_{\mathcal{I},(d-2)/2}$ attains its maximum on $[0,4]$ on the points that have the largest distance to the midpoint of $\mathcal{I}$. Hence $\max_{x\in[0,4]}\mathcal{T}_{\mathcal{I},(d-2)/2}(x)\leq\max\{\mathcal{T}_{\mathcal{I},(d-2)/2}(0),\mathcal{T}_{\mathcal{I},(d-2)/2}(4)\}=\mathrel{\mathop{\mathchar 58\relax}}C_{\max}$.

But then we get the brute upper bound that for every $I\subseteq[0,4]$ and $a\in\mathbb{S}^{n-1}$:

[TABLE]

And therefore for $W_{\text{test}}=\mathbb{E}_{\mu}[\mathcal{T}_{\mathcal{I},(d-2)/2}(\|X-a_{j}\|^{2})+1]^{-1}\left(\mathcal{T}_{\mathcal{I},(d-2)/2}(\|X-a_{j}\|^{2})+1\right)$ we have:

[TABLE]

due to

[TABLE]

The last estimate holds since the Chebyshev polynomial is nonnegative outside of $\mathcal{I}$.

With this bound we get a new complexity-truncated SOS optimisation programme

[TABLE]

By adding the new constraint, the feasible space gets compact with respect to the norm $\|W\|_{\mathcal{U}}\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}_{\mathcal{U}}[|W|]$ on $\mathbb{R}[X]_{\leq d-2}$ and therefore there exists an optimiser $W^{\ast}$ of $(P_{\mu,v,C})$.