{'fraction_non_alphanumeric': 0.036143651297439745, 'fraction_numerical': 0.00902148614076596, 'mean_word_length': 4.345260127493317, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Statistical mechanics of low-rank tensor decomposition

Often, large, high dimensional datasets collected across multiple modalities can be organized as a higher order tensor. Low-rank tensor decomposition then arises as a powerful and widely used tool to discover simple low dimensional structures underlying such data. However, we currently lack a theoretical understanding of the algorithmic behavior of low-rank tensor decompositions. We derive Bayesian approximate message passing (AMP) algorithms for recovering arbitrarily shaped low-rank tensors buried within noise, and we employ dynamic mean field theory to precisely characterize their performance. Our theory reveals the existence of phase transitions between easy, hard and impossible inference regimes, and displays an excellent match with simulations. Moreover, it reveals several qualitative surprises compared to the behavior of symmetric, cubic tensor decomposition. Finally, we compare our AMP algorithm to the most commonly used algorithm, alternating least squares (ALS), and demonstrate that AMP significantly outperforms ALS in the presence of noise.

AMP algorithms have a long history dating back to early work on the statistical physics of perceptron learning [7,8] (see [9] for a review). The term AMP was coined by Donoho, Maleki and Montanari in their work on compressed sensing [10] (see also [11,12,13,14,15,16] for replica approaches to compressed sensing and high dimensional regression). AMP approximates belief propagation in graphical models and a rigorous analysis of AMP was carried out in [17]. For a rank-one matrix estimation problem, AMP was first introduced and analyzed in [18]. This framework has been extended in a beautiful body of work by Krzakla and Zdeborova and collaborators to various low-rank matrix factorization problems in [19,20,21,22]). Also, recently low-rank tensor decomposition through AMP was studied in [21], but their analysis was limited to symmetric tensors which are then necessarily cubic in shape.
However, tensors that occur naturally in the wild are almost never cubic in shape, nor are they symmetric. The reason is that the p different modes of an order p tensor correspond to measurements across very different modalities, resulting in very different numbers of dimensions across modes, yielding highly irregularly shaped, non-cubic tensors with no symmetry properties. For example in EEG studies 3 different tensor modes could correspond to time, spatial scale, and electrodes [23]. In fMRI studies the modes could span channels, time, and patients [24]. In neurophysiological measurements they could span neurons, time, and conditions [25] or neurons, time, and trials [3]. In studies of visual cortex, modes could span neurons, time and stimuli [26].
Thus, given that tensors in the wild are almost never cubic, nor symmetric, to bridge the gap between theory and experiment, we go beyond prior work to derive and analyze Bayes optimal AMP algorithms for arbitrarily shaped high order and low rank tensor decomposition with different priors for different tensor modes, reflecting their different measurement types. We find that the low-rank decomposition problem admits two phase transitions separating three qualitatively different inference regimes: (1) the easy regime at low noise where AMP works, (2) the hard regime at intermediate noise where AMP fails but the ground truth tensor is still possible to recover, if not in a computationally tractable manner, and (3) the impossible regime at high noise where it is believed no algorithm can recover the ground-truth low rank tensor.
From a theoretical perspective, our analysis reveals several surprises relative to the analysis of symmetric cubic tensors in [21]. First, for symmetric tensors, it was shown that the easy inference regime cannot exist, unless the prior over the low rank factor has non-zero mean. In contrast, for non-symmetric tensors, one tensor mode can have zero mean without destroying the existence of the easy regime, as long as the other modes have non-zero mean. Furthermore, we find that in the space of all possible tensor shapes, the hard regime has the largest width along the noise axis when the shape is cubic, thereby indicating that tensor shape can have a strong effect on inference performance, and that cubic tensors have highly non-generic properties in the space of all possible tensor shapes.
Before continuing, we note some connections to the statistical mechanics literature. Indeed, AMP is closely equivalent to the TAP equations and the cavity method [27,28] in glassy spin systems. Furthermore, the posterior distribution of noisy tensor factorization is equivalent to p-spin magnetic systems [29], as we show below in section 2.2. For Bayes-optimal inference, the phase space of the problem is reduced to the Nishimori line [30]. This ensures that the system does not exhibit replicasymmetry breaking. Working in the Bayes-optimal setting thus significantly simplifies the statistical analysis of the model. Furthermore, it allows theoretical insights into the inference phase-transitions, as we shall see below. In practice, for many applications the prior or underlying rank of the tensors are not known a-priori. The algorithms we present here can also be applied in a non Bayes-optimal setting, where the parametric from of the prior can not be determined. In that case, the theoretical asymptotics we describe here may not hold. However, approximate Bayesian-optimal settings can be recovered through parameter learning using expectation-maximization algorithms [31]. We discuss these consequences in section 4. Importantly,the connection to the statistical physics of magnetic systems allows the adaptation of many tools and intuitions developed extensively in the past few decades, see e.g. [32]. We discuss more connections to statistical mechanics as we proceed below.

Low rank decomposition using approximate message passing
In the following we define the low-rank tensor decomposition problem and present a derivation of AMP algorithms designed to solve this problem, as well as a dynamical mean field theory analysis of their performance. A full account of the derivations can be found in the supplementary material.

Low-rank tensor decomposition
Consider a general tensor Y of order-p, whose components are given by a set of p indices, Y i1,i2,...,ip . Each index i α is associated with a specific mode of the tensor. The dimension of the mode α is N α so the index i α ranges from 1, . . . , N α . If N α = N for all α then the tensor is said to be cubic. Otherwise we define N as the geometric mean of all dimensions N = ( p α N α ) 1/p , and denote n α ≡ N α /N so that p α n α = 1. We employ the shorthand notation Y i1,i2,...,ip ≡ Y a , where a = {i 1 , . . . , i p } is a set of p numbers indicating a specific element of Y . A rank-1 tensor of order-p is the outer product of p vectors (order-1 tensors) ⊗ 1≤α≤p x α ,where x α ∈ R Nα . A rank-r tensor of order-p has a special structure that allows it to be decomposed into a sum of r rank-1 tensors, each of order-p. The goal of the rank decomposition is to find all x ρ α ∈ R Nα , for α = 1, ..p, and ρ = 1, . . . , r, given a tensor Y of order-p and rank-r. In the following, we will use x αi ∈ R r to denote the vector of values at each entry of the tensor, spanning the r rank-1 components. In a low-rank decomposition it is assumed that r < N . In noisy low-rank decomposition, individual elements Y a are noisy measurements of a low-rank tensor [ Figure 1.A]. A comprehensive review on tensor decomposition can be found in [4].
We state the problem of low-rank noisy tensor decomposition as follows: Given a rank-r tensor we would like to find all the underlying factors x ρ αi . We note that we have used the shorthand notation i = i α to refer to the index i α which ranges from 1 to N α , i.e. the dimensionality of mode α of the tensor. Now consider a noisy measurement of the rank-r tensor w given by where ǫ is a random noise tensor of the same shape as w whose elements are distributed i.i.d according to a standard normal distribution, yielding a total noise variance ∆ ∼ O(1) [ Fig. 1.A]. The underlying factors x ρ αi are sampled i.i.d from a prior distribution P α (x), that may vary between the modes α. This model is a generalization of the spiked-tensor models studied in [21,33].
We study the problem in the thermodynamic limit where N → ∞ while r,n α ∼ O(1). In that limit, the mean-field theory we derive below becomes exact. The achievable performance in the decomposition problem depends on the signal-to-noise ratio (SNR) between the underlying low rank tensor (the signal) and the noise variance ∆. In eq. (1) we have scaled the SNR (signal variance divided by noise variance) with N so that the SNR is proportional to the ratio between the O(N ) unknowns and the N p measurements, making the inference problem neither trivially easy nor always impossible. From a statistical physics perspective, this same scaling ensures that the posterior distribution over the factors given the data corresponds to a Boltzmann distribution whose Hamiltonian has extensive energy proportional to N , which is necessary for nontrivial phase transitions to occur.

Tensor decomposition as a Bayesian inference problem
In Bayesian inference, one wants to compute properties of the posterior distribution Here P out (Y a |w a ) is an element-wise output channel that introduce independent noise into individual measurements. For additive white Gaussian noise, the output channel in eq. (3) is given Y Y 112 x 21 x 11 x 32 x 2j x 11 x 12 x 13 x 14 x 21 x 22 x 23 x 24 x 31 x 32 by where g(·) is a quadratic cost function. The denominator Z(Y, w), in 3 is a normalization factor, or the partition function in statistical physics. In a Bayes-optimal setting, the priors P α (x), as well as the rank r and the noise ∆ are known.
The channel universality property [19] states that for low-rank decomposition problems, r ≪ N , any output channel is equivalent to simple additive white Gaussian noise, as defined in eq. (2). Briefly, the output channel can be developed as a power series in w a . For low-rank estimation problems we have w a ≪ 1 [eq. (1)], and we can keep only the leading terms in the expansion. One can show that the remaining terms are equivalent to random Gaussian noise, with variance equal to the inverse of the Fisher information of the channel [See supplementary material for further details]. Thus, nonadditive and non-Gaussian measurement noise at the level of individual elements, can be replaced with an effective additive Gaussian noise, making the theory developed here much more generally applicable to diverse noise scenarios.
The motivation behind the analysis below, is the observation that the posterior (3), with the quadratic cost function (4) is equivalent to a Boltzmann distribution of a magnetic system at equilibrium, where x αi ∈ R r can be though of as the r-dimensional vectors of a spherical (xy)-spin model [29].

Approximate message passing on factor graphs
To solve the problem of low-rank decomposition we frame the problem as a graphical model with an underlying bipartite factor graph. The variable nodes in the graph represent the rN α n α unknowns x ρ αi and the N p factor nodes correspond to the measurements Y a . The edges in the graph are between factor node Y a and the variable nodes in the neighbourhood ∂a [ Figure 1.B]. More precisely, for each factor node a = {i 1 , i 2 , ..., i p }, the set of variable nodes in the neighbourhood ∂a are precisely {x 1i1 , x 2i2 , ..., x pip }, where each x αiα ∈ R r . Again, in the following we will use the shorthand notation x αi for x αiα . The state of a variable node is defined as the marginal probability distribution η αi (x) for each of the r components of the vectors x αi ∈ R r . The estimatorsx αi ∈ R r for the values of the factors x αi are given by the means of each of the marginal distributions η αi (x).
In the approximate message passing framework, the state of each node (also known as a 'belief'), η αi (x) is transmitted to all other variable nodes via its adjacent factor nodes [ Fig. 1.C]. The state of each node is then updated by marginalizing over all the incoming messages, weighted by the cost function and observations in the factor nodes they passed on the way in: Here P α (x) is the prior for each factor x αi associated with mode α, and Z αi =´dxη αi (x) is the partition function for normalization. The first product in (5) spans all factor nodes adjacent to variable node αi. The second product is over all variable nodes adjacent to each of the factor nodes, excluding the target node αi. The trace T r x βj denotes the marginalization of the cost function g(y a , w a ) over all incoming distributions.
The mean of the marginalized posterior at node αi is given bŷ and its covariance isσ Eq. (5) defines an iterative process for updating the beliefs in the network. In what follows, we use mean-field arguments to derive iterative equations for the means and covariances of the these beliefs in (6)- (7). This is possible given the assumption that incoming messages into each node are probabilistically independent. Independence is a good assumption when short loops in the underlying graphical model can be neglected. One way this can occur is if the factor graph is sparse [34,35]. Such graphs can be approximated by a directed acyclic graph; in statistical physics this is known as the Bethe approximation [36]. Alternatively, in low-rank tensor decomposition, the statistical independence of incoming messages originates from weak pairwise interactions that scale as w ∼ N −(p−1)/2 . Loops correspond to higher order terms interaction terms, which become negligible in the thermodynamic limit [17,32].
Exploiting these weak interactions we construct an accurate mean-field theory for AMP. Each node αi receives N p−2 messages from every node βj with β = α, through all the factor nodes that are connected to both nodes,{y b |b ∈ ∂αi ∪ ∂βj} [ Fig. 1.D]. Under the independence assumption of incoming messages, we can use the central limit theorem to express the state of node αj in (5) as where A −1 β u βj and A −1 β are the mean and covariance of the local incoming messages respectively. The distribution is normalized by the partition function The mean and covariance of the distribution, eq. (6) and (7) are the moments of the partition function Finally, by expanding g(Y a , w a ) in eq. (5) to quadratic order in w, and averaging over the posterior, one can find a self consistent equation for A αi and u αi in terms of x αi and Y [see supplemental material for details].

AMP algorithms
Using equations (10), and the self-consistent equations for A αi and u αi , we construct an iterative algorithm whose dynamics converges to the solution of the self-consistent equations [ see supple-mental material for details]. The resulting update equations for the parameters are The second term on the RHS of (11), is given by This term originates from the the exclusion the target node αi from the product in equations (5) and (8). In statistical physics it corresponds to an Onsager reaction term due to the removal of the node yielding a cavity field [37]. In the above, the notations ⊙, ⊙ denote component-wise multiplication between two, and multiple tensors respectively.
Note that in the derivation of the iterative update equation above, we have implicitly used the assumption that we are in the Bayes-optimal regime which simplifies eq. (11)- (14) [see supplementary material for details]. The AMP algorithms can be derived without the assumption of Bayes-optimality, resulting in a slightly more complicated set of algorithms [See supplementary material for details]. However, further analytic analysis, which is the focus of this current work, and the derivation of the dynamic mean-field theory which we present below is applicable in the Bayes-optimal regime, were there is no replica-symmetry breaking, and the estimators are self-averaging. Once the update equations converge, the estimates for the factors x αi and their covariances are given by the fixed point value of equations (13) and (14) respectively. A statistical treatment for the convergence in typical settings is presented in the following section.

Dynamic mean-field theory
To study the performance of the algorithm defined by eq. (11)-(14), we use another mean-field approximation that estimates the evolution of the inference error. As before, the mean-field becomes exact in the thermodynamic limit. We begin by defining order parameters that measure the correlation of the estimatorsx t αi with the ground truth values x αi for each mode α of the tensor Technically, the algorithm is permutation invariant, so one should not expect the high correlation values to necessarily appear on the diagonal of M t α . In the following, we derive an update equation for M t α , which will describe the performance of the algorithm across iterations. An important property of Bayes-optimal inference is that there is no statistical difference between functions operating on the ground truth values, or on values sampled uniformly from the posterior distribution. In statistical physics this property is known as one of the Nishimori conditions [30]. These conditions allow us to derive a simple equation for the update of the order parameter (16). For example, from (12) one easily finds that in Bayes-optimal settings A t α =M t α . Furthermore, averaging the expression for u αi over the posterior, we find that [ supplemental material] Similarly, the covariance matrix of u αi under the posterior is Finally, using eq. (13) for the estimation ofx αi , and the definition of M t α in (16) we find a dynamical equation for the evolution of the order parameters M t α : where f α ≡ ∂ ∂u log Z α (A, u) is the estimation ofx t+1 αi from (13). The average in (20) is over the prior P α (x) and over the standard Gaussian variables z ∈ R r . The average over z represents fluctuations in the meanM t α x αi in (17), due to the covarianceM t α in (19). Finally, the performance of the algorithm is given by the fixed point of (20), As we will see below, the inference error can be calculated from the fixed point order parameters M * α in a straightforward manner.

Phase transitions in generic low-rank tensor decomposition
The dynamics of M t α depend on the SNR via the noise level ∆. To study this dependence, we solve equations (20) and (21) with specific priors. Below we present the solution of using Gaussian priors. In the supplementary material we also solve for Bernoulli and Gauss-Bernoulli distributions, and discuss mixed cases where each mode of the tensor is sampled from a different prior. Given our choice of scaling in (1) For simplicity we focus on the case of rank r = 1, where the order parameters M t α in (16) become scalars, which we denote m t α .

Solution with Gaussian priors
We study the case where x αi are sampled from normal distributions with mode-dependent mean and variance P α (x) ∼ N (µ α , σ 2 α ). The mean-field update equation (20) can be written as wherem t α ≡ ∆ −1 n α β =α m β , as in (18). We define the average inference error for all modes where m * α is the fixed point of eq. (22). Though we focus here on the r = 1 case for simplicity, the theory is equally applicable to higher-rank tensors.
Solutions to the theory in (22) and (23) are plotted in Fig. 2.A together with numerical simulations of the algorithm (11)-(14) for order-3 tensors generated randomly according to (2). The theory and simulations match perfectly. The AMP dynamics for general tensor decomposition is qualitatively similar to that of rank-1 symmetric matrix and tensor decompositions, despite the fact that such symmetric objects possess only one mode. As a consequence, the space of order parameters for these two problems is only one-dimensional; in contrast for the general case we consider here, it is p-dimensional. Indeed, the p = 3 order parameters are all simultaneously and correctly predicted by our theory.
For low levels of noise, the iterative dynamics converge to a stable fixed point of (22) with low MSE. As the noise increases beyond a bifurcation point ∆ alg , a second stable fixed point emerges with As either of the two nonzero means increases, ∆ dyn increases with them, reflecting an increase in the regime of noise level ∆ over which inference is easy. Importantly, the transition ∆ dyn is finite even when only one prior has non zero mean. (C) Same as (B) but for ∆ alg . Again, as either mean increases, ∆ alg increases also, reflecting a delay in the onset of the impossible regime as the noise level ∆ increases. The algorithmic phase transition is finite when at most one prior has zero mean. (D) Lower and higher transition points ∆ Alg (blue) and ∆ Dyn (orange) as a function of tensor shape. The ratios between the mode dimensions are n α = {1, n x , 1/n x }. The width of the bi-stable or hard inference regime is widest at the cubic point where n x = 1.
m * α ≪ 1 and M SE ≈1. Above this point AMP may not converge to the true factors. The basin of attraction of the two stable fixed points are separated by a p − 1 dimensional sub-manifold in the p-dimensional order parameter space of m α . If the initial values x 0 αi have sufficiently high overlap with the true factors x αi , then the AMP dynamics will converge to the low error fixed point; we refer to this as the informative initialization, as it requires prior knowledge about the true structure. For uninformative initializations, the dynamics will converge to the high error fixed point almost surely in the thermodynamic limit.
At a higher level of noise, ∆ dyn , another pitchfork bifurcation occurs and the high error fixed point becomes the only stable point. With noise levels ∆ above ∆ dyn , the dynamic mean field equations will always converge to a high error fixed point. In this regime AMP cannot overcome the high noise and inference is impossible.
From eq. (22), it can be easily checked that if the prior means µ α = 0, ∀α then the high error fixed point with m α = 0 is stable for any finite ∆. This implies that ∆ alg = 0, and there is no easy regime for inference, so AM P with uninformed initialization will never find the true solution. This difficulty was previously noted for the low-rank decomposition of symmetric tensors [21], and it was further shown there that the prior must be non-zero for the existence of an easy inference regime. However, for general tensors there is higher flexibility; one mode α can have a zero mean without destroying the existence of an easy regime. To show this we solved (22), with different prior means for different modes and we plot the phase boundaries in Fig. 2.B-C. For the p = 3 case, ∆ dyn is finite even if two of the priors have zero mean. Interestingly, the algorithmic transition ∆ alg is finite if at most one prior is has zero mean. Thus, the general tensor decomposition case is qualitatively different than the symmetric case in that an easy regime can exist even when a tensor mode has zero mean.

Non-cubic tensors
The shape of the tensor, defined by the different mode dimensions n α , has an interesting effect on the phase transition boundaries, which can be studied using (22). In figure 2.D the two transitions, ∆ alg and ∆ dyn are plotted as a function of the shape of the tensor. Over the space of all possible tensor shapes, the boundary between the hard and impossible regimes, ∆ dyn is maximized, or pushed furthest to the right in Fig. 2.D, when the shape takes the special cubic form where all dimensions are equal n α = 1, ∀α. This diminished size of the impossible regime at the cubic point can be understood by noting the cubic tensor has the highest ratio between the number of observed data points N p and the number of unknowns rN α n α .
Interestingly the algorithmic transition is lowest at this point. This means that although the ratio of observations to unknowns is the highest, algorithms may not converge, as the width of the hard  8 10 , 10 8 regime is maximized. To explain this observation, we note that in (18), the noise can be rescaled independently in each mode by defining ∆ → ∆ α = ∆/n α . It follows that for non-cubic tensors the worst case effective noise across modes will be necessarily higher than in the cubic case. As a consequence, moving from cubic to non-cubic tensors lowers the minimum noise level ∆ alg at which the uninformative solution is stable, thereby extending the hard regime to the left in Fig. 2.D.

Bayesian AMP compared to maximum a-posteriori (MAP) methods
We now compare the performance of AMP to one of the most commonly used algorithms in practice, namely alternating least squares (ALS) [4]. ALS is motivated by the observation that optimizing one mode while holding the rest fixed is a simple least-squares subproblem [6,5]. Typically, ALS performs well at low noise levels, but here we explore how well it compares to AMP at high noise levels, in the scaling regime defined by defined by (1) and (2), where inference can be non-trivial.
In Fig. 3 we compare the performance of ALS with that of AMP on the same underlying large (N = 500) tensors with varying amounts of noise. First, we note that that ALS does not exhibit a sharp phase transition, but rather a smooth cross-over between solvable and unsolvable regimes. Second, the robustness of ALS to noise is much lower than that of AMP. This difference is more substantial as the size of the tensors, N , is increased [data not shown].
One can understand the difference in performance by noting that ALS is like a MAP estimator, while Bayesian AMP attempts to find the minimal mean square error (MMSE) solution. AMP does so by marginalizing probabilities at every node. Thus AMP is expected to produce better inferences when the posterior distribution is rough and dominated by noise. From a statistical physics perspective, ALS is a zero-temperature method, and so it is subject to replica symmetry breaking. AMP on the other hand is Bayes-optimal and thus operates at the Nishimori temperature [30]. At this temperature the system does not exhibit replica symmetry breaking, and the true global ground state can be found in the easy regime, when ∆ < ∆ alg .

Summary
In summary, our work partially bridges the gap between theory and practice by creating new AMP algorithms that can flexibly assign different priors to different modes of a high-order tensor, thereby enabling AMP to handle arbitrarily shaped high order tensors that actually occur in the wild. Moreover, our theoretical analysis reveals interesting new phenomena governing how irregular tensor shapes can strongly affect inference performance and the positions of phase boundaries, and highlights the special, non-generic properties of cubic tensors. Finally, we hope the superior performance of our flexible AMP algorithms relative to ALS will promote the adoption of AMP in the wild. Code to reproduce all simulations presented in this paper is available at https://github.com/gangulilab/tensorAMP.
In the following sections, we derive the approximate message passing (AMP) algorithms for arbitrary low-rank tensors. The derivation follows similar lines as the derivation in [20,32] for the matrix p = 2 case, only here the model is generalized for higher modes p>2. We then derive dynamical mean field theory (also known as state evolution) and find the phase transition of the inference problem using non-linear analysis on the recursive dynamical equations of the order parameters. Lastly, we explicitly solve the equations for several simple examples of mode-3 tensors with a mixture of different prior distributions. For notational simplicity, we start the derivation assuming all modes are of the same dimension, N , which we assume to be in the thermodynamic limit, N → ∞. Then, we will generalize for noncubic tensors. We emphasize that even in the cubic case the tensor is non-symmetric, and each mode is independent and is iid drawn from a prior distribution, which is potentially different for each mode.
A Message passing on factorized graph

A.1 Factor graph for tensor decomposition -notations
We consider a given low-rank tensor Here the vectors x 0 αi ∈ R r denote the ground-truth values to the estimation problem. Each entry of the tensor is denoted with a lower-case latin letter {a, b, c, ...}. The notation stands for the set of p indices that define that tensor element, However, we only have access to noisy measurements of the ground-truth vectors, denoted by where ǫ a is a random tesnsor, whose elements are i.i.d. gaussians with zero mean and unit variance. We assume no covariance between two measurements, E(ǫ a ǫ b ) = 0 ∀a = b.
The goal of the low rank decomposition is to find the estimatorsx αi ∈ R r that minimize the mean square errorx = arg min To solve the Bayesian inference problem, using message passing, we frame it as a bipartite graphical model. Each of the variable nodes corresponds to an estimator x αi [See figure 1.b in the main text]. We use the notation ∂a to denote all the neighboring nodes to a and the notation ∂a\iα to denote all the neighboring variable nodes adjacent to a, excluding the node αi. The cardinality of the set of all factor points is |{a}| = N P .
Each variable point on the graph is connected to N p−1 factor nodes. The set of neighboring factor nodes that is connected to the variable node αi is denoted as

A.2 Weakly connected graph
An underlying assumption in belief propagation and message passing algorithms is that the incoming messages into each node are statistically independent. It can be achieved, for example in sparse graphs, where at each node the graph can be approximately considered as a tree (directed acyclic graph), without recurring loops. In the physics literature such approximation is often referred to as Bethe Lattice. In the current model this is possible due to the scaling of individual elements, w ∼ N (1−p)/2 , as defined in Eq. (1) in the main text. Since correlations in the messages are due to loops in the underlying graph, that pass through several nodes, we neglect them when the interactions are sufficiently weak [38].

A.3 Message passing
We start by defining two different types of messages (beliefs): one for messages outgoing from a variable node into a factor node η. Messages going from factor nodes into variable nodes are denoted byη. Messages are the marginal probabilities at each node, measuring the posterior probability density of the estimator at that source node. A message outgoing from the variable node αi to a factor node a can be written in terms of the product of messages originating from all its connected nodes excluding a, The denominator Z αi→a is a normalization factor Incoming messages into variable nodes are obtained by marginalizing over distribution over all the messages. A message outgoing from a factor node into a variable node is given bỹ where The normalization factor in the denominator of eq. (30) is given by The cost function g(·) at the exponent can be expanded as a power series in Ñ where S b and R b are the first and second derivative of the cost function g(Y, w) evaluated at Y b and w b = 0:

A.3.1 Belief propagation
The mean values of outgoing messages from variable node αi into factor node a, are obtained by integrating over the marginal probabilities η αi→a : Note that we have used the transpose of the vector x T αi , which will become useful for the notation below. Their covariance matrix is equal tô Using the first-and second-order statistics, we can write explicit expressions for the moments of w b appearing in the expansion (32) above. The first moment reads Similarly, the second moment, w 2 b , is given by Introducing the explicit moments back into eq. (32), the incoming messages into variable nodes are given bỹ Since we are interested in the marginals in the variable nodes, we can replace this result in the expression for messages outgoing from a variable node (28), we obtain Note that g(Y b , 0) is a constant and can be absorbed into the normalization function. We define the two order parameters and Using the order parameters we rewrite equation (40) as The normalization, or partition function Z αi→a , can be written in terms of the order parameters u T αi→a and A βj→b as Finally, the moments of the local variables x αi with distribution η αi→a (x αi ) can be found directly from the partition functions Z αi→a by standard derivations. The mean is given bŷ and the covariance matrices are

A.4 AMP algorithms
The mean-field equations, describing the equilibrium of the local estimators can be used to iteratively into an algorithm by iteratively calculating the statistics of the messages given their estimators using eq. (41) and (42) and then reevaluating the estimatorsx and σ using eq. (45) and (46). Defining the upper-script t denoting the time step of the algorithm iteration, we can write the iterative equations as A.5 Approximate message passing -local mean-field approximation for the messages In the equations above (47)-(51), the number of overall messages (and thus calculations) scale with the number of edges in the factorized graph, i.e. as O(N P ). However, the dependence of each message on the state of the target node is weak. Therefore, the values of u t αi→a and A t αi→a are very close to their mean, when marginalized over all target nodes a. The local deviations about that mean scale as N (1−p)/2 . For that reason, we can consider the statistics of all outgoing messages from each node (i.e., average over all the adjacent edges), and assume small fluctuations due to the state of the targets. This procedure is essentially performing mean-field approximation at every node. The result will be the AMP equations which scale with the number of variable nodes P N, rather than with the number of edges in the graph. In physics, this analogous to the cavity method.
To apply this reasoning to the equations, we define the order parameters A αi and u αi , which explicitly exclude the dependence of the target node: The difference between the non-directed and the directed messages is the component that depends on the target node S a . For the mean-messages, the correction terms scale as O(N (1−p)/2 ), and is given by For the fluctuations in the local messages about their mean, the correction term scales as and we will be neglecting it.
To transform the equations for the local messages statistics, to use only the target-agnostic variables, x t βj and σ t βj , we calculate the difference between the two mean values We develop the second term on the RHS to linear order in the small parameter of the difference δu αi→a , and note that the leading order cancel with the first term in the RHS above, yielding Using eq. (54) and (57), we can write an expression for the node-average local messages, Expanding the product of the βj factors, and keeping terms up to linear order in the small difference δx, we obtain The first correction for the above, involves the quadratic terms in the expansion of (58). The mixed terms involve the values at time t and at time t − 1, which originate from the expansion of δx in (57). The mixed term is given by This expression, which couples the dynamical variablex into its previous time step is an Onsager response term. It reflects the changes to the fields of the nodes surrounding the node αi due to the activity of the node αi in the previous time step.
Importantly, up until this point, we have not yet used the assumption of Bayes optimality, nor have we used the Nishimori identities that follow the Bayes-optimal assumption. Consequently, algorithms based on the approximate message-passing above should be general and do not require the Bayes-optimal assumption. In the following section, we consider simplification due to the Bayesoptimal assumption. Beyond simplification of the mathematical expressions, it will allow us to systematically derive a dynamical mean-field theory for the errors in section B.

A.6 Simplifications for Bayes-optimal settings
The covariance matrix in the Bayes-optimal case can be much simplified. First, one can show, using the Nishimori identities at the equilibrium, that Here the angular brackets denote averaging over the posterior. The posterior variance of S b is given by the Fisher information of the output channel E post S 2 b = 1 ∆ . Furthermore, in Bayes-optimal setting all samples from the equilibrium ensemble are similar, the sum over nodes becomes selfaveraging and so Using the above simplifications, the covariance matrix of the messages can be written as Importantly, the covariance matrix does not depend on the specific node i, but it does depend on the mode of the tensor α. This is a significant difference from the algorithms for symmetric-tensor decomposition, for which A t was uniform for all nodes in the factor graph [21].
The Onsager term, which appears in the iterative mean-field equations for the local mean of the messages can also be simplified under the Bayes-optimal setting. Using some algebra, the Onsager correction becomes Finally, we write the simplified AMP equations as where the Onsager term is given by and the partition function reads There are two parameters in these equations (apart from the prior distributions P α (x). One is the Fisher information of the output channel, which is a global parameter that can tune the global dynamics. The other S b , which is the Fisher score of the entry at Y b . The last one is what yields the structure in the solution of the estimators.

B Dynamic mean field theory (state evolution)
In the previous section we have derived the AMP algorithm for general tensors, and show their simplified form in the case of the Bayes-optimal assumption, where the priors are known, and the system follow Nishimori identities at equilibrium. These algorithms follow the iterative evolution of the estimators in each of the variable nodes in the factor graph. In order to analytically study the performance of the algorithm, we want to know how the mean error reduces from one iteration of the algorithm to the next. To do that, we derive a dynamical mean-field theory (also known as state-evolution). As mentioned above, following the Bayes-optimal assumption, the estimators are self-averaging; thus a mean-field description of the error is a good measure for the typical evolution of any given system.
We define an order parameter that measures the overlap between each of the underlying vectors of estimatorsx t α ∈ R p×r and the ground truth values x 0 α ∈ R p×r . The overlap matrix is defined as In total, there are p matrices of dimensions r × r, each for each mode of the tensor. In the Bayesoptimal regime, the ground-truth values can be replaced with any typical sample from the posterior distribution. Thus, in Bayes-optimal inference, M t α is also the typical covariance matrix of the estimators It follows that under the Bayes-optimality condition, M t α is a symmetric matrix. To study the typical dynamics of the algorithm using the mean overlap, we derive yet another meanfield approximation, now on the spatial degrees of freedom -i.e., the nodes. Given the self-averaging property of the nodes under the Bayes-optimal setting, and using the central-limit theorem, we need to find the first two moments of the distribution of the local values u iα (note that A α are already node-independent). Following the usual procedure of mean-field theory, we then close the equations self-consistently using the overlap parameter M t α . Using the definition of u αi from eq. (52), we average over the posterior P out : Similar to the the approximation carried above for the AMP algorithms, we develop the posterior probability about w = 0, and keep only the leading terms, Carrying the integration in eq. (75), the leading order will vanisĥ which is the consequence of the Nishimori identities. In other words, in a Bayes-optimal setting, and when the interactions are weak, then the average value of the messages when averaged over the entire graph vanish to leading order. Intuitively, since the underlying graph is isotropic, we expect that the dynamics will be similar at every node on average.
Performing the integration on the next, quadratic, term in (75) we get Note that the original tensor components, denoted by w a are the ground-truth in the context of the inference problem, and we can write Replacing this into the expression for the expectation above we get Finally we can write In a similar manner, we can calculate the covariance matrix of the mean-messages Keeping the leading order after the expansion of the distribution P out for small w we get In the Bayes-optimal setting this is equal to While the mean-message u αi varies from node to node, the mean covariance (not to be confused with the covariance of the mean calculated above), A t α , is node-independent, as we have established in the previous section. In the Bayes-optimal setting, where the Nishimori identities hold, it is equal to Using the definition of the mean overlap in eq. (73), and eq. (68), we write where and with the partition function Replacing the average over all nodes i in (86) with the expectation, we write an iterative update equation for the order parameter M t α , Here, z ∈ R r are random variables with standard normal distribution. The expectation in the RHS of eq. (89) is over two random variables: First are expected values for the underlying ground-truth x 0 α , which follows the prior distribution P α ; The second is of a standard gaussian variable z, which represent the node-to-node fluctuations in the local mean-messages, with mean 1 Mean square error The real quantity of interest is the mean square error (MSE) of the estimate. This can be easily obtained from the mean overlap at any time-step of the algorithm using.
Finally, the expected error of the AMP algorithms, once it has converged is given by

C Convergence of the AMP algorithms
Approximate message passing, and belief-propagation algorithms in general are known to have convergence issues (see for example [39,40,41,42,43]). A typical naive implementation of the algorithms will reduce the overall mean square error of the estimator, M SE t . However, at some point, M SE t will start increasing and may diverge to large deviations from the ground-truth values or will oscillate about some fixed value. Loosely speaking, the step size of the iterative update equations (68) and (69) is too big, and the algorithm may 'overshoot' the MMSE estimator. One possible way to correct this behavior (see e.g., [42] and reference therein) is to reduce the step size. Since the differential change to x t and σ t is proportional to derivatives of the partition function in (68) and (69), a good normalization scheme could use an energy estimation of the configuration at time-step t. To do this, one can evaluate the Bethe free energy at every time step [42,43,20]. However, since this report does not focus on possible implementations of the algorithms, it is sufficient to use a simpler -and potentially less efficient -scheme, using fixed step-size reduction, or damping.
To implement the fixed damping algorithm, eq (68) and (69) can be rewritten aŝ where 0 ≤ λ < 1 is the damping coefficient that controls the effective step size and the speed of convergence. In this simple implementation, the level of damping is a control parameter of the algorithm. A more sophisticated approach would use adaptive dampingλ t , where the effective step size decreases as the Bethe free energy of the configuration {x t α } decreases [39,20].

D Noncubic tensors
In the above, we have assumed that the dimensionality of all p modes is N , implying that the underlying tensor is cubic (i.e., all modes have the same length). To study how the shape of the tensors influence the AMP algorithm and the performance, we allow for the different modes to have different dimensionality N α . Importantly, we assume that all modes are in the thermodynamic regime, i.e., N α → ∞ α = {1, ..., p}. Furthermore, we assume all modes scale in a similar way. This is done by defining N α = n α N where all n α = O(1) and α n α = 1. The thermodynamic limit is then understood by taking N → ∞.
First we note that the scaling of the tensor elements does not change with this choice of scaling, However, the algorithms have no symmetry with respect to the dimensionality of the different modes in this case. This broken symmetry is in the iterative mean-field equations for the local mean messages u αi , in eq. (66) , which now is scaled by proportion of the dimensionality respective mode: The other mean-field equations of the algorithms are left unchanged.
Correction to the dynamic mean-filed equations In order to make the necessary changes to the dynamic mean-field theory in section B, we redefine the mean overlap with the appropriate scaling, which now depends on the mode α, Using the rescaled overlap, we re-derive the iterative dynamic mean-field equations, following the same steps as in section B: Substituting the expression for w, Finally we arrive at Note that the only difference between this result and eq. (81) is the factor n α . It follows that in order to generalize the dynamic mean-field theory to noncubic tensors, we simply need to replace M α → n α M α throughout the results of section B. The final equations are presented in the main text.
In order to theoretically evaluate the performance of the AMP algorithms given different tensors, we explicitly derive the dynamic mean-field equation of the overlap, and the error, for some common priors P α (x). These derivations closely resemble the analysis done in [20] for the p = 2 case, only here we allow different mixture of prior and different sizes for the modes, and consider arbitrary order p. In the following we will use rank r = 1 tensors, where the estimators x, the ground truth x 0 , and the the overlaps M , which we will denote here as m, are all scalar values. The same analysis hods with multivariate calculation, when r ≥ 2.
In order to theoretically evaluate the performance of the AMP algorithms given different tensors, we explicitly derive the dynamic mean-field equations of the overlap, and the error, for some common priors P α (x). These derivations closely resemble the analysis done in [20] for the p = 2 case, only here we allow a different mixture of prior and different sizes for the modes, and consider arbitrary order p. In the following, we will use rank r = 1 tensors, where the estimators x, the ground truth x 0 , and the overlaps M -which we will denote here as m -are all scalar values. The same analysis holds with multivariate calculation when r ≥ 2.

E.1 Gaussian prior
The first, and perhaps most common, choice for prior is a normal distribution of x α , with variance σ 2 α and and mean µ α , We use that prior to explicitly calculate the update rule, The nominator of (101) can be written aŝ where b = uσ 2 α and a = Aσ 2 α .Completing the quadratic form, we have Similar treatment in performed on the denominator. It is straight forward to see that the function in (101) reduces to Next, we want to use this functional form in the dynamic mean-field eq. (89). Denotem t α = 1 ∆ ⊙ β =α m t β , then we have A t =m t α and u t α =m t α x 0 α + m t α z then we want to compute Averaging over the distribution of the ground-truth values P (x 0 ), Finally the dynamic mean-field iterative equation on the mean overlap are given by In the case of zero mean priors, µ α = 0, the equation is reduced to We note that if all modes are Gaussian with zero means, then the solution M α = 0 ∀α is a stable fixed point of the dynamics, implying that if we start from random initial conditions, that are uncorrelated with the true values, the algorithms will not converge. A numerical analysis of eq. (106) for order p = 3 tensors is presented in the main text.
Consider the case of µ α = µ and σ α = σ, with all priors are similar. From the structure of (106) we find that in the fixed point M * which is what would be expected from the symmetry of the problem. Note however that unlike the derivation in [21], the underlying tensor in non-symmetric.
Noncubic tensors. If we have different population sizes, then we have a ratio between the order parameters, and replace m α with n α m α .

E.2 Bernoulli distribution
For many applications, it is expected that some of the modes in the underlying low rank tensors are sparse, meaning they contribute information to only a small subset of the measurements. A simple way of modeling such data is using the Bernoulli distribution, P α (x) = ρδ(x − 1) + (1 − ρ)δ(x).
As in the derivation of the Gaussian priors in the previous subsection, we compute the function (101). The nominator is equal tô with first derivative equal to In the bayes optimal case we have A α =m α and u αi =m α x 0 i + √m α z ,where we have defined m α ≡ 1 ∆ β =α m β . In the expression in the exponent of (114) we have Next we integrate over the prior and ground-truth to get m t+1 α = ρE z f α m α ,m α + m α z = In the sparse case, where ρ ≪ 1this can be simplified further Note that in a complete overlap we have m = ρ so em α/2 = 1/ρ and In instances where all of the modes have similar statistics, then we would have 1 2∆ ρ p−1 = − log ρ ⇒ ∆ = ρ p−1 2 |log ρ| .
Here, we can expect that for ∆ ∼ ρ p−1 / |log ρ|, where p is the order of the tensor, we will have high overlap with zero error. However, in the case of non-symmetric tensors, not all directions have to be sparse, and may have different distributions. In that case the noise scale as ∆ ∼ ρp −1 / |log ρ|, wherep is the number of sparse modes in the underlying tensor.

E.3 Gauss-Bernoulli
The next logical step is to combine the continuous irregularity of the Gaussian distribution and the sparse nature of the Bernoulli distribution. The Gauss-Bernoulli distribution is given by For brevity we will use zero mean µ = 0 and unit variance σ 2 = 1, and note that the results can be easily rescaled. The update function is given by .
Using some algebra we get

E.4 Mixed priors
In the case of general asymmetric tensors, we can construct a tensor using different priors for the different modes. It is particularly useful in real applications, as different modes of the tensors can originate from entirely different sources. Consider for example an order-3 tensor holding neural firing rate data r itk . The index i marks the neuron recorded; index t is the time bin within a single trial, and k is the trial index. If we believe that the data originates from the low-dimensional dynamical system, we would want to write the tensor as where D is the dimensions of the dynamical system, and ∆ is the noise of a single measurement.We may ask how should we design an experiment so that low-rank decomposition of the recorded data would be possible. In this case, we would assert different priors to the different modes. The mode x t represent the D dimensional dynamical system. We could assume for example that is generated by some Gaussian process, thus follows Gaussian statistics. The mode u i represents the projections of the low-dimensional dynamical system onto the set measured neurons. It may be a valid assumption that only a fraction of the neurons responds in coherence with the underlying dynamics; a Gauss-Bernoulli distribution will be suitable for this mode. Lastly, the trial modulus mode v k can have Gaussian distribution about some mean with small variance, suggesting small trial-to-trial modulations.
To solve the dynamic mean field theory for this case, and find the boundaries of the inference we would use the appropriate equation for each of the modes. For example, for two Gaussian distributions and one Gauss-Bernoulli, we would have