A statistical mechanics approach to de-biasing and uncertainty estimation in LASSO for random measurements
Takashi Takahashi Yoshiyuki Kabashima


A statistical mechanics approach to de-biasing and uncertainty estimation in LASSO for random measurements]
A statistical mechanics approach to de-biasing and uncertainty estimation in LASSO for random measurements


takahashi.t.cc@m.titech.ac.jp
Department of Mathematical and Computing Science 

Tokyo Institute of Technology

2-12-1, Ookayama, Meguro-ku, Tokyo, Japan


In high-dimensional statistical inference in which the number of parameters to be estimated is larger than that of the holding data, regularized linear estimation techniques are widely used. 
These techniques have, however, some drawbacks. 
First, estimators are biased in the sense that their absolute values are shrunk toward zero because of the regularization effect. 
Second, their statistical properties are difficult to characterize as they are given as numerical solutions to certain optimization problems.
In this manuscript, we tackle such problems concerning LASSO, which is a widely used method for sparse linear estimation, when the measurement matrix is regarded as a sample from a rotationally invariant ensemble.
We develop a new computationally feasible scheme to construct a de-biased estimator with a confidence interval and conduct hypothesis testing for the null hypothesis that a certain parameter vanishes.
It is numerically confirmed that the proposed method successfully de-biases the LASSO estimator and constructs confidence intervals and p-values by experiments for noisy linear measurements.



[
    Takashi Takahashi and Yoshiyuki Kabashima
    
=============================================








§ INTRODUCTION
  

Estimating high-dimensional unknown variables from a limited number of data precisely and reliably is an important task in statistics, machine learning, signal processing, and so on.
    For instance, such demands arise in compressed sensing <cit.> and genomics <cit.>.  
	 Since, in these problems, the number of parameters often far surpasses that of observed data, it is 0 intuitively clear that some sparsity assumptions on the parameters are necessary to reasonably estimate them.
     Therefore, one needs to simultaneously solve two problems: variable selection, which seeks relevant (or non-zero) parameters for the data generation process, and parameter estimation.
     In the past few decades, a number of methods have been developed to tackle such problems.
     One of the most successful approaches is the least absolute shrinkage and selection operator (LASSO) <cit.> method for high-dimensional linear regression problems in which the estimator is obtained by minimizing the L_1 norm regularized likelihood function.
     As LASSO estimators can be easily obtained by versatile algorithms for convex optimization <cit.> and have appealing consistency properties <cit.>, they have received considerable attention.
     

			Specifically, let us consider the linear measurement model:
            
    y_i = a_i^⊤x_0 + ξ_i,  ξ_i∼_i.i.d𝒩(0, σ^2),   i = 1,2,...,M,

         where x_0∈ℝ^N and a_i∈ℝ^N are the parameter (signal) and measurement vectors, respectively, σ^2 ∈ℝ is a parameter that describes the strength of the measurement noise, and 𝒩(μ, σ^2) is the normal distribution with mean μ and variance σ^2.
         Notation ⊤ means the operation of matrix/vector transpose. 
         In matrix notation, this model is expressed as 
         
    y = Ax_0 + ξ,

         where a_i^⊤ corresponds to the i'-th row of the matrix A∈ℝ^M× N.
         A is called the observation or measurement matrix by cases.
         The objective of high-dimensional linear regression is to find the parameter vector x_0, where the number of measurements M is smaller than that of the parameter N.
         Note that in this high-dimensional setting, one cannot obtain a true solution with simple linear algebra because A^⊤ A is not invertible; by contrast, in the classical setting where M > N, the unique unbiased estimator is easily obtained as x_classical = (A^⊤ A)^-1A^⊤y by using the least squares method.
     		To achieve this aim, LASSO seeks an estimator by solving an optimization problem that imposes sparsity via an L_1 penalty:　
            
    x^LASSO(y, A;λ) ≡min_x[1/2y - Ax_2^2 + λx_1],

            where λ is a hyperparameter that controls the strength of the regularization.
            This convex optimization problem can be solved efficiently by using various versatile algorithms. 
            Although LASSO might be seen as simple heuristics, it has an appealing consistency property: in a certain sparsity condition on the true parameter x_0 and an appropriate control of the regularization strength λ, the LASSO solution and x_0 are consistent in the sense that x^ LASSO - x_0 _2^2/N vanishes as the measurement ratio γ≡ M/N tends to infinity.
            For a more comprehensive review of LASSO in the context of high-dimensional settings, see <cit.>.

     

			Unfortunately, LASSO also has some drawbacks.
        First, the LASSO solution is biased as long as λ>0 is finite.
        The amplitude of the LASSO estimator x^ LASSO is shrunk toward zero by the regularization term and its absolute value is typically smaller than that of the true parameter x_0 even in an ideal sparsity assumption. 
			Second, no explicit form of the distribution is available for the estimator, as it is just expressed as a numerical solution of (<ref>). 
        Consequently, one can neither construct confidence intervals nor perform hypothesis testing for the null hypothesis that a certain element of the parameter vanishes. These bottlenecks are considered to be problematic in real applications in which the statistical reliability of the estimation result should be assessed.
        This situation is different from the one of classical statistics in which one can analytically obtain an unbiased estimator and its distribution. 
            

        To resolve the problems stated above, in this study, we develop a new scheme for de-biasing and uncertainty estimation in the LASSO estimation in the case that the observation matrix A is generated from rotationally invariant random matrix ensembles, which are concretely defined in the next section. 
        The uncertainty addressed in this study concerns the randomness that arises from the random observation matrix A and measurement noise ξ.
        Our approach is based on a careful observation of the replica analysis of LASSO and an advanced mean-field method known as expectation consistent approximation or the adaptive Thouless–Anderson–Palmer (TAP) approach <cit.> developed in machine learning <cit.> and statistical mechanics. 
        We numerically show that the proposed algorithm effectively de-biases the LASSO estimator and estimates its uncertainty.
                
  The rest of this manuscript is organized as follows.
  In section <ref>, we explain the problem setting.
  In section <ref>, we describe the result of the replica analysis of LASSO and its physical implications. 
  Then, the design of our scheme is introduced.
  The derivation of the free energy density is in <ref>.
  In section <ref>, the proposed scheme is numerically tested by experiments for noisy linear measurements using various matrix ensembles.
  The last section provides a summary.
    



§ PROBLEM SETTING
 

	

 §.§ Model specification

   In this study, we focus on random design models of (2), in which A is a random matrix and the true parameter vector x_0 is sparse in the sense that the number of its non-zero components is limited to ϱ N (0≤ϱ < 1). 
   More precisely for A, we assume that for eigenvalue decomposition A^⊤ A = O D O^⊤, O can be regarded as a random sample from the uniform distribution of the N× N orthogonal matrices and the empirical eigenvalue distribution ∑_i=1^N δ(λ - λ_i) /N, where {λ_i}_i are the eigenvalues of A^⊤ A, converges to a certain distribution ρ(λ) in the limit N→∞ with probability one.
    
    

 §.§ De-biasing and uncertainty estimation in LASSO

    Let x^LASSO(y, A;λ) be the LASSO estimator for the given y, A, and λ.
	We are interested in the two problems associated with x^LASSO(y, A; λ).
    The first problem is that the LASSO estimator is biased. 
    In other words, |𝔼[x^LASSO_i]_A, ξ-x_0,i|, (i=1,2,...,M) remains finite for λ > 0 because of the shrinkage effect caused by the regularization term λx_1.
    The second is that the LASSO estimator does not have an explicit form of the distribution. As a consequence, one can neither construct a confidence interval nor compute a p-value to conduct hypothesis testing for the null hypothesis that a certain parameter vanishes. 
    
	 In response to the aforementioned problems, we construct the following quantities.
    The first quantity is the de-biased estimators {x_i^ debiased}_i that have confidence intervals {ℐ_i(α_ sig)≡[x_i^ debiased - L_i(α_ sig), x_i^ debiased + U_i(α_ sig)]}_i with significance α_ sig.  
    The term de-biased means that this estimator coincides with the true parameter on average: 𝔼[x_i^ debiased]_A,ξ = x_0,i. 
    The second quantity is the p-values to test whether the LASSO estimator is zero or not.
    We are interested in hypothesis testing with the null hypothesis H_0,i: x_0,i = 0.
    The confidence intervals concerning the de-biased estimators and hypothesis testing via p-values assess the uncertainty in LASSO.
	
	In the past few years, several researchers have been working on the issue closely related to that stated here <cit.>. 
   These studies discuss de-biasing and hypothesis testing in high-dimensional statistics for a fixed observation matrix where the randomness comes from the measurement noise, under tight sparsity assumptions on a true sparse signal, which corresponds to the ϱ→ 0 limit in the current setting.
   In contrast to these studies, we concentrate on the case that the randomness comes from both the random observation matrix and the measurement noise without an explicit sparsity assumption on the true parameter keeping ϱ∼ O(1).
   




§ A STATISTICAL MECHANICS APPROACH
 

	
    

 §.§ Replica analysis for general rotationally invariant random design matrices and its physical implications

    To investigate how the LASSO solution depends on the true solution, observation matrix, and measurement noise, we first evaluate the free energy density corresponding to the LASSO Hamiltonian H(x) ≡y - Ax_2^2/2 + λx_1 at a zero-temperature limit:
    
    f(λ) ≡ -lim_β→∞lim_N→∞1/Nβ𝔼[ln Z(y,A;λ)]_A,ξ,

     where β is the inverse temperature and Z is the partition function:
     
    Z(A, y;λ) = ∫exp(-β/2y - Ax_2^2 - βλx_1)dx.

    We take the limit N→∞ with γ=M/N∼ O(1) fixed.
    In the zero-temperature limit β→∞, the Boltzmann distribution e^-β H(x)/Z is dominated by the configurations of the LASSO solution (<ref>).
    Hence, one can evaluate how the LASSO estimator depends on x_0, A, ξ by analyzing the macroscopic behavior of the typical free energy density (<ref>) using statistical mechanics.
    
   Since the Hamiltonian defined above has a mean-field nature in the sense that all the variables are weakly connected, the free energy density (<ref>) can be evaluated by using the replica method:
   
    f =  extr_χ,χ,Q,Q, m, m[
         G'(-χ;J)(Q-2m + ϱ -χσ^2) + γ/2σ^2
        - QQ/2  + χχ/2 + mm  .
    .+ lim_N→∞1/N∑_i=1^N ∫min_x_i{-Q/2x_i^2 + (mx_0,i + √(χ)z_i)x_i - λ|x_i| }  Dz_i 
        ],

   where extr_χ,χ,Q,Q, m, mℱ (χ,χ,Q,Q, m, m) denotes the extremization of the function ℱ with respect to its arguments and G^'(x;J) is the derivative of G(x;J) with respect to x.
   We have defined ∫ (...) Dz, J, G(x) as follows:
   
    ∫ (...) Dz≡∫ (...) e^-z^2/2/√(2π)dz, 
    
           J≡ A^⊤ A, 
    
           G(x;J) ≡ extr_z [-∫ρ_J (s)ln|z - s|ds  + zx/2] - 1/2ln x - 1/2,

   where ρ_J(s) is the asymptotic eigenvalue distribution of J.
   The derivative of the function G(x;J) has the following form:
   
    G^'(x;J) = 1/2(z(x) - 1/x),

   where z(x) is implicitly determined by the extremal condition of (<ref>):
   
    x = 𝒮_J(z(x)) ≡∫ρ_J(λ)/z(x)-λdλ.

   The transformation 𝒮_J that appears in (<ref>) is called the Stieltjes transformation of ρ_J.
   The introduced function G is connected to the R-transform ℛ_J(·) of the asymptotic eigenvalue distribution of J in studies of free probability theory <cit.>: G(x;J)=∫_0^xℛ_J(t) dt.
   <ref> provides a brief derivation of the free energy density (<ref>).
   
   The connection between the free energy density (<ref>) and macroscopic observables is as follows. 
   At the extremum, Q, m, and χ correspond to the macroscopic physical observables: Q=⟨ |x|^2⟩_A,ξ/N, m=⟨x_0^⊤x⟩_A,ξ/N, and χ = β⟨|x|^2⟩ - |⟨x⟩|^2_A,ξ/N. Each of these corresponds to the self-overlap, the overlap between the LASSO solutions and true solutions, and the macroscopic susceptibility. The notation ⟨ ... ⟩ represents the Boltzmann average in the zero-temperature limit: ⟨...⟩≡lim_β→∞∫(...)e^-β H(x)dx/Z. In addition, from direct calculation, one can show the following relationships between the free energy density, regularization term, and residual sum of squares:
   
    f    = γ/2RSS + r̅, 
    r̅   ≡⟨1/N∑_i=1^N |x_i | ⟩_A,ξ = 
         χχ + mm - QQ, 
    RSS   =  RSS_A,ξ≡⟨1/My - Ax_2^2⟩_A,ξ,
       = 2/γ[G'(-χ; J ) (Q- 2m + ϱ - χσ^2) + γ/2σ^2 - 1/2χχ],

   where r̅ and RSS represent the per-element average of the regularization term and residual sum of squares, respectively.
   By using the relationships (<ref>) and (<ref>) and the extremal condition concerning Q, m, χ, the conjugate fields Q, m, χ can be represented via the macroscopic physical variables:
   
    χ = γ G”(-χ; J)/G'(-χ; J) - G”(-χ; J)χRSS  + -G”(-χ; J)γ + 2(G'(-χ; J))^2/G'(-χ; J) - G”(-χ; J)χσ^2,  
    Q = m = 2G'(-χ;J).

   Here, χ, G^'(-χ;J) and G^''(-χ;J) are given as follows:
   
    χ = -𝒮_J(z(-χ)) ,
    
       	G^'(-χ;J) = 1/2(z(-χ) + 1/χ), 
    
        G^''(-χ;J) = 1/2(z^'(-χ) + 1/χ^2),

   where z^'(-χ) is obtained from the derivative of equation (<ref>):
   
    z^'(-χ) = -[∫ρ_J(λ)/(z(-χ)-λ)^2dλ]^-1.

   
   The minimization problem in equation (<ref>) corresponds to the effective single body problem, which determines the value of the local magnetization ⟨ x_i⟩. 
   Splitting into effective single body problems from the original multi-body estimation problem is called the decoupling principle in the literature on information theory <cit.>.
   A comparison with the TAP/cavity analysis indicates that h_i≡mx_0,i + √(χ)z_i and m correspond to the local field and Onsager reaction coefficient, respectively <cit.>. 
   Here, z_i∼𝒩(0,1) effectively represents the randomness that comes from the observation matrix and measurement noise.
   Figure <ref> schematically shows the distribution of the local fields and how the local field determines the LASSO solution.
   Each local field is distributed according to the normal distribution 𝒩(mx_0,i, χ) and the LASSO solution is obtained by acting the soft-thresholding operator ST_λ, Q on it: 
   
    x_i^LASSO =  ST_λ, Q(h_i)≡h_i - λsgn(h_i)/QΘ(|h_i | - λ),

   where Θ(z) is Heaviside's step function.
   
   
  
  
  
	The LASSO solution takes a non-zero value if the amplitude of the corresponding local field is larger than λ. Conversely, if and only if it is smaller than λ, the LASSO solution is exactly zero.
    Hereafter, we call the non-zero and zero components of the LASSO solution the active and inactive components, respectively.
    
    The above observations indicate that once the local fields and m, χ are estimated from the LASSO solutions, one can construct an intended p-value P_i as 
    
    P_i ≡ 2{1- Φ(h_i/√(χ))} ,
    Φ(x) ≡∫_-∞^x e^-t^2/2/√(2π)dt,

    and an unbiased estimator as
    
    x_i^ debiased≡h_i/Q,

    with a confidence interval 
    
    ℐ_i(α_ sig) = [h_i/Q - Φ^-1(1-α_ sig/2)√(χ)/Q,  h_i/Q + Φ^-1(1-α_ sig/2)√(χ)/Q ],

    of significance α_ sig. 
    These are the key observations for the design of our scheme. 
   
   
    
    

 §.§ Adaptive TAP approach to constructing local fields and their variances from LASSO solutions

    

  §.§.§ Derivation of the adaptive TAP equations:

    To derive the relation between the LASSO solution x^LASSO(=⟨x⟩) and the local fields, let us consider Gibbs free energy:
  
    G(m) ≡ extr_h[
            	h^⊤m - 1/βln{
                	e^-β/2y - Ax_2^2 + βh^⊤x -βλx_1dx}].

    The average ⟨x⟩ is determined as the global minimizer of G(m):
    ⟨x⟩ = min_mG(m). 
    Once the above Gibbs free energy is exactly calculated, the extremal conditions of h and m generally associate the average ⟨x⟩ and local field <cit.>.     
    However, the evaluation of equation (<ref>) is computationally difficult in general.
    To overcome this difficulty, we take the following expectation consistent or the adaptive TAP approach <cit.>.
    First, we define an alternative Gibbs free energy:
    
    G(m, Q) ≡ extr_h, Λ[h^⊤m - N/2Λ Q -1/βln{∫ e^-β/2y - Ax_2^2 +βh^⊤x -β/2Λx_2^2 -βλx_1 dx}],

    which provides the constraints on the first and macroscopic second moments so that ⟨x⟩, ⟨ |x|^2 ⟩/N = min_m, Q G(m, Q). 

	Unfortunately, equation (<ref>) is also difficult to evaluate directly.
    The adaptive TAP approach resorts this calculation to the following approximation:
	
    G(m, Q) ≃ϕ_ada(m, Q) ≡ϕ (m, Q; l = 0) 
            + ϕ^G(m, Q; l=1) 
            - ϕ^G(m, Q; l=0), 
    ϕ (m, Q; l) ≡ extr_h, Λ{h^⊤m - N/2Λ Q 
                    - 1/βln∫ e^-β l/2y - Ax_2^2 + βh^⊤x - β/2Λx_2^2 -βλx_1  dx}, 
    ϕ^G(m, Q; l) ≡ extr_h^G, Λ^G{h_G^⊤m - N/2Λ_GQ 
                -1/βln∫ e^-β l /2y - A x_2^2+ βh_G^⊤x -β/2Λ_Gx_2^2  dx},

	where ϕ(m, Q;l=0), ϕ^ G(m, Q;l=1), and ϕ^ G(m, Q;l=0) are the free energies for the modified distributions: the first term is a factorized distribution but contains the original non-Gaussian prior factor e^-βλx_1, while the second and third terms are the global and factorized multivariate Gaussian distribution that replaces the prior factor e^-βλx_1 with a Gaussian factor e^-βΛ_ Gx_2^2/2.
    In contrast to the original form of Gibbs free energy (<ref>), adaptive TAP free energy (<ref>) can be easily calculated as it is composed of only integration over the multivariate Gaussian and factorized distributions.
    The evaluation of the integrals and extremal conditions in the second and third terms of equation (<ref>) provides the following expression of ϕ_ ada:
       
    ϕ_ ada(m, Q) =  extr_h, Λ[
                  1/2y - Am_2^2- NΛ Q/2
                  - N/β G(-χ;J) . 
    .
                  + h^⊤m - 1/2Λ∑_i=1^N (|h_i|- λ)^2 Θ(|h_i|-λ)
               ],

    where χ≡β(Q-q), q≡∑_i m_i^2/N.
    It has been shown <cit.> that the above free energy ϕ_ada(m, Q) is asymptotically consistent with replica theory in the sense that lim_β→∞, N→∞ extr_m, Qϕ_ ada(m, Q)_A,ξ/N = f_A,ξ when A is a sample from the rotationally invariant ensemble.
	Thus, the extremal condition on h, Λ, m, Q and linear response argument give the intended TAP/cavity equations, which connect the local field and LASSO estimator for the current matrix ensembles for β→∞,N→∞:
    
    h = Λm + A^⊤(y - Am), 
    
                m_i = h_i - λ sgn(h_i)/ΛΘ(|h_i| - λ), 
    Λ = 2G'(-χ), 
    χ = 1/NΛ∑_i=1^NΘ(|h_i| - λ) = ϱ_ active/Λ,

    where ϱ_ active≡∑_i=1^N Θ(|h_i| -λ)/N=|{i|x^ LASSO_i≠ 0}|/N is the active component density of the LASSO solution (<ref>).
    
    

  §.§.§ General construction procedure of the de-biased estimator, confidence interval, and p-value:
 
    In summary, once the LASSO estimator x^ LASSO(y,A;λ) is obtained for a set of (y, A; λ) by using versatile algorithms for the optimization problem (<ref>) such as least-angle regression <cit.>, coordinate descent <cit.>, and approximate message passing <cit.>, one can estimate the local fields h(y, A; λ), de-biased estimator x^ debiased(y, A; λ), confidence interval {ℐ_i(α_ sig)}_i, and p-value P_i as follows.
    We emphasize here that there is no need to use the derived TAP equation to obtain a LASSO estimator.
    
    First, the active component density ϱ_ active is calculated from the LASSO solution:
    
    ϱ_ active(y, A;λ) = 1/N| {i | x^ LASSO_i(y,A;λ)≠ 0}|.

    Second, z(-χ) is obtained by combining equations (<ref>), (<ref>), (<ref>), and (<ref>):
    z(-χ) is obtained as the solution of
    
    z = 1-ϱ_ active/𝒮_J(z).

	This equation is solved analytically or numerically depending on the cases. 
    Note that this equation is easily solved by using a simple iteration algorithm even if an analytical expression is not obtained.
    Then, z^'(-χ), χ , G^'(-χ;J), G^''(-χ;J), the Onsager coefficient Q=Λ, the local field h(y,A;λ), the de-biased estimator x^ debiased(y,A;λ), the residual sum of squares, and the variance of the local field χ are obtained by subsequently substituting the obtained values into equations (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>):
    
    z^'(-χ) = -[∫ρ_J(λ)/(z(-χ)-λ)^2dλ]^-1,
    χ = -𝒮_J(z(-χ)), 
    
        	   	G^'(-χ;J) = 1/2(z(-χ) + 1/χ), 
    
        G^''(-χ;J) = 1/2(z^'(-χ) + 1/χ^2),
    Q = Λ = z(-χ) + 1/χ, 
    h(y,A;λ) = Qx^ LASSO(y,A;λ) + A^⊤(y-Ax^ LASSO(y,A;λ)), 
    x^ debiased(y,A;λ) = h(y,A;λ)/Q,
     RSS = 1/My - Ax^ LASSO(y,A;λ) _2^2,
    χ = γ G”(-χ; J)/G'(-χ; J) - G”(-χ; J)χRSS  + -G”(-χ; J)γ + 2(G'(-χ; J))^2/G'(-χ; J) - G”(-χ; J)χσ^2.

  	Finally, the de-biased estimator's confidence interval {ℐ_i(α_ sig)}_i and p-value { P_i}_i are obtained based on equations (<ref>)–(<ref>).

	Note that a consistent estimator of the error variance σ^2 should be needed when it is unknown.





§ NUMERICAL EXPERIMENT
  
     
    

 §.§ Settings

    We perform numerical experiments to assess the usefulness of the proposed scheme.
    For this, we artificially generate the true sparse parameter x_0, observation matrix A, and measurement noise ξ.
    The true sparse parameter x_0 is generated from the Bernoulli–Gauss distribution: x_0,i∼_ i.i.d. (1-ϱ)δ(x_0,i) + ϱ𝒩(0,1) for i=1,2,...,N.
    The measurement noise is distributed according to the Gaussian distribution ξ∼𝒩(0_M, σ^2 I_M).
    For the random observation matrix ensembles, the following ensembles are considered.
    
	    
  * The random i.i.d. Gaussian ensemble in which all entries of A are i.i.d. Gaussian variables with mean 0 and variance 1/N. For this ensemble, the asymptotic eigenvalue distribution is given as the Marchenko-Pastur distribution <cit.>: 
        
    ρ(s) = (1-γ)δ(s) + γ/2π√((λ_+ - s)(s-λ_-))/s𝕀_[λ_-, λ_+](s), 
    λ_± = (1±√(γ))^2 , 
    𝕀_S(x) = 1   if x∈ S
    
          0   otherwise 
    .

         Then, the form of G^'(-χ;J), G^''(-χ;J), χ and Q are given as follows:
        
    G^'(-χ;J) = γ/21/1+χ, 
    
                G^''(-χ;J) = γ/21/(1+χ)^2, 
    χ = ϱ_ active/γ - ϱ_ active,
    Q = γ - ϱ_ active.

        By substituting the above expressions of G^', G^'' into (<ref>), one can show that χ does not depend on σ^2. 
        This is the characteristic property of this ensemble.
        Generally, χ depends on the measurement noise σ^2.
        
        
  * The row-orthogonal ensemble <cit.> constructed by randomly selecting M rows from a randomly generated N× N orthogonal matrix. For this ensemble, the asymptotic eigenvalue distribution is given as ρ(s) = (1-γ)δ(s) + γδ(s-1). In this case, the form of G^'(-χ;J), G^''(-χ;J), χ and Q are given as follows:
        
    G^'(-χ;J) = 1/2(z(-χ) + 1/χ),
    
                G^''(-χ;J) = 1/2(z^'(-χ) + 1/χ^2),
    χ = ρ_ A(1-ϱ_ active)/γ - ϱ_ active,
    Q = γ - ϱ_ active/1-ϱ_ active,

          where 
          
    z(-χ) = - 1- χ + √((χ+1)^2- 4γχ)/2 χ,
    
                z^'(-χ) = - 1- 2 γχ + χ + √((χ+1)^2- 4 γχ)/2 χ^2√((χ + 1)^2 - 4 γχ).
 

        
  * The random discrete cosine transform (DCT) ensemble in which A is constructed by randomly selecting M rows from N× N DCT matrix.
       While this ensemble shares the same eigenvalue distribution as the row-orthogonal one, it is much more relevant for practical purposes, as the computational cost for observation and inference can be significantly reduced by using the fast Fourier transform technique. 
        In addition, although the rotationally invariant assumption on O does not hold, this ensemble is also compatible with the current adaptive TAP scheme, as pointed by <cit.>.
        
        
  * The geometric setup <cit.> in which A is constructed as A=UΣ V^⊤, where U∈ℝ^M× M and V∈ℝ^N× N are random samples from the uniform distribution of orthogonal matrices, and Σ∈ℝ^M× N is a diagonal matrix whose (i,i)-th element is given by ν_i∝τ^i-1 for i=1,2,...,M. 
        The parameter τ∈ (0,1] is chosen so that the given value of the peak-to-average eigenvalue ratio 
        
    κ≡ν_1^2/M^-1∑_i=1^Mν_i^2

        is met and the singular values are scaled to satisfy the power constraint 1 = 1/N∑_i=1^M ν_i^2.
        The asymptotic eigenvalue distribution is given as 
        
    ρ(s)=(1-γ)δ(s) + γ/η s𝕀_(Be^-η, B](s),

        where η and B are related to the peak-to-average ratio κ:
        
    κ = η/1-e^-η, 
    
                B = κ/γ.

        
        In this case, the explicit form of G^' ,G^'' cannot be obtained.
        Thus, it should be evaluated numerically.
        To achieve this aim, we conduct the procedure explained in section <ref>, using the expression of the Stieltjes transform and z^'(-χ):
        
    χ = -𝒮_J(z(-χ))=-∫ρ(s)/z(-χ)-sdλ = -1/z(-χ)[1-α/ηlnz(-χ)-B/z(-χ)-Be^-η], 
     
              z^'(-χ) = z(-χ)^2/-1+γ/ηlnz(-χ)-B/z(-χ)-Be^-η -z(-χ)/(z(-χ)-Be^-η)(z(-χ)-B).

        
    
    
    We mainly use the random i.i.d. Gaussian ensemble and random DCT ensemble for the numerical experiments. The geometric setup is only used in section <ref>. We do not use the original row-orthogonal setup.
    
	Once a tuple of (x_0, A, ξ) is generated, we calculate x^ LASSO, h, χ, χ and Q = Λ, x^ debiased by using the procedure explained in section <ref>.
        
 	   
   
	To estimate the error variance σ^2 needed in the random DCT case, we use the naive cross-validation-based estimator:
    
    σ̂^2(y, A;λ̂) ≡1/M - N ϱ_ activey - Ax^ LASSO(y, A;λ̂)_2^2,

    where λ̂ is selected by K-fold cross-validation.
	In <cit.>, it is empirically shown that this estimator robustly estimates the error variance, more so than its competitors.
    
    We use N_ s = 1000 different sets of pairs (A, ξ) for fixed x_0 to evaluate the statistical properties of the observables.
    We set ϱ=0.1, γ=0.5, σ^2=0.02, and K=40, except for the geometric setup. 
    In the geometric setup, we set ϱ =0.1, γ=0.8, σ^2=0.02, and κ=8.

	
    

 §.§ Results

    	

  §.§.§ Distribution of the local fields and de-biased estimators:

        First, we examine the statistical properties of the local fields and de-biased estimators.
        Figure <ref> plots the sample quantiles of {(h_i - Qx_0,i)/√(χ)}_i versus the theoretical quantiles of the standard normal distribution for one configuration of (x_0, A,ξ).
        It is clear that all the points are close to the line with unit slope and zero intercept.
        Further, Figure <ref> plots the average values of the x^ LASSO and x^ debiased versus the true parameter x_0.
        In contrast to the LASSO estimators, which are shrunk toward zero by the regularization term, x^ debiased efficiently reduces the LASSO estimator's bias. 
        The average is taken over N_s realizations of (A,ξ).
        These results validate our theoretical predictions on the local fields and de-biased estimators.
        Figure <ref> plots the constructed de-biased estimators and their 95% confidence intervals.
        We show only the first 80 components for the sake of clarity.
        
        
    
    
     Although Figures <ref>–<ref> show the results for one value of λ, 
     the same results are obtained for a wide range of λ.
     The means of {h_i - Qx_0,i}_i and {x^ debiased - x_0,i}_i are zero in both the i.i.d. Gaussian and the random DCT cases (Figure <ref> (a) and (b)).
     Further, the variances of { h_i-Q x_0,i}_i and {x^ debiased_i -x_0,i}_i agree with their estimates of χ and χ/Q^2, respectively for the whole range of the weight of the L_1 regularizer λ (Figure <ref> (c) and (d)).
    
    
     
     

  §.§.§ Hypothesis testing:

    An important advantage of the proposed scheme over LASSO is that it provides a hypothesis testing method with a null hypothesis that a certain parameter vanishes.
    Although LASSO provides a parameter selection rule that selects an active component set 𝒜(y, A;λ) as 𝒜(y, A;λ)={i | x^ LASSO_i(y, A;λ)≠0 }, it cannot measure the statistical significance for finding an active component.
    
     
    

     
    
    
    Specifically, we are interested in testing an individual null hypothesis H_0,i: x_0,i=0 versus the alternative hypothesis H_1,i : x_0,i≠ 0, assigning a p-value of P_i for these tests.
    To this end, we evaluate the p-value of { P_i} by using equation (<ref>) for a two-tailed test.
    Then, the decision rule is to reject the null hypothesis H_0,i if the observed p-value P_i is lower than α̃_ sig and to accept the alternative hypothesis otherwise:
    
    T_i(y, A;λ) = 1   if P_i ≤α̃_ sig (reject)
    
          0   otherwise (accept)
    ,

    where α̃_ sig is the significance level.
    We use T as a rejection flag.
    This procedure ensures that the type I error probability or the FPR is α̃_ sig.
    Here, the FPR is the probability of falsely rejecting the null hypothesis H_0,i:
    
    FPR≡| {i | T_i =1 and x_0,i=0}| /| {i | x_0,i=0}|.

    Indeed, Figure <ref> shows that the significance level α̃_ sig and empirical true positive rate (TPR) are in excellent agreement.
    
    Further, we examine the relation between the FPR and TPR or the statistical power achieved by LASSO and our hypothesis testing procedure.
    Here, the TPR is the probability that the test correctly rejects the null hypothesis H_0,i:
    
    TPR≡| {i | T_i =1 and x_0,i≠ 0}| /| {i | x_0,i≠ 0}|.

    Note that although we can control the FPR by varying the significance level α̃_ sig, the TPR cannot be controlled. 
    Thus, a performance measure of the variable selection procedure by hypothesis testing can be given as the TPR for each value of the FPR.
    We evaluate the performance of hypothesis testing by using the ROC curve, which plots the TPR versus the FPR as an implicit function of α̃_ sig.
    We examine the TPR and FPR by varying the significance level α̃_ sig for each regularization parameter λ.
    For comparison purposes, we also plot the ROC curve for LASSO.
    For LASSO, the TPR and FPR are examined by changing the regularization parameter λ.
    
    
    
 
    
    Figure <ref> summarizes the results averaged over N_s configurations of (A,ξ).
    It is observed that for some values of λ around which the variance of the de-biased estimator is minimized, our testing procedure performs slightly better than LASSO in the sense that the TPR of the testing method is slightly larger than that of LASSO's one for some values of the FPR.
    In the case of LASSO, when the measurement ratio γ is sufficiently small, the TPR and FPR do not coincide with (1, 1) for finite λ>0, as the consistency property does not hold in such a situation and the number of active components of the LASSO estimator is always smaller than min(N,M) <cit.>.
    On the contrary, as our hypothesis testing procedure always approaches the point (1,1) from (0,0), we can examine the TPR for all the values of the FPR ∈[0,1].
    The superiority of the TPR comes from the fact that we are using the knowledge of the ensemble of the observation matrix.
    Further, as the hypothesis testing procedure controls the FPR and TPR by varying the significance α_ sig but not λ, one does not suffer from the shrinkage effect in the variable selection procedure.
    This is another advantage over variable selection by LASSO.
    These observations show the utility of our hypothesis testing procedure.

    
    

		

  §.§.§ Hyperparameter selection via confidence interval minimization:

        The issue of hyperparameter selection is noteworthy here.
        As LASSO has the hyperparameter λ that controls the strength of the regularization, one should choose a value of λ based on some criteria.
        As shown in Figure <ref>, the estimated variance of the de-biased estimator χ/Q^2 has a minimum value at some λ>0.
        At this point, one can estimate x_0 with the highest conviction in the sense that the confidence interval has the smallest width.
        It is therefore expected that the estimated variance of the de-biased estimator itself serves as a hyperparameter selection criterion.
        Indeed, in the i.i.d. Gaussian and row-orthogonal/random DCT cases, one can analytically show that minimizing the confidence interval is the equivalent to the minimization of the leave-one-out cross-validation error 𝒞:
        
    χ/Q^2 = 1/γ𝒞    the i.i.d. Gaussian case, 
    1-γ/γ𝒞 + σ^2     the row-orthogonal or the random DCT cases.

        In other words, the leave-one-out cross-validation error and width of the confidence intervals are related with the linear transformation in these cases (Figure <ref>).

        Here, the leave-one-out cross-validation error is a widely used hyperparameter selection criterion that evaluates prediction performance, defined as follows:
        
    𝒞(y, A;λ) = 1/M∑_i = 1^M1/2 y_i -  a_i^⊤x^ LASSO(y_\ i, A_\ i;λ)_2^2,

        where the symbol \ i denotes the absence of the i-th component (e.g., a_\ i=(a_1,...,a_i-1,a_i+1,...,a_N)^⊤) and each term in the summation evaluates the fitness to the i-th data when the true signal is inferred from the other data.
        In the settings considered here, the above leave-one-out cross-validation error is expressed as follows <cit.>:
        
    𝒞 = (1-ϱ_ active/γ)^-2 RSS = (1-2χ G^'(-χ;J)/γ)^-2 RSS.

        By substituting the expression of the leave-one-out cross-validation error (<ref>) into equation (<ref>), the relations (<ref>) are obtained. 


        To investigate the validity of the above observation that the confidence interval minimization and leave-one-out cross-validation error minimization provide the same λ, we test the geometric setup case in which χ/Q^2 is not expressed as a linear function of 𝒞.
        Figure <ref> compares the variance of {x^ debiased - x_0,i}_i with the leave-one-out cross-validation error (<ref>).
        Surprisingly, the minimization of these two quantities seems to provide the same value of λ, although they do not have a functional relation as (<ref>).

        From the above observations, we speculate that the minimization of the confidence interval proposed here and the minimization of the leave-one-out cross-validation error yields the same value of λ for LASSO in general rotationally invariant observation matrices, but further investigation in this direction is still needed.







§ SUMMARY

	We developed a new computationally feasible scheme for de-biasing and uncertainty estimation in LASSO in the case of rotationally invariant observation matrix ensembles and validated the proposed scheme by using numerical experiments.
    We focused on the development of a de-biased estimator that has a confidence interval and hypothesis testing scheme for the null hypothesis that a certain parameter vanishes.
	The numerical experiments showed that the proposed method efficiently constructed de-biased estimators with confidence intervals and p-values for the intended hypothesis testing. 
    We revealed that the proposed hypothesis testing slightly improved the variable selection performance in the sense that the TPR of the testing method achieves a slightly larger value than that of the LASSO's one for some values of the FPR.
    Further, we examined the utility of the estimator of the confidence interval as a criterion for determining the hyperparameter. Surprisingly, minimizing the width of the confidence interval was equivalent to the minimization of the leave-one-out cross-validation error in our investigation.

    
    Although we only focused on LASSO for linear models, future work could include an extension to other sparse regression methods such as the elastic net <cit.> as well as generalized linear models.




§ ACKNOWLEDGEMENT

Support by JSPS KAKENHI Nos. 25120013 and 17H00764 (YK) is acknowledged. 






§ DERIVATION OF THE FREE ENERGY DENSITY
  
	To take the average that appears in (<ref>), we use the replica method <cit.> based on the identity for n∈ℝ:
    
    f = -lim_β→∞lim_N→∞1/β Nlim_n→ 0Z^n_A,ξ/n.

    In the replica method, we first take the average of the n-th power of the partition function over the randomness of A, ξ for the positive integer n∈ℕ, and then analytically continue the obtained expression to real n∈ℝ to take the limit n→0, exchanging the order of the limits.
    
    For the general matrix ensembles considered here, it is convenient to first take the average over ξ.
    By taking this average, we obtain the following expression under the replica symmetric ansatz:
    
    Z^n_A,ξ  = ∫exp( 1/2Tr J L )_A  e^S dQdqdm,

      where L, u_a, and S are defined as follows:
      
    L ≡β^2σ^2/ 1 + β nσ^2(∑_a u_a)(∑_a u_a)^⊤
              -β∑_a u_a u_a^⊤, 
    u_a ≡x_a - x_0, 
     e^S ≡∫∏_a=1^nδ(NQ - x_a^⊤x_a)
          δ(Nm - x_a^⊤x_0) 
    ×∏_1≤ a<b≤ nδ(Nq - x_a^⊤x_b)exp{-Nγ/2β nσ^2 -βλ∑_ax_a_1 }dx,

    where u_a and x_a are the a-th replica's variable.
    In <cit.>, it was shown that under the rotational invariance assumption on the random matrix J=A^⊤ A for eigenvalue decomposition J=ODO^⊤ considered in this study, the average over A that appears in equation (<ref>) is evaluated by using the eigenvalues {s_i}_i of L/N for sufficiently large N:
    
    exp( 1/2Tr J L )_A = exp{N ∑_iG(s_i; J) },

    where G(x;J) is the function defined in (<ref>).
    Under the replica symmetric ansatz, L/N has three types of eigenvalues: s_1=βΔ Q - β n (q- 2m + ϱ) + n β^2 σ^2 Δ Q, s_2=-βΔ Q, and s_3 =0. 
    The number of degeneracy is 1, n-1, and N-n, respectively.
    Thus, we obtain the following expression up to the leading order in n:
    
    e^1/2Tr J L _A 
            = exp[-Nnβ{-G(-βΔ Q; J)/β..
    ..+ G'(-βΔ Q; J)(q-2m + ϱ - βΔ Q σ^2) }].

    
    On the contrary, by using the Fourier transform of the delta function and Hubbard–Stratonovich transform: e^B^2/2A=∫ e^-Ax^2/2 + Bx√(A/2π)dx for A>0,B∈ℝ, the factor e^S is given as follows:
    
    e^S = ∫exp[Nn{γσ^2/2 + 
               qq/2 + QQ/2 - mm.. 
    .. +1/N∑_i=1^N ∫lnϕ(x_0,i, z_i, Q,q, m;β, λ)Dz_i
               }]dQdqdm,
    ϕ(x_0,i, z_i, Q,q, m;β, λ) =
             ∫exp{-Q + q/2x_i^2 + (mx_0,i + √(q)z_i)x_i - βλ |x_i| }dx_i.

    For β→∞, the relevant variables scale as β(Q-q)=χ, Q + q = βQ, m=βm, and q=β^2χ of order unity to ensure an appropriate limit f exists.
    Finally, by combining equations (<ref>–<ref>) and evaluating the integrals by adopting the saddle point method, we obtain equation (<ref>) for β, N→∞.
    



§ REFERENCES

100

D2006 Donoho D L, 2006 Compressed sensing IEEE Trans. Inform. Theory 52 1289.
CT2006 Candes E J and Tao T, 2006 Near-optimal signal recovery from random projections: Universal encoding strategies?, IEEE Trans. Inform. Theory 52 5406.


Petal2010 Peng J, Zhu J, Bergamaschi A, Han W, Noh J-Y, Pollack J R, and Wang P, 2010 Regularized multivariate regression for identifying master predictors with application to integrative genomics study of breast cancer, Ann. Appl. Stat. 4(1) 53–77.


T1996 Thibshirani, R, 1996 Regression shrinkage and selection via the Lasso, J. R. Stat. Soc. B 58 267–288.


CT2005 Candes E J and Tao T, 2005 Decoding by linear programming, IEEE Trans. Inform. Theory 51 4203.


BRT2009 Bickej PJ, R Yitov and Tsybakov A B, 2009 Simultaneous analysis of lasso and dantzig selector, Ann. Stat. 37 1705–1732.


DMM2009 Donoho D L, Maleki A and Montanari A, 2009 Message-passing algorithms for compressed sensing, Proc. Natl. Acad. Sci. 106, 18914.

KWT2009 Kabashima Y, Wadayama T and Tanaka T, 2009 A typical reconstruction limit for compressed sensing based on L_p-norm minimization, J. Stat. Mech. L09003.

BS2011 Geer S V and Bühlmann P, 2011 Statistics for High-Dimensional Data: Methods,
Theory and Applications, Springer Series in Statistics. Springer, Berlin,
Heidelberg.


OW Opper M and Winther O, 2001 Tractable approximations for probabilistic models: The adaptive Thouless-Anderson-Palmer mean field approach, Phys. Rev. Lett 86 3695.
OW2Opper M and Winther O, 2001 Adaptive and self-averaging Thouless-Anderson-Palmer mean-field theory for probabilistic modeling, Phy. Rev. E 64 056131.
OW3Opper M and Winther O, 2005 Expectation consistent approximate inference, J. Mach. Learn. Resear. 6 2177.
MMinka P T, 2001 Expectation propagation for approximate Bayesian inference, in Proceedings UAI-2001, 362.




 

JM2014 Javanmard A and Montanari A, 2014 Confidence intervals and hypothesis testing for high-dimensional regression, J. Mach. Learn. Resear., 15, 2869–2909.
GBRD2014Geer S V, Bühlmann P, Ritov Y and Dezeure R, 2014 On asymptotically optimal confidence regions and tests for high-dimensional models, Ann. Stat., 42: 1166–1202.
ZZ2014Zhang C and Zhang S, 2014 Confidence intervals for low dimensional parameters in high dimensional linear models, J. R. Stat. Soc. Ser. B. Stat. Methodol, 76 217–242. MR3153940.


V1986 Voiculescu D, 1986 Addition of certain non-commuting random variables, J. Funct. Anal. 66 323.


GV2006 Guo D and Verdú S, 2006 Randomly spread CDMA: Asymptotics via statistical physics, IEEE Trans. Inform. Theory, 51 1983.
RFG2009 Rangan S, Fletcher A K and Goyal V K, 2009 Asymptotics Analysis of MAP Estimation via the Replica Method and Compressed Sensing, Advances in Neural Information Processing Systems 22 (Curran Associates, Inc., Red Hook), pp. 1545–1553.


P1982 Plefka T, 1982 Convergence condition of the TAP equation for the infinite-ranged Ising spin glass model, J. Phys. A: Math. Gen., 15, 1971–1978.


OS2001 Opper M and Saad D, editors, 2001 Advanced Mean-field Methods: Theory and Practice. MIT
Press.


KV2014 Kabashima Y and Vehkaperä M, 2014 Signal recovery using expectation consistent approximation for linear observations, in Proc. IEEE Int.
Symp. Inf. Theory, Jun./Jul. pp. 226–230.


RTF2013 Reid S, Tibshirani R and Friedman J, 2013 A study of error variance estimation in lasso regression. Preprint. arXiv:1311.5274.


CO2018 Çakmak B and Opper M, 2018 Expectation propagation for approximate inference: Free probability framework. preprint. arXiv:1801.05411.


EHJT2004 Efron B, et al., 2004 Least angle regression, Ann. Stat., 32, 407–499.
FHT2010 Friedman J, Hastie T, and Tibshirani R, 2010 Regularization paths for generalized linear models via coordinate descent, J. Stat. Softw. 33, 1.
R2011 Rangan S, 2011 Generalized approximate message passing for estimation with random linear mixing. Information Theory Proceedings (ISIT), 2011 IEEE International Symposium on. IEEE.


MP1967 Marčenko V A and Pastur L A, 1967 Distribution of eigenvalues for some sets of random matrices. Mathematics of the USSR-Sbornik 1.4 (1967): 457.


VKC2016 Vehkaperä, M, Kabashima Y, and Chatterjee S, 2016. Analysis of regularized LS reconstruction and random matrix ensembles in compressed sensing, IEEE Trans. Inform. Theory 62.4 2100-2124.
RFSK2017 Rangan S, et al., 2017 Inference for generalized linear models via alternating directions and Bethe free energy minimization, IEEE Trans. Inform. Theory, 63, 676–697.


MPV1987 Mézard M, Parisi G and Virasoro M, 1987 Spin Glass Theory and Beyond: An Introduction to the Replica Method and Its Applications, World Scientific,
Singapore.


OK2016Obuchi T and  Kabashima Y, 2016 Cross validation in LASSO and its acceleration, J. Stat. Mech. Theory Exp., 5, 053304.
RM2018 Rad K R and Maleki A, 2018 A scalable estimate of the extra-sample prediction error via approximate leave-one-out. arXiv preprint arXiv:1801.10243.


ZH2005 Zou H and Hastie T, 2005 Regularization and variable selection via the elastic net, J. R. Stat. Soc. B, 67, 301.


MPR1994 Marinari E, Parisi G and Ritort F, 1994 Replica field theory for deterministic models: II. A non-random spin glass with glassy behavior, J. Phys. A: Math. Gen, 27, 7647.

