Estimator of prediction error based on approximate message passing for penalized linear regression

We propose an estimator of prediction error using an approximate message passing (AMP) algorithm that can be applied to a broad range of sparse penalties. Following Stein’s lemma, the estimator of the generalized degrees of freedom, which is a key quantity for the construction of the estimator of the prediction error, is calculated at the AMP fixed point. The resulting form of the AMP-based estimator does not depend on the penalty function, and its value can be further improved by considering the correlation between predictors. The proposed estimator is asymptotically unbiased when the components of the predictors and response variables are independently generated according to a Gaussian distribution. We examine the behavior of the estimator for real data under nonconvex sparse penalties, where Akaike’s information criterion does not correspond to an unbiased estimator of the prediction error. The model selected by the proposed estimator is close to that which minimizes the true prediction error.


Introduction
In recent decades, variable selection using sparse penalties, referred to here as sparse estimation, has become an attractive estimation scheme [1,2,3].The sparse estimation is mathematically formulated as the minimization of the estimating function associated with the sparse penalties.In this paper, we concentrate on the linear regression problem with an arbitrary sparse regularization.Let y ∈ R M and A ∈ R M ×N be a response vector and predictor matrix, respectively, where each column of A corresponds to a predictor.We denote the regression coefficient to be estimated as x ∈ R N , and the problem considered here can be formulated as min where J(x; η) = selection problem, the variables are estimated by a maximum likelihood method under the constraint on the support of the variables used in the model.In sparse estimation, the support is determined by the regularization parameter η; hence, the adjustment of η is regarded as the selection of a model, and is crucial in the sparse estimation.
In general, the regularization parameter is determined based on the model selection criteria.One of the criteria is the prediction performance, which is measured by the prediction error defined for the test data generated according to the true distribution.Following this criterion, the model with the regularization parameter that minimizes the prediction error provides an appropriate description of the data.However, the exact computation of the prediction error is impossible, because we do not know the true distribution.Therefore, the model selection is generally implemented using estimators of the prediction error, rather than the prediction error itself.For the maximum likelihood estimation, Akaike's information criterion (AIC) is an unbiased estimator of the prediction error when the true generative model is included in the set of models under consideration, and is widely used for model selection when concentrating on the prediction of unknown data [4].However, for sparse estimation, AIC is generally not an unbiased estimator of the prediction error.In such cases, a C p -type criterion, which evaluates the prediction error using generalized degrees of freedom (GDF) [5,6], is useful when the variance of the data is known.However, the exact computation of GDF also requires the true probability distribution of the data, which means that an estimator must be constructed.In sparse estimation, GDF and its unbiased estimator can be derived for certain regularizations.A well-known result obtained using the least absolute shrinkage and selection operator (LASSO) is that AIC corresponds to the unbiased estimator of the prediction error [7], a nontrivial fact that has attracted the interest of statisticians.Recently, sparse estimations using nonconvex regularizations have been studied with the aim of achieving a more compact representation of the data than LASSO [8,9].Hence, model selection methods that can be applied to such regularizations are required.
In this paper, we propose a numerical method for the construction of estimators for GDF and prediction error based on approximate message passing (AMP).The derived form of GDF does not depend on the regularization, although the property of the AMP fixed point is regularization-dependent.The estimator is asymptotically unbiased when the predictor matrix and response variables are independently and identically distributed (i.i.d.) Gaussian, which is indicated by showing the correspondence between AMP and the replica method [10].In this case, we can confirm that AIC is the unbiased estimator for the 1 penalty, which is consistent with an earlier study [7].We apply our estimator to real data and show that it adequately emulates model selection according to the prediction error.
The remainder of this paper is organized as follows.In Sec. 2, we summarize the background materials that support the following discussion.The problem setting considered in this paper is explained in Sec. 3, and the AMP algorithm for the penalized linear regression problem is introduced in Sec. 4. The estimator of GDF using the AMP fixed point is derived in Sec. 5, and its asymptotic property for Gaussian random data and predictor are studied in Sec.6 using the replica method.In Sec. 7, the proposed estimator is applied to real data, and the improvement it offers is quantified in Sec. 8. Finally, Sec. 9 presents the conclusions to this study and a discussion of the results.

Prediction error and generalized degrees of freedom
We denote the estimated regression coefficient obtained by solving (1) and its fit to the response variable y as x(y) and ŷ = Ax(y), respectively.The prediction performance of the model under the given response variable is measured by the prediction error where z ∈ R M represents test data whose statistical property is equivalent to that of y.
The value of η that minimizes the prediction error gives an appropriate description of the data y.We aim to approximate the prediction error by constructing its estimator.The simplest estimator is the training error, which is defined by but this underestimates the prediction error because of the correlation between y and ŷ.We consider the case in which the mean and covariance matrix of the response vector y are given by where T denotes the matrix transpose, and I M is the M -dimensional identity matrix.In this case, the prediction error and the training error satisfy the relationship where df is the generalized degrees of freedom (GDF) defined by [11] df = cov(y, ŷ(y)) and cov(y, ŷ . The expression ( 5) is known as the C p criterion for model selection [5,6].

Stein's lemma
The covariance (6) is not observable in the sense that the expectation with respect to the true distribution is required.Hence, estimators of df, denoted by df(y), are used instead of df for approximating the prediction error.An idea for the construction of the estimator appears in Stein's unbiased risk estimate (SURE), known as Stein's lemma [12].Following Stein's lemma, the component-wise covariance between the response variable and its fit is given by Equation (7) means that an unbiased estimator of GDF is given by Note that ( 7) is exact when the response variables are i.i.d.Gaussian random variables.Equation ( 8) is an observable quantity, and hence one of the estimators of the prediction error can be defined by Equation ( 9) is an unbiased estimator of the prediction error when the response variables are i.i.d.Gaussian, but the unbiasedness is not generally guaranteed for other response variables.
For the maximum likelihood method, the unbiased estimator of GDF is given by the number of parameters used in the model, which is the AIC given by [4] when the true generative model is included in the assumed models.As defined here, AIC has been multiplied by σ 2 y /M compared with its conventional definition, but this does not affect the following discussion.For penalized regression problems, when the fit of the response variable is represented as ŷ = Hy using a matrix H, it can be shown that GDF corresponds to TrH [13], which is widely used in ridge regression.Furthermore, many studies have derived SURE in problems where (8) is easily computed [14,15,16].

GDF for sparse penalties
For sparse penalties, the derivation of GDF and the construction of its estimator is not straightforward, because the penalty term is not differentiable at the origin.For the 1 penalty, the GDF is equal to the expectation of the density of non-zero components [7]; and hence || x(y)|| 0 /M is an unbiased estimator of GDF.It means that AIC corresponds to an unbiased estimator of the prediction error.The extension of this result for general problem settings has been extensively discussed [17,18].However, such a simple relationship does not generally hold, and efficient computation methods for (8) are required.Computation techniques for df have been developed using cross-validation and the parametric bootstrap method, but these techniques are computationally expensive [6,19].In this paper, we propose a method to construct the estimator of the prediction error based on Stein's lemma using AMP.In principle, the method can be applied for any J(x; η), and the extra computational cost of obtaining the fixed points is less than in other methods.

Sparse penalties considered in this paper
The penalty functions considered in this paper are 1 , the smoothly clipped absolute deviation (SCAD), and the minimax concave penalty (MCP).Under these regularizations, we explain the properties of the estimator using the one-dimensional problem θ(w) = arg min where w is the training sample.The ordinary least square (OLS) estimator, which minimizes the squared error between θ and w, is θOLS (w) = w for any w.The penalty term changes the behavior of the estimator from that of the OLS estimator.The penalties considered here give analytical solutions of (12).As explained in Sec. 4, the estimates given by AMP reduce to the one-body problem (12) with effective σ 2 w and w.We summarize the solution of (12) in the following subsections.

1 penalty
The 1 penalty is given by where λ is the regularization parameter determined by the model selection criterion.The solution of (12) under the 1 penalty is given by θ where w = w/σ 2 w and The behavior of the estimator at λ = 1 is shown in Figure 1(a).Across the whole region of w, the 1 estimator is shrunk from the ordinary least square (OLS) estimator, denoted by the dashed line.

MCP
The MCP is defined by [9] where η = {λ, a}.The estimator of (12) under MCP is given by θ(w; where Figure 1(c) shows the behavior of the MCP estimator at a = 3.The MCP estimator behaves like the OLS estimator when |w| > aλσ −2 , and is connected from zero to the OLS estimator in the region aλσ −2 ≥ |w| > λ.We call this the transient region of MCP.

Problem setting
Problem (1) discussed in this paper is equivalent to searching the ground state, but here we extend the problem to finite temperatures by introducing the inverse temperature β.We construct the posterior distribution for finite β, and extract the ground state corresponding to the maximum a posteriori (MAP) estimator in the limit β → ∞.
The likelihood we consider here is the Gaussian distribution given by The distribution (25) depends on the parameter x only through the form of Ax; hence, we hereafter denote (25) as P l (y|u(A, x)), where u(A, x) = Ax.The prior distribution of x is given by the penalty function and the inverse temperature as and hence the posterior distribution of x is given by following Bayes' formula, where Z β (y) is the normalization constant.The estimate for the solution of ( 1) is expressed as where • β denotes the expectation according to (27) at β.As β → ∞, (27) becomes the uniform distribution over the minimizers of (1).Hence, when problem (1) has a unique minimizer, it corresponds to the estimate (28).The exact computation of the expectation in (28) is intractable.Thus, we resort to the approximated message passing (AMP) algorithm to approximate (28), a technique that is widely used in signal processing and statistical physics [20,21,22,23].

Approximate message passing
In this section, we briefly summarize AMP for penalized regression problems.The explanation here only covers the aspects required to reach the result shown in Sec.5; a detailed derivation is given in [24].In the derivation of AMP, we introduce the following assumptions for sufficiently large M and N .

A2: The correlation between predictors is negligible: the off-diagonal components of
The message passing algorithm defined for the probability distribution P (x|y, A) is expressed by 2M N messages that are propagated from factors to variables mµ→i and variables to factors m i→µ according to mµ→i Here, V(µ) and F(i) represent the variable nodes and factor nodes connected to the µ-th factor node and i-th variable node, respectively, and \i denotes that i is not included.Using the messages, the marginal distribution is given by Updating the messages in the form of the probability distribution is intractable, and we represent the distributions using the mean and rescaled variance, which is the variance multiplied by β, assuming that they can be defined: Further, we define the mean and rescaled variance with respect to the marginal distribution as where the i-th component of the estimate (28), denoted by xi , corresponds to a i at β → ∞.Following the calculation under assumptions A1 and A2, the marginal distribution is derived as [24] m where and Ẑ(Σ 2 i , R i ) is the normalization constant required to satisfy dxM(x; Σ 2 , R) = 1.Equation (36) suggests that the product of messages from factors connected to i, ν∈F (i) m(y ν ), corresponds to a Gaussian distribution with mean R i and variance Σ 2 i .The functions g out and g out are given by where we define A µi a i→µ (42) and the function Ξ(y µ |ω µ , V µ ) contains the dependence on the likelihood as Under assumption A1, the following relationship holds [24]: which leads to

Specific form for penalized linear regression
Substituting the output function for linear regression given by into Eqs.( 40) and (41), we obtain the specific form of g out and g out as Using these equations, the mean and variance parameters are given by In summary, the AMP algorithm consists of (52)-(57).
where we define and note that At β → ∞, the integral of (58) can be computed by the saddle point method as which corresponds to the one-dimensional problem (12) with the replacements σ 2 w → Σ 2 and w → R.Under the penalties considered here, namely, 1 , SCAD, and MCP penalties, these functions are computed analytically, and have the form where R = R/Σ 2 and the subscript r denotes dependence on the regularization: r ∈ { 1 , SCAD, MCP}.

AMP-based estimator
At the AMP fixed point, Stein's lemma implies that the covariance between the response variable and its fit can be expressed as The dominant dependence on y ν of the estimate a i is induced through the mean parameter R i .Hence, using (60), Substituting the definition of ω µ , (42), into (55), we obtain where the second term can be ignored under the assumption A1: Hence, we have where a k→κ is approximated by a k and ( 65) is used to derive (68) from (67).Substituting (68) into (65), the second term of (68) can be ignored, and we obtain This means that df is the AMP expression of the unbiased estimator of GDF, and we define the corresponding estimator of the prediction error as

(y). (71)
These expressions do not depend on the form of the function J(x; η), hence they can be applied to any regularization J(x; η), such as elastic-net.In deriving (70), note that assumptions A1 and A2 are required, and hence its unbiasedness is not generally guaranteed.

Asymptotic behavior of the AMP-based estimator for Gaussian predictor matrix
When the components of the predictor matrix and data are i.i.d.according to a Gaussian distribution, the asymptotic property of the AMP fixed point can be described by the replica method, where N → ∞ and M → ∞ while α = M/N remains O(1).Here, we concentrate on the case y µ ∼ N (0, σ 2 y ) and A µi ∼ N (0, M −1 ), but the discussion can be extended to the case in which the generative model of y contains a "true" sparse expression x 0 such as y = Ax 0 + ξ, where ξ ∈ R M is a noise vector.It is known that the Gaussian i.i.d.predictor matrix is a full-rank matrix with probability 1; the rank of A is equal to min(N, M ).

Replica analysis
The basis of replica analysis is the derivation of the free energy density, which is defined by The free energy density and physical quantities derived from it are averaged over the predictor matrix and response variable, a procedure that is less common in the context of statistics.This averaging is implemented with the purpose of describing the typical behavior of the problem under consideration.
After the calculation under the replica symmetry assumption, the free energy density is derived as [24] f = extr where extr Q,χ, Q, χ denotes extremization with respect to the variables {Q, χ, Q, χ}.The function ξ is given by where √ χz is the random field that effectively represents the randomness of the problem introduced by y and A, and Dz = dz exp(−z 2 /2)/ √ 2π.Equation ( 75) is equivalent to the one-dimensional problem (12) with the correspondence σ 2 w → Q−1 and w → √ χz.The solution of (75), denoted by x * (z; Q, χ), is given by and is statistically equivalent to the solution of the original problem (1).Therefore, the expected density of nonzero components in the estimate is given by ρ where I(x) is the indicator function, which takes the value 1 when x is satisfied, otherwise it is zero.At the extremum of (73), the variables Q, χ, Q, χ are given by Note that the functional form of the parameters χ and Q does not depend on the regularization, but the values of χ and Q are regularization-dependent.

GDF for Gaussian random predictors
As discussed in [10], the expression of GDF for Gaussian random predictors using the replica method is given by for any regularization.The specific form of this expression depends on the regularization.We summarize the form of GDF for three regularizations.

1 penalty
For 1 penalty, the saddle point equation of χ is given by where ρ = erfc(λ/ √ 2 χ) and Substituting ( 83) and ( 81) into (82), GDF for 1 regularization is given by df Equation ( 85) is the number of parameters divided by M , and coincides with the wellknown results of GDF for LASSO [7], where AIC is the unbiased estimator of the prediction error.

SCAD
The saddle point equation of χ for SCAD regularization is given by where Here, γ S corresponds to the fraction of components of the regression coefficients that are in the transient region of SCAD (Sec.2.4.2).Equation (86) leads to the following expression for GDF: The coefficient of γ S in the second term corresponds to the rescaled variance of the estimate in the transient region.From (89), the difference between the prediction error and AIC is given by This result means that AIC underestimates the prediction error when SCAD regularization is used, and the transient region contributes to the increase in GDF.

MCP
The saddle point equation of χ for MCP is given by where and γ M is the fraction of the regression coefficient in the transient region of MCP (Sec.2.4.3).Equation (91) leads to where the coefficient of γ M corresponds to the rescaled variance of the estimate in the transient region.Hence, the difference between the prediction error and AIC is given by As with SCAD, the transient region contributes to the increase in GDF.

Unbiasedness of the AMP-based estimator
When the predictor matrix is homogeneous in the sense that the variance of each component is M −1 , the variable V ν for any ν in AMP converges to V , which is defined by [24,25] V (y, The Gaussian random predictor considered here corresponds to this case, and the estimator of GDF can be simplified as df The variable V (y, A) fluctuates depending on y and A, but its variance vanishes in the limit as N → ∞ and M → ∞ with N/M → α: which is a consequence of the density evolution equation [24,25].Further, the density evolution equation indicates that E y,A [V (y, A)] in AMP corresponds to χ in replica method.Therefore, from (82), (y, A)] = df (99) holds for sufficiently large system size.Equation (99) means that the estimator df (y, A) is asymptotically unbiased.
Figure 2 and Figure 3 show df at the AMP fixed point and GDF calculated by the replica method with α = 0.5 and σ 2 y = 1 for SCAD and MCP, respectively.In both figures, a = 3.7 and (a) and (b) represent the dependence on λ and ρ/α of df (1 ) .In AMP, we set N = 200, M = 100, and the result is the average over 1000 randomly generated y and A. The vertical dashed lines represent the de Almeida-Thouless (AT) condition [24].AMP does not converge on the left of the dashed line, and the diagonal dashed line of gradient 1 represents the bias term of AIC, df = ρ/α.For all parameter regions, the AMP result is consistent with the replica method, and the difference between df and the bias term of AIC is correctly described.

Application to real data
We applied our AMP-based estimators to the "Communities and Crime Unnormalized Data Set", which is available from the UCI machine learning repository.We used N = 70 of the 125 predictors in the original data, and redefined them to reduce the similarity (correlation) between the predictors.After preprocessing, the maximum value of the correlation between predictors was suppressed to be less than 0.685 (absolute value).However, this was insufficient to guarantee that assumption A2 would be satisfied.In addition, a part of the coefficients had larger values for which O(1/ √ M ) did not hold, and so assumption A1 was also violated.Among the 2215 response variables in the original data, we used M = 302 response variables for which the corresponding row vectors in the predictor matrix were not missing.The original data contained 18 kinds of response variables.Here, we demonstrate the model selection based on the proposed estimator using the response variable of "the number of murders per 100K population".
The simulation condition was as follows.The predictors and output variables were standardized.First, we computed the OLS estimate xOLS = A + y, where A + is the pseudo-inverse matrix of A. The signal x 0 estimated in this simulation was generated from the OLS estimate.We denote the set of indices of the largest K components of xOLS (in absolute value) as S K , and set and consider a synthetic model where each component of ξ ∈ R M is an i.i.d.Gaussian random variable with mean 0 and variance 1.
Figure 4 shows the λ-dependence of the prediction error and our proposed estimator ˆ (1) pre under the 1 penalty for K = 7 and K = 14.To evaluate the prediction error, we prepared 1000 test samples and imitated the expectation according to the generative process using the sample average.The value of the estimator was averaged over 1000 training samples to describe its typical behavior, and the averaged value was compared with the prediction error.In the case of the 1 penalty, the behavior of ˆ pre is similar to the prediction error, and it can select the model that minimizes the prediction error.
However, in the case of SCAD and MCP, the proposed estimators do not provide such a good match to the prediction error.Figure 5 shows the λ-dependence of the prediction error, the estimators of prediction error, and AIC at K = 7 for SCAD and MCP.The regularization parameter a is set to 3.7 in both cases.The prediction error attains a minimum at λ = 0.8 for SCAD and λ = 0.9 for MCP, but AIC cannot detect these minima.AIC tends to deviate from the prediction error as λ decreases, and reaches a minimum at smaller values of λ.In other words, the model selected based on AIC has more non-zero components than that selected by the minimization of the prediction error.In both cases, ˆ pre is slightly larger than AIC.The difference between AIC and ˆ (1) pre is interpreted as the contribution of the estimates in the transient region, referring to (90) and (95) for a Gaussian random predictor and response.However, under assumptions A1 and A2, the AMP-based estimator of GDF underestimates its true value.We introduce a heuristic method to modify the proposed estimator in the following section.

Correction of the estimator taking into account correlation between predictors
The estimator (70) is an increasing function of V µ (> 0), and hence the underestimation of the prediction error is caused by that of V µ .We correct the value of the rescaled variance {v i } included in V µ by considering the off-diagonal elements of A T A, after the convergence of AMP.
From (57) and (60), the variable v i corresponds to the variation of x i with respect to R i /Σ 2 .As shown in (61), R i /Σ 2 is the effective field of the single-body minimization problem of x i .The assumptions introduced for the derivation of AMP convert the original problem to the single-body problem, ignoring the correlation between the predictors.We go back to the original problem with respect to the support set at the fixed point of AMP, to re-calculate v and V taking into account the correlation between the predictors.The corrected v and V are denoted by ṽ and Ṽ , respectively.We denote the support indices at the fixed point of AMP by K, and denote the vector of the variables x j (j ∈ K) as x K ∈ R |K| , where |K| is the size of the support set K and the i-th component of K is denoted by K i .Furthermore, the submatrix of A consists of i(∈ K)-th columns and is denoted by A K .We consider the original minimization problem with respect to x K under an infinitesimal external field h ∈ R |K| as [7,26] The minimizer is the solution of the following equation where the i-th component of J ∈ R |K| , denoted by J i , is given by We correct v by implementing the variation with respect to the external field in the original problem.The infinitesimal external field helps the procedure as Therefore, the corrected variable Ṽ is expressed as Algorithm 1: AMP-based estimator of the prediction error(ˆ pre (y)) 1) Get the fixed point of AMP : Update variables according to (52)-(57) until they converge.Calculate the training error train (y) at the fixed point.2) Find K : Set the support of a at the fixed point of AMP as K, and construct a submatrix A K .3) Get x K : Solve (102) and set the solution as x K .4) Calculate ṽ : Calculate (104) for all components in the support.5) Get estimator of df : Calculate Ṽ as (105) and put it into (106), then obtain df (2) (y).6) Get estimator of the prediction error : According to (107), calculate the estimator of the prediction error ˆ pre (y).
Using Ṽµ , we define the estimator of GDF as df and the corresponding estimator of the prediction error as ˆ (2)  pre (y) ≡ train (y) + 2σ 2 y df

(y). (107)
Figure 6 is the pseudocode for the calculation of the AMP-based estimator of the prediction error.As a by-product of this correction, ṽi can be calculated using other algorithms that do not update v i , such as coordinate descent (CD) and Alternating Direction Method of Multipliers (ADMM).It is known that the convergence of AMP is slow when the predictor matrix is not i.i.d.Gaussian.In such cases, the combination of the AMP-based estimator with CD or ADMM is useful, as these techniques are less influenced by the properties of the predictor matrix than AMP.
We summarize the specific form of ṽi for SCAD and MCP.In the case of SCAD, the vector J is given by where Ψ ∈ R |K|×|K| and Φ ∈ R |K|×|K| are diagonal matrices whose components are given by The solution of (102) is given by and setting U In the case of MCP, the vector J is given by where Φ is the diagonal matrix whose components are We examined the performance of the modified AMP-based estimator ˆ pre in terms of the prediction error for the synthetic model (100) under SCAD and MCP.The results are shown in Figure 5.The estimator ˆ (2) pre captures the qualitative behavior of the prediction error better than AIC and ˆ (1) pre .In particular, ˆ pre attains a minimum at a value of λ that is close to that which minimizes the prediction error.Thus, it can effectively emulate the prediction error in the model selection process.Figures 7 and 8 show the models selected according to the minimum prediction error given by ˆ (2) pre and AIC for SCAD and MCP, respectively.The model selected based on the minimization of ˆ (2) pre is closer to that selected by the prediction error than the model given by AIC.These results indicate the advantage of the estimator ˆ (2) pre over AIC when it is necessary to select a model based on the minimization of the prediction error under nonconvex penalties.
Matlab codes for calculating the AMP-based estimator and the data used in this simulation are provided on our webpage.

Summary and Conclusion
We have proposed an AMP-based estimator of GDF and the prediction error for penalized linear regression problems.The proposed GDF estimator is given by the variable V µ in AMP for any regularization.The asymptotic property of the estimator df (1) for a Gaussian random predictor and response variable is consistent with the result derived by the replica method, and the asymptotic unbiasedness is mathematically guaranteed.The prediction error estimator using df (1) could be improved when applied to real data whose predictors are correlated.We corrected the estimator by considering pre and AIC, respectively.The minima of the prediction error, ˆ pre , and AIC are obtained with λ = 0.8, a = 3.7, λ = 1, a = 3.2, and λ = 0.45, a = 3.7, respectively.
the correlation between predictors, and thus constructed df (2) .We demonstrated the model selection process using the proposed estimators for nonconvex sparse penalties with real data, where AIC does not correspond to the unbiased estimator of the prediction error.The proposed estimator selects models that are close to the model that minimizes the prediction error.
Our proposed method can determine the regularization parameter, particularly for regularizations in which the unbiased estimator of the prediction error cannot be analytically derived.However, the correspondence between the proposed estimator and the prediction error is not mathematically guaranteed, except when the predictor matrix is i.i.d.Gaussian.The validity of our proposed method should be examined under several scenarios, especially for nonconvex penalties.For the real data used in this paper, qualitatively similar results are obtained for other response variables that the model selected based on our proposed estimators is close to the minimum of the prediction error.The generalization of our discussion to correlated predictors and rank-deficient predictors is required for more mathematically rigorous discussions [27].Furthermore, in practical usage, an efficient method of finding the minimum of the estimated prediction error should be developed.

Figure 2 .
Figure 2. GDF for SCAD calculated by replica method (solid line) and its estimator df(1) in AMP (circles) for α = 0.5, ρ = 0.1, and σ 2 y = 1: (a) as a function of λ and (b) as a function of ρ/α.The dashed vertical lines represent the AT condition, and the diagonal dashed line in (b) represents the bias term of AIC.

Figure 3 .
Figure 3. GDF for MCP calculated by the replica method (solid line) and its estimator df(1) in AMP (circles) for α = 0.5, ρ = 0.1, and σ 2 y = 1: (a) as a function of λ and (b) as a function of ρ/α.The dashed vertical lines represent the AT condition, and the diagonal dashed line in (b) represents the bias term of AIC.

Figure 4 .
Figure 4. λ-dependence of the prediction error and its estimator for "Communities and Crime Unnormalized Data Set" under 1 penalty with (a) K = 7 and (b) K = 14.

Figure 6 .
Figure 6.Pseudocode for the calculation of AMP-based estimator of the prediction error, ˆ

Figure 7 .
Figure 7.Comparison of models for SCAD.Models with the minimum prediction error are denoted by bars, and the circles in (a) and (b) represent the model given by the minimum of ˆ

Figure 8 .
Figure 8.Comparison of models for MCP.Models with the minimum prediction error are denoted by bars, and the circles in (a) and (b) represent the model given by the minimum of ˆ