Bilinear Generalized Approximate Message Passing

I Introduction

and we likewise assume that the likelihood function of $\boldsymbol{Z}$ is known and separable, i.e.,

Recently, various special cases of this problem have gained the intense interest of the research community, e.g.,

Matrix Completion: In this problem, one observes a few (possibly noise-corrupted) entries of a low-rank matrix and the goal is to infer the missing entries. In our framework, $\boldsymbol{Z}\!=\!\boldsymbol{AX}$ would represent the complete low-rank matrix (with tall $\boldsymbol{A}$ and wide $\boldsymbol{X}$) and $p_{\textsf{y}_{ml}|\textsf{z}_{ml}}$ the observation mechanism, which would be (partially) informative about $\textsf{z}_{ml}$ at the observed entries $(m,l)\in\Omega$ and non-informative at the missing entries $(m,l)\notin\Omega$.

Robust PCA: Here, the objective is to recover a low-rank matrix (or its principal components) observed in the presence of noise and sparse outliers. In our framework, $\boldsymbol{Z}\!=\!\boldsymbol{AX}$ could again represent the low-rank matrix, and $p_{\textsf{y}_{ml}|\textsf{z}_{ml}}$ the noise-and-outlier-corrupted observation mechanism. Alternatively, $\boldsymbol{X}$ could also capture the outliers, as described in the sequel.

Dictionary Learning: Here, the objective is to learn a dictionary $\boldsymbol{A}$ for which there exists a sparse data representation $\boldsymbol{X}$ such that $\boldsymbol{AX}$ closely matches the observed data $\boldsymbol{Y}$. In our framework, $\{p_{\textsf{x}_{nl}}\}$ would be chosen to induce sparsity, $\boldsymbol{Z}=\boldsymbol{AX}$ would represent the noiseless observations, and $\{p_{\textsf{y}_{ml}|\textsf{z}_{ml}}\}$ would model the (possibly noisy) observation mechanism.

In the context of CS, the AMP framework yields algorithms with remarkable properties: i) solution trajectories that, in the large-system limit (i.e., as $M,N\rightarrow\infty$ with $M/N$ fixed, under iid sub-Gaussian $\boldsymbol{A}$) are governed by a state-evolution whose fixed points—when unique—yield the true posterior means and ii) a low implementation complexity (i.e., dominated by one multiplication with $\boldsymbol{A}$ and $\boldsymbol{A}^{\textsf{T}}$ per iteration, and relatively few iterations) . Thus, a natural question is whether the AMP framework can be successfully applied to the generalized bilinear problem described earlier.

In this manuscript, we propose an AMP-based approach to generalized bilinear inference that we henceforth refer to as Bilinear Generalized AMP (BiG-AMP), and we uncover special cases under which the general approach can be simplified. In addition, we propose an adaptive damping mechanism, an expectation-maximization (EM)-based method of tuning the parameters of $p_{\textsf{a}_{mn}}$, $p_{\textsf{x}_{nl}}$, and $p_{\textsf{y}_{ml}|\textsf{z}_{ml}}$ (in case they are unknown), and methods to select the rank $N$ (in case it is unknown). In the case that $p_{\textsf{a}_{mn}}$, $p_{\textsf{x}_{nl}}$, and/or $p_{\textsf{y}_{ml}|\textsf{z}_{ml}}$ are completely unknown, they can be modeled as Gaussian-mixtures with mean/variance/weight parameters learned via EM . Finally, we present a detailed numerical investigation of BiG-AMP applied to matrix completion, robust PCA, and dictionary learning. Our empirical results show that BiG-AMP yields an excellent combination of estimation accuracy and runtime when compared to existing state-of-the-art algorithms for each application.

Although the AMP methodology is itself restricted to separable known pdfs (1)-(3), the results of Part II suggest that this limitation is not an issue for many practical problems of interest. However, in problems where the separability assumption is too constraining, it can be relaxed through the use of hidden (coupling) variables, as originally proposed in the context of “turbo-AMP” and applied to BiG-AMP in . Due to space limitations, however, this approach will not be discussed here. Finally, although we focus on real-valued random variables, all of the methodology described in this work can be easily extended to circularly symmetric complex-valued random variables.

We now discuss related work. One possibility of applying AMP methods to matrix completion was suggested by Montanari in [31, Sec. 9.7.3] but the approach described there differs from BiG-AMP in that it was i) constructed from a factor graph with vector-valued variables and ii) restricted to the (incomplete) additive white Gaussian noise (AWGN) observation model. Moreover, no concrete algorithm nor performance evaluation was reported. Since we first reported on BiG-AMP in , Rangan and Fletcher proposed an AMP-based approach for the estimation of rank-one matrices from AWGN-corrupted observations, and showed that it can be characterized by a state evolution in the large-system limit. More recently, Krzakala, Mézard, and Zdeborová proposed an AMP-based approach to blind calibration and dictionary learning in AWGN that bears similarity to a special case of BiG-AMP, and derived a state-evolution using the cavity method. Their method, however, was not numerically successful in solving dictionary learning problems . The BiG-AMP algorithm that we derive here is a generalization of those in in that it handles generalized bilinear observations rather than AWGN-corrupted ones. Moreover, our work is the first to detail adaptive damping, parameter tuning, and rank-selection mechanisms for AMP based bilinear inference, and it is the first to present an in-depth numerical investigation involving both synthetic and real-world datasets. An application/extension of the BiG-AMP algorithm described here to hyperspectral unmixing (an instance of non-negative matrix factorization) was recently proposed in .

The remainder of the document is organized as follows. Section II derives the BiG-AMP algorithm, and Sec. III presents several special-case simplifications of BiG-AMP. Section IV describes the adaptive damping mechanism, and Sec. V the EM-based tuning of prior parameters and selection of rank $N$. Application-specific issues and numerical results demonstrating the efficacy of our approach for matrix completion, robust PCA, and dictionary learning, are discussed in Sections VI–VIII, respectively, and concluding remarks are offered in Sec. IX.

II Bilinear Generalized AMP

For the statistical model (1)-(3), the posterior distribution is

where (4) employs Bayes’ rule and $\propto$ denotes equality up to a constant scale factor.

The posterior distribution can be represented with a factor graph, as depicted in Fig. 1. There, the factors of $p_{\textsf{{{X}}},\textsf{{{A}}}|\textsf{{{Y}}}}$ from (6) are represented by “factor nodes” that appear as black boxes, and the random variables are represented by “variable nodes” that appear as white circles. Each variable node is connected to every factor node in which that variable appears. The observed data $\{y_{ml}\}$ are treated as parameters of the $p_{\textsf{y}_{ml}|\textsf{z}_{ml}}$ factor nodes in the middle of the graph, and not as random variables. The structure of Fig. 1 becomes intuitive when recalling that $\textsf{{{Z}}}=\textsf{{{A}}}\textsf{{{X}}}$ implies $\textsf{z}_{ml}=\sum_{n=1}^{N}\textsf{a}_{mn}\textsf{x}_{nl}$.

II-B Loopy Belief Propagation

In this work, we aim to compute minimum mean-squared error (MMSE) estimates of $\boldsymbol{X}$ and $\boldsymbol{A}$, i.e., the meansAnother worthwhile objective could be to compute the joint MAP estimate $\arg\max_{\boldsymbol{X},\boldsymbol{A}}p_{\textsf{{{X}}},\textsf{{{A}}}|\textsf{{{Y}}}}(\boldsymbol{X},\boldsymbol{A}\,|\,\boldsymbol{Y})$; we leave this to future work. of the marginal posteriors $p_{\textsf{x}_{nl}|\textsf{{{Y}}}}(\cdot\,|\,\boldsymbol{Y})$ and $p_{\textsf{a}_{mn}|\textsf{{{Y}}}}(\cdot\,|\,\boldsymbol{Y})$, for all pairs $(n,l)$ and $(m,n)$. Although exact computation of these quantities is generally prohibitive, they can be efficiently approximated using loopy belief propagation (LBP) .

In LBP, beliefs about the random variables (in the form of pdfs or log pdfs) are propagated among the nodes of the factor graph until they converge. The standard way to compute these beliefs, known as the sum-product algorithm (SPA) , stipulates that the belief emitted by a variable node along a given edge of the graph is computed as the product of the incoming beliefs from all other edges, whereas the belief emitted by a factor node along a given edge is computed as the integral of the product of the factor associated with that node and the incoming beliefs on all other edges. The product of all beliefs impinging on a given variable node yields the posterior pdf for that variable. In cases where the factor graph has no loops, exact marginal posteriors result from two (i.e., forward and backward) passes of the SPA . For loopy factor graphs, exact inference is in general NP hard and so LBP does not guarantee correct posteriors. That said, LBP has shown state-of-the-art performance in many applications, such as inference on Markov random fields , turbo decoding , LDPC decoding , multiuser detection , and compressive sensing .

In high-dimensional inference problems, exact implementation of the SPA is impractical, motivating approximations of the SPA. A notable example is the generalized approximate message passing (GAMP) algorithm, developed in to solve the generalized CS problem, which exploits the “blessings of dimensionality” that arise when $\boldsymbol{A}$ is a sufficiently large and dense and which was rigorously analyzed in . In the sequel, we derive an algorithm for the generalized bilinear inference BiG-AMP algorithm that employs GAMP-like approximations to the SPA on the factor graph in Fig. 1. As we shall see, the approximations are primarily based on central-limit-theorem (CLT) and Taylor-series arguments.

II-C Sum-product Algorithm

Applying the SPA to the factor graph in Fig. 1, we arrive at the following update rules for the four messages in Table I.

where const is an arbitrary constant (w.r.t $x_{nl}$ in (7) and (8), and w.r.t $a_{mn}$ in (9) and (10)). In the sequel, we denote the mean and variance of the pdf $\frac{1}{C}\exp(\Delta_{m{\scriptscriptstyle\leftarrow}nl}^{\textsf{x}}(t,.))$ by $\widehat{x}_{m,nl}(t)$ and $\nu^{x}_{m,nl}(t)$, respectively, and we denote the mean and variance of $\frac{1}{C}\exp(\Delta_{l{\scriptscriptstyle\leftarrow}mn}^{\textsf{a}}(t,.))$ by $\widehat{a}_{l,mn}(t)$ and $\nu^{a}_{l,mn}(t)$. For the log-posteriors, the SPA implies

and we denote the mean and variance of $\frac{1}{C}\exp(\Delta_{nl}^{\textsf{x}}(t,.))$ by $\widehat{x}_{nl}(t)$ and $\nu^{x}_{nl}(t)$, and the mean and variance of $\frac{1}{C}\exp(\Delta_{mn}^{\textsf{a}}(t,.))$ by $\widehat{a}_{mn}(t)$ and $\nu^{a}_{mn}(t)$.

II-D Approximated Factor-to-Variable Messages

We now apply AMP approximations to the SPA updates (7)-(12). As we shall see, the approximations are based primarily on central-limit-theorem (CLT) and Taylor-series arguments that become exact in the large-system limit, where $M,L,N\to\infty$ with fixed ratios $M/N$ and $L/N$. (Due to the use of finite $M,L,N$ in practice, we still regard them as approximations.) In particular, our derivation will neglect terms that vanish relative to others as $N\rightarrow\infty$, which requires that we establish certain scaling conventions. First, we assume w.l.o.gOther scalings on $\operatorname{E}\{\textsf{z}_{ml}^{2}\}$, $\operatorname{E}\{\textsf{x}_{nl}^{2}\}$, and $\operatorname{E}\{\textsf{a}_{mn}^{2}\}$ could be used as long as they are consistent with the relationship $\textsf{z}_{ml}=\sum_{n=1}^{N}\textsf{a}_{mn}\textsf{x}_{nl}$. that $\operatorname{E}\{\textsf{z}_{ml}^{2}\}$ and $\operatorname{E}\{\textsf{x}_{nl}^{2}\}$ scale as $O(1)$, i.e., that the magnitudes of these elements do not change as $N\rightarrow\infty$. In this case, the relationship $\textsf{z}_{ml}=\sum_{n=1}^{N}\textsf{a}_{mn}\textsf{x}_{nl}$ implies that $\operatorname{E}\{\textsf{a}_{mn}^{2}\}$ must scale as $O(1/N)$. These scalings are assumed to hold for random variables $\textsf{z}_{ml}$, $\textsf{a}_{mn}$, and $\textsf{x}_{ml}$ distributed according to the prior pdfs, according to the pdfs corresponding to the SPA messages (7)-(10), and according to the pdfs corresponding to the SPA posterior approximations (11)-(12). These assumptions lead straightforwardly to the scalings of $\widehat{z}_{ml}(t)$, $\nu^{z}_{ml}(t)$, $\widehat{x}_{m,nl}(t)$, $\nu^{x}_{m,nl}(t)$, $\widehat{x}_{nl}(t)$, $\nu^{x}_{nl}(t)$, $\widehat{a}_{l,mn}(t)$, $\nu^{a}_{l,mn}(t)$, $\widehat{a}_{mn}(t)$, and $\nu^{a}_{mn}(t)$ specified in Table II. Furthermore, because $\Delta_{m{\scriptscriptstyle\rightarrow}nl}^{\textsf{x}}(t,\cdot)$ and $\Delta_{nl}^{\textsf{x}}(t,\cdot)$ differ by only one term out of $M$, it is reasonable to assume that the corresponding difference in means $\widehat{x}_{m,nl}(t)-\widehat{x}_{nl}(t)$ and variances $\nu^{x}_{m,nl}(t)-\nu^{x}_{nl}(t)$ are both $O(1/\sqrt{N})$, which then implies that $\widehat{x}^{2}_{m,nl}(t)-\widehat{x}^{2}_{nl}(t)$ is also $O(1/\sqrt{N})$. Similarly, because $\Delta_{l{\scriptscriptstyle\rightarrow}mn}^{\textsf{a}}(t,\cdot)$ and $\Delta_{mn}^{\textsf{a}}(t,\cdot)$ differ by only one term out of $N$, where $\widehat{a}_{l,mn}(t)$ and $\widehat{a}_{mn}(t)$ are $O(1/\sqrt{N})$, it is reasonable to assume that $\widehat{a}_{l,mn}(t)-\widehat{a}_{mn}(t)$ is $O(1/N)$ and that both $\nu^{a}_{l,mn}(t)-\nu^{a}_{mn}(t)$ and $\widehat{a}^{2}_{l,mn}(t)-\widehat{a}^{2}_{mn}(t)$ are $O(1/N^{3/2})$. The remaining entries in Table II will be explained below.

We start by approximating the message $\Delta_{m{\scriptscriptstyle\rightarrow}nl}^{\textsf{x}}(t,.)$. Expanding (7), we find

For large $N$, the CLT motivates the treatment of $\textsf{z}_{ml}$, the random variable associated with the $z_{ml}$ identified in (13), conditioned on $\textsf{x}_{nl}=x_{nl}$, as Gaussian and thus completely characterized by a (conditional) mean and variance. Defining the zero-mean r.v.s $\widetilde{\textsf{a}}_{l,mn}\triangleq\textsf{a}_{mn}-\widehat{a}_{l,mn}(t)$ and $\widetilde{\textsf{x}}_{m,nl}=\textsf{x}_{nl}-\widehat{x}_{m,nl}(t)$, where $\textsf{a}_{mn}\sim\frac{1}{C}\exp(\Delta_{l{\scriptscriptstyle\leftarrow}mn}^{\textsf{a}}(t,\cdot))$ and $\textsf{x}_{nl}\sim\frac{1}{C}\exp(\Delta_{m{\scriptscriptstyle\leftarrow}nl}^{\textsf{x}}(t,\cdot))$, we can write

after which it is straightforward to see that

With this conditional-Gaussian approximation, (13) becomes

Unlike the original SPA message (7), the approximation (20) requires only a single integration. Still, additional simplifications are possible. First, notice that $\widehat{p}_{n,ml}(t)$ and $\nu^{p}_{n,ml}(t)$ differ from the corresponding $n$-invariant quantities

by one term. In the sequel, we will assume that $\widehat{p}_{ml}(t)$ and $\nu^{p}_{ml}(t)$ are $O(1)$ since these quantities can be recognized as the mean and variance, respectively, of an estimate of $\textsf{z}_{ml}$, which is $O(1)$. Writing the $H_{ml}$ term in (20) using (22)-(23),

where in (25) we used the facts that $\widehat{a}_{l,mn}(t)(\widehat{x}_{nl}(t)-\widehat{x}_{m,nl}(t))$ and $\nu^{a}_{l,mn}(t)(\widehat{x}_{m,nl}^{2}(t)-\widehat{x}_{nl}^{2}(t)))-\widehat{a}_{l,mn}^{2}(t)\nu^{x}_{m,nl}(t)-\nu^{a}_{l,mn}(t)\nu^{x}_{m,nl}(t)$ are both $O(1/N)$.

Rewriting (20) using a Taylor series expansion in $x_{nl}$ about the point $\widehat{x}_{nl}(t)$, we get

where $H_{ml}^{\prime}$ and $H_{mn}^{\prime\prime}$ are the first two derivatives of $H_{mn}$ w.r.t its first argument and $\dot{H}_{ml}$ is the first derivative w.r.t its second argument. Note that, in (26) and elsewhere, the higher-order terms in the Taylor’s expansion are written solely in terms of their scaling dependence on $N$, which is what will eventually allow us to neglect these terms (in the large-system limit).

We now approximate (26) by dropping terms that vanish, relative to the second-to-last term in (26), as $N\rightarrow\infty$. Since this second-to-last term is $O(1/N)$ due to the scalings of $\widehat{a}^{2}_{l,mn}(t)$, $\widehat{p}_{ml}(t)$, and $\nu^{p}_{ml}(t)$, we drop terms that are of order $O(1/N^{3/2})$, such as the final term. We also replace $\nu^{a}_{l,mn}(t)$ with $\nu^{a}_{mn}(t)$, and $\widehat{a}_{l,mn}^{2}(t)$ with $\widehat{a}_{mn}^{2}(t)$, since in both cases the difference is $O(1/N^{3/2})$. Finally, we drop the $O(1/N)$ terms inside the $H_{ml}$ derivatives, which can be justified by taking a Taylor series expansion of these derivatives with respect to the $O(1/N)$ perturbations and verifying that the higher-order terms in this latter expansion are $O(1/N^{3/2})$. All of these approximations are analogous to those made in previous AMP derivations, e.g., , , and .

Applying these approximations to (26) and absorbing $x_{nl}$-invariant terms into the const term, we obtain

Note that (27) is essentially a Gaussian approximation to the pdf $\frac{1}{C}\exp(\Delta_{m{\scriptscriptstyle\rightarrow}nl}^{\textsf{x}}(t,.))$.

computed according to the (conditional) pdf

where here C=\int_{z}p_{\textsf{y}_{ml}|\textsf{z}_{ml}}(y_{ml}\,|\,z)\mathcal{N}\big{(}z;\widehat{p}_{ml}(t),\nu^{p}_{ml}(t)\big{)}. In fact, (35) is BiG-AMP’s iteration-$t$ approximation to the true marginal posterior $p_{\textsf{z}_{ml}|\textsf{{{Y}}}}(\cdot|\boldsymbol{Y})$. We note that (35) can also be interpreted as the (exact) posterior pdf for $\textsf{z}_{ml}$ given the likelihood $p_{\textsf{y}_{ml}|\textsf{z}_{ml}}(y_{ml}|\cdot)$ from (3) and the prior \textsf{z}_{ml}\sim\mathcal{N}\big{(}\widehat{p}_{ml}(t),\nu^{p}_{ml}(t)\big{)} that is implicitly assumed by iteration-$t$ BiG-AMP.

Since $\textsf{{{Z}}}^{\textsf{T}}=\textsf{{{X}}}^{\textsf{T}}\textsf{{{A}}}^{\textsf{T}}$, the derivation of the BiG-AMP approximation of $\Delta_{l{\scriptscriptstyle\rightarrow}mn}^{\textsf{a}}(t,.)$ closely follows the derivation for $\Delta_{m{\scriptscriptstyle\rightarrow}nl}^{\textsf{x}}(t,.)$. In particular, it starts with (similar to (13))

where again the CLT motivates the treatment of $\textsf{z}_{ml}$, conditioned on $\textsf{a}_{mn}=a_{mn}$, as Gaussian. Eventually we arrive at the Taylor-series approximation (similar to (27))

II-E Approximated Variable-to-Factor Messages

We now turn to approximating the messages flowing from the variable nodes to the factor nodes. Starting with (8) and plugging in (27) we obtain

Since $\widehat{a}^{2}_{mn}(t)$ and $\nu^{a}_{mn}(t)$ are $O(1/N)$, and recalling $\widehat{s}^{2}_{ml}(t)$ and $\nu^{s}_{ml}(t)$ are $O(1)$, we take $\nu^{r}_{m,nl}(t)$ to be $O(1)$. Meanwhile, since $\widehat{r}_{m,nl}(t)$ is an estimate of $\textsf{x}_{nl}$, we reason that it is $O(1)$.

The mean and variance of the pdf associated with the $\Delta_{m{\scriptscriptstyle\leftarrow}nl}^{\textsf{x}}(t\!+\!1,.)$ approximation in (40) are

where here C=\int_{x}p_{\textsf{x}_{nl}}(x)\mathcal{N}\big{(}x;\widehat{r}_{m,nl}(t),\nu^{r}_{m,nl}(t)\big{)} and $g^{\prime}_{\textsf{x}_{n}l}$ denotes the derivative of $g_{\textsf{x}_{nl}}$ with respect to the first argument. The fact that (43) and (44) are related through a derivative was shown in .

We now derive approximations of $\widehat{x}_{m,nl}(t)$ and $\nu^{x}_{m,nl}(t)$ that avoid the dependence on the destination node $m$. For this, we introduce $m$-invariant versions of $\widehat{r}_{m,nl}(t)$ and $\nu^{r}_{m,nl}(t)$:

Comparing (45)-(46) with (41)-(42) and applying previously established scalings from Table II reveals that $\nu^{r}_{m,nl}(t)-\nu^{r}_{nl}(t)$ is $O(1/N)$ and that $\widehat{r}_{m,nl}(t)=\widehat{r}_{nl}(t)-\nu^{r}_{nl}(t)\widehat{a}_{mn}(t)\widehat{s}_{ml}(t)+O(1/N)$, so that (43) implies

Above, (48) follows from taking Taylor series expansions around each of the $O(1/N)$ perturbations in (47); (II-E) follows from a Taylor series expansion in the first argument of (48) about the point $\widehat{r}_{nl}(t)$; and (50) follows by neglecting the $O(1/N)$ term (which vanishes relative to the others in the large-system limit) and applying the definitions

which match (43)-(44) sans the $m$ dependence. Note that (50) confirms that the difference $\widehat{x}_{m,nl}(t)-\widehat{x}_{nl}(t)$ is $O(1/\sqrt{N})$, as was assumed at the start of the BiG-AMP derivation. Likewise, taking Taylor series expansions of $g^{\prime}_{\textsf{x}_{nl}}$ in (44) about the point $\widehat{r}_{nl}(t)$ in the first argument and about the point $\nu^{r}_{nl}(t)$ in the second argument and then comparing the result with (52) confirms that $\nu^{x}_{m,nl}(t)-\nu^{x}_{nl}(t)$ is $O(1/\sqrt{N})$.

We then repeat the above procedure to derive an approximation to $\Delta_{l{\scriptscriptstyle\leftarrow}mn}^{\textsf{a}}(t\!+\!1,.)$ analogous to (40), whose corresponding mean is then further approximated as

Arguments analogous to the discussion following (42) justify the remaining scalings in Table II.

II-F Closing the Loop

The penultimate step in the derivation of BiG-AMP is to approximate earlier steps that use $\widehat{a}_{l,mn}(t)$ and $\widehat{x}_{m,nl}(t)$ in place of $\widehat{a}_{mn}(t)$ and $\widehat{x}_{nl}(t)$. For this, we start by plugging (50) and (53) into (22), which yieldsRecall that the error of the approximation in (50) is $O(1/N)$ and the error in (53) is $O(1/N^{3/2})$.

where, for (60), we used $\widehat{a}_{mn}^{2}(t)$ in place of $\widehat{a}_{mn}(t)\widehat{a}_{mn}(t\!-\!1)$, used $\widehat{x}_{nl}^{2}(t)$ in place of $\widehat{x}_{nl}(t)\widehat{x}_{nl}(t\!-\!1)$, and neglected terms that are $O(1/\sqrt{N})$, since they vanish relative to the remaining $O(1)$ terms in the large-system limit.

Next we plug (50), (53), $\nu^{x}_{m,nl}(t)=\nu^{x}_{nl}(t)+O(1/\sqrt{N})$, and $\nu^{a}_{l,mn}(t)=\nu^{a}_{mn}(t)+O(1/N^{3/2})$ into (23), giving

where (62) retains only the $O(1)$ terms from (61).

Similarly, we plug (53) into (46) and (50) into (58) to obtain

where the approximations involve the use of $\widehat{s}^{2}_{ml}(t)$ in place of $\widehat{s}_{ml}(t)\widehat{s}_{ml}(t-1)$, of $\widehat{a}_{mn}(t)$ in place of $\widehat{a}_{mn}(t-1)$, of $\widehat{x}_{nl}(t)$ in place of $\widehat{x}_{nl}(t-1)$, and the dropping of terms that vanish in the large-system limit. Finally, we make the approximations

by neglecting the $\widehat{s}^{2}_{ml}(t)-\nu^{s}_{ml}(t)$ terms in (45) and (57), as explained in Appendix B.

II-G Approximated Posteriors

The final step in the BiG-AMP derivation is to approximate the SPA posterior log-pdfs in (11) and (12). Plugging (27) and (37) into those expressions, we get

using steps similar to (40). The associated pdfs are

for C_{x}\triangleq\int_{x}p_{\textsf{x}_{nl}}(x)\mathcal{N}\big{(}x;\widehat{r}_{nl}(t),\nu^{r}_{nl}(t)\big{)} and C_{a}\triangleq\int_{a}p_{\textsf{a}_{mn}}(a)\mathcal{N}\big{(}a;\widehat{q}_{mn}(t),\nu^{q}_{mn}(t)\big{)}, which are iteration-$t$ BiG-AMP’s approximations to the true marginal posteriors $p_{\textsf{x}_{nl}|\textsf{{{Y}}}}(x_{nl}\,|\,\boldsymbol{Y})$ and $p_{\textsf{a}_{mn}|\textsf{{{Y}}}}(a_{mn}\,|\,\boldsymbol{Y})$, respectively.

Note that $\widehat{x}_{nl}(t\!+\!1)$ and $\nu^{x}_{nl}(t\!+\!1)$ from (51)-(52) are the mean and variance, respectively, of the posterior pdf in (69). Note also that (69) can be interpreted as the (exact) posterior pdf of $\textsf{x}_{nl}$ given the observation $\textsf{r}_{nl}=\widehat{r}_{nl}(t)$ under the prior model $\textsf{x}_{nl}\sim p_{\textsf{x}_{nl}}$ and the likelihood model p_{\textsf{r}_{nl}|\textsf{x}_{nl}}(\widehat{r}_{nl}(t)\,|\,x_{nl};\nu^{r}_{nl}(t))=\mathcal{N}\big{(}\widehat{r}_{nl}(t);x_{nl},\nu^{r}_{nl}(t)\big{)} implicitly assumed by iteration-$t$ BiG-AMP. Analogous statements can be made about the posterior pdf of $\textsf{a}_{mn}$ in (70).

This completes the derivation of BiG-AMP.

II-H Algorithm Summary

The BiG-AMP algorithm derived in Sections II-C to II-G is summarized in Table III. There, we have included a maximum number of iterations, $T_{\textrm{max}}$ and a stopping condition (R17) based on the (normalized) change in the residual and a user-defined parameter $\tau_{\textrm{BiG-AMP}}$. We have also written the algorithm in a more general form that allows the use of complex-valued quantities [note the complex conjugates in (R10) and (R12)], in which case $\mathcal{N}$ in (D1)-(D3) would be circular complex Gaussian. For ease of interpretation, Table III does not include the important damping modifications that will be detailed in Sec. IV-A. Suggestions for the initializations in (I2) will be given in the sequel.

We note that BiG-AMP avoids the use of SVD or QR decompositions, lending itself to simple and potentially parallel implementations. Its complexity order is dominatedThe computations in steps (R4)-(R8) are $O(ML)$, while the remainder of the algorithm is $O(MN+NL)$. Thus, as $N$ grows, the matrix multiplies dominate the complexity. by ten matrix multiplications per iteration [in steps (R1)-(R3) and (R9)-(R12)], each requiring $MNL$ multiplications, although simplifications will be discussed in Sec. III.

The steps in Table III can be interpreted as follows. (R1)-(R2) compute a “plug-in” estimate $\boldsymbol{\overline{P}}$ of the matrix product $\textsf{{{Z}}}\!=\!\textsf{{{A}}}\textsf{{{X}}}$ and a corresponding set of element-wise variances $\{\overline{\nu}^{p}_{ml}\}$. (R3)-(R4) then apply “Onsager” correction (see and for discussions in the contexts of AMP and GAMP, respectively) to obtain the corresponding quantities $\boldsymbol{\widehat{P}}$ and $\{\nu^{p}_{ml}\}$. Using these quantities, (R5)-(R6) compute the (approximate) marginal posterior means $\boldsymbol{\widehat{Z}}$ and variances $\{\nu^{z}_{ml}\}$ of Z. Steps (R7)-(R8) then use these posterior moments to compute the scaled residual $\boldsymbol{\widehat{S}}$ and a set of inverse-residual-variances $\{\nu^{s}_{ml}\}$. This interpretation becomes clear in the case of AWGN observations with noise variance $\nu^{w}$, where

Steps (R9)-(R10) then use the residual terms $\boldsymbol{\hat{S}}$ and $\{\nu^{s}_{ml}\}$ to compute $\boldsymbol{\hat{R}}$ and $\{\nu^{r}_{nl}\}$, where $\widehat{r}_{nl}$ can be interpreted as a $\nu^{r}_{nl}$-variance-AWGN corrupted observation of the true $\textsf{x}_{nl}$. Similarly, (R11)-(R12) compute $\boldsymbol{\hat{Q}}$ and $\{\nu^{q}_{mn}\}$, where $\widehat{q}_{mn}$ can be interpreted as a $\nu^{q}_{mn}$-variance-AWGN corrupted observation of the true $\textsf{a}_{mn}$. Finally, (R13)-(R14) merge these AWGN-corrupted observations with the priors $\{p_{\textsf{x}_{nl}}\}$ to produce the posterior means $\boldsymbol{\hat{X}}$ and variances $\{\nu^{x}_{nl}\}$; (R15)-(R16) do the same for the $\textsf{a}_{mn}$ quantities.

The BiG-AMP algorithm in Table III is a direct (although non-trivial) extension of the GAMP algorithm for compressive sensing , which estimates X assuming perfectly known A, and even stronger similarities to the A-uncertain GAMP from , which estimates X assuming knowledge of the marginal means and variances of unknown random A, but which makes no attempt to estimate A itself. In Sec. III-B, a simplified version of BiG-AMP will be developed that is similar to the Bayesian-AMP algorithm for compressive sensing.

III BiG-AMP Simplifications

We now describe simplifications of the BiG-AMP algorithm from Table III that result from additional approximations and from the use of specific priors $p_{\textsf{y}_{ml}|\textsf{z}_{ml}}$, $p_{\textsf{x}_{nl}}$, and $p_{\textsf{a}_{mn}}$ that arise in practical applications of interest.

The BiG-AMP algorithm in Table III stores and processes a number of element-wise variance terms whose values vary across the elements (e.g., $\nu^{x}_{nl}$ can vary across $n$ and $l$). The use of scalar variances (i.e., uniform across $m,n,l$) significantly reduces the memory and complexity of the algorithm.

To derive scalar-variance BiG-AMP, we first assume $\forall n,l:\nu^{x}_{nl}(t)\approx\nu^{x}(t)\triangleq\frac{1}{NL}\sum_{n=1}^{N}\sum_{l=1}^{L}\nu^{x}_{nl}(t)$ and $\forall m,n:\nu^{a}_{mn}(t)\approx\nu^{a}(t)\triangleq\frac{1}{MN}\sum_{m=1}^{M}\sum_{n=1}^{N}\nu^{a}_{mn}(t)$, so from (R1)

Note that using (74) in place of (R1) avoids two matrix multiplies. Plugging these approximations into (R3) gives

which, when used in place of (R3), avoids another matrix multiply. Even with the above scalar-variance approximations, $\{\nu^{s}_{ml}(t)\}$ from (R5) are not guaranteed to be equal (except in special cases like AWGN $p_{\textsf{y}_{ml}|\textsf{z}_{ml}}$). Still, they can be approximated as such using $\nu^{s}(t)\triangleq\frac{1}{ML}\sum_{m=1}^{M}\sum_{l=1}^{L}\nu^{s}_{ml}(t)$, in which case

Using (76) in place of (R9) and (77) in place of (R11) avoids two matrix multiplies and $NL\!+\!MN\!-\!2$ scalar divisions, and furthermore allows (R10) and (R12) to be implemented as

saving two more matrix multiplies, and leaving a total of only three matrix multiplies per iteration.

III-B Possibly Incomplete AWGN Observations

We now consider a particular observation model wherein the elements of $\boldsymbol{Z}\!=\!\boldsymbol{AX}$ are AWGN-corrupted at a subset of indices $\Omega\subset(1\dots M)\!\times\!(1\dots L)$ and unobserved at the remaining indices, noting that the standard AWGN model (71) is the special case where $|\Omega|\!=\!ML$. This “possibly incomplete AWGN” (PIAWGN) model arises in a number of important applications, such as matrix completion and dictionary learning.

We can state the PIAWGN model probabilistically as

When the PIAWGN model is combined with the scalar-variance approximations from Sec. III-A, BiG-AMP simplifies considerably. To see this, we start by using $\nu^{p}(t)$ from (75) in place of $\nu^{p}_{ml}(t)$ in (81)-(82), resulting in

since $P_{\Omega}$ is a projection operator, and using (R4) and (83).

Scalar-variance BiG-AMP under PIAWGN observations is summarized in Table IV. Note that the residual matrix $\boldsymbol{\widehat{U}}(t)\triangleq P_{\Omega}(\boldsymbol{Y}-\boldsymbol{\widehat{A}}(t)\boldsymbol{\widehat{X}}(t))$ needs to be computed and stored only at the observed entries $(m,l)\in\Omega$, leading to significant savingsSimilar computational savings also occur with incomplete non-Gaussian observations. when the observations are highly incomplete (i.e., $|\Omega|\ll ML$). The same is true for the Onsager-corrected residual, $\boldsymbol{\widehat{V}}(t)$. Thus, the algorithm in Table IV involves only three (partial) matrix multiplies [in steps (R3p), (R8p), and (R10p), respectively], each of which can be computed using only $N|\Omega|$ scalar multiplies.

We note that Krzakala, Mézard, and Zdeborová recently proposed an AMP-based approach to blind calibration and dictionary learning that bears close similarityThe approach in does not compute (or use) $\nu^{p}(t)$ as given in lines (R4p)-(R5p) of Table IV, but rather uses an empirical average of the squared Onsager-corrected residual in place of our $\nu^{p}(t)+\nu^{w}$ throughout their algorithm. to BiG-AMP under the special case of AWGN-corrupted observations (i.e., $|\Omega|=ML$) and scalar variances. Their derivation differs significantly from that in Sec. II due to the many simplifications offered by this special case.

III-C Zero-mean iid Gaussian Priors on A and X

In this section we will investigate the simplifications that result in the case that both $p_{\textsf{a}_{mn}}$ and $p_{\textsf{x}_{nl}}$ are zero-mean iid Gaussian, i.e.,

which, as will be discussed later, is appropriate for matrix completion. In this case, straightforward calculations reveal that $\operatorname{E}\{\textsf{x}_{nl}\,|\,\textsf{r}_{nl}\!=\!\widehat{r}_{nl};\nu^{r}_{nl}\}=\widehat{r}_{nl}\nu^{x}_{0}/(\nu^{r}_{nl}+\nu^{x}_{0})$ and $\operatorname{var}\{\textsf{x}_{nl}\,|\,\textsf{r}_{nl}\!=\!\widehat{r}_{nl};\nu^{r}_{nl}\})=\nu^{x}_{0}\nu^{r}_{nl}/(\nu^{r}_{nl}+\nu^{x}_{0})$ and, similarly, that $\operatorname{E}\{\textsf{a}_{mn}\,|\,\textsf{q}_{mn}\!=\!\widehat{q}_{mn};\nu^{q}_{mn}\}=\widehat{q}_{mn}\nu^{a}_{0}/(\nu^{q}_{mn}+\nu^{a}_{0})$ and $\operatorname{var}\{\textsf{a}_{mn}\,|\,\textsf{q}_{mn}\!=\!\widehat{q}_{mn},\nu^{q}_{mn}\}=\nu^{a}_{0}\nu^{q}_{mn}/(\nu^{q}_{mn}+\nu^{a}_{0})$. Combining these iid Gaussian simplifications with the scalar-variance simplifications from Sec. III-A yields an algorithm whose computational cost is dominated by three matrix multiplies per iteration, each with a cost of $MNL$ scalar multiplies. The precise number of multiplies it consumes depends on the assumed likelihood model that determines steps (R7g)-(R8g).

Additionally incorporating the PIAWGN observations from Sec. III-B reduces the cost of the three matrix multiplies to only $N|\Omega|$ scalar multiplies each, and yields the “BiG-AMP-Lite” algorithm summarized in Table V, consuming $(3N+5)|\Omega|+3(MN+NL)+29$ multiplies per iteration.

IV Adaptive Damping

The approximations made in the BiG-AMP derivation presented in Sec. II were well-justified in the large system limit, i.e., the case where $M,N,L\rightarrow\infty$ with fixed $\frac{M}{N}$ and $\frac{L}{N}$. In practical applications, however, these dimensions (especially $N$) are finite, and hence the algorithm presented in Sec. II may diverge. In case of compressive sensing, the use of “damping” with GAMP yields provable convergence guarantees with arbitrary matrices . Here, we propose to incorporate damping into BiG-AMP. Moreover, we propose to adapt the damping of these variables to ensure that a particular cost criterion decreases monotonically (or near-monotonically), as described in the sequel. The specific damping strategy that we adopt is similar to that described in and coded in .

In BiG-AMP, the iteration-$t$ damping factor $\beta(t)\in(0,1]$ is used to slow the evolution of certain variables, namely $\overline{\nu}^{p}_{ml}$, $\nu^{p}_{ml}$, $\nu^{s}_{ml}$, $\widehat{s}_{ml}$, $\widehat{x}_{nl}$, and $\widehat{a}_{mn}$. To do this, steps (R1), (R3), (R7), and (R8) in Table III are replaced with

and the following are inserted between (R8) and (R9):

The newly defined state variables $\overline{x}_{nl}(t)$ and $\overline{a}_{mn}(t)$ are then used in place of $\widehat{x}_{nl}(t)$ and $\widehat{a}_{mn}(t)$ in steps (R9)-(R12) [but not (R1)-(R2)] of Table III. A similar approach can be used for the algorithm in Table IV (with the damping applied to $\boldsymbol{\widehat{V}}(t)$ instead of $\boldsymbol{\widehat{S}}(t)$) and those in Table V. Notice that, when $\beta(t)\!=\!1$, the damping has no effect, whereas when $\beta(t)\!=\!0$, all quantities become frozen in $t$.

IV-B Adaptive Damping

The idea behind adaptive damping is to monitor a chosen cost criterion $J(t)$ and decrease $\beta(t)$ when the cost has not decreased sufficientlyThe following adaptation procedure is borrowed from GAMPmatlab, where it has been established to work well in the context of GAMP-based compressive sensing. When the current cost $J(t)$ is not smaller than the largest cost in the most recent stepWindow iterations, then the “step” is deemed unsuccessful, the damping factor $\beta(t)$ is reduced by the factor stepDec, and the step is attempted again. These attempts continue until either the cost criterion decreases or the damping factor reaches stepMin, at which point the step is considered successful, or the iteration count exceeds $T_{\max}$ or the damping factor reaches stepTol, at which point the algorithm terminates. When a step is deemed successful, the damping factor is increased by the factor stepInc, up to the allowed maximum value stepMax. relative to $\{J(\tau)\}_{\tau=t-1-T}^{t-1}$ for some “step window” $T\geq 0$. This mechanism allows the cost criterion to increase over short intervals of $T$ iterations and in this sense is similar to the procedure used by SpaRSA . We now describe how the cost criterion $J(t)$ is constructed, building on ideas in .

Notice that, for fixed observations $\boldsymbol{Y}$, the joint posterior pdf solves the (trivial) KL-divergence minimization problem

The factorized form (5) of the posterior allows us to write

To judge whether a given time-$t$ BiG-AMP approximation “$b_{\textsf{{{X}}},\textsf{{{A}}}}(t)$” of the joint posterior $p_{\textsf{{{X}}},\textsf{{{A}}}|\textsf{{{Y}}}}$ is better than the previous approximation $b_{\textsf{{{X}}},\textsf{{{A}}}}(t\!-\!1)$, one could in principle plug the posterior approximation expressions (69)-(70) into (103) and then check whether $J(b_{\textsf{{{X}}},\textsf{{{A}}}}(t))<J(b_{\textsf{{{X}}},\textsf{{{A}}}}(t\!-\!1))$. But, since the expectation in (103) is difficult to evaluate, we approximate the cost (103) by using, in place of AX, an independent Gaussian matrixThe GAMP work uses a similar approximation. whose component means and variances are matched to those of AX. Taking the joint BiG-AMP posterior approximation $b_{\textsf{{{X}}},\textsf{{{A}}}}(t)$ to be the product of the marginals from (69)-(70), the resulting component means and variances are

In this way, the approximate iteration-$t$ cost becomes

Intuitively, the first term in (108) penalizes the deviation between the (BiG-AMP approximated) posterior and the assumed prior on X, the second penalizes the deviation between the (BiG-AMP approximated) posterior and the assumed prior on A, and the third term rewards highly likely estimates Z.

V Parameter Tuning and Rank Selection

Recall that BiG-AMP requires the specification of priors $p_{\textsf{{{X}}}}(\boldsymbol{X})=\prod_{n,l}p_{\textsf{x}_{nl}}(x_{nl})$, $p_{\textsf{{{A}}}}(\boldsymbol{A})=\prod_{m,n}p_{\textsf{a}_{mn}}(a_{mn})$, and $p_{\textsf{{{Y}}}|\textsf{{{Z}}}}(\boldsymbol{Y}|\boldsymbol{Z})=\prod_{m,l}p_{\textsf{y}_{ml}|\textsf{z}_{ml}}(y_{ml}|z_{ml})$. In practice, although one may know appropriate families for these distributions, the exact parameters that govern them are generally unknown. For example, one may have good reason to believe apriori that the observations are AWGN corrupted, justifying the choice $p_{\textsf{y}_{ml}|\textsf{z}_{ml}}(y_{ml}|z_{ml})=\mathcal{N}(y_{ml};z_{ml},\nu^{w})$, but the noise variance $\nu^{w}$ may be unknown. In this section, we outline a methodology that takes a given set of BiG-AMP parameterized priors $\{p_{\textsf{x}_{nl}}(\cdot;\boldsymbol{\theta}),p_{\textsf{a}_{mn}}(\cdot;\boldsymbol{\theta}),p_{\textsf{y}_{ml}|\textsf{z}_{ml}}(y_{ml}|\cdot;\boldsymbol{\theta})\}_{\forall m,n,l}$ and tunes the parameter vector $\boldsymbol{\theta}$ using an expectation-maximization (EM) based approach, with the goal of maximizing the likelihood, i.e., finding $\boldsymbol{\hat{\theta}}\triangleq\operatorname*{arg\,max}_{\boldsymbol{\theta}}p_{\textsf{{{Y}}}}(\boldsymbol{Y};\boldsymbol{\theta})$. The approach presented here can be considered as a generalization of the GAMP-based work to BiG-AMP.

Taking X, A, and Z to be the hidden variables, the EM recursion can be written as

As a concrete example, consider updating the noise variance $\nu^{w}$ under the PIAWGN model (80). Equation (109) suggests

where the true marginal posterior $p_{\textsf{z}_{ml}|\textsf{{{Y}}}}(\cdot|\boldsymbol{Y})$ is replaced with the most recent BiG-AMP approximation $p_{\textsf{z}_{ml}|\textsf{p}_{ml}}(\cdot|\widehat{p}_{ml}(T_{\max});\nu^{p}_{ml}(T_{\max}),\boldsymbol{\widehat{\theta}}^{k})$, where “most recent” is with respect to both EM and BiG-AMP iterations. Zeroing the derivative of the sum in (110) with respect to $\nu^{w}$,

where $\widehat{z}_{ml}(t)$ and $\nu^{z}_{ml}(t)$ are the BiG-AMP approximated posterior mean and variance from (33)-(34).

The overall procedure can be summarized as follows. From a suitable initialization $\boldsymbol{\hat{\theta}}^{0}$, BiG-AMP is run using the priors $\{p_{\textsf{x}_{nl}}(\cdot;\boldsymbol{\hat{\theta}}^{0}),p_{\textsf{a}_{mn}}(\cdot;\boldsymbol{\hat{\theta}}^{0}),p_{\textsf{y}_{ml}|\textsf{z}_{ml}}(y_{ml}|\cdot;\boldsymbol{\hat{\theta}}^{0})\}_{\forall m,n,l}$ and iterated to completion, yielding approximate marginal posteriors on $\{\textsf{x}_{nl},\textsf{a}_{mn},\textsf{z}_{ml}\}_{\forall m,n,l}$. These posteriors are used in (109) to update the parameters $\boldsymbol{\theta}$ one element at a time, yielding $\boldsymbol{\hat{\theta}}^{1}$. BiG-AMP is then run using the priors $\{p_{\textsf{x}_{nl}}(\cdot;\boldsymbol{\hat{\theta}}^{1}),p_{\textsf{a}_{mn}}(\cdot;\boldsymbol{\hat{\theta}}^{1}),p_{\textsf{y}_{ml}|\textsf{z}_{ml}}(y_{ml}|\cdot;\boldsymbol{\hat{\theta}}^{1})\}_{\forall m,n,l}$, and so on. A detailed discussion in the context of GAMP, along with explicit update equations for the parameters of Bernoulli-Gaussian-mixture pdfs, can be found in .

V-B Rank Selection

BiG-AMP and EM-BiG-AMP, as described up to this point, require the specification of the rank $N$, i.e., the number of columns in A (and rows in X) in the matrix factorization $\textsf{{{Z}}}=\textsf{{{A}}}\textsf{{{X}}}$. Since, in many applications, the best choice of $N$ is difficult to specify in advance, we now describe two procedures to estimate $N$ from the data $\boldsymbol{Y}$, building on well-known rank-selection procedures.

have been developed, such as the Bayesian Information Criterion (BIC) and Akaike’s Information Criterion (AIC) . In (112), $\boldsymbol{\hat{\Theta}}_{N}$ is the ML estimate of $\mathsfbf{\Theta}_{N}$ under $\boldsymbol{Y}$, and $\eta(\cdot)$ is a penalty function that depends on the effective number of scalar parameters $N_{\text{eff}}$ estimated under model $\mathcal{H}_{N}$ (which depends on $N$) and possibly on the number of scalar parameters $|\Omega|$ that make up the observation $\boldsymbol{Y}$.

Applying this methodology to EM-BiG-AMP, where $p_{\textsf{{{Y}}}|\mathsfbf{\Theta}_{N}}(\boldsymbol{Y}\,|\,\boldsymbol{\Theta}_{N})=p_{\textsf{{{Y}}}|\textsf{{{Z}}}}(\boldsymbol{Y}\,|\,\boldsymbol{A}_{N}\boldsymbol{X}_{N};\boldsymbol{\theta})$, we obtain the rank-selection rule

Since $N_{\text{eff}}$ depends on the application (e.g., matrix completion, robust PCA, dictionary learning), detailed descriptions of $\eta(\cdot)$ are postponed to .

To perform the maximization over $N$ in (113), we start with a small hypothesis $N_{1}$ and run EM-BiG-AMP to completion, generating the (approximate) MMSE estimates $\boldsymbol{\hat{A}}_{N_{1}},\boldsymbol{\hat{X}}_{N_{1}}$ and ML estimate $\boldsymbol{\hat{\theta}}$, which are then used to evaluateSince we compute approximate MMSE estimates rather than ML estimates, we are in fact evaluating a lower bound on the penalized log-likelihood. the penalized log-likelihood in (113). The $N$ hypothesis is then increased by a fixed value (i.e., $N_{2}=N_{1}+\texttt{rankStep}$), initializations of $(\boldsymbol{A}_{N_{2}},\boldsymbol{X}_{N_{2}},\boldsymbol{\theta})$ are chosen based on the previously computed $(\boldsymbol{\hat{A}}_{N_{1}},\boldsymbol{\hat{X}}_{N_{1}},\boldsymbol{\hat{\theta}})$, and EM-BiG-AMP is run to completion, yielding estimates $(\boldsymbol{\hat{A}}_{N_{2}},\boldsymbol{\hat{X}}_{N_{2}},\boldsymbol{\hat{\theta}})$ with which the penalized likelihood is again evaluated. This process continues until either the value of the penalized log-likelihood decreases, in which case $\widehat{N}$ is set at the previous (i.e., maximizing) hypothesis of $N$, or the maximum-allowed rank $\overline{N}$ is reached.

V-B2 Rank contraction

We now describe an alternative rank-selection procedure that is appropriate when Z has a “cliff” in its singular value profile and which is reminiscent of that used in LMaFit . In this approach, EM-BiG-AMP is initially configured to use the maximum-allowed rank, i.e., $N=\overline{N}$. After the first EM iteration, the singular values $\{\sigma_{n}\}$ of the estimate $\boldsymbol{\hat{X}}$ and the corresponding pairwise ratios $R_{n}=\sigma_{n}/\sigma_{n+1}$ are computed,In some cases the singular values of $\boldsymbol{\hat{A}}$ could be used instead. from which a candidate rank estimate $\widehat{N}=\operatorname*{arg\,max}_{n}R_{n}$ is identified, corresponding to the largest gap in successive singular values. However, this candidate is accepted only if this maximizing ratio exceeds the average ratio by the user-specified parameter $\tau_{\text{MOS}}$ (e.g., $\tau_{\text{MOS}}=5$), i.e., if

and if $\widehat{N}/\overline{N}$ is sufficiently small. Increasing $\tau_{\text{MOS}}$ makes the approach less prone to selecting an erroneous rank during the first few iterations, but making the value too large prevents the algorithm from detecting small gaps between the singular values. If $\widehat{N}$ is accepted, then the matrices A and X are pruned to size $\widehat{N}$ and EM-BiG-AMP is run to convergence. If not, EM-BiG-AMP is run for one more iteration, after which a new candidate $\widehat{N}$ is identified and checked for acceptance, and so on.

In many cases, a rank candidate is accepted after a small number of iterations, and thus only a few SVDs need be computed. This procedure has the advantage of running EM-BiG-AMP to convergence only once, rather than several times under different hypothesized ranks. However, when the singular values of Z decay smoothly, this procedure can mis-estimate the rank, as discussed in .

VI Matrix Completion

In this and the next two sections, we detail the application of BiG-AMP to the problems of matrix completion (MC), robust principle components analysis (RPCA), and dictionary learning (DL), respectively. For each application, we discuss the BiG-AMP’s choice of matrix representation, priors, likelihood, initialization, adaptive damping, EM-driven parameter learning, and rank-selection. Also, for each application, we provide an extensive empirical study comparing BiG-AMP to state-of-the-art solvers on both synthetic and real-world datasets. These results demonstrate that BiG-AMP yields excellent reconstruction performance (often best in class) while maintaining competitive runtimes. For each application of BiG-AMP discussed in the sequel, we recommend numerical settings for necessary parameter values, as well as initialization strategies when appropriate. Although we cannot guarantee that our recommendations are universally optimal, they worked well for the range of problems considered in this paper, and we conjecture that they offer a useful starting point for further experimentation. Nevertheless, modifications may be appropriate when applying BiG-AMP outside the range of problems considered here. Our BiG-AMP Matlab code can be found as part of the GAMPmatlab package at https://sourceforge.net/projects/gampmatlab/, including examples of BiG-AMP applied to the MC, RPCA, and DL problems.

As in several existing Bayesian approaches to matrix completion (e.g., ), we choose Gaussian priors for the factors A and X. Although EM-BiG-AMP readily supports the use of priors with row- and/or column-dependent parameters, we focus on simple iid priors of the form

where the mean and variance in (116) can be tuned using EM-BiG-AMP, as described in the sequel, and where the variance in (115) is fixed to avoid a scaling ambiguity between A and X. Section VI-F demonstrates that this simple approach is effective in attacking several MC problems of interest. Assuming the observation noise to be additive and Gaussian, we then choose the PIAWGN model from (80) for the likelihood $p_{\textsf{{{Y}}}|\textsf{{{Z}}}}$ given by

Note that, by using (115)-(116) with $\widehat{x}_{0}=0$ and the scalar-variance approximation from Sec. III-A, the BiG-AMP algorithm from Table III reduces to the simpler BiG-AMP-Lite algorithm from Table V with $\nu^{a}_{0}=1$.

VI-B Initialization

In most cases we advocate initializing the BiG-AMP quantities $\boldsymbol{\hat{X}}(1)$ and $\boldsymbol{\hat{A}}(1)$ using random draws from the priors $p_{\textsf{{{X}}}}$ and $p_{\textsf{{{A}}}}$, although setting either $\boldsymbol{\hat{X}}(1)$ or $\boldsymbol{\hat{A}}(1)$ at zero also seems to perform well in the MC application. Although it is also possible to use SVD-based initializations of $\boldsymbol{\hat{X}}(1)$ and $\boldsymbol{\hat{A}}(1)$ (i.e., for SVD $\boldsymbol{Y}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{D}^{T}$, set $\boldsymbol{\hat{A}}(1)=\boldsymbol{U}\boldsymbol{\Sigma}^{1/2}$ and $\boldsymbol{\hat{X}}(1)=\boldsymbol{\Sigma}^{1/2}\boldsymbol{D}^{T}$) as done in LMaFit and VSBL , experiments suggest that the extra computation required is rarely worthwhile for BiG-AMP.

As for the initializations $\nu_{nl}^{x}(1)$ and $\nu_{mn}^{a}(1)$, we advocate setting them at $10$ times the prior variances in (115)-(116), which has the effect of weighting the measurements $\boldsymbol{Y}$ more than the priors $p_{\textsf{{{X}}}},p_{\textsf{{{A}}}}$ during the first few iterations.

VI-C Adaptive damping

For the assumed likelihood (117) and priors (115)-(116), the adaptive-damping cost criterion $\widehat{J}(t)$ described in Sec. IV-B reduces to

To derive (118), one can start with the first term in (108) and leverage the Gaussianity of the approximated posterior on $\textsf{x}_{nl}$:

which then directly yields the first term in (118). The second term in (118) follows using a similar procedure, and the third and fourth terms follow directly from the PIAWGN model.

In the noise free setting (i.e., $\nu^{w}\rightarrow 0$), the third term in (118) dominates, omitting the need to compute the other terms.

VI-D EM-BiG-AMP

For the likelihood (117) and priors (115)-(116), the distributional parameters $\boldsymbol{\theta}=[\nu^{w},\widehat{x}_{0},\nu_{0}^{x}]^{T}$ can be tuned using the EM approach from Sec. V-A.For the first EM iteration, we recommend initializing BiG-AMP using $\nu_{nl}^{x}(1)=\nu^{x}_{0}$, $\widehat{x}_{nl}(1)=\widehat{x}_{0}$, $\nu_{mn}^{a}(1)=1$, and $\widehat{a}_{mn}(1)$ drawn randomly from $p_{\textsf{a}_{mn}}$. After the first iteration, we recommend warm-starting BiG-AMP using the values from the previous EM iteration. To initialize $\boldsymbol{\theta}$ for EM-BiG-AMP, we adapt the procedure outlined in to our matrix-completion problem, giving the EM initializations $\widehat{x}_{0}=0$ and

where $\textrm{SNR}^{0}$ is an initial estimate of the signal-to-noise ratio that, in the absence of other knowledge, can be set at $100$.

VI-E Rank selection

For MC rank-selection under the penalized log-likelihood strategy (113), we recommend using the small sample corrected AIC (AICc) penalty $\eta(N)=2\frac{|\Omega|}{|\Omega|-N_{\text{eff}}-1}N_{\text{eff}}$. For the MC problem, $N_{\text{eff}}=\mathsf{df}+3$, where $\mathsf{df}\triangleq N(M\!+\!L\!-\!N)$ counts the degrees-of-freedom in a rank-$N$ real-valued $M\times L$ matrix and the three additional parameters come from $\boldsymbol{\theta}$. Based on the PIAWGN likelihood (117) and the standard form of the ML estimate of $\nu^{w}$ (see, e.g., [55, eq. (7)]), the update rule (113) becomes

We note that a similar rule (but based on BIC rather than AICc) was used for rank-selection in .

MC rank selection can also be performed using the rank contraction scheme described in Sec. V-B2. We recommend choosing the maximum rank $\overline{N}$ to be the largest value such that $\overline{N}(M+L-\overline{N})<|\Omega|$ and setting $\tau_{\text{MOS}}=1.5$. Since the first EM iteration runs BiG-AMP with the large value $N=\overline{N}$, we suggest limiting the number of allowed BiG-AMP iterations during this first EM iteration to $\texttt{nitFirstEM}=50$. In many cases, the rank learning procedure will correctly reduce the rank after these first few iterations, reducing the added computational cost of the rank selection procedure.

VI-F Matrix Completion Experiments

We now present the results of experiments used to ascertain the performance of BiG-AMP relative to existing state-of-the-art algorithms for matrix completion. For these experiments, we considered IALM , a nuclear-norm based convex-optimization method; LMaFit , a non-convex optimization-based approach using non-linear successive over-relaxation; GROUSE , which performs gradient descent on the Grassmanian manifold; Matrix-ALPS , a greedy hard-thresholding approach; and VSBL , a variational Bayes approach. In general, we configured BiG-AMP as described in Sec. VIUnless otherwise noted, we used the BiG-AMP parameters $T_{\max}=1500$ (see Sec. II-H for descriptions) and the adaptive damping parameters $\texttt{stepInc}=1.1$, $\texttt{stepDec}=0.5$, $\texttt{stepMin}=0.05$, $\texttt{stepMax}=0.5$, $\texttt{stepWindow}=1$, and $\beta(1)=\texttt{stepMin}$. (See Sec. IV-B for descriptions). and made our best attempt to configure the competing algorithms for maximum performance. That said, the different experiments that we ran required somewhat different parameter settings, as we detail in the sequel.

Defining “successful” matrix completion as NMSE $<-100$ dB, Fig. 2 shows the success rate of each algorithm over a grid of sampling ratios $\delta\triangleq\frac{|\Omega|}{ML}$ and ranks $N$. As a reference, the solid line superimposed on each subplot delineates the problem feasibility boundary, i.e., the values of $(\delta,N)$ yielding $|\Omega|=\mathsf{df}$, where $\mathsf{df}=N(M+L-N)$ is the degrees-of-freedom in a rank-$N$ real-valued $M\times L$ matrix; successful recovery above this line is impossible by any method.

Figure 2 shows that each algorithm exhibits a sharp phase-transition separating near-certain success from near-certain failure. There we see that BiG-AMP yields the best PTC. Moreover, BiG-AMP’s PTC is near optimal in the sense of coming very close to the feasibility boundary for all tested $\delta$ and $N$. In addition, Fig. 2 shows that BiG-AMP-Lite yields the second-best PTC, which matches that of BiG-AMP except under very low sampling rates (e.g., $\delta<0.03$). Recall that the only difference between the two algorithms is that BiG-AMP-Lite uses the scalar-variance simplification from Sec. III-A.

Figure 3 plots median runtimeThe reported runtimes do not include the computations used for initialization nor those used for runtime evaluation. to NMSE $=-100$ dB versus rank $N$ for several sampling ratios $\delta$, uncovering orders-of-magnitude differences among algorithms. For most values of $\delta$ and $N$, LMaFit was the fastest algorithm and BiG-AMP-Lite was the second fastest, although BiG-AMP-Lite was faster than LMaFit at small $\delta$ and relatively large $N$, while BiG-AMP-Lite was slower than GROUSE at large $\delta$ and very small $N$. In all cases, BiG-AMP-Lite was faster than IALM and VSBL, with several orders-of-magnitude difference at high rank. Meanwhile, EM-BiG-AMP was about $3$ to $5$ times slower than BiG-AMP-Lite. Although none of the algorithm implementations were fully optimized, we believe that the reported runtimes are insightful, especially with regard to the scaling of runtime with rank $N$.

VI-F2 Approximately low-rank matrices

As in , we first tried to recover $\boldsymbol{Z}$ from the noiseless incomplete observations $\{z_{ml}\}_{(m,l)\in\Omega}$, with $\Omega$ chosen uniformly at random. Figure 4 shows the performance of several algorithms that are able to learn the underlying rank: LMaFit,LMaFit was run under the settings provided in their source code for this example. VSBL,VSBL was allowed at most $100$ iterations and run with DIMRED_THR $=10^{3}$, UPDATE_BETA $=1$, and tolerance $=10^{-8}$. and EM-BiG-AMP under the penalized log-likelihood rank selection strategy from Sec. V-B1.Rank-selection rule (113) was used with up to $5$ EM iterations for each rank hypothesis $N$, a minimum of $30$ and maximum of $100$ BiG-AMP iterations for each EM iteration, and a BiG-AMP tolerance of $10^{-8}$. All three algorithms were allowed a maximum rank of $\overline{N}=30$. The figure shows that the NMSE performance of BiG-AMP and LMaFit are similar, although BiG-AMP tends to find solutions with lower rank but comparable NMSE at low sampling ratios $\delta$. For this noiseless experiment, VSBL consistently estimates ranks that are too low, leading to inferior NMSEs.

Next, we examined noisy matrix completion by constructing the matrix $\boldsymbol{Z}=\boldsymbol{U\Sigma V}^{\textsf{T}}$ as above but then corrupting the measurements with AWGN. Figure 5 shows NMSE and estimated rank versus the measurement signal-to-noise ratio (SNR) $\sum_{(m,l)\in\Omega}|z_{ml}|^{2}/\sum_{(m,l)\in\Omega}|y_{ml}-z_{ml}|^{2}$ at a sampling rate of $\delta=0.2$. There we see that, for SNRs $<50$ dB, EM-BiG-AMP and VSBL offer similarly good NMSE performance and nearly identical rank estimates, whereas LMaFit overestimates the rank and thus performs worse in NMSE. Meanwhile, for SNRs $>50$ dB, EM-BiG-AMP and LMaFit offer similarly good NMSE performance and nearly identical rank estimates, whereas VSBL underestimates the rank and thus performs worse in NMSE. Thus, in these examples, EM-BiG-AMP is the only algorithm to successfully estimate the rank across the full SNR range.

VI-F3 Image completion

We now compare the performance of several matrix-completion algorithms for the task of reconstructing an image from a subset of its pixels. For this, we repeated the experiment in the Matrix-ALPS paper , where the $512\times 512$ boat image was reconstructed from $35\%$ of its pixels sampled uniformly at random. Figure 6 shows the complete (full-rank) image, the images reconstructed by several matrix-completion algorithmsAll algorithms were run with a convergence tolerance of $10^{-4}$. VSBL was run with $\beta$ hand-tuned to maximize performance, as the adaptive version did not converge on this example. GROUSE was run with $\texttt{maxCycles}=600$ and $\texttt{step\_size}=0.1$. Matrix-ALPS II with QR was run under default parameters and $300$ allowed iterations. Other settings are similar to earlier experiments. under a fixed rank of $N=40$, and the NMSE-minimizing rank-$40$ approximation of the complete image, computed using an SVD. In all cases, the sample mean of the observations was subtracted prior to processing and then added back to the estimated images, since this approach generally improved performance. Figure 6 also lists the median reconstruction NMSE over $10$ sampling-index realizations $\Omega$. From these results, it is apparent that EM-BiG-AMP provides the best NMSE, which is only $3$ dB from that of the NMSE-optimal rank-$40$ approximation.

VI-F4 Collaborative Filtering

In our final experiment, we investigate the performance of several matrix-completion algorithms on the task of collaborative filtering. For this, we repeated an experiment from the VSBL paper that used the MovieLens 100k dataset, which contains ratings $\{z_{ml}\}_{(m,l)\in\mathcal{R}}$, where $|\mathcal{R}|=100\,000$ and $z_{ml}\in\{1,2,3,4,5\}$, from $M=943$ users about $L=1682$ movies. The algorithms were provided with a randomly chosen training subset $\{z_{ml}\}_{(m,l)\in\Omega}$ of the ratings (i.e., $\Omega\subset\mathcal{R}$) from which they estimated the unseen ratings $\{\widehat{z}_{ml}\}_{(m,l)\in\mathcal{R}\setminus\Omega}$. Performance was then assessed by computing the Normalized Mean Absolute Error (NMAE)