# Prediction Errors for Penalized Regressions based on Generalized Approximate
Message Passing

Ayaka Sakata1,2,3 1Department of Statistical Modeling, The Institute of
Statistical Mathematics, Midori-cho, Tachikawa, Tokyo 190-8562, Japan
2Department of Statistical Science, The Graduate University for Advanced
Science (SOKENDAI), Hayama-cho, Kanagawa 240-0193, Japan  3JST PRESTO, 4-1-8
Honcho, Kawaguchi, Saitama 332-0012, Japan<EMAIL_ADDRESS>
###### Abstract

We discuss the prediction accuracy of assumed statistical models in terms of
prediction errors for the generalized linear model and penalized maximum
likelihood methods. We derive the forms of estimators for the prediction
errors, such as $C_{p}$ criterion, information criteria, and leave-one-out
cross validation (LOOCV) error, using the generalized approximate message
passing (GAMP) algorithm and replica method. These estimators coincide with
each other when the number of model parameters is sufficiently small; however,
there is a discrepancy between them in particular in the parameter region
where the number of model parameters is larger than the data dimension. In
this paper, we review the prediction errors and corresponding estimators, and
discuss their differences. In the framework of GAMP, we show that the
information criteria can be expressed by using the variance of the estimates.
Further, we demonstrate how to approach LOOCV error from the information
criteria by utilizing the expression provided by GAMP.

##  1 Introduction

Over the last decade, there has been a considerable increase in the use of
statistical methods and related machine learning techniques in various
research fields. This trend, which can be called the ‘data science revolution’
[1], is driven by qualitative and quantitative changes in data, increase in
computational power, and development of algorithms. Further, physics has also
been part of this trend; for example, the reconstruction of black hole images
in the event horizon telescope collaboration using sparse estimation is a
recent cooperative achievement between astrophysics and statistics [2].

Natural science is the quest for physical laws as the ground truth, and hence,
statistical inference in natural science aims at the inference of the truth.
However, all data do not necessarily possess the ground truth. After
industrialization, data from non-natural science fields have become
diversified [3], and the statistical inference for such non-nature data can be
considered practical science for the human society. For example, industrial
data are used for the reasonable determination, control, and creation of
values. In such cases, it is irrelevant if there are true physical laws, or if
our statistical models are different from these laws. The statistical model is
considered acceptable if practical operations based on this model do not
result in any problems.

Once the purpose of statistical inference becomes people-oriented, the
goodness of a statistical model becomes subjective, as the justification by
nature is impossible. Mathematical schemes are required for the construction
of statistical models to make statistical inference a scientific procedure.
The Akaike information criterion (AIC) [4] and subsequent cross validation
(CV) method [5] are schemes where the statistical model is objectively
selected using data among the set of models prepared subjectively. In these
methods, the goodness of the model is quantified by its predictive
performance, i.e. based on whether the model built by the given data can
describe the unknown data, which are expected to be observed in future. The
key idea for the mathematical formalism of the prediction is that the
predictive performance of the model is defined by the distributions rather
than the value of a certain statistic. This idea led to the geometrical
understanding of the maximum likelihood method and associated model selection
in terms of information theory as a projection of the real to the model space
[6].

Conventional statistics focuses on the construction of statistical models that
take into account empirical knowledge on the system, sometimes it requires
experts’ wisdom. In contrast to such approach, deep learning, which is a
central topic in the recent research on statistics and associated machine
learning, provides a framework to construct the generative model with less
manual tuning. The development of deep learning does not indicate the decline
of the traditional statistical inference. There have been attempts to
understand the properties of deep learning using statistical methodologies
such as understanding the double descent phenomena in an overparametrized
region using ridge regression [7], and characterizing the flat minima [8] in
deep neural networks using information criteria [9]. The new paradigm offered
by deep learning gives conventional statistics a role as a simplified system
for understanding the picture behind complicated inference problems. This
relationship will resonate with physicists, who are trying to extract
essentials of physical phenomena using idealized theoretical frameworks or
simplified experimental systems. In addition, there exists problems wherein it
is difficult to acquire large amounts of datasets because of the operating
cost and technical difficulties. In such cases, conventional statistical
modeling remains an important approach.

In this paper, we review model quantifications in terms of prediction
performance, and explain the physical consideration for statistical
methodologies. We introduce a statistical-physics-based approach for the
mathematical analysis of the prediction accuracy of the model under
consideration. The mathematical formalism for statistical inference has
similarities with statistical physics, in particular with the theory of spin-
glass systems, as the mathematical treatment of the randomness between spins’
interactions can be regarded as that of data [10, 11, 12]. Utilizing the
similarity, statistical physics-based methods have been applied to the
analyses of the inference problems as a mathematical tool [13, 14]. The
statistical physics methods are effective even in parameter regions where
asymptotic normality does not hold; this is key in classical statistics. This
paper indicates the complementarity between traditional statistical methods
and statistical physics methods.

We consider regression problems on the generalized linear model under penalty
functions. We focus on the estimator of the generalization gap, which is the
difference between the prediction and training errors [15, 16]. As estimators
for the prediction errors and corresponding generalization gaps, we analyze
the $C_{p}$ criterion [17], Widely Applicable Information Criterion (WAIC)
[18] or Takeuchi Information Criterion (TIC) [19] and CV method. We compare
these estimators in the same setting of the penalized maximum likelihood
method, i.e. the zero-temperature limit of the Bayesian inference in the
terminology of the statistical physics. The generalization gaps are quantified
by the generalized message passing (GAMP) algorithm. We show that the
generalization gap obtained using the $C_{p}$ criterion can be expressed by
the variance of the estimated parameters. This relationship resembles the
fluctuation-response relationship in statistical physics, in the sense that
the response of the error is proportional to the variance of the estimate at
the equilibrium [20]. We also show that both $C_{p}$ and WAIC or TIC are
insufficient for describing the parameter region where the number of parameter
is sufficiently large. We correct TIC towards CV by removing the assumption
introduced for the derivation of TIC. The resultant form reproduces the CV
error with high accuracy even in the non-asymptotic limit in statistics, and
one time estimation is sufficient under full data to obtain the form.

The remainder of the paper is organized as follows. Sec. 2 provides a brief
overview of the prediction errors and their estimators. In sec. 3, we explain
the problem setting discussed in this paper. We explain the GAMP algorithm in
sec. 4, and we apply the algorithm for the calculation of the prediction error
and corresponding generalization gap. In sec. 5, we present the difference
between the information criteria and CV error, and we explain how to approach
the CV error from TIC by utilizing the expression obtained using GAMP. We
derive the analytical form of the prediction error for random predictors where
each component is independently and identically distributed according to
Gaussian distribution in sec. 6. Finally, sec. 7 presents the summary and the
conclusion.

##  2 Prediction errors and their estimators

We focus on a parametric model $f(z|\bm{\theta})$ with parameter
$\bm{\theta}\in\bm{\Theta}\subset\mathbb{R}^{M}$ where $\bm{\Theta}$
represents a parameter space. The model parameter $\bm{\theta}$ is estimated
under the model $f(\cdot|\bm{\theta})$ and training data
$\bm{y}\in\mathbb{R}^{M}$ to describe the data generative process. Here, we
consider regression problems wherein the parameter $\bm{\theta}$ is described
by $\bm{\theta}(\bm{F},\bm{x})=\bm{F}\bm{x}/\sqrt{N}$, where
$\bm{F}\in\mathbb{R}^{M\times N}$ and $\bm{x}\in\mathbb{R}^{N}$ are the
predictor matrix and the regression coefficient to be estimated, respectively.
Therefore, parameter space $\Theta$ is specified by the distributions of the
predictors and regression coefficients. The data set consists of $\bm{y}$ and
$\bm{F}$, hence we denote ${\cal D}=\\{\bm{y},\bm{F}\\}$. In the following
discussion, we assume that each component of the training data
$y_{\mu}~{}(\mu=1,\cdots,M)$ is independently and identically generated
according to the predetermined rule. For mathematical convenience, we denote
the true rule as the distribution $q(y_{\mu}|\bm{F})$. Furthermore, we
consider a case in which the columns of the predictors are standardized such
that $\mathbb{E}_{\bm{F}}[\sum_{\mu=1}^{M}F_{\mu i}^{2}]=M$ for any $i$. Let
us prepare a set of candidate models ${\cal
M}=\\{f_{1}(z|\bm{\theta}(\bm{F},\widehat{\bm{x}}_{1})),\cdots,f_{m}(z|\bm{\theta}(\bm{F},\widehat{\bm{x}}_{m}))\\}$
for the true distribution, where $\widehat{\bm{x}}_{k}({\cal D})$ represents
the estimated regression coefficient under the $k$-th model $f_{k}$ using
training data ${\cal D}$. The model selection is the problem of adopting a
model based on a certain criterion from the set of candidate models.

We consider the penalized maximum likelihood method to construct the estimate
$\widehat{\bm{x}}({\cal D})$ as

$\displaystyle\widehat{\bm{x}}({\cal
D})=\mathop{\mathrm{argmax}}_{\bm{x}}\varphi({\cal D},\bm{x}),$ (2.1)

where

$\displaystyle\varphi({\cal D},\bm{x})=\ln
f(\bm{y}|\bm{\theta}(\bm{F},\bm{x}))-h(\bm{x})$ (2.2)

and $f(\bm{y}|\bm{\theta}(\bm{F},\bm{x}))$ and $h(\bm{x})$ are the likelihood
and penalty terms, respectively. In this paper, we focus on separable
likelihoods

$\displaystyle\ln f(\bm{y}|\bm{\theta}(\bm{F},\bm{x}))=\sum_{\mu=1}^{M}\ln
f(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x})),$ (2.3)

where $\mathbf{f}_{\mu}\in\mathbb{R}^{1\times N}$ represents the $\mu$-th row
vector of the predictor matrix, and
$\theta_{\mu}(\mathbf{f}_{\mu},\bm{x})=\frac{1}{\sqrt{N}}\mathbf{f}_{\mu}\bm{x}$.
Here, $h(\bm{x})$ contains the smoothing or regularization parameter
introduced to avoid overfitting [15, 21]. Models with varying regularization
parameter constitute the model space; the determination of the regularization
parameter corresponds to a model-selection problem.

In general, a model is constructed to represent the given training data sample
${\cal D}$ accurately. The accuracy of representing the training data can be
quantified by the training error defined by

$\displaystyle\mathrm{err}_{\mathrm{train}}({\cal
D})=-\frac{1}{M}\sum_{\mu=1}^{M}\ln
f(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\widehat{\bm{x}}({\cal D}))).$ (2.4)

The adequacy of the models with the parameter
$\bm{\theta}(\bm{F},\widehat{\bm{x}}({\cal D}))$ is quantified by the
prediction error, which indicates the goodness of the model for the
description of the test data not used for training. The high prediction error
implies the overfitting, and therefore, a model with a small prediction error
is required in practice.

###  2.1 Extra-sample and in-sample prediction errors

Figure 1: Definitions of extra-sample and in-sample prediction errors. A test
sample whose position is not in the training data is used for the evaluation
of the extra-sample prediction error. Test samples whose positions are the
same as that of the training sample are used for the in-sample prediction
error.

There are two definitions of the prediction error: extra-sample prediction
error and in-sample prediction error, which are considering different
generative processes of the test data used for evaluating the prediction error
[15]. Their definitions are given as follows.

######  Definition 1 (Extra-sample prediction error).

Let us denote ${\bf f}_{\mathrm{new}}\in\mathbb{R}^{1\times N}$ as a predictor
value that is not included in the training sample ${\cal D}$. The extra-sample
prediction error for the penalized maximum likelihood method is defined by

$\displaystyle\mathrm{err}_{\mathrm{pre}}^{(\mathrm{ex})}({\cal
D})=-\mathbb{E}_{z_{\mathrm{new}},{\bf f}_{\mathrm{new}}}\left[\ln
f\left(z_{\mathrm{new}}\Big{|}\frac{1}{\sqrt{N}}{\bf
f}_{\mathrm{new}}\widehat{\bm{x}}({\cal D})\right)\right],$ (2.5)

where $z_{\mathrm{new}}\in\mathbb{R}$ represents the test data associated with
the predictor ${\bf f}_{\mathrm{new}}$, and $\mathbb{E}_{z_{\mathrm{new}},{\bf
f}_{\mathrm{new}}}[\cdot]$ denotes the expectation with respect to
$z_{\mathrm{new}}$ and ${\bf f}_{\mathrm{new}}$.

######  Definition 2 (In-sample prediction error).

The in-sample prediction error for the penalized maximum likelihood method is
given by

$\displaystyle\mathrm{err}_{\mathrm{pre}}^{(\mathrm{in})}({\cal
D})=-\frac{1}{M}\sum_{\mu=1}^{M}\mathbb{E}_{z_{\mu}\sim q(z|{\bf
f}_{\mu})}[\ln{f(z_{\mu}|\theta_{\mu}({\bf f}_{\mu},\widehat{\bm{x}}({\cal
D})))}],$ (2.6)

where $z_{\mu}$ represents the new response variable (test data) observed at
$\mu$-th training point.

Fig.1 schematically shows the difference between the test data for the extra-
sample and in-sample prediction errors. The terms ‘extra’ and ‘in’ indicate
whether the position of the test data overlaps the training set. The test data
is on the location not covered by training samples in the case of the extra-
sample error, hence the extra-sample error describes the prediction error for
varying the predictors. The in-sample prediction error considers the test data
generated by repeated observation on the same position. The in-sample
prediction error sometimes does not match the actual problem settings;
however, it yields a convenient comparison method between models and leads to
an understanding of the model complexity as explained later. In the parameter
region where the conventional statistics focused on, the tendencies of the two
prediction errors are considered to be similar with each other [15].

The exact values of both prediction errors cannot be obtained because we do
not know the true generative process of the data. Instead, estimators for the
prediction errors are used for evaluating the prediction accuracy.

###  2.2 Cross-validation error

The cross-validation error is an estimator for the extra-sample prediction
error [5].

######  Definition 3 (Leave-one-out cross validation error).

Let us denote a leave-one-out (LOO) sample ${\cal D}_{\setminus\mu}$ where the
$\mu$-th data sample ${\cal D}_{\mu}\equiv\\{y_{\mu},\mathbf{f}_{\mu}\\}$ is
discarded from the data set ${\cal D}$. The estimate under the LOO sample,
which is denoted by $\widehat{\bm{x}}({\cal D}_{\setminus\mu})$, is defined as

$\displaystyle\widehat{\bm{x}}({\cal
D}_{\setminus\mu})=\mathop{\mathrm{argmin}}_{\bm{x}}\sum_{\nu\in\bm{M}\setminus\mu}\ln
f(y_{\nu}|\theta_{\nu}(\mathbf{f}_{\nu},\widehat{\bm{x}}({\cal
D}_{\nu})))-h(\bm{x}),$ (2.7)

where $\bm{M}=\\{1,2,\ldots,M\\}$ and $\bm{M}\setminus\mu$ denotes $\bm{M}$
except $\mu$. The leave-one-out cross validation (LOOCV) error is defined as

$\displaystyle\mathrm{err}_{\mathrm{LOOCV}}({\cal
D})=-\frac{1}{M}\sum_{\mu=1}^{M}\ln
f\left(y_{\mu}\Big{|}\frac{1}{\sqrt{N}}{\bf f}_{\mu}\widehat{\bm{x}}({\cal
D}_{\setminus\mu})\right).$ (2.8)

Intuitively, for a sufficiently large $M$, the empirical mean with respect to
the data sample is assumed to converge to the expectation value with respect
to the true distribution. Therefore, the LOOCV error is considered an adequate
estimator of the extra-sample prediction error.

####  2.2.1 Special form of LOOCV error for the linear estimation rule

For calculating the LOOCV error, we need to compute the estimate $M$ times
under $M$-LOO-samples ${\cal D}_{\setminus\mu}$ for $\mu\in\bm{M}$. There is a
long history of analytical and numerical attempts to reduce the computational
effort involved in evaluating the LOOCV error. Here, we focus our discussion
on the linear estimation rule.

######  Definition 4 (Linear estimation rule).

In the linear estimation rule, the output $\bm{y}$ is estimated as
$\widehat{\bm{y}}=\bm{H}\bm{y}$ on the assumed model with an arbitrary matrix
$\bm{H}\in\mathbb{R}^{M\times M}$.

The ordinary least square (OLS) method is one of the One of the linear
estimation rules, wherein the regression coefficient is given as
$\widehat{\bm{x}}({\cal
D})=\frac{1}{\sqrt{N}}(\frac{1}{N}\bm{F}^{\top}\bm{F})^{-1}\bm{F}^{\top}\bm{y}$
at $M>N$, then $\widehat{\bm{y}}=\bm{H}\bm{y}$ with
$\bm{H}=\frac{1}{N}\bm{F}(\bm{F}^{\top}\bm{F})^{-1}\bm{F}^{\top}$. Here, the
estimated output corresponds to the projection of the actual output onto the
column space of the predictor matrix, and the matrix $\bm{H}$ for the OLS case
is known as the hat matrix [22].

In the linear estimation methods, the computation of the LOOCV error is
reduced to a one-time calculation of the estimate under full training data
${\cal D}$ and is expressed using the matrix $\bm{H}$. The expression for the
Gaussian likelihood shown in the following theorem is known as the predicted
residual sum of squares (PRESS) statistics [23], and was derived by [24, 25].

######  Theorem 1.

For the linear estimation rule where $\widehat{\bm{y}}=\bm{H}\bm{y}$ holds,
the LOOCV error under the Gaussian likelihood is given as

$\displaystyle\mathrm{err}_{\mathrm{LOOCV}}({\cal
D})=\frac{1}{2M}\sum_{\mu=1}^{M}\left(\frac{y_{\mu}-\widehat{y}_{\mu}({\cal
D})}{1-H_{\mu\mu}}\right)^{2}.$ (2.9)

###### Proof.

We prepare a data set ${\cal
D}_{\setminus\mu^{+}}=\\{\bm{y}_{\setminus\mu^{+}},\bm{F}\\}$, where
$\bm{y}_{\setminus\mu^{+}}=\\{y_{1},\cdots,y_{\mu-1},\widehat{y}_{\mu}({\cal
D}_{\setminus\mu}),y_{\mu+1},\cdots,y_{M}\\}$ and $\widehat{y}_{\mu}({\cal
D}_{\setminus\mu})=\mathbf{f}_{\mu}\widehat{\bm{x}}({\cal
D}_{\setminus\mu})/\sqrt{N}$, as shown in Fig.2. The estimate under the data
${\cal D}_{\setminus\mu^{+}}$ is defined by

$\displaystyle\widehat{\bm{x}}({\cal D}_{\setminus\mu^{+}})$
$\displaystyle=\mathop{\mathrm{argmax}}_{\bm{x}}\Big{\\{}\sum_{\nu\in\bm{M}\setminus\mu}\ln
f(y_{\nu}|\theta_{\nu}(\mathbf{f}_{\nu},\bm{x}))$ (2.10) $\displaystyle\hskip
85.35826pt+\ln f(\widehat{y}_{\mu}({\cal
D}_{\setminus\mu})|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x}))-h(\bm{x})\Big{\\}}.$

We can show that $\widehat{\bm{x}}({\cal
D}_{\setminus\mu})=\widehat{\bm{x}}({\cal D}_{\setminus\mu^{+}})$ holds for
any $\mu\in\bm{M}$. By definition, the following inequality holds for any
$\bm{x}$; i.e.

$\displaystyle\sum_{\nu\in\bm{M}\setminus\mu}\ln
f(y_{\nu}|\theta_{\nu}(\widehat{\bm{x}}(\mathbf{f}_{\nu},{\cal
D}_{\setminus\mu}))$ $\displaystyle-h(\widehat{\bm{x}}({\cal
D}_{\setminus\mu}))$ (2.11)
$\displaystyle\geq\sum_{\nu\in\bm{M}\setminus\mu}\\!\\!\ln
f(y_{\nu}|\theta_{\nu}(\mathbf{f}_{\nu},\bm{x}))-h(\bm{x}).$

Further, it is natural to consider that the following relationship holds for
any $\theta$

$\displaystyle\ln f(y_{\mu}=\theta_{\mu}|\theta_{\mu})\geq\ln
f(y_{\mu}=\theta_{\mu}|\theta).$ (2.12)

Combining (2.11) with (2.12), we obtain

$\displaystyle\varphi({\cal D}_{\setminus\mu^{+}},\bm{x})$
$\displaystyle=\sum_{\nu\in\bm{M}\setminus\mu}\ln
f(y_{\nu}|\theta_{\nu}(\mathbf{f}_{\nu},\bm{x}))+\ln
f(\theta_{\mu}(\mathbf{f}_{\mu},\widehat{\bm{x}}({\cal
D}_{\setminus\mu}))|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x}))-h(\bm{x})$ (2.13)
$\displaystyle\leq\sum_{\nu\in\bm{M}\setminus\mu}\ln
f(y_{\nu}|\theta_{\nu}(\mathbf{f}_{\nu},\widehat{\bm{x}}({\cal
D}_{\setminus\mu})))-h(\widehat{\bm{x}}({\cal D}_{\setminus\mu}))$
$\displaystyle\hskip 85.35826pt+\ln
f(\theta_{\mu}(\mathbf{f}_{\mu},\widehat{\bm{x}}({\cal
D}_{\setminus\mu}))|\theta_{\mu}(\mathbf{f}_{\mu},\widehat{\bm{x}}({\cal
D}_{\setminus\mu})))$ $\displaystyle=\varphi({\cal
D}_{\setminus\mu^{+}},\widehat{\bm{x}}({\cal D}_{\setminus\mu}))$

hence $\widehat{\bm{x}}({\cal D}_{\setminus\mu^{+}})=\widehat{\bm{x}}({\cal
D}_{\setminus\mu})$ is one of the solutions of (2.10), and for the linear
estimator, the following relationship holds:

$\displaystyle\frac{1}{\sqrt{N}}\sum_{i=1}^{N}F_{\nu i}\widehat{x}_{i}({\cal
D}_{\setminus\mu^{+}})=\sum_{\eta\in\bm{M}\setminus\mu}H_{\nu\eta}y_{\eta}+H_{\nu\mu}\widehat{y}_{\mu}({\cal
D}_{\setminus\mu}).$ (2.14)

Meanwhile, by utilizing the relationship $\widehat{\bm{x}}({\cal
D}_{\setminus\mu^{+}})=\widehat{\bm{x}}({\cal D}_{\setminus\mu})$, we obtain

$\displaystyle\widehat{y}_{\mu}({\cal D}_{\setminus\mu})$
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}F_{\mu i}\widehat{x}_{i}({\cal
D}_{\setminus\mu})=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}F_{\mu
i}\widehat{x}_{i}({\cal D}_{\setminus\mu^{+}})$ (2.15)
$\displaystyle=H_{\mu\mu}\widehat{y}_{\mu}({\cal
D}_{\setminus\mu})+\sum_{\eta\in\bm{M}\setminus\mu}H_{\mu\eta}y_{\eta}$
$\displaystyle=H_{\mu\mu}\widehat{y}_{\mu}({\cal
D}_{\setminus\mu})+\widehat{y}_{\mu}({\cal D})-H_{\mu\mu}y_{\mu}.$

For the Gaussian likelihood, the LOOCV error corresponds to
$\frac{1}{2M}\sum_{\mu=1}^{M}(y_{\mu}-\widehat{y}_{\mu}({\cal
D}_{\setminus\mu}))$, and using (2.15), we obtain (2.9). ∎

Eq.(2.15) is valid for the linear estimation rule, wherein one can compute the
LOOCV error by substituting $\hat{y}_{\mu}({\cal D}_{\setminus\mu})$ with the
likelihood under consideration.

Figure 2: Schematic of data samples ${\cal D}_{\setminus\mu}$ and ${\cal
D}_{\setminus\mu^{+}}$. The data set ${\cal D}_{\setminus\mu}$ comprises the
input–output relationship except the $\mu$-th sample denoted by $\varnothing$.
The point $\square$ of ${\cal D}_{\setminus\mu^{+}}$ is at
$\\{\mathbf{f}_{\mu},\theta_{\mu}(\widehat{\bm{x}}({\cal
D}_{\setminus\mu}))\\}$, where $\widehat{\bm{x}}({\cal D}_{\setminus\mu})$
denotes the estimate under ${\cal D}_{\setminus\mu}$, and ${\cal
D}_{\setminus\mu^{+}}=\\{{\cal D}_{\setminus\mu}\cup\square\\}$.

####  2.2.2 Sherman–Morrison Formula for the LOOCV error

The same expression as (2.9) can be obtained using the Sherman–Morrison
formula [26] for linear estimators, as shown in [27]. We define the predictor
matrix without the $\nu$-th row vector as $\bm{F}_{\setminus\nu}$. The
Sherman–Morrison formula is given as

$\displaystyle\left(\frac{1}{N}{\bm{F}_{\setminus\nu}}^{\top}\bm{F}_{\setminus\nu}\right)^{-1}$
$\displaystyle\\!\\!\\!=\\!\left(\frac{1}{N}\bm{F}^{\top}\bm{F}\right)^{-1}\\!\\!\\!+\\!\frac{\displaystyle\left(\frac{1}{N}\bm{F}^{\top}\bm{F}\right)^{-1}\\!\\!\\!\frac{1}{N}{\bf
f}_{\nu}^{\top}{\bf
f}_{\nu}\left(\frac{1}{N}\bm{F}^{\top}\bm{F}\right)^{-1}}{1-\displaystyle\frac{1}{N}{\bf
f}_{\nu}\left(\frac{1}{N}\bm{F}^{\top}\bm{F}\right)^{-1}{\bf
f}_{\nu}^{\top}}.$ (2.16)

The estimates under the LOO sample ${\cal D}_{\setminus\mu}$ is given by
$\widehat{\bm{x}}({\cal
D}_{\setminus\mu})=\frac{1}{\sqrt{N}}(\frac{1}{N}\bm{F}_{\setminus\mu}^{\top}\bm{F}_{\setminus\mu})^{-1}\bm{F}_{\setminus\mu}^{\top}$,
and hence, $\widehat{\bm{y}}({\cal
D}_{\setminus\mu})=\bm{H}^{\setminus\mu}\bm{y}$, where
$\bm{H}^{\setminus\mu}=\frac{1}{N}\bm{F}(\frac{1}{N}\bm{F}_{\setminus\mu}^{\top}\bm{F}_{\setminus\mu})^{-1}\bm{F}_{\setminus\mu}^{\top}$.
Multiplying $\bm{F}/\sqrt{N}$ and $\bm{F}_{\setminus\mu}^{\top}/\sqrt{N}$ from
the left and right of both sides of (2.16), respectively, we get

$\displaystyle\widehat{y}_{\mu}({\cal
D}_{\setminus\mu})=\widehat{y}_{\mu}({\cal
D})-H_{\mu\mu}y_{\mu}+\frac{H_{\mu\mu}(\widehat{y}_{\mu}({\cal
D})-H_{\mu\mu}y_{\mu})}{1-H_{\mu\mu}},$ (2.17)

Then, transforming (2.17), we obtain the same expression as (2.9).

###  2.3  Information criteria as estimators of the in-sample prediction error

Information criteria are known as the estimators of the in-sample prediction
error. These criteria were originally introduced to evaluate the quality of
the statistical model constructed with the maximum likelihood method using the
Kullback–Leibler (KL) divergence defined by

$\displaystyle{\rm KL}(q:f)=\mathbb{E}_{\bm{z}\sim
q(\bm{z})}\left[\ln{q(\bm{z})}\right]-\mathbb{E}_{\bm{z}\sim
q(\bm{z})}[\ln{f(\bm{z}|\bm{\theta}(\bm{F},\widehat{\bm{x}}_{\rm ML}({\cal
D}))}],$ (2.18)

where $\widehat{\bm{x}}_{\rm ML}({\cal D})$ represents the maximum likelihood
estimator under the training sample ${\cal D}$. A statistical model with
smaller KL divergence to the true distribution is expected to be appropriate
for the expression of given data. The second term of KL divergence (2.18)
corresponds to the in-sample prediction error under the ML estimator, and the
first term is a constant. Therefore, the minimization of KL divergence
corresponds to that of the in-sample prediction error.

Takeuchi information criterion (TIC) is an estimator of the in-sample
prediction error and in particular under the maximum likelihood method, TIC
corresponds to an unbiased estimator of the in-sample prediction error [19].

######  Definition 5 (Takeuchi information criterion).

The TIC is defined as

$\displaystyle\mathrm{TIC}({\cal
D})=2\left\\{\mathrm{err}_{\mathrm{train}}({\cal
D})+\frac{1}{M}\mathrm{Tr}\left(\widehat{\cal I}(\widehat{\bm{x}}({\cal
D}))\widehat{\cal J}^{-1}(\widehat{\bm{x}}({\cal D}))\right)\right\\},$ (2.19)

where $\widehat{\cal I}$ and $\widehat{\cal J}$ represent the estimator of the
Fisher information matrix $\overline{\cal I}$ and the expected Hessian
$\overline{\cal J}$ given as

$\displaystyle\overline{\cal I}(\widehat{\bm{x}})$
$\displaystyle=\frac{1}{M}\mathbb{E}_{\bm{z}\sim
q(\bm{z})}\left[\frac{\partial\ln
f(\bm{z}|\bm{\theta}(\bm{F},\bm{x}))}{\partial\bm{x}}\frac{\partial\ln
f(\bm{z}|\bm{\theta}(\bm{F},\bm{x}))}{\partial\bm{x}^{\top}}\Big{|}_{\bm{x}=\widehat{\bm{x}}}\right]$
(2.20) $\displaystyle\overline{\cal J}(\widehat{\bm{x}})$
$\displaystyle=-\frac{1}{M}\mathbb{E}_{\bm{z}\sim
q(\bm{z})}\left[\frac{\partial^{2}\ln
f(\bm{z}|\bm{\theta}(\bm{F},\bm{x}))}{\partial\bm{x}\partial\bm{x}^{\top}}\Big{|}_{\bm{x}=\widehat{\bm{x}}}\right],$
(2.21)

respectively.

The factor 2 in the definition of TIC comes from the convention that TIC for
the Gaussian likelihood corresponds to the estimator for the squared error. As
the estimators of Fisher information matrix and the expected Hessian, sample-
mean of the (non-centered) covariance of gradient ${\cal I}(\bm{x})$ and that
of the negative Hessian ${\cal J}(\bm{x})$ are used, where

$\displaystyle{\cal I}_{ij}(\widehat{\bm{x}})$
$\displaystyle=\frac{1}{M}\sum_{\mu=1}^{M}\\!\frac{\partial}{\partial
x_{i}}\\!\ln\\!f(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}}\\!\frac{\partial}{\partial
x_{j}}\\!\ln\\!f(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}},$
(2.22) $\displaystyle{\cal J}_{ij}(\widehat{\bm{x}})$
$\displaystyle=-\frac{1}{M}\frac{\partial^{2}}{\partial x_{i}\partial
x_{j}}\sum_{\mu=1}^{M}f(y_{\mu}|\theta_{\mu}(\bm{F},\bm{x})\Big{|}_{\bm{x}=\widehat{\bm{x}}}.$
(2.23)

with the expectation that they converge to $\overline{\cal I}$ and
$\overline{\cal J}$ from the law of large numbers.

For the special case where the model space contains the true model, where a
regression coefficient $\bm{x}^{*}$ exists such that
$q(\bm{z})=f(\bm{z}|\frac{1}{\sqrt{N}}\bm{F}\bm{x}^{*})$, it is shown that TIC
is reduced to Akaike’s information criterion (AIC) [4], which is defined as

$\displaystyle\mathrm{AIC}({\cal
D})=2\left(\mathrm{err}_{\mathrm{train}}({\cal D})+\frac{N}{M}\right).$ (2.24)

This is the consequence of $\overline{\cal I}\overline{\cal
J}^{-1}=\bm{I}_{N}$ that holds in the special case [28], where $\bm{I}_{N}$
represents the $N$-dimensional identity matrix.

The difference between the prediction and training error is termed the
generalization gap. Based on the correspondence between the generalization gap
obtained by AIC and the dimension of the parameters, the generalization gap is
considered for expressing the model complexity. The model can express various
distributions with an increase in the complexity of the model, and therefore,
there is a decrease in the training error. However, the models with large
parameters are prone to overfitting, which hampers the prediction of unknown
data, because these models can fit the randomness in given training data. The
form of AIC is given by the training error penalized by the number of model
parameters. Therefore, model selection based on the minimization of AIC
represents a tradeoff between fitting accuracy to training data and model
complexity.

We note that the meaning of AIC as an unbiased estimator of the in-sample
prediction error holds in the case where the the maximum likelihood estimator
is considered, and the model candidates include true distribution. In the case
where the latter condition does not hold, TIC is theoretically adequate for
evaluating the in-sample prediction error. This does not mean, however, that
the AIC should not be used, although theoretical justification is difficult.
In fact, AIC is, at times, superior to TIC in terms of low computational cost
and the small variance. To evaluate TIC, we need to compute the inverse of
Hessian, whose computational cost is $O(N^{3})$, and the fluctuation of the
generalization gap of TIC tends to be large, in particular when eigenvalues
close to zero exist [29].

###  2.4 Asymptotic equivalence between the LOOCV error and information
criterion

The LOOCV error has an asymptotic equivalence to the information criterion, as
proved in [30] with the following assumptions:

A1

The differences in log-likelihoods under full sample and LOO sample can be
expressed by the first order of the difference between the corresponding
estimates, under the consideration that the difference $\widehat{x}_{i}({\cal
D}_{\setminus\mu})-\widehat{x}_{i}({\cal D})$ is sufficiently small for any
$i\in\bm{N}$ and $\mu\in\bm{M}$,

Here, we set $\bm{N}=\\{1,\ldots,N\\}$. Under assumption A1, we expand the
difference in log-likelihoods for ${\cal D}$ and ${\cal D}_{\setminus\mu}$ up
to the first order of $\widehat{\bm{x}}({\cal D})-\widehat{\bm{x}}({\cal
D}_{\setminus\mu})$ as

$\displaystyle\ln f(y_{\mu}$
$\displaystyle|\theta_{\mu}(\mathbf{f}_{\mu},\widehat{\bm{x}}({\cal D})))-\ln
f(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\widehat{\bm{x}}({\cal
D}_{\setminus\mu})))$ (2.25)
$\displaystyle=\sum_{i=1}^{N}\frac{\partial}{\partial x_{i}}\ln
f(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\mu})}(\widehat{x}_{i}({\cal D})-\widehat{x}_{i}({\cal
D}_{\setminus\mu})).$

Following the procedure explained in A under assumption A1, we obtain the
relationship

$\displaystyle\sum_{j=1}^{N}{\cal J}_{ij}^{\varphi}(\widehat{\bm{x}}({\cal
D}_{\setminus\nu}))(\widehat{x}_{i}({\cal D})-\widehat{x}_{i}({\cal
D}_{\setminus\nu}))=\frac{1}{M}\frac{\partial}{\partial x_{i}}\ln
f(y_{\nu}|\theta_{\nu}(\mathbf{f}_{\nu},\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})},$ (2.26)

where

$\displaystyle{\cal J}_{ij}^{\varphi}(\widehat{\bm{x}})$
$\displaystyle=-\frac{1}{M}\frac{\partial^{2}}{\partial x_{i}\partial
x_{j}}\varphi({\cal D},\bm{x})\Big{|}_{\bm{x}=\widehat{\bm{x}}}.$ (2.27)

Here, we introduce the assumption:

A2

${\cal J}^{\varphi}_{ij}(\widehat{\bm{x}}({\cal D}_{\setminus\nu}))={\cal
J}^{\varphi}_{ij}(\widehat{\bm{x}}({\cal D}))$ holds for any $\nu\in\bm{M}$
and $i,j\in\bm{N}$.

Under assumption A2, we obtain

$\displaystyle\widehat{\bm{x}}({\cal D})-\widehat{\bm{x}}({\cal
D}_{\setminus\nu})=\frac{1}{M}({\cal J}^{\varphi}(\widehat{\bm{x}}({\cal
D})))^{-1}\frac{\partial}{\partial\bm{x}}\ln
f(y_{\nu}|\theta_{\nu}(\mathbf{f}_{\nu},\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})},$ (2.28)

where we implicitly assume that the Hessian (2.27) is invertible.

We introduce another assumption as follows:

A3

${\cal I}^{\setminus\cdot}(\\{\widehat{\bm{x}}({\cal
D}_{\setminus\nu})\\})={\cal I}(\widehat{\bm{x}}({\cal D}))$ holds, where

$\displaystyle{\cal I}_{ij}^{\setminus\cdot}(\\{\widehat{\bm{x}}({\cal
D}_{\setminus\nu})\\})=\frac{1}{M}\sum_{\mu=1}^{M}\frac{\partial}{\partial
x_{i}}\ln
f(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\mu})}$ $\displaystyle\hskip
113.81102pt\times\frac{\partial}{\partial x_{j}}\ln
f(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\mu})}.$ (2.29)

Under assumption A3, substituting (2.28) into (2.25) yields

$\displaystyle\mathrm{err}_{\mathrm{LOOCV}}({\cal
D})=\mathrm{err}_{\mathrm{train}}({\cal D})+\frac{1}{M}\mathrm{Tr}\left({\cal
I}(\widehat{\bm{x}}({\cal D}))({\cal J}^{\varphi}(\widehat{\bm{x}}({\cal
D})))^{-1}\right),$ (2.30)

where $\mathrm{Tr}(\cdot)$ denotes the summation of the diagonal components.
The expression (2.30) is equivalent to TIC (2.19) for the ML estimator where
${\cal J}^{\varphi}={\cal J}$.

The equivalence between LOOCV error and TIC under the assumptions A1-A3 is
called asymptotic equivalence in the terminology of statistics. We mention the
case where these assumptions are valid. The assumption A1 holds when $N$ is
sufficiently small. At $N\to\infty$, the higher order terms with respect to
the difference $\widehat{x}_{i}({\cal D})-\widehat{x}_{i}({\cal
D}_{\setminus\nu})$ have a contribution onto the difference in the log-
likelihoods, as shown in A.1, and the correction of (2.26) is required. The
assumptions A2 and A3 holds when $M$ is sufficiently large compared with $N$.
In addition, for the derivation of (2.28) under the assumption A2, the matrix
${\cal J}^{\varphi}$ needs to be invertible. This condition is violated in the
singular models such as mixture models, matrix factorization, and deep neural
networks. For such problems, TIC cannot be applied in the original form, but
it is suggested that TIC discarding zero-eigenmodes well characterize the
generalization property of the deep neural network [9].

###  2.5 Mallows’ $C_{p}$ and model complexity in linear models

The model selection criteria explained thus far are for models constructed
using the maximum likelihood method. A general expression of the estimator for
the in-sample generalization error is derived from another model selection
criterion, that is, Mallows’ $C_{p}$ [17]. We consider a linear estimation
rule where $\widehat{\bm{y}}=\bm{H}\bm{y}$ holds with a matrix $\bm{H}$, and
Gaussian likelihood with variance $\sigma^{2}$.

######  Definition 6 ($C_{p}$ criterion and degrees of freedom).

$C_{p}$ criterion is an estimator for the in-sample prediction error defined
by

$\displaystyle{c_{p}}(\bm{y})=2\sigma^{2}{\rm err}_{\rm
train}(\bm{y})+2{\sigma^{2}}\mathrm{df}.$ (2.31)

The term df is called degrees of freedom (DF) [31, 15] defined by

$\displaystyle\mathrm{df}=\frac{1}{M}\mathrm{Tr}(\bm{H}).$ (2.32)

Under the known variance assumption, $C_{p}$ criterion divided by $\sigma^{2}$
is reduced to AIC. To demonstrate this, let us apply a singular value
decomposition to $\bm{F}$ as $\bm{F}=\bm{U}\bm{D}\bm{V}^{\top}$, where
$\bm{U}\in\mathbb{R}^{M\times N}$ and $\bm{V}\in\mathbb{R}^{N\times N}$
represent the left and right singular vectors, respectively, and $\bm{D}$
represents the $N\times N$ dimensional diagonal matrix whose $(i,i)$-component
corresponds to the $i$-th singular value $d_{i}$. We obtain
$\mathrm{Tr}(\bm{H})=\sum_{i=1}^{N}d_{i}^{2}/d_{i}^{2}=N$, which indicates the
equivalence of the $C_{p}$ criterion and AIC.

$C_{p}$ can be extended to general estimators $\widehat{\bm{x}}({\cal D})$ in
regression problems. The following theorem shown in [16] is the basis for the
extension:

######  Theorem 2.

Let us consider the case where the following conditions hold:

  * •

The data generative process is given by Gaussian distribution $\bm{y}\sim{\cal
N}(\bm{\mu},\sigma^{2}\bm{I}_{M})$, where $\bm{I}_{M}$ represents the
$M$-dimensional identity matrix.

  * •

The likelihood is given by Gaussian distribution as

$\displaystyle
f(\bm{y}|\bm{\theta})\propto\exp\left(-\frac{||\bm{y}-\bm{\theta}||_{2}^{2}}{2\sigma^{2}}\right)$
(2.33)

  * •

The variance $\sigma^{2}$ is known.

The following estimator

$\displaystyle\widehat{\mathrm{Err}}({\cal
D})=2\sigma^{2}\mathrm{err}_{\mathrm{train}}({\cal
D})+\frac{2}{M}\widehat{\mathrm{Cov}}(\bm{y},\widehat{\bm{y}}({\cal D}))$
(2.34)

is an unbiased estimator of the in-sample prediction error with the
multiplying factor $2\sigma^{2}$ when
$\widehat{\mathrm{Cov}}(\bm{y},\widehat{\bm{y}}({\cal D}))$ is an unbiased
estimator of the covariance
$\mathrm{Cov}_{\bm{y}}(\bm{y},\widehat{\bm{y}}({\cal
D}))=\mathbb{E}_{\bm{y}\sim q(\bm{y})}[\bm{y}^{\top}\widehat{\bm{y}}({\cal
D})]-\mathbb{E}_{\bm{y}\sim q(\bm{y})}[\bm{y}]^{\top}\mathbb{E}_{\bm{y}\sim
q(\bm{y})}[\widehat{\bm{y}}({\cal D})]$.

###### Proof.

In the conditions considered in Theorem 2, the in-sample prediction and
training errors are respectively given by

$\displaystyle\mathrm{err}_{\mathrm{pre}}^{(\mathrm{in})}({\cal D})$
$\displaystyle=\frac{1}{2\sigma^{2}M}\mathbb{E}_{\bm{z}\sim
q(\bm{z})}[||\bm{z}-\bm{\theta}(\bm{F},\widehat{\bm{x}}({\cal
D}))||_{2}^{2}],$ (2.35) $\displaystyle\mathrm{err}_{\mathrm{train}}({\cal
D})$
$\displaystyle=\frac{1}{2\sigma^{2}M}||\bm{y}-\bm{\theta}(\bm{F},\widehat{\bm{x}}({\cal
D}))||_{2}^{2}.$ (2.36)

Using the relationship $\mathbb{E}_{\bm{z}\sim
q(\bm{z})}[\bm{z}]=\mathbb{E}_{\bm{y}\sim q(\bm{y})}[\bm{y}]$, the
generalization gap is given by the form known as covariance penalty [16]

$\displaystyle\Delta^{(\mathrm{in})}=\frac{1}{\sigma^{2}M}\mathrm{Cov}_{\bm{y}}\left(\bm{y},\widehat{\bm{y}}({\cal
D})\right).$ (2.37)

Therefore, we obtain

$\displaystyle 2\sigma^{2}\mathrm{err}_{\mathrm{pre}}^{(\mathrm{in})}({\cal
D})=2\sigma^{2}\mathrm{err}_{\mathrm{train}}({\cal
D})+\frac{2}{M}\mathrm{Cov}_{\bm{y}}\left(\bm{y},\widehat{\bm{y}}({\cal
D})\right),$ (2.38)

and this expression indicates that (2.34) is an unbiased estimator of the in-
sample prediction error. ∎

For constructing the estimator for the covariance penalty, the following
identity that holds when $\bm{y}\sim{\cal N}(\bm{\mu},\sigma^{2}\bm{I}_{M})$
is convenient, known as Stein’s lemma [32];

$\displaystyle\frac{1}{\sigma^{2}M}\mathrm{Cov}_{\bm{y}}\left(\bm{y},\widehat{\bm{y}}({\cal
D})\right)=\frac{1}{M}\sum_{\mu=1}^{M}\mathbb{E}_{\bm{y}}\left[\frac{\partial\widehat{y}_{\mu}({\cal
D})}{\partial y_{\mu}}\right].$ (2.39)

For the linear estimation rule, we obtain $\partial\widehat{y}_{\mu}({\cal
D})/\partial y_{\mu}=H_{\mu\mu}$; therefore, DF is an unbiased estimator of
the covariance penalty. We can obtain the estimator of the generalization gap
even when the analytical form of $\widehat{\bm{y}}$ cannot be derived using
the relationship (2.39). Stein’s unbiased risk estimate (SURE) of the mean
squared error [32] is an estimator of
$2\sigma^{2}\mathrm{err}_{\mathrm{pre}}^{(\mathrm{in})}$ for the Gaussian
likelihood applicable to general estimators of $\widehat{\bm{y}}$.

######  Definition 7 (SURE and generalized degrees of freedom).

SURE is an estimator of the in-sample prediction error defined as

$\displaystyle\mathrm{SURE}({\cal
D})=2\sigma^{2}\mathrm{err}_{\mathrm{train}}({\cal
D})+2\sigma^{2}\mathrm{gdf}({\cal D}),$ (2.40)

where $\mathrm{gdf}(\bm{y})$ is termed the generalized degree of freedom (GDF)
[33], which is defined as

$\displaystyle\mathrm{gdf}({\cal
D})=\frac{1}{M}\sum_{\mu=1}^{M}\frac{\partial\widehat{y}_{\mu}({\cal
D})}{\partial y_{\mu}}.$ (2.41)

GDF coincides with DF for the linear estimator, as we have already observed.
The expression (2.41) is applicable to various models; however, the exact
correspondence between GDF and the number of estimated variables in the model
does not hold in general. The functional form of $\widehat{y}_{\mu}({\cal D})$
depends on the true generative process of the data, assumed model, and
estimation method; therefore, GDF depends on these factors.

####  2.5.1 Shrinkage estimators reduce model complexity

A well-known tendency of GDF is that it decreases under shrinking estimators.
As an example, let us consider ridge regression

$\displaystyle\widehat{\bm{x}}({\cal
D})=\mathop{\mathrm{argmax}}_{\bm{x}}\left\\{-\frac{1}{2}\Big{|}\Big{|}\bm{y}-\frac{1}{\sqrt{N}}\bm{F}\bm{x}\Big{|}\Big{|}_{2}^{2}-\frac{\lambda}{2}||\bm{x}||_{2}^{2}\right\\}.$
(2.42)

The estimate $\widehat{\bm{x}}$ and the expression of data using the estimate
$\widehat{\bm{y}}$ are provided by

$\displaystyle\widehat{\bm{x}}$
$\displaystyle=\frac{1}{\sqrt{N}}\left(\frac{1}{N}\bm{F}^{\top}\bm{F}+\lambda\bm{I}\right)^{-1}\bm{F}^{\top}\bm{y}$
(2.43) $\displaystyle\widehat{\bm{y}}$ $\displaystyle=\bm{H}\bm{y}$ (2.44)

where
$\bm{H}=\frac{1}{N}\bm{F}\left(\frac{1}{N}\bm{F}^{\top}\bm{F}+\lambda\bm{I}\right)^{-1}\bm{F}^{\top}$.
Applying the singular value decomposition
$\frac{1}{\sqrt{N}}\bm{F}=\bm{U}\bm{D}\bm{V}^{\top}$, DF is given by

$\displaystyle\mathrm{df}=\frac{1}{M}\sum_{i=1}^{N}\frac{d_{i}^{2}}{d_{i}^{2}+\lambda},$
(2.45)

where $d_{i}$ represents the $i$-th singular value. Therefore, DF decrease
with an increase in the regularization parameter $\lambda$.

However, in case of the $\ell_{1}$-penalized maximum likelihood method, GDF
corresponds to the number of parameters in the model although the estimator is
shrinking from ordinary least square estimator [34, 35, 36]. This reason is
considered that the selection of the non-zero components increase model
complexity, which is an increase balanced with the decrease of the model
complexity attributed to shrinkage.

###  2.6 Importance sampling cross validation and functional variance

We introduce another expression of the LOOCV error in the framework of the
Bayesian inference; however, we consider the penalized maximum likelihood
method in this paper. In our problem setting, posterior distribution is
defined by

$\displaystyle p_{\beta}(\bm{x}|{\cal
D})=\frac{1}{Z_{\beta}}\exp(\beta\varphi({\cal D},\bm{x})),$ (2.46)

where $Z_{\beta}=\int d\bm{x}\exp(\beta\varphi({\cal D},\bm{x}))$ and $\beta$
is the inverse temperature. In the usual Bayesian inference, $\beta=1$ is
considered. This parameter connects the Bayesian inference and the penalized
maximum likelihood method as

$\displaystyle\widehat{x_{i}}({\cal D})=\lim_{\beta\to\infty}\int
d\bm{x}x_{i}p_{\beta}(\bm{x}|{\cal D}).$ (2.47)

We introduce the importance sampling cross validation (ISCV) error as follows.

######  Definition 8 (Importance sampling cross-validation error).

Under the posterior distribution $p_{\beta}(\bm{x}|{\cal D})$, ISCV error is
defined by

$\displaystyle\mathrm{err}_{\mathrm{ISCV}}({\cal
D};\beta)=\frac{1}{M}\sum_{\mu=1}^{M}\ln\int
d\bm{x}\frac{1}{f^{\beta}(y_{\mu}|\theta_{\mu}(\bm{x},\mathbf{f}_{\mu}))}p_{\beta}(\bm{x}|{\cal
D}).$ (2.48)

######  Theorem 3.

ISCV error is equivalent to the LOOCV error for the Bayesian inference at any
$\beta$ [37, 38].

###### Proof.

We start with the LOOCV error for the Bayesian inference

$\displaystyle\mathrm{err}_{\mathrm{LOOCV}}({\cal
D};\beta)=-\frac{1}{M}\sum_{\mu=1}^{M}\ln\int
d\bm{x}f(y_{\mu}|\theta_{\mu}(\bm{x},\mathbf{f}_{\mu}))p_{\beta}(\bm{x}|{\cal
D}_{\setminus\mu}),$ (2.49)

which is reduced to (2.8) at $\beta\to\infty$. Setting
$\phi(\bm{x})=\exp(-h(\bm{x}))$, (2.49) is transformed as

$\displaystyle\mathrm{err}_{\mathrm{LOOCV}}({\cal D};\beta)$
$\displaystyle=\frac{1}{M}\sum_{\mu=1}^{M}\ln\frac{\int
d\bm{x}\prod_{\nu\in\bm{M}\setminus\mu}f^{\beta}(y_{\nu}|\theta_{\nu}(\mathbf{f}_{\nu},\bm{x}))\phi^{\beta}(\bm{x})}{\int
d\bm{x}\prod_{\nu\in\bm{M}}f^{\beta}(y_{\nu}|\theta_{\nu}(\mathbf{f}_{\nu},\bm{x}))\phi^{\beta}(\bm{x})}$
(2.50) $\displaystyle=\frac{1}{M}\sum_{\mu=1}^{M}\ln\int
d\bm{x}\frac{1}{f^{\beta}(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x}))}p_{\beta}(\bm{x}|{\cal
D}),$

and the R.H.S. corresponds to ISCV error. ∎

For the sake of generality, we define the posterior expectation of the
$\eta$-powered likelihood as

$\displaystyle{\cal
T}(\eta)=\lim_{\beta\to\infty}\frac{1}{M\beta}\sum_{\mu=1}^{M}\ln\int
d\bm{x}p_{\beta}(\bm{x}|{\cal
D})f^{\eta\beta}(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x})).$ (2.51)

The training and ISCV errors correspond to $-{\cal T}(1)$ and ${\cal T}(-1)$,
respectively. Applying Maclaurin expansion to ${\cal T}(\eta)$ as

$\displaystyle{\cal
T}(\eta)=\sum_{k=0}^{\infty}\frac{1}{k!}\frac{\partial^{k}{\cal
T}(\eta)}{\partial\eta^{k}}\Big{|}_{\eta=0}\eta^{k},$ (2.52)

and setting $\eta=\pm 1$, we obtain

$\displaystyle\mathrm{err}_{\mathrm{train}}({\cal D})$
$\displaystyle=-\frac{\partial{\cal
T}(\eta)}{\partial\eta}\Big{|}_{\eta=0}-\frac{1}{2}\frac{\partial^{2}{\cal
T}(\eta)}{\partial\eta^{2}}\Big{|}_{\alpha=0}-\frac{1}{6}\frac{\partial^{3}{\cal
T}(\eta)}{\partial\eta^{3}}\Big{|}_{\alpha=0}+O(\eta^{4})$ (2.53)
$\displaystyle\mathrm{err}_{\mathrm{ISCV}}({\cal D})$
$\displaystyle=-\frac{\partial{\cal
T}(\eta)}{\partial\eta}\Big{|}_{\eta=0}+\frac{1}{2}\frac{\partial^{2}{\cal
T}(\eta)}{\partial\eta^{2}}\Big{|}_{\alpha=0}-\frac{1}{6}\frac{\partial^{3}{\cal
T}(\eta)}{\partial\eta^{3}}\Big{|}_{\alpha=0}+O(\eta^{4}),$ (2.54)

where ${\cal T}(0)=0$. Thus, the following relationship holds:

$\displaystyle\mathrm{err}_{\mathrm{ISCV}}({\cal
D})=\mathrm{err}_{\mathrm{train}}({\cal D})+\frac{\partial^{2}{\cal
T}(\eta)}{\partial\eta^{2}}\Big{|}_{\eta=0}+O(\eta^{4}).$ (2.55)

The second derivative of ${\cal T}$ is given by

$\displaystyle\frac{\partial^{2}}{\partial\eta^{2}}{\cal T}(\eta)$
$\displaystyle=\beta\sum_{\mu=1}^{M}\Big{[}\frac{\mathbb{E}_{\mathrm{post}}\left[\left(\ln
f(y_{\mu}|\theta_{\mu})\right)^{2}\exp(\eta\beta\ln
f(y_{\mu}|\theta_{\mu}))\right]}{\mathbb{E}_{\mathrm{post}}\left[\exp(\eta\beta\ln
f(y_{\mu}|\theta_{\mu}))\right]}$ (2.56)
$\displaystyle-\left\\{\frac{\mathbb{E}_{\mathrm{post}}\left[\ln
f(y_{\mu}|\theta_{\mu})\exp(\eta\beta\ln
f(y_{\mu}|\theta_{\mu}))\right]}{\mathbb{E}_{\mathrm{post}}\left[\exp(\eta\beta\ln
f(y_{\mu}|\theta_{\mu}))\right]}\right\\}^{2}\Big{]},$

which is called functional variance at $\eta=0$ [18].

######  Definition 9 (Functional variance and widely applicable information
criterion).

The functional variance (FV) is defined by the variance of the component-wise
log-likelihood as

$\displaystyle\mathrm{FV}({\cal
D})=\frac{\beta}{M}\sum_{\mu=1}^{M}\Big{\\{}\int d\bm{x}$ $\displaystyle
p_{\beta}(\bm{x}|{\cal D})\left(\ln
f(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x}))\right)^{2}$ (2.57)
$\displaystyle-\left(\int d\bm{x}p_{\beta}(\bm{x}|{\cal D})\ln
f(y_{\mu}|{\theta}_{\mu}(\mathbf{f}_{\mu},\bm{x}))\right)^{2}\Big{\\}}.$

The estimator of the ISCV using the functional variance is known as the widely
applicable information criterion (WAIC)

$\displaystyle\mathrm{WAIC}({\cal D})=\mathrm{err}_{\mathrm{train}}({\cal
D})+\mathrm{FV}({\cal D}).$ (2.58)

The terms of $O(\eta^{3})$ have a negligible contribution in the asymptotic
limit, and hence, WAIC is asymptotically equivalent to the LOOCV or ISCV error
[39]. However, for cases with a large number of parameters, the convergence of
the series expansion (2.52) should be considered, which immediately concerns
the adequacy of WAIC.

For the understanding of the FV, we introduce the following form of FV;

$\displaystyle\mathrm{FV}({\cal
D})=\frac{1}{M}\sum_{\mu=1}^{M}\frac{\partial^{2}}{\partial\gamma_{\mu}^{2}}\ln
Z_{\beta,\bm{\gamma}}({\cal D})\Big{|}_{\bm{\gamma}=1},$ (2.59)

where

$\displaystyle Z_{\beta,\bm{\gamma}}({\cal D})=\int
d\bm{x}f^{\beta\gamma_{\mu}}(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x}))\phi^{\beta}(\bm{x})$
(2.60)

is the extended partition function with an auxiliary variable
$\bm{\gamma}\in\mathbb{R}^{M}$ [40]. This form is useful for the derivation of
FV in the framework of GAMP later.

####  2.6.1 Equivalence of TIC and WAIC

The derivations of TIC and WAIC from the LOOCV error differ from each other.
TIC is obtained by removing one observation and assuming that the difference
in estimates due to the remove is sufficiently small. Meanwhile, WAIC
considers component-wise variance of the log-likelihood. However, it can be
shown that TIC and WAIC are equivalent to each other, as explained in [41],
when the Laplace approximation is valid. By expanding the log-likelihood
around the MAP estimate $\widehat{\bm{x}}({\cal D})$, we obtain

$\displaystyle\ln f(y_{\mu}|\theta_{\mu}$
$\displaystyle(\mathbf{f}_{\mu},\bm{x}))=\ln
f(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\widehat{\bm{x}}({\cal D})))$ (2.61)
$\displaystyle+\sum_{i=1}^{N}\frac{\partial}{\partial x_{i}}\ln
f(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D})}(x_{i}-\widehat{x}_{i}({\cal D})).$

The function to be minimized is expanded as

$\displaystyle\varphi({\cal D},\bm{x})$ $\displaystyle=\varphi({\cal
D},\widehat{\bm{x}}({\cal D}))$ (2.62)
$\displaystyle+\frac{1}{2}\sum_{i,j}\frac{\partial^{2}}{\partial x_{i}\partial
x_{j}}\varphi({\cal D},\bm{x})\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D})}(x_{i}-\widehat{x}_{i}({\cal D}))(x_{j}-\widehat{x}_{j}({\cal D})),$

where we consider that the gradient at $\widehat{\bm{x}}$ is zero. Hence, we
obtain the Laplace approximation of the posterior distribution as

$\displaystyle p_{\beta}(\bm{x}|{\cal D})\propto\exp\left(-\frac{\beta
M}{2}(\bm{x}-\widehat{\bm{x}}({\cal D})){\cal
J}^{\varphi}(\widehat{\bm{x}}({\cal D}))(\bm{x}-\widehat{\bm{x}}({\cal
D}))^{\mathrm{T}}\right),$ (2.63)

where ${\cal J}^{\varphi}(\widehat{\bm{x}})$ is given by (2.23). Using (2.61)
and (2.63), we obtain

$\displaystyle\mathrm{FV}({\cal D})$
$\displaystyle=\frac{\beta}{M}\sum_{\mu=1}^{M}\sum_{ij}\int
d\bm{x}p_{\beta}(\bm{x}|{\cal D})\frac{\partial}{\partial x_{i}}\ln
f(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D})}$ (2.64) $\displaystyle\hskip 28.45274pt\times\frac{\partial}{\partial
x_{j}}\ln
f(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D})}(x_{i}-\widehat{x}_{i}({\cal D}))(x_{j}-\widehat{x}_{j}({\cal D}))$
$\displaystyle=\frac{1}{M}\mathrm{Tr}\left({\cal I}(\widehat{\bm{x}}({\cal
D})){\cal J}^{-1}(\widehat{\bm{x}}({\cal D}))\right),$

which corresponds to the bias term of TIC.

##  3 Our purpose and problem setting

In the conventional problem settings in statistics, $M\to\infty$ with
sufficiently small $N$ are considered, hence the difference between the in-
sample prediction error and the extra-sample prediction error has not been
focused on. Recently, characteristic phenomena known as double descent
phenomena, coined by [42], wherein the extra-sample prediction error exhibits
two minima: one each in the underparameterized and overparametrized regions,
has been studied [7]. These phenomena were reported in [43, 44, 45] for linear
models, and have since been reported [46, 47] for neural networks, attracting
considerable attention. For a broader understanding of these phenomena, some
studies are attempted on the random features model [48, 49, 50, 51]. In
general, the modern neural network considers the overparametrized region, and
some studies attempted to extend the conventional statistical methods to the
overparametrized region [9]. Following such direction, we focus on the
differences between the estimator of the prediction errors, namely, $C_{p}$,
TIC or WAIC, and LOOCV error, in the parameter region $M\sim O(N)$, including
the regime where the number of parameters is larger than the data dimension.
In the subsequent sections, we derive their forms using the generalized
approximate message passing (GAMP) algorithm [52, 53], and discuss their
behavior in the $M\sim O(N)$ region.

Hereafter, we denote $\alpha=M/N$. For the dense regression coefficient,
$\alpha$ corresponds to the ratio of data dimension to the number of
parameters, and $\alpha\gg 1$ is the parameter region where the conventional
studies of statistics focused on. We focus on the finite $\alpha$ region where
$\alpha\sim O(1)$. In this paper, we consider the generalized linear model
(GLM), and demonstrate the results for the penalized linear regression and
penalized logistic regression. The forms shown in this study can be applied to
general exponential family distributions.

We here summarize our contribution shown in the following sections.

  * •

We derive the forms of the estimators for the generalization gaps using GAMP,
which is applicable to any $\alpha$ region as long as the local stability
condition of GAMP is not violated.

  * •

We show that GDF and FV are related to the variance of the estimated parameter
in the framework of GAMP.

  * •

Using the expression by GAMP, it is shown that the convergence condition of
the series expansion (2.52), which is the basis of WAIC, is related to GDF in
case of the Gaussian likelihood.

  * •

The resulting forms of the estimators for the in-sample prediction error have
differences with the extra-sample prediction error in the finite $\alpha$
region. We derive the differences in the finite $\alpha$ region using the
expression by GAMP.

###  3.1 Problem setting

We consider GLM as the likelihood where the assumed model is provided by the
natural exponential family distribution;

$\displaystyle\ln
f(\bm{y}|\bm{\theta})=\bm{\theta}^{\top}\bm{y}-a(\bm{\theta})+b(\bm{y}),$
(3.1)

where $a(\bm{\theta})=\sum_{\mu=1}^{M}a(\theta_{\mu})$ and
$b(\bm{y})=\sum_{\mu=1}^{M}b(y_{\mu})$. In the natural exponential family,
$\bm{\theta}$ relates to the statistics of the distribution as

$\displaystyle\mathbb{E}_{{\cal L}(\bm{\theta})}[y_{\mu}]$
$\displaystyle=\frac{\partial}{\partial\theta_{\mu}}a(\theta_{\mu}),~{}~{}~{}\mathbb{V}_{{\cal
L}(\bm{\theta})}[y_{\mu}]$
$\displaystyle=\frac{\partial^{2}}{\partial^{2}\theta_{\mu}}a(\theta_{\mu}),$
(3.2)

where $\mathbb{E}_{{\cal L}(\bm{\theta})}[\cdot]$ and $\mathbb{V}_{{\cal
L}(\bm{\theta})}[\cdot]$ denote the mean and variance according to the
likelihood $f$ with the parameter $\bm{\theta}$. In GLM, the model parameter
and the mean of the distribution are related to each other through link
function $\ell$ as

$\displaystyle\ell(\mathbb{E}_{{\cal
L}(\bm{\theta})}[y_{\mu}])=\theta_{\mu},~{}~{}~{}\ell^{-1}(\theta_{\mu})=\mathbb{E}_{{\cal
L}(\bm{\theta})}[y_{\mu}].$ (3.3)

The distribution is specified by the functions $a(\cdot)$ and $b(\cdot)$. For
example, the linear regression with Gaussian likelihood corresponds to

$\displaystyle a(\bm{\theta})$
$\displaystyle=\frac{1}{2}||\bm{\theta}||_{2}^{2},~{}~{}~{}b(\bm{y})$
$\displaystyle=-\frac{1}{2}||\bm{y}||_{2}^{2},$ (3.4)

and the link function is given by the identity function $\ell(\theta)=\theta$.
Other representative GLMs are presented in B. The generalization gap of the
in-sample prediction error for GLM is written in the same form as that of
(2.37):

$\displaystyle\mathrm{err}_{\mathrm{train}}({\cal D})$
$\displaystyle=-\frac{1}{M}\ln
f(\bm{y}|\bm{\theta}(\bm{F},\widehat{\bm{x}}({\cal D})))$ (3.5)
$\displaystyle=\frac{1}{M}\left\\{-\bm{\theta}(\bm{F},\widehat{\bm{x}}({\cal
D}))^{\top}\bm{y}+a(\bm{\theta}(\bm{F},\widehat{\bm{x}}({\cal
D})))-b(\bm{y})\right\\}$ $\displaystyle{\rm err}_{\rm
pre}^{(\mathrm{in})}({\cal D})$
$\displaystyle=-\frac{1}{M}\mathbb{E}_{\bm{z}\sim q(\bm{z})}[\ln
f(\bm{z}|\bm{\theta}(\bm{F},\widehat{\bm{x}}({\cal D})))]$ (3.6)
$\displaystyle=\frac{1}{M}\Big{\\{}-\bm{\theta}(\bm{F},\widehat{\bm{x}}({\cal
D}))^{\top}\mathbb{E}_{\bm{z}\sim
q(\bm{z})}[\bm{z}]+a(\bm{\theta}(\bm{F},\widehat{\bm{x}}({\cal D})))$
$\displaystyle\hskip 142.26378pt-\mathbb{E}_{\bm{z}\sim
q(\bm{z})}[b(\bm{z})]\Big{\\}},$

hence

$\displaystyle\Delta^{(\mathrm{in})}({\cal
D})=\frac{1}{M}\mathrm{Cov}_{\bm{y}}(\bm{y},|\bm{\theta}(\bm{F},\widehat{\bm{x}}({\cal
D}))).$ (3.7)

For the penalty term, we restrict our discussion to the separable penalties,
where $h(\bm{x})=\sum_{i=1}^{N}h(x_{i})$. We focus on the elastic net penalty

$\displaystyle h(x)=\lambda_{1}|x|+\frac{\lambda_{2}}{2}x^{2},$ (3.8)

where $\lambda_{1}$ and $\lambda_{2}$ are regularization parameters for
controlling the intensity of the penalty. The elastic net penalty at
$\lambda_{1}=0$ and $\lambda_{2}=0$ correspond to $\ell_{2}$ and $\ell_{1}$
norm penalty, respectively.

##  4 GAMP for GLM

Message passing is an algorithm used to approximately calculate the marginal
posterior distribution. The algorithm is based on the graphical representation
of the conditional independency of the variables according to the probability
distribution [54, 55]. Graphical models provide visual support to draw
inferences and construct efficient algorithms. In practice, message passing
shows efficient performances as a decoding algorithm for low density parity
check (LDPC) code [56], code division multiple access (CDMA) [57], and others.

Message passing with a Gaussian approximation is known as approximate message
passing (AMP) [58] that reduces the computational cost required for the
inference. The AMP was proposed in the study on compressed sensing and it was
extended to a generalized approximate message passing (GAMP) [52] to deal with
general separable likelihood and prior. The AMP has a similarity between the
cavity approach [59] and Thouless–Anderson–Palmer (TAP) equation [60] in
statistical physics for spin glasses. This relationship induces a statistical-
physics based study on the information science and machine learning [61] such
as compressed sensing [14], regression [62], and matrix factorization [63].
The details of the derivation of the GAMP for the regression problem of the
separable likelihoods and priors are shown in C.

The message passing procedure provides the exact marginal distribution when
the graph consists of tree structures. Our problem setting corresponds to the
case where the graph has loops, hence the marginal distribution calculated by
the message passing is an approximation. More precisely, the resulting
marginal posterior distributions are locally consistent marginals, and do not
necessarily constitute an overall posterior distribution [61]. Furthermore, it
is known that GAMP with random initialization only provides an accurate
approximation of the posterior distribution only in the case in which the
global minima of the free energy is the closest minima to zero overlap [64,
65]. Although a lack of the exactness and limitations exist in accurate
approximation, GAMP is useful for calculating the generalization gap. In the
estimation of the in-sample prediction error, analytical expression of the
degrees of freedom for general regression problems are limited, and specific
expressions are studied for LASSO [34], group lasso, fused lasso [66], and
minimax concave penalty [67]. Numerically, the bootstrap method is useful for
evaluating the generalization gap [16]; however, its computational cost is
expensive, and it often leads unstable estimates. Extension of the path
seeking algorithm [68] is another numerical approach that is applied for
estimating the generalization gap [69], which requires additional matrix
operation. Compared with these numerical approaches, GAMP has several
advantages as an estimation method of the generalization gap for the in-sample
prediction error. First, GAMP provides a general expression of the
generalization gap for the in-sample prediction error. This expression leads
to a unified interpretation of the prediction as a relationship between the
fluctuation and response, as discussed in sec.4.4. Second, GAMP does not
require additional computation for evaluating the generalization gap for the
in-sample prediction error, as the generalization gap can be expressed using
variables that constitute GAMP.

In addition to the advantages of evaluating the in-sample prediction error, we
can obtain LOOCV error using GAMP by assuming that the cavity distribution in
GAMP corresponds to the posterior distribution under the LOO sample. While the
verification of the assumption is required, it greatly reduces the
computational cost required to estimate the LOOCV error, as discussed in
sec.4.3.3.

###  4.1 General form of message passing algorithm

Figure 3: Graphical representation of regression problem and messages defined
on the edges

The graphical representation of the regression problems considered here is
shown in Figure 3. Factor nodes $\square$, variable nodes $\bigcirc$, and the
edges represent the data $\bm{y}$, regression coefficient $\bm{x}$, and
predictor matrix $\bm{F}$, respectively. The message passing procedure is
described by following two messages defined on the edges between the factor
nodes and the variable nodes:

  * •

${\xi}_{i\to\mu}(x_{i})$: Message from the variable $x_{i}$ to data $y_{\mu}$
(output message)

  * •

$\widetilde{\xi}_{\mu\to i}(x_{i})$: Message from data $y_{\mu}$ to the
variable $x_{i}$ (input message)

Intuitively, output and input messages are the probability distributions
before and after the observation of $\mu$-th data. Hence, in particular, the
output message ${\xi}_{i\to\mu}(x_{i})$ is considered to describe the
probability distribution under the LOO sample ${\cal D}_{\setminus\mu}$ [70].
We term the graphical model without $\mu$-th data as $\mu$-cavity system.

For general graphs, these messages need to be recursively updated to obtain
the posterior distribution. Here, we place the superscript $(t)$ on the
messages and related quantities to indicate that their values are the results
of the $t$-times update. The messages at the $t$-th and $(t+1)$-th step are
related as

$\displaystyle\widetilde{\xi}^{(t)}_{\mu\to i}(x_{i})$
$\displaystyle\propto\int d\bm{x}_{\backslash
i}f^{\beta}(y_{\mu}|\bm{\theta})\prod_{j\in\bm{N}\setminus
i}{\xi}^{(t)}_{j\to\mu}(x_{j})$ (4.1)
$\displaystyle{\xi}_{i\to\mu}^{(t+1)}(x_{i})$
$\displaystyle\propto\phi^{\beta}(x_{i})\prod_{\nu\in\bm{M}\setminus\mu}\widetilde{\xi}_{\nu\to
i}^{(t)}(x_{i}),$ (4.2)

where $\bm{N}=\\{1,\cdots,N\\}$, $\bm{M}=\\{1,\cdots,M\\}$, and
$\bm{N}\setminus i$ denote $\bm{N}$ except $i$. Using the input messages, the
marginal distribution is given by

$\displaystyle
p^{(t)}_{\beta,i}(x_{i})\propto\phi^{\beta}(x_{i})\prod_{\nu\in\bm{M}}\widetilde{\xi}^{(t)}_{\nu\to
i}(x_{i}),$ (4.3)

and the $i$-th component of the penalized maximum likelihood estimate is given
by

$\displaystyle\widehat{x}^{(t)}_{i}$ $\displaystyle=\lim_{\beta\to\infty}\int
d{x_{i}}~{}x_{i}p_{\beta,i}^{(t)}(x_{i}).$ (4.4)

###  4.2 Approximate message passing

The updating rule (4.1) and (4.2) consists of $2MN$ messages. For
simplification, we introduce the following assumptions for sufficiently large
$M$ and $N$.

a1:

All components of the predictor matrix are $O(1)$.

a2:

The correlation between predictors is negligible: the off-diagonal components
of $\frac{1}{N}\bm{F}^{\top}\bm{F}$ are $O(N^{-1/2})$.

Under assumptions a1 and a2, the input messages $\\{\widetilde{\xi}_{\mu\to
i}\\}$ can be considered Gaussian distributions, as shown in C. Therefore, it
is sufficient to consider the means and variances of the messages; we denote
them as

$\displaystyle\widehat{x}_{i\to\mu}^{(t)}$ $\displaystyle=\int
dx_{i}x_{i}\xi_{i\to\mu}^{(t)}(x_{i})$ (4.5) $\displaystyle s_{i\to\mu}^{(t)}$
$\displaystyle=\beta\int
dx_{i}(x_{i}-\widehat{x}_{i\to\mu}^{(t)})^{2}\xi_{i\to\mu}^{(t)}(x_{i})$ (4.6)
$\displaystyle\widetilde{x}_{\mu\to i}^{(t)}$ $\displaystyle=\int
dx_{i}x_{i}\widehat{\xi}_{\mu\to i}^{(t)}(x_{i})$ (4.7)
$\displaystyle\widetilde{s}_{\mu\to i}^{(t)}$ $\displaystyle=\beta\int
dx_{i}(x_{i}-\widetilde{x}_{\mu\to i}^{(t)})^{2}\widetilde{\xi}_{\mu\to
i}^{(t)}(x_{i}).$ (4.8)

An algorithm for updating these means and variances under the arbitrary
likelihood and prior distribution is called generalized approximate message
passing (GAMP). In GAMP, the marginalized distribution $p_{\beta,i}(x_{i})$ is
represented by

$\displaystyle p_{\beta,i}^{(t)}(x_{i})={\cal
P}_{\beta}(x_{i};\Sigma_{i}^{(t)},\mathrm{m}_{i}^{(t)}),$ (4.9)

where ${\cal P}_{\beta}$ is the distribution defined by

$\displaystyle{\cal P}_{\beta}(x;\Sigma,\mathrm{m})=\frac{1}{{\cal
Z}}\phi^{\beta}(x)\exp\left\\{-\frac{\beta(x-\mathrm{m})^{2}}{2\Sigma}\right\\},$
(4.10)

with the normalization constant ${\cal Z}$. We define the mean and variance of
(4.10) as $\mathbb{M}_{\beta}(\Sigma,\mathrm{m})$ and
$\beta^{-1}\mathbb{V}_{\beta}(\Sigma,\mathrm{m})$, respectively, and denote

$\displaystyle\widehat{x}_{i}^{(t)}$ $\displaystyle=\int dx_{i}x_{i}{\cal
P}_{\beta}(x;\Sigma_{i}^{(t)},\mathrm{m}_{i}^{(t)})=\mathbb{M}_{\beta}(\Sigma_{i}^{(t)},\mathrm{m}_{i}^{(t)})$
(4.11) $\displaystyle s_{i}^{(t)}$ $\displaystyle=\beta\int
dx_{i}(x_{i}-\widehat{x}_{i}^{(t)})^{2}{\cal
P}_{\beta}(x;\Sigma_{i}^{(t)},\mathrm{m}_{i}^{(t)})=\mathbb{V}_{\beta}(\Sigma_{i}^{(t)},\mathrm{m}_{i}^{(t)}),$
(4.12)

where $\widehat{x}_{i}^{(t)}$ corresponds to the penalized maximum likelihood
estimate of step $t$ at $\beta\to\infty$, which limit will be taken later.
Here, the following relationship holds.

$\displaystyle\mathbb{V}_{\beta}(\Sigma,\mathrm{m})=\beta^{-1}\Sigma\frac{\partial}{\partial\mathrm{m}}\mathbb{M}_{\beta}(\Sigma,\mathrm{m}).$
(4.13)

The form of (4.10) indicates that the posterior distribution is given by the
prior distribution for $x_{i}$ and the Gaussian distribution with mean
$\mathrm{m}_{i}$ and variance $\Sigma_{i}$, which represents a contribution
from other variables. Here, we set

$\displaystyle\widehat{\theta}_{\mu}^{\setminus\mu(t)}$
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}F_{\mu
i}\widehat{x}_{i\to\mu}^{(t)},$ (4.14)
$\displaystyle{s_{\theta}}_{\mu}^{\setminus\mu(t)}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}F_{\mu i}^{2}s_{i\to\mu}^{(t)},$
(4.15)

where (4.14) corresponds to the estimated model parameter in the $\mu$-cavity
system. After some calculations shown in C, we derive $\mathrm{m}_{i}$ and
$\Sigma_{i}$ as

$\displaystyle\mathrm{m}_{i}^{(t)}$
$\displaystyle={\Sigma}_{i}^{(t)}\sum_{\mu=1}^{M}\Big{\\{}\frac{1}{\sqrt{N}}F_{\mu
i}\beta^{-1}g_{\mathrm{out}}({s_{\theta}}_{\mu}^{\setminus\mu(t)},\widehat{\theta}_{\mu}^{\setminus\mu(t)},y_{\mu})$
(4.16) $\displaystyle\hskip 85.35826pt-\frac{1}{N}F_{\mu
i}^{2}\widehat{x}_{i}^{(t)}\beta^{-1}\partial_{\widehat{\theta}}g_{\mathrm{out}}({s_{\theta}}_{\mu}^{\setminus\mu(t)},\widehat{\theta}_{\mu}^{\setminus\mu(t)},y_{\mu})\Big{\\}}$
$\displaystyle\frac{1}{{\Sigma}_{i}^{(t)}}$
$\displaystyle=-\frac{1}{N}\sum_{\mu=1}^{M}F^{2}_{\mu
i}\beta^{-1}\partial_{\widehat{\theta}}g_{\mathrm{out}}({s_{\theta}}_{\mu}^{\setminus\mu(t)},\widehat{\theta}_{\mu}^{\setminus\mu(t)},y_{\mu}),$
(4.17)

where

$\displaystyle
g_{\mathrm{out}}({s_{\theta}}_{\mu}^{\setminus\mu(t)},\widehat{\theta}_{\mu}^{\setminus\mu(t)},y_{\mu})$
$\displaystyle=\Big{\langle}\frac{\beta(\theta_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu(t)})}{{s_{\theta}}_{\mu}^{\setminus\mu(t)}}\Big{\rangle}_{\theta_{\mu}}$
(4.18)
$\displaystyle\partial_{\widehat{\theta}}g_{\mathrm{out}}({s_{\theta}}_{\mu}^{\setminus\mu(t)},\widehat{\theta}_{\mu}^{\setminus\mu(t)},y_{\mu})$
$\displaystyle=\frac{\partial}{\partial\widehat{\theta}_{\mu}^{\setminus\mu(t)}}g_{\mathrm{out}}({s_{\theta}}_{\mu}^{\setminus\mu(t)},\widehat{\theta}_{\mu}^{\setminus\mu(t)},y_{\mu}).$
(4.19)

The expectation $\langle\cdot\rangle_{\theta_{\mu}}$ in (4.18) denotes the
expectation with respect to $\theta_{\mu}$ based on the distribution

$\displaystyle
p_{\mathrm{c},\mu}^{(t)}(\theta_{\mu}|y_{\mu})=\displaystyle\frac{f^{\beta}(y_{\mu}|\theta_{\mu})\psi^{\beta}_{\mathrm{c}}(\theta_{\mu}|{s_{\theta}}_{\mu}^{\setminus\mu(t)},\widehat{\theta}_{\mu}^{\setminus\mu(t)})}{\int
d\theta_{\mu}f^{\beta}(y_{\mu}|\theta_{\mu})\psi^{\beta}_{\mathrm{c}}(\theta_{\mu}|{s_{\theta}}_{\mu}^{\setminus\mu(t)},\widehat{\theta}_{\mu}^{\setminus\mu^{(t)}})},$
(4.20)

where

$\displaystyle\psi_{\mathrm{c}}(\theta_{\mu}|{s_{\theta}}_{\mu}^{\setminus\mu},\widehat{\theta}^{\setminus\mu}_{\mu})=\exp\left(-\frac{(\theta_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu})^{2}}{2{s_{\theta}}_{\mu}^{\setminus\mu}}\right).$
(4.21)

Defining the following quantities

$\displaystyle\widehat{\theta}^{(t)}_{\mu}$
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}F_{\mu i}\widehat{x}_{i}^{(t)}$
(4.22) $\displaystyle{s_{\theta}}_{\mu}^{(t)}$
$\displaystyle=\frac{1}{N}\sum_{j=1}^{N}F_{\mu j}^{2}s_{j}^{(t)},$ (4.23)

where (4.22) corresponds to the estimated model parameter, their values are
related to those for cavity distribution as

$\displaystyle\widehat{\theta}_{\mu}^{(t)}$
$\displaystyle=\widehat{\theta}_{\mu}^{\setminus\mu(t)}+\beta^{-1}(g_{\mathrm{out}}^{(t-1)})_{\mu}{s_{\theta}}_{\mu}^{\setminus\mu(t)}$
(4.24) $\displaystyle{s_{\theta}}_{\mu}^{(t)}$
$\displaystyle={s_{\theta}}_{\mu}^{\setminus\mu(t)},$ (4.25)

where the second term of (4.24) is known as Onsager reaction term [60].

###  4.3 $\beta\to\infty$ limit

At $\beta\to\infty$, the saddle point method can be used to calculate the
integrals with respect to $x$ and $\theta$ based on the distributions (4.9)
and (4.20). First, the mean of the distribution (4.9) denoted by $\widehat{x}$
is given as the solution of the maximization problem

$\displaystyle\widehat{x}_{i}=\mathop{\mathrm{argmax}}_{x}f_{\xi}(x;\mathrm{m}_{i},\Sigma_{i})$
(4.26)

where

$\displaystyle
f_{\xi}(x;\mathrm{m},\Sigma)=\ln\phi(x)-\frac{(x-\mathrm{m})^{2}}{2\Sigma}.$
(4.27)

The solution depends on the prior distribution. Here, we show the solution
under the elastic net penalty as

$\displaystyle\widehat{x}_{i}$
$\displaystyle=\frac{\mathrm{m}_{i}-\mathrm{sgn}(\mathrm{m}_{i})\lambda_{1}\Sigma_{i}}{\lambda_{2}\Sigma_{i}+1}\mathbb{I}(|\mathrm{m}_{i}|>\lambda_{1}\Sigma_{i}),$
(4.28) $\displaystyle s_{i}$
$\displaystyle=\frac{\Sigma_{i}}{\lambda_{2}\Sigma_{i}+1}\mathbb{I}(|\mathrm{m}_{i}|>\lambda_{1}\Sigma_{i}),$
(4.29)

where $\mathbb{I}(a)$ represents the indicate function returns 1 if $a$ is
true, otherwise it is 0. Setting $\lambda_{1}=0$ and $\lambda_{2}=0$, we
obtain the solution under the $\ell_{2}$ penalty and $\ell_{1}$ penalty,
respectively.

Next, the mean (4.20) at $\beta\to\infty$, denoted by $\theta_{\mu}^{*}$, is
given by

$\displaystyle\theta_{\mu}^{*}$
$\displaystyle=\mathop{\mathrm{argmax}}_{\theta}\widehat{f}_{\xi}(\theta;y_{\mu},\widehat{\theta}_{\mu}^{\setminus\mu},{s_{\theta}}_{\mu}^{\setminus\mu}),$
(4.30)

with the function $\widehat{f}_{\xi}$ defined as

$\displaystyle\widehat{f}_{\xi}(\theta;y,\widehat{\theta},s)=\ln
f(y|\theta)-\frac{(\theta-\widehat{\theta})^{2}}{2s}.$ (4.31)

The maximizer $\theta_{\mu}^{*}$ in (4.30) satisfies

$\displaystyle\frac{\theta_{\mu}^{*}-\widehat{\theta}_{\mu}^{\setminus\mu}}{{s_{\theta}}_{\mu}^{\setminus\mu}}=y_{\mu}-a^{\prime}(\theta_{\mu}^{*}),$
(4.32)

where
$a^{\prime}(\theta_{\mu}^{*})=\frac{\partial}{\partial\theta}a(\theta)|_{\theta=\theta_{\mu}^{*}}$
corresponds to the mean of $y_{\mu}$ expressed by the assumed model, i.e.
$a^{\prime}(\theta_{\mu}^{*})=\mathbb{E}_{{\cal
L}(\bm{\theta}^{*})}[y_{\mu}]=\ell^{-1}(\theta_{\mu}^{*})$. From (4.32),
$\partial_{\hat{\theta}}\theta_{\mu}^{*}\equiv\partial\theta_{\mu}^{*}/\partial\widehat{\theta}_{\mu}^{\setminus\mu}$
satisfies

$\displaystyle\frac{\partial_{\hat{\theta}}\theta_{\mu}^{*}-1}{{s_{\theta}}_{\mu}}=-\partial_{\widehat{\theta}}\theta_{\mu}^{*}a^{\prime\prime}(\theta_{\mu}^{*}),$
(4.33)

where $a^{\prime\prime}(\theta_{\mu}^{*})$ corresponds to the variance of
$y_{\mu}$ under the assumed model, i.e.
$a^{\prime\prime}(\theta_{\mu}^{*})=\mathbb{V}_{{\cal
L}(\bm{\theta}^{*})}[y_{\mu}]$. Solving (4.33), we obtain

$\displaystyle\frac{\partial_{\widehat{\theta}}\theta_{\mu}^{*}-1}{{s_{\theta}}_{\mu}}=-\frac{\mathbb{V}_{{\cal
L}(\theta^{*}_{\mu})}[y_{\mu}]}{1+{s_{\theta}}_{\mu}\mathbb{V}_{{\cal
L}(\theta^{*}_{\mu})}[y_{\mu}]}.$ (4.34)

For convenience, we introduce rescaled variables

$\displaystyle\beta^{-1}(g_{\mathrm{out}})_{\mu}$
$\displaystyle=(\check{g}_{\mathrm{out}})_{\mu},~{}~{}~{}\beta^{-1}(\partial_{\hat{\theta}}g_{\mathrm{out}})_{\mu}$
$\displaystyle=(\check{\partial}_{\theta}g_{\mathrm{out}})_{\mu},$ (4.35)

and they are given by

$\displaystyle(\check{g}_{\rm
out})_{\mu}=\frac{\theta_{\mu}^{*}-\widehat{\theta}_{\mu}^{\setminus\mu}}{{s_{\theta}}_{\mu}},~{}~{}~{}(\check{\partial}_{\hat{\theta}}g_{\rm
out})_{\mu}=\frac{\partial_{\hat{\theta}}\theta_{\mu}^{*}-1}{{s_{\theta}}_{\mu}},$
(4.36)

at $\beta\to\infty$. Based on scaling with respect to $\beta$, (4.6), (4.8)
and (4.35), the variables in GAMP stay finite even at the limit
$\beta\to\infty$. This scaling can be understood through Laplace approximation
of the posterior distribution around the MAP estimates. From eqs.(4.32) and
(4.33), we obtain

$\displaystyle(\check{g}_{\mathrm{out}})_{\mu}$
$\displaystyle=y_{\mu}-\mathbb{E}_{{\cal L}}(\theta_{\mu}^{*})[y_{\mu}],$
(4.37)
$\displaystyle(\check{\partial}_{\widehat{\theta}}g_{\mathrm{out}})_{\mu}$
$\displaystyle=-\frac{\mathbb{V}_{{\cal
L}(\theta^{*}_{\mu})}[y_{\mu}]}{1+{s_{\theta}}_{\mu}\mathbb{V}_{{\cal
L}(\theta^{*}_{\mu})}[y_{\mu}]}.$ (4.38)

The solution $\theta^{*}_{\mu}$ of the maximization problem (4.30) is
equivalent to the estimate for the parameter $\widehat{\theta}_{\mu}$ defined
by (4.22), which can be confirmed by substituting (4.37) to (4.24).

The analytical form of the estimate $\bm{\theta}^{*}$ and related quantities
depend on the likelihood. For Gaussian likelihood, the estimate is given by

$\displaystyle\theta_{\mu}^{*}=\frac{{s_{\theta}}_{\mu}y_{\mu}}{1+{s_{\theta}}_{\mu}}+\frac{\widehat{\theta}_{\mu}^{\setminus\mu}}{1+{s_{\theta}}_{\mu}}.$
(4.39)

Further, $\mathbb{E}_{{\cal L}}(\theta_{\mu}^{*})[y_{\mu}]=\theta_{\mu}^{*}$
and $\mathbb{V}_{{\cal L}}(\theta_{\mu}^{*})[y_{\mu}]=1$, and hence,

$\displaystyle(\check{g}_{\mathrm{out}})_{\mu}$
$\displaystyle=\frac{y_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu}}{1+{s_{\theta}}_{\mu}},~{}~{}~{}(\check{\partial}_{\widehat{\theta}}g_{\mathrm{out}})_{\mu}$
$\displaystyle=-\frac{1}{1+{s_{\theta}}_{\mu}}.$ (4.40)

We summarise the GAMP solutions for representative distributions in the
exponential family in B. The procedure of GAMP for general likelihoods and
priors is summarised in Algorithm 1, where $\circ$ and $\oslash$ denote the
component-wise product and division, respectively, and $\bm{0}_{N\times 1}$
denotes the $N$-dimensional vector whose components are zero.

Algorithm 1 GAMP for penalized regression problems

1:${\cal D}=\\{\bm{y},\bm{F}\\}$

2:$\widehat{\bm{x}}$

3:$\widehat{x}^{(-1)}\leftarrow\bm{0}_{N\times 1}$,
$\bm{g}^{(-1)}_{\mathrm{out}}\leftarrow\bm{0}_{M\times 1}$

4:$\widehat{\bm{x}}^{(0)},~{}\bm{s}^{(0)}\leftarrow$ Initial values from
$\mathbb{R}^{N}$

5:for $t=0\,\ldots\,T_{\max}$ do

6:
$\widehat{\bm{\theta}}^{(t)}\leftarrow\frac{1}{\sqrt{N}}\bm{F}\widehat{\bm{x}}^{(t)}-\frac{1}{N}\bm{g}_{\mathrm{out}}^{(t-1)}\circ\left((\bm{F}\circ\bm{F})\bm{s}^{(t)}\right)$

7: $\bm{s}_{\theta}^{(t)}\leftarrow\frac{1}{N}(\bm{F}\circ\bm{F})\bm{s}^{(t)}$

8:
$\bm{g}_{\mathrm{out}}^{(t)}\leftarrow\bm{g}_{\mathrm{out}}(\bm{s}_{\theta}^{(t)},\widehat{\bm{\theta}}^{(t)};\bm{y})$

9:
$\partial_{\hat{\theta}}\bm{g}_{\mathrm{out}}^{(t)}\leftarrow\partial_{\widehat{\theta}}\bm{g}_{\mathrm{out}}(\bm{s}_{\theta}^{(t)},\widehat{\bm{\theta}}^{(t)};\bm{y})$

10: $\bm{\Sigma}^{(t)}\leftarrow-\bm{1}_{N\times
1}\oslash\left(\frac{1}{N}(\bm{F}\circ\bm{F})^{\top}\partial_{\hat{\theta}}\bm{g}_{\mathrm{out}}^{(t)}\right)$

11:
$\bm{\mathrm{m}}^{(t)}\leftarrow\bm{\Sigma}^{(t)}\circ\left(\frac{1}{\sqrt{N}}\bm{F}^{\top}\bm{g}_{\mathrm{out}}^{(t)}-\frac{1}{N}\widehat{\bm{x}}^{(t)}\circ((\bm{F}\circ\bm{F})^{\top}\partial_{\hat{\theta}}\bm{g}_{\mathrm{out}}^{(t)})\right)$

12:
$\widehat{\bm{x}}^{(t+1)}\leftarrow\mathbb{M}_{\beta}(\bm{\Sigma}^{(t)},\bm{\mathrm{m}}^{(t)})$

13:
$\bm{s}^{(t+1)}\leftarrow\mathbb{V}_{\beta}(\bm{\Sigma}^{(t)},\bm{\mathrm{m}}^{(t)})$

14:end for

####  4.3.1 Neglecting a part of site-dependency for random predictors

For random predictors with i.i.d. random entries of zero mean and variance 1,
${s_{\theta}}_{\mu}$ converges to the averaged value of $s_{i}$

$\displaystyle{s_{\theta}}_{\mu}=\frac{1}{N}\sum_{i=1}^{N}F_{\mu
i}^{2}s_{i}\to s_{\theta}=\frac{1}{N}\sum_{i=1}^{N}s_{i}$ (4.41)

at sufficiently large $N$ because of the law of large numbers. Further, the
site dependency of $\Sigma_{i}$ can be ignored as

$\displaystyle\frac{1}{{\Sigma}^{(t)}}$
$\displaystyle=-\frac{\alpha}{M}\sum_{\mu=1}^{M}\check{\partial}_{\widehat{\theta}}g_{\mathrm{out}}({s_{\theta}}_{\mu}^{\setminus\mu(t)},\widehat{\theta}_{\mu}^{\setminus\mu(t)},y_{\mu}).$
(4.42)

####  4.3.2 Comparison of GAMP with the linear algebraic expression

Now, we compare the expression of the LOOCV error obtained using GAMP and that
obtained from the linear algebraic form given by (2.9) for the Gaussian
likelihood. In GAMP, we denote the model expression of data under the
$\mu$-cavity system as
$\widehat{y}_{\mu}^{\setminus\mu}=\frac{1}{\sqrt{N}}F_{\mu
i}\widehat{x}_{i\to\mu}$; this corresponds to
$\widehat{\theta}_{\mu}^{\setminus\mu}$. Hence, we obtain

$\displaystyle
y_{\mu}-\widehat{y}_{\mu}^{\setminus\mu}=(y_{\mu}-\widehat{y}_{\mu})(1+{s_{\theta}}_{\mu}).$
(4.43)

Meanwhile, from (2.15), we obtain

$\displaystyle y_{\mu}-\widehat{y}_{\mu}({\cal
D}_{\setminus\mu})=\frac{y_{\mu}-\widehat{y}_{\mu}({\cal D})}{1-H_{\mu\mu}}.$
(4.44)

Assuming the equivalence between $\widehat{y}_{\mu}^{\setminus\mu}$ in GAMP
and $\widehat{y}_{\mu}({\cal D}_{\setminus\mu})$ in the exact calculation, the
correspondence between ${s_{\theta}}_{\mu}$ in GAMP and the matrix $\bm{H}$ is
given as

$\displaystyle
H_{\nu\nu}\leftrightarrow\frac{{s_{\theta}}_{\nu}}{1+{s_{\theta}}_{\nu}}$
(4.45)

for any $\nu\in\bm{M}$. For computing the hat matrix in the linear estimator,
we need to calculate the matrix inverse. An advantage of GAMP is that there is
no need to compute the inverse of the matrix to obtain the estimate, although
a recursive update is required. Further, the R.H.S of (4.45) can be calculated
for any estimator, even when the linear estimation rule is not considered and
the matrix $\bm{H}$ is not defined; hence, (4.45) can be regarded as the
extension of the matrix $\bm{H}$ for other estimators.

####  4.3.3 Correspondence between $\mu$-cavity system and LOO sample

Here, we discuss the correspondence between $\widehat{x}_{i\to\mu}$ in GAMP
and exact $\widehat{x}_{i}({\cal D}_{\mu})$. When these estimates are
equivalent to each other at any $i$ and $\mu$, the parameters
$\widehat{\theta}_{\mu}^{\setminus\mu}$ and $\widehat{\theta}({\cal
D}_{\setminus\mu})$ coincide. In the case where the graph structure consists
of tree structures, the exactness of the output message $\xi_{i\to\mu}(x_{i})$
as the posterior distribution under the LOO sample ${\cal D}_{\setminus\mu}$
can be mathematically proved [61], hence,
$\widehat{x}_{i\to\mu}=\widehat{x}({\cal D}_{\setminus\mu})$ holds. In the
current problem setting, the graph contains loops, and the correspondence does
not hold in general. Hereafter, we use $\widehat{\theta}_{\mu}^{\setminus\mu}$
as an estimated value of $\widehat{\theta}({\cal D}_{\setminus\mu})$, and the
related generalization gap of LOOCV error given by the estimated value is
denoted as $\widehat{\Delta}_{\mathrm{LOOCV}}$.

####  4.3.4 Local stability condition of GAMP

The GAMP algorithm sometimes fails to converge depending on the model
parameters $N$, $M$, and regularization parameters. The instability condition
can be derived by considering the linear stability analysis around the GAMP’s
fixed point.

We consider the case that $\widehat{\bm{x}}$ is at a fixed point with the
corresponding values of $\widehat{\theta}_{\mu}^{\setminus\mu}$ for
$\mu\in\bm{M}$. When $\widehat{\bm{x}}$ is slightly deviated from the value at
the fixed point, the deviation propagates to
$\widehat{\theta}_{\mu}^{\setminus\mu}$ by GAMP, and back to
$\widehat{\bm{x}}$. If estimates under the small deviance are moved away from
the fixed point by the GAMP update, the fixed point can be regarded as locally
unstable.

We consider the difference in $\widehat{\theta}_{\mu}^{\setminus\mu}$ induced
by the small deviance in $\widehat{\bm{x}}$ as

$\displaystyle\delta\widehat{\theta}_{\mu}^{\setminus\mu}\equiv\widehat{\theta}_{\mu}^{\setminus\mu}(\\{\widehat{x}_{i}+\delta\widehat{x}_{i}\\})-\widehat{\theta}_{\mu}^{\setminus\mu}(\\{\widehat{x}_{i}\\})$
$\displaystyle=\sum_{i=1}^{N}\frac{\partial\widehat{\theta}_{\mu}^{\setminus\mu}}{\partial\widehat{x}_{i}}\delta\widehat{x}_{i}$
(4.46) $\displaystyle=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}F_{\mu
i}\delta\widehat{x}_{i},$

where $\widehat{\theta}_{\mu}^{\setminus\mu}(\widehat{\bm{x}})$ denotes the
value of $\widehat{\theta}_{\mu}^{\setminus\mu}$ under $\widehat{\bm{x}}$. We
restrict our discussion to the predictors whose entries are i.i.d. with mean 0
and variance 1. For such predictors, the summation of L.H.S. in (4.46) is
reduced to zero at sufficiently large $N$, hence we consider the squared sum
as

$\displaystyle{\delta\widehat{\theta}_{\mu}^{\setminus\mu}}^{2}\simeq\frac{1}{N}\sum_{i=1}^{N}\delta\widehat{x}_{i}^{2}.$
(4.47)

Meanwhile, we consider the deviation of $\widehat{\bm{x}}$ attributed to the
slight deviance of $\theta_{\mu}^{\setminus\mu}$ defined by

$\displaystyle\delta\widehat{x}_{i}\equiv\widehat{x}_{i}(\\{\widehat{\theta}_{\mu}^{\setminus\mu}+\delta\widehat{\theta}_{\mu}^{\setminus\mu}\\})-\widehat{x}_{i}(\\{\widehat{\theta}_{\mu}^{\setminus\mu}\\})=\sum_{\mu=1}^{M}\frac{\partial\widehat{x}_{i}}{\partial\widehat{\theta}_{\mu}^{\setminus\mu}}\delta\widehat{\theta}_{\mu}^{\setminus\mu},$
(4.48)

where $\widehat{x}_{i}(\\{\widehat{\theta}_{\mu}^{\setminus\mu}\\})$ denotes
the estimates under $\\{\widehat{\theta}_{\mu}^{\setminus\mu}\\}$. The
deviation of $\widehat{\theta}_{\mu}^{\setminus\mu}$ propagates to
$\widehat{x}_{i}$ through $\mathrm{m}_{i}$, and hence, we obtain

$\displaystyle\frac{\partial\widehat{x}_{i}}{\partial\widehat{\theta}_{\mu}^{\setminus\mu}}$
$\displaystyle=\frac{\partial\mathrm{m}_{i}}{\partial\widehat{\theta}_{\mu}^{\setminus\mu}}\frac{\partial\widehat{x}_{i}}{\partial\mathrm{m}_{i}}$
(4.49) $\displaystyle=\Sigma_{i}\left\\{\frac{1}{\sqrt{N}}F_{\mu
i}(\check{\partial}_{\hat{\theta}}{g}_{\mathrm{out}})_{\mu}+O(N^{-1})\right\\}s_{i}\Sigma_{i}^{-1},$

where we use the relationship (4.13). For the random preditors, the summation
of (4.49) with respect to $\mu$ is expected to converge to zero for a
sufficiently large $M$, and hence, we consider the squared summation under the
assumptions as

$\displaystyle\delta\widehat{x}_{i}^{2}\simeq\frac{s_{i}^{2}}{N}\sum_{\mu=1}^{M}F_{\mu
i}^{2}(\check{\partial}_{\hat{\theta}}{g}_{\mathrm{out}})_{\mu}^{2}{\delta\widehat{\theta}_{\mu}^{\setminus\mu}}^{2}\simeq\frac{\alpha
s_{i}^{2}}{M}\sum_{\nu=1}^{M}(\check{\partial}_{\hat{\theta}}{g}_{\mathrm{out}})_{\nu}^{2}{\delta\widehat{\theta}_{\mu}^{\setminus\mu}}^{2}.$
(4.50)

Therefore, setting
$\overline{\delta\widehat{x}^{2}}\equiv\frac{1}{N}\sum_{i=1}^{N}\delta\widehat{x}_{i}^{2}$,
we obtain the recursive relationship of $\overline{\delta\widehat{x}^{2}}$ as

$\displaystyle\overline{\delta\widehat{x}^{2}}=\alpha\overline{s^{2}}\overline{(\check{\partial}{g}_{\mathrm{out}})^{2}}\overline{\delta\widehat{x}^{2}},$
(4.51)

where we set $\overline{s^{2}}\equiv\frac{1}{N}\sum_{i=1}^{N}s_{i}^{2}$ and
$\overline{(\check{\partial}{g}_{\mathrm{out}})^{2}}\equiv\frac{1}{M}\sum_{\nu=1}^{M}(\check{\partial}_{\hat{\theta}}{g}_{\mathrm{out}})_{\nu}^{2}$.
Thus, the solution $\overline{\delta\widehat{x}^{2}}=0$ is unstable when

$\displaystyle\alpha\overline{s^{2}}\overline{(\check{\partial}{g}_{\mathrm{out}})^{2}}>1,$
(4.52)

which indicates the local instability of the fixed point. This condition
corresponds the instability of the replica symmetric ansatz in the replica
analysis shown in sec.6.4.

###  4.4 Expression of GDF by GAMP

We start with the form (4.24) for the estimated parameter
$\widehat{\theta}_{\mu}$. By definition,
$\widehat{\theta}_{\mu}^{\setminus\mu}$ and $s_{\theta_{\mu}}^{\setminus\mu}$
are defined in the $\mu$-cavity system, and hence, their values are
independent from $y_{\mu}$ when the graphical model does not have any loops
including $y_{\mu}$. In fact, when the graph contains loops,
$\widehat{\theta}_{\mu}^{\setminus\mu}$ and $s_{\theta_{\mu}}^{\setminus\mu}$
are dependent on $y_{\mu}$, but we assume that their dependency on $y_{\mu}$
is negligible. Under this consideration, the dependency on $y_{\mu}$ only
appears in the Onsager reaction term, i.e. the term $g_{\mathrm{out}}$, and
hence,

$\displaystyle\frac{\partial\widehat{\theta}_{\mu}}{\partial
y_{\mu}}={s_{\theta}}_{\mu}\frac{\partial(\check{g}_{\mathrm{out}})_{\mu}}{\partial
y_{\mu}}.$ (4.53)

From (4.18), in the case of GLM,

$\displaystyle\frac{\partial(\check{g}_{\mathrm{out}})_{\mu}}{\partial
y_{\mu}}$
$\displaystyle=\frac{\beta}{{s_{\theta}}_{\mu}}\left\\{\left\langle(\theta_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu})\\!\left(\theta_{\mu}+b^{\prime}(y_{\mu})\right)\right\rangle_{\theta_{\mu}}\\!\\!\\!\\!-\\!\left\langle\theta_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu}\right\rangle_{\theta_{\mu}}\\!\\!\\!\\!\left\langle\theta_{\mu}+b^{\prime}(y_{\mu})\right\rangle_{\theta_{\mu}}\\!\right\\}$
(4.54) $\displaystyle=\frac{{\chi_{\theta}}_{\mu}}{{s_{\theta}}_{\mu}},$

where we define the variance of the model parameter as

$\displaystyle{\chi_{\theta}}_{\mu}=\beta\left(\langle\theta_{\mu}^{2}\rangle_{\theta_{\mu}}-\langle\theta_{\mu}\rangle_{\theta_{\mu}}^{2}\right).$
(4.55)

Hence, we obtain

$\displaystyle\frac{\partial\widehat{\theta}_{\mu}}{\partial
y_{\mu}}={\chi_{\theta}}_{\mu},$ (4.56)

which holds for the general value of $\beta$, and implies that the GAMP
evaluates GDF and the generalization gap for the in-sample error as [36]

$\displaystyle\mathrm{gdf}^{(\mathrm{AMP})}({\cal
D})=\frac{1}{M}\sum_{\mu=1}^{M}{\chi_{\theta}}_{\mu}$ (4.57)

Eq. (4.57) indicates the proportionality between the response of the learning
error and the variance of the estimated parameter at equilibrium, which can be
regarded as a fluctuation-response relationship described by GAMP [20].

The form of ${\chi_{\theta}}_{\mu}$ can be derived at the $\beta\to\infty$
limit, when the following second order approximation is valid. Expanding the
function
$f_{\xi}(\theta;y_{\mu},\widehat{\theta}_{\mu}^{\setminus\mu},{s_{\theta}}_{\mu})$
around ${\theta}_{\mu}^{*}$ up to the second order, we obtain

$\displaystyle
f_{\xi}(\theta_{\mu};y_{\mu},\widehat{\theta}_{\mu}^{\setminus\mu},{s_{\theta}}_{\mu})$
$\displaystyle=f_{\xi}(\theta^{*}_{\mu};y_{\mu},\widehat{\theta}_{\mu}^{\setminus\mu},{s_{\theta}}_{\mu})$
(4.58)
$\displaystyle+\frac{(\theta_{\mu}-\theta_{\mu}^{*})^{2}}{2}\frac{\partial^{2}}{\partial\theta^{2}}f_{\xi}(\theta;y_{\mu},\widehat{\theta}_{\mu}^{\setminus\mu},{s_{\theta}}_{\mu})\Big{|}_{\theta=\theta_{\mu}^{*}}.$

Hence, in case of GLM, we obtain [53]

$\displaystyle{\chi_{\theta}}_{\mu}=\frac{{s_{\theta}}_{\mu}}{1+{s_{\theta}}_{\mu}a^{\prime\prime}(\theta_{\mu}^{*})}.$
(4.59)

For the Gaussian likelihood, where $a^{\prime\prime}(\theta_{\mu}^{*})=1$, GDF
is given by

$\displaystyle\mathrm{gdf}^{(\mathrm{AMP})}({\cal
D})=\frac{1}{M}\sum_{\mu=1}^{M}\frac{{s_{\theta}}_{\mu}}{1+{s_{\theta}}_{\mu}}.$
(4.60)

Considering the correspondence between the linear estimation rule and GAMP
(4.45), the expression of GDF by GAMP can be regarded as DF using the extended
hat matrix (4.45) for the Gaussian likelihood.

###  4.5 Expression of the functional variance by GAMP

From the relationship (2.59), FV is expressed using the effective distribution
by GAMP (4.20) as

$\displaystyle\mathrm{FV}^{(\mathrm{AMP})}({\cal D})$
$\displaystyle=\frac{\beta}{M}\sum_{\mu=1}^{M}\left\\{\left\langle\left(\ln
f(y_{\mu}|\theta_{\mu})\right)^{2}\right\rangle_{\theta_{\mu}}-\left\langle\ln
f(y_{\mu}|\theta_{\mu})\right\rangle_{\theta_{\mu}}^{2}\right\\}$ (4.61)
$\displaystyle=\frac{1}{\beta
M}\sum_{\mu=1}^{M}\frac{\partial^{2}}{\partial\gamma^{2}}\ln\int
d\theta_{\mu}f^{\gamma\beta}(y_{\mu}|\theta_{\mu})\psi_{c}^{\beta}(\theta_{\mu}|{s_{\theta}}_{\mu},\widehat{\theta}_{\mu}^{\setminus\mu})\Big{|}_{\gamma=1}$

where $\gamma$ is an auxiliary variable. At $\beta\to\infty$, FV is given by

$\displaystyle\mathrm{FV}^{(\mathrm{AMP})}({\cal
D})=\frac{1}{M}\sum_{\mu=1}^{M}\frac{\partial^{2}}{\partial\gamma^{2}}{\cal
T}_{\gamma}(\gamma)\Big{|}_{\gamma=1}$ (4.62)

where

$\displaystyle{\cal T}_{\gamma}(\gamma)$ $\displaystyle=\gamma\ln
f(y_{\mu}|\theta_{\mu}^{*}(\gamma))-\frac{(\theta_{\mu}^{*}(\gamma)-\widehat{\theta}_{\mu})^{2}}{2{{s_{\theta}}_{\mu}}}$
(4.63) $\displaystyle\theta_{\mu}^{*}(\gamma)$
$\displaystyle=\mathop{\mathrm{argmax}}_{\theta}\left\\{\gamma\ln
f(y_{\mu}|\theta)-\frac{(\theta-\widehat{\theta}_{\mu}^{\setminus\mu})^{2}}{2{s_{\theta}}_{\mu}}\right\\}.$
(4.64)

In the case of GLM, the following relationship holds for any $\gamma$ and
$\mu\in\bm{M}$:

$\displaystyle\gamma(y_{\mu}-a^{\prime}(\theta_{\mu}^{*}(\gamma)))-\frac{\theta_{\mu}^{*}(\gamma)-\widehat{\theta}_{\mu}^{\setminus\mu}}{{s_{\theta}}_{\mu}}=0,$
(4.65)

which is reduced to (4.32) at $\gamma=1$, and we denote
$\theta_{\mu}^{*}(\gamma=1)=\theta_{\mu}^{*}$. Further, by differentiating
both sides of (4.65) with respect to $\gamma$, we obtain

$\displaystyle y_{\mu}-a^{\prime}(\theta_{\mu}^{*}(\gamma))-\gamma
a^{\prime\prime}(\theta_{\mu}^{*}(\gamma))\frac{\partial\theta_{\mu}^{*}(\gamma)}{\partial\gamma}=\frac{1}{{s_{\theta}}_{\mu}}\frac{\partial\theta_{\mu}^{*}(\gamma)}{\partial\gamma}.$
(4.66)

Hence, at $\gamma=1$,

$\displaystyle\frac{\partial\theta_{\mu}^{*}(\gamma)}{\partial\gamma}\Big{|}_{\gamma=1}=(\check{g}_{\mathrm{out}})_{\mu}\frac{{s_{\theta}}_{\mu}}{1+{s_{\theta}}_{\mu}a^{\prime\prime}(\theta_{\mu}^{*})},$
(4.67)

where we used the relationship (4.38).

Utilizing (4.65), (4.62) is transformed to

$\displaystyle\mathrm{FV}^{(\mathrm{AMP})}({\cal D})$ $\displaystyle\hskip
28.45274pt=2(\check{g}_{\mathrm{out}})_{\mu}\frac{\partial\theta_{\mu}^{*}(\gamma)}{\partial\gamma}\Big{|}_{\gamma=1}-\left(a^{\prime\prime}(\theta_{\mu}^{*})+\frac{1}{{s_{\theta}}_{\mu}}\right)\left(\frac{\partial\theta_{\mu}^{*}(\gamma)}{\partial\gamma}\Big{|}_{\gamma=1}\right)^{2},$
(4.68)

and substituting (4.67) to (4.68), we get

$\displaystyle\mathrm{FV}^{(\mathrm{AMP})}({\cal
D})=\frac{1}{M}\sum_{\mu=1}^{M}(\check{g}_{\mathrm{out}})_{\mu}^{2}\frac{{s_{\theta}}_{\mu}}{1+{s_{\theta}}_{\mu}a^{\prime\prime}(\theta_{\mu}^{*})}.$
(4.69)

Comparing FV with $\sigma^{2}\times\mathrm{gdf}$, where both are estimators
for the generalization gap with respect to the in-sample prediction error, the
strength of the noise $\sigma^{2}$, which is assumed to be known in $C_{p}$
criterion, is replaced with $(\check{g}_{\mathrm{out}})_{\mu}^{2}$ component-
wise, in case of FV. In the actual usage of the $C_{p}$ criterion, an
estimator is used instead of the variance, when its value is not known in
advance [16]. FV can be regarded as the $C_{p}$ criterion defined with the
estimator of the variance. When the mean of the model
$a^{\prime}(\theta_{\mu}^{*})$ is exactly equal to the data $y_{\mu}$, FV is
reduced to zero as $(\check{g}_{\mathrm{out}})_{\mu}^{2}$ is the squared error
between them given in (4.37).

####  4.5.1 Convergence of the series expansion

We consider the convergence condition of the series expansion of (2.52) for
the derivation of WAIC using the expression obtained using GAMP.

######  Theorem 4.

For the Gaussian likelihood, the radius of convergence for the the series
expansion with respect to $\eta$ is given by

$\displaystyle\eta\lim_{n\to\infty}\left|\frac{\sum_{\mu}(\check{g}_{\mathrm{out}})^{2}_{\mu}\left(\displaystyle\frac{{s_{\theta}}_{\mu}}{1+{s_{\theta}}_{\mu}}\right)^{n}}{\sum_{\mu}(\check{g}_{\mathrm{out}})^{2}_{\mu}\left(\displaystyle\frac{{s_{\theta}}_{\mu}}{1+{s_{\theta}}_{\mu}}\right)^{n-1}}\right|<1.$
(4.70)

When the components of the predictors are i.i.d. random variables of mean zero
and variance 1, the radius of the convergence for the series expansion is
given by

$\displaystyle\eta\left|\frac{{s_{\theta}}}{1+{s_{\theta}}}\right|<1.$ (4.71)

###### Proof.

We prove the theorem in an inductive manner. The first-order term of (2.52) is
given by

$\displaystyle\frac{\partial{\cal
T}_{\gamma}(\gamma)}{\partial\gamma}=\frac{1}{M}\sum_{\mu=1}^{M}\ln
f(y_{\mu}|\theta_{\mu}^{*}(\gamma)).$ (4.72)

We start our proof from the second-order term obtained as

$\displaystyle\frac{\partial^{2}{\cal
T}_{\gamma}(\gamma)}{\partial\gamma^{2}}=\frac{1}{M}\sum_{\mu=1}^{M}(\check{g}_{\mathrm{out}}(\gamma))^{2}_{\mu}\frac{{s_{\theta}}_{\mu}}{1+\gamma{s_{\theta}}_{\mu}},$
(4.73)

which corresponds to $\mathrm{FV}^{(\mathrm{AMP})}$ at $\gamma=1$, and we
define

$\displaystyle(\check{g}_{\mathrm{out}}(\gamma))_{\mu}=y_{\mu}-\theta_{\mu}^{*}(\gamma)=\frac{\theta_{\mu}^{*}(\gamma)-\widehat{\theta}_{\mu}^{\setminus\mu}}{\gamma{s_{\theta}}_{\mu}},$
(4.74)

where
$(\check{g}_{\mathrm{out}}(\gamma=1))_{\mu}=(\check{g}_{\mathrm{out}})_{\mu}$
For $k=3$, we obtain

$\displaystyle\frac{\partial^{3}{\cal
T}_{\gamma}(\gamma)}{\partial\gamma^{3}}=-\frac{3}{M}\sum_{\mu=1}^{M}(\check{g}_{\mathrm{out}}(\gamma))^{2}_{\mu}\left(\frac{{s_{\theta}}_{\mu}}{1+\gamma{s_{\theta}}_{\mu}}\right)^{2}.$
(4.75)

For general $k$, we can show that if the following relationship holds

$\displaystyle\frac{\partial^{k}{\cal
T}(\gamma)}{\partial\gamma^{k}}=\frac{k!}{2M}(-1)^{k}\sum_{\mu=1}^{M}(\check{g}_{\mathrm{out}}(\gamma))^{2}_{\mu}\left(\frac{{s_{\theta}}_{\mu}}{1+\gamma{s_{\theta}}_{\mu}}\right)^{k-1},$
(4.76)

then

$\displaystyle\frac{\partial^{k+1}{\cal
T}(\gamma)}{\partial\gamma^{k+1}}=\frac{(k+1)!}{2M}(-1)^{k+1}\sum_{\mu=1}^{M}(\check{g}_{\mathrm{out}}(\gamma))^{2}_{\mu}\left(\frac{{s_{\theta}}_{\mu}}{1+\gamma{s_{\theta}}_{\mu}}\right)^{k}.$
(4.77)

Therefore, (4.76) holds for any $k\geq 2$, and we obtain

$\displaystyle{\cal T}(\eta)=\frac{\eta}{M}\sum_{\mu=1}^{M}\ln
f(y_{\mu}|\theta_{\mu}^{*})+\frac{1}{2M}\sum_{k=2}^{\infty}(-\eta)^{k}\sum_{\mu=1}^{M}(\check{g}_{\mathrm{out}})^{2}_{\mu}\left(\frac{{s_{\theta}}_{\mu}}{1+{s_{\theta}}_{\mu}}\right)^{k-1}.$
(4.78)

This expression indicates that the radius of the convergence is given by
(4.70). For random i.i.d. predictors, ${s_{\theta}}_{\mu}$ is reduced
${s_{\theta}}$ for any $\mu$, and the radius of the convergence is given by
(4.71). ∎

For the Gaussian likelihood with random i.i.d. predictor, the radius of the
convergence of the series expansion is given by GDF, as indicated in (4.71).

##  5 Approaching LOOCV error from FV at finite $\alpha$

###  5.1 Discrepancy between the FV and LOOCV errors in the finite $\alpha$
region

We show the prediction errors and their estimators for the random predictor
$\bm{F}$, whose components are independently and identically generated by
Gaussian distribution with mean 0 and variance 1. In the numerical simulation,
we set the true regression coefficient $\bm{x}^{(0)}\in\mathbb{R}^{N}$ and the
model parameter as $\bm{\theta}^{(0)}=\frac{1}{\sqrt{N}}\bm{F}\bm{x}^{(0)}$.

####  5.1.1 Ridge regression under Gaussian noise

Figure 4: Prediction errors for the ridge regression at $M=200$ for (a)
$\lambda=1,\sigma=1$ and (b) $\lambda=0.01,\sigma=1$.
$\mathrm{err}_{\mathrm{LOOCV}}$, SURE, and WAIC calculated by GAMP are denoted
by $\circ$, $\square$, and $*$, respectively. The value of each point is
averaged over 100-times numerical simulations, and the bar on the point
indicates standard error. The solid and dashed lines represent the analytical
results obtained using the replica method for the extra-sample prediction
error $\mathrm{err}_{\mathrm{pre}}^{(\mathrm{ex})}$ and in-sample prediction
error $\mathrm{err}_{\mathrm{pre}}^{(\mathrm{in})}$, respectively.

Figure 5: Generalization gaps for ridge regression at $M=200$ for (a)
$\lambda=1,\sigma=1$ and (b) $\lambda=0.01,\sigma=1$.
$\Delta_{\mathrm{LOOCV}}$, GDF, and FV calculated by GAMP are denoted by
$\circ$, $\square$, and $*$, respectively. The value of each point is averaged
over 100-times numerical simulations, and the bar on the point indicates
standard error. The solid and dashed lines denote the analytical results
obtained using the replica method for the generalization gap of the extra-
sample prediction error $\Delta^{(\mathrm{ex})}$ and the in-sample prediction
error $\Delta^{(\mathrm{in})}$, respectively.

There is no need to introduce message passing to obtain the estimate
$\widehat{\bm{x}}({\cal D})$ for the ridge regression (2.42) because the ridge
regression can be solved analytically. We apply our method for evaluating the
FV and LOOCV error to ridge regression to clarify its adequacy. In the
numerical simulation for the ridge regression, we independently generate the
components of $\bm{x}^{(0)}$ according to the Gaussian distribution with mean
0 and variance 1, and we generate the data as $y_{\mu}\sim{\cal
N}(\theta_{\mu}^{(0)},\sigma^{2})$. Figs.5 and 5 show prediction errors and
generalization gaps for the ridge regression at $M=200$, where $N$ runs
$[1,3000]$; therefore, the region $\alpha\in[0.667,200]$ is shown. Parameters
for Figs.5 and 5 are set as (a) $\lambda=1$ and $\sigma=1$, and (b)
$\lambda=0.01$ and $\sigma=1$, respectively. The LOOCV error,
$\mathrm{err}_{\mathrm{LOOCV}}$, and corresponding generalization gap
$\Delta_{\mathrm{LOOCV}}$ are computed by (2.9), and SURE, WAIC, and their
corresponding generalization gaps, namely GDF and FV, are calculated by GAMP.
The value of each point is averaged over 100 realization of
$\\{\bm{y},\bm{F},\bm{x}^{(0)}\\}$. Theoretical results obtained using the
replica method (see sec. 6) for the extra-sample and in-sample prediction
errors are shown by solid and dashed lines, respectively. As shown in the
figure, the GDF by GAMP works well for all region of $\alpha$ as an estimator
of the in-sample prediction error. When the value of $\alpha$ is sufficiently
large, all prediction errors and all generalization gaps coincide with each
other; this means all estimators can work as appropriate estimates of the
generalization gap of the extra-sample prediction error. However, for small
$\alpha$, the differences between the extra-sample prediction error and their
estimates increase, and the cusp around $\alpha=1$ for a sufficiently small
$\lambda$, which is known as the double descent phenomena [7], cannot be
described by WAIC and SURE, as shown in Fig.5 (b) and Fig.5 (b).

####  5.1.2 Linear regression with elastic net penalty under Gaussian noise

Figure 6: Prediction errors for linear regression with elastic net penalty at
$M=200$, $N=400$, $\lambda_{2}=0.01$, and $\sigma=1$ for (a) $\rho=0.1$ and
(b) $\rho=1$; $\lambda_{1}$ is controlled to set
$\widetilde{\alpha}=M/||\widehat{\bm{x}}||_{0}$ at each point. The LOOCV error
by the CD algorithm, SURE, and WAIC calculated by GAMP are denoted by $\circ$,
$\square$, and $*$, respectively. The value of each point is averaged over
100-times numerical simulations, and in particular the bar on the point
indicates the standard error. The solid and dashed lines represent the
analytical results using the replica method for the extra-sample and in-sample
prediction errors, respectively.

,

Figure 7: Generalization gaps for the linear regression with elastic net
penalty at $M=200$, $N=400$, $\lambda_{2}=0.01$ and $\sigma=1$ for (a)
$\rho=0.1$ and (b) $\rho=1$, where $\lambda_{1}$ is controlled to set
$\widetilde{\alpha}=M/||\widehat{\bm{x}}||_{0}$ at each point.
$\Delta_{\mathrm{LOOCV}}$ by CD algorithm, GDF, and FV calculated by GAMP are
denoted by $\circ$, $\square$, and $*$, respectively. The value of each point
is averaged over a 100-times numerical simulation, and the bar on the point
indicates the standard error. The solid and dashed lines denote the analytical
results obtained using the replica method for the generalization gap of the
extra-sample and in-sample prediction errors, respectively.

We show the results of GAMP for the linear regression with elastic net penalty
as an example where the estimate cannot be described by the hat matrix. The
estimates have zero components induced by the regularization parameter
$\lambda_{1}$ of (3.8). Hence, we prepare sparse coefficient $\bm{x}^{(0)}$ in
the numerical simulation. We identically and independently generate the
components in $\bm{x}^{(0)}\in\mathbb{R}^{N}$ according to the Gaussian
distribution with mean zero and variance 1, and then, we choose
$N(1-\rho)$-components and set them to zero, where $\rho\in[0,1]$. Using this
sparse coefficient $\bm{x}^{(0)}$, we set
$\bm{\theta}^{(0)}=\frac{1}{\sqrt{N}}\bm{F}\bm{x}^{(0)}$, and the components
of data $\bm{y}$ are generated as $y_{\mu}\sim{\cal
N}(\theta_{\mu}^{(0)},\sigma^{2})$. When $\lambda_{1}>0$, the estimate
$\widehat{\bm{x}}({\cal D})$ has zero components, hence the number of the
parameter is given by $||\bm{x}||_{0}$, where $||\cdot||_{0}$ is the
$\ell_{0}$ norm. Therefore, we set $\tilde{\alpha}=M/||\bm{x}||_{0}$
distinguishing from $\alpha=M/N$, and the parameter regime where the number of
the model parameters is larger than data dimension corresponds to
$\tilde{\alpha}>1$.

For the elastic net regularization, we cannot derive the exact form of the
estimate as with the ridge regression. Therefore, we compare our result with
the LOOCV error calculated by the coordinate descent (CD) algorithm [35]
considered to lead to reasonable solutions. For the details of the CD
algorithm, see D.

Figs. 7 and 7 show the prediction errors and corresponding generalization gaps
at $M=200$, $N=400$, $\sigma=1$, and $\lambda_{2}=0.01$, respectively, for
linear regression with the elastic net regularization under Gaussian noise. In
the figs., the regularization parameter $\lambda_{1}$ runs from $[0,2]$ for
(a) $\rho=0.1$ and $[0,3]$ for (b) $\rho=1$. The solid and dashed lines
represent the analytical result obtained using the replica method for the
extra-sample and in-sample prediction errors, respectively. As with the ridge
regression, the FV tends to underestimate the prediction error and the
generalization gap for the small $\widetilde{\alpha}$ region.

###  5.2 From FV to the LOOCV error

We discuss the discrepancy between FV and $\Delta_{\mathrm{LOOCV}}$. Here, we
start with FV and not GDF because FV corresponds to the first-order
approximation of the generalization gap of LOOCV error, as WIC is equivalent
to TIC. We approach LOOCV error by removing assumptions imposed in the
derivation of TIC, one by one.

####  5.2.1 Differences in metrics

In deriving TIC, we imposed assumption A3, which corresponds to the ignorance
of the difference between $a^{\prime}(\bm{\theta}(\bm{F},\bm{x}({\cal
D}_{\setminus\mu})))$ and $a^{\prime}(\bm{\theta}(\bm{F},\bm{x}({\cal D})))$
in the case of GLM. This assumption is supposed to be inadequate in the region
where $N$ is sufficiently large because differences between the estimates
under the full sample and that under the LOO sample are accumulated by
$N$-components. We consider the correction of the first-order approximation of
the LOOCV error in the framework of GAMP by removing assumption A3. We denote
FV without assumption A3 as cFV(corrected FV), which is defined by

$\displaystyle\mathrm{cFV}({\cal D})=\frac{1}{M}\mathrm{Tr}\left({\cal
J}^{\varphi}(\widehat{\bm{x}}({\cal D}))^{-1}{\cal
I}^{\setminus\cdot}(\\{\widehat{\bm{x}}({\cal D}_{\setminus\mu})\\})\right).$
(5.1)

Using the expression by GAMP, we obtain

$\displaystyle{\cal I}_{ij}^{\setminus\cdot}(\\{\widehat{\bm{x}}({\cal
D}_{\setminus\mu})\\})$ $\displaystyle=\frac{1}{MN}\sum_{\mu=1}^{M}F_{\mu
i}F_{\mu j}(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu}))^{2},$
(5.2)

and hence $\mathrm{cFV}$ is given by

$\displaystyle\mathrm{cFV}^{(\mathrm{GAMP})}$ $\displaystyle=$
$\displaystyle=\frac{1}{M}\sum_{\mu=1}^{M}(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu}))^{2}{\chi_{\theta}}_{\mu}.$
(5.3)

From (4.32), we know the relationship between $\widehat{\theta}_{\mu}$ and
$\widehat{\theta}^{\setminus\mu}_{\mu}$.

Figure 8: Generalization gaps for ridge regression at $M=200$ and $\sigma=1$
for (a) $\lambda=1$ and (b) $\sigma=0.01$. FV ($*$), cFV ($\triangle$), and
cFV+F$\chi_{2}$ ($\triangledown$) by GAMP and $\Delta^{(\mathrm{ex})}$
obtained using the replica method. The values of the each points are averaged
over 100 data samples ${\cal D}$, and the bars on the points represent
standard error.

Figure 9: Generalization gaps for the linear regression with elastic net
penalty at $M=200$, $N=400$, $\sigma=1$, and $\lambda_{2}=0.01$ for (a)
$\rho=0.1$ and (b) $\rho=1$. The values of $\lambda_{1}$ are controlled to
satisfy $\widetilde{\alpha}=M/||\widehat{\bm{x}}||_{0}$, and FV ($*$), cFV
($\triangle$), and cFV+F$\chi_{2}$ ($\triangledown$) by GAMP and
$\Delta^{(\mathrm{ex})}$ by replica method are also shown. The values of each
points are averaged over 100 data samples ${\cal D}$, and the bars on the
points denote the standard error.

In Figs. 8 and 9, the behavior of the cFV is shown for ridge regression and
linear regression with the elastic net penalty, respectively. The setting of
the numerical simulation for Figs. 8 and 9 are the same as that explained in
the previous subsection. cFV overestimates $\Delta_{\mathrm{LOOCV}}$ for a
sufficiently small $\alpha$; however, the cFV qualitatively approaches the
behavior of $\Delta_{\mathrm{LOOCV}}$; i.e. cFV has a cusp around $\alpha=1$
in Fig.8(b), and the minimum at $\alpha<1$ in Figure 9(b), as with
$\Delta^{(\mathrm{ex})}$.

####  5.2.2 Higher order expansion of cross validation error

We relax assumption A1 for further correction, and consider higher-order terms
with respect to $\widehat{x}_{i}({\cal
D}_{\setminus\mu})-\widehat{x}_{i}({\cal D})$. When we remove the assumption
A1, the expression (2.28) and (2.25) needs to be corrected. For example, when
we consider the second order of $\widehat{x}_{i}({\cal
D}_{\setminus\mu})-\widehat{x}_{i}({\cal D})$, its value should be determined
to satisfy

$\displaystyle\frac{1}{M}\sum_{k=1}^{N}\frac{\partial^{2}}{\partial
x_{k}\partial x_{j}}\varphi({\cal
D},\bm{x})\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})}(\widehat{x}_{k}({\cal
D}_{\setminus\nu})-\widehat{x}_{k}({\cal D}))$
$\displaystyle+\frac{1}{2M}\sum_{k\ell}\frac{\partial^{3}}{\partial
x_{\ell}\partial x_{k}\partial x_{j}}\varphi({\cal
D},\bm{x})\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})}(\widehat{x}_{k}({\cal
D}_{\setminus\nu})-\widehat{x}_{k}({\cal D}))(\widehat{x}_{\ell}({\cal
D}_{\setminus\nu})-\widehat{x}_{\ell}({\cal D}))$
$\displaystyle=\frac{1}{M}\frac{\partial}{\partial x_{j}}\ln
f(y_{\nu}|\theta_{\nu}(\mathrm{f}_{\nu},\bm{x})),$ (5.4)

as shown in A.1, and the value differs from that given by (2.28) in general.

For simplicity, we impose an assumption in the penalty term as

A4

Higher-order derivative of the penalty term $h(x)$ than the second order is
zero for any $x$.

Penalties considered in this paper, i.e. elastic net penalty, satisfies the
condition A4 for any regularization parameter $\lambda_{1}$ and $\lambda_{2}$.
Higher order derivatives of $\varphi$ than the second order is ignorable when
we focus on the Gaussian likelihood with penalty satisfying A4, and hence,
(5.4) is reduced to (2.28). Therefore, when removing the assumption A1 under
the assumpgion A4, the correction of (2.28) does not need to be considered,
and only that of (2.25) needs to be considered. By continuing the expansion of
(2.25) up to the second order, and substituting the expression (2.28), we
obtain the second order term of the generalization gap denoted by
$\mathrm{F}\chi_{2}$ as

$\displaystyle\mathrm{F}\chi_{2}$
$\displaystyle=-\frac{1}{2M}\sum_{\mu=1}^{M}a^{\prime\prime}(\widehat{\theta}_{\mu}^{\setminus\mu})\left(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu})\right)^{2}{\chi_{\theta}}_{\mu}^{2},$
(5.5)

where we used the relationship

$\displaystyle\beta\left(\langle x_{i}x_{j}\rangle_{\beta}-\langle
x_{i}\rangle_{\beta}\langle
x_{j}\rangle_{\beta}\right)=\frac{1}{M}\left({{\cal
J}^{\varphi}}^{-1}\right)_{ij},$ (5.6)

(4.55), and (4.59). In the case of Gaussian likelihood with penalties
satisfying A4, we obtain the GAMP-based estimator of the generalization gap
for $\mathrm{err}_{\mathrm{LOOCV}}$ as

$\displaystyle\widehat{\Delta}_{\mathrm{LOOCV}}=\mathrm{cFV}+\mathrm{F}\chi_{2}.$
(5.7)

Similarly, the $k$-th order term denoted by $\mathrm{F}\chi_{k}$ under the
assumptions A2 and A4 is given by

$\displaystyle\mathrm{F}\chi_{k}$
$\displaystyle=-\frac{1}{k!M}\sum_{\mu=1}^{M}a^{\prime(k)}(\widehat{\theta}_{\mu}^{\setminus\mu})\left(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu})\right)^{k}{\chi_{\theta}}_{\mu}^{k},$
(5.8)

where $a^{\prime(k)}(\widehat{\theta}_{\mu}^{\setminus\mu})$ denotes the
$k$-th order derivative. For the Gaussian likelihood, $\mathrm{F}\chi_{k}$ for
$k\geq 3$ are zero as $a^{\prime(k)}(\theta)$ for $k\geq 3$ are zero for any
$\theta$.

The behaviour of $\widehat{\Delta}_{\mathrm{LOOCV}}$ (5.7) by GAMP is shown in
Figs.8 and 9 for the ridge regression and linear regression with elastic net
penalty, respectively. As the second-order term is negative,
$\widehat{\Delta}_{\mathrm{LOOCV}}$ decreases from cFV and well matches with
the generalization gap for the extra-sample prediction error.

In the discussion so far, we did not mention assumption A2. For Gaussian
likelihood, ${\cal J}(\bm{x})$ do not depend on the argument, and therefore,
we do not need to consider the adequacy of assumption A2.

####  5.2.3 Correction of the Hessian

We consider the removal of assumption A1 and A2 for general likelihoods of
GLM. For the case, correction of (2.28) and (5.8) are required. As shown in
A.1, (2.28) is corrected by removing assumptions A1 and A2 as follows:

$\displaystyle\bm{x}({\cal D}_{\setminus\nu})-\bm{x}({\cal D})=({\cal
K}^{\setminus\nu})^{-1}\frac{1}{M}\frac{\partial}{\partial\bm{x}}\ln
f(y_{\nu}|\theta_{\nu}(\mathbf{f}_{\nu},\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})},$ (5.9)

where ${\cal K}^{\setminus\mu}$ for $\mu\in\bm{M}$ is defined by

$\displaystyle{\cal K}^{\setminus\nu}=-{\cal
J}+\frac{d_{\nu}\mathbf{f}_{\nu}^{\top}\mathbf{f}_{\nu}}{NM}$ (5.10)

with

$\displaystyle
d_{\nu}=\frac{a^{\prime}(\widehat{\theta}_{\nu})-a^{\prime}(\widehat{\theta}_{\nu}^{\setminus\nu})}{\widehat{\theta}_{\nu}-\widehat{\theta}_{\nu}^{\setminus\nu}}-a^{\prime\prime}(\widehat{\theta}_{\nu}).$
(5.11)

When the difference between $\widehat{\theta}_{\nu}$ and
$\widehat{\theta}_{\nu}^{\setminus\nu}$ is sufficiently small, $d_{\nu}$
corresponds to
$a^{\prime\prime}(\widehat{\theta}_{\nu}^{\setminus\nu})-a^{\prime\prime}(\widehat{\theta}_{\nu})$.
Then, for Gaussian likelihood where $a^{\prime\prime}(\theta)=1$ for any
$\theta$, ${\cal K}^{\setminus\mu}={\cal J}$ holds for any $\mu\in\bm{M}$.
Using the matrix ${\cal K}^{\setminus\mu}$, the $k$-th order term of the
generalization gap is corrected as

$\displaystyle\mathrm{F}\chi_{k}=-\frac{1}{k!M}\sum_{\mu=1}^{M}a^{\prime(k)}(\widehat{\theta}_{\mu}^{\setminus\mu})\left(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu})\right)^{k}\left({\chi_{\theta}}_{\mu}^{\setminus\mu}\right)^{k},$
(5.12)

where we define

$\displaystyle{\chi_{\theta}}_{\mu}^{\setminus\mu}=\frac{\beta}{N}\sum_{i,j}F_{\mu
i}F_{\mu j}({{\cal K}^{\setminus\mu}}^{-1})_{ij}.$ (5.13)

As explained in B.4, the $k$-th derivative of $a(\theta)$ corresponds to the
$k$-th cumulant of the distribution $f(y|\theta)$ in case of the exponential
family distribution. We denote the $k$-th cumulant of the distribution
$f(\cdot|\theta)$ as $\kappa_{k}[f(\cdot|\theta)]$, and obtain the following
expression:

$\displaystyle\sum_{k=2}^{\infty}\mathrm{F}\chi_{k}$
$\displaystyle=-\frac{1}{M}\sum_{\mu=1}^{M}\sum_{k=2}^{\infty}\frac{1}{k!}\kappa_{k}[f(\cdot|\widehat{\theta}_{\mu}^{\setminus\mu})]\left\\{{\chi_{\theta}}_{\mu}^{\setminus\mu}\left(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu})\right)\right\\}^{k}$
(5.14) $\displaystyle=-\frac{1}{M}\sum_{\mu=1}^{M}{\cal
C}_{f(\cdot|\widehat{\theta}_{\mu}^{\setminus\mu})}\left({\chi_{\theta}}_{\mu}^{\setminus\mu}(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu}))\right)$
$\displaystyle\hskip
42.67912pt+\frac{1}{M}\sum_{\mu=1}^{M}a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu}){\chi_{\theta}}_{\mu}^{\setminus\mu}\left(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu})\right),$

where ${\cal C}_{f(\cdot|\theta)}(t)$ is the cumulant generating function;

$\displaystyle{\cal C}_{f(\cdot|\theta)}(t)$ $\displaystyle=\ln\int
dy\exp(ty)f(y|\theta)=a(\theta+t)-a(\theta).$ (5.15)

In summarize, we obtain the following form of the estimator of the
generalization gap for LOOCV error:

$\displaystyle\widehat{\Delta}_{\mathrm{LOOCV}}=\mathrm{F}\chi\mathrm{s},$
(5.16)

where

$\displaystyle\mathrm{F}\chi\mathrm{s}$
$\displaystyle=\mathrm{cFV}-\frac{1}{M}\sum_{\mu=1}^{M}\Big{\\{}a\left(\widehat{\theta}_{\mu}^{\setminus\mu}+{\chi_{\theta}}_{\mu}^{\setminus\mu}(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu}))\right)$
(5.17) $\displaystyle\hskip
85.35826pt-a(\widehat{\theta}_{\mu}^{\setminus\mu})-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu}){\chi_{\theta}}_{\mu}^{\setminus\mu}\left(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu})\right)\Big{\\}}.$

The remaining task is the computation of ${\cal K}^{\setminus\mu}$ and
associated ${{\chi}_{\theta}}_{\mu}^{\setminus\mu}$ for $\mu\in\bm{M}$. We
skip this $M$-times computation using the Sherman–Morrison formula (2.16).
Using the formula, we obtain

$\displaystyle{\chi_{\theta}}_{\mu}^{\setminus\mu}=\frac{{\chi_{\theta}}_{\mu}}{1+d_{\mu}{\chi_{\theta}}_{\mu}},$
(5.18)

hence we do not need the computation of the inverse of matrix for each
$\mu\in\bm{M}$.

Figure 10: (a) Prediction errors and (b) generalization gaps for logistic
regression with $\ell_{2}$ penalty at $N=200$, $\lambda=0.01$, and $\sigma=1$.
LOOCV error ($\mathrm{err}_{\mathrm{LOOCV}}$) calculated by CD algorithm,
SURE, WAIC, and $\mathrm{err}_{\mathrm{train}}+\mathrm{F}\chi\mathrm{s}$
calculated by GAMP are denoted by $\bigcirc$, $\square$, $*$ and
$\triangledown$ in (a), respectively, and corresponding generalization gaps
are shown in (b) with the same markers. The solid line denotes the
analytically derived results obtained using the replica method.

In Fig. 10, the GAMP-based estimators of (a) prediction error and (b)
generalization gap are shown for the logistic regression with the $\ell_{2}$
penalty at $N=200$ and $\lambda=0.01$, where $a^{\prime(n)}(\theta)$ for
$n\geq 3$ have finite values. We set the true value of the regression
coefficient $\bm{x}^{(0)}$ and generate the components according to the
Gaussian distribution with mean 0 and variance 1. The data
$\bm{y}\in[0,1]^{M}$ is generated as
$y_{\mu}=\mathbb{I}(\theta_{\mu}^{(0)}+\varepsilon_{\mu})$, where
$\bm{\theta}^{(0)}=\frac{1}{\sqrt{N}}\bm{F}\bm{x}^{(0)}$ and
$\varepsilon_{\mu}\sim{\cal N}(0,\sigma^{2})$. For comparison, we show the
estimates obtained by CD algorithm denoted by $\bigcirc$, and the analytical
results of the extra-sample prediction error obtained using the replica method
by solid line. CD algorithm for the GLM including logistic regression is
explained in D.2. As shown in Fig.10, the GAMP-based generalization gap given
by (5.16) shows good agreement with the results by CD algorithm and replica
method for any $\alpha$.

We summarize the procedure for obtaining the estimator for the generalization
gap of the LOOCV error in Algorithm 2, where the notation $\mathrm{GAMP}({\cal
D})$ in the algorithm denotes to obtain the estimates under given data ${\cal
D}$ by running the GAMP of Algorithm 1.

Algorithm 2 Generalization gap for CV by GAMP

1:${\cal D}=\\{\bm{y}$, $\bm{F}\\}$, and
$\\{\widehat{\bm{x}},\widehat{\bm{\theta}},\bm{s}_{\theta}\\}\leftarrow\mathrm{GAMP}({\cal
D})$

2:$\widehat{\Delta}_{\mathrm{LOOCV}}$

3:for $\mu=1,\ldots,M$ do

4:
$\widehat{\theta}_{\mu}^{\setminus\mu}\leftarrow\widehat{\theta}_{\mu}-{s_{\theta}}_{\mu}(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}))$

5:
${\chi_{\theta}}_{\mu}\leftarrow\displaystyle\frac{{s_{\theta}}_{\mu}}{1+{s_{\theta}}_{\mu}a^{\prime\prime}(\widehat{\theta}_{\mu})}$

6:
$d_{\mu}\leftarrow\displaystyle\frac{a^{\prime}(\widehat{\theta}_{\mu})-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu})}{\widehat{\theta}_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu}}-a^{\prime\prime}(\widehat{\theta}_{\mu})$

7:
${\chi_{\theta}}_{\mu}^{\setminus\mu}\leftarrow\displaystyle\frac{{\chi_{\theta}}_{\mu}}{1+d_{\mu}{\chi_{\theta}}_{\mu}}$

8:end for

9:$\mathrm{cFV}\leftarrow\displaystyle\frac{1}{M}\sum_{\mu=1}^{M}(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu}))^{2}{\chi_{\theta}}_{\mu}^{\setminus\mu}$

10:$\mathrm{F}\chi\leftarrow-\displaystyle\frac{1}{M}\sum_{\mu=1}^{M}\left\\{a\left(\widehat{\theta}_{\mu}^{\setminus\mu}+(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu})){\chi_{\theta}}_{\mu}^{\setminus\mu}\right)-a(\widehat{\theta}_{\mu}^{\setminus\mu})-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu})(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu})){\chi_{\theta}}_{\mu}^{\setminus\mu}\right\\}$

11:$\widehat{\Delta}_{\mathrm{LOOCV}}\leftarrow\mathrm{cFV}+\mathrm{F}\chi$

####  5.2.4 Calculation of $\mathrm{cFV}$ and $\mathrm{F}\chi\mathrm{s}$
without GAMP

The problem of GAMP is that it does not converge when the stability condition
(4.52) is not satisfied. In the parameter region where GAMP does not converge,
there could be a local minima of the exponential order with respect to $N$, as
indicated by the correspondence with the de Almeida–Thouless instability [71].
However, even in such cases, one can obtain the GAMP-inspired estimator
without using GAMP at each minimum obtained by the arbitrary algorithm that
can converge to a local minimum. The characteristic variables of GAMP are
$s_{i}$, ${s_{\theta}}_{\mu}$, and
$\partial_{\widehat{\theta}}g_{\mathrm{out}}$, in the sense that general
algorithms do not contain these variables. If we obtain $s_{i}$ numerically,
the other variables ${s_{\theta}}_{\mu}$, ${\chi_{\theta}}_{\mu}$, and
${\chi_{\theta}}_{\mu}^{\setminus\mu}$ can be obtained, and further,
$\partial_{\widehat{\theta}}g_{\mathrm{out}}$ can also be obtained. Thus, to
mimic the calculation of $s_{i}$ is important to extend the GAMP-based
estimator for other algorithms. This calculation of $s_{i}$ is resolved using
the relationship (5.6).

Algorithm 3 GAMP-inspired generalization gap without GAMP

1:${\cal D}=\\{\bm{y}$, $\bm{F}\\}$, and
$\\{\widehat{\bm{x}},\widehat{\bm{\theta}}\\}\leftarrow\mathrm{arbitrary~{}algorithm}({\cal
D})$

2:$\widehat{\Delta}_{\mathrm{LOOCV}}$

3:$\bm{D}\leftarrow\mathrm{diag}(a^{\prime\prime}(\widehat{\bm{\theta}}))$

4:$\bm{\chi}\leftarrow\displaystyle\left(\frac{1}{N}\bm{F}^{\top}\bm{D}\bm{F}+\frac{\partial^{2}}{\partial\bm{x}\partial\bm{x}^{\top}}h(\bm{x})\right)^{-1}$
$\triangleright$ For sparse estimators, use $\bm{\chi}^{{\cal L}}$.

5:for $i=1,...,N$ do

6: $s_{i}\leftarrow\chi_{ii}$ $\triangleright$ For sparse estimators,
$s_{i}\leftarrow\chi^{\cal L}_{ii}$ for $i\in{\cal L}$, otherwise
$s_{i}\leftarrow 0$.

7:end for

8:for $\mu=1,\ldots,M$ do

9: ${s_{\theta}}_{\mu}\leftarrow\displaystyle\frac{1}{N}\sum_{i=1}^{N}F_{\mu
i}^{2}s_{i}$

10:
${\chi_{\theta}}_{\mu}\leftarrow\displaystyle\frac{{s_{\theta}}_{\mu}}{1+{s_{\theta}}_{\mu}a^{\prime\prime}(\theta_{\mu})}$

11:
$\widehat{\theta}_{\mu}^{\setminus\mu}\leftarrow\widehat{\theta}_{\mu}-{s_{\theta}}_{\mu}(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}))$

12:
$d_{\mu}\leftarrow\displaystyle\frac{a^{\prime}(\widehat{\theta}_{\mu})-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu})}{\widehat{\theta}_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu}}-a^{\prime\prime}(\widehat{\theta}_{\mu})$

13:
${\chi_{\theta}}_{\mu}^{\setminus\mu}\leftarrow\displaystyle\frac{{\chi_{\theta}}_{\mu}}{1+d_{\mu}{\chi_{\theta}}_{\mu}}$

14:end for

15:$\mathrm{cFV}\leftarrow\displaystyle\frac{1}{M}\sum_{\mu=1}^{M}(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu}))^{2}{\chi_{\theta}}_{\mu}^{\setminus\mu}$

16:$\mathrm{F}\chi\leftarrow-\displaystyle\frac{1}{M}\sum_{\mu=1}^{M}\left\\{a\left(\widehat{\theta}_{\mu}^{\setminus\mu}+(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu})){\chi_{\theta}}_{\mu}^{\setminus\mu}\right)-a(\widehat{\theta}_{\mu}^{\setminus\mu})-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu})(y_{\mu}-a^{\prime}(\widehat{\theta}_{\mu}^{\setminus\mu})){\chi_{\theta}}_{\mu}^{\setminus\mu}\right\\}$

17:$\widehat{\Delta}_{\mathrm{LOOCV}}\leftarrow\mathrm{cFV}+\mathrm{F}\chi$

In Algorithm 3, we showed the procedure to obtain the GAMP-based estimator
without using GAMP. Here, we consider that the algorithm returns the estimates
$\widehat{\bm{x}}$ and $\widehat{\bm{\theta}}$ corresponding to the input data
${\cal D}$. The calculation of $s_{i}$ is implemented in the 2nd line of the
Algorithm 3, where $\chi_{ii}$ corresponds to $s_{i}$. In the case of sparse
estimators which contains zero components, such as the regression under the
$\ell_{1}$ or the elastic net regularization, it is natural to set $s_{i}=0$
for $i$ with $\widehat{x}_{i}=0$, assuming that the set of zero components do
not change against the small perturbation on data. We denote the support set
as ${\cal L}$, where $x_{i}\neq 0$ for $i\in{\cal L}$; otherwise, $x_{i}=0$.
For the calculation of $s_{i}$ for $i\in{\cal L}$, we define the submatrix of
the predictor $\bm{F}$ as $\bm{F}^{{\cal L}}$ that comprises row vectors in
${\cal L}$ [72]. We define the Hessian matrix for the support set as

$\displaystyle{\cal H}^{\cal L}=\frac{\partial^{2}}{\partial{\bm{x}^{\cal
L}}\partial{\bm{x}^{\cal L}}^{\top}}\left\\{\ln
f\left(\bm{y}\Big{|}\frac{1}{\sqrt{N}}\bm{F}^{{\cal L}}\bm{x}^{\cal
L}\right)-\phi(\bm{x}^{\cal L})\right\\}$ (5.19)

and the inverse $\bm{\chi}^{\cal L}=({\cal H}^{\cal L})^{-1}$; we set

$\displaystyle s_{i}=\left\\{\begin{array}[]{ll}\chi_{ii}^{\cal L}&i\in{\cal
L}\\\ 0&\mathrm{otherwise}.\end{array}\right.$ (5.22)

###  5.3 Correction of cavity-variance for correlated predictors

The correlation between predictors is neglected in GAMP, as shown in
assumption a1 and a2, and therefore, there is no mathematical support for the
accuracy of the estimators for the non-i.i.d. predictors. Here, we consider
two non-i.i.d. predictors and apply the GAMP-based calculation of the
estimators for the generalization gap.

  * •

Correlated Gaussian predictor: The row vectors of the predictor matrix are
generated by ${\cal N}(\bm{0},{\cal V})$ with a covariance matrix ${\cal
V}\in\mathbb{R}^{N\times N}$. We set the component of ${\cal V}$ as ${\cal
V}_{ij}=\sigma_{d}^{-|i-j|}$ with $\sigma_{d}>0$; i.e., the power law decay in
the distance.

  * •

Rank deficient predictor: We generate i.i.d. Gaussian predictor
$\bm{G}\in\mathbb{R}^{M\times N}$ and apply singular value decomposition as
$\bm{G}=\bm{U}\bm{S}\bm{V}^{\top}$, where $\bm{U}\in\mathbb{R}^{M\times r}$
and $\bm{V}\in\mathbb{R}^{r\times N}$ are the set of left and right singular
vectors, respectively, and $r=\min(M,N)$. The matrix
$\bm{S}\in\mathbb{R}^{r\times r}$ is a diagonal matrix whose $(i,i)$ component
is the $i$-th singular value. We choose singular values with fraction
$r\times(1-\rho_{F})$, and set them as 0. We denote the matrix of singular
values with the zero components as $\widetilde{\bm{S}}$, and construct the
rank-deficient predictor matrix by $\bm{F}=\bm{U}\widetilde{\bm{S}}\bm{V}$
with the normalization $\sum_{\mu=1}^{M}F_{\mu i}^{2}=M$ for any $i$.

Under these predictor matrix, we consider the same generating rule of the data
$\bm{y}$ as with the previous sections.

Fig. 11 shows the generalization gaps at $N=200$, $\lambda=0.01$ and
$\sigma=1$ for (a) correlated Gaussian predictors at $\sigma_{d}=0.5$ and (b)
rank-deficient predictors at $\rho_{F}=0.9$ under the ridge regression. As the
likelihoods are Gaussian distribution, $\mathrm{cFV}+\mathrm{F}\chi_{2}$ is
theoretically sufficient for the description of the generalization gap, but it
overestimates the generalization gap $\Delta_{\mathrm{LOOCV}}$ in the small
$\alpha$ region. Fig. 12 shows the generalization gap for the linear
regression elastic net penalty at $N=400$ and $\lambda_{2}=0.01$. As with Fig.
11, the GAMP-based estimator overestimates the actual value of the
$\Delta_{\mathrm{LOOCV}}$.

Figure 11: Generalization gaps under the ridge regression at $M=200$,
$\sigma=1$, and $\lambda=0.01$. $\Delta_{\mathrm{LOOCV}}$ in the linear
algebraic form is denoted by a solid line, and FV ($*$), cFV+F$\chi_{2}$
($\triangledown$), and cFV+cF$\chi_{2}$ with cavity variance ($\diamond$)
calculated by GAMP are shown. (a) and (b) are for the correlated Gaussian
predictor with $\sigma_{d}=0.5$, and rank-deficient predictor with
$\rho_{F}=0.9$, respectively. Each data point is averaged over 100 samples of
${\cal D}$, and the bars on the points represent standard errors.

Figure 12: Generalization gaps under linear regression with elastic net
regularization at $M=200$, $N=400$, $\lambda_{2}=0.01$, $\sigma=1$, and
$\rho=1$. LOOCV obtained using the CD algorithm is denoted by the solid line;
FV ($*$), cFV+F$\chi_{2}$ ($\triangledown$), and cFV+cF$\chi_{2}$ with cavity
variance ($\diamond$) calculated by GAMP are shown. The GAMP-inspired
generalization gap using the estimates given by CD algorithm is denoted by
pentagons. (a) and (b) are for the correlated Gaussian predictor with
$\sigma_{d}=0.5$, and the rank-deficient predictor with $\rho_{F}=0.9$,
respectively, and $\lambda_{1}$ is controlled to set
$\widetilde{\alpha}=M/||\widehat{\bm{x}}||_{0}$. Each data point is averaged
over 100 samples of ${\cal D}$, and bars on the points represent the standard
errors.

A part of the discrepancies between $\widehat{\Delta}_{\mathrm{LOOCV}}$ and
$\Delta_{\mathrm{LOOCV}}$ attributed to the non-i.i.d.-ness in predictor
matrices is improved by considering the differences between the variables
${s_{\theta}}_{\mu}$ and ${s_{\theta}}_{\mu}^{\setminus\mu}$, which we call
the cavity variance, and it cannot be determined in the framework of GAMP. As
shown in (4.25), we ignore their difference as it scales as $O(N^{-1})$ for
the random predictors satisfying assumptions a1 and a2; however, for the
correlated predictors, this is not the case. GAMP do not have closed equations
to determine ${s_{\theta}}^{\setminus\mu}$, and therefore, we heuristically
derive the form of ${s_{\theta}}_{\mu}^{\setminus\mu}$ using the relationship
(5.6).

By definition,

$\displaystyle{s_{\theta}}_{\mu}^{\setminus\mu}=\frac{\beta}{N}\sum_{ij}F_{\mu
i}F_{\mu j}\left(\langle x_{i}x_{j}\rangle_{\beta}^{\setminus\mu}-\langle
x_{i}\rangle_{\beta}^{\setminus\mu}\langle
x_{j}\rangle_{\beta}^{\setminus\mu}\right),$ (5.23)

where $\langle\cdot\rangle_{\beta}^{\setminus\mu}$ denotes the posterior
average under the LOO sample ${\cal D}_{\setminus\mu}$,

$\displaystyle p_{\beta}(\bm{x}|{\cal
D}_{\setminus\mu})\propto\prod_{\nu\in\bm{M}\neq\mu}f^{\beta}(y_{\nu}|\bm{\theta}(\bm{x}))\phi^{\beta}(\bm{x}).$
(5.24)

As with (5.6), we obtain

$\displaystyle\beta\left(\langle
x_{i}x_{j}\rangle^{\setminus\mu}_{\beta}-\langle
x_{i}\rangle^{\setminus\mu}_{\beta}\langle
x_{j}\rangle^{\setminus\mu}_{\beta}\right)=\frac{1}{M}({{\cal
J}^{\varphi\setminus\mu}(\widehat{\bm{x}}({\cal