D}_{\setminus\mu}))}^{-1})_{ij},$ (5.25)

where

$\displaystyle{\cal J}_{ij}^{\varphi\setminus\mu}(\widehat{\bm{x}}({\cal
D}_{\setminus\mu}))=-\frac{1}{M}\frac{\partial^{2}}{\partial x_{i}\partial
x_{j}}\left\\{\sum_{\nu\in\bm{M}\setminus\mu}\ln
f(y_{\nu}|\bm{\theta}(\bm{x}))-h(\bm{x})\right\\}\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\mu})}.$ (5.26)

We avoid the $\mu$-times computation of ${{\cal
J}^{\varphi\setminus\mu}}^{-1}$ for $\mu\in\bm{M}$ by utilizing the
Sherman–Morrison formula (2.16). Setting
$\widetilde{\bm{F}}=\sqrt{\bm{D}(\bm{\theta}^{\setminus\nu})}\bm{F}$, where
$\bm{D}(\bm{\theta}^{\setminus\nu})$ is the diagonal matrix whose $(\mu,\mu)$
component is $\sqrt{a^{\prime\prime}(\theta_{\mu}^{\setminus\nu})}$, we obtain

$\displaystyle{s_{\theta}}_{\nu}^{\setminus\nu}$
$\displaystyle={\chi_{\theta}}_{\nu}+\frac{\frac{1}{N}||\widetilde{\textbf{f}}_{\nu}(\frac{1}{N}\widetilde{\bm{F}}^{\top}\widetilde{\bm{F}}+\lambda\bm{I}_{N})^{-1}\textbf{f}_{\nu}^{\top}||_{2}^{2}}{1-\frac{1}{N}\widetilde{\textbf{f}}_{\nu}(\frac{1}{N}\widetilde{\bm{F}}^{\top}\widetilde{\bm{F}}+\lambda\bm{I}_{N})^{-1}\widetilde{\textbf{f}}_{\nu}^{\top}}$
(5.27)
$\displaystyle=\frac{{\chi_{\theta}}_{\nu}}{1-a^{\prime\prime}(\theta_{\nu}^{\setminus\nu}){\chi_{\theta}}_{\nu}}.$

When we put ${\chi_{\theta}}_{\nu}$ derived GAMP, (4.59), into (5.27), we
obtain the equality ${s_{\theta}}_{\nu}^{\setminus\nu}={s_{\theta}}_{\nu}$,
and hence, the computation of the cavity variance should be based on Algorithm
3.

Figure 13: Generalization gaps under the logistic regression with $\ell_{2}$
regularization at $M=200$, $\lambda=0.01$, and $\sigma=1$. LOOCV obtained by
the CD algorithm is denoted by a solid line; FV ($*$), F$\chi$s
($\triangledown$), and F$\chi$s 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. Each data point is averaged over
100 samples of ${\cal D}$, and the bars on the points represent standard
errors.

The behaviours of the estimators for the generalization gap considering the
cavity variance at $N=200$ are shown in Figs. 11, 12, and 13 for the ridge
regression, linear regression with elastic net penalty, and logistic
regression with $\ell_{2}$ penalty, respectively. In these examples, the
accuracy of the estimators are significantly improved by considering the
cavity variance.

The procedure to obtain the generalization gap considering the cavity variance
is summarised in Algorithm 4.

Algorithm 4 GAMP-inspired generalization gap without GAMP (with cavity
variance)

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}(\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 the sparse estimators, use $\bm{\chi}^{{\cal L}}$.

5:for $i=1,\ldots,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}\chi_{ii}$

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

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

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

13:
$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})$

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

15:end for

16:$\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}$

17:$\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\\}$

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

##  6 Replica method

For a part of predictor matrices, the prediction error and generalization gap
can be analytically derived by using replica method. In this section, we show
the analysis for the random predictor matrix $\bm{F}$ whose components are
independently and identically distributed according to the Gaussian
distribution. In such cases, it is shown that the typical trajectory of the
time evolution of the GAMP, termed as state evolution, corresponds to the
saddle point equations of the replica method [14]. The rigorousness of the
state evolution for GAMP has been proved under i.i.d. Gaussian data [73], and
correlated Gaussian data [51].

In the replica analysis, we derive the expectation value of the errors with
respect to ${\cal D}=\\{\bm{y},\bm{F}\\}$, to understand the typical property
under the given ${\cal D}$ [64]. The basis for the analysis is the free energy
density, which is defined as

$\displaystyle\Psi=-\lim_{\beta\to\infty}\frac{1}{N\beta}\mathbb{E}_{{\cal
D}}[\ln Z_{\beta}({\cal D})].$ (6.1)

The training error relates to the free energy density at zero temperature, and
as discussed in sec.2.6, the expression of FV by the replica method can be
derived from the partition function by introducing auxiliary variables.

In the analysis, we consider the same setting as the numerical simulations by
GAMP summarized in sec.5, where the data $\bm{y}$ is generated by a certain
rule governed by the ‘true’ parameter
$\bm{\theta}^{(0)}=\bm{F}\bm{x}^{(0)}/\sqrt{N}$ with a ‘true’ regression
coefficient $x_{i}^{(0)}$. Here, we consider that the true regression
coefficient is generated by Bernoulli–Gaussian distribution

$\displaystyle\phi^{(0)}(x^{(0)})=(1-\rho)\delta(x^{(0)})+\frac{\rho}{\sqrt{2\pi\sigma_{x}^{2}}}\exp\left(-\frac{(x^{(0)})^{2}}{2\sigma_{x}^{2}}\right),$
(6.2)

where $\rho\in[0,1]$ is the fraction of nonzero components in the true
coefficient. We assume that the model parameter is subjected to additive noise
$\bm{\varepsilon}\in\mathbb{R}^{M}$; hence we denote $\bm{y}={\cal
G}(\bm{\theta}^{(0)}+\bm{\varepsilon})$, where ${\cal G}$ is an arbitrary
function that describe the data generative process. The expectation of the
free energy density with respect to $\bm{y},\bm{\varepsilon},\bm{F}$ and
$\bm{x}^{(0)}$ is implemented by replica method based on the following
identity

$\displaystyle\mathbb{E}[\ln Z_{\beta}({\cal D})]=\lim_{n\to
0}\frac{\partial}{\partial n}\mathbb{E}_{\cal D}[Z_{\beta}^{n}({\cal D})].$
(6.3)

We focus on the $N\to\infty$ and $M\to\infty$ limit which keep $M/N=\alpha\sim
O(1)$. For the sake of convenience, we define free entropy density for the
finite $\beta$ and finite $n$ as

$\displaystyle\Phi_{\beta}(n)=\frac{1}{N}\ln\mathbb{E}_{{\cal
D}}[Z_{\beta}^{n}({\cal D})].$ (6.4)

Assuming $n$ is a positive integer, we express the $n$-th power by introducing
the $n$-replicated system as

$\displaystyle\Phi_{\beta}(n)$ $\displaystyle=\frac{1}{N}\ln\int
d\bm{\varepsilon}P_{\varepsilon}(\bm{\varepsilon})\int
d\bm{F}P_{F}(\bm{F})\int d\bm{x}^{(0)}\phi^{(0)}(\bm{x}^{(0)})$ (6.5)
$\displaystyle\times\prod_{a=1}^{n}\int d\bm{x}^{(a)}f^{\beta}\left({\cal
G}(\bm{\theta}^{(0)}+\bm{\varepsilon})|\bm{\theta}^{(a)}\right)\phi^{\beta}(\bm{x}^{(a)})$

where we set $\bm{\theta}^{(a)}=\frac{1}{\sqrt{N}}\bm{Fx}^{(a)}$. For
implementing the integrals, we introduce the subshells with respect to the
macroscopic variables $\\{q^{(ab)}\\}$ and $\\{m^{(a)}\\}$ by inserting
trivial identities

$\displaystyle 1$ $\displaystyle=\int
dq^{(ab)}\delta\left(q^{(ab)}-\frac{1}{N}\sum_{i=1}^{N}x_{i}^{(a)}x_{i}^{(b)}\right)$
(6.6) $\displaystyle=\int
dq^{(ab)}d\widehat{q}^{(ab)}\exp\left(-Nq^{(ab)}\widehat{q}^{(ab)}+\widehat{q}^{(ab)}\sum_{i=1}^{N}x_{i}^{(a)}x_{i}^{(b)}\right)$
$\displaystyle 1$ $\displaystyle=\int
dm^{(a)}\delta\left(m^{(a)}-\frac{1}{N}\sum_{i=1}^{N}x_{i}^{(0)}x_{i}^{(a)}\right)$
(6.7) $\displaystyle=\int
dm^{(a)}d\widehat{m}^{(a)}\exp\left(-Nm^{(a)}\widehat{m}^{(a)}+\widehat{m}^{(a)}\sum_{i=1}^{N}x_{i}^{(0)}x_{i}^{(a)}\right)$

for all pairs of $(a,b)~{}(a,b\in\bm{n})$ and $a\in\bm{n}$, where
$\bm{n}=\\{1,\ldots,n\\}$, and $\\{\widehat{q}^{(ab)}\\}$ and
$\\{\widehat{m}^{(a)}\\}$ are conjugate variables for the integral expressions
of the delta functions. The expectation with respect to $\bm{F}$, whose
components are independently and identically distributed according to Gaussian
distribution with mean 0 and variance 1, within the subshells lead to

$\displaystyle\mathbb{E}_{\bm{F}|q^{(ab)}}[\theta^{(a)}_{\mu}\theta^{(b)}_{\mu}]$
$\displaystyle=q^{(ab)}$ (6.8)
$\displaystyle\mathbb{E}_{\bm{F}|m^{(a)}}[\theta^{(a)}_{\mu}\theta^{(0)}_{\mu}]$
$\displaystyle=m^{(a)}$ (6.9)
$\displaystyle\mathbb{E}_{\bm{F}}[\theta^{(0)}_{\mu}\theta^{(0)}_{\mu}]$
$\displaystyle=\rho\sigma_{x}^{2}.$ (6.10)

Thus, in the subshell, the joint distribution of
$\widetilde{\bm{\theta}}_{\mu}\equiv\\{\theta_{\mu}^{(0)},\theta_{\mu}^{(1)},\cdots\theta_{\mu}^{(n)}\\}$
is given by

$\displaystyle P_{\theta}(\tilde{\bm{\theta}}_{\mu}|{\cal
Q})=\frac{1}{\sqrt{2\pi|{\cal
Q}|}}\exp\left(-\frac{1}{2}\widetilde{\bm{\theta}}_{\mu}^{\top}{\cal
Q}^{-1}\widetilde{\bm{\theta}}_{\mu}\right)$ (6.11)

for any $\mu\in\bm{M}$, where ${\cal Q}\in\mathbb{R}^{(n+1)\times(n+1)}$ is
given by

$\displaystyle{\cal Q}_{ab}=\left\\{\begin{array}[]{c c}q^{(ab)}&a\neq
0~{}\mbox{and}~{}b\neq 0\\\ m^{(a)}&a\neq 0~{}\mbox{and}~{}b=0\\\
m^{(b)}&b=0~{}\mbox{and}~{}b\neq 0\\\
\rho\sigma_{x}^{2}&a=b=0\end{array}\right..$ (6.16)

Using the joint distribution with respect to $\widetilde{\bm{\theta}}_{\mu}$,
we obtain

$\displaystyle\Phi_{\beta}$ $\displaystyle(n)=\frac{1}{N}\ln\int\prod_{a\leq
b}dq^{(ab)}d\widehat{q}^{(ab)}\prod_{a=1}^{n}dm^{(a)}d\widehat{m}^{(a)}$
$\displaystyle\times\exp\left(-N\sum_{a\leq
b}q^{(ab)}\widehat{q}^{(ab)}-N\sum_{a=1}^{n}m^{(a)}\hat{m}^{(a)}\right)$
$\displaystyle\times\left[\int d\widetilde{\bm{\theta}}\int d\varepsilon
P_{\varepsilon}(\varepsilon)P_{\theta}(\widetilde{\bm{\theta}}|{\cal
Q})\prod_{a=1}^{n}f\left({\cal
G}(\theta^{(0)}+\varepsilon)|\theta^{(a)}\right)\right]^{M}$
$\displaystyle\times\left[\int\\!\\!d\widetilde{\bm{x}}\\!\prod_{a=1}^{n}\\!\phi(x^{(a)})\phi^{(0)}(x^{(0)})\\!\exp\\!\\!\left(\sum_{a\leq
b}\widehat{q}^{(ab)}x^{(a)}x^{(b)}\\!+\\!\\!\sum_{a=1}^{n}\widehat{m}^{(a)}x^{(0)}x^{(a)}\\!\\!\right)\\!\right]^{N},$

where $\widetilde{\bm{x}}=\\{x^{(0)},x^{(1)},\cdots x^{(n)}\\}$ and
$P(\varepsilon)$ represents the distribution of the noise $\varepsilon$. By
applying the saddle point equation, we obtain

$\displaystyle\Phi_{\beta}(n)=\mathop{\mathrm{extr}}_{\\{q^{(ab)}\\},\\{m^{(a)}\\},\\{\widehat{q}^{(ab)}\\},\\{\widehat{m}^{(a)}\\}}\left\\{\Phi_{\beta}^{(C)}(n)+\alpha\Phi^{(E)}_{\beta}(n)+\Phi^{(S)}_{\beta}(n)\right\\},$
(6.18)

where

$\displaystyle\Phi_{\beta}^{(C)}(n)$ $\displaystyle=-\sum_{a\leq
b}q^{(ab)}\widehat{q}^{(ab)}-\sum_{a=1}^{n}m^{(a)}\widehat{m}^{(a)}$ (6.19)
$\displaystyle\Phi_{\beta}^{(E)}(n)$ $\displaystyle=\ln\int d\varepsilon
P_{\varepsilon}(\varepsilon)\int
d\widetilde{\bm{\theta}}P_{\theta}(\widetilde{\bm{\theta}}|{\cal
Q})\prod_{a=1}^{n}f^{\beta}\left({\cal
G}(\theta^{(0)}+\varepsilon)|\theta^{(a)}\right)$ (6.20)
$\displaystyle\Phi_{\beta}^{(S)}(n)$ $\displaystyle=\ln\int
d\widetilde{\bm{x}}\phi^{(0)}(x^{(0)})\phi^{\beta}(\widetilde{\bm{x}}^{(a)})$
(6.21) $\displaystyle\hskip 56.9055pt\times\exp\left(\sum_{a\leq
b}\widehat{q}^{(ab)}x^{(a)}x^{(b)}+\sum_{a=1}^{n}\widehat{m}^{(a)}x^{(0)}x^{(a)}\right).$

###  6.1 Replica symmetric assumption

For further calculations, we introduce replica symmetric (RS) assumption as

$\displaystyle(q^{(ab)},\widehat{q}^{(ab)})$
$\displaystyle=\left\\{\begin{array}[]{ll}\left(Q,-\displaystyle\frac{\widehat{Q}}{2}\right)&{\rm
if}~{}a=b\\\ (q,\widehat{q})&{\rm if}~{}a\neq b\\\ \end{array},\right.$ (6.24)
$\displaystyle(m^{(a)},\widehat{m}^{(a)})$
$\displaystyle=(m,\widehat{m})~{}{\rm for~{}any}~{}a.$ (6.25)

RS assumption restricts the extremum of (6.18) to be in the space where (6.24)
and (6.25) hold. Therefore, it is always necessary to examine the validity of
the solution under the RS assumptions. The RS formula at zero-temperature
derived here has been proven to be exact for convex likelihoods in [65], and
that for Bayes-optimal inference was proven to be rigorous in [74]. The RS
assumption is known as to be equivalent with the assumption introduced in the
cavity approach [59].

For the integral of $\theta$ in (6.20) under the RS assumption, it is
convenient to introduce the transformation using independent Gaussian random
variables $w^{(a)}~{}(a\in\bm{n})$, $z$, and $v$ as

$\displaystyle\theta^{(a)}=\sqrt{Q-q}w^{(a)}+\sqrt{q}z,~{}(a\neq 0)$ (6.26)
$\displaystyle\theta^{(0)}=\sqrt{\rho\sigma_{x}^{2}-\frac{m^{2}}{q}}v+\frac{m}{\sqrt{q}}z,$
(6.27)

where this transformation preserves the relationships (6.8)–(6.10). We denote
the functions $\Phi_{\beta}^{(C)}(n)$, $\Phi_{\beta}^{(E)}(n)$, and
$\Phi_{\beta}^{(S)}(n)$ under the RS assumption as
$\Phi_{\beta}^{(C,\mathrm{RS})}(n)$, $\Phi_{\beta}^{(E,\mathrm{RS})}(n)$, and
$\Phi_{\beta}^{(S,\mathrm{RS})}(n)$, respectively. Using them, we obtain

$\displaystyle\Psi_{\beta}^{(C,\mathrm{RS})}$ $\displaystyle\equiv-\lim_{n\to
0}\frac{\partial}{\partial n}\Phi_{\beta}^{(C,\mathrm{RS})}(n)$ (6.28)
$\displaystyle=\widehat{m}m-\frac{\widehat{Q}Q+\widehat{q}q}{2}$
$\displaystyle\Psi_{\beta}^{(E,\mathrm{RS})}$ $\displaystyle\equiv-\lim_{n\to
0}\frac{\partial}{\partial n}\Phi_{\beta}^{(E,\mathrm{RS})}(n)$ (6.29)
$\displaystyle={-\int d\varepsilon P_{\varepsilon}(\varepsilon)\int Dz\int
Dv\ln\widehat{\xi}^{(\mathrm{RS})}_{\beta}(v,z,\varepsilon)}$
$\displaystyle\Psi_{\beta}^{(S,\mathrm{RS})}$ $\displaystyle\equiv-\lim_{n\to
0}\frac{\partial}{\partial n}\Phi_{\beta}^{(S,\mathrm{RS})}(n)$ (6.30)
$\displaystyle=-\int dx^{(0)}\phi^{(0)}(x^{(0)})\int
Dz\ln\xi^{(\mathrm{RS})}_{\beta}(z,x^{(0)}),$

and the free energy under the RS assumption is given by

$\displaystyle\Psi^{(\mathrm{RS})}=\lim_{\beta\to\infty}\frac{1}{\beta}\left\\{\Psi_{\beta}^{(C,\mathrm{RS})}+\alpha\Psi_{\beta}^{(E,\mathrm{RS})}+\Psi_{\beta}^{(S,\mathrm{RS})}\right\\},$
(6.31)

where we set

$\displaystyle\widehat{\xi}^{(\mathrm{RS})}_{\beta}(v,z)={\int Dw}$
$\displaystyle\hskip 28.45274pt{\times f^{\beta}\left({\cal
G}\left(\sqrt{\sigma_{0}^{2}-\frac{m^{2}}{q}}v+\frac{m}{\sqrt{q}}z+\varepsilon\right)\Big{|}\sqrt{Q-q}w+\sqrt{q}z\right)}$
(6.32) $\displaystyle\xi^{(\mathrm{RS})}_{\beta}(z,x^{(0)})=\int
dx\phi^{\beta}(x)\exp\left(-\frac{\widehat{Q}+\hat{q}}{2}x^{2}+(\sqrt{\hat{q}}z+\hat{m}x^{(0)})x\right),$
(6.33)

and $\sigma_{0}^{2}=\rho\sigma_{x}^{2}$.

####  6.1.1 $\beta\to\infty$

At $\beta\to\infty$, we implement some integrals by introducing the saddle
point method. We scale the variables as $\chi=\beta(Q-q)\sim O(1)$ and
$\widehat{Q}+\widehat{q}=\beta\widehat{\Theta}$,
$\widehat{q}=\beta^{2}\widehat{\chi}$, and $\widehat{m}=\beta\widehat{\mu}$.
Then, we obtain

$\displaystyle\lim_{\beta\to\infty}\frac{1}{\beta}\Psi_{\beta}^{(C,\mathrm{RS})}=\widehat{\mu}\mu-\\!\frac{\widehat{\Theta}Q-\widehat{\chi}\chi}{2}.$
(6.34)

Next, we introduce variable $u=\sqrt{Q-q}w+\sqrt{q}z$ and replace $w$ with it
as $w=(u-\sqrt{q}z)/\sqrt{Q-q}$. Further, by applying the transformation from
$(v,z)$ to $(\nu,\zeta)$, which keeps $dvdz=d\nu d\zeta$, as

$\displaystyle\sqrt{\sigma_{0}^{2}-\frac{m^{2}}{Q}}v+\frac{m}{\sqrt{Q}}z$
$\displaystyle=\sigma_{0}\nu$ (6.35)
$\displaystyle-\frac{m}{\sqrt{Q}}v+\sqrt{\sigma_{0}^{2}-\frac{m^{2}}{Q}}z$
$\displaystyle=\sigma_{0}\zeta,$ (6.36)

hence $\sqrt{Q}z=(m\nu+\sqrt{Q\sigma_{0}^{2}-m^{2}}\zeta)/\sigma_{0}$. With
these variables, we obtain

$\displaystyle\lim_{\beta\to\infty}\frac{1}{\beta}\ln\widehat{\xi}^{(\mathrm{RS})}_{\beta}$
$\displaystyle=\widehat{f}_{\xi}^{(\mathrm{RS})*}(\nu,\zeta,\varepsilon),$
(6.37)

where
$\widehat{f}_{\xi}^{(\mathrm{RS})*}(\nu,\zeta,\varepsilon)=\widehat{f}_{\xi}^{(\mathrm{RS})}(u_{\mathrm{RS}}^{*}(\nu,\zeta,\varepsilon);\nu,\zeta,\varepsilon)$
with

$\displaystyle\widehat{f}_{\xi}^{(\mathrm{RS})}(u;\nu,\zeta,\varepsilon)$
$\displaystyle=-\frac{(u-(m\nu+\sqrt{Q\sigma_{0}^{2}-m^{2}}\zeta)(\sigma_{0})^{-1})^{2}}{2\chi}+\ln
f\left({\cal G}\left(\sigma_{0}\nu+\varepsilon\right)|u\right)$ (6.38)
$\displaystyle
u_{\mathrm{RS}}^{*}(\nu,\zeta,\varepsilon)=\mathop{\mathrm{argmax}}_{u}\widehat{f}_{\xi}(u;\nu,\zeta,\varepsilon)$
(6.39)

In the case of GLM, $u^{*}(\nu,\zeta)$ should satisfy

$\displaystyle\frac{u_{\mathrm{RS}}^{*}-(m\nu+\sqrt{Q\sigma_{0}^{2}-m^{2}}\zeta)\sigma_{0}^{-1}}{\chi}={\cal
G}\left(\sigma_{0}\nu+\varepsilon\right)-a^{\prime}(u^{*}_{\mathrm{RS}}).$
(6.40)

For $\Psi_{\beta}^{(S,\mathrm{RS})}$, under the scaling of the parameters at
$\beta\to\infty$, we obtain

$\displaystyle\lim_{\beta\to\infty}\frac{1}{\beta}\ln\xi^{(\mathrm{RS})}_{\beta}(z,x^{(0)})$
$\displaystyle=f^{(\mathrm{RS})*}_{\xi}(z,x^{(0)})$ (6.41)

where
$f^{(\mathrm{RS})*}_{\xi}(z,x^{(0)})=f^{(\mathrm{RS})}_{\xi}(x_{\mathrm{RS}}^{*}(z,x^{(0)});z,x^{(0)})$
with

$\displaystyle f^{(\mathrm{RS})}_{\xi}(x;z,x^{(0)})$
$\displaystyle=-\frac{\widehat{\Theta}}{2}x^{2}+(\sqrt{\widehat{\chi}}z+\widehat{\mu}x^{(0)})x+\ln\phi(x)$
(6.42) $\displaystyle x_{\mathrm{RS}}^{*}(z,x^{(0)})$
$\displaystyle=\mathop{\mathrm{argmax}}_{x}f^{(\mathrm{RS})}_{\xi}(x;z,x^{(0)}).$
(6.43)

In summary, the free energy density under RS assumption is given by

$\displaystyle\Psi^{(\mathrm{RS})}$
$\displaystyle=\mathop{\mathrm{extr}}_{\Omega,\widehat{\Omega}}\Big{[}\widehat{\mu}m-\\!\frac{\widehat{\Theta}Q-\widehat{\chi}\chi}{2}-\int
dx^{(0)}\phi^{(0)}(x^{(0)})\int Dzf_{\xi}^{(\mathrm{RS})*}(z,x^{(0)})$ (6.44)
$\displaystyle\hskip 28.45274pt-\alpha{\int D\zeta D\nu d\varepsilon
P_{\varepsilon}(\varepsilon)\widehat{f}_{\xi}^{(\mathrm{RS})*}(\nu,\zeta,\varepsilon)}\Big{]}$

where $\Omega=\\{\mu,Q,\chi\\}$,
$\widehat{\Omega}=\\{\widehat{\mu},\widehat{\Theta},\widehat{\chi}\\}$, and
$\mathrm{extr}_{\Omega,\widehat{\Omega}}$ denotes the extremization with
respect to $\Omega$ and $\widehat{\Omega}$.

Hereafter, we consider the Gaussian noise with mean 0 and variance
$\sigma^{2}$. In this case, we can implement the integral with respect to
$\varepsilon$, and the resulting form of (6.38) is given as

$\displaystyle\widehat{f}_{\xi}^{(\mathrm{RS})}(u;\nu,\zeta)$
$\displaystyle=-\frac{(u-(m\nu+\sqrt{Q\sigma_{T}^{2}-m^{2}}\zeta)(\sigma_{T})^{-1})^{2}}{2\chi}+\ln
f\left({\cal G}\left(\sigma_{T}\nu\right)|u\right)$ (6.45)

where $\sigma_{T}^{2}=\sigma_{0}^{2}+\sigma^{2}$. The corresponding maximizer
is given as
$\widehat{f}_{\xi}^{(\mathrm{RS})*}(\nu,\zeta)=\widehat{f}_{\xi}^{(\mathrm{RS})}(u_{\mathrm{RS}}^{*}(\nu,\zeta),\nu,\zeta)$
with
$u_{\mathrm{RS}}^{*}=\mathrm{argmax}_{u}\widehat{f}_{\xi}^{(\mathrm{RS})}(u;\nu,\zeta)$.

####  6.1.2 Saddle point equations and simplification of the free energy
density

From the extremization conditions of the free energy density, we obtain the
saddle point equations for macroscopic quantities as

$\displaystyle Q$ $\displaystyle=\int dx^{(0)}\phi^{(0)}(x^{(0)})\int
Dz(x_{\mathrm{RS}}^{*}(z,x^{(0)}))^{2}$ (6.46) $\displaystyle\chi$
$\displaystyle=\frac{1}{\sqrt{\widehat{\chi}}}\int
dx^{(0)}\phi^{(0)}(x^{(0)})\int Dz\frac{\partial
x_{\mathrm{RS}}^{*}(z,x^{(0)})}{\partial z}$ (6.47) $\displaystyle m$
$\displaystyle=\int dx^{(0)}\phi^{(0)}(x^{(0)})\int
Dzx_{\mathrm{RS}}^{*}(z,x^{(0)})x^{(0)},$ (6.48)

and for conjugate variables

$\displaystyle\widehat{\Theta}$ $\displaystyle=\alpha\int D\zeta
D\nu~{}\frac{a^{\prime\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta))}{1+\chi
a^{\prime\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta))}$ (6.49)
$\displaystyle\widehat{\chi}$ $\displaystyle=\alpha\int D\zeta D\nu~{}({\cal
G}(\sigma_{T}\nu,\varepsilon)-a^{\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta)))^{2}$
(6.50) $\displaystyle\widehat{\mu}$ $\displaystyle=\alpha\int D\zeta
D\nu\frac{{\cal G}^{\prime}(\sigma_{T}\nu)}{1+\chi
a^{\prime\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta))}.$ (6.51)

The derivations of (6.49)-(6.51) are shown in E.1.

Using the saddle point equations, we can simplify
$\Psi_{\beta}^{(S,\mathrm{RS})}$ and $\Psi_{\beta}^{(E,\mathrm{RS})}$ as

$\displaystyle\frac{1}{\beta}\Psi_{\beta}^{(S,\mathrm{RS})}$
$\displaystyle=\frac{\widehat{\Theta}Q}{2}-\widehat{\chi}\chi-\widehat{\mu}m-\int
dx^{(0)}\phi^{(0)}(x^{(0)})\int Dz\ln\phi(x_{\mathrm{RS}}^{*}(z,x^{(0)}))$
(6.52) $\displaystyle\frac{1}{\beta}\Psi_{\beta}^{(E,\mathrm{RS})}$
$\displaystyle=\frac{\chi\widehat{\chi}}{2}-\alpha\int D\zeta D\nu\ln f({\cal
G}(\sigma_{T}\nu)|u_{\mathrm{RS}}^{*}(\nu,\zeta)),$ (6.53)

and the free energy density as

$\displaystyle\Psi^{\mathrm{RS}}$ $\displaystyle=-\alpha\int D\zeta D\nu\ln
f({\cal G}(\sigma_{T}\nu)|u_{\mathrm{RS}}^{*}(\nu,\zeta))$ (6.54)
$\displaystyle\hskip 56.9055pt-\int dx^{(0)}\phi^{(0)}(x^{(0)})\int
Dz\ln\phi(x_{\mathrm{RS}}^{*}(z,x^{(0)})).$

The first and the second terms of (6.54) correspond to the expectation of the
training error and regularization term, respectively, hence we can obtain the
training error from the free energy density (6.54), by discarding the
contribution from the penalty term.

####  6.1.3 Relationship with GAMP

The calculation by the replica analysis has correspondence with GAMP. The
equations for determining $u^{*}_{\mathrm{RS}}$ and $x^{*}_{\mathrm{RS}}$ in
replica method, given by (6.39) and (6.43), respectively, have the equivalent
form with (4.30) and (4.26), which appear in GAMP. In table 1, we summarize
the correspondence of the quantities appear in the replica method and GAMP.
The difference between GAMP and the replica method is that GAMP is defined for
one realization of data, and meanwhile replica method considers the typical
property with respect to the data distribution. For example, as shown in table
1, the estimate $\widehat{x}_{i}$ in GAMP is defined on given data ${\cal D}$,
and has site-dependency. In the replica method, the estimate
$x^{*}_{\mathrm{RS}}$ is defined on random variables $z$ and $x^{(0)}$. These
random variables effectively represent the randomness induced by
$\bm{\varepsilon}$, $\bm{F}$, and $\bm{x}^{(0)}$, and the estimates in GAMP
and the replica method are expected to be statistically equivalent.

GAMP | Replica (RS)
---|---
$\widehat{x}_{i}$ | $x_{\mathrm{RS}}^{*}(z,x^{(0)})$
$s$ | $\chi$
$\theta^{*}_{\mu}$ | $u_{\mathrm{RS}}^{*}(\nu,\zeta)$
${s_{\theta}}$ | $\chi$
$\widehat{\theta}_{\mu}^{\setminus\mu}$ | $(m\nu-\sqrt{Q\sigma_{T}^{2}-m^{2}}\zeta)\sigma_{T}^{-1}$
$\Sigma$ | $\widehat{\Theta}^{-1}$
$\mathrm{m}_{i}\Sigma^{-1}$ | $\sqrt{\widehat{\chi}}z+\widehat{\mu}x^{(0)}$
$\frac{1}{M}\sum_{\mu=1}^{M}(\check{g}_{\mathrm{out}})_{\mu}^{2}$ | $\widehat{\chi}$
$-\frac{1}{M}\sum_{\mu=1}^{M}(\check{\partial}_{\hat{\theta}}g_{\mathrm{out}})_{\mu}$ | $\widehat{\Theta}$
Table 1: Correspondence between GAMP and replica method under RS assumption.

###  6.2 Analytical form of extra-sample prediction error

The training error is naturally derived from the free energy density as shown
in (6.54), but we need additional computation for deriving the extra-sample
prediction error. Following the definition of the extra-sample prediction
error, we consider new predictor $\mathbf{f}_{\mathrm{new}}$ and define

$\displaystyle\theta^{0}_{\mathrm{new}}$
$\displaystyle=\frac{1}{\sqrt{N}}\mathbf{f}_{\mathrm{new}}\bm{x}^{(0)}$ (6.55)
$\displaystyle\theta_{\mathrm{new}}$
$\displaystyle=\frac{1}{\sqrt{N}}\mathbf{f}_{\mathrm{new}}\widehat{\bm{x}}({\cal
D}),$ (6.56)

and set the test data as $y_{\mathrm{new}}={\cal
G}(\theta_{\mathrm{new}}^{0}+\varepsilon)$. The expectation value of the
extra-sample prediction error is defined by

$\displaystyle\mathbb{E}_{{\cal
D},\mathbf{f}_{\mathrm{new}},\varepsilon}[\mathrm{err}_{\mathrm{pre}}({\cal
D})]=-\mathbb{E}_{{\cal D},\mathbf{f}_{\mathrm{new}},\varepsilon}[\ln
f(y_{\mathrm{new}}|\theta_{\mathrm{new}}(\mathbf{f}_{\mathrm{new}},\widehat{\bm{x}}({\cal
D})))].$ (6.57)

Here, we consider that the each component of the new predictor
$\mathbf{f}_{\mathrm{new}}$ is distributed according to the Gaussian
distribution with mean 0 and variance 1. For the expectation in the framework
of replica method under RS assumption, we introduce Gaussian random variables
$v$ and $w$ with mean zero and variance 1 and set

$\displaystyle\theta^{0}_{\mathrm{new}}$
$\displaystyle=\sqrt{\sigma_{x}^{2}-\frac{m^{2}}{Q}}v+\frac{m}{\sqrt{Q}}w$
(6.58) $\displaystyle\theta_{\mathrm{new}}$ $\displaystyle=\sqrt{Q}w,$ (6.59)

to satisfy the conditions

$\displaystyle\mathbb{E}_{\mathbf{f}_{\mathrm{new}}}[\theta^{0}_{\mathrm{new}}\theta_{\mathrm{new}}]$
$\displaystyle=m$ (6.60)
$\displaystyle\mathbb{E}_{\mathbf{f}_{\mathrm{new}}}[\theta^{0}_{\mathrm{new}}\theta_{\mathrm{new}}^{0}]$
$\displaystyle=\rho\sigma_{x}^{2}$ (6.61)
$\displaystyle\mathbb{E}_{\mathbf{f}_{\mathrm{new}}}[\theta_{\mathrm{new}}\theta_{\mathrm{new}}]$
$\displaystyle=Q.$ (6.62)

Using (6.55) and (6.56), the extra-sample prediction error for GLM is given by

$\displaystyle\mathbb{E}_{{\cal
D}}[\mathrm{err}_{\mathrm{pre}}^{(\mathrm{in})}({\cal D})]$
$\displaystyle=-\int DvDw\Big{\\{}{\cal
G}\left(\sqrt{\sigma_{T}^{2}-\frac{m^{2}}{Q}}v+\frac{m}{\sqrt{Q}}w\right)\sqrt{Q}w$
(6.63) $\displaystyle\hskip 28.45274pt-a(\sqrt{Q}w)+b\left({\cal
G}\left(\sqrt{\sigma_{T}^{2}-\frac{m^{2}}{Q}}v+\frac{m}{\sqrt{Q}}w\right)\right)\Big{\\}}.$

For the linear regression, we obtain

$\displaystyle\mathbb{E}_{{\cal
D}}\left[\mathrm{err}_{\mathrm{pre}}^{(\mathrm{ex})}({\cal D})\right]$
$\displaystyle=\frac{Q-2m+\sigma_{T}^{2}}{2},$ (6.64)

and for the logistic regression, we obtain

$\displaystyle\mathbb{E}_{{\cal
D}}\left[\mathrm{err}_{\mathrm{pre}}^{(\mathrm{ex})}({\cal D})\right]$
$\displaystyle=-\int Dw\left\\{\frac{\sqrt{Q}}{2}\frac{\partial}{\partial
w}\mathrm{erfc}\left(-\frac{\frac{m}{\sqrt{Q}}w}{\sqrt{2(\sigma_{T}^{2}-\frac{m^{2}}{Q})}}\right)-a(\sqrt{Q}w)\right\\}$
(6.65) $\displaystyle=-\frac{m}{\sqrt{2\pi\sigma_{T}^{2}}}+\int
Dwa(\sqrt{Q}w).$

The behaviors of the extra-sample prediction error and corresponding
generalization gap by replica method are shown in Figs. 5-10. The results by
replica method coincides with that by GAMP.

###  6.3 Analytical forms of estimators for the generalization gap

####  6.3.1 GDF by replica method under RS assumption

The data $y$ is replaced with $\sigma_{T}\nu$ in the replica analysis, hence
GDF is defined by the derivative of $u_{\mathrm{RS}}^{*}(\nu)$ with respect to
$\sigma_{T}\nu$. From the solution of $u_{\mathrm{RS}}^{*}(\nu)$ (6.39), we
obtain

$\displaystyle\frac{\partial
u_{\mathrm{RS}}^{*}(\nu,\zeta)}{\partial(\sigma_{T}\nu)}=\frac{\chi}{1+\chi
a^{\prime\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta))},$ (6.66)

and hence the expectation of GDF by replica method is derived as

$\displaystyle\mathbb{E}_{{\cal D}}[\mathrm{gdf}]=\int D\nu
D\zeta\frac{\chi}{1+\chi a^{\prime\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta))}.$
(6.67)

This form is consistent with GDF derived by GAMP (4.57): ${s_{\theta}}$ in
GAMP corresponds to $\chi$ in the replica analysis, and the summation with
respect to $\theta_{\mu}^{*}$ is replaced with the integral with respect to
$u_{\mathrm{RS}}^{*}(\nu,\zeta)$, as shown in Table 1. The behaviors of GDF
and corresponding in-sample prediction errors are shown in Figs.5-10. The GDF
by replica method coincident with the expectation value of GDF calculated by
GAMP.

####  6.3.2 FV by replica method under RS assumption

For deriving FV in the framework of the replica method, we introduce the
auxiliary variable $\gamma$, following the similar discussion with GAMP, as

$\displaystyle{\cal T}_{\gamma}^{(\mathrm{RS})}(\gamma)$
$\displaystyle=\gamma\ln f({\cal G}(\sigma_{T}\nu)|u^{*}(\gamma))$ (6.68)
$\displaystyle\hskip
28.45274pt-\frac{(u_{\gamma}^{*}(\gamma)-(m\nu+\\!\\!\sqrt{Q\sigma_{x}^{2}-m^{2}}\zeta)\sigma_{T}^{-1})^{2}}{2\chi}$
$\displaystyle u_{\gamma}^{*}(\gamma,\nu,\zeta)$
$\displaystyle=\mathop{\mathrm{argmax}}_{u}\Big{\\{}\gamma\ln f({\cal
G}(\sigma_{T}\nu)|u)$ (6.69) $\displaystyle\hskip
28.45274pt-\frac{(u-(m\nu+\\!\\!\sqrt{Q\sigma_{x}^{2}-m^{2}}\zeta)\sigma_{T}^{-1})^{2}}{2\chi}\Big{\\}}.$

Using these functions, the expectation of FV by the replica method is given by

$\displaystyle\mathbb{E}_{{\cal D}}[\mathrm{FV}]$ $\displaystyle=\int D\nu
D\zeta\frac{\partial^{2}}{\partial\gamma^{2}}{\cal
T}_{\gamma}^{(\mathrm{RS)}}(\gamma)\Big{|}_{\gamma=1}$ (6.70)
$\displaystyle=\int D\nu D\zeta\frac{({\cal
G}(\sigma_{T}\nu)-u_{\mathrm{RS}}^{*}(\nu,\zeta))^{2}\chi}{1+\chi
a^{\prime\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta))}.$

Eq.(6.70) can be interpreted as the expectation value of FV under GAMP from
the correspondence shown in Tabel 1.

###  6.4 Instability of RS solution

The RS solutions discussed so far lose local stability when subjected to
perturbations that break the symmetry between replicas in a certain parameter
region. This phenomenon is known as de Almeida-Thouless (AT) instability [71]
and occurs when the following conditions hold;

$\displaystyle\alpha\int Dzdx^{(0)}\phi^{(0)}(x^{(0)})\left(\frac{\partial
x_{\mathrm{RS}}^{*}(z,x^{(0)})}{\partial\sqrt{\widehat{\chi}}z}\right)^{2}$
$\displaystyle\times\int D\nu
D\zeta\left(\frac{\partial}{\partial\sqrt{Q-\frac{m^{2}}{\sigma_{T}^{2}}}\zeta}\left({\cal
G}\left(\sigma_{T}\nu\right)-a^{\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta))\right)\right)^{2}>1.$
(6.71)

The derivation of the AT instability condition is explained in E.3, by
introducing the solution under the 1 step replica symmetry breaking
assumption. Based on the correspondence between AMP and the RS saddle point
discussed so far and summarized in Table 1, the AT instability condition
(6.71) is coincident with the local instability condition of the GAMP fixed
point. Similar correspondence has been shown in the analysis for CDMA [57].

###  6.5 Specific forms of the saddle point equations

####  6.5.1 Gaussian likelihood

The saddle point equations for the conjugate variables can be simplified in
case of the Gaussian likelihood. The maximizer $u^{*}_{\mathrm{RS}}$ of (6.39)
is given by

$\displaystyle
u_{\mathrm{RS}}^{*}(\nu,\zeta)=\frac{\chi\sigma_{T}\nu+(m\nu+\\!\\!\sqrt{Q\sigma_{T}^{2}-m^{2}}\zeta)\sigma_{T}^{-1}}{1+\chi}$
(6.72)

hence we obtain

$\displaystyle\int D\zeta D\nu\ln\widehat{\xi}_{\beta}^{(\mathrm{RS})}$
$\displaystyle=\int D\zeta
D\nu\left(-\frac{(\sigma_{T}\nu-(m\nu+\\!\\!\sqrt{Q\sigma_{T}^{2}-m^{2}}\zeta)\sigma_{T}^{-1})^{2}}{2(1+\chi)}\right)$
(6.73) $\displaystyle=-\frac{Q-2m+\sigma_{T}^{2}}{2(1+\chi)}.$

Using (6.75), we obtain the saddle point equations

$\displaystyle\widehat{\Theta}$ $\displaystyle=\frac{\alpha}{1+\chi}$ (6.74)
$\displaystyle\widehat{\chi}$
$\displaystyle=\frac{\alpha(Q-2m+\sigma_{T}^{2})}{(1+\chi)^{2}}$ (6.75)
$\displaystyle\widehat{m}$ $\displaystyle=\frac{\alpha}{1+\chi}.$ (6.76)

Further, the AT instability condition is given by

$\displaystyle\frac{\alpha}{(1+\chi)^{2}}\int
Dzdx^{(0)}\phi^{(0)}(x^{(0)})\left(\frac{\partial
x_{\mathrm{RS}}^{*}(z,x^{(0)})}{\partial\sqrt{\widehat{\chi}}z}\right)^{2}>1.$
(6.77)

####  6.5.2 Elastic net penalty

The estimate $x_{\mathrm{RS}}^{*}$ (6.43) for the elastic net penalty is given
by

$\displaystyle
x_{\mathrm{RS}}^{*}=\frac{\sqrt{\widehat{\chi}}z+\widehat{\mu}x^{(0)}-\mathrm{sgn}(\sqrt{\widehat{\chi}}z+\widehat{\mu}x^{(0)})\lambda_{1}}{\widehat{\Theta}+\lambda_{2}}\mathbb{I}\left(|\sqrt{\widehat{\chi}}z+\widehat{\mu}x^{(0)}|>\lambda_{1}\right),$
(6.78)

hence the saddle point equations under the elastic net penalty are given by

$\displaystyle Q$
$\displaystyle=\frac{1}{(\widehat{\Theta}+\lambda_{2})^{2}}\int
d\sigma_{\rho}P_{\rho}(\sigma_{\rho})$ (6.79) $\displaystyle\hskip
28.45274pt\times\left\\{\sigma_{\rho}^{2}\left(1+2\Theta(\sigma_{\rho})^{2}\right)\mathrm{erfc}\left(\Theta(\sigma_{\rho})\right)-\frac{2\sigma_{\rho}}{\sqrt{\pi}}\Theta(\sigma_{\rho})\right\\}$
$\displaystyle m$
$\displaystyle=\frac{\rho\sigma_{x}\widehat{\mu}}{\widehat{\Theta}+\lambda_{2}}\mathrm{erfc}\left(\frac{\lambda_{1}}{\sqrt{2(\widehat{\chi}+\widehat{\mu}^{2}\sigma_{x}^{2})}}\right)$
(6.80) $\displaystyle\chi$
$\displaystyle=\frac{\widehat{\rho}}{\widehat{\Theta}+\lambda_{2}},$ (6.81)

where

$\displaystyle\Theta(\sigma_{\rho})$
$\displaystyle=\frac{\lambda_{1}}{\sqrt{2}\sigma_{\rho}},$ (6.82)
$\displaystyle\widehat{\rho}$ $\displaystyle=\int
d\sigma_{\rho}P_{\rho}(\sigma_{\rho})\mathrm{erfc}\left(\Theta(\sigma_{\rho})\right)$
(6.83) $\displaystyle P_{\rho}(\sigma_{\rho})$
$\displaystyle=(1-\rho)\delta(\sigma_{\rho}-\sqrt{\widehat{\chi}})+\rho\delta(\sigma_{\rho}-\sqrt{\widehat{\chi}+\widehat{\mu}^{2}\sigma_{x}^{2}}),$
(6.84)

and $\widehat{\rho}$ corresponds $\mathbb{E}_{{\cal
D}}[||\widehat{\bm{x}}({\cal D})||_{0}]/N$, namely expected value of the
fraction of non-zero component in the estimate. For the Gaussian likelihood,
the AT instability condition for the elastic net penalty is given by

$\displaystyle\frac{\alpha}{\widehat{\rho}}\left(\frac{\chi}{1+\chi}\right)^{2}>1.$
(6.85)

As discussed in sec.6.3.1, the term $\chi/(1+\chi)$ corresponds to the
expectation value of GDF.

####  6.5.3 $\ell_{2}$ penalty

For the $\ell_{2}$ penalty at $\rho=1$, the saddle point equation is given by

$\displaystyle Q$
$\displaystyle=\frac{\widehat{\chi}+\widehat{\mu}^{2}\sigma_{x}^{2}}{(\widehat{\Theta}+\lambda)^{2}}$
(6.86) $\displaystyle m$
$\displaystyle=\frac{\widehat{\mu}\sigma_{x}^{2}}{\widehat{\Theta}+\lambda}$
(6.87) $\displaystyle\chi$ $\displaystyle=\frac{1}{\widehat{\Theta}+\lambda},$
(6.88)

where we set $\lambda_{1}=0$, $\widehat{\rho}=1$ and $\lambda_{2}=\lambda$ of
the saddle point equations for the elastic net penalty (6.79)-(6.81). For the
Gaussian likelihood, the saddle point equation of $\chi$ can be solved by
using (6.88) and (6.75). In particular for $\lambda=0$ case, we can solve the
variable $\chi$ as

$\displaystyle\chi=\frac{1}{\alpha-1}.$ (6.89)

This expression indicates that $\chi\to\infty$ as $\alpha$ approaches 1 from
above. Using (6.89), we obtain

$\displaystyle\mathbb{E}_{{\cal
D}}[\mathrm{gdf}]=\frac{\chi}{1+\chi}=\frac{1}{\alpha},$ (6.90)

which corresponds to the GF or GDF using the hat matrix.

For the Gaussian likelihood with $\ell_{2}$ penalty of $\lambda=0$ case, we
term ridgeless limit, the following theorem with respect to the AT instability
and the series expansion, which is the basis for the derivation of WAIC as
(2.52), holds.

######  Theorem 5.

In the ridgeless limit, the instability of the RS solution and the divergence
of the series expansion of ${\cal T}(\eta)$ at $\eta=1$ arise at $\alpha<1$.

###### Proof.

From (6.85) and (6.90), the AT instability condition for the ridgeless limit
is given by

$\displaystyle\frac{1}{\alpha}>1.$ (6.91)

Therefore, AT instability always appears in $\alpha<1$ for $\lambda=0$.

Putting the expression (6.90) to the convergence condition of the series
expansion for the derivation of WAIC given by (4.71) with the correspondence
between $\chi$ and $s_{\theta}$, we obtain

$\displaystyle\frac{\eta}{\alpha}>1,$ (6.92)

hence the expansion at $\eta=1$ is inappropriate when ${\alpha}^{-1}>1$.
Therefore, in the ridgeless limit, the AT instability corresponds to the
inadequacy of the series expansion for deriving WAIC. ∎

####  6.5.4 $\ell_{1}$ penalty

In case of $\ell_{1}$ penalty, which corresponds to $\lambda_{2}=0$ of the
elastic net penalty, we obtain

$\displaystyle\chi=\frac{\widehat{\rho}}{\widehat{\Theta}}.$ (6.93)

For the Gaussian likelihood with $\ell_{1}$ penalty termed LASSO, we obtain
the expression using (6.74),

$\displaystyle\chi=\frac{\widehat{\rho}}{\alpha-\widehat{\rho}},$ (6.94)

and hence

$\displaystyle\mathbb{E}_{{\cal
D}}[\mathrm{gdf}]=\frac{\chi}{1+\chi}=\frac{\widehat{\rho}}{\alpha}.$ (6.95)

Eq. (6.95) indicates that the following theorem hold.

######  Theorem 6.

For LASSO, the number of the parameters divided by $M$, namely
$||\widehat{\bm{x}}({\cal D})||_{0}/M$, is an unbiased estimator of the
covariance penalty (2.37). Hence, the $C_{p}$ criterion is equivalent to AIC,
and they are unbiased estimators of the in-sample prediction error of LASSO
[34].

As with the ridgeless limit, we can show that following theorem for LASSO.

######  Theorem 7.

For LASSO, the instability of the RS solution and the divergence of the series
expansion of ${\cal T}(\eta)$ at $\eta=1$ arise at $\alpha<\widehat{\rho}$.

###### Proof.

From (6.85) and (6.95), AT instability arises when

$\displaystyle\frac{\widehat{\rho}}{\alpha}>1.$ (6.96)

Substituting (6.95) to the radius of convergence of the series expansion
(2.52), we obtain

$\displaystyle\frac{\eta\widehat{\rho}}{\alpha}>1.$ (6.97)

hence the invalidity of the series expansion at $\eta=1$ and the instability
of the RS solution arise at $\alpha<\widehat{\rho}$. ∎

##  7 Summary and Conclusion

We discussed prediction errors, generalization gaps and their estimators in
GLM. We derived the form of $C_{p}$ and WAIC using GAMP, and we discussed the
relationship between them. We showed that GDF and FV, which are the estimators
of the generalization gap, can be expressed by the variance of the model
parameter in the framework of GAMP.

We focused on the differences between the LOOCV error and WAIC or TIC at the
finite $\alpha$ region. The first-order expansion of the LOOCV error with
respect to the difference between the estimates under the full and LOO samples
is insufficient for describing the LOOCV error for the finite $\alpha$. Thus,
we derived the difference between the LOOCV error and WAIC within the
framework of GAMP by considering the higher order terms; we showed that the
difference can be described by the cumulant generating function of the
likelihood. The resultant form can be calculated under full data samples,
without considering LOO samples. Further, we extended GAMP-based estimators of
the prediction errors to deal with correlated predictors and showed that these
estimators exhibit a good agreement with LOOCV errors calculated by the linear
algebraic form or coordinate descent algorithm.

We show that GAMP or replica method describes the prediction accuracy of the
model by using the variance of the parameter at equilibrium, which
relationship can be regarded as an exemplification of fluctuation-response
relationship in the prediction based modeling [20]. The concept of the
prediction-based modeling is compatible with that in statistical physics. The
derivation of such a relationship in this study is restricted to the Gaussian
predictors, whose components independently and identically distributed
according to the Gaussian distribution. However, this setting is not as
restrictive as it appears. In fact, recent studies have shown Gaussian
universality for GLMs with a large class of features, wherein their asymptotic
behaviors can be computed to leading order under a simpler model with Gaussian
features [65, 75, 76, 77]. Utilizing this correspondence, the statistical
physics-based method is expected to be a valuable tool for analyzing the
statistical methods to find limitation in the conventional statistical theory,
analyze the $\alpha\sim O(1)$ region, and construct efficient algorithms.

The author thanks Yukito Iba for discussion on the information criteria as a
fluctuation-response relationship and teaching me regarding reference [40].
The author also thanks Yoshiyuki Kabashima for the helpful discussions and
comments. This work was partially supported by JST PRESTO Grant Number
JPMJPR19M2, Japan.

##  Appendix A Expansion of the leave-one-out cross validation error

The maximizer of $\varphi({\cal D},\bm{x})$ denoted by $\widehat{\bm{x}}({\cal
D})$ satisfies

$\displaystyle\frac{\partial}{\partial x_{j}}\varphi({\cal
D},\bm{x})\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal D})}=0$ (A.1)

for any $j\in\bm{N}$. Meanwhile, the maximizer of
$\sum_{\mu\in\bm{M}\setminus\nu}\ln
f(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x}))-h(\bm{x})$, denoted by
$\widehat{\bm{x}}({\cal D}_{\setminus\nu})$, satisfies

$\displaystyle\frac{\partial}{\partial
x_{j}}\left\\{\sum_{\mu\in\bm{M}\setminus\nu}\ln
f(y_{\mu}|\theta_{\mu}(\mathbf{f}_{\mu},\bm{x}))-h(\bm{x})\right\\}\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})}=0$ (A.2)

for any $j\in\bm{N}$. Hence, the following relationship holds:

$\displaystyle\frac{\partial}{\partial x_{j}}\varphi({\cal
D},\bm{x})\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})}=\frac{\partial}{\partial x_{j}}\ln
f(y_{\nu}|\theta_{\nu}(\mathbf{f}_{\nu},\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})}.$ (A.3)

Discarding (A.1) from (A.3), we obtain

$\displaystyle\frac{\partial\varphi({\cal D},\bm{x})}{\partial
x_{j}}\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})}\\!\\!\\!\\!-\frac{\partial\varphi({\cal
D},\bm{x})}{\partial x_{j}}\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D})}\\!\\!\\!\\!\\!=\frac{\partial}{\partial x_{j}}\ln
f(y_{\nu}|\theta_{\nu}(\mathbf{f}_{\nu},\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})}.$ (A.4)

Under the assumption A1 that the difference $\widehat{\bm{x}}({\cal
D})-\widehat{\bm{x}}({\cal D}_{\setminus\nu})$ is sufficiently small, we
expand the L.H.S. up to the first order as

$\displaystyle\frac{\partial}{\partial x_{j}}\varphi({\cal
D},\bm{x})\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})}-\frac{\partial}{\partial x_{j}}\varphi({\cal
D},\bm{x})\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal D})}$
$\displaystyle=\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})).$ (A.5)

Therefore, for any $j\in\bm{N}$ and $\nu\in\bm{M}$,

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

where we introduce the coefficient $1/M$ for convenience and define

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

With assumption A2, we assume that ${\cal J}_{kj}(\widehat{\bm{x}}({\cal
D}_{\setminus\nu}))={\cal J}_{kj}(\widehat{\bm{x}}({\cal D}))$ holds for any
$\nu\in\bm{M}$ and $i,j\in\bm{N}$, and we obtain

$\displaystyle\widehat{\bm{x}}({\cal D}_{\setminus\nu})-\widehat{\bm{x}}({\cal
D})=-({\cal J}^{\varphi}(\widehat{\bm{x}}({\cal
D}))^{-1}\frac{1}{M}\frac{\partial}{\partial\bm{x}}\ln
f(y_{\nu}|\bm{\theta}(\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})}.$ (A.8)

###  A.1 Further expansion

We continue the expansion of (A.4) for higher orders. For simplicity, we
consider that the case where the higher order derivative of the penalty term
$h(x)$ is greater than the second order is zero as assumption A4. Under this
assumption, we obtain

$\displaystyle\frac{\partial}{\partial x_{j}}\varphi({\cal
D},\bm{x})\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})}-\frac{\partial}{\partial x_{j}}\varphi({\cal
D},\bm{x})\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal D})}$
$\displaystyle=\sum_{k=1}^{N}(\widehat{x}_{k}({\cal
D}_{\setminus\nu})-\widehat{x}_{k}({\cal
D}))\Big{\\{}\frac{\partial^{2}\varphi({\cal D},\bm{x})}{\partial
x_{k}\partial x_{j}}\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})}$
$\displaystyle+\sum_{n=2}^{\infty}\frac{(-1)^{n}}{n!}\\!\\!\sum_{i_{1}=1}^{N}\\!\\!\cdots\\!\\!\\!\\!\sum_{i_{n-1}=1}^{N}\\!\\!\frac{\partial^{n+1}\varphi({\cal
D},\bm{x})}{\partial x_{i_{1}}\\!\cdots\partial x_{i_{n-1}}\partial
x_{k}\partial x_{j}}\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})}\\!\\!\prod_{\ell=1}^{n-1}(\widehat{x}_{i_{\ell}}({\cal
D}_{\setminus\nu})\\!-\\!\widehat{x}_{i_{\ell}}({\cal D}))\\!\Big{\\}}$
$\displaystyle=\sum_{k=1}^{N}(\widehat{x}_{k}({\cal
D}_{\setminus\nu})\\!-\\!\widehat{x}_{k}({\cal
D}))\left\\{\frac{\partial^{2}\varphi}{\partial x_{k}\partial
x_{j}}\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal
D}_{\setminus\nu})}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-\\!\frac{F_{\nu
k}F_{\nu
j}}{N}\\!\sum_{n=2}^{\infty}\\!\frac{(\widehat{\theta}_{\nu}\\!-\\!\widehat{\theta}_{\nu}^{\setminus\nu})^{n}}{n!}a^{\prime(n+1)}(\widehat{\theta}_{\nu}^{\setminus\nu})\\!\right\\}$
$\displaystyle=\sum_{k=1}^{N}(\widehat{x}_{k}({\cal
D}_{\setminus\nu})\\!-\\!\widehat{x}_{k}({\cal
D}))\left\\{\\!\frac{\partial^{2}\varphi}{\partial x_{k}\partial
x_{j}}\Big{|}_{\bm{x}=\widehat{\bm{x}}({\cal D}_{\setminus\nu})}-\frac{F_{\nu
k}F_{\nu
j}}{N}\\!\left(\\!\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}^{\setminus\nu})\\!\right)\\!\right\\},$

where the arguments of $\varphi$ are abbreviated, and
$\widehat{\theta}_{\mu}^{\setminus\nu}\sim\widehat{\theta}_{\mu}$ holds for
$\mu\neq\nu$. Note that in (A.1), we consider all possible combinations of
$n$-th order terms for any integer, that is, we are implicitly assuming that
$N\to\infty$. By defining the matrix

$\displaystyle{\cal K}^{\setminus\nu}={\cal
J}^{\varphi}(\widehat{\bm{x}}({\cal
D}_{\setminus\nu}))-\frac{\mathbf{f}_{\nu}^{\top}\mathbf{f}_{\nu}}{NM}\left(\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})\right)$
(A.9)

for $\nu\in\bm{M}$, we obtain

$\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})}.$ (A.10)

For the Gaussian likelihood, ${\cal K}^{\setminus\nu}={\cal J}$ holds for any
$\nu\in\bm{M}$.

##  Appendix B Generalized linear models

###  B.1 Logistic regression with Bernoulli distribution

The logistic regression on the Bernoulli distribution is defined for binary
data $y_{\mu}\in\\{0,1\\}$ as

$\displaystyle\ln
f(\bm{y}|\bm{\pi})=\sum_{\mu=1}^{M}\left\\{y_{\mu}\log\pi_{\mu}+(1-y_{\mu})\log(1-\pi_{\mu})\right\\}$
(B.1)

where the Bernoulli parameter $\pi_{\mu}$ provides the mean of $y_{\mu}$, and
the probability that $y_{\mu}$ is 1. In the logistic regression, the model
parameter $\bm{\theta}=\frac{1}{\sqrt{N}}\bm{F}\bm{x}$ is connected to the
Bernoulli parameter as

$\displaystyle\pi_{\mu}=\ell^{-1}({\theta}_{\mu})=\frac{1}{1+\exp(-\theta_{\mu})},$
(B.2)

using the link function $\ell(\cdot)$. Therefore,
$\ell(\pi_{\mu})=\ln\left(\pi_{\mu}/(1-\pi_{\mu})\right)$. Using (B.2), the
distribution (B.1) is transformed into the form of the exponential family
distribution

$\displaystyle\ln
f(\bm{y}|\bm{\theta})=\bm{\theta}^{\top}\bm{y}-\sum_{\mu=1}^{M}\log(1+e^{\theta_{\mu}}),$
(B.3)

where functions $a$ and $b$ respectively correspond to

$\displaystyle
a(\bm{\theta})=\log\prod_{\mu=1}^{M}(1+e^{\theta_{\mu}}),~{}~{}~{}b(\bm{y})=0.$
(B.4)

The maximization problem in GAMP (4.30) for the logistic regression is given
by

$\displaystyle\theta_{\mu}^{*}=\mathop{\mathrm{argmax}}_{\theta_{\mu}}\left\\{-\theta_{\mu}(1-y_{\mu})-\log(1+e^{-\theta_{\mu}})-\frac{(\theta_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu})^{2}}{2{s_{\theta}}_{\mu}}\right\\},$
(B.5)

and the estimate $\theta_{\mu}^{*}$ satisfies

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

The estimated mean of $y_{\mu}$ is given by
$\pi_{\mu}^{*}\equiv(1+e^{-\theta_{\mu}^{*}})^{-1}$, and the variance is
$\pi_{\mu}^{*}(1-\pi_{\mu}^{*})$; therefore, we get

$\displaystyle(\check{g}_{\rm out})_{\mu}$
$\displaystyle=y_{\mu}-\pi_{\mu}^{*},~{}~{}~{}(\partial_{\widehat{\theta}}\check{g}_{\rm
out})_{\mu}$
$\displaystyle=\frac{\pi_{\mu}^{*}(1-\pi_{\mu}^{*})}{1+{s_{\theta}}_{\mu}\pi_{\mu}^{*}(1-\pi_{\mu}^{*})}.$
(B.7)

###  B.2 Poisson regression

Poisson distribution is defined for $y_{\mu}\in\mathbb{N}$, where $\mathbb{N}$
denotes the set of natural numbers including zero as

$\displaystyle\ln
f(\bm{y}|\bm{\lambda})=\sum_{\mu=1}^{M}\left\\{y_{\mu}\ln\lambda_{\mu}-\lambda_{\mu}-\log
y_{\mu}!\right\\},$ (B.8)

where $\lambda_{\mu}$ denotes the mean of $y_{\mu}$. In Poisson regression,
the model parameter $\bm{\theta}$ is connected to the mean as

$\displaystyle\lambda_{\mu}=\ell^{-1}(\theta_{\mu})=\exp(\theta_{\mu}),$ (B.9)

and hence, $\ell(\lambda_{\mu})=\log(\lambda_{\mu})$. This link function
reduces Poisson distribution to the exponential family’s form as

$\displaystyle
a(\bm{\theta})=\sum_{\mu=1}^{M}\exp(\theta_{\mu}),~{}~{}~{}b(\bm{y})=-\sum_{\mu=1}^{M}\log
y_{\mu}!.$ (B.10)

The maximization problem in GAMP (4.30) for the Poisson regression is given by

$\displaystyle\theta_{\mu}^{*}=\max_{\theta}\left\\{y_{\mu}\theta_{\mu}-\exp(\theta_{\mu})-\frac{(u_{\mu}-\widehat{\theta}^{\setminus\mu}_{\mu})^{2}}{2{s_{\theta}}_{\mu}}\right\\},$
(B.11)

and the solution $\theta^{*}_{\mu}$ satisfies the condition

$\displaystyle\frac{\theta_{\mu}^{*}-\widehat{\theta}^{\setminus\mu}_{\mu}}{{s_{\theta}}_{\mu}}=y_{\mu}-\lambda_{\mu}^{*},$
(B.12)

where $\lambda_{\mu}^{*}\equiv\exp(u_{\mu}^{*})$ is the estimated average.
Therefore, we obtain

$\displaystyle(\check{g}_{\mathrm{out}})_{\mu}$
$\displaystyle=y_{\mu}-\lambda_{\mu}^{*},~{}~{}~{}(\check{\partial}_{\widehat{\theta}}{g}_{\mathrm{out}})_{\mu}$
$\displaystyle=\frac{\lambda_{\mu}^{*}}{1+{s_{\theta}}_{\mu}\lambda_{\mu}^{*}}.$
(B.13)

###  B.3 Exponential regression

Exponential distribution with mean $\lambda_{\mu}$ is defined for
$y_{\mu}\in[0,\infty)$ as

$\displaystyle\ln
f(\bm{y}|\bm{\eta})=\sum_{\mu=1}^{M}\left\\{-\frac{y_{\mu}}{\eta_{\mu}}-\log(\eta_{\mu})\right\\}.$
(B.14)

In the exponential regression, the model parameter $\bm{\theta}$ relates to
the mean as

$\displaystyle-\eta_{\mu}=\ell^{-1}(\theta_{\mu}),$ (B.15)

and hence, $\ell(\theta_{\mu})=-\eta_{\mu}^{-1}$, where we introduce the
factor $-1$ for convenience. Using the parameter $\bm{\theta}$, the
exponential distribution is rewritten as an exponential family with functions
$a$ and $b$ given by

$\displaystyle a(\bm{\theta})$
$\displaystyle=-\sum_{\mu=1}^{M}\log(-\theta_{\mu}),~{}~{}~{}b(\bm{y})$
$\displaystyle=0.$ (B.16)

The maximization problem in GAMP (4.30) for exponential regression is given by

$\displaystyle\theta_{\mu}^{*}=\max_{\theta}\left\\{y_{\mu}\theta_{\mu}+\ln(-\theta_{\mu})-\frac{(\theta_{\mu}-\widehat{\theta}^{\setminus\mu}_{\mu})^{2}}{2{s_{\theta}}_{\mu}}\right\\},$
(B.17)

and the solution satisfies

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

The estimated mean and variance corresponds to
$\eta_{\mu}^{*}\equiv-(\theta_{\mu}^{*})^{-1}$ and $(\theta_{\mu}^{*})^{-2}$,
respectively; hence we obtain

$\displaystyle(\check{g}_{\rm
out})_{\mu}=y_{\mu}-\eta_{\mu}^{*},~{}~{}~{}(\partial_{\widehat{\theta}}\check{g}_{\rm
out})_{\mu}=\frac{(\theta_{\mu}^{*})^{-2}}{1+{s_{\theta}}_{\mu}(\theta_{\mu}^{*})^{-2}}.$
(B.19)

###  B.4 Cumulant generating function for exponential family distribution

The moment-generating function for the exponential family distribution is
given by

$\displaystyle{\cal M}_{{\cal L}(\theta)}(t)$ $\displaystyle=\int
dy\exp(ty)f(y|\theta)$ (B.20)
$\displaystyle=\frac{\exp(a(\theta+t))}{\exp(a(\theta))}\int
dy\exp((t+\theta)y-a(t+\theta)+b(y)).$

Hence, the cumulant generating function ${\cal C}_{{\cal
L}(\theta)}(t)=\ln{\cal M}_{{\cal L}(\theta)}(t)$ is given by

$\displaystyle{\cal C}_{{\cal L}(\theta)}(t)=a(\theta+t)-a(\theta).$ (B.21)

The $k$-th cumulant for the distribution $f(y|\theta)$ denoted by
$\kappa_{k}[f(\cdot|\theta)]$ is given by

$\displaystyle\kappa_{k}[f(\cdot|\theta)]=\frac{\partial^{k}}{\partial
t^{k}}{\cal
C}_{f(\cdot|\theta)}(t)\Big{|}_{t=0}=\frac{\partial^{k}}{\partial\theta^{k}}a(\theta),$
(B.22)

hence the derivative of $a$ corresponds to the cumulant.

##  Appendix C Derivation of the updating rule of AMP

###  C.1 General form of belief propagation

The likelihood depends on the variable $\bm{x}$ through the parameter
$\bm{\theta}=\frac{1}{\sqrt{N}}\bm{F}\bm{x}$. Considering the input message to
$i$, we focus on
$\theta^{-i}_{\mu}\equiv{\theta}_{\mu}-\frac{1}{\sqrt{N}}F_{\mu i}x_{i}$ and
obtain the probability distribution of $\theta_{\mu}^{-i}$ using the messages.
For a variable transformation, we introduce the identity for all
$\mu\in\bm{M}$ as

$1=\int
d\theta_{\mu}^{-i}\delta\left(\theta_{\mu}^{-i}-\frac{1}{\sqrt{N}}\sum_{j\neq
i}^{N}F_{\mu j}x_{j}\right).$ (C.1)

The distribution of $\theta_{\mu}^{-i}$ at step $t$, defined as
$\mathrm{Pr}[\theta_{\mu}^{-i}=\theta]=W^{-i^{(t)}}_{\mu}(\theta)$, is given
by

$\displaystyle W_{\mu}^{-i^{(t)}}(\theta)$ $\displaystyle=\int
d\bm{x}_{\setminus i}\delta\left(\theta-\frac{1}{\sqrt{N}}\sum_{j\neq i}F_{\mu
j}x_{j}\right)\prod_{j\in\bm{N}\backslash i}\xi^{(t)}_{j\to\mu}(x_{j})$ (C.2)
$\displaystyle=\frac{1}{2\pi}\int
d\tilde{\theta}e^{-\sqrt{-1}\tilde{\theta}\theta}\prod_{j\in\bm{N}\setminus
i}\int dx_{j}\exp\left(\sqrt{-1}\tilde{\theta}\frac{F_{\mu
j}x_{j}}{\sqrt{N}}\right)\xi^{(t)}_{j\to\mu}(x_{j})$
$\displaystyle\simeq\frac{1}{2\pi}\int
d\tilde{\theta}e^{-\sqrt{-1}\tilde{\theta}\theta}\exp\left(\sqrt{-1}\tilde{\theta}\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}}-\frac{\widetilde{\theta}^{2}}{2}\beta^{-1}s_{\theta_{\mu}}^{\backslash\mu,-i^{(t)}}\right)$
$\displaystyle=\frac{1}{\sqrt{2\pi
s^{\backslash\mu,-i(t)}_{\theta_{\mu}}}}\exp\left(-\frac{\beta(\theta-\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}})^{2}}{2s_{\theta_{\mu}}^{\backslash\mu,-i^{(t)}}}\right),$

where the integral representation of the delta function with the variable
$\tilde{\theta}$ is introduced, and we define

$\displaystyle\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}}$
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{j\neq i}{F}_{\mu
j}\widehat{x}^{(t)}_{j\to\mu}$ (C.3) $\displaystyle
s^{\setminus\mu,-i^{(t)}}_{\theta_{\mu}}$
$\displaystyle=\frac{1}{N}\sum_{j\neq i}{F}^{2}_{\mu j}{s}^{(t)}_{j\to\mu}.$
(C.4)

Using the distribution $W_{\mu}^{-i^{(t)}}(\theta)$, the integral in the input
message is operated as

$\displaystyle\widetilde{\xi}^{(t)}_{\mu\to i}(x_{i})$
$\displaystyle\propto\int
d\theta_{\mu}d\theta_{\mu}^{-i}f^{\beta}(y_{\mu}|\theta_{\mu})\delta\left(\theta_{\mu}-\theta_{\mu}^{-i}-\frac{F_{\mu
i}x_{i}}{\sqrt{N}}\right)W_{\mu}^{-i^{(t)}}(\theta_{\mu}^{-i})$ (C.5)
$\displaystyle=\int
d\theta_{\mu}f^{\beta}(y_{\mu}|\theta_{\mu})\exp\left(-\frac{\beta(\theta_{\mu}-\frac{F_{\mu
i}x_{i}}{\sqrt{N}}-\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}})^{2}}{2s^{\setminus\mu,-i^{(t)}}_{\theta_{\mu}}}\right)$
$\displaystyle=\int
d\theta_{\mu}f^{\beta}(y_{\mu}|\theta_{\mu})\exp\left(-\frac{\beta(\theta_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}})^{2}}{2s^{\setminus\mu,-i^{(t)}}_{\theta_{\mu}}}\right)$
$\displaystyle\times\Big{\\{}1+\frac{\beta F_{\mu
i}x_{i}(\theta_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}})}{\sqrt{N}s^{\setminus\mu,-i^{(t)}}_{\theta_{\mu}}}+\frac{\beta^{2}F^{2}_{\mu
i}x^{2}_{i}(\theta_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}})^{2}}{2N(s^{\setminus\mu,-i^{(t)}}_{\theta_{\mu}})^{2}}$
$\displaystyle\hskip 28.45274pt-\frac{\beta F^{2}_{\mu
i}x^{2}_{i}}{2Ns^{\setminus\mu,-i^{(t)}}_{\theta_{\mu}}}+O(N^{-3/2})\Big{\\}}.$

Dividing by a constant does not change the proportionality, and therefore, we
divide both sides of (C.5) by
$\Xi_{\beta}(y_{\mu};s^{\setminus\mu,-i^{(t)}}_{\theta_{\mu}},\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}})$
defined as

$\displaystyle\Xi_{\beta}(y;{s_{\theta}},\widehat{\theta})$
$\displaystyle=\sqrt{\displaystyle\frac{\beta}{2\pi s_{\theta}}}\int d\theta
f_{\beta}(y|\theta)\exp\left(-\frac{\beta(\theta-\widehat{\theta})^{2}}{2s_{\theta}}\right).$
(C.6)

Further, we define the effective posterior distribution after the observation
of $y$ as

$\displaystyle
p_{\mathrm{c},\mu}(\theta_{\mu};{s_{\theta}}_{\mu},\widehat{\theta}_{\mu}^{\setminus\mu},y_{\mu})$
$\displaystyle\hskip
14.22636pt=\frac{1}{\Xi_{\beta}(y_{\mu};{s_{\theta}}_{\mu},\widehat{\theta}_{\mu}^{\setminus\mu})}\sqrt{\displaystyle\frac{\beta}{2\pi
s_{\theta}}}f_{\beta}(y_{\mu}|\theta_{\mu})\exp\left(-\frac{\beta(\theta_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu})^{2}}{2{s_{\theta}}_{\mu}}\right),$
(C.7)

and

$\displaystyle
g_{\mathrm{out}}({s_{\theta}}_{\mu},\widehat{\theta}_{\mu}^{\setminus\mu},y_{\mu})$
$\displaystyle=\frac{\partial}{\partial\hat{\theta}_{\mu}^{\setminus\mu}}\ln\Xi_{\beta}(y_{\mu};{s_{\theta}}_{\mu},\widehat{\theta}_{\mu}^{\setminus\mu})=\Big{\langle}\frac{\beta(\theta_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu})}{{s_{\theta}}_{\mu}}\Big{\rangle}_{\theta_{\mu}}$
(C.8)
$\displaystyle\partial_{\hat{\theta}}g_{\mathrm{out}}({s_{\theta}}_{\mu},\widehat{\theta}_{\mu}^{\setminus\mu},y_{\mu})$
$\displaystyle=\frac{\partial}{\partial\hat{\theta}_{\mu}^{\setminus\mu}}g_{\mathrm{out}}({s_{\theta}}_{\mu},\widehat{\theta}_{\mu}^{\setminus\mu},y_{\mu}),$
(C.9)

where $\langle\cdot\rangle_{\theta}$ denotes the expectation with respect to
$\theta$ according to the distribution
$p_{\mathrm{c},\mu}(\theta;{s_{\theta}},\widehat{\theta},y)$. The following
relationship exists:

$\displaystyle\partial_{\hat{\theta}}g_{\mathrm{out}}({s_{\theta}}_{\mu},\widehat{\theta}_{\mu}^{\setminus\mu},y_{\mu})$
$\displaystyle\\!=\\!-\frac{\beta}{{s_{\theta}}_{\mu}}+\\!\left\langle\\!\\!\left(\displaystyle\frac{\beta(\theta_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu})}{{s_{\theta}}_{\mu}}\right)^{2}\right\rangle_{\theta_{\mu}}\\!\\!\\!\\!\\!-g_{\mathrm{out}}^{2}({s_{\theta}}_{\mu},\widehat{\theta}_{\mu}^{\setminus\mu},y_{\mu}).$
(C.10)

Using Eqs.(C.8)-(C.10) and ignoring $O(N^{-3/2})$ terms, we obtain

$\displaystyle\widetilde{\xi}_{\mu\to i}^{(t)}(x_{i})$ $\displaystyle\propto
1+\frac{F_{\mu
i}x_{i}}{\sqrt{N}}g_{\mathrm{out}}(s_{\theta_{\mu}}^{\setminus\mu,-i^{(t)}},\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}},y_{\mu})$
(C.11) $\displaystyle+\frac{X^{2}_{i\ell}F^{2}_{\mu
i}}{2N}\Big{\langle}\frac{\beta^{2}(\theta_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}})^{2}}{(s^{\setminus\mu,-i^{(t)}}_{\theta_{\mu}})^{2}}\Big{\rangle}_{\theta_{\mu}}-\frac{\beta
x^{2}_{i}}{2Ns^{\setminus\mu,-i^{(t)}}_{\theta_{\mu}}}F^{2}_{\mu i}$
$\displaystyle=1+\frac{F_{\mu
i}x_{i}}{\sqrt{N}}g_{\mathrm{out}}(s_{\theta_{\mu}}^{\setminus\mu,-i^{(t)}},\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}},y_{\mu})$
$\displaystyle+\frac{F^{2}_{\mu
i}x^{2}_{i}}{2N}g_{\mathrm{out}}^{2}(s_{\theta_{\mu}}^{\setminus\mu,-i^{(t)}},\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}},y_{\mu})+\frac{x^{2}_{i}F^{2}_{\mu
i}}{2N}\Big{\langle}\frac{\beta^{2}(\theta_{\mu}-\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}})^{2}}{(s^{\setminus\mu,-i^{(t)}}_{\theta_{\mu}})^{2}}\Big{\rangle}_{\theta_{\mu}}$
$\displaystyle-\frac{\beta
x^{2}_{i}}{2Ns^{\setminus\mu,-i^{(t)}}_{\theta_{\mu}}}F^{2}_{\mu
i}-\frac{F^{2}_{\mu
i}x^{2}_{i}}{2N}g_{\mathrm{out}}^{2}(s_{\theta_{\mu}}^{\setminus\mu,-i^{(t)}},\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}},y_{\mu}),$

hence

$\displaystyle\widetilde{\xi}^{(t)}_{\mu\to i}(x_{i})$
$\displaystyle\propto\exp\left\\{-\frac{\beta x_{i}^{2}}{2N}A^{(t)}_{\mu\to
i}+x_{i}\frac{\beta B^{(t)}_{\mu\to i}}{\sqrt{N}}\right\\},$ (C.12)

where

$\displaystyle A_{\mu\to i}^{(t)}$ $\displaystyle=-F^{2}_{\mu
i}\beta^{-1}\partial_{\widehat{\theta}}g_{\mathrm{out}}(s_{\theta_{\mu}}^{\setminus\mu,-i^{(t)}},\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}},y_{\mu})$
(C.13) $\displaystyle B_{\mu\to i}^{(t)}$ $\displaystyle=F_{\mu
i}\beta^{-1}g_{\mathrm{out}}(s_{\theta_{\mu}}^{\setminus\mu,-i^{(t)}},\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}},y_{\mu}).$
(C.14)

Eq. (C.12) is a Gaussian distribution whose mean $\widetilde{x}^{(t)}_{\mu\to
i}$ and the rescaled variance $\widetilde{s}^{(t)}_{\mu\to i}$ are
respectively given by

$\displaystyle\widetilde{x}^{(t)}_{\mu\to i}$
$\displaystyle=\frac{\sqrt{N}B^{(t)}_{\mu\to i}}{A^{(t)}_{\mu\to
i}},~{}~{}~{}\widetilde{s}^{(t)}_{\mu\to i}$
$\displaystyle=\frac{N}{A^{(t)}_{\mu\to i}}.$ (C.15)

We immediately obtain the following expressions of the output messages using
Gaussian approximation (C.12):

$\displaystyle{\xi}^{(t+1)}_{i\to\mu}(x_{i})$
$\displaystyle\propto\phi^{\beta}(x_{i})\exp\left\\{-\sum_{\nu\in\bm{M}\setminus\mu}\frac{\beta}{2\widetilde{s}^{(t)}_{\nu\to
i}}\left(x_{i}^{2}-2x_{i}\widetilde{x}_{\nu\to i}^{(t)}\right)\right\\}.$
(C.16)

We define the following form of the distribution to express the output
messages and marginal distribution in the unified form:

$\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\\},$
(C.17)

where ${\cal Z}$ represents the normalization constant, and denotes its mean
and variance as $\mathbb{M}_{\beta}(\Sigma,\mathrm{m})$ and
$\beta^{-1}\mathbb{V}_{\beta}(\Sigma,\mathrm{m})$, respectively. Transforming
(C.16), the mean and rescaled variance of the output message, denoted by
$\widehat{x}_{i\to\nu}$ and ${s}_{i\to\nu}$, respectively, are given by

$\displaystyle\widehat{x}_{i\to\nu}^{(t+1)}=\mathbb{M}_{\beta}\left({\Sigma}^{(t)}_{i\to\nu},\mathrm{m}^{(t)}_{i\to\nu}\right),~{}~{}~{}{s}^{(t+1)}_{i\to\nu}=\beta\mathbb{V}_{\beta}\left({\Sigma}^{(t)}_{i\to\nu},\mathrm{m}^{(t)}_{i\to\nu}\right),$
(C.18)

where

$\displaystyle{\Sigma}^{(t)}_{i\to\mu}=\frac{1}{\sum_{\nu\in\bm{M}\backslash\mu}(\widetilde{s}^{(t)}_{\nu\to
i})^{-1}},~{}~{}~{}\mathrm{m}^{(t)}_{i\to\mu}=\frac{\sum_{\nu\in\bm{M}\backslash\mu}(\widetilde{s}^{(t)}_{\nu\to
i})^{-1}\widetilde{x}^{(t)}_{\nu\to
i}}{\sum_{\nu\in\bm{M}\backslash\mu}(\widetilde{s}^{(t)}_{\nu\to i})^{-1}}.$
(C.19)

Similarly, marginal distributions at time $t$ is given by

$\displaystyle p_{\beta,i}^{(t)}(x_{i})$ $\displaystyle={\cal
P}_{\beta}(x_{i};{\Sigma}_{i}^{(t)},\mathrm{m}_{i}^{(t)}),$ (C.20)

where

$\displaystyle{\Sigma}^{(t)}_{i}=\frac{1}{\sum_{\nu=1}^{M}(\widetilde{s}^{(t)}_{\nu\to
i})^{-1}},~{}~{}~{}\mathrm{m}^{(t)}_{i}=\frac{\sum_{\nu=1}^{M}(\widetilde{s}^{(t)}_{\nu\to
i})^{-1}\widetilde{x}^{(t)}_{\nu\to
i}}{\sum_{\nu=1}^{M}(\widetilde{s}^{(t)}_{\nu\to i})^{-1}},$ (C.21)

and the estimate and rescaled variance are given by

$\displaystyle\widehat{x}_{i}^{(t)}=\mathbb{M}_{\beta}({\Sigma}_{i}^{(t)},\mathrm{m}_{i}^{(t)}),~{}~{}~{}s_{i}^{(t)}=\beta\mathbb{V}_{\beta}({\Sigma}_{i}^{(t)},\mathrm{m}_{i}^{(t)}).$
(C.22)

###  C.2 TAP form of GAMP

An iteration step of BP comprises the updating of $2MN$ messages following
(C.13)-(C.15), (C.18), and (C.19) for $i=1,\cdots,N$ and $\mu=1,\cdots,M$.
Computational cost per iteration is $O(M^{2}\times N+M\times N^{2})$, and
therefore, the application of BP to large systems are not effective. However,
BP equations can be rewritten in terms of $M+N$ messages instead of $2MN$
messages under the assumption that matrix $\bm{F}$ is not sparse and all
elements scale as $O(1)$. The BP with reduced messages are known as
Thouless–Anderson–Palmer (TAP) equations [60] in the study of the statistical
physics for spin glasses. The TAP form is equivalent to the BP equations at a
sufficiently large $N$, and this is called the GAMP algorithm [58]. The TAP
form results in a significant reduction of computational complexity to $O(MN)$
per iteration.

First, we ignore $\frac{1}{N}F_{\mu i}^{2}s_{i\to\mu}\sim O(N^{-1})$ for a
sufficiently large system; hence
${s}_{\theta_{\mu}}^{\setminus\mu,-i^{(t)}}={s}^{\setminus\mu^{(t)}}_{\theta_{\mu}}$.
Further, using the relationship

$\displaystyle\widehat{\theta}_{\mu}^{\setminus\mu^{(t)}}$
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{j=1}^{N}F_{\mu
j}\widehat{x}_{j\to\mu}^{(t)}=\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}}+\frac{1}{\sqrt{N}}F_{\mu
i}\widehat{x}_{i\to\mu}^{(t)},$ (C.23)

and applying Taylor expansion with respect to $O(N^{-1/2})$ terms, we obtain

$\displaystyle
g_{\mathrm{out}}(s_{\theta_{\mu}}^{\setminus\mu,-i^{(t)}},\widehat{\theta}_{\mu}^{\setminus\mu,-i^{(t)}},y_{\mu})$
$\displaystyle=g_{\mathrm{out}}(s_{\theta_{\mu}}^{\setminus\mu^{(t)}},\widehat{\theta}_{\mu}^{\setminus\mu^{(t)}},y_{\mu})$
$\displaystyle-\partial_{\hat{\theta}}g_{\mathrm{out}}({s_{\theta}}^{\setminus\mu^{(t)}}_{\mu},\widehat{\theta}_{\mu}^{\setminus\mu^{(t)}},y_{\mu})\frac{F_{\mu
i}\widehat{x}_{i\to\mu}^{(t)}}{\sqrt{N}}+O(N^{-1}).$

For the sake of simplicity, we denote
$(g_{\mathrm{out}}^{(t)})_{\mu}=g_{\mathrm{out}}(s_{\theta_{\mu}}^{\setminus\mu^{(t)}},\widehat{\theta}_{\mu}^{\setminus\mu^{(t)}},y_{\mu})$
and
$(\partial_{\hat{\theta}}g_{\mathrm{out}}^{(t)})_{\mu}=\partial_{\hat{\theta}}g_{\mathrm{out}}(s_{\theta_{\mu}}^{\setminus\mu^{(t)}},\widehat{\theta}_{\mu}^{\setminus\mu^{(t)}},y_{\mu})$.

From (C.2), we obtain

$\displaystyle\mathbb{M}(\Sigma^{(t)}_{i},\mathrm{m}^{(t)}_{i})$
$\displaystyle=\mathbb{M}(\Sigma^{(t)}_{i\to\mu},\mathrm{m}^{(t)}_{i\to\mu})$
$\displaystyle+\frac{\partial}{\partial
m^{(t)}_{i\to\mu}}\mathbb{M}(\Sigma^{(t)}_{i\to\mu},\mathrm{m}^{(t)}_{i\to\mu})\frac{1}{\sqrt{N}}\Sigma^{(t)}_{i\to\mu}B_{\mu\to
i}^{(t)}+O(N^{-1})$
$\displaystyle=\mathbb{M}(\Sigma^{(t)}_{i\to\mu},\mathrm{m}^{(t)}_{i\to\mu})+\beta\mathbb{V}(\Sigma^{(t)}_{i},\mathrm{m}^{(t)}_{i})\frac{F_{\mu
i}}{\sqrt{N}}\beta^{-1}(g_{\mathrm{out}}^{(t)})_{\mu}+O(N^{-1}),$

where we use the relationship (4.13) and
${\Sigma}_{i}^{(t)}={\Sigma}_{i\to\mu}^{(t)}+O(N^{-1})$. Immediately, we get

$\displaystyle\widehat{x}_{i\to\mu}^{(t+1)}=\widehat{x}_{i}^{(t+1)}-\frac{1}{\sqrt{N}}{s}_{i}^{(t+1)}F_{\mu
i}\beta^{-1}g_{\mathrm{out}}(s_{\theta_{\mu}}^{\setminus\mu^{(t)}},\widehat{\theta}_{\mu}^{\setminus\mu^{(t)}},y_{\mu}),$
(C.26)

and $\widehat{x}_{i}^{(t+1)}$ and $s_{i}^{(t+1)}$ are obtained from
$\mathbb{M}_{\beta}(\Sigma_{i}^{(t)},\mathrm{m}_{i}^{(t)})$ and
$\beta\mathbb{V}_{\beta}(\Sigma_{i}^{(t)},\mathrm{m}_{i}^{(t)})$ where

$\displaystyle\mathrm{m}_{i}^{(t)}$
$\displaystyle={\Sigma}_{i}^{(t)}\sum_{\mu=1}^{M}\left\\{\frac{1}{\sqrt{N}}F_{\mu
i}\beta^{-1}(g_{\mathrm{out}}^{(t)})_{\mu}-\frac{1}{N}F_{\mu
i}^{2}\widehat{x}_{i}^{(t)}\beta^{-1}(\partial_{\widehat{\theta}}g_{\mathrm{out}}^{(t)})_{\mu}\right\\}$
(C.27) $\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}}^{(t)})_{\mu}.$ (C.28)

Further, at a sufficiently large $N$, the relationship between the mean and
variance of the parameter $\theta_{\mu}$ in the $\mu$-cavity and full system
are presented by

$\displaystyle\widehat{\theta}_{\mu}^{\setminus\mu^{(t)}}$
$\displaystyle=\widehat{\theta}_{\mu}^{(t)}-\beta^{-1}(g_{\mathrm{out}}^{(t-1)})_{\mu}{s_{\theta}}_{\mu}^{\setminus\mu^{(t)}}$
(C.29) $\displaystyle{s_{\theta}}_{\mu}^{\setminus\mu^{(t)}}$
$\displaystyle={s_{\theta}}_{\mu}^{(t)}=\frac{1}{N}\sum_{j=1}^{N}F^{2}_{\mu
j}{s}^{(t)}_{j}.$ (C.30)

##  Appendix D Coordinate descent algorithm

Coordinate descent (CD) is an algorithm that can optimise each variable
successively where other variables are fixed. The CD algorithm is efficient
for penalised regression problems, and partially for non-convex optimization
problems [78, 79].

###  D.1 Gaussian likelihood with penalty

We focus on the variable $x_{i}$ and optimise it under the fixed variables
$x_{j}=\widehat{x}_{j}$ where $j\in\bm{N}\setminus i$. For Gaussian
likelihood, we rewrite the function to be minimised as

$\displaystyle\varphi$
$\displaystyle=-\frac{1}{2}\sum_{\mu=1}^{M}\\!\left\\{-2\frac{y_{\mu}}{\sqrt{N}}F_{\mu
i}x_{i}\\!+\\!\frac{1}{N}F_{\mu i}^{2}x_{i}^{2}\\!+\\!\frac{2}{\sqrt{N}}F_{\mu
i}x_{i}\widehat{\theta}_{\mu}^{-i}\\!\right\\}\\!-\\!h(x_{i})\\!+\\!\mathrm{const.},$
(D.1)

where we set $\widehat{\theta}_{\mu}^{-i}=\frac{1}{\sqrt{N}}\sum_{j\neq
i}F_{\mu j}\widehat{x}_{j}$. In the case of the elastic net penalty, the
maximiser of (D.1) is given by

$\displaystyle\widehat{x}_{i}=\frac{b_{i}-\lambda_{1}\mathrm{sgn}(b_{i})}{\frac{1}{N}\sum_{\mu=1}^{M}F_{\mu
i}^{2}+\lambda_{2}}\mathbb{I}\left(|b_{i}|>\lambda_{1}\right)$ (D.2)

where

$\displaystyle
b_{i}=\frac{1}{\sqrt{N}}\sum_{\mu=1}^{M}(y_{\mu}-\widehat{\theta}_{\mu}^{-i})F_{\mu
i}.$ (D.3)

One can obtain the estimate by repeating this procedure for $i\in\bm{N}$
sufficient time steps.

###  D.2 Coordinate descent for GLM

For general likelihood, the proximal-Newton iterative approach can be combined
with the CD algorithm where the likelihood is repeatedly approximated by a
quadratic function [80, 35]. We obtain the following by expanding the log-
likelihood up to the second order around the estimate at step $t$, denoted by
$\widehat{\bm{x}}^{(t)}$:

$\displaystyle\ln f(\bm{y}|\bm{\theta}(\bm{F},\bm{x}))$
$\displaystyle\simeq\ln
f(\bm{y}|\bm{\theta}(\bm{F},\widehat{\bm{x}}^{(t)}))+\sum_{i=1}^{N}\frac{\partial}{\partial
x_{i}}\ln
f(\bm{y}|\bm{\theta}(\bm{F},\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}^{(t)}}(x_{i}-\widehat{x}_{i}^{(t)})$
$\displaystyle+\frac{1}{2}\sum_{ij}\frac{\partial^{2}}{\partial x_{i}\partial
x_{j}}\ln
f(\bm{y}|\bm{\theta}(\bm{x}))\Big{|}_{\bm{x}=\widehat{\bm{x}}^{(t)}}(x_{i}-\widehat{x}_{i}^{(t)})(x_{j}-\widehat{x}_{j}^{(t)})$
$\displaystyle=\ln
f(\bm{y}|\bm{\theta}(\bm{F},\widehat{\bm{x}}^{(t)}))+\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{\mu=1}^{M}F_{\mu
i}(y_{\mu}-a^{\prime}(\theta_{\mu}^{(t)}))(x_{i}-\widehat{x}_{i}^{(t)})$
$\displaystyle-\frac{1}{2N}\sum_{ij}\sum_{\mu=1}^{M}F_{\mu i}F_{\mu
j}a^{\prime\prime}(\theta_{\mu}^{(t)})(x_{i}-\widehat{x}_{i}^{(t)})(x_{j}-\widehat{x}_{j}^{(t)}).$
(D.4)

Hence, the penalized likelihood is approximated by

$\displaystyle\widetilde{\varphi}^{(t)}(\bm{x})$
$\displaystyle=-\frac{1}{2}(\bm{x}-\widehat{\bm{x}}^{(t)}-{\bm{W}^{(t)}}^{-1}\bm{R}^{(t)})^{\top}\bm{W}(\bm{x}-\widehat{\bm{x}}^{(t)}-{\bm{W}^{(t)}}^{-1}\bm{R}^{(t)})$
(D.5) $\displaystyle\hskip 28.45274pt-h(\bm{x})+\mathrm{const.},$

where

$\displaystyle\bm{R}^{(t)}$
$\displaystyle=\frac{1}{\sqrt{N}}\bm{F}^{\top}(\bm{y}-a^{\prime}(\bm{\theta}(\widehat{\bm{x}}^{(t)})))$
(D.6) $\displaystyle\bm{W}^{(t)}$
$\displaystyle=\frac{1}{N}\bm{F}^{\top}\bm{D}(\bm{\theta}(\widehat{\bm{x}}^{(t)}))\bm{F},$
(D.7)

and $\bm{D}(\bm{\theta})\in\mathbb{R}^{M\times M}$ is the diagonal matrix
whose $(\mu,\mu)$ component corresponds to
$a^{\prime\prime}(\theta_{\mu}(\widehat{\bm{x}}^{(t)}))$. We solve (D.5) using
the CD algorithm. When the penalty is $\ell_{2}$ norm, the solution is given
by

$\displaystyle\hat{\bm{x}}^{(t+1)}=(\bm{W}^{(t)}+\lambda\bm{I}_{N})^{-1}\bm{W}^{(t)}(\widehat{\bm{x}}^{(t)}+{\bm{W}^{(t)}}^{-1}\bm{R}^{(t)}).$
(D.8)

##  Appendix E Replica method

###  E.1 Derivation of the saddle point equations for conjugate variables

From the extremization condition of the free energy density, we obtain saddle
point equations as

$\displaystyle\widehat{\Theta}$ $\displaystyle=-2\alpha\int D\zeta
D\nu\frac{\partial}{\partial Q}\widehat{f}_{\xi}^{(\mathrm{RS})*}(\nu,\zeta)$
(E.1) $\displaystyle\widehat{\chi}$ $\displaystyle=2\alpha\int D\zeta
D\nu\frac{\partial}{\partial\chi}\widehat{f}_{\xi}^{(\mathrm{RS})*}(\nu,\zeta)$
(E.2) $\displaystyle\widehat{\mu}$ $\displaystyle=\alpha\int D\zeta
D\nu\frac{\partial}{\partial m}\widehat{f}_{\xi}^{(\mathrm{RS})*}(\nu,\zeta),$
(E.3)

where

$\displaystyle\widehat{f}_{\xi}^{(\mathrm{RS})*}(\nu,\zeta)$
$\displaystyle=-\frac{1}{2\chi}\left(u_{\mathrm{RS}}^{*}(\nu,\zeta)-\frac{m\nu+\sqrt{Q\sigma_{T}^{2}-m^{2}}\zeta}{\sigma_{T}}\right)^{2}$
(E.4) $\displaystyle\hskip 85.35826pt+\ln f\left({\cal
G}\left(\sigma_{T}\nu\right)|u_{\mathrm{RS}}^{*}(\nu,\zeta)\right)$

and in case of GLM, $u_{\mathrm{RS}}^{*}(\nu,\zeta)$ satisfies

$\displaystyle\frac{u_{\mathrm{RS}}^{*}(\nu,\zeta)-(m\nu+\sqrt{Q\sigma_{T}^{2}-m^{2}}\zeta)\sigma_{T}^{-1}}{\chi}={\cal
G}\left(\sigma_{T}\nu\right)-a^{\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta)).$
(E.5)

For the derivation of the saddle point equations, it is convenient to derive
the expression of $\partial u_{\mathrm{RS}}^{*}(\nu,\zeta)/\partial\zeta$ and
$\partial u_{\mathrm{RS}}^{*}(\nu,\zeta)/\partial\nu$. From (E.5), $\partial
u_{\mathrm{RS}}^{*}(\nu,\zeta)/\partial\zeta$ satisfies

$\displaystyle\frac{1}{\chi}\left(\frac{\partial
u_{\mathrm{RS}}^{*}(\nu,\zeta)}{\partial\zeta}-\sqrt{\frac{Q\sigma_{T}^{2}-m^{2}}{\sigma_{T}^{2}}}\right)=-\frac{\partial
u_{\mathrm{RS}}^{*}(\nu,\zeta)}{\partial\zeta}a^{\prime\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta)).$
(E.6)

Solving this, we obtain

$\displaystyle\frac{\partial
u_{\mathrm{RS}}^{*}(\nu,\zeta)}{\partial\zeta}=\sqrt{\frac{Q\sigma_{T}^{2}-m^{2}}{\sigma_{T}^{2}}}\frac{1}{1+\chi
a^{\prime\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta))}.$ (E.7)

Similarly, $\partial u^{*}(\nu,\zeta)/\partial\nu$ satisfies

$\displaystyle\frac{1}{\chi}\left(\frac{\partial
u_{\mathrm{RS}}^{*}(\nu,\zeta)}{\partial\nu}-\frac{m}{\sigma_{T}}\right)=\frac{\partial{\cal
G}(\sigma_{T}\nu,\varepsilon)}{\partial\nu}-\frac{\partial
u_{\mathrm{RS}}^{*}(\nu,\zeta)}{\partial\nu}a^{\prime\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta)),$
(E.8)

and solving this, we obtain

$\displaystyle\frac{\partial
u_{\mathrm{RS}}^{*}(\nu,\zeta)}{\partial\nu}=\frac{1}{1+\chi
a^{\prime\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta))}\left(\chi\frac{\partial{\cal
G}(\sigma_{T}\nu,\varepsilon)}{\partial\nu}+\frac{m}{\sigma_{T}}\right).$
(E.9)

Applying (E.7) to the saddle point equation of $\widehat{\Theta}$,

$\displaystyle\int D\zeta\frac{\partial}{\partial Q}$
$\displaystyle\widehat{f}_{\xi}^{(\mathrm{RS})*}(\nu,\zeta)$ (E.10)
$\displaystyle=\int
D\zeta\frac{\zeta\sigma_{T}}{2\chi\sqrt{Q\sigma_{T}^{2}-m^{2}}}\left(u_{\mathrm{RS}}^{*}(\nu,\zeta)-\frac{m\nu+\sqrt{Q\sigma_{T}^{2}-m^{2}}\zeta}{\sigma_{T}}\right)$
$\displaystyle=\frac{\sigma_{T}}{2\chi\sqrt{Q\sigma_{T}^{2}-m^{2}}}\int
D\zeta\left(\frac{\partial
u_{\mathrm{RS}}^{*}(\nu,\zeta)}{\partial\zeta}-\frac{\sqrt{Q\sigma_{T}^{2}-m^{2}}}{\sigma_{T}}\right)$
$\displaystyle=\frac{1}{2\chi}\int D\zeta\left(\frac{1}{1+\chi
a^{\prime\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta))}-1\right),$

hence we obtain

$\displaystyle\widehat{\Theta}=\alpha\int D\zeta
D\nu\frac{a^{\prime\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta))}{1+\chi
a^{\prime\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta)))}.$ (E.11)

Next, for the saddle point equation of $\widehat{\chi}$,

$\displaystyle\frac{\partial}{\partial\chi}\widehat{f}_{\xi}^{(\mathrm{RS})*}(\nu,\zeta)=\frac{1}{\chi^{2}}\left(u_{\mathrm{RS}}^{*}(\nu,\zeta)-\frac{m\nu+\sqrt{Q\sigma_{T}^{2}-m^{2}}\zeta}{\sigma_{T}}\right)^{2},$
(E.12)

and using (E.5), we obtain

$\displaystyle\widehat{\chi}=\alpha\int D\zeta D\nu\left({\cal
G}\left(\sigma_{T}\nu\right)-a^{\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta))\right)^{2}.$
(E.13)

For the saddle point equation of $\widehat{\mu}$, we obtain

$\displaystyle\int D\zeta D\nu\frac{\partial}{\partial
m}\widehat{f}_{\xi}^{(\mathrm{RS})*}(\nu,\zeta)$
$\displaystyle=\frac{1}{\chi\sigma_{T}}\int D\zeta
D\nu\left(\nu-\frac{m\zeta}{\sqrt{Q\sigma_{T}^{2}-m^{2}}}\right)\left(u_{\mathrm{RS}}^{*}-\frac{m\nu+\sqrt{Q\sigma_{T}^{2}-m^{2}}\zeta}{\sigma_{T}}\right)$
$\displaystyle=\frac{1}{\chi\sigma_{x}}\int D\zeta D\nu\Big{\\{}\frac{\partial
u_{\mathrm{RS}}^{*}(\nu,\zeta)}{\partial\nu}-\frac{m}{\sigma_{T}}$
$\displaystyle\hskip
42.67912pt-\frac{m}{\sqrt{Q\sigma_{T}^{2}-m^{2}}}\left(\frac{\partial
u_{\mathrm{RS}}^{*}(\nu,\zeta)}{\partial\zeta}-\frac{\sqrt{Q\sigma_{T}^{2}-m^{2}}}{\sigma_{T}}\right)\Big{\\}},$
(E.14)

and substituting (E.7) and (E.9) into (E.14), the following equation is
obtained

$\displaystyle\widehat{\mu}=\alpha\int D\zeta D\nu\frac{1}{1+\chi
a^{\prime\prime}(u_{\mathrm{RS}}^{*}(\nu,\zeta))}\frac{\partial{\cal
G}(\sigma_{T}\nu)}{\partial(\sigma_{T}\nu)}.$ (E.15)

###  E.2 Simplified calculation for Gaussian noise

When the likelihood under consideration is the Gaussian distribution, the
calculation can be simplified. We define

$\displaystyle u^{(a)}_{\mu}=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}F_{\mu
i}(x_{i}^{(0)}-x_{i}^{(a)})+\varepsilon_{\mu}$ (E.16)

for $\mu\in\bm{M}$ and $a=1,\ldots,n$. These variables satisfy

$\displaystyle\mathbb{E}_{\bm{y}|q^{(ab)}}\left[u_{\mu}^{(a)}u_{\mu}^{(b)}\right]=q^{(ab)}-m^{(a)}-m^{(b)}+\sigma_{T}^{2}$
(E.17)

under the subshell, hence their joint distribution is given by

$\displaystyle P_{\bm{u}}(\widetilde{\bm{u}}|{\cal
V})=\frac{1}{\sqrt{(2\pi)^{n}|{\cal
V}|}}\exp\left(-\frac{1}{2}\widetilde{\bm{u}}{\cal
V}^{-1}\widetilde{\bm{u}}^{\top}\right).$ (E.18)

Combining with the likelihood, we obtain

$\displaystyle\Phi_{\beta}^{(E)}(n)=-\frac{1}{2}\ln|\bm{I}_{n}+\beta{\cal
V}|.$ (E.19)

Under the RS assumption, the diagonal components and offdiagonal components of
${\cal V}$, denoted by ${\cal V}_{\mathrm{d}}$ and ${\cal V}_{\mathrm{nd}}$,
respectively, are given by

$\displaystyle{\cal V}_{\mathrm{d}}$ $\displaystyle=Q-2m+\sigma_{T}^{2}$
(E.20) $\displaystyle{\cal V}_{\mathrm{nd}}$
$\displaystyle=q-2m+\sigma_{T}^{2}.$ (E.21)

Calculating the matrix determinant, we obtain

$\displaystyle\Psi_{\beta}^{(E,\mathrm{RS})}(n)\\!=\\!-\frac{\alpha}{2}\left[\ln\left\\{1+\frac{n\beta(q-2m+\sigma_{T}^{2})}{1+\beta(Q-q)}\right\\}+n\ln\left\\{1+\beta(Q-q)\right\\}\\!\right].$
(E.22)

Taking $n\to 0$ limit and $\beta\to\infty$ limit, the free energy density for
the Gaussian likelihood is derived as

$\displaystyle\Psi=-\widehat{\mu}m+\frac{\widehat{\Theta}Q-\widehat{\chi}\chi}{2}-\frac{\alpha(Q-2m+\sigma_{T}^{2})}{2(1+\chi)}+\frac{\widehat{\chi}+\widehat{\mu}^{2}\sigma_{x}^{2}}{2(\widehat{\Theta}+\lambda)}.$
(E.23)

###  E.3 1RSB assumption and AT instability

The local stability of the RS solution against a symmetry breaking
perturbation is examined by constructing the one-step replica symmetry
breaking (1RSB) solution. We introduce the 1RSB assumption with respect to the
saddle point as

$\displaystyle(q^{(ab)},\widehat{q}^{(ab)})$
$\displaystyle=\left\\{\begin{array}[]{ll}\left(Q,-\displaystyle\frac{\widehat{Q}}{2}\right)&{\rm
if}~{}a=b\\\ (q_{1},\widehat{q}_{1})&{\rm if}~{}a\neq b,~{}B(a)=B(b)\\\
(q_{0},\widehat{q}_{0})&{\rm if}~{}B(a)\neq B(b)\\\ \end{array},\right.$
(E.27) $\displaystyle(m^{(a)},\widehat{m}^{(a)})$
$\displaystyle=(m,\widehat{m})~{}{\rm for~{}any}~{}a,$ (E.28)

where $n$ replicas are separated into $n/\ell$ blocks that contain $\ell$
replicas. The index of the block that includes the $a$-th replica is denoted
by $B(a)$. For the integral of $\theta$ under the 1RSB assumption, we set

$\displaystyle\theta^{(a)}$
$\displaystyle=\sqrt{Q-q_{1}}w^{(a)}+\sqrt{\Delta_{q}}z_{1}^{(B(a))}+\sqrt{q_{0}}z_{0}$
(E.29) $\displaystyle\theta^{(0)}$
$\displaystyle=\sqrt{\sigma_{T}-\frac{m^{2}}{q_{0}}}v+\frac{m}{\sqrt{q_{0}}}z_{0},$
(E.30)

where $w^{(a)}~{}(a\in\bm{n})$, $z_{1}^{(B)}~{}(B\in\\{1,\ldots,n/\ell\\})$,
$v$, and $z_{0}$ are the independent Gaussian random variables with mean 0 and
variance 1. Utilizing these expressions, we obtain

$\displaystyle\Psi_{\beta}^{(C,\mathrm{1RSB})}$
$\displaystyle=-\widehat{m}m+\frac{\widehat{Q}Q+\widehat{q}_{1}q_{1}}{2}-\frac{\ell(\widehat{q}_{1}q_{1}-\widehat{q}_{0}q_{0})}{2}$
(E.31) $\displaystyle\Psi_{\beta}^{(E,\mathrm{1RSB})}$
$\displaystyle=\frac{1}{\ell}\int
Dz_{0}Dv\ln\widehat{\xi}_{\beta}^{(\mathrm{RSB})}(v,z_{0})$ (E.32)
$\displaystyle\Psi_{\beta}^{(S,\mathrm{1RSB})}$
$\displaystyle=\frac{1}{\ell}\int
Dz_{0}dx^{(0)}\phi^{(0)}(x^{(0)})\ln\xi_{\beta}^{(\mathrm{RSB})}(z_{0},x^{(0)}),$
(E.33)

where we set

$\displaystyle\widehat{\xi}_{\beta}^{(\mathrm{RSB})}(v,z_{0})=\int
Dz_{1}\left\\{\\!\int Dwf^{\beta}\\!\\!\left({\cal
G}\left(\\!\sqrt{\sigma_{T}^{2}\\!-\\!\frac{m^{2}}{q_{0}}}v\\!+\\!\frac{m}{\sqrt{q_{0}}}z\\!\right)\\!\Big{|}u\\!\right)\\!\\!\right\\}^{\ell}$
(E.34) $\displaystyle\xi_{\beta}^{(\mathrm{RSB})}(z_{0},x^{(0)})=\int Dz_{1}$
(E.35) $\displaystyle\times\left\\{\int
dx\phi^{\beta}(x)\exp\left(-\frac{\widehat{Q}+\widehat{q}_{1}}{2}x^{2}+(\sqrt{\widehat{q}_{1}-\widehat{q}_{0}}z_{1}+\sqrt{\hat{q}_{0}}z_{0}+\hat{m}x^{(0)})x\right)\right\\}^{\ell},$

and $u=\sqrt{Q-q_{1}}w+\sqrt{q_{1}-q_{0}}z_{1}+\sqrt{q_{0}}z_{0}$. Using $u$,
we apply the expression
$w=(u-\sqrt{q_{1}-q_{0}}z_{1}-\sqrt{q_{0}}z_{0})/\sqrt{Q-q_{1}}$, and further,
we apply following transformation from $(v,z)$ to $(\nu,\zeta)$ that keeps
$dvdz=d\nu d\zeta$ as

$\displaystyle\sqrt{\sigma_{T}^{2}-\frac{m^{2}}{q_{0}}}v+\frac{m}{\sqrt{q_{0}}}z$
$\displaystyle=\sigma_{T}\nu$ (E.36)
$\displaystyle-\frac{m}{\sqrt{q_{0}}}v+\sqrt{\sigma_{T}^{2}-\frac{m^{2}}{q_{0}}}z$
$\displaystyle=\sigma_{T}\zeta,$ (E.37)

hence
$\sqrt{q_{0}}z=(m\nu+\sqrt{q_{0}\sigma_{T}^{2}-m^{2}}\zeta)\sigma_{T}^{-1}$.
At $\beta\to\infty$, the variables scale as $\beta(Q-q_{1})=\chi$,
$\beta\widehat{\Theta}=\widehat{Q}+\widehat{q}_{1}$,
$\beta^{2}\widehat{\chi}_{1}=\widehat{q}_{1}$,
$\beta^{2}\widehat{\chi}_{0}=\widehat{q}_{0}$,
$\beta\widehat{\mu}=\widehat{m}$, and $\widetilde{\ell}=\beta\ell$. With these
scaled variables, we obtain

$\displaystyle\widehat{\xi}_{\beta}^{\mathrm{(RSB)}}(\nu,\zeta)=\int
Dz_{1}\exp\left(\widetilde{\ell}\widehat{f}_{\xi}^{\mathrm{(RSB)}*}(z_{1},\nu,\zeta)\right),$
(E.38)

where we set $\Delta_{q}=q_{1}-q_{0}$,
$\widehat{f}_{\xi}^{\mathrm{(RSB)}*}(z_{1},\nu,\zeta)=\widehat{f}_{\xi}^{(\mathrm{RSB})}(u^{*}_{\mathrm{(RSB)}};z_{1},\nu,\zeta)$
with

$\displaystyle\widehat{f}_{\xi}^{(\mathrm{RSB})}(u;z_{1},\nu,\zeta)$
$\displaystyle=\frac{-(u-\sqrt{\Delta_{q}}z_{1}-(m\nu+\sqrt{q_{0}\sigma_{T}^{2}-m^{2}}\zeta)\sigma_{T}^{-1})^{2}}{2\chi}$
(E.39) $\displaystyle\hskip 142.26378pt+\ln f\left({\cal
G}\left(\sigma_{T}\nu\right)|u\right).$ $\displaystyle
u_{\mathrm{RSB}}^{*}(z_{1},\nu,\zeta)$
$\displaystyle=\mathop{\mathrm{argmax}}_{u}\widehat{f}^{\mathrm{(RSB)}}_{\xi}(u;z_{1},\nu,\zeta).$
(E.40)

In case of GLM, $u^{*}_{\mathrm{RSB}}$ satisfies

$\displaystyle\frac{u_{\mathrm{RSB}}^{*}-\sqrt{\Delta_{q}}z_{1}-(m\nu+\sqrt{q_{0}\sigma_{T}^{2}-m^{2}}\zeta)\sigma_{T}^{-1}}{\chi}={\cal
G}\left(\sigma_{T}\nu\right)-a^{\prime}(u^{*}_{\mathrm{RSB}}).$ (E.41)

Next, by applying the saddle point method to $\xi_{\beta}^{\mathrm{RSB}}$, we
obtain

$\displaystyle\xi_{\beta}^{\mathrm{(RSB)}}(z_{0},x^{(0)})=\int
Dz_{1}\exp\left(\widetilde{\ell}f_{\xi}^{\mathrm{(RSB)}*}(z_{1},z_{0},x^{(0)})\right),$
(E.42)

where we set $\widehat{\Delta}_{q}=\widehat{q}_{1}-\widehat{q}_{0}$, and
$f_{\xi}^{*}(z_{1},z_{0},x^{(0)})=f_{\xi}^{(\mathrm{RSB})}(x^{*}_{\mathrm{RSB}};z_{1},z_{0},x^{(0)})$
with

$\displaystyle f_{\xi}^{(\mathrm{RSB})}(x;z_{1},z_{0},x^{(0)})$
$\displaystyle=-\frac{\hat{\Theta}}{2}{x}^{2}+\\!(\sqrt{\hat{\Delta}_{q}}z_{1}+\\!\sqrt{\hat{\chi}_{0}}z+\\!\hat{\mu}x^{(0)})x+\\!\ln\phi(x)$
(E.43) $\displaystyle x_{\mathrm{RSB}}^{*}(z_{1},z_{0},x^{(0)})$
$\displaystyle=\mathop{\mathrm{argmax}}_{x}f_{\xi}^{(\mathrm{RSB})}(x;z_{1},z_{0},x^{(0)}).$
(E.44)

In summary, the free energy density under the 1RSB assumption is given by

$\displaystyle\Psi^{(\mathrm{1RSB})}$
$\displaystyle=\mathop{\mathrm{extr}}_{\Omega,\widehat{\Omega}}\Big{[}\widehat{\mu}m-\frac{\widehat{\Theta}Q-\widehat{\chi}_{1}\chi}{2}+\frac{\widetilde{\ell}(\widehat{\chi}_{1}Q-\widehat{\chi}_{0}q_{0})}{2}$
(E.45) $\displaystyle-\frac{1}{\widetilde{\ell}}\int
dx^{(0)}\phi^{(0)}(x^{(0)})Dz_{0}\ln\xi_{\beta}^{(\mathrm{RSB})}(z_{0},x^{(0)})$
$\displaystyle-\frac{\alpha}{\widetilde{\ell}}\int D\nu
D\zeta\ln\widehat{\xi}_{\beta}^{(\mathrm{RSB})}(\nu,\zeta)\Big{]},$

where $\Omega=\\{Q,\chi,q_{0},m\\}$ and
$\widehat{\Omega}=\\{\widehat{\Theta},\widehat{\chi}_{1},\widehat{\chi}_{0},\widehat{\mu}\\}$.
The saddle point equations are derived as

$\displaystyle Q$ $\displaystyle=\int
dx^{(0)}\phi^{(0)}Dz_{0}\langle{x^{*}_{\mathrm{RSB}}}^{2}\rangle_{\mathrm{RSB}}$
(E.46) $\displaystyle\chi$
$\displaystyle=\frac{1}{\sqrt{\widehat{\Delta}_{q}}}\int
dx^{(0)}\phi^{(0)}Dz_{0}\left\langle\frac{\partial
x^{*}_{\mathrm{RSB}}}{\partial z_{1}}\right\rangle_{\mathrm{RSB}}$ (E.47)
$\displaystyle q_{0}$ $\displaystyle=\int
Dz_{0}dx^{(0)}\phi^{(0)}(x^{(0)})\langle{x^{*}_{\mathrm{RSB}}}\rangle_{\mathrm{RSB}}^{2}$
(E.48) $\displaystyle m$ $\displaystyle=\int
Dz_{0}dx^{(0)}\phi^{(0)}(x^{(0)})x^{(0)}\langle{x^{*}_{\mathrm{RSB}}}\rangle_{\mathrm{RSB}}$
(E.49) $\displaystyle\widehat{\Theta}$ $\displaystyle=\alpha\int D\nu
D\zeta\left\langle\frac{a^{\prime\prime}(u^{*}_{\mathrm{RSB}})}{1+\chi
a^{\prime\prime}(u^{*}_{\mathrm{RSB}})}\right\rangle_{\widehat{\mathrm{RSB}}}$
(E.50) $\displaystyle\widehat{\chi}_{1}$ $\displaystyle=\alpha\int D\nu
D\zeta\left\langle\left({\cal
G}\left(\sigma_{T}\nu\right)-a^{\prime}(u^{*}_{\mathrm{RSB}})\right)^{2}\right\rangle_{\widehat{\mathrm{RSB}}}$
(E.51) $\displaystyle\widehat{\chi}_{0}$ $\displaystyle=\alpha\int D\nu
D\zeta\left\langle{\cal
G}\left(\sigma_{T}\nu\right)-a^{\prime}(u^{*}_{\mathrm{RSB}})\right\rangle^{2}_{\widehat{\mathrm{RSB}}}$
(E.52) $\displaystyle\widehat{\mu}$ $\displaystyle=\alpha\int D\nu
D\zeta\frac{\partial{\cal
G}(\sigma_{T}\nu)}{\partial(\sigma_{T}\nu)}\left\langle\frac{1}{1+\chi
a^{\prime\prime}(u^{*}_{\mathrm{RSB}})}\right\rangle_{\widehat{\mathrm{RSB}}},$
(E.53)

where we define

$\displaystyle\langle\varphi(z_{1})\rangle_{\mathrm{RSB}}$
$\displaystyle=\frac{\int
Dz_{1}\varphi(z_{1})\exp(\widetilde{\ell}f_{\xi}^{(\mathrm{RSB})*})}{\int
Dz_{1}\exp(\widetilde{\ell}f_{\xi}^{(\mathrm{RSB})*})}$ (E.54)
$\displaystyle\langle\varphi(z_{1})\rangle_{\widehat{\mathrm{RSB}}}$
$\displaystyle=\frac{\int
Dz_{1}\varphi(z_{1})\exp(\widetilde{\ell}\widehat{f}_{\xi}^{(\mathrm{RSB})*})}{\int
Dz_{1}\exp(\widetilde{\ell}\widehat{f}_{\xi}^{(\mathrm{RSB})*})}.$ (E.55)

When $\Delta_{q}=\widehat{\Delta}_{q}=0$, the 1RSB solution is reduced to the
RS solution. Hence, the stability of the RS solution is examined by the
stability analysis of the 1RSB solution of
$\Delta_{q}=\widehat{\Delta}_{q}=0$. We apply Taylor expansion to them by
assuming that they are sufficiently small.

For the expansion of $\Delta_{q}$ around $\widehat{\Delta}_{q}=0$ and that of
$\widehat{\Delta}_{q}$ around $\Delta_{q}$, $x^{*}_{\mathrm{RSB}}$ and
$u^{*}_{\mathrm{RSB}}$ for finite $\beta$ are useful. We denote them as

$\displaystyle\langle x\rangle=\frac{\int dxx\exp(-\beta
f_{\xi}^{(\mathrm{RSB})}(x;z_{1},z_{0},x^{(0)}))}{\int dx\exp(-\beta
f_{\xi}^{(\mathrm{RSB})}(x;z_{1},z_{0},x^{(0)}))}$ (E.56)
$\displaystyle\langle u\rangle=\frac{\int
duu\exp(-\beta\widehat{f}_{\xi}^{(\mathrm{RSB})}(u;z_{1},\nu,\zeta))}{\int
du\exp(-\beta\widehat{f}_{\xi}^{(\mathrm{RSB})}(u;z_{1},\nu,\zeta))},$ (E.57)

respectively. Using eqs. (E.56) and (E.57), we obtain

$\displaystyle\Delta_{q}$ $\displaystyle=\int
Dzdx^{(0)}\phi^{(0)}(x^{(0)})\beta^{2}\left(\langle x^{2}\rangle-\langle
x\rangle^{2}\right)^{2}\widehat{\Delta}_{q}$ (E.58)
$\displaystyle\widehat{\Delta}_{q}$ $\displaystyle=\alpha\int D\nu
D\zeta\beta^{2}\left(\left\langle({\cal
G}\left(\sigma_{T}\nu\right)-a^{\prime}(u))^{2}\right\rangle-\langle{\cal
G}\left(\sigma_{T}\nu\right)-a^{\prime}(u)\rangle^{2}\right)^{2}{\Delta}_{q}.$
(E.59)

These expression corresponds to

$\displaystyle\Delta_{q}$ $\displaystyle=\int
Dzdx^{(0)}\phi^{(0)}(x^{(0)})\left(\frac{\partial
x^{*}_{\mathrm{RSB}}}{\partial\sqrt{\widehat{\chi}_{0}}z}\right)^{2}\widehat{\Delta}_{q}$
(E.60) $\displaystyle\widehat{\Delta}_{q}$ $\displaystyle=\alpha\int D\nu
D\zeta\left(\frac{\partial}{\partial\sqrt{q_{0}-\frac{m^{2}}{\sigma_{T}^{2}}}\zeta}\left({\cal
G}\left(\sigma_{T}\nu\right)-a^{\prime}(u^{*}_{\mathrm{RSB}})\right)\right)^{2}{\Delta}_{q},$
(E.61)

hence the instability condition of $\Delta_{q}=\widehat{\Delta}_{q}=0$ is
given by (6.71).

## References

## References

  * [1] Taddy M 2019 The Economics of Artificial Intelligence: An Agenda (University of Chicago Press) chap The Technological Elements of Artificial Intelligence, pp 61–87
  * [2] The Event Horizon Telescope Collaboration et al 2019 The Astrophysical Journal Letters 875 L1
  * [3] Best M and Neuhauser D 2006 Qual Saf Health Care 15 142–143
  * [4] Akaike H 1973 Information theory and an extension of the maximum likelihood principle Second International Symposium on Information Theory (Alademiai Kiado, Budapest) pp 267–281
  * [5] Stone M 1974 Journal of the Royal Statistical Society, Series B 36 111–147
  * [6] Amari S 2016 Information Geometry and Its Applications (Springer)
  * [7] Hastie T, Montanari A, Rosset S and Tibshirani R J 2022 The Annals of Statistics 50 949–986
  * [8] Hochreiter S and Schmidhuber J 1997 Neural Computation 9 1–42
  * [9] Thomas V, Pedregosa F, van MerriÃ«nboer B, Mangazol P A, Bengio Y and Le Roux N 2020 On the interplay between noise and curvature and its effect on optimization and generalization International Conference on Artificial Intelligence and Statistics (PMLR) pp 3503–3513
  * [10] Ruján P 1993 Phys. Rev. Lett. 70 2968–2971
  * [11] Iba Y 1999 Journal of Physics A: Mathematical and General 32 3875
  * [12] Nishimori H 2001 Statistical physics of spin glasses and information processing: an introduction 111 (Clarendon Press)
  * [13] Kabashima Y, Wadayama T and Tanaka T 2009 Journal of Statistical Mechanics: Theory and Experiment 2009 L09003
  * [14] Krzakala F, Mézard M, Sausset F, Sun Y and Zdeborová L 2012 Journal of Statistical Mechanics: Theory and Experiment 2012 P08009
  * [15] Hastie T, Tibshirani R and Friedman J 2009 The Elements of Statistical Learning (Springer Science & Business Media)
  * [16] Efron B 2004 J. Amer. Statist. Assoc. 99 619–632
  * [17] Mallows C 1973 Technometrics 661–675
  * [18] Watanabe S 2009 Algebraic geometry and statistical learning theory 25 (Cambridge university press)
  * [19] Konishi S and Kitagawa G 2008 Information Criteria and Statistical Modeling (Springer)
  * [20] Watanabe S 2021 Jpn J Stat Data Sci 4 1–19
  * [21] Goodd I J and Gaskins R A 1971 Biometrika 58 255–277
  * [22] Hoaglin D C and Welsch R E 1978 The American Statistician 32 17–22
  * [23] Allen D M 1974 technometrics 16 125–127
  * [24] Cook R D 1977 Technometrics 19 15–18
  * [25] Cook R D 1979 Journal of the American Statistical Association 74 169–174
  * [26] Sherman J and Morrison W J 1950 The Annals of Mathematical Statistics 21 124–127
  * [27] Seber G A and Lee A J 2003 Linear regression analysis (John Wiley & Sons)
  * [28] Cox D R and Hinkley D V 1979 Theoretical statistics (CRC Press)
  * [29] Ripley B D 2007 Pattern recognition and neural networks (Cambridge university press)
  * [30] Stone M 1977 Journal of the Royal Statistical Society, Series B 39 44–47
  * [31] Hastie T and Tibshirani R 1990 Generalized Adaptive Models (Chapman and Hall (London))
  * [32] Stein C 1981 Annal. Stat. 9 1135–1151
  * [33] Ye J 1998 Journal of the American Statistical Association 93 120–131
  * [34] Zou H, Hastie T and Tibshirani R 2007 Annal. Stat. 35 2173–2192
  * [35] Hastie T, Tibshirani R and Wainwright M 2015 Statistical Learning with Sparsity: The Lasso and Generalizations (Chapman & Hall/CRC)
  * [36] Sakata A 2016  2016 123302
  * [37] Gelfand A E, Dey D K and Chang H 1992 Model determination using predictive distributions with implementation via sampling-based methods Tech. rep. Stanford Univ CA Dept of Statistics
  * [38] Gelman A, Hwang J and Vehtari A 2014 Stat. Comput. 24 997–1016
  * [39] Watanabe S 2010 Journal of Machine Learning Research 11 3571–3594
  * [40] Iba Y and Yano K Posterior covariance information criterion for arbitrary loss functions arXiv:2206.05887
  * [41] Watanabe S 2010 IEICE transactions on fundamentals of electronics, communications and computer sciences 93 617–626
  * [42] Belkin M, Hsu D, Ma S and Mandal S 2019 Proceedings of the National Academy of Sciences 116 15849–15854
  * [43] Opper M, Kinzel W, Kleinz J and Nehl R 1990 Journal of Physics A: Mathematical and General 23 L581
  * [44] Krogh A and Hertz J A 1992 Journal of Physics A: Mathematical and General 25 1135
  * [45] Opper M and Kinzel W 1996 Statistical mechanics of generalization Models of neural networks III (Springer) pp 151–209
  * [46] Advani M S, Saxe A M and Sompolinsky H 2020 Neural Networks 132 428–446
  * [47] Geiger M, Spigler S, d’Ascoli S, Sagun L, Baity-Jesi M, Biroli G and Wyart M 2019 Physical Review E 100 012115
  * [48] Mei S and Montanari A The generalization error of random features regression: Precise asymptotics and double descent curve arXiv:1908.05355
  * [49] Gerace F, Loureiro B, Krzakala F, Mézard M and Zdeborová L 2020 Generalisation error in learning with random features and the hidden manifold model International Conference on Machine Learning (PMLR) pp 3452–3462
  * [50] d’Ascoli S, Refinetti M, Biroli G and Krzakala F 2020 Double trouble in double descent: Bias and variance (s) in the lazy regime International Conference on Machine Learning (PMLR) pp 2280–2290
  * [51] Loureiro B, Gerbelot C, Refinetti M, Sicuro G and Krzakala F 2022 Fluctuations, bias, variance & ensemble of learners: Exact asymptotics for convex losses in high-dimension Proceedings of the 39th International Conference on Machine Learning (Proceedings of Machine Learning Research vol 162) ed Chaudhuri K, Jegelka S, Song L, Szepesvari C, Niu G and Sabato S (PMLR) pp 14283–14314
  * [52] Rangan S 2011 Generalized approximate message passing for estimation with random linear mixing Information Theory Proceedings (ISIT), 2011 IEEE International Symposium on (IEEE) pp 2168–2172
  * [53] Sakata A 2018  2018 063404
  * [54] Pearl J 1988 Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference (Morgan Kaufmann)
  * [55] Koller D and Friedman N 2009 Probabilistic Graphical Models: Principles and Techniques (The MIT Press)
  * [56] Richardson T J and Urbanke R L 2001 IEEE Transactions on Information Theory 47 599–618
  * [57] Kabashima Y 2003 J. Phys. A: Math. Gen. 36 11111–11121
  * [58] Donoho D L, Maleki A and Montanari A 2009 Proc. Natl. Acad. Sci. U.S.A 106 18914–18919
  * [59] Mézard M, Parisi G and Virasoro M 1987 Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications (World Scientific Publishing Co Inc)
  * [60] Thouless D J, Anderson P W and Palmer R G 1977 Philosophical Magazine 35 593–601
  * [61] Mezard M and Montanari A 2009 Information, physics, and computation (Oxford University Press)
  * [62] Gerbelot C, Abbara A and Krzakala F 2020 Asymptotic errors for high-dimensional convex penalized linear regression beyond gaussian matrices Proceedings of Machine Learning Research (PMLR) pp 1682–1713
  * [63] Kabashima Y, Krzakala F, Mézard M, Sakata A and Zdeborová L 2016 IEEE Transactions on Information Theory 62 4228–4265
  * [64] Zdeborová L and Krzakala F 2016 Advances in Physics 65 453–552
  * [65] Loureiro B, Gerbelot C, Cui H, Goldt S, Krzakala F, Mezard M and Zdeborová L 2021 Advances in Neural Information Processing Systems 34 18137–18151
  * [66] Kato K 2009 Journal of Multivariate Analysis 100 1338–1352
  * [67] Mazumder R, Friedman J H and Hastie T 2011 Journal of the American Statistical Association 106 1125–1138
  * [68] Friedman J H 2012 International Journal of Forecasting 28 722–738
  * [69] Hirose K, Tateishi S and Konishi S 2013 Computational Statistics & Data Analysis 59 28–40
  * [70] Opper M and Winther O 2655–2684 Neural Computation 12 2000
  * [71] De Almeida J and Thouless D J 1978 Journal of Physics A: Mathematical and General 11 983
  * [72] Obuchi T and Kabashima Y 2016 Journal of Statistical Mechanics: Theory and Experiment 2016 053304 URL https://doi.org/10.1088/1742-5468/2016/05/053304
  * [73] Javanmard A and Montanari A 2013 Information and Inference: A Journal of the IMA 2 115–144
  * [74] Barbier J, Krzakala F, Macris N, Miolane L and Zdeborová L 2019 Proceedings of the National Academy of Sciences 116 5451–5460
  * [75] Goldt S, Loureiro B, Reeves G, Krzakala F, Mézard M and Zdeborová L 2022 The gaussian equivalence of generative models for learning with shallow neural networks Mathematical and Scientific Machine Learning (PMLR) pp 426–471
  * [76] Montanari A and Saeed B N 2022 Universality of empirical risk minimization Conference on Learning Theory (PMLR) pp 4310–4312
  * [77] Gerace F, Krzakala F, Loureiro B, Stephan L and Zdeborová L 2022 arXiv preprint arXiv:2205.13303
  * [78] Breheny P and Huang J 2011 The annals of applied statistics 5 232
  * [79] Sakata A and Xu Y 2018  2018 033404
  * [80] Lee J D, Sun Y and Saunders M A 2014 SIAM Journal on Optimization 24 1420–1443