On the Universality of Noiseless Linear Estimation with Respect to the Measurement Matrix

TL;DR

This paper demonstrates that the universality of noiseless linear estimation extends beyond random i.i.d. measurement matrices to structured matrices, including the optimal Bayesian reconstruction, using message passing methods.

Contribution

It shows that universality applies to the Bayes-optimal MMSE error and structured matrices, broadening the understanding of noiseless linear estimation.

Findings

Universality extends to Bayes-optimal MMSE error.

Universality applies to a range of structured matrices.

Message passing methods reveal these universal properties.

Abstract

In a noiseless linear estimation problem, one aims to reconstruct a vector x* from the knowledge of its linear projections y=Phi x*. There have been many theoretical works concentrating on the case where the matrix Phi is a random i.i.d. one, but a number of heuristic evidence suggests that many of these results are universal and extend well beyond this restricted case. Here we revisit this problematic through the prism of development of message passing methods, and consider not only the universality of the l1 transition, as previously addressed, but also the one of the optimal Bayesian reconstruction. We observed that the universality extends to the Bayes-optimal minimum mean-squared (MMSE) error, and to a range of structured matrices.

Equations

$$\mathbf{y}=\Phi {\mathbf{x}}^{\ast },$$

$$x_i\sim (1-\rho )\delta (x)+\rho \mathcal{N}(0,1).$$

$$\Phi =U\Sigma V$$

$$\mathbf{y}=\Phi \mathbf{x}=U\Sigma V\mathbf{x}.$$

$$\Sigma^{inv}=\left[\begin{array}[]{c}\tilde{\Sigma}^{-1}\\ \hline\cr 0\end{array}\right]\text{ such that }\Sigma^{inv}\Sigma=\left[\begin{array}[]{c|c}I_{m}&0\\ \hline\cr 0&0\end{array}\right].$$

$$\tilde{\mathbf{y}}={\Sigma }^{inv}{U}^T\mathbf{y}=\tilde{V}\mathbf{x}$$

$$\Sigma^{inv}=\left[\begin{array}[]{c}\tilde{\Sigma}\\ \hline\cr 0\end{array}\right]$$

$$\tilde{\mathbf{y}}={\Sigma }^{inv}{U}^T\mathbf{y}=V\mathbf{x}.$$

$${\hat{\mathbf{x}}}^{t+1}={\eta }_t({\Phi }^T{\hat{\mathbf{x}}}^t)$$

$${\mathbf{z}}^t=\mathbf{y}-\Phi {\hat{\mathbf{x}}}^t+\frac{1}{\alpha }{\mathbf{z}}^{t-1}⟨{\eta }_{t-1}^{\prime }({\Phi }^T{\mathbf{z}}^{t-1}+{\hat{\mathbf{x}}}^{t-1})⟩.$$

$${\sigma }_{t+1}^2=\Psi ({\sigma }_t^2),\Psi ({\sigma }^2)=\mathbb{E}\left[{\left({\eta }_t(X+\frac{\sigma }{\sqrt{\alpha }}Z)-X\right)}^2\right],$$

$$U_0{\Sigma }_0{\tilde{\Sigma }}^{-1}{U}^T\mathbf{y}$$

$${\mathbf{y}}^{\prime }$$

$${\mathbf{u}}_{\ell }^{t+1}=\frac{{{\hat{\mathbf{x}}}_{\mathbf{l}}}^t}{⟨{\mathrm{Var}}_{\ell }^t(\mathbf{x})⟩}-{\mathbf{u}}_r^t,{\rho }_{\ell }^{t+1}=\frac{1}{⟨{\mathrm{Var}}_{\ell }^t(\mathbf{x})⟩}-{\rho }_r^t,$$

$${\mathbf{u}}_r^{t+1}=\frac{{{\hat{\mathbf{x}}}_{\mathbf{r}}}^t}{⟨{\mathrm{Var}}_r(\mathbf{x})⟩}-{\mathbf{u}}_{\ell }^t,{\rho }_r^{t+1}=\frac{1}{⟨{\mathrm{Var}}_r^t(\mathbf{x})⟩}-{\rho }_{\ell }^t,$$

$${{\hat{\mathbf{x}}}_{\mathbf{l}}}^t$$

$$⟨{\mathrm{Var}}_{\ell }^t(\mathbf{x})⟩$$

$${({\hat{x}}_r)}_j$$

$$⟨{\mathrm{Var}}_r^t(\mathbf{x})⟩$$

$${\sigma }^t=⟨{\mathrm{Var}}_r^t(\mathbf{x})⟩\mathrm{and}{ϵ}^{\mathrm{t}}=⟨{\mathrm{Var}}_{\mathrm{l}}^{\mathrm{t}}(\mathbf{x})⟩,$$

$${\sigma }^t({\rho }_l^t)$$

$$ϵ({\rho }_r^t)$$

$$\underset{\Delta \to 0}{\lim }ϵ({\rho }_r^t)=\frac{1-\alpha }{{\rho }_r^t}.$$

$$Y_{jk}=\sqrt{\frac{2}{n}}{ϵ}_k\cos \left(\frac{\pi (2j+1)k}{2n}\right),$$

$$\left[\begin{array}{cc}H& H\\ H& -H\end{array}\right]$$

$$W_2=\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]\text{ and }{W}_{2k}=\left[\begin{array}{c}{H}_k\otimes [1,-1]\\ I_k\otimes [1,1]\end{array}\right]$$

Full text

1 Introduction

The problem of recovering a signal through the knowledge of its linear projections is ubiquitous in modern information theory, statistics and machine learning. In particular, many applications require to reconstruct an unknown $n-$dimensional signal vector ${\bf x^{*}}$ from the linear projections

[TABLE]

where ${\bf y}$ is a m-dimensional vector, and $\Phi$ is a $m\times n$ random matrix. For instance, if ${\bf x^{*}}$ is sparse, this task of estimating the signal from its linear random projections is at the roots of compressed sensing [1]. A fundamental question in the field is how much the algorithmic and the information theoretic performance depends on the choice of the random matrix $\Phi$.

In the present letter, we concentrate on the noiseless and asymptotic, large $n$, regime with a fixed value $\alpha\!=\!m/n$. We consider $\bf x^{*}$ to be $k$-sparse, i.e. to have only $k$ non-zero values, and we shall work in the limit where $n\to\infty$, $k\to\infty$, and a finite value of $\rho\!=\!k/n$. In such setting, a classical result is the following: for random matrices $\Phi$ with independent standard Gaussian entries, the (convex) reconstruction with $\ell_{1}$ penalty displays a precisely determined phase transition. For a certain region in the ($\alpha,\rho$)-phase diagram, it typically finds back the vector $\bf x^{*}$, being the sparsest solution, whereas outside that region, it typically fails. The obundary between these two regions is called the Donoho-Tanner line [2]. It has been shown empirically that the very same phase transition location seems to hold for a wider range of random matrix ensembles, see e.g. [3, 4], suggesting a large universality of the Donoho-Tanner phase transitions. Another line of work showed that the convex $\ell_{1}$ reconstruction problem can be treated through conic geometry, and the success probability of signal recovery only depends on a geometric number characterizing a subcone (statistical dimension or Gaussian width) [5, 6].

Here we investigate the universality of the phase transition not only for the $\ell_{1}$ transition, but also to the performance of the optimal Bayesian reconstruction. We analyze this question through the prism of information theory, message passing methods, and random matrix theory. We shall see that the universality indeed extends to a more generic set of properties than the $\ell_{1}$ transition, such as the minimum mean-squared (MMSE) error or the easy-hard phase transition for optimal Bayesian learning, and empirically to structured matrices such as the one appearing in [7, 8].

We note that investigation of universality are very common to physics problems, and understanding how large is the class of model for which a given result applied is a very fundamental question. The message-passing-based algorithm that we investigate in this paper to demonstrate the universality also has their origin in pysics works, such as [9].

2 A short review of results for i.i.d. random matrices

A first well-understood case of universality holds for random matrices $\Phi$ where all the elements are generated i.i.d. from a well-behaved distribution -with zero mean and unit variance- which all exhibit the same transitions as Gaussian random matrices. This is known for multiple retrieval problems:

2.1 $\ell_{1}$ recovery

Consider for instance the Donoho-Tanner line[2] that regulates the $\ell_{1}$ recovery. Thanks to the approximate message passing solver (see below) that has been shown to be universal with respect to all i.i.d. distributions with finite moments [10, 11], we know that the Donoho-Tanner phase transition is the same for all such random matrices.

2.2 Information theoretic optimal reconstruction

There has been a considerable amount of work in the information theory community on the computation of the mutual information and on the MMSE for problems such as (1) with Gaussian matrices. In particular, following the replica method from statistical physics (the Tanaka formula [12]), a heuristic formula has been postulated in different situations, see e.g. [13, 14, 15, 16]. This heuristic replica result has been recently rigorously proven in a series of papers [17, 17, 18]. In a more recent proof [19], it has been shown, again, that the formula is not specific to Gaussian i.i.d. matrices, but that any matrix with i.i.d elements of unit variance and zero mean leads to the same exact result for the mutual information and the MMSE.

2.3 Hard phase for Bayesian decoders

A third interesting point is to ask about tractacle decoders that aim to perform the optimal Bayesian estimation, i.e. with a perfect prior knowledge on the distribution of $\bf x^{*}$. For simplicity, consider for instance the case where each element of $\bf x^{*}$ has been sampled from a Gauss-Bernoulli distribution:

[TABLE]

In this case, the best known solver is again AMP, using a Bayesian decoder (instead of the soft thresholding function for $\ell_{1}$ recovery) [20, 21, 14, 15]. Interestingly, it shares with the $\ell_{1}$ recovery a similar phase transition: for a certain region in the ($\alpha,\rho$) plane it typically finds back the vector $\bf x^{*}$, whereas outside that region it fails. We shall denote the limit between these regions the "Bayesian hard-phase" transition. The "Bayesian hard-phase" line, that has been precisely computed in [14, 15] is always better than the Donoho-Tanner line (as it should, since it exploits additional information). Once more, the universality of AMP shows that this phase transition is not restricted to Gaussian matrices, but extends as well to all (well normalized) i.i.d. matrices.

The fact that these three properties (the $\ell_{1}$, the hard-phase line, as well as the MMSE) are universal for all i.i.d. matrices makes the case for Gaussian computations, as done in theoretical computation, stronger. We shall see that this universality extends well beyond these simple cases.

3 Random rotationally invariant matrices

Moving away from the well-known i.i.d. examples, we start by considering a much larger set of random matrices defined through their singular value decomposition (SVD): any real matrix $\Phi$ can be decomposed into ${\Phi=U\Sigma V}$, with $U$ and $V$ orthogonal matrices, and $\Sigma$’s elements being $\Phi$’s singular values. We shall look at the left rotationally invariant random matrix ensemble: these are matrices $\Phi$ that can be written as

[TABLE]

with an arbitrary rotation matrix $U$ and singular values $\Sigma$, but where the matrix $V$ has been randomly (and independently of $\Sigma$ and $U$) generated from the Haar measure (that is, uniformly from all possible rotations).

When the singular values are different from zero, it is straightforward to justify the universality property for matrices from this subclass. We start by the definition of the problem: we wish to find $\bm{x}$ such that

[TABLE]

If $m\leq n$, then $\Sigma$ is written as \Sigma=\left[\begin{array}[]{c|c}\tilde{\Sigma}&0\\ \end{array}\right] and we define

[TABLE]

Multiplying (2) on both sides by $U^{T}$, and then by $\Sigma^{inv}$; one reaches

[TABLE]

where $\tilde{V}$ is a $m\times n$ matrix composed of the first $m$ lines of $V$. If instead $m>n$, $\Sigma$ is written as

[TABLE]

and we define \Sigma^{inv}=\left[\begin{array}[]{c|c}\tilde{\Sigma}^{-1}&0\\ \end{array}\right] such that $\Sigma^{inv}\Sigma=I_{n}$. Multiplying (2) by $U^{T}$ then $\Sigma^{inv}$, we obtain

[TABLE]

In both cases, we thus see that the problem has been transformed —in a constructive way— into a standard linear system with the sensing matrix $\tilde{V}$ when $m\leq n$ being a (sub-sampled) random rotation one, or sensing matrix $V$ when $m>n$. This shows that all rotationally invariant matrices, which satisfy $U$ and $\Sigma$’s independence on $V$, can be transformed the same way and are in the same universality class as far as noiseless linear recovery is concerned, i.e. they will display the same phase transitions.

Since Gaussian i.i.d. matrices belong among random rotationally invariant matrices (in this case $\Sigma$ follows the Marcenko-Pastur law [22]) this means that all the information theoretic rigorous results (such as phase transitions and MMSE value) with zero noise for random Gaussian matrices applies verbatim to all rotationally invariant ensemble, as long as the SVD’s matrices $U$ and $\Sigma$ are independent of $V$. This is a very strong universality, that applies to the three cases (1, 2, 3) from sec. 2. Note that the universality of the Donoho-Tanner line with rotationally invariant matrices was already hinted by the replica method [23].

Notice, however, that the above construction depends crucially on the fact that we consider here noiseless measurements. It would not work if an additional Gaussian noise were added in eq. (1): in this case, the transformation would make the i.i.d. Gaussian noise a correlated one. Indeed, the replica formula for noisy measurements underlines that the MMSE depends on the precise set of matrices in noisy reconstruction [24, 13] (this formula is not yet fully rigorous, but see [25] for a proof in a restricted setting). Any differences, however, must go to zero in the noiseless limit.

4 Approximate Message Passing

Having discussed the universality with respect to random rotationally invariant matrices, we now wish to discuss its effect on specific solvers, concretely the message passing algorithms.

4.1 AMP

We first consider the original approximate message passing (AMP) [26] to compute the phase transition between the phase where the algorithm reconstructs $\bf x^{*}$ perfectly, and the one where reconstruction may be possible but is not achieved by the algorithm. AMP is an iterative algorihm that follows:

[TABLE]

where $t$ is the iteration index, ${\bf x}^{t}$ is the current estimate of ${\bf x^{*}}$, ${\bf z}^{t}$ the current residual, $\langle\cdot\rangle$ is an averaged sum of components, and $\eta_{t}$ is a prior-dependent threshold function applied component-wise (the soft thresholding for $\ell_{1}$, or the Bayesian decoder [14, 15]).

One of the most interesting features of AMP is that, if $\Phi$ is a Gaussian i.i.d. matrix, its mean squared error (MSE) $\sigma_{t}$ can be tracked accurately by the state evolution formalism [26, 10, 11]. State evolution is a relatively simple recursive equation:

[TABLE]

where the expectation is with respect to independent random variables $Z\sim{\cal N}(0,1)$ and $X$, whose distribution coincides with the empirical distribution of the entries of $x^{*}$. Analyzing the evolution of this equation for the $\ell_{1}$ decoder yields the Donoho-Tanner line[26], while using the Bayesian decoder it yields the hard-phase line for Bayesian decoding [14].

It would be interesting to use AMP for rotationally invariant matrices. In order to do this, we follow the construction of sec. 3: starting from equation (3) we then multiply by $\Sigma_{0}$, a $m\times m$ diagonal matrix with singular values sampled from Marcenko-Pastur law (singular values of a Gaussian i.i.d. matrix 111The singular values of a Gaussian matrix are correlated, so in fact we may want to generate $\Sigma_{0}$ by first generating a random Gaussian matrix, and then calculating its singular values.), and $U_{0}$ a $m\times m$ Haar-generated orthogonal matrix, thus ensuring that $\Sigma_{0}$ and $U_{0}$ are generated independently of $V$:

[TABLE]

After this transformation, $\Phi^{\prime}=U_{0}\Sigma_{0}\tilde{V}$ is a random matrix that belongs to an ensemble very close to the Gaussian i.i.d. matrices ensemble. In fact, a recent work showed that AMP applied to a Gaussian matrix follows the same state evolution as matrices such as $\Phi^{\prime}$ where $U_{0},\tilde{V}$ are uniform orthogonal matrices and $\Sigma_{0}$ diagonal’s elements are singular values sampled from the Marcenko-Pastur law [27]. Combining this result with the matrix transformation, we have thus constructively mapped the noiseless reconstruction problem back to the well-understood noiseless compressed sensing case for a Gaussian i.i.d. matrix, where we can safely apply the algorithm, and its state evolution. In the section 5.2, we apply this matrix transformation for numerical experiments using AMP.

4.2 Vector-AMP

While the transformation trick allows to make AMP work with random rotationally invariant matrices, another alternative is to work directly with a dedicated solver. To this means, different but related approaches were proposed [24, 28], in particular, using the general expectation-propagation (EP) [29, 30] scheme. Ma and Ping proposed a variation of EP called OAMP [31] specially adapted to rotation matrices. Rangan, Schniter and Fletcher introduced a similar approach called VAMP [32] and proved that it follows state evolution equations corresponding to the fixed point of the replica potential [24, 13, 25]. The multi-layer AMP algorithm of[33] also display the same fixed point.

We shall concentrate here on the VAMP (Vector-AMP) approach, and for a moment, put back a small additional random Gaussian i.i.d. noise of variance $\Delta$ in the measurement in eq. (1) as it is needed for stating the algorithm. VAMP then consists in the following fixed-point iteration:

[TABLE]

where we denote by $\mathbb{E}_{\ell,r}^{t}$ the expectation w.r.t. the tilted distributions ${\tilde{Q}_{\ell,r}^{t}(\bm{x})\propto P_{\ell,r}(\bm{x})Q_{\ell,r}^{t}(\bm{x})}$, and by $\operatorname{Var}_{r,\ell}^{t}(\bm{x})$ the variance of these distributions. Here, we have defined $Q_{l,r}({\bf x})=e^{-\frac{1}{2}\rho_{l,r}\bm{x}^{T}\bm{x}+\bm{u}_{l,t}^{T}\bm{x}}$, $P_{l}(\bm{x})\propto e^{-||{\bf y}-\Phi{\bf x}||_{2}^{2}/2\Delta}$ and $P_{r}(\bm{x})$ is the prior used in the algorithm (i.e. the Laplace prior for the $\ell_{1}$ model, or the actual distribution of the signal for Bayesian reconstruction). In particular

[TABLE]

where, as for AMP, we define the denoiser that yields the estimates of $x$ by $z(u,\rho)=\int dxP_{r}(x)e^{-\frac{1}{2}\rho x^{2}+ux}$,

[TABLE]

Again, the performance of the recursion can be analyzed rigorously through the state evolution [32]. For simplicity, let us concentrate on the Bayes optimal case in which case the state evolution can be closed on the variables (see [32]):

[TABLE]

by writing

[TABLE]

where the expectation is above the distribution of the singular values $\Sigma$ of the matrix $\Phi$, and where we recognize the Stieljes transform $S_{X}(r)={\mathbb{E}}\left[1/{X-r}\right]$.

Though this transform, we see that the performance depends crucially on the distribution of eigenvalues. Let us now go back on the noiseless limit when $\Delta\to 0$ and analyze how the universality shows up. Consider again the Stieltjes transform: out of the $n$ singular values of the $n\times n$ matrix $\bm{\Phi}^{T}\bm{\Phi}$, we shall have $(1-\alpha)n$ of them to be zero (assuming $\alpha<1$) while the rest are positive (since $m<n$). In this case, the limit $r\to 0$ of the Stieltjes transform will behave as $S_{X}(r)\approx-(1-\alpha)/r$ so that

[TABLE]

Again, we see that all the complicated dependence on the spectrum of the matrix $\Phi$ has been eliminated. This is a direct, alternative, proof that VAMP will also yield universal results in the zero noise limit for the Bayesian reconstruction. Given that VAMP has the same fixed point as the replica mutual information [13, 25], this argument applies to the replica prediction for the MMSE as well.

5 Structured matrices

We now move to very structured matrices, in order to test the universality as well as the quality and the prediction of the state evolution out of its comfort zone. In order to do so, we have considered different matrix ensembles:

5.1 Tested ensembled of matrices

Discrete cosine transform matrices

The first ensemble we consider consists in Fourier-like matrices. A $n\times n$ discrete cosine transform (DCT) matrix $Y$ is defined by:

[TABLE]

where $j,k\in\llbracket 0,n-1\rrbracket$, $\epsilon_{0}=1/\sqrt{2}$, $\epsilon_{i}=1$ for $i=1,...,n-1$. We used a sub-sampled version of these matrices in which we picked some rows randomly.

Hadamard matrices

A natural variant of DCT is given by the Hadamard matrices. $H$ is a $n\times n$ Hadamard matrix if its entries are $\pm 1$ and its rows are pairwise orthogonal, i.e. $HH^{T}=nI_{n}$. For every integer $k$, there exists a Hadamard matrix $H_{k}$ of size $2^{k}$. These can be created with Sylvester’s construction: Let $H$ be a Hadamard matrix of order $n$. Then the partitioned matrix

[TABLE]

is a Hadamard matrix of order $2n$.

Random features maps

Finally, we wanted to consider here random features maps (RFM) as encountered in nonlinear regression problems. In such settings, a random features matrix $\Phi=f(WX)$ is obtained from the raw data matrix $X$ by means of a random projection matrix $W$ and a pointwise nonlinear activation $f$. Kernel regression models, nonlinear in the original data $X$, can then be approximately but efficiently solved by the linear estimation problem (1), with an appropriate choice for $f$ and the $W$-distribution [34]. Such matrices, that can be seen as the ouput of a neuron with random weights, have been investigated in particular in the context of neural networks [7, 8]. Indeed, in neural networks configurations with random weights play an important role as they define the initial loss landscape. They are also fundamental in the random kitchen sinks algorithm in machine learning [34] and it is thus of interest to test our understanding of linear reconstructions with AMP and VAMP in this case.

In what follows we will test random features matrices where both $W$ and $X$ are random Gaussian i.i.d. matrices.

5.2 Numerical results

We provide the codes used to generate the data on github in the repo http://sphinxteam/Universality-CS-2019. To generate Figure 1 and 2, we ran VAMP $50$ times on $50\times 50$ points spanning the ($\alpha$,$\rho$)-space, and computed the average mean-squared error (MSE) between the signal ${\bf x^{*}}$ and the reconstructed configuration x. The MSE is represented with a color bar (white means perfect reconstruction). For a DCT and a Hadamard matrix, we observe a phase transition in the Bayes-optimal case that matches the theoretical transition for Gaussian i.i.d. matrices. We also ran VAMP for the $\ell_{1}$ reconstruction problem. Averaging on 20 executions (or 50 for small $\alpha$ where finite-size effects are more important), we recover again a phase transition matching the theoretical Donoho-Tanner line for Gaussian i.i.d. matrices [3]. Besides, we compared the MSE obtained by VAMP at each point of the phase diagram for different matrices. In figures 3 and 4, we plot the MSE averaged on 20 executions for $\rho$ fixed and $\alpha$ ranging between 0 and 1. We get the same error in reconstruction for all matrices, following the MSE for Gaussian i.i.d. matrix for $\rho=0.25$, $0.5$ and $0.75$. We also checked that AMP, provided one uses the trick eq. (7), reproduce these results as well: indeed the two algorithms returned extremely similar results.

5.3 Discussion

Figures of the previous section perfectly illustrate our main point: the universality in noiseless compressed sensing is not limited to the $\ell_{1}$-type reconstruction as in [3, 4], but extends to other quantities and estimators, such as the hard-phase line in Bayesian reconstruction, and the MMSE. Besides, it is not limited to random orthogonal matrices, but empirically extends to Fourier-type matrices and to the random features maps currently studied in machine learning. It is an open question to extend the proof of state evolution to these challenging matrices. It would be interesting to find a good creterion to identify which matrices satisfy this universality and which do not; this is something that we are yet unable to predict in advance. An example of structured matrices that do not seem to follow these universal phase transitions is given by Haar wavelet matrices, which can be defined recursively by:

[TABLE]

where $I_{k}$ is the identity matrix of size $k$ and $\otimes$ is the Kronecker product. In fact, VAMP even fails to converge for these matrices. Investigating this behavior is an interesting direction of research.

Acknowledgment

We thank Andre Manoel and Galen Reeves for useful discussions. We acknowledge funding from the ERC under the European Union’s Horizon 2020 Research and Innovation Program Grant Agreement 714608-SMiLe; and from the French National Research Agency (ANR) grant PAIL.