Average case analysis of Lasso under ultra-sparse conditions

TL;DR

This paper provides an average-case analysis of Lasso in ultra-sparse linear models using a novel replica method approach, offering insights into support recovery and performance without restrictive assumptions.

Contribution

It introduces a new analytical framework for Lasso's performance in ultra-sparse regimes, extending previous results to more general settings and noise conditions.

Findings

Provides a lower bound on sample complexity for support recovery

Generalizes previous bounds to non-Gaussian noise

Supports analysis with extensive numerical experiments

Abstract

We analyze the performance of the least absolute shrinkage and selection operator (Lasso) for the linear model when the number of regressors $N$ grows larger keeping the true support size $d$ finite, i.e., the ultra-sparse case. The result is based on a novel treatment of the non-rigorous replica method in statistical physics, which has been applied only to problem settings where $N$ ,$d$ and the number of observations $M$ tend to infinity at the same rate. Our analysis makes it possible to assess the average performance of Lasso with Gaussian sensing matrices without assumptions on the scaling of $N$ and $M$, the noise distribution, and the profile of the true signal. Under mild conditions on the noise distribution, the analysis also offers a lower bound on the sample complexity necessary for partial and perfect support recovery when $M$ diverges as $M = O(\log N)$. The obtained bound for perfect support recovery is a generalization of that given in previous literature, which only considers the case of Gaussian noise and diverging $d$. Extensive numerical experiments strongly support our analysis.

Tables

Table 1. Table 1: Values of α C subscript 𝛼 𝐶 \alpha_{C} evaluated from figure 2.

${𝒙}_S^0$	$(\lambda ,{\sigma }^2)$	${\alpha }_C$
$[1,1,1]$	$(0.5,0.0)$	6.00
$[1,1,1]$	$(0.5,0.5)$	10.0
$[1/3,2/3,1]$	$(0.5,0.0)$	4.89
$[1/3,2/3,1]$	$(0.5,0.5)$	8.89

Equations

$$\mathbf{y}=\mathbf{A}{\mathbf{x}}^0+\mathbf{\xi },$$

$${\hat{\mathbf{x}}}_{\lambda }(\mathbf{A},\mathbf{y}):=\underset{\mathbf{x}}{\arg \min }\text{quantity}(\frac{1}{2}\text{norm}{\mathbf{A}\mathbf{x}-\mathbf{y}}^2+M\lambda \text{norm}{\mathbf{x}}_1),$$

$${\mathbb{E}}_{\mathbf{A},\mathbf{y}}(\cdots )$$

$$=\int \text{differential}\mathbf{y}\text{differential}\mathbf{A}\text{differential}\mathbf{\xi }{p}_{\xi }(\mathbf{\xi })(\cdots )\frac{e^{-\frac{1}{2}\mathrm{Tr}{\mathbf{A}}^{\mathsf{T}}\mathbf{A}}}{(2\pi {)}^{(NM/2)}}\delta (\mathbf{y}-\mathbf{A}{\mathbf{x}}^0-\mathbf{\xi }),$$

$$\mathrm{FP}(\mathbf{A},\mathbf{y})$$

$$\mathrm{TP}(\mathbf{A},\mathbf{y})$$

$$\begin{array}{c}{P}_{\beta }(\mathbf{x}|\mathbf{A},\mathbf{y})\\ :=Z_{\beta }^{-1}\text{quantity}(\mathbf{A},\mathbf{y})\exp \text{quantity}(-\frac{\beta }{2}\text{norm}{\mathbf{A}\mathbf{x}-\mathbf{y}}^2-\beta M\lambda \text{norm}{\mathbf{x}}_1),\end{array}$$

$$\mathcal{F}=-\underset{\beta \to \infty }{\lim }{\beta }^{-1}{\mathbb{E}}_{\mathbf{A},\mathbf{y}}\log Z_{\beta }(\mathbf{A},\mathbf{y}).$$

$${\mathbb{E}}_{\mathbf{A},\mathbf{y}}\log Z_{\beta }(\mathbf{A},\mathbf{y})=\underset{n\to +0}{\lim }\frac{{\mathbb{E}}_{\mathbf{A},\mathbf{y}}{Z}_{\beta }^n(\mathbf{A},\mathbf{y})-1}{n}.$$

$$\begin{array}{c}{\mathbb{E}}_{\mathbf{A},\mathbf{y}}{Z}_{\beta }^n(\mathbf{A},\mathbf{y})={\mathbb{E}}_{\mathbf{A},\mathbf{y}}\int \prod _{a=1}^n\text{differential}{\mathbf{x}}^a\\ \exp \text{quantity}(-\frac{\beta }{2}\sum _{a=1}^n\text{norm}{\mathbf{A}{\mathbf{x}}^a-\mathbf{y}}^2-\beta M\lambda \sum _{a=1}^n\text{norm}{{\mathbf{x}}^a}_1).\end{array}$$

$$\sum _{i\notin S}{\Delta }_i^a{\Delta }_i^b:=\{\begin{array}{ll}Q& a=b\\ Q-\chi /\beta & \text{otherwise}\end{array},$$

$${\mathbb{E}}_{{\mathbf{A}}_S,\mathbf{y}}\int \prod _{a=1}^n\text{differential}{\mathbf{\Delta }}^a\int \text{differential}Q\text{differential}\chi e^{-\beta M\lambda {\sum }_{a=1}^n\text{norm}{{\mathbf{\Delta }}^a}_1}\mathcal{I}\mathcal{L},$$

$$\mathcal{I}:=\prod _{a=1}^n\delta \text{quantity}(Q-\sum _{i\notin S}({\Delta }_i^a{)}^2)\prod _{a<b}\delta \text{quantity}(Q-\frac{\chi }{\beta }-\sum _{i\notin S}{\Delta }_i^a{\Delta }_i^b),$$

$$\begin{array}{c}\mathcal{L}:=\int D\mathbf{z}\text{quantity}(\int \text{differential}{\mathbf{x}}_S{e}^{-M\beta G({\mathbf{x}}_S;\mathbf{z})}{)}^n,\\ G({\mathbf{x}}_S;\mathbf{z}):=\frac{\text{norm}{{\mathbf{A}}_S{\mathbf{x}}_S+\sqrt{Q}\mathbf{z}-\mathbf{y}}^2}{2M(1+\chi )}+\lambda \text{norm}{{\mathbf{x}}_S}_1.\end{array}$$

$$\begin{array}{c}\mathcal{I}={\int }_{-\mathrm{i}\infty }^{+\mathrm{i}\infty }\text{differential}\hat{Q}\text{differential}\hat{\chi }{e}^{\frac{Mn\beta }{2}\text{quantity}(Q\hat{Q}+(n-1)\chi \hat{\chi }-n\beta Q\hat{\chi })}\\ \times \int D\hat{\mathbf{z}}{e}^{-\frac{M\beta \hat{Q}}{2}{\sum }_{a=1}^n\text{norm}{{\mathbf{\Delta }}^a}^2+\beta \sqrt{M\hat{\chi }}{\mathbf{z}}^{\mathsf{T}}{\mathbf{\Delta }}^a+o(\beta )}.\end{array}$$

$$\begin{gathered}\mathcal{F}=\mathbb{E}_{\mathbf{A}_{S},\bm{y}}\ \operatorname*{{\rm Extr}}_{\Theta}\Bigg{\{}-\frac{Q\hat{Q}-\chi\hat{\chi}}{2}\\ -\frac{\tilde{N}}{2\hat{Q}}\quantity[(\Lambda+\hat{\chi}){\rm erfc}\quantity(\sqrt{\frac{\Lambda}{2\hat{\chi}}})-\sqrt{\frac{2\Lambda\hat{\chi}}{\pi}}e^{-\Lambda/2\hat{\chi}}]\\ +\int D\bm{z}\min_{\bm{x}_{S}}G(\bm{x}_{S};\bm{z})\Bigg{\}}.\end{gathered}$$

$$Q$$

$$\chi$$

$$\hat{Q}$$

$$=\frac{M-\int D\mathbf{z}\text{norm}{{\hat{\mathbf{x}}}_{(1+\chi )\lambda }({\mathbf{A}}_S,\sqrt{Q}\mathbf{z}+\mathbf{y})}_0}{1+\chi },$$

$$\hat{\chi }$$

$$\begin{array}{c}({\hat{\mathbf{x}}}_{\lambda }(\mathbf{A},\mathbf{y}){)}_{i\notin S}\sim g_{\lambda }(\hat{Q},\sqrt{\hat{\chi }}{z}_i)\\ =\underset{x}{\min }\text{quantity}(\frac{\hat{Q}}{2}{x}^2-\sqrt{\hat{\chi }}{z}_{i}x+M\lambda \text{absolutevalue}x),\end{array}$$

$$\begin{array}{c}⟨\Psi (\mathbf{x})⟩:=\underset{\beta \to \infty }{\lim }{\mathbb{E}}_{\mathbf{A},\mathbf{y}}\int \text{differential}\mathbf{x}{P}_{\beta }(\mathbf{x}|\mathbf{A},\mathbf{y})\Psi (\mathbf{x})\\ =-\underset{\beta \to \infty }{\lim }\underset{h\to 0}{\lim }\text{partialderivative}h{\beta }^{-1}{\mathbb{E}}_{\mathbf{A},\mathbf{y}}\log Z_{\beta }(\mathbf{A},\mathbf{y};h\Psi ),\end{array}$$

$$\begin{array}{c}{Z}_{\beta }(\mathbf{A},\mathbf{y};h\Psi ):=\int \text{differential}\mathbf{x}{e}^{-\frac{\beta }{2}\text{norm}{\mathbf{A}\mathbf{x}-\mathbf{y}}^2-\beta M\lambda \text{norm}{\mathbf{x}}_1-\beta h\Psi (\mathbf{x})}.\end{array}$$

$$⟨\sum _{i\notin S}\psi (x_i)⟩=\tilde{N}{\mathbb{E}}_{{\mathbf{A}}_S,\mathbf{y}}\int Dz\psi (g_{\lambda }(\hat{Q},\sqrt{\hat{\chi }}z)),$$

$$⟨\Psi ({\mathbf{x}}_S)⟩={\mathbb{E}}_{\mathrm{eff}}\Psi ({\hat{\mathbf{x}}}_{(1+\chi )\lambda }({\mathbf{A}}_S,\sqrt{Q}\mathbf{z}+\mathbf{y})),$$

$$⟨\mathrm{TP}⟩$$

$$⟨\mathrm{FP}⟩$$

$$⟨{ϵ}_x⟩$$

$$s_{\lambda }^{(M)}:=\frac{1}{M}\text{norm}{{\mathbf{\gamma }}_{\lambda }(\mathbf{y})-\mathbf{y}}^2$$

$${\Gamma }^{(M)}:=\frac{1}{M}\int \text{differential}\mathbf{\xi }{p}_{\xi }(\mathbf{\xi })\text{norm}{\mathbf{\xi }}^2<C$$

$$\alpha (1+ϵ)>{\alpha }_C=\frac{{\bar{s}}_{\lambda }}{2{\lambda }^2},$$

$$\alpha (1+ϵ)>2\text{quantity}(d+\frac{\Gamma }{{\lambda }^2}),$$

$$\mathrm{Pr}\text{quantity}[{\lambda }_{\max }(\mathbf{B})\ge \text{quantity}(\sqrt{M}+\sqrt{d}+t{)}^2]\le e^{-\frac{t^2}{2}}$$

Full text

Our analysis makes it possible to assess the average performance of Lasso with Gaussian sensing matrices without assumptions on the scaling of $N$ and $M$, the noise distribution, and the profile of the true signal. Under mild conditions on the noise distribution, the analysis also offers a lower bound on the sample complexity necessary for partial and perfect support recovery when $M$ diverges as $M=O(\log N)$. The obtained bound for perfect support recovery is a generalization of that given in previous literature, which only considers the case of Gaussian noise and diverging $d$. Extensive numerical experiments strongly support our analysis.

1 Introduction

An important objective of high dimensional statistics is to extract information in situations where the signal’s dimension $N$ is overwhelmingly large compared to the accumulated sample size $M$. It is crucial to incorporate prior knowledge on the signal structure to reduce the signal space dimension for reliable estimation. A particularly common assumption is sparsity, which postulates that the true signal has few nonzero entries. Exploiting this property allows one to obtain robust and interpretable results specifying the few relevant variables explaining the retrieved data (Donoho, , 2006).

For instance, consider the sparse linear regression problem where measurements $\bm{y}\in\mathbb{R}^{M}$ of the signal $\bm{x}^{0}\in\mathbb{R}^{N}$ with $d$ non-zero components are given by the linear model

[TABLE]

where $\mathbf{A}\in\mathbb{R}^{M\times N}$ is the sensing matrix, and $\bm{\xi}\in\mathbb{R}^{M}$ is the noise vector distributed according to $p_{\xi}(\bm{\xi})$. The most fundamental yet popular sparse signal estimation method is the least absolute shrinkage and selection operator (Lasso) (Tibshirani, , 1996), which offers the estimator by solving the following convex program:

[TABLE]

where $\lambda$ is a regularization parameter. Since its introduction, this simple $\ell_{1}$-regularization scheme has been successfully adapted as a backbone technique for solving a wide variety of sparse estimation problems. A particularly interesting question to ask is if one can make any guarantees on the performance of Lasso under general scalings of $(N,M,d)$, its dependence on $\lambda$, and statistical properties of the noise and true signal.

A sheer amount of research has been devoted to assessing the performance of Lasso. Traditionally, research based on the irrepresentability condition (Meinshausen and Bühlmann, , 2006; Zhao and Yu, , 2006) has been popular in establishing guarantees in terms of support recovery of the sparse signal (Wainwright, 2009b, ; Dossal et al., , 2012; Meinshausen and Bühlmann, , 2006; Zhang and Huang, , 2008; Candès and Plan, , 2009; Zhao and Yu, , 2006). A different approach based on approximate message-passing (AMP) theory (Donoho et al., , 2009), and the heuristical replica method (Mézard et al., , 1986) from statistical physics has focused on assessing the sharp, asymptotic properties of Lasso in the large $N$ and $M$ limit under random sensing matrix designs. Despite the previous works, the understanding of the Lasso estimator is still limited. Analysis based on the irrepresentability condition often offers only scaling guarantees with respect to $(N,M,d)$, or statements with strong assumptions on the regularization parameter. Besides, the AMP/replica-based analysis has been only limited to linear sparsity, i.e. $d/N=O(1)$ and $M/N=O(1)$ as $N\to\infty$, which may be somewhat unrealistic compared to real-world situations.

1.1 Contributions

In this work, we complement the drawbacks in both the irrepresentability condition approach and AMP / replica approach by theoretically analyzing the average performance of Lasso when $d=O(1)$, i.e. the ultra-sparse case (Donoho et al., , 1992; Bhadra et al., , 2017), which is a more typical situation in certain applications such as materials informatics (Ghiringhelli et al., , 2015; Kim et al., , 2016; Pilania et al., , 2016). Moreover, our result offers a necessary condition for support recovery in the limit $N,M\to\infty$. Specifically, our contributions are summarized as follows:

• •

We provide a new way to apply the replica method in the ultra-sparsity regime. This is done by explicitly handling the correlations and finite-size effects acting on the active set ${\rm supp}(\bm{x}^{0})=\quantity{i\ |\ x_{i}^{0},\neq 0\ i=1,\cdots,N}$, which is otherwise ignored in conventional analysis (Section 2.1, Claim 1).

• •

Using this enhanced replica method, we precisely evaluate the average property of Lasso under ultra-sparsity and standard Gaussian matrix design, i.e. each element of $\mathbf{A}$ is i.i.d. according to a standard Gaussian distribution. This provides an extension to previous results derived from the AMP theory and the replica method, where linear sparsity is necessary for the analysis (Section 2.2, Claim 2).

• •

We derive a necessary condition for partial support recovery ${\rm supp}(\hat{\mathbf{x}}_{\lambda}(\mathbf{A},\mathbf{y}))\subseteq{\rm supp}(\bm{x}^{0})$ under some mild conditions (Assumption 1). Specifically, the number of false positives, and subsequently the model misselection probability vanishes only if $M>\alpha_{C}\log N$ for $N\to\infty$. This constant $\alpha_{C}$ is determined by the mean prediction error of an oracle (Section 2.3, Claim 3, 4).

• •

In addition to partial support recovery, the analysis also provides a necessary condition for perfect support recovery ${\rm supp}(\hat{\mathbf{x}}_{\lambda}(\mathbf{A},\mathbf{y}))={\rm supp}(\bm{x}^{0})$, which generalizes the sample complexity bound given by Wainwright, 2009b for i.i.d. Gaussian noise distributions in the limit $d\to\infty$ to more general noise distributions under constant $d$ (Section 2.3, Claim 5).

• •

We demonstrate that our theory agrees well with experiment by conducting extensive numerical simulations (Section 3).

Note that all of the results are derived from the enhanced replica method, which is yet to be proven rigorously; hence the statements are presented as claims.

1.2 Related Work

Irrepresentability Condition.

As aforementioned, the irrepresentability condition, first introduced by Meinshausen and Bühlmann, (2006) and Zhao and Yu, (2006), has been an important cornerstone, as it establishes a sufficient condition for perfect support recovery. This condition indicates whether the covariates, i.e. the columns of $\mathbf{A}$, are linearly independent enough to be distinguishable from one another, and hence variable selection is relatively feasible. It has been revealed that Lasso is an “optimal” support estimator in the sublinear regime $d=o(N)$, i.e. Lasso has its success/failure threshold for sample complexity in the same order as the informational-theoretical one (Fletcher et al., , 2009; Wainwright, 2009a, ). However, little is known about the constants involved in these conditions. Wainwright, 2009b provided necessary and sufficient conditions for perfect support recovery under random Gaussian matrices for diverging $d$. This is a simple and explicit bound which depends on the regularization parameter and intensity of the noise, which is restricted to i.i.d. Gaussian. Focusing on the case $d=O(1)$, Dossal et al., (2012) derived sufficient conditions for partial and perfect support recovery under deterministic noise, whose bound is similar to the one given in Wainwright, 2009b .

AMP theory.

A particular line of work has aimed in assessing the properties of Lasso under general random matrix designs via careful analysis of the dynamical behavior of the AMP algorithm (Kabashima, , 2003; Donoho et al., , 2009; Takahashi and Kabashima, , 2022), whose convergence point coincides with (2) in the large $N$ limit. Rather than establishing inequality bounds or conditions, the objective is to establish sharp results on the Lasso for a random instance of $(\mathbf{A},\bm{y})$. Although analysis is limited to linear sparsity regime, powerful and precise results have been proven rigorously under this framework (Bayati and Montanari, , 2012). For instance, Su et al., (2017) and Wang et al., (2020) determine the possible rate of false positives and true positives achievable under certain settings, which can be obtained by solving a small set of nonlinear equations. Nevertheless, the analysis does not give insight on support recovery, since this is impossible in the linear sparsity regime (Fletcher et al., , 2009; Wainwright, 2009a, ).

Replica method.

Results similar to those from AMP theory have also been derived by using the non-rigorous replica method in statistical mechanics. Unlike AMP theory, which is based on a convergence analysis of a particular algorithm, the replica method aims at directly calculating the average over $(\mathbf{A},\bm{y})$ of a cumulant generating function for some probability distribution, i.e. of the form $K_{\phi}(t)=\mathbb{E}_{\mathbf{A},\bm{y}}\ \log\int\differential\bm{x}\ e^{t\phi(\bm{x})}p(\bm{x}|\mathbf{A},\bm{y})$. This calculation is often encountered in the field of statistics, where one is interested in the average behavior of a statistical model. While lacking a complete proof, this method has been successful in predicting the average performance of machine learning and optimization methods under general random designs in the linear sparsity regime (Vehkaperä et al., , 2016; Zdeborovà and Krzakala, , 2016). In fact, under certain assumptions, the average predictions given by the replica method have been proven to be consistent with the asymptotic results obtained from AMP theory and other rigorous methods (Stojnic, , 2013; Thrampoulidis et al., , 2018). Similar to AMP theory, however, reliable adaptations of this method outside linear sparsity are still open problems. Previous research such as Abbara et al., (2020), Meng et al., 2021a and Meng et al., 2021b analyzed the performance of sparse Ising model selection using a variation of the replica method. However, this was accomplished through a series of ansatzes which are generally difficult to justify theoretically.

1.3 Preliminaries

Here we summarize the notations used in this paper. The expression $\norm{\cdot}$ denotes the $\ell_{2}$ norm. The active set $S$ is defined as the support of the $d$-sparse true signal $\bm{x}^{0}$, $S:={\rm supp}(\bm{x}^{0})=\{i\ |\ x^{0}_{i}\neq 0,\ i=1,\cdots,N\}$. Define $\tilde{N}:=N-d$, the size of the inactive set. The matrix $\mathbf{A}_{S}$ denotes the submatrix constructed by concatenating the columns of $\mathbf{A}$ with indices in $S$. The vector $\bm{x}^{0}_{S}$ denotes the subvector of $\bm{x}^{0}$ with indices in $S$. For simplicity, $\bm{x}^{0}$ is assumed to be a deterministic, although this can be extended to random signals trivially. The expression $\mathbb{E}_{\mathbf{A},\bm{y}}$ denotes the average over the joint probability with respect to the pair $(\mathbf{A},\bm{y})$, i.e.

[TABLE]

where $\delta(\cdot)$ denotes the Dirac delta function. The definition of $\mathbb{E}_{\mathbf{A}_{S},\bm{y}}$ follows straightforwardly from the above. Also, define $D\bm{z}$ as the standard Gaussian measure, $D\bm{z}=\differential\bm{z}e^{-\norm{\bm{z}}^{2}/2}/(2\pi)^{n/2}$ for $\bm{z}\in\mathbb{R}^{n}$. Given $(\mathbf{A},\bm{x}^{0},\bm{y})$ and regularization parameter $\lambda$, the oracle Lasso estimator is defined as $\hat{\bm{x}}_{\lambda}(\mathbf{A}_{S},\bm{y})$, which is the Lasso estimator with the true support identified beforehand. It is also convenient to define the oracle Lasso fit, defined by $\bm{\gamma}_{\lambda}(\bm{y}):=\mathbf{A}_{S}\hat{\bm{x}}_{\lambda}(\mathbf{A}_{S},\bm{y})$, with its dependence on $\mathbf{A}_{S}$ suppressed for convenience.

Given configuration $(\mathbf{A},\bm{x}^{0},\bm{y})$, and regularization parameter $\lambda$, the number of false positives $\rm FP$ and the number of true positives $\rm TP$ of the lasso estimator is defined as

[TABLE]

where $S^{C}$ denotes the complement of set $S$ from $\quantity{1,\cdots,N}$. Without confusion, the dependence on $(\mathbf{A}_{S},\bm{y})$ is suppressed for convenience.

We say that an event $A$ holds with asymptotically high probability (w.a.h.p.) if there exists a constant $c>0$ such that ${\Pr}[A]>1-O(N^{-c})$. We also say that $A$ holds with probability approaching one (w.p.a.1) if ${\Pr}[A]>1-o(1)$ as $N\to\infty$.

2 Replica analysis

Define the Boltzmann distribution as

[TABLE]

where $Z_{\beta}(\mathbf{A},\bm{y})$ is the normalization constant. Note that in the limit $\beta\to\infty$, (5) converges to a point-wise distribution concentrated on the Lasso estimator $\hat{\bm{x}}_{\lambda}(\mathbf{A},\bm{y})$. The main objective of our analysis is to calculate the average of the logarithm of $Z_{\beta}\quantity(\mathbf{A},\bm{y})$ over the random variables $(\mathbf{A},\bm{y})$ in the limit $\beta\to\infty$, which is called the free energy or the cumulant generating function

[TABLE]

The properties of $\hat{\bm{x}}_{\lambda}(\mathbf{A},\bm{y})$ averaged over the population of $(\mathbf{A},\bm{y})$ can then be assessed by taking appropriate derivatives of $\mathcal{F}$.

Although (6) is difficult to calculate straightforwardly, this can be resolved by using the replica method (Mézard and Montanari, , 2009; Mézard et al., , 1986), which is based on the following equality

[TABLE]

Instead of handling the cumbersome $\log$ expression in (6) directly, one calculates the average of the $n$-th power of $Z_{\beta}$ for $n\in\mathbb{N}$, analytically continues this expression to $n\in\mathbb{R}$, and finally takes the limit $n\to+0$. Based on this replica ”trick”, it suffices to calculate

[TABLE]

up to the first order of $n$ to take the $n\to+0$ limit in the right hand side of (7).

2.1 Outline of the derivation

Here, we only give a brief outline of the derivation; for details, see Supplementary Materials. Rewriting $\Delta^{a}_{i}:=x^{a}_{i}\ (i\notin S)$, it is convenient to introduce the auxillary variable $h_{\mu}^{a}\equiv\sum_{i\notin S}A_{\mu i}\Delta_{i}^{a}\ (\mu=1,\cdots,M)$, which accounts for the effect from the variables not in the true support in each replica $a$. A crucial observation is that $\quantity{A_{\mu i}}_{1\leq\mu\leq M,\\ i\notin S}$ is statistically independent from $(\mathbf{A}_{S}$, $\bm{y})$, which allows the average to be taken individually.

By taking the average over the Gaussian variables $\quantity{A_{\mu i}}_{1\leq\mu\leq M,\\ i\notin S}$ first, we find that $h_{\mu}^{a}$ is Gaussian with zero mean and covariance $\mathbb{E}h_{\mu}^{a}h_{\nu}^{b}=\delta_{\mu\nu}\sum_{i\notin S}\Delta_{i}^{a}\Delta_{i}^{b}$. By assuming the replica symmetric (RS) ansatz (Mézard et al., , 1986)

[TABLE]

the integral for the replicated vectors $\quantity{\bm{\Delta}^{a}}_{a=1}^{n}$ over the whole $\mathbb{R}^{\tilde{N}\times n}$ space is restricted to a subspace satisfying the constraints (9). More explicitly, one can rewrite (8) as

[TABLE]

where $\mathcal{I}$ corresponds to the contribution from the RS constraint: i.e.

[TABLE]

and $\mathcal{L}$ is the contribution from the second line of (8), albeit simplified as a result of replica symmetry:

[TABLE]

By using the Fourier representation of the delta function, (11) can be further rewritten as

[TABLE]

Using this expression, the integral with respect to $\{\bm{\Delta}^{a}\}_{1\leq a\leq n}$ in (10) can be calculated analytically. Performing the saddle point approximation for large $M$ to the integrals with respect to $(Q,\hat{Q},\chi,\hat{\chi})$, and finally taking the limit $\beta\to\infty$ after $n\to+0$ in (7) yields the following expression for $\mathcal{F}$.

Claim 1**.**

The free energy is given by

[TABLE]

Here, $\Lambda:=(M\lambda)^{2}$, ${\rm erfc}$ is the complementary error function ${\rm erfc}(x):=2/\sqrt{\pi}\int_{x}^{\infty}\differential te^{-t^{2}}$, and $\operatorname*{{\rm Extr}}$ refers to the extremum condition with respect to $\Theta:=(Q,\hat{Q},\chi,\hat{\chi})$, which are random variables dependent on $(\mathbf{A}_{S},\bm{y})$.

Straightforward calculation shows that the extremum conditions are given by

[TABLE]

where the second equality in (2.1) is from Theorem 1 in Tibshirani and Taylor, (2012). Note that the dependence of $\Theta$ on $(\mathbf{A}_{S},\bm{y})$ is not explicitly written for sake of simplicity. This evaluation of $\mathcal{F}$ reduces the high-dimensional integral over $\mathbf{A}$ and $\bm{y}$ to an average over a four-dimensional extremum problem involving a $M$–dimensional integral with respect to $\bm{z}$, which can be numerically computed via iterative substitution and Monte Carlo sampling over $(\bm{A}_{S},\bm{y})$ and $\bm{z}$.

It is interesting to compare our replica analysis in the large $N$ and $M$ limit to the ones considering linear sparsity (Kabashima et al., , 2009; Vehkaperä et al., , 2016). In linear sparsity, the lasso estimator’s statistical property can effectively be described by a population of $N$ decoupled, independent scalar estimators under Gaussian noise with identical intensity as $N\to\infty$. This is often referred to as the decoupling principle in information theory; see Guo and Verdú, (2005) and Bayati and Montanari, (2011) for details. In the ultra-sparse case, the elements of the Lasso estimator in the active set, consisting of $d=O(1)$ terms, cannot be expected to decouple, as finite-size effects of non-Gaussian and correlated nature are expected to be significant to describe its profile. This is why a $d-$body optimization procedure and the average with respect to $(\mathbf{A}_{S},\mathbf{y})$ appears explicitly in (14). On the other hand, the decoupling principle is implicitly employed for the $\tilde{N}$ non-active variables conditioned on $(\mathbf{A}_{S},\bm{y})$. More explicitly, for each configuration of $(\mathbf{A}_{S},\bm{y})$, each element of the non-active Lasso estimator is statistically equivalent to

[TABLE]

where $z_{i}$ are i.i.d. according to $\mathcal{N}(0,1)$. Note that the decoupling principle, rigorously proven under AMP theory, does not necessarily need $N$ and $M$ to diverge at the same rate (Rush and Venkataramanan, , 2018).

2.2 Performance assessment of Lasso

The free energy allows convenient evaluation of averages of certain functions of the estimator. More explicitly, for a function $\Psi:\mathbb{R}^{N}\to\mathbb{R}$, its average with respect to the Boltzmann distribution (5) and $(\mathbf{A},\bm{y})$ is given by

[TABLE]

where

[TABLE]

For a class of functions $\Psi$, the above can be calculated trivially, which we state in the following claim:

Claim 2** (Average with respect to active and inactive sets).**

For arbitrary functions $\psi:\mathbb{R}\to\mathbb{R}$ and $\Psi:\mathbb{R}^{d}\to\mathbb{R}$, we have

[TABLE]

and

[TABLE]

where $\mathbb{E}_{\rm eff}:=\mathbb{E}_{\mathbf{A}_{S},\bm{y}}\int D\bm{z}$, and $(Q,\hat{Q},\hat{\chi},\chi)$ is given by the solution of the extremum conditions (15)–(18) for each $(\mathbf{A}_{S},\bm{y})$. In particular, performance measures such as the average of true positives (${\mathrm{TP}}$), false positives (${\mathrm{FP}}$) and $\ell_{2}$ error $\epsilon_{x}:=\norm{\bm{x}_{\lambda}(\mathbf{A},\bm{y})-\bm{x}^{0}}^{2}$ is given by

[TABLE]

2.3 Necessary condition for support recovery

A particular topic of interest is partial support recovery, and the minimum number of samples $M$ necessary for the false positives to vanish in the limit $N\to\infty$. Although the fixed point equations (15) –(18) do not admit a closed form solution, a necessary condition in terms of the sample complexity can be derived under the following mild conditions:

Assumption 1**.**

**

1. A:

(Uniqueness of fixed point) The solutions of the fixed point equations (15)–(18) are unique and satisfy $(Q,\hat{Q},\chi,\hat{\chi})\in(0,\infty)^{4}$. 2. B:

(Concentration of the oracle Lasso estimator) The random variable

[TABLE]

*has finite mean $\bar{s}_{\lambda}^{(M)}$ and variance converging to zero. *** 3. C:

(Bounded variance of noise distribution) The distribution $p_{\xi}$ satisfies

[TABLE]

for some constant $C$.

Claim 3** (Necessary sample complexity for asymptotically zero false positives).**

*Let $M$ diverge with $N$ with scaling $M=\alpha\log N\ (\alpha>0)$. Under Claim 1 and Assumption 1, if there exists a constant $c>0$ such that $\left\langle{{\mathrm{FP}}}\right\rangle<O(N^{-c})$ in the limit $N\to\infty$, then *

[TABLE]

holds for any constant $\epsilon>0$, where $\displaystyle\bar{s}_{\lambda}=\lim_{M\to\infty}\bar{s}_{\lambda}^{(M)}.$

The proof is postponed to Section 4. From this claim, the necessary sample complexity for partial support recovery follows immediately:

Claim 4** (Necessary sample complexity for partial support recovery).**

Under the settings in Claim 3, if ${\rm supp}(\hat{\bm{x}}_{\lambda}(\mathbf{A},\bm{y}))\subseteq{\rm supp}(\bm{x}^{0})$ w.a.h.p., then $\alpha>\alpha_{C}$.

By definition, $\bar{s}_{\lambda}$ is the prediction error of the oracle, which is given the sensing submatrix $\mathbf{A}_{S}$ and observation vector $\bm{y}$. This is reminiscent of the primal-dual witness construction in Wainwright, 2009b , where sufficient conditions for asymptotically zero FPs are derived by solving the oracle Lasso first, and observing whether the oracle solution concatenated with $N-d$ zero elements is a unique solution of the original Lasso problem (2).

Furthermore, the necessary condition for perfect support recovery can also be derived using Claim 3.

Claim 5** (Necessary sample complexity for perfect support recovery).**

Under the settings in Claim 3, suppose ${\rm supp}(\hat{\bm{x}}_{\lambda}(\mathbf{A},\bm{y}))={\rm supp}(\bm{x}^{0})$ holds w.a.h.p. Then

[TABLE]

holds for any constant $\epsilon>0$, where $\displaystyle\Gamma=\lim_{M\to\infty}\Gamma^{(M)}$.

Note that in the special case of Gaussian noise with variance $\sigma^{2}$, we have $\Gamma=\sigma^{2}$, which extends the result of Wainwright, 2009b , Theorem 4 to the case $d=O(1)$. Moreover, our result can be applied to any noise distribution satisfying Assumption 1.C.

3 Numerical experiments

3.1 Non-asymptotic results

To verify the derived results based on Claim 1, numerical experiments were conducted. For simplicity, we consider the case where the active set has size $d=3$ with $\bm{x}^{0}_{S}=\bm{1}_{3}$, and $\bm{\xi}$ is generated from a Gaussian distribution with variance $\sigma^{2}$. Here, the value of $d$ is taken to be small enough such that finite-size effects are nonignorable. The values of $\left\langle{{\mathrm{TP}}}\right\rangle,\left\langle{{\mathrm{FP}}}\right\rangle$ and $\epsilon_{x}$ obtained from our replica predictions (24)-(26) are compared with the average over $10^{4}$ experimental runs. The average with respect to $(\bm{A}_{S},\bm{y})$ for obtaining the replica prediction was approximated using a Monte Carlo procedure over $10^{6}$ samples.

Figure 1 shows that all three values from theory and experiment are in good agreement for parameters $(\lambda,\sigma^{2})=(0.5,0.0)$ and $(0.5,0.5)$.

3.2 Asymptotic results

Claims 3 and 4 are also verified via numerical experiments; see Supplementary Materials for numerical experiments on Claim 5. In order to access the critical point $\alpha_{C}$ in (27), Monte Carlo experiments were conducted to evaluate $s_{\lambda}^{(M)}$ for different values of $M$. Figure 2 shows the value of $s_{\lambda}^{(M)}$ at $(\lambda,\sigma^{2})=(0.5,0.0)$ and $(0.5,0.5)$ for both $\bm{x}_{S}^{0}=\bm{1}_{3}$ and $\bm{x}_{S}^{0}=[\frac{1}{3},\frac{2}{3},1]$. From its asympototic behavior, $\alpha_{C}$ can be evaluated as the values given in Table 1. Interestingly, for the case $\bm{x}_{S}^{0}=\bm{1}_{3}$, $s_{\lambda}$ approaches $6$ and $10$ for $\sigma^{2}=0$ and $0.5$ respectively, which is equivalent to $2(d+\Gamma/\lambda^{2})$ given in Claim 5.

Figure 3 shows the average number of FP and partial support recovery probability over 10,000 experimental runs for $\alpha$ in the vicinity of the numerically evaluated $\alpha_{C}$ for different values of $N$. We observe that for $\alpha<\alpha_{C}$, the average FP is consistently nondecreasing with respect to $N$, while partial support recovery probability is consistently nonincreasing with respect to $N$.

4 Proofs

4.1 Proof of Claim 3

The following lemmas will be useful in the proof.

Lemma 1** (Lemma 1, Dossal et al., (2012)).**

There is a finite increasing sequence $(\lambda_{t})_{t\leq K}$ with $\lambda_{0}=0$ such that for all $t<K$, the sign and support of $\hat{\bm{x}}_{\lambda}(\mathbf{A}_{S},\bm{y})$ are constant on each interval $(\lambda_{t},\lambda_{t+1})$.

Lemma 2** (Lemma 1, Tibshirani and Taylor, (2012)).**

The Lasso fit is 1-Lipschitz continuous with respect to $\ell_{2}$ norm.

Lemma 3** (Theorem II.13, Davidson and Szarek, (2001)).**

Let $\mathbf{A}\in\mathbb{R}^{M\times d}$ be a random matrix with i.i.d standard Gaussian entries. The largest and smallest eigenvalue of $\mathbf{B}=\mathbf{A}^{\mathsf{T}}\mathbf{A}$ satisfy

[TABLE]

for $t>0$ and

[TABLE]

for $0<t<\sqrt{M}-\sqrt{d}$.

We now prove Claim 3. Define

[TABLE]

Let us evaluate the difference between $s^{(M)}_{(1+\chi)\lambda,Q}$ and $s^{(M)}_{\lambda,0}=s^{(M)}_{\lambda}$ when $\left\langle{{\mathrm{FP}}}\right\rangle<O(N^{-c})$. Using the Cauchy Schwartz inequality and symmetry $\bm{\gamma}_{{\lambda}}({\bm{y}})=-\bm{\gamma}_{{\lambda}}({-\bm{y}})$,

[TABLE]

The triangle inequality and Lemma 2 implies that

[TABLE]

and similarily,

[TABLE]

To derive a bound for the last term in (4.1) and (4.1), Lemma 1 is employed. Let the support and sign of $\hat{\bm{x}}_{\lambda}(\mathbf{A}_{S},\bm{y})$ be constant in intervals $(M\lambda_{t},M\lambda_{t+1})\ (t=0,\cdots,K-1)$, where $\lambda=\lambda_{0}<\cdots<\lambda_{K}=(1+\chi)\lambda$. Let the support set in interval $(\lambda_{t},\lambda_{t+1})$ be given by $I_{t}$, and define $\bm{s}_{t}\in\quantity{-1,0,1}^{\absolutevalue{I_{t}}}$ be the sign vector of $\hat{\bm{x}}_{\lambda^{\prime}}(\mathbf{A}_{S},\bm{y})\ (\lambda^{\prime}\in(\lambda_{t},\lambda_{t+1}))$ restricted to $I_{t}$. From the KKT conditions, the Lasso fit is expressed as

[TABLE]

where $\mathbf{M}^{+}$ denotes the pseudoinverse of matrix $\mathbf{M}$. We deduce

[TABLE]

Lemma 3, with the inclusion principle $\rho((\mathbf{A}_{SI_{t}}^{\mathsf{T}}\mathbf{A}_{SI_{t}})^{-1})\leq\rho((\mathbf{A}_{S}^{\mathsf{T}}\mathbf{A}_{S})^{-1})$ implies that w.a.h.p., $\norm{\bm{\gamma}_{{(1+\chi)\lambda}}({\bm{y}})-\bm{\gamma}_{{\lambda}}({\bm{y}})}\leq 2\chi\sqrt{d\lambda^{2}M}$. The relations (31) – (4.1), and inequality $\int D\bm{z}\norm{\bm{z}}=\sqrt{2}\Gamma((M+1)/2)/\Gamma(M/2)<\sqrt{M}$ then leads to the following holding w.a.h.p.

[TABLE]

We now use the following lemma which shows that $Q$ and $\chi$ are negligible almost surely.

Lemma 4**.**

Under the assumptions of Claim 3, $\chi<N^{-c/2}$ and $Q<N^{-c/4}$ holds w.a.h.p.

The proof is given in Supplementary Materials. Since $\norm{\bm{y}}<\norm{\mathbf{A}_{S}\bm{x}^{0}}+\norm{\bm{\xi}}$ is bounded by $M^{2}$ w.p.a.1 from Lemma 3 and Assumption 1.C, the right hand side of eq. (4.1) is of $O(N^{-c/8})$ w.p.a.1.

We therefore have

[TABLE]

On the other hand, the extremum conditions (16) and (2.1) imply that $\hat{\chi}$ is always bounded.

Lemma 5**.**

Suppose the extremum conditions (15)–(18) are satisfied. Then, the variable $\hat{\chi}$ satsfies

[TABLE]

Combined with (35), for sufficiently large $M$

[TABLE]

holds for arbitrary constant $\epsilon>0$. This implies that $\frac{1}{2}\alpha\lambda^{2}(1+\epsilon)$ must be larger than the median of $s_{\lambda}^{(M)}$. Now, the difference between the median and average is no larger than one standard deviation, which is negligible from Assumption 1.B. This yields the statement of the claim in the limit $M\to\infty$.

4.2 Proof of Claim 4

From Theorem 6 in Osborne et al., (2000), the number of false positives is bounded by $\min(M,N)$. Hence, we have $\left\langle{{\mathrm{FP}}}\right\rangle<M\times{\Pr}\quantity[{\mathrm{FP}}\neq 0]=O(N^{-c})$ for some $c>0$. The statement of Claim 4 then follows from Claim 3.

4.3 Proof of Claim 5

From Claim 3 and 4, it suffices to show that

[TABLE]

The KKT conditions imply that w.p.a.1, $\mathbf{A}_{S}\hat{\bm{x}}=\mathbf{A}_{S}\mathbf{A}_{S}^{+}(\bm{y}-M\lambda(\mathbf{A}_{S}^{+})^{\mathsf{T}}\bm{s}),$ where we abbreviated $\hat{\bm{x}}:=\hat{\bm{x}}_{\lambda}(\mathbf{A}_{S},\bm{y}),$ and $\bm{s}={\rm sgn}\quantity(\hat{\bm{x}})$. Therefore, $\bm{y}-\mathbf{A}_{S}\hat{\bm{x}}$ can be decomposed into a sum of two linearly independent vectors

[TABLE]

where $\bm{v}:=M\lambda\mathbf{A}_{S}(\mathbf{A}_{S}^{\mathsf{T}}\mathbf{A}_{S})^{-1}\bm{s},$ $\bm{v}_{\perp}:=\mathcal{P}_{\ker(\mathbf{A}_{S})}(\bm{y})=\mathcal{P}_{\ker(\mathbf{A}_{S})}(\bm{\xi})$, and $\mathcal{P}_{\ker(\mathbf{A}_{S})}$ is the projection onto the kernel of $\mathbf{A}_{S}$. The average of the squared norm of $\bm{v}$ can be evaluated as

[TABLE]

where the last inequality follows from Jensen’s inequality and $\mathbb{E}_{\rm eff}\lambda_{\min}(\bm{A}_{S}^{\mathsf{T}}\bm{A}_{S})\geq(\sqrt{M}-\sqrt{d})^{2}$ (Davidson and Szarek, , 2001).

To obtain a lower bound on the squared norm of $\bm{v}_{\perp}$, fix the vector $\bm{\xi}$. Noticing that entries of $\mathbf{A}_{S}^{\mathsf{T}}\bm{\xi}/\norm{\bm{\xi}}^{2}$ are i.i.d. standard Gaussian, the tail bound for $\chi^{2}$–random variables (Laurent and Massart, , 2000) implies that for some constant $C>0$,

[TABLE]

Using this inequality, (30) with $t=\sqrt{2\log M}$ and the union bound, we have that

[TABLE]

Equation (38) immediately follows from (40) and (4.3), which completes the proof.

5 Conclusion

In this paper, we provided an analysis based on an enhanced replica method for assessing the average performance of the Lasso estimator under ultra-sparse conditions. Besides, we deduced conditions necessary for support recovery which are derived from the oracle Lasso estimator. Numerical experiments strongly support the validity of our analysis.

The methodological novelty originates from an observation of finite-size effects and correlations within the active set, which is implicitly assumed to be negligible in the conventional replica analysis. We anticipate that this framework is applicable to analysis of other machine learning or optimization problems where finite-size effects are nonnegligible. Extending this method further to more general sensing matrix ensembles is also another exciting direction for future work.

Acknowledgements

This work was partially supported by JSPS KAKENHI Grant Nos. 22J21581 (KO), 21K21310 (TT), 17H00764, 19H01812, 22H05117 (YK) and JST CREST Grant Number JPMJCR1912 (YK).