Role of Bootstrap Averaging in Generalized Approximate Message Passing

Consider estimating an unknown signal $\mathbf{\boldsymbol{w}}_{0}\in\mathbb{R}^{N}$ from a random linear measurement $\mathbf{\boldsymbol{y}}\in\mathbb{R}^{M}$ in the form,

In (1), $X\in\mathbb{R}^{M\times N}$ is a known measurement matrix whose elements are independent and identically distributed (i.i.d.) standard Gaussian variables, and $\mathbf{\boldsymbol{\epsilon}}\sim\mathcal{N}(\mathbf{\boldsymbol{0}},\Delta I_{M})$ is the measurement noise. We also assume that each element of the unknown signal $\mathbf{\boldsymbol{w}}_{0}$ is i.i.d., according to a distribution $q_{0}$.

Generalized approximate message passing (GAMP) [1, 2] is a computationally efficient algorithm for solving this problem. A striking feature of GAMP is its applicability to various hypothesis testing [3]. Specifically, GAMP can provide an unbiased estimator $\hat{\mathbf{\boldsymbol{r}}}^{\rm G}\sim\mathcal{N}(\mathbf{\boldsymbol{w}}_{0},\hat{s}^{2}I_{M})$ [2] in a high-dimensional asymptotic setting with $M,N\to\infty,M/N\to\alpha\in(0,\infty)$, where the variance $\hat{s}^{2}$ depends on the quality of the measurement $\mathbf{\boldsymbol{y}}$ and denoising function used in GAMP. This unbiased estimator has been used to test the significance of estimated signals [4, 5, 3, 6].

The statistical power of these tests depends on the variance $\hat{s}^{2}$, with a lower variance leading to a higher statistical power. However, reducing the variance $\hat{s}^{2}$ is not a trivial task. Replacing the denoising function used in GAMP with a powerful one based on nonconvex regularization, for example, can worsen the convergence of GAMP [7] or, even if convergence is achieved, the improvement may be insignificant [8]. This study aims to find an alternative way to reduce the variance without nonconvex regularization.

To reduce the variance, we use the bootstrap averaging [9] of computational statistics (commonly known as ensemble learning in machine learning [10, 11]). Specifically, we consider averaging the unbiased estimators of GAMP for multiple bootstrap samples with arbitrary size $M\mu_{B},\,\mu_{B}\in(0,\infty]$. For denoising function of GAMP, we consider the one for the elastic net [12].

However, the efficient computation and theoretical analysis of an averaged unbiased estimator remains unresolved. To resolve this problem, we use AMP with resampling (AMPR) [13, 14] that has been proposed for computing bootstrap statistics of the elastic net estimator by running a variant of GAMP once. We will argue that AMPR is actually computing the bootstrap average of the unbiased estimators of GAMP. That is, the averaged unbiased estimator can be computed efficiently, and its variance can be analyzed using the state evolution (SE) of AMPR, which has been developed to analyze the performance of AMPR in [14, 13]. We then conduct a thorough numerical analysis of the SE of the AMPR to investigate when bootstrap averaging reduces the variance of the unbiased estimator, and what phenomena occur when optimizing the bootstrap sample size and the hyperparameter of the elastic net regularization.

The findings of this study are summarized as follows:

• •

  As mentioned above, the averaged unbiased estimator is obtained by AMPR. Furthermore, its variance can be estimated using the output of AMPR without knowing the actual signal $\mathbf{\boldsymbol{w}}_{0}$. Thus, we can minimize the variance by adjusting the bootstrap sample size and the hyperparameters of the elastic net (see Section III and IV)

• •

  The variance of the averaged unbiased estimator can be reduced via bootstrapping, especially when the true data generation process is inconsistent with the sparsity assumption of the regularization and the data size is insufficient (see Section V-B)

• •

  When $\mathbf{\boldsymbol{w}}_{0}$ is less sparse, and the data size is small, a phase transition occurs. This phase transition indicates the existence of the region where the value of the regularization parameter is infinitesimally small, and the number of unique data points in each bootstrap sample is less than the dimension of $\mathbf{\boldsymbol{w}}_{0}$. That is, in this region, the ensemble average of the unbiased estimators of GAMP for the elastic norm minimization problem (also known as the minimum norm interpolation in machine learning [15, 16, 17]) yields the best averaged unbiased estimator (see Section V-C)

Here, we explain the meaning of the SE of AMPR. Subsequently, we argue that $\hat{\mathbf{\boldsymbol{r}}}_{t}$ of AMPR is the bootstrap average of the unbiased estimators of GAMP.

The first point is the meaning of $\xi$ and $\eta$ that appear in SE. Propositions 1 and 2 indicate that, once the AMPR reaches its fixed point at sufficiently large $T_{\rm it}$, for any functions $\phi,\psi:\mathbb{R}\to\mathbb{R}$ such that the following expression is well defined, the following holds in the asymptotic setting $M,N\to\infty,M/N\to\alpha\in(0,\infty)$:

where the variables $\hat{Q},\hat{\chi},\hat{v}$ are the outputs of SE at a fixed point. In (11), one can interpret that $\xi$ and $\eta$ effectively represents the randomness originating from $(X,\mathbf{\boldsymbol{y}})$ and $\mathbf{\boldsymbol{c}}$, respectively.

The next point is the relationship between AMPR and GAMP. For this, we define $\bar{q}_{t}\equiv\mathcal{E}_{t}+v_{t}$ and $\bar{\chi}_{t}=\hat{\chi}_{t}+\hat{v}_{t}$. Then, the evolution of $(\hat{Q}_{t},\bar{\chi}_{t},\bar{q}_{t},\chi_{t})$ is equivalent to the SE of GAMP for computing $\hat{\mathbf{\boldsymbol{w}}}^{\ast}$ for one realization of $\mathbf{\boldsymbol{c}}$ [1, 2]. Recall that $\hat{\mathbf{\boldsymbol{w}}}^{\ast}$ is specified in (8). In addition, $\bar{\chi}_{t}/\hat{Q}_{t}^{2}\,(=(\mathcal{E}_{t}+\hat{v}_{t})/\hat{Q}_{t}^{2})$ represents the variance of the unbiased estimate of GAMP $\hat{\mathbf{\boldsymbol{r}}}_{t}^{\rm G}(\mathbf{\boldsymbol{c}})$ at each iteration. These observations indicate that the unbiased estimator of GAMP at each iteration $t$ can be decomposed into two Gaussian variables:

where $\mathbf{\boldsymbol{\xi}}_{0},\mathbf{\boldsymbol{\xi}},\mathbf{\boldsymbol{\eta}}\sim\mathcal{N}(0,I_{N})$, and $\mathbf{\boldsymbol{\xi}}$ and $\mathbf{\boldsymbol{\eta}}$ represent the randomness originating from $(X,\mathbf{\boldsymbol{y}})$ and $\mathbf{\boldsymbol{c}}$, respectively. We can interpret that AMPR computes the bootstrap statistics by decomposing $\mathbf{\boldsymbol{\xi}}_{0}$ into $\mathbf{\boldsymbol{\xi}}$ and $\mathbf{\boldsymbol{\eta}}$, and averaging over $\mathbf{\boldsymbol{\eta}}$. This yields the unbiased estimator $\hat{\mathbf{\boldsymbol{r}}}_{t}$ as in (10) that is the bootstrap average of $\hat{\mathbf{\boldsymbol{r}}}^{\rm G}(\mathbf{\boldsymbol{cka}})$. We verify this point numerically in section V-A (see Fig. 1 and 2).

In the sequel, by numerically minimizing the variance using SE, we investigate when bootstrapping reduces the variance of the unbiased estimator and the phenomena that occur when optimizing the bootstrap sample size $\mu_{B}$ and the hyperparameter of the elastic net regularization $(\lambda,\gamma)$.

For this, we searched for the optimal parameter $(\mu_{B}^{\star},\lambda^{\star},\gamma^{\star})$ that yielded the minimum variance using the SE of AMPR and the Nelder-Mead algorithm in the Optim.jl library [18]. We obtained the fixed point of the SE by iterating the SE a sufficient number of times. For comparison, the same optimization was performed for the non-bootstrap case. For the signal distribution $q_{0}$, we consider the Gauss-Bernoulli model: $q_{0}=\rho\delta_{0}+(1-\rho)\mathcal{N}(0,1)$, with sparsity $\rho\in(0,1)$.

In this section, we denote the outputs of the AMPR or GAMP at fixed points by unindexed variables.

In this study, we investigated the behavior of the bootstrap-averaged unbiased estimator of GAMP using AMPR and its SE. We found that the bootstrap averaging procedure can effectively reduce the variance of the unbiased estimator when the actual data generation process is inconsistent with the sparsity assumption of the regularization and the data size is insufficient. We also found a phase transition where the regularization strength drops to infinitesimally small values by decreasing the measurement ratio $\alpha$ or increasing $\rho$.

Although increasing the variance of weak learners is key to the success of ensemble learning [22, 23, 24], the phase transition to an ensemble of interpolators may be unexpected. Investigating whether similar phase transitions occur in other more sophisticated machine-learning models, such as neural networks, would be an interesting future direction.

We also remark that the landscape in $(\mu_{B},\lambda)$-space may not be convex in general. Thus, the uniqueness of the optimizer should be investigated in the future, paying attention to the relationship with the implicit regularization [25, 26, 27].

On the technical side, the key to this study was the precise performance characterization of the averaged estimator by AMPR, which was developed by combining GAMP and the replica method of statistical physics [28, 29]. Such a combination of approximate inference algorithms and the replica method has been used to develop approximate computation algorithms [13, 30, 14, 31, 32, 33, 34] and has not been applied to precise performance analysis of ensemble methods. It would be an interesting direction to try similar performance analysis for other bootstrap methods or ensemble learning.

This study was supported by JSPS KAKENHI Grant Number 21K21310 and 17H00764.