Statistical and computational phase transitions in spiked tensor estimation

TL;DR

This paper rigorously analyzes phase transitions in spiked tensor estimation, showing when Bayesian methods and AMP algorithms succeed or fail, revealing a wider 'easy' region than previously known.

Contribution

It provides a rigorous computation of mutual information and MMSE, and characterizes the algorithmic and information-theoretic phase transitions in tensor factorization.

Findings

AMP achieves MMSE in a large parameter region

Existence of a 'hard' region where AMP fails

Factorization is algorithmically 'easy' over a wider range

Abstract

We consider tensor factorizations using a generative model and a Bayesian approach. We compute rigorously the mutual information, the Minimal Mean Squared Error (MMSE), and unveil information-theoretic phase transitions. In addition, we study the performance of Approximate Message Passing (AMP) and show that it achieves the MMSE for a large set of parameters, and that factorization is algorithmically "easy" in a much wider region than previously believed. It exists, however, a "hard" region where AMP fails to reach the MMSE and we conjecture that no polynomial algorithm will improve on AMP.

Figures

Tables

Table 1. Table 1: Examples of the information-theoretic Δ IT subscript Δ IT \Delta_{\rm IT} and algorithmic Δ Alg subscript Δ Alg \Delta_{\rm Alg} thresholds for order- p 𝑝 p tensor decomposition for different priors on the factors. For the Gaussian case Δ IT ​ p ​ log ⁡ ( p ) subscript Δ IT 𝑝 𝑝 \Delta_{\rm IT}p\log(p) converges to 1 at large p 𝑝 p . For the Bernoulli case the rescaling in power of ρ 𝜌 \rho is for convenience to present quantities of order one, we did not check if it describes the large p 𝑝 p limit.

	Gaussian $𝒩(0,1)$		Rademacher		Bernoulli $\rho =0.1$		3 clusters
	${\Delta }_{\mathrm{IT}}p\log \left(p\right)$	${\Delta }_{\mathrm{Alg}}$	${\Delta }_{\mathrm{IT}}$	${\Delta }_{\mathrm{Alg}}$	${\Delta }_{\mathrm{IT}}{\rho }^{-p}$	${\Delta }_{\mathrm{Alg}}{\rho }^{-2p+2}$	$\frac{{\Delta }_{\mathrm{IT}}}{{\Delta }_{\mathrm{Alg}}}$	$\frac{{\Delta }_{\mathrm{Alg}}{r}^{2p-2}}{p-1}$
$2$	$2\log 2$	$1$	$1$	$1$	$-$	$-$	$1$	$1$
$3$	$0.754$	$0$	$0.2828$	$0$	$0.577$	$3.738$	$1$	$1$
$4$	$0.701$	$0$	$0.1902$	$0$	$0.398$	$6.017$	$1.18$	$1$
$5$	$0.685$	$0$	$0.1473$	$0$	$0.311$	$8.251$	$1.62$	$1$
$10$	$0.677$	$0$	$0.07216$	$0$	$0.154$	$19.30$	$6.59$	$1$

Equations

$$Y=\frac{\sqrt{(p-1)!}}{N^{\frac{p-1}{2}}}\sum _{k=1}^r(X_k^0{)}^{\otimes p}+V,$$

$$W_{i_1,i_2,\dots ,i_p}^0=\frac{\sqrt{(p-1)!}}{N^{\frac{p-1}{2}}}\sum _{k=1}^r{x}_{i_1,k}^0{x}_{i_2,k}^0\cdots x_{i_p,k}^0.$$

$$P(Y|X^0)=\prod _{i_1<i_2<\cdots <i_p}{P}_{\mathrm{out}}\left(Y_{i_1,i_2,\dots ,i_p}|W_{i_1,i_2,\dots ,i_p}^0\right).$$

$$P(X|Y)=\frac{1}{{\mathcal{Z}}_N}\prod _{i=1}^N{P}_X(x_i)\prod _{i_1<i_2<\cdots <i_p}{P}_{\mathrm{out}}\left(Y_{i_1,i_2,\dots ,i_p}|W_{i_1,i_2,\dots ,i_p}\right),$$

$$S$$

$$\frac{1}{\Delta }$$

$$\mathcal{M}(x)=\frac{1}{{\mathcal{Z}}_X(A,B)}{P}_X(x)e^{B^{\top }x-\frac{x^{\top }Ax}{2}},$$

$$u\cdot v=\frac{1}{N}\sum _{j=1}^n{u}_j{v}_j^{\top }\in {\mathbb{R}}^{r\times r}.$$

$$B_i^t$$

$$A^t$$

$${\hat{x}}_i^{t+1}$$

$${\sigma }_i^{t+1}$$

$$\displaystyle M^{t+1}=\mathbb{E}_{Z,x_{0}}\left[f_{\rm in}\left(\widehat{M}^{t},\widehat{M}^{t}x_{0}+\big{(}\widehat{M}^{t}\big{)}^{1/2}Z\right)x_{0}^{\top}\right],$$

$${\hat{M}}^t={(M^t)}^{\circ (p-1)}/\Delta ,$$

$${\mathrm{MSE}}_{\mathrm{AMP}}=\mathrm{Tr}\left[{\Sigma }_X-M_{\mathrm{AMP}}\right].$$

$${\text{MMSE}}_N=\underset{\hat{\theta }}{\inf }\left\{\frac{1}{N}\mathbb{E}\left[{\Vert X^0-\hat{\theta }(Y)\Vert }^2\right]\right\}=\frac{1}{N}\mathbb{E}\left[{\Vert X^0-\mathbb{E}[X^0|Y]\Vert }^2\right],$$

$$\frac{1}{N}I(X^0;Y)\underset{N\to \infty }{\to }\frac{1}{2p\Delta }\sum _{k,k^{\prime }=1}^r({\Sigma }_X{)}_{k,k^{\prime }}^p-\underset{M\in S_r^{+}}{\sup }{\varphi }_{\text{RS}}(M),$$

$${\varphi }_{\mathrm{RS}}(M)={\mathbb{E}}_{Z,x_0}\left[\log {\mathcal{Z}}_X\left(\hat{M},\hat{M}{x}_0+{\left(\hat{M}\right)}^{1/2}Z\right)\right]-\frac{p-1}{2p\Delta }\sum _{k,k^{\prime }=1}^r{M}_{k{k}^{\prime }}^p$$

$$\mathrm{MMSE}=\underset{N\to \infty }{\lim }{\text{MMSE}}_N=\mathrm{Tr}\left[{\Sigma }_X-M^{\ast }\right].$$

$$M^{t+1}=\frac{1}{\Delta }{\Sigma }_X\left[{(M^t)}^{\circ (p-1)}\right]{\Sigma }_X.$$

$$Y_{i_1,\dots ,i_p}=\frac{\sqrt{(p-1)!}}{N^{(p-1)/2}}{x}_{i_1}^0\dots x_{i_p}^0+V_{i_1,\dots ,i_p}\text{for}1\le i_1<\cdots <i_p\le N,$$

$$H_N(X)={\Delta }^{-1}\sum _{i_1<\cdots <i_p}\frac{\sqrt{(p-1)!}}{N^{(p-1)/2}}{Y}_{i_1,\dots ,i_p}{x}_{i_1}\dots x_{i_p}-\frac{(p-1)!}{2N^{p-1}}{\left(x_{i_1}\dots x_{i_p}\right)}^2,$$

$$dP(X^0=X|Y)=\frac{1}{{\mathcal{Z}}_N}d{P}_X(X)e^{H_N(X)},$$

$$I(X^0;Y)=\frac{N}{2p\Delta }{\Sigma }_X^p-\mathbb{E}[\log {\mathcal{Z}}_N]+O(1).$$

$${\varphi }_{\mathrm{RS}}:m\ge 0\mapsto \mathbb{E}\left[\log \int d{P}_X(x)\exp \left(\sqrt{\frac{m^{p-1}}{\Delta }}Zx+\frac{m^{p-1}}{\Delta }x{x}^0-\frac{m^{p-1}}{2\Delta }{x}^2\right)\right]-\frac{p-1}{2p\Delta }{m}^p,$$

$$F_N\equiv \frac{1}{N}\mathbb{E}\left[\log {\mathcal{Z}}_N\right]\underset{N\to \infty }{\to }\underset{m\ge 0}{\sup }{\varphi }_{\mathrm{RS}}(m)\equiv F_{\mathrm{RS}}(\Delta ).$$

$${\text{T-MMSE}}_N(\Delta )=\underset{\hat{\theta }}{\inf }\left\{\frac{p!}{N^p}\sum _{i_1<\cdots <i_p}{\left(x_{i_1}^0\dots x_{i_p}^0-\hat{\theta }(Y{)}_{i_1\dots i_p}\right)}^2\right\},$$

$$\frac{\partial F_{N}}{\partial\lambda}=\frac{N(N-1)\dots(N-p+1)}{2pN^{p}}\big{(}\Sigma_{X}^{p}-\text{T-MMSE}_{N}(\Delta)\big{)}.$$

$${\text{T-MMSE}}_N\underset{N\to \infty }{\to }{\Sigma }_X^p-m^{\ast }(\Delta {)}^p.$$

$${\Delta }_{\mathrm{IT}}=\sup \left\{\Delta >0|m^{\ast }(\Delta )>{\mathbb{E}}_{P_X}[x{]}^2\right\}.$$

$$\mathbb{E}\left|{\left(\frac{1}{N}\sum _{i=1}^N{x}_i^0{x}_i\right)}^p-m^{\ast }(\Delta {)}^p\right|\underset{N\to \infty }{\to }0,$$

$${\text{MMSE}}_N\underset{N\to \infty }{\to }{\Sigma }_X-m^{\ast }(\Delta ).$$

$$\mathbb{E}⟨f(Y,X^{(1)},\dots ,X^{(k)})⟩=\mathbb{E}⟨f(Y,X^{(1)},\dots ,X^{(k-1)},X)⟩.$$

Full text

††​​​​​​​​​​​​​This paper was presented at the IEEE International Symposium on Information Theory (ISIT) 2017 in Aachen, Germany.

$\dagger$ Institut de Physique Théorique, CNRS & CEA & Université Paris-Saclay, Saclay, France.

$\star$ Laboratoire de Physique Statistique, CNRS & Université Pierre et Marie Curie & École Normale Supérieure & PSL Université, Paris, France.

$\Diamond$ Département d’Informatique de l’ENS, École Normale Supérieure & CNRS & PSL Research University & Inria, Paris, France.

This study inscribes into the line of research on low-rank tensor decomposition, a problem with many applications ranging from signal processing to machine learning[2, 6, 25]. We consider the model of [24] where the tensor is a noisy version of a $r$-dimensional randomly generated spike and analyze the Bayes-optimal inference of the spike, compute the associated mutual information and the minimum mean-squared error (MMSE). We also investigate whether the MMSE is achievable with some known efficient algorithms, and most particularly by approximate message passing (AMP).

1 The spiked tensor model

One observes an order-$p$ tensor $Y\in\bigotimes^{p}\mathbb{R}^{N}$ created as

[TABLE]

where $X^{0}_{1},\dots,X^{0}_{r}\in\mathbb{R}^{N}$ are $r$ unknown vectors to be inferred from $Y$, and $V\!\in\bigotimes^{p}\mathbb{R}^{N}$ is a symmetric tensor accounting for noise. We denote by $X$ the $N\!\times\!r$ matrix that collects the $r$ vectors $X_{k}$. The observed tensor $Y$ can thus be seen as a rank $r$ perturbation of a random symmetric tensor $V$. Consider now the setting where the $X^{0}$ is generated at random from a known prior distribution. The core question considered in this paper is: What is the best possible reconstruction of $X^{0}$ one can hope for?

In fact, we can look at even more general noise than just additive one as in (1). Denote for $i=1,\dots,N$, $x_{i}=(x_{i,1},\dots,x_{i,r})\in\mathbb{R}^{r}$ the $r$-dimensional vector created by aggregating the $i^{\rm th}$ coordinates of the $r$ vectors $X_{k}$. Assume that for $1\leq i\leq N$ the $x_{i}^{0}$ are generated independently from a probability distribution $P_{X}$ over $\mathbb{R}^{r}$. We denote, for $(i_{1},i_{2},\dots,i_{p})\in\{1,\dots,N\}^{p}$

[TABLE]

For simplicity, we will assume to only observe the extra-diagonal elements of $Y$, i.e. the coefficients $Y_{i_{1},i_{2},\dots,i_{p}}$ for $1\leq i_{1}<\dots<i_{p}\leq N$. The case where all coefficients are observed can be directly deduced from this case. The observed tensor $Y$ is generated from $W^{0}$ using a noisy component-wise output channel $P_{\rm out}$ so that

[TABLE]

The simplest situation corresponds to eq. (1) with additive white Gaussian noise (AWGN), i.e. $P_{\rm out}(\,\cdot\,|\,w)={\cal N}(w,\Delta)$.

Given the above generative model and assuming that both the prior distribution $P_{X}$ and the output channel $P_{\rm out}$ are known we can write the Bayes-optimal estimator of $X^{0}$ as marginalization of the following posterior distribution

[TABLE]

where ${\cal Z}_{N}$ is a normalization constant depending of the observed tensor $Y$, $W_{i_{1},i_{2},\dots i_{p}}$ is defined analogously to (2) (with $X$ instead of $X^{0}$).

We will study this tensor estimation problem in the limit where the dimension $N\!\rightarrow\!\infty$ while the rank $r$ remains constant. The factor $N^{\frac{p-1}{2}}$ is here to ensure that information-theoretically the inference problem is neither trivially hard nor trivially easy when one deals with signals such that $\|x_{i}\|$ and the noise magnitude are of order 1. The factor $\sqrt{(p-1)!}$ is used for convenient rescaling of the signal-to-noise ratio.

2 Related work and summary of results

Recently there have been numerous works on the matrix ($p=2$) version of the above setting. In particular an explicit single-letter characterization of the mutual information and of the optimal Bayesian reconstruction error have been rigorously established [12, 9, 13, 4, 14]. A large part of these results rely on the approximate message passing (AMP) algorithm. For the rank-one matrix estimation problems AMP has been introduced by [23], who also derived the state evolution (SE) formula to analyze its performance, generalizing techniques developed by [5]. In [16, 15, 17] the generalization to larger rank, and general output channel, was considered. Following the theorem proven in [9, 4, 14], we know that indeed AMP is Bayes-optimal and achieves the minimum mean-squared error (MMSE) for a large set of parameters of the problem. There, however, might exist a region denoted as hard, where this is not the case, and polynomial algorithms improving on AMP are not known.

In comparison, there has been much less work on Bayesian low-rank tensor estimation. In statistical physics, the measure (4) was considered for $Y$ with random i.i.d. components. For a Gaussian $P_{X}$, it is called the spherical $p$-spin glass [8], while for Rademacher $P_{X}$ it is the Ising $p$-spin glass [19]. AMP for tensor estimation is actually equivalent to the so-called Thouless-Anderson-Palmer equations in spin glass theory [27, 7, 12]. In the context of tensor PCA these equations have been studied by Richard and Montanari [24] for the maximum likelihood estimation. Interestingly, they showed that the hard phase was particularly large in the tensor estimation case and that, with side information (such that for each component $x_{i,k}$ we have its direct noisy observation), the estimation problem becomes easier. However, such a kind of component-wise side information is very strong and rarely available in applications. The tight statistical limits for the present tensor model were also studied in [12] for the special case of the Rademacher (Ising) prior. For more generic priors only upper and lower bounds are known rigorously [22].

Summary of results: In this contribution, we aim to bridge the gap between what is known for the general $r,P_{x},P_{\rm out}$ Bayesian estimation for low-rank matrices and what is known for low-rank tensors. We present the following contributions:

1. (A)

The AMP algorithm and its state evolution analysis for the Bayes-optimal tensor estimation, see sections 3 and 4. 2. (B)

The so-called channel universality result that allows us to map any generic channel $P_{\rm out}$ on a model with additive Gaussian noise, see section 3. 3. (C)

Rigorous formula for the asymptotic mutual information and the MMSE, thus generalizing the matrix results of [4, 14], see section 5. 4. (D)

The identification of statistical and computational phase transitions. In fact, we show that as soon as the effect of a non-zero-mean prior is taken properly into account, the hard region shrinks considerably, making the tensor decomposition problem much easier than hitherto believed, at least for algorithms that do take the prior information into account. Having a reliable prior information on the distribution of $x_{i}$ (not on each of the components as in [24]) is rather realistic in applications, for instance when constraints of negativity or membership to clusters are imposed. This is presented in sections 4 and 6.

3 AMP algorithm & channel universality

We discuss in this section the Approximate Message Passing (AMP) algorithm for the Bayesian version of the problem. This is a relatively straightforward generalization of what has been done for the low-rank matrix estimation in e.g. [23, 18, 17], i.e. $p=2$ case of the present setting. In general, AMP is derived from belief propagation by taking into account that every variable in the corresponding graphical model has a large number of neighbors. Since the incoming messages are considered independent one can use the central limit theorem and represent each message as a Gaussian with a given mean $\hat{x}_{i}\in\mathbb{R}^{r}$ and covariance $\sigma_{i}\in\mathbb{R}^{r\times r}$.

A crucial property, called channel universality, that the tensor-AMP shares with the low-rank matrix estimation, allows to drastically simplify the problem of tensor estimation with generic output channel $P_{\rm out}$. The justification of this property follows closely the low-rank matrix estimation case, and we refer the reader to [15, 13, 17]. First, we define the Fisher score tensor $S$ associated to the output channel $P_{\rm out}$ and its Fisher information $\Delta$ as

[TABLE]

where it is understood in (5) that the function y\mapsto\frac{\partial\log{P_{\rm out}(y,w)}}{\partial w}\big{|}_{w=0} acts component-wise on $Y$. Informally speaking, the channel universality property states that the mutual information of the problem defined by the output channel $P_{\rm out}$ is the same as the one of a AWGN (1) with variance $\Delta$, and that the AMP algorithm written for the Bayes-optimal inference of low-rank tensors then depends on the data tensor $Y$ and the output channel $P_{\rm out}$ only trough the tensor $S$ and the effective noise $\Delta$.

AMP involves an auxiliary function that depends explicitly on the prior as follows. Define the probability distribution

[TABLE]

where ${\cal Z}_{X}(A,B)$ is a normalization factor. Then AMP uses the function ${f_{\rm in}}(A,B)\in\mathbb{R}^{r},A\in\mathbb{R}^{r\times r},B\in\mathbb{R}^{r}$ defined by the expectation $f_{\rm in}(A,B)={\mathbb{E}}_{{\cal M}(x)}[x]$ as well as the covariance matrix $\partial_{B}{f_{\rm in}}(A,B)$. We shall denote the overlap of $u=(u_{1},\dots u_{N}),v=(v_{1},\dots v_{N})\in(\mathbb{R}^{r})^{N}$ by

[TABLE]

AMP is then written as an iterative update procedure on the estimates of the posterior means and co-variances $\hat{x}_{i}$ and $\sigma_{i}$ that uses auxiliary variables $B_{i}\in\mathbb{R}^{r}$ and $A\in\mathbb{R}^{r\times r}$:

[TABLE]

where $\circ$ denotes a component-wise (Hadamard) product of matrices, and $x^{\circ p}$ the corresponding component-wise power.

4 Theoretical analysis

4.1 State evolution of AMP

The evolution of the AMP algorithm in the limit of large systems $N\to\infty$ can be tracked via a low-dimensional set of equations called the state evolution (SE). For maximum-likelihood estimation the state evolution have been used in [24]. Its heuristic derivation for the present case of general rank $r$, prior $P_{X}$, and output $P_{\rm out}$ follows line by line the matrix estimation case detailed in [17].

For the Bayes-optimal inference, SE is written in terms of an order parameter $M^{t}\in{\mathbb{R}}^{r\times r}$ describing the overlap between $\hat{x}^{t}$ (the AMP estimator at iteration $t$) and the ground truth $x^{0}$ defined as $M^{t}=\hat{x}^{t}\cdot x^{0}$, and reads

[TABLE]

where $Z\sim\mathcal{N}(0,I_{r})$ and $x_{0}\sim P_{X}$ are independent. $M^{\circ(n)}$ is again the $n$-th Hadamard power of a matrix $M$.

We shall not present a rigorous proof of the SE for tensor estimation and rely instead on standard arguments from statistical physics. The performance of the AMP algorithm can be understood by initializing the SE at $M^{t=0}=0$. Or when $M=0$ is a fixed point of SE we initialize as $M^{t=0}=\epsilon$, an infinitesimally small number (accounting for the fact that a random initialization of AMP will —due to finite size fluctuations— be infinitesimally correlated with the ground truth). We denote $M_{\rm AMP}$ the fixed point of the state evolution resulting from iterations of (12-13) from this initialization. The mean-squared error achieved by tensor-AMP is then

[TABLE]

where $\Sigma_{X}=\mathbb{E}_{x}[xx^{\top}]$. When $P_{X}$ has zero mean, this is the covariance matrix of $P_{X}$.

4.2 Information-theoretically optimal inference

Our next goal is to analyze the performance of (possibly intractable) Bayes-optimal inference that evaluates the marginals of the posterior probability distribution (4). The error achieved by this procedure will be denoted the minimum mean-squared error (MMSE) and is formally defined as

[TABLE]

where the infimum is taken over all measurable functions $\hat{\theta}$ of $Y$. In order to compute the MMSE it is instrumental to compute the mutual information $I(X^{0};Y)$. This quantity is related to the free energy from statistical physics (see section 5 and [13]). To compute the limit of such quantities, one traditionally applies the replica method stemming from statistical physics [19]. We take advantage of the fact that for the Bayes-optimal inference the so-called replica symmetric version of this method yields the correct free energy [28]. The replica method yields

[TABLE]

where $\widehat{M}={M}^{\circ(p-1)}/{\Delta}$, ${\cal Z}_{X}(A,B)$ is defined in eq. (7), $x_{0}\sim P_{X}$ and $Z\sim\mathcal{N}(0,I_{r})$ are independent random variables. $S_{r}^{+}$ denotes the set of $r\times r$ symmetric positive semi-definite matrices. In section 5 we prove this result for the rank-one case ($r=1$).

The replica free energy (16) not only provides the limit of the mutual information $I(X^{0};Y)$, but thanks to an “I-MMSE Theorem” (similar to [11]) it yields the value of the MMSE for tensor estimation, see sec. 5. Denoting $M^{*}={\rm argmax}_{M}\phi_{\rm RS}(M)$ we get

[TABLE]

We proved (17) rigorously, but only in the rank-one case and for odd values of $p$, see again sec. 5. Notice that when $r\geq 2$ the estimation problem is symmetric under permutations of the $r$ columns of $X^{0}$: (17) is not expected to be true without further assumptions.

4.3 Statistical and computational trade-off

By evaluation of the derivative of (16) with respect to $M$ one can check that critical points of (16) are fixed points of the state evolution equations (12-13) allowing all the results to be read of the curve $\phi_{\rm RS}(M)$: The global maximum of (16) gives the MMSE while the (possibly local) maximum reached by iteration of (12-13) from the uninformative initialization yields the ${\rm MSE}_{\rm AMP}$.

We now discuss the interplay between the MMSE and ${\rm MSE}_{\rm AMP}$. The working hypothesis in this paper is that AMP yields lowest MSE among known polynomial algorithms. Depending on the parameters of model (4), i.e. the order of the tensor $p$, rank $r$, prior distribution $P_{X}$, and output channel $P_{\rm out}$ that appears in the SE only via its Fisher information $\Delta$, we can distinguish between two cases: the easy phase where asymptotically AMP is Bayes optimal so that ${\rm MMSE}={\rm MSE}_{\rm AMP}$, and the hard phase where ${\rm MMSE}<{\rm MSE}_{\rm AMP}$.

Given both the ${\rm MMSE}$ and ${\rm MSE}_{\rm AMP}$ are non-decreasing in $\Delta$ we denote the borders of the hard phase (when it exists) as follows: Information theoretic threshold $\Delta_{\rm IT}$ as the (limsup of the) highest $\Delta$ for which ${\rm MMSE}<{\rm MSE}_{\rm AMP}$. Algorithmic threshold $\Delta_{\rm Alg}$ as the (liminf of the) lowest $\Delta$ for which ${\rm MMSE}<{\rm MSE}_{\rm AMP}$. Another threshold used in this paper is that of a critical value $\Delta_{c}$ defined as smallest $\Delta$ such that for $\Delta>\Delta_{c}$ one has $M_{\rm AMP}=M^{*}(\Delta=+\infty)$ (the estimate one can do when the noise is infinite), and for $\Delta<\Delta_{c}$ one has $M_{\rm AMP}>M^{*}(\Delta=+\infty)$. Note that from the definition we must have $\Delta_{c}\geq\Delta_{\rm Alg}$. In cases where the hard phase does not exist, but $\Delta_{c}<\infty$ we will consider that $\Delta_{c}=\Delta_{\rm IT}=\Delta_{\rm Alg}$.

Existing results on maximum likelihood estimation [24] suggest that for tensor decomposition $p\geq 3$ we have $\Delta_{\rm Alg}=\Delta_{c}=0$ in the limit $N\to 0$ considered in this paper. This means that the spiked model of low-rank tensor decomposition is algorithmically very hard, compared to matrix $p=2$ case. The authors of [24] give a good account on how $\Delta$ needs to scale with $N$ for known polynomial algorithms to work.

For the Bayes-optimal estimation the situation seems at first sight similar. Indeed, whenever the prior $P_{X}$ has a zero mean, for $p\geq 3$ we get $\Delta_{\rm Alg}=\Delta_{c}=0$ and the hard phase is consequently huge. This can be seen as follows. Indeed if the mean of the prior $P_{X}$ is zero then the state evolution equations (12-13) have a fixed point $M=0$. Expanding the state evolution around this fixed point we find

[TABLE]

Whenever $p\geq 3$ the fixed point $M=0$ is stable for all $\Delta>0$. Hence $\Delta_{\rm Alg}=\Delta_{c}=0$ for priors of zero mean.

A closer look, however, shows that the situation is not so pessimistic. Indeed, as soon as the mean of the prior $P_{X}$ is non-zero, $M=0$ is no longer a fixed point of the state evolution and once we solve the state evolution equations we observe either $\Delta_{\rm Alg}>0$ (with AMP performing optimally for $\Delta<\Delta_{\rm Alg}$) or the hard phase is completely absent and AMP has information-theoretically optimal performance for all $\Delta$. We give examples of such priors in section 6.

5 Rigorous results

We present in this section rigorous results for the rank-one case ($r=1$). As mentioned above, the universality property [15, 13] reduces the computation of the mutual information to the case of additive white Gaussian noise.

Consider a probability distribution $P_{X}$ over $\mathbb{R}$ that admits a finite second moment $\Sigma_{X}$. The observation model (1) reduces in the rank-one case to

[TABLE]

where X^{0}=(x_{1}^{0},\dots,x_{N}^{0})\overset{\text{\tiny i.i.d.}}{\mathrel{\raisebox{-2.0pt}{\sim}}}P_{X} and (V_{i_{1},\dots,i_{p}})_{i_{1}<\dots<i_{p}}\overset{\text{\tiny i.i.d.}}{\mathrel{\raisebox{-2.0pt}{\sim}}}\mathcal{N}(0,\Delta) are independent. We define the Hamiltonian

[TABLE]

for $X=(x_{1},\dots,x_{N})\in\mathbb{R}^{N}$. We also write $dP_{X}(X)=\prod_{i=1}^{N}dP_{X}(x_{i})$. The posterior distribution (4) of $X^{0}$ given $Y$ reads then:

[TABLE]

where ${\mathcal{Z}}_{N}$ is the appropriate normalizing factor. Then the free energy is defined as (minus) the logarithm of ${\mathcal{Z}}_{N}$ of the Boltzmann probability divided by $N$ and averaged over $Y$. This is of particular interest since it is related to the mutual information (see [13]):

[TABLE]

In the rank-one case, the expression (16) of $\phi_{\rm RS}$ simplifies, so that we will use in this section

[TABLE]

where $\mathbb{E}$ is the expectation with respect to the independent random variables $x_{0}\sim P_{X}$ and $Z\sim\mathcal{N}(0,1)$. The proof of (15) reduces then to the following Theorem.

Theorem 1** (Replica-symmetric formula for the free energy).**

Let $P_{X}$ be a probability distribution over  $\mathbb{R}$, with finite second moment. Then, for all $\Delta>0$

[TABLE]

We now define the tensor-MMSE, $\text{T-MMSE}_{N}$ by

[TABLE]

where the infimum is taken over all measurable functions $\hat{\theta}$ of the observations $Y$.

Let us write $\lambda=\frac{1}{\Delta}$. Using an “I-MMSE Theorem” (see [11]) and the fact that the tensor MMSE is achieved by the posterior mean of $(X^{0})^{\otimes p}$ given $Y$, it is not difficult to verify that

[TABLE]

The arguments are the same than in the matrix ($p=2$) case, see [14] Corollary 17. $\text{T-MMSE}(\Delta)$ increases with the noise level $\Delta$, so that $\frac{\partial}{\partial\lambda}F_{N}$ is a non-decreasing function of $\lambda$. $F_{N}$ is thus a convex function of $\lambda$, and so is $F_{\rm RS}$ its pointwise limit. Consequently, $\frac{\partial}{\partial\lambda}F_{N}\to\frac{\partial}{\partial\lambda}F_{\rm RS}$ at all values of $\lambda$ at which $F_{\rm RS}$ is differentiable, that is for almost every $\Delta>0$. For these values of $\Delta$, one can also verify that the maximizer $m^{*}$ of $\phi_{\rm RS}$ is unique: we refer to [14] for a detailed proof in the matrix case $p=2$. We thus obtain the following theorem:

Theorem 2** (Tensor-MMSE).**

For almost every $\Delta>0$, $\phi_{\rm RS}$ admits a unique maximizer $m^{*}(\Delta)$ over $\mathbb{R}_{+}$ and

[TABLE]

The information-theoretic threshold $\Delta_{\rm IT}$ is the maximal value of $\Delta$ such that ${\lim\text{T-MMSE}_{N}<\Sigma_{X}^{p}-\mathbb{E}_{P_{X}}[x]^{2p}}$ (which is the asymptotic performance achieved by random guess). We obtain thus the precise location of the information-theoretic threshold:

[TABLE]

Let $X=(x_{1},\dots,x_{N})$ be a sample from the posterior (4), independently of everything else. An extension of Theorem 2 of [14] (that was derived for priors $P_{X}$ with bounded support) to the tensor case, gives that for almost every $\Delta>0$,

[TABLE]

i.e. the $p^{\text{th}}$-power of the overlap $X\cdot X^{0}$ concentrates around $m^{*}$. This leads to

Theorem 3** (Vector-MMSE for odd $p$).**

Suppose that $P_{X}$ has a bounded support. If $p$ is odd, then for almost every $\Delta>0$

[TABLE]

Before showing how (23) implies Theorem 3 we need to introduce a fundamental property of Bayesian inference: the Nishimori identity.

Proposition 1** (Nishimori identity).**

Let $(X,Y)$ be a couple of random variables on a polish space. Let $k\geq 1$ and let $X^{(1)},\dots,X^{(k)}$ be $k$ i.i.d. samples (given $Y$) from the distribution $P(X=\cdot\,|\,Y)$, independently of every other random variables. Let us denote $\langle\cdot\rangle$ the expectation with respect to $P(X=\cdot\,|\,Y)$ and $\mathbb{E}$ the expectation with respect to $(X,Y)$. Then, for all continuous bounded function $f$

[TABLE]

Proof.

It is equivalent to sample the couple $(X,Y)$ according to its joint distribution or to sample first $Y$ according to its marginal distribution and then to sample $X$ conditionally to $Y$ from its conditional distribution $P(X=\cdot\,|\,Y)$. Thus the $(k+1)$-tuple $(Y,X^{(1)},\dots,X^{(k)})$ is equal in law to $(Y,X^{(1)},\dots,X^{(k-1)},X)$. ∎

We will now use Proposition 1 to prove Theorem 3.

Proof of Theorem 3.

Let $\langle\cdot\rangle$ denote the expectation with respect to the posterior distribution $P(X^{0}=\cdot\ |\ Y)$, and let $X$ be a sample from this distribution, independently of everything else. The best estimator of $X^{0}$ in term of mean-squared error is the posterior mean $\langle X\rangle=(\langle x_{1}\rangle,\dots,\langle x_{N}\rangle)$. Therefore

[TABLE]

where $X^{\prime}$ is another sample from $\langle\cdot\rangle$, independently of everything else. We apply now the Nishimori identity (Proposition 1) to obtain $\mathbb{E}\langle X\cdot X^{\prime}\rangle=\mathbb{E}\langle X^{0}\cdot X\rangle$. This gives

[TABLE]

We then deduce from (23) that $\mathbb{E}\langle X\cdot X^{0}\rangle\xrightarrow[N\to\infty]{}m^{*}$, because $p$ is here supposed to be odd. This concludes the proof. ∎

We will now prove Theorem 1. For the matrix case ($p=2$), this has been proved in [13, 4, 14] and we explain here how this can be adapted to the case $p\geq 2$. To prove the limit (22), one shows successively an upper bound on $\limsup F_{N}$ and the matching lower bound on $\liminf F_{N}$. As shown in [14] (Section 6.2.2) one only need to prove Theorem 1 for input distributions $P_{X}$ with finite support $S$. We now assume to be in this situation.

5.1 Adding a small perturbation

One of the key ingredient of the proof is the introduction of a small perturbation of our model, that takes the form of a small amount of side information. This kind of techniques are frequently used for the study of spin glasses, where these small perturbations forces the Gibbs measure to verify some crucial identities, see [21]. In our context of Bayesian inference, we will see that small quantities of side information “breaks” the correlations of the signal variables under the posterior distribution.

Let us fix $\epsilon\in[0,1]$, and suppose we have access to the additional information, for $1\leq i\leq N$

[TABLE]

where L_{i}\overset{\text{\tiny i.i.d.}}{\mathrel{\raisebox{-2.0pt}{\sim}}}\text{Ber}(\epsilon) and $*$ is a symbol that does not belong to $\mathbb{R}$. The posterior distribution of $X$ is now

[TABLE]

where $\mathcal{Z}_{N,\epsilon}$ is the appropriate normalization constant. For $X=(x_{1},\dots,x_{N})\in\mathbb{R}^{N}$ we will use the notation

[TABLE]

$\bar{X}$ is thus obtained by replacing the coordinates of $X$ that are revealed by $Y^{\prime}$ by their revealed values. The notation $\bar{X}$ will allow us to obtain a convenient expression for the free energy of the perturbed model

[TABLE]

The next lemma shows that the perturbation does not change the free energy up to the order of $\epsilon$. Recall that we supposed the support $S$ of $P_{X}$ to be finite, so we can find a constant $K$ such that $S\subset[-K,K]$.

Lemma 1**.**

For all $n\geq 1$ and all $\epsilon,\epsilon^{\prime}\in[0,1]$, we have

[TABLE]

Lemma 1 follows from a direct adaptation of Proposition 23 from [14] to the tensor case. Consequently, if we suppose $\epsilon\sim\mathcal{U}([0,1])$ and define $\epsilon_{N}=N^{-1/2}\epsilon$ and L_{i}\overset{\text{\tiny i.i.d.}}{\mathrel{\raisebox{-2.0pt}{\sim}}}\text{Ber}(\epsilon_{N}), independently of everything else, we have

[TABLE]

where $\mathbb{E}_{\epsilon}$ denotes the expectation with respect to $\epsilon$ only. It remains therefore to compute the limit of the free energy under a small perturbation. As shown in [20], the perturbation (24) forces the correlations to vanish asymptotically.

Lemma 2** (Lemma 3.1 from [20]).**

[TABLE]

Let us write $\langle\cdot\rangle$ the expectation with respect to $P(X=\cdot\ |\ Y,Y^{\prime})$, and let $X^{(1)},X^{(2)}$ be two independents samples from $P(X=\cdot\ |\ Y,Y^{\prime})$, independently of everything else. We define $Q=\langle X^{(1)}\cdot X^{(2)}\rangle$. Notice that $Q$ is a non-negative random variable. As a consequence of Lemma 2, the overlaps under the posterior distribution concentrates around $Q$:

Lemma 3**.**

[TABLE]

where $\mathbb{E}$ denotes the expectation with respect all random variables.

Lemma 3 follows from the arguments of section 4.4 from [14].

The arguments presented in this section are robust and apply to a large class of Hamiltonians. In particular, we will be able to apply in the sequel Lemmas 1 and 3 to other Hamiltonians and posterior distributions (and corresponding free energies).

5.2 Guerra’s interpolation scheme

The lower bound is obtained by extending the bound derived for $p=2$ in [13], using a Guerra-type interpolation [10] as was already done for tensors by Korada and Macris in [12] (who consider tensors in the special case of Rademacher $P_{X}$).

Lemma 4**.**

[TABLE]

Proof.

We use a Guerra-type interpolation [10]: Let $0\leq t\leq 1$ and $m\in\mathbb{R}_{+}$. We suppose to observe $Y$ and $\tilde{Y}$ given by

[TABLE]

where the variables $V_{i_{1},\dots,i_{p}}$ and $\tilde{V}_{j}$ are i.i.d. $\mathcal{N}(0,\Delta)$ random variables. We define the interpolating Hamiltonian

[TABLE]

Then, the posterior distribution of $X^{0}$ given $Y$ and $\tilde{Y}$ reads

[TABLE]

where $\mathcal{Z}_{N,t}$ is the appropriate normalization. Let $\psi_{N}(t)=\frac{1}{N}\mathbb{E}[\log\mathcal{Z}_{N,t}]$ be the corresponding free energy. Notice that

[TABLE]

Let $\langle\cdot\rangle_{t}$ denote the expectation with respect to the posterior (28) and let $X$ be a sampled from (28), independently of everything else.

Using Gaussian integration by parts and the Nishimori identity of Proposition 1 one can show (see [13, 14]) that for all $0\leq t\leq 1$

[TABLE]

By convexity of the function $a\mapsto a^{p}$ on $\mathbb{R}_{+}$ we have, for all $a,b\geq 0$: $a^{p}-pb^{p-1}a\geq(1-p)b^{p}$. We would like to use this inequality with $a=X\cdot X^{0}$ and $b=m$ to obtain that $\psi_{N}^{\prime}(t)\geq\frac{(1-p)m^{p}}{2\Delta p}$. This would conclude the proof of the lower bound because

[TABLE]

However, we do not know that $X\cdot X^{0}\geq 0$ almost surely. To bypass this issue we can add, as in sec. 5.1, a small perturbation (24) that forces $X\cdot X^{0}$ concentrates around a non-negative value (Lemma 3), without affecting the “interpolating free energy” $\psi_{N}(t)$ in the $N\to\infty$ limit, see (26). The arguments are the same than in sec. 5.1, so we omit the details and the rewriting of the previous calculations with the perturbation term. This concludes the proof. ∎

5.3 Proving the upper-bound: Aizenman-Sims-Starr scheme

We are now going to show how the arguments of [14] for the upper bound —using cavity computations with an Aizenman-Sims-Starr approach [1]— can be extended to the tensor case.

Lemma 5**.**

[TABLE]

Proof.

We are going to compare the system with $N$ variables to the system with $N+1$ variables. Define $A_{N}=\mathbb{E}[\log{\mathcal{Z}}_{N+1}]-\mathbb{E}[\log{\mathcal{Z}}_{N}]$. Consequently, $F_{N}=\frac{1}{N}\sum_{k=0}^{N-1}A_{k}$ and $\limsup F_{N}\leq\limsup A_{N}$.

We are thus going to upper-bound $A_{N}$. Let $X\in S^{N}$ be the $N$-first variables and $\sigma\in S$ the $(N+1)^{\text{th}}$ variable. We decompose $H_{N+1}(X,\sigma)=H_{N}^{\prime}(X)+\sigma z(X)+\sigma^{2}s(X)$ where

[TABLE]

[TABLE]

One can also decompose $H_{N}(X)=H_{N}^{\prime}(X)+y(X)$ in law, where

[TABLE]

In the above definition, the $V^{\prime}$ are i.i.d. $\mathcal{N}(0,\Delta)$ random variables, independent of everything else, and $r_{n}=o(N^{-p})$. If we denote by $\langle\cdot\rangle^{\prime}$ the Gibbs measure on $S^{N}$ corresponding to the Hamiltonian $\log P_{X}+H_{N}^{\prime}$ we can rewrite

[TABLE]

where $X$ is a sample from $\langle\cdot\rangle^{\prime}$, independently of everything else. $A_{N}$ is thus a difference of two terms that will correspond exactly to the terms of (21). As in sec. 5.1, we can show that under a small perturbation of the system, the overlap $X\cdot X^{0}$ with the planted configuration concentrates around a non-negative value $Q^{\prime}$. This leads to simplifications in (29):

[TABLE]

For a precise derivation of (30), the reader is invited to report to the matrix case (see [14], sec. 4.6), since there is no major difference with the tensor case on this point. The arguments presented there are commonly used in the study of spin glasses and are the analog of cavity computations in the SK model developed in [26], sec. 1.5. This concludes the proof. ∎

6 Examples of phase transitions

We used the state evolution eqs. (12-13), and the free energy (16), to compute the values of the thresholds $\Delta_{c}$, $\Delta_{\rm IT}$ and $\Delta_{\rm Alg}$ for several examples of the prior distributions: Gaussian (spherical spins), $P_{X}(x)={\cal N}\left(\mu,1\right)$; Rademacher (Ising spins), $P_{X}(x)=\frac{1}{2}\left[\delta(x-1)+\delta(x+1)\right]$; Bernoulli (sub-tensor localization), $P_{X}(x)=\rho\delta(x-1)+(1-\rho)\delta(x)$; and clustering (tensor stochastic block model), $P_{X}(x)=\frac{1}{r}\sum\limits_{k=1}^{r}\delta(x-\vec{e}_{k})$, where $\vec{e}_{k}\in\mathbb{R}^{r}$ is a vector with a 1 at coordinate $k$ and 0 elsewhere. Examples of values of the thresholds for the above priors are given in Table 1. For the zero mean Gaussian and the Rademacher prior our results for $\Delta_{\rm IT}$ indeed agree with those presented in [12, 22]. Central and right part of Fig. 1 present the thresholds for the Gaussian and Bernoulli prior as a function of the mean $\mu$ and density $\rho$, respectively. Left part of Fig. 1 illustrates that indeed the fixed points of the state evolution agree with the fixed points of the AMP algorithm.

6.1 Results for Gaussian prior

In this section we detail the analysis of the state evolution for rank $r=1$ Gaussian prior of mean $\mu$ and variance 1.

[TABLE]

Using (12) one gets for the SE equation

[TABLE]

where $M$ is a scalar, and $\Delta$ is the inverse Fisher information of the output channel. It turns out that as soon as $p\geq 3$ the SE equation exhibits multiple stable fixed points.

For the zero mean $\mu=0$ case one gets

[TABLE]

Here the fixed point $M=0$ is stable whatever the noise $\Delta>0$ and therefore AMP will not achieve performance better than random guessing for any $\Delta>0$. Ref. [24] studies the scaling of $\Delta$ with $N$ for which AMP and other algorithms succeed.

For positive mean $\mu>0$, however, the AMP algorithm is able to recover the signal for values of $\Delta<\Delta_{\rm Alg}$ with

[TABLE]

where we defined a new threshold $\Delta_{\rm Dyn}$ as the smallest such that for $\Delta>\Delta_{\rm Dyn}$ the state evolution has a unique fixed point. We know of no analytical formula for $\Delta_{\rm IT}$ and for Figure 1 we computed it numerically. The tri-critical point where all these curve meet is located at

[TABLE]

Using the above expressions we derive that

[TABLE]

We can also compute the limit of the $\Delta_{\rm IT}(\mu=0,p)$ as $p\rightarrow\infty$ and get

[TABLE]

This scaling agrees with the large $p$ behavior derived in [24] and [22].

6.2 Results for clustering prior

An interesting example of the prior for rank $r$ tensor estimation is

[TABLE]

This describes a model of $r$ non-overlapping clusters. Due to the channel universality, this prior also describes the stochastic block model on dense hyper-graphs as considered for sparse hyper-graph in e.g. [3]. This model was considered in detail for $p=2$ in [17].

The above clustering prior has non-zero mean, and it also exhibits the transition $\Delta_{c}$ from a phase where recovery of clusters better than chance is not possible, to a phase where it is.

To analyze the SE equations we first notice that the stable fixed point will be of the form

[TABLE]

where $I_{r}$ is the identity matrix and $J_{r}$ is a matrix filled with ones. $b=0$ means that the estimate of the marginals does not carry any information. $b=1$ means perfect reconstruction. The state evolution now becomes

[TABLE]

where ${\cal M}_{r}$ is a function that was defined and studied in [15]. Its Taylor expansion is

[TABLE]

We further notice that $b=0$ is always a fixed point of (43). By expanding (43) to first order one gets

[TABLE]

This fixed point therefore becomes unstable when

[TABLE]

By analyzing eq. (43) further we can prove that $\forall x\in\mathbb{R}^{+}$

[TABLE]

[TABLE]

$m$ is a fixed point of (43) when $m=m(x)$ and $\Delta=\Delta(x)$. Rather than finding the fixed point iteratively, the above equations allow us to draw all the fixed point of (43), be it stable or unstable We have that $m(x)$ is a stable fixed point of (43) if and only if

[TABLE]

The next question is whether there is a first or second order phase transition at $\Delta_{c}$. To answer this, one needs to analyze whether the fixed point close to $b=0$ is stable or unstable. For this we to compute $\frac{\partial\Delta(x)}{\partial x}$ at $x=0$ to get using (44) that

[TABLE]

Therefore if $-2p-r+pr>0$ there will be no stable fixed point close to $b=0$ and the system must have a first order phase transition (discontinuity in the ${\rm MSE}_{\rm AMP}$) at $\Delta_{c}=\Delta_{\rm Alg}$.

For two clusters $r=2$, there is a second order phase transition at $\Delta_{c}$ for all $p\geq 2$. However, analyzing the state evolution numerically we observed that for $p\geq 5$ there is a discontinuity later at some $\Delta_{\rm Alg}<\Delta_{c}$. For three and more clusters $r\geq 3$ we always have $\Delta_{\rm Alg}=\Delta_{c}$, and for $-2p-r+pr\leq 0$ we have not detected any other discontinuities. Values of $\Delta_{\rm IT}$ for three clusters and some values of $p$ are given in Table 1.

Acknowledgment

This work has been supported by the ERC under the European Union’s FP7 Grant Agreement 307087-SPARCS.