A Generalization Gap Estimation for Overparameterized Models via the Langevin Functional Variance

Abstract

This article discusses the estimation of the generalization gap, the difference between generalization performance and training performance, for overparameterized models including neural networks. We first show that a functional variance, a key concept in defining a widely-applicable information criterion, characterizes the generalization gap even in overparameterized settings where a conventional theory cannot be applied. As the computational cost of the functional variance is expensive for the overparameterized models, we propose an efficient approximation of the function variance, the Langevin approximation of the functional variance (Langevin FV). This method leverages only the first-order gradient of the squared loss function, without referencing the second-order gradient; this ensures that the computation is efficient and the implementation is consistent with gradient-based optimization algorithms. We demonstrate the Langevin FV numerically by estimating the generalization gaps of overparameterized linear regression and nonlinear neural network models, containing more than a thousand of parameters therein. Supplementary materials for this article are available online.

KEYWORD:

• Bayesian inference
• Langevin dynamics
• Neural networks
• Overparameterization

1 Introduction

The great success of deep neural networks (LeCun, Bengio, and Hinton 2015; Goodfellow, Bengio, and Courville 2016) has been reported in many applied fields, such as natural language processing, image processing, and the natural sciences. This has now altered our understanding of generalization in data science. Generalization is our model’s ability to adapt properly to new data drawn from the same distribution as that of the training data. A classical discipline for generalization is the bias-variance dilemma, that is: the richness of models reduces the modeling bias but suffers from fitting spurious patterns, leading to poor generalization performance. However, the modern practice of deep neural networks yields counterexamples to this discipline, that is: deep neural networks are capable of exactly fitting the training data (known as interpolation) through the use of an overabundance of parameters of the models (known as over-parameterization) with accurately predicted test data (Chiyuan et al. 2017; Belkin, Hsu, and Mitra 2018). Several such surprising phenomena have been reported: the double descent phenomenon (Belkin, Hsu, and Mitra 2018; Mario et al. 2020; Nakada and Imaizumi 2021; Hastie et al. 2022), and the multiple descent phenomenon (Adlam and Pennington 2020; d’ Ascoli, Sagun, and Biroli 2020). These interesting but controversial phenomena give us a chance to rethink generalization in modern data science.

From a practical point of view, estimating the generalization performance helps us understand what is actually happening. A naive estimator of the generalization performance is the training performance. However, real-case studies using deep neural networks (Azulay and Weiss 2019; Zhang et al. 2021) suggest that there is a gap between the generalization performance and the training performance, a generalization gap. Double descent phenomena also imply the existence of non-negligible generalization gap. These empirical evidences advocate the importance of accurately estimating the generalization gap.

Statistical tools for estimating the generalization gap have been reported in the literature on information criteria or optimism estimation. The well-known Akaike information criterion (AIC; Akaike 1974) estimates the generalization log-loss of the plug-in predictive density using a maximum likelihood estimator. The Takeuchi information criterion (TIC; Takeuchi 1976) is a modification of AIC to deal with model misspecification. The regularization information criterion (RIC; Shibata 1989) is an extension of TIC to accommodate maximum penalized likelihood estimation. Moody (1992) and Murata, Yoshizawa, and Amari (1994) generalized RIC to an arbitrary loss and advocated its use in nonlinear systems, such as neural networks. Mallows’ C_p (Mallows 1973) and Stein’s unbiased risk estimates (SURE; Stein 1981) offer an elegant estimation scheme of the generalization gap using the covariance between a predictor and its outcome in Gaussian models with the ${\mathcal{l}}_2$ loss. Ramani, Blu, and Unser (2008) proposed an efficient Monte Carlo sampling-based method (Monte Carlo SURE) to estimate the covariance. Recently, new methods associated with Bayesian learning have been developed. The deviance information criterion (DIC; Spiegelhalter et al. 2002) and widely-applicable information criterion (WAIC; Watanabe 2010) offer computationally efficient devices that estimate the generalization loss of Bayesian learning. In particular, Watanabe (2010) showed that in a statistical model with a fixed dimension, the generalization gap of Bayesian learning is asymptotically equal to the functional variance (FV), that is, the posterior variance of the log-likelihood.

However, studies on the data-driven measurement of the generalization gap in the overparameterized regime has been relatively scarce. Gao and Jojic (2016) extended Monte Carlo SURE to analyze the generalization gap in deep neural networks. Thomas et al. (2020) modified TIC for the aforementioned purpose. These studies empirically investigated the gaps in common datasets, such as MNIST and CIFAR-10, and successfully developed useful tools to understand the generalization gap in the overparameterized regime. However, these approaches lacked theoretical guarantees in overparameterized models although they provided guarantees in statistical models with a fixed dimension. This results in unclarified theoretical applicability of these approaches. Furthermore, the computational costs of these approaches related to memory and speed are relatively high. This is because Gao and Jojic (2016) requires training several times, and Thomas et al. (2020) needs the second-order gradient (Hessian) of the loss function.

In this article, we focus on yet another tool to measure the generalization gap, the Functional Variance (FV) developed in Watanabe (2010). We present the following contributions:

• Theoretical applicability: We prove that FV is asymptotically unbiased to the generalization gap of Bayesian learning for overparameterized linear regression models.

• Computational efficiency: We propose a computationally-efficient approximation of FV, Langevin FV (LFV), by leveraging only the first-order gradient of the loss function.

We show the theoretical applicability of FV as a measure of the generalization gap in overparameterized regime. Our theory employs the overparameterized linear regression model (see Belkin, Hsu, and Xu 2020; Bartlett et al. 2020; Hastie et al. 2022), but we should mention that this linear model can be regarded as a linear approximation of nonlinear overparameterized models including deep neural networks and the linearization has been widely employed to analyze behaviors of deep neural networks (see, Jacot, Gabriel, and Hongler 2018; Arora et al. 2019; Chizat, Oyallon, and Bach 2019; Yang and Littewin 2021). One of such theories is known as neural tangent kernel (Jacot, Gabriel, and Hongler 2018), that is, an approximation of neural networks in a small region of the parameters around an estimate or an initial value. The condition of the neural tangent kernel approximation has been investigated in Chizat, Oyallon, and Bach (2019) and known to be applicable to any architecture of deep neural networks in Yang and Littewin (2021).

Although FV is theoretically favorable, it is defined by using the full posterior covariance and the full Bayesian inference in the overparameterized regime is often prohibited because of the computational burden. For efficiently computing the posterior, we employ the Langevin dynamics (see Risken 1996), which sequentially adds a normal random perturbation to each update of the gradient descent optimization and obtains the stationary distribution approximating the posterior distribution (Cheng et al. 2018). This approach has several merits from the computational perspective. For example, the implementation is easy and consistent with the gradient-based optimization that is the de-facto standard in deep neural network applications. Since the effect of linearization in theory remains unclear and controversial in practice (see Lee et al. 2020), we confirm the applicability of the developed method to regression using neural networks with real datasets in Section 4.3.

The rest of the article is organized as follows. Section 2 provides our theoretical results on FV, Section 3 proposes Langevin FV, which is an efficient implementation of FV using Langevin dynamics, Section 4 demonstrates the numerical experiments, and finally, Section 5 concludes this article. Supplementary material gives proofs of all theorems, and descriptions of source codes to reproduce experimental results.

2 Theoretical Results for the Functional Variance

In this section, we present theoretical results for the functional variance in overparameterized linear regression models. Throughout this article, ${\mathit{I}}_n$ denotes the n × n identity matrix for any $n\in \mathbb{N}$, and $\Vert \mathit{\beta }{\Vert }_2={({\beta }_1^2+{\beta }_2^2+\cdots +{\beta }_p^2)}^{1/2}$ and $\Vert \mathit{\beta }{\Vert }_{\infty }={\max }_{i=1,2,\dots ,p}\left|{\beta }_i\right|$ for any vector $\mathit{\beta }={({\beta }_1,{\beta }_2,\dots ,{\beta }_p)}^{\top }\in {\mathbb{R}}^p$.

2.1 Problem Setting

We begin by introducing the problem setting for the theory. We consider a linear regression model $$\mathit{y}=\mathit{X}{\mathit{\beta }}_0+\mathit{\epsilon },$$where $\mathit{y}={(y_1,\dots ,y_n)}^{\top }\in {\mathbb{R}}^n$ is a vector of observed outcomes, $\mathit{X}={({\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_n)}^{\top }\in {\mathbb{R}}^{n\times p}$ is an n × p random design matrix, ${\mathit{\beta }}_0\in {\mathbb{R}}^p$ is an unknown coefficient vector, and $\mathit{\epsilon }={({\epsilon }_1,\dots ,{\epsilon }_n)}^{\top }\in {\mathbb{R}}^n$ is a vector of iid error terms with mean zero and variance $0<{\sigma }_0^2<\infty$. Our interest is the overparameterized situation, where the number of regressors p is larger than the sample size n: $$n\le p.$$

We take the quasi-Bayesian approach on the vector $\mathit{\beta }$. Working under the quasi-likelihood $f(y_i\mid {\mathit{x}}_i,\mathit{\beta })$ with a Gaussian distribution with mean ${\mathit{x}}_i^{\top }\mathit{\beta }$ and variance ${\sigma }_0^2$, and a Gaussian prior $N_p(\mathit{0},\{{\sigma }_0^2/(\alpha n)\}{\mathit{I}}_p)$ on $\mathit{\epsilon }$ with mean $\mathit{0}$ and covariance matrix $\{{\sigma }_0^2/(\alpha n)\}{\mathit{I}}_p$, we obtain the quasi-posterior distribution of $\mathit{\beta }$: (1) $$\begin{array}{c}{\Pi }_{\alpha }(\mathrm{d}\mathit{\beta }\mid \mathit{y},\mathit{X})\\ \hspace{1em}={(2\pi )}^{-p/2}{(\text{det}{\mathit{Q}}_{\alpha })}^{-1/2}{\mathrm{e}}^{-{(\mathit{\beta }-{\widehat{\mathit{\beta }}}_{\alpha })}^{\top }{\mathit{Q}}_{\alpha }^{-1}(\mathit{\beta }-{\widehat{\mathit{\beta }}}_{\alpha })/2}\mathrm{d}\mathit{\beta },\end{array}$$(1) where ${\widehat{\mathit{\beta }}}_{\alpha }$ is the maximum a posterior estimate (2) $$\begin{array}{c}{\widehat{\mathit{\beta }}}_{\alpha }:={\left(n^{-1}{\mathit{X}}^{\top }\mathit{X}+\alpha {\mathit{I}}_p\right)}^{-1}{n}^{-1}{\mathit{X}}^{\top }\mathit{y}\hspace{1em}\\ \hspace{1em}=\underset{\mathit{\beta }\in {\mathbb{R}}^d}{\arg \min }\left\{{\mathcal{l}}_{\alpha }(\mathit{\beta }) := n^{-1}\Vert \mathit{y}-\mathit{X\beta }{\Vert }_2^2+\alpha \Vert \mathit{\beta }{\Vert }_2^2\right\}\end{array}$$(2) and ${\mathit{Q}}_{\alpha }$ is the matrix (3) $${\mathit{Q}}_{\alpha }:=n^{-1}{\sigma }_0^2{\left(n^{-1}{\mathit{X}}^{\top }\mathit{X}+\alpha {\mathit{I}}_p\right)}^{-1}.$$(3)

We now analyze the Gibbs generalization gap (4) $$\begin{array}{c}\Delta (\alpha ;\mathit{X}):=\frac{1}{2{\sigma }_0^2}\{{\mathbb{E}}_{\mathit{y},{\mathit{y}}^{*}}\left[{\mathbb{E}}_{\mathit{\beta }}\left[\Vert {\mathit{y}}^{*}-\mathit{X\beta }{\Vert }_2^2\right]\right]\\ \hspace{1em}-{\mathbb{E}}_{\mathit{y}}\left[{\mathbb{E}}_{\mathit{\beta }}\left[\Vert \mathit{y}-\mathit{X\beta }{\Vert }_2^2\right]\right]\},\end{array}$$(4) where ${\mathit{y}}^{*}$ is an independent copy of $\mathit{y}$ with given $\mathit{X}$, ${\mathbb{E}}_{\mathit{\beta }}[·]$ is the expectation with respect to the quasi-posterior distribution (1), and ${\mathbb{E}}_{\mathit{y}}[·]$ (and ${\mathbb{E}}_{\mathit{y},{\mathit{y}}^{*}}[·]$) is the expectation with respect to $\mathit{y}$ (and $\mathit{y},{\mathit{y}}^{*}$, respectively,) with given $\mathit{X}$. The Gibbs generalization gap considers the generalization gap for one stochastic sample from the quasi-posterior distribution as estimates of the parameters. To estimate the Gibbs generalization gap from the current observations, we focus on the functional variance (5) $$\text{FV}(\alpha ;\mathit{X}):=\sum _{i=1}^n{\mathbb{V}}_{\mathit{\beta }}[\log \hspace{0.17em}f(y_i\mid {\mathit{x}}_i,\mathit{\beta })],$$(5) where ${\mathbb{V}}_{\mathit{\beta }}$ is the variance with respect to the quasi-posterior distribution (1).

2.2 Asymptotic Unbiasedness of the Functional Variance

Here we present our theoretical findings on the functional variance in the overparameterized regime. The following theorem shows that the functional variance $\text{FV}(\alpha ;\mathit{X})$ is an asymptotically unbiased estimator of the Gibbs generalization gap $\Delta (\alpha ;\mathit{X})$ for the overparameterized linear regression with Gaussian covariates. Its proof is given in Supplement A.

Theorem 1.

Let $n,p\in \mathbb{N}$, and let $\Sigma$ be a p × p nonzero and nonnegative definite matrix. Let ${\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_n\in {\mathbb{R}}^p$ be iid random vectors from the Gaussian distribution with a mean zero and a covariance matrix $\Sigma$. There exists an absolute constant C > 0, such that we have, for any $\epsilon ,\alpha >0$, $${\mathbb{P}}_{\mathit{X}}\left(\left|{\mathbb{E}}_{\mathit{y}}\right[\text{FV}(\alpha ;\mathit{X})]-\Delta (\alpha ;\mathit{X})|>\epsilon \right)$$ $$\hspace{1em}\le C\left(\frac{{\xi }^2}{{\alpha }^2}\frac{1}{n\epsilon }+\frac{{\xi }^4}{{\alpha }^4}\frac{1}{n\epsilon }+\frac{{\xi }^4{b}^4}{{\alpha }^2{\sigma }_0^4}\frac{1}{n{\epsilon }^2}+\frac{1}{n}\right),$$where $\xi :=\text{tr}\left\{\Sigma \right\},b:=p^{1/2}\Vert {\mathit{\beta }}_0{\Vert }_{\infty }$, and ${\mathbb{P}}_{\mathit{X}}$ denotes the probability with respect to $\mathit{X}$. Furthermore, under the conditions (C1) $\underset{p\in \mathbb{N}}{\sup }\xi <\infty$ and (C2) $\underset{p\in \mathbb{N}}{\sup }b<\infty$, we have $${\mathbb{P}}_{\mathit{X}}(|{\mathbb{E}}_{\mathit{y}}[\text{FV}(\alpha ;\mathit{X})]-\Delta (\alpha ;\mathit{X})|\le C\prime a_n/\sqrt{n})\to 1$$ for some $C\prime >0$ not depending on n and p with an arbitrary slowly increasing sequence a_n.

Theorem 1 implies that FV successfully estimates the Gibbs generalization gap in the overparameterized regime, which supports the use of FV, even in the overparameterized regime. Interestingly, its rate of convergence $(a_n/\sqrt{n})$ is not affected by the dimension p but affected by the value of $\underset{p\in \mathbb{N}}{\sup }\xi$. As $\underset{p\in \mathbb{N}}{\sup }\xi$ gets larger, the decay of the difference becomes slower. Furthermore, the result does not restrict the true distribution of the additive error to the Gaussian distribution; in contrast, some classical theories such as SURE require the error term to either be Gaussian or be of related distributions.

Let us mention the conditions in the theorem. Condition (C1) indicates the trace boundedness of ${\mathit{X}}^{\top }\mathit{X}$ (divided by n) because $\Sigma =\mathbb{E}[{\mathit{X}}^{\top }\mathit{X}]/n$. Both shallow and deep neural network models satisfy this condition in general settings (see, e.g., Karakida, Akaho, and Amari 2019). Condition (C2) together with Condition (C1) indicates that $n^{-1}\Vert \mathit{\mu }{\Vert }_2^2\le n^{-1}\Vert \mathit{X}{\Vert }_{\mathrm{F}}^2\Vert {\mathit{\beta }}_0{\Vert }_2^2$ is upper-bounded by some constant C > 0 with probability approaching 1. It implies that the average of the entries in the outcome expectation $\mathit{\mu }$ is upper-bounded with high probability.

The main ingredient of the proof of Theorem 1 is the explicit identity of the difference between the generalization gap and FV. See Lemma 2. Its proof is presented in Supplement B. For two matrices $\mathit{A}=(a_{ij}),\mathit{B}=(b_{ij})$ of the same size, $\mathit{A}°\mathit{B}=(a_{ij}{b}_{ij})$ denotes the Hadamard element-wise product, and let ${\mathit{A}}^{\otimes 2}:=\mathit{A}{\mathit{A}}^{\top }$.

Lemma 2.

We have (6) $$\Delta (\alpha ;\mathit{X})=\text{tr}\left\{{\mathit{H}}_{\alpha }\right\}\hspace{1em}\text{and}$$(6) (7) $$\begin{array}{c}{\mathbb{E}}_{\mathit{y}}[\text{FV}(\alpha ;\mathit{X})]-\Delta (\alpha ;\mathit{X})=-\frac{3}{2}\text{tr}\{{\mathit{H}}_{\alpha }°{\mathit{H}}_{\alpha }\}+\text{tr}\{{\mathit{H}}_{\alpha }°{\mathit{H}}_{\alpha }^2\}\\ \hspace{1em}+\frac{1}{{\sigma }_0^2}\text{tr}\{{\mathit{H}}_{\alpha }°{(({\mathit{I}}_n-{\mathit{H}}_{\alpha })(\mathit{X}{\mathit{\beta }}_0))}^{\otimes 2}\},\end{array}$$(7)

where ${\mathit{H}}_{\alpha }$ is the regularized hat matrix (8) $${\mathit{H}}_{\alpha }:=n^{-1}\mathit{X}{\left(n^{-1}{\mathit{X}}^{\top }\mathit{X}+\alpha {\mathit{I}}_p\right)}^{-1}{\mathit{X}}^{\top }.$$(8)

Lemma 2 highlights the role of the regularized hat matrix ${\mathit{H}}_{\alpha }$ in evaluating the generalization gap and the residual. Same as the fixed-dimension theory, the hat matrix controls the magnitude of the generalization gap. With singular values $s_1\ge s_2\ge \cdots \ge s_n\ge 0$ of the matrix $\mathit{X}$, under the conditions in Theorem 1, we have $$\Delta (\alpha ;\mathit{X})=\text{tr}\left\{{\mathit{H}}_{\alpha }\right\}=\sum _{i=1}^n\frac{s_i^2/n}{(s_i^2/n)+\alpha }\to \exists {\Delta }_{\infty }(\alpha )\in [0,\infty )$$ $$(n\to \infty )$$as $\text{tr}\left\{{\mathit{H}}_{\alpha }\right\}\le {\alpha }^{-1}{n}^{-1}({\sum }_{i=1}^n{s}_i^2)<\infty$ with the probability approaching 1 as $n\to \infty$. Consider a simple case where $\mathit{X}$ has a fixed intrinsic dimension $p_{*}\ll n$ with growing n, p, namely $s_1\ge s_2\ge \cdots \ge s_{p_{*}}>0=s_{p_{*}+1}=s_{p_{*}+2}=\cdots =s_n$. Then, FV approaches the fixed intrinsic dimension $p_{*}$ (as $\alpha \to 0$), which coincides with the degrees of freedom (Mallows 1973; Ye 1998) and AIC (Akaike 1974) for linear regression equipped with $p_{*}$ regressors.

Remark 1

(Sample-wise and joint log-likelihoods). One may think the use of the posterior covariance of the joint log-likelihood $$\mathrm{J}-\text{FV}(\alpha ;\mathit{X}):={\mathbb{V}}_{\mathit{\beta }}[\log \hspace{0.17em}f(\mathit{y}\mid \mathit{X},\mathit{\beta })]$$instead of FV (5) using sample-wise log-likelihood. In this case, we have the following expression of the difference between the J-FV and the generalization gap. Its proof is in Supplement C.

Proposition 1.

For any $\alpha >0$, we have $${\mathbb{E}}_{\mathit{y}}[\mathrm{J}-\text{FV}(\alpha ;\mathit{X})]-\Delta (\alpha ;\mathit{X})$$ $$\hspace{1em}=-\frac{3}{2}\text{tr}\left\{{\mathit{H}}_{\alpha }^2\right\}+\text{tr}\left\{{\mathit{H}}_{\alpha }^3\right\}+\frac{1}{{\sigma }_0^2}\text{tr}\left\{{\mathit{H}}_{\alpha }{((\mathit{I}-{\mathit{H}}_{\alpha })(\mathit{X}{\mathit{\beta }}_0))}^{\otimes 2}\right\}.$$

This exhibits an interesting correspondence between FV and J-FV (see Lemma 2 and Proposition 1); replacing the Hadamard products in ${\mathbb{E}}_{\mathit{y}}[\text{FV}(\alpha ;\mathit{X})]$ with simple matrix products yields ${\mathbb{E}}_{\mathit{y}}[\mathrm{J}-\text{FV}(\alpha ;\mathit{X})]$. Since the simple matrix products in ${\mathbb{E}}_{\mathit{y}}[\mathrm{J}-\text{FV}(\alpha ;\mathit{X})]$ do not vanish as $n\to \infty$, J-FV is not an asymptotically unbiased estimator of the generalization gap.

Remark 2

(Extension of the theorem). Theorem 1 is further extended in Theorem 3, which proves the convergence of FV under milder conditions. Its proof is shown in Supplement A.

Theorem 3.

Let $n,p\in \mathbb{N}$ and assume the setting in Section 2.1. Let $\mathit{X}\in {\mathbb{R}}^{n\times p}$ be a random matrix and let ${\mathbb{S}}^{n-1}=\{\mathit{u}\in {\mathbb{R}}^d\mid \Vert \mathit{u}{\Vert }_2=1\}$ be a $(n-1)$-sphere. Let $q_1,q_2,\dots ,q_n$ be marginal probability densities of the left-singular vectors ${\mathit{u}}_1,{\mathit{u}}_2,\dots ,{\mathit{u}}_n\in {\mathbb{S}}^{n-1}$ of $\mathit{X}$. Write $$b:=p^{1/2}\Vert {\mathit{\beta }}_0{\Vert }_{\infty },\hspace{1em}\psi :={\max }_{i=1,2,\dots ,n}{\max }_{\mathit{u}\in {\mathbb{S}}^{n-1}}\{q_i(\mathit{u})\}{\int }_{{\mathbb{S}}^{n-1}}\mathrm{d}\mathit{v},$$ $$\eta :=\frac{1}{n}\text{tr}\left\{{\mathit{X}}^{\top }\mathit{X}\right\}.$$

Then, there exists an absolute constant C > 0 such that, for any $\epsilon ,\alpha >0$, we have (9) $$\begin{array}{c}{\mathbb{P}}_{\mathit{X}}\left(\left|{\mathbb{E}}_{\mathit{y}}\right[\text{FV}(\alpha ;\mathit{X})]-\Delta (\alpha ;\mathit{X})|>\epsilon \mid \eta \right)\\ \hspace{1em}\le C\left(\frac{\psi {\eta }^2}{{\alpha }^2}\frac{1}{n\epsilon }+\frac{\psi {\eta }^4}{{\alpha }^4}\frac{1}{n\epsilon }+\frac{\psi {\eta }^4{b}^4}{{\alpha }^2{\sigma }_0^4}\frac{1}{n{\epsilon }^2}\right),\end{array}$$(9) where ${\mathbb{P}}_{\mathit{X}}(·\mid \eta )$ denotes the conditional probability of $\mathit{X}$ given η.

Theorem 1 (Gaussian covariates) follows from Theorem 3 with ψ = 1; see Proposition 7.1 of Eaton (1989), which proves that the left-singular vectors of the Gaussian design matrix follow a uniform distribution over the unit $(n-1)$-sphere ${\mathbb{S}}^{n-1}$, that is, $q_i(\mathit{u})={\mathbb{1}}_{{\mathbb{S}}^{n-1}}(\mathit{u})/{\int }_{{\mathbb{S}}^{n-1}}\mathrm{d}\mathit{v}$ with ${\mathbb{1}}_{{\mathbb{S}}^{n-1}}(·)$ the indicator function of ${\mathbb{S}}^{n-1}$.

Remark 3

(Definition of the generalization gap). Regarding the definition of the generalization gap $\Delta (\alpha ;\mathit{X})$, we may consider another design matrix ${\mathit{X}}^{*}$ for generating ${\mathit{y}}^{*}$. However, the covariate difference is less effective in the generalization gap (4) since we focus on the prediction of the conditional random variable $\mathit{y}\mid \mathit{X}$ and not the covariate $\mathit{X}$. For theoretical simplicity, we employed a single design matrix $\mathit{X}$ to generate both outcomes $\mathit{y},{\mathit{y}}^{*}$.

3 Langevin Functional Variance

Despite FV being theoretically attractive as proved in Theorem 3, there exist two computational difficulties:

• (D1) generating samples from the quasi-posterior (1) requires the p × p matrix ${\mathit{Q}}_{\alpha }$, which is computationally intensive for overparameterized models $n\ll p$, and

• (D2) computing the quasi-posterior (1) is inconsistent with gradient-based optimization approaches, such as the stochastic gradient descent (SGD), which are often used in optimizing overparameterized models.

To resolve these difficulties (D1) and (D2), we consider a Langevin approximation of the quasi-posterior in Section 3.1, and propose Langevin FV (LFV) for linear models in Section 3.2; we then extend it to nonlinear models in Section 3.3. We compare FV and the proposed LFV with the existing estimators based on TIC (Takeuchi 1976) in Section 3.4.

3.1 Langevin Approximation of the Quasi-Posterior

A key idea of our algorithm is approximating the quasi-posterior (1) via a Langevin process. Starting with an estimator ${\mathit{\gamma }}^{(1)}:={\widehat{\mathit{\beta }}}_{\alpha }$, we stochastically update ${\mathit{\gamma }}^{(t)}$ by (10) $${\mathit{\gamma }}^{(t+1)}={\mathit{\gamma }}^{(t)}-\frac{1}{4}\delta {\kappa }_n\frac{\partial {\mathcal{l}}_{\alpha }({\mathit{\gamma }}^{(t)})}{\partial \mathit{\gamma }}+{\delta }^{1/2}{\mathit{e}}^{(t)}\hspace{1em}(t=1,2,\dots ),$$(10) where ${\kappa }_n:=n/{\sigma }_0^2$ and $\delta >0$ are user-specified parameters, let $\{{\mathit{e}}^{(1)},{\mathit{e}}^{(2)},\dots \}$ denote a sequence of iid standard normal random vectors, and ${\mathcal{l}}_{\alpha }$ is defined in (2). Then, the distribution of the Langevin process (10) approximates the quasi-posterior as discussed in the following. First, the Langevin process (10) is a discretization of the Ornstein-Uhlenbeck process $$\mathrm{d}{\tilde{\mathit{\gamma }}}_{\tau }=-\frac{1}{2}{\mathit{Q}}_{\alpha }^{-1}({\tilde{\mathit{\gamma }}}_{\tau }-{\widehat{\mathit{\beta }}}_{\alpha })\mathrm{d}\tau +\mathrm{d}{\mathit{e}}_{\tau }$$ equipped with a Wiener process ${\left\{{\mathit{e}}_{\tau }\right\}}_{\tau \ge 0}$, that is, ${\mathit{e}}_{\tau }-{\mathit{e}}_{\tau \prime }\sim N_p(\mathit{0},(\tau -\tau \prime ){\mathit{I}}_p)$ for any $\tau >\tau \prime \ge 0$. The distribution of ${\tilde{\mathit{\gamma }}}_{\tau }$ with the initialization ${\tilde{\mathit{\gamma }}}_0:={\widehat{\mathit{\beta }}}_{\alpha }$ coincides with the quasi-posterior (11) $$\mathcal{D}({\tilde{\mathit{\gamma }}}_{\tau })=N_p({\widehat{\mathit{\beta }}}_{\alpha },{\mathit{Q}}_{\alpha })\hspace{1em}(\tau >0).$$(11)

See Risken (1996) p. 156 for an example of the distribution (11). Next, with n and p fixed, for any sufficiently small $\varphi >0$, Theorem 2 of Cheng et al. (2018) evaluates the 1-Wasserstein distance between the distributions of the Ornstein-Uhlenbeck and the Langevin processes as (12) $$W_1(\mathcal{D}({\tilde{\mathit{\gamma }}}_{t\delta }),\mathcal{D}({\mathit{\gamma }}^{(t)}))\le \varphi$$(12) for some $\delta =O({\varphi }^2/p)$ and $t\asymp p/{\varphi }^2$; see Proposition 2 in Supplement D for details. Thus, the relations (11) and (12) imply that the distribution $\mathcal{D}({\mathit{\gamma }}^{(t)})$ approximates the quasi-posterior $N_p({\widehat{\mathit{\beta }}}_{\alpha },{\mathit{Q}}_{\alpha })$.

This Langevin approximation resolves difficulties (D1) and (D2) since

• (S1) the Langevin process (10) is computed with only the first order gradient of the squared loss function ${\mathcal{l}}_{\alpha }(\mathit{\gamma })$ defined in (2), meaning that the large matrix ${\mathit{Q}}_{\alpha }$ is not explicitly computed, and

• (S2) the Langevin process (10) is in the form of the gradient descent up to the normal noise ${\delta }^{1/2}{\mathit{e}}^{(t)}$; this process can be implemented consistently with gradient-based algorithms, which are often used to optimize overparameterized models.

Remark 4

(The other approach for the posterior approximation). Besides the Langevin dynamics, we can use the other approaches for approximating the posterior distribution. One such approach is to use the stochastic noise in stochastic gradient descent (Sato and Nakagawa 2014; Mandt, Hoffman, and Blei 2017). Though powerful, this approach relies on the normal approximation in the stochastic gradient and the normal approximation is hardly expected in the overparameterized regime, which results in the computational instability.

3.2 Langevin Functional Variance

Using the Langevin approximation (10), we propose LFV: for a time step $T\in \mathbb{N}$, let ${\left\{{\mathit{\gamma }}^{(t)}\right\}}_{t=1,2,\dots ,T}$ be the samples from the Langevin process. Let ${\widehat{\mathbb{V}}}_{\mathit{\gamma }}[·]$ denotes an empirical variance with respect to ${\left\{{\mathit{\gamma }}^{(t)}\right\}}_{t=1,2,\dots ,T}$. Let ${\mu }_i^{(t)}:={\mathit{x}}_i^{\top }{\mathit{\gamma }}^{(t)}$. Then, we define LFV as (13) $$\begin{array}{c}\text{LFV}(\alpha ;\mathit{X}):=\sum _{i=1}^n{\widehat{\mathbb{V}}}_{\mathit{\gamma }}[\log \hspace{0.17em}f(y_i\mid {\mathit{x}}_i,\mathit{\gamma })]\\ =\sum _{i=1}^n\frac{1}{T}\sum _{t=1}^T\{\frac{1}{2{\sigma }_0^2}{\left(y_i-{\mu }_i^{(t)}\right)}^2\\ \hspace{1em}-{\frac{1}{T}\sum _{t\prime =1}^T\frac{1}{2{\sigma }_0^2}{\left(y_i-{\mu }_i^{(t\prime )}\right)}^2\}}^2.\end{array}$$(13)

As shown in (12), the distribution of the Langevin process approaches the quasi-posterior: LFV is expected to approximate FV for a sufficiently large $T\in \mathbb{N}$. Thus, LFV is expected to inherit the nature of FV, such as the asymptotic unbiasedness when estimating the generalization gap.

The distributional approximation of the high-dimensional vector $\mathit{\gamma }\in {\mathbb{R}}^p$ faces the curse of dimensionality. However, FV is a statistic taking a value in $\mathbb{R}$; empirically, FV and LFV approximate the generalization gap $\Delta (\alpha ;\mathit{X})$ with a reasonable number of samples $T\in \mathbb{N}$, even if the dimension p is relatively large. See numerical experiments in Section 4.1.

Lastly, let us mention the computational time of FV and LFV. For each iteration, the number of operations to obtain the gradient in LFV is $O(p^2)$, while that to compute the covariance matrix in FV is $O(p^3)$: therefore, the computational time is expected to be reduced, as also demonstrated for linear models in Remark 5.

3.3 Application to Nonlinear Models

This section provides a procedure in applying LFV (13) to a nonlinear overparameterized model $g_{\mathit{\theta }}({\mathit{z}}_i)$ equipped with a parameter vector $\mathit{\theta }$ and an input vector ${\mathit{z}}_i$. The procedure is as follows: first, by replacing ${\mathcal{l}}_{\alpha }(\mathit{\beta })$ in the Langevin process (10) with $${\rho }_{\alpha }(\mathit{\theta }):=\frac{1}{n}\sum _{i=1}^n{(y_i-g_{\mathit{\theta }}({\mathit{z}}_i))}^2+\alpha \Vert \mathit{\theta }{\Vert }_2^2,$$the update (10) yields a sequence ${\left\{{\mathit{\gamma }}^{(t)}\right\}}_{t=1,2,\dots }$. Next, by substituting ${\mu }_i^{(t)}:=g_{{\mathit{\gamma }}^{(t)}}({\mathit{z}}_i)$ in (13), we obtain LFV for the nonlinear model $g_{\mathit{\theta }}$. We here emphasize that, while the exact form of the full quasi-posterior for nonlinear models is difficult to obtain and thus computing FV for nonlinear models is almost prohibited, LFV can be applied to even such nonlinear models.

The gradient of the overparameterized model $g_{\mathit{\theta }}$ is compatible with the gradient of its linear approximation. This application works if the linear approximation of nonlinear neural network is reasonable around the initial estimate. Again, we remark that this assumption is important in the neural tangent kernel literature (Jacot, Gabriel, and Hongler 2018; Arora et al. 2019).

3.4 Comparison to the Existing Methods

Here, we compare FV and LFV to existing methods for the estimation of the generalization gap using TIC (Takeuchi 1976). TIC estimates the generalization gap of regular statistical models by using a quantity (14) $$\text{tr}\left\{\widehat{\mathit{F}}{\widehat{\mathit{G}}}^{-1}\right\}$$(14) equipped with p × p matrices $$\widehat{\mathit{F}}:=\sum _{i=1}^n\frac{\partial \hspace{0.17em}\log \hspace{0.17em}f(y_i\mid {\mathit{x}}_i,\mathit{\beta })}{\partial \mathit{\beta }}\frac{\partial \hspace{0.17em}\log \hspace{0.17em}f(y_i\mid {\mathit{x}}_i,\mathit{\beta })}{\partial {\mathit{\beta }}^{\top }}{|}_{\mathit{\beta }=\widehat{\mathit{\beta }}}\text{and}$$ $$\widehat{\mathit{G}}:=\sum _{i=1}^n\frac{{\partial }^2\hspace{0.17em}\log \hspace{0.17em}f(y_i\mid {\mathit{x}}_i,\mathit{\beta })}{\partial \mathit{\beta }\partial {\mathit{\beta }}^{\top }}{|}_{\mathit{\beta }=\widehat{\mathit{\beta }}},$$ where $\widehat{\mathit{\beta }}$ is a maximum likelihood estimate. However, TIC cannot be applied to singular statistical models whose Hessian matrix $\widehat{\mathit{G}}$ degenerates (i.e., the inverse of $\widehat{\mathit{G}}$ does not exist), such as overparameterized models. To overcome this limitation of TIC, Thomas et al. (2020) replace the inverse matrix ${\widehat{\mathit{G}}}^{-1}$ in (14) with the κ-generalized inverse matrix ${\widehat{\mathit{G}}}_{\kappa }^{+}$; ${\widehat{\mathit{G}}}_{\kappa }^{+}$ is defined as $\mathit{V}{\Sigma }_{\kappa }^{+}{\mathit{V}}^{\top }$, where $\mathit{V}$ is a matrix of the eigenvectors of $\widehat{\mathit{G}}$, and ${\Sigma }_{\kappa }^{+}$ is the diagonal matrix of which the jth diagonal component is $${\lambda }_j({\widehat{\mathbf{G}}}_{\kappa }^{+}):=\{\begin{array}{cc}1/{\lambda }_j(\widehat{\mathbf{G}})\hfill & ({\lambda }_j(\widehat{\mathbf{G}})>\kappa )\hfill \\ 0\hfill & ({\lambda }_j(\widehat{\mathbf{G}})\le \kappa )\hfill \end{array},\hspace{1em}(j=1,2,\dots ,p)$$ with $\{{\lambda }_j(\widehat{\mathit{G}}):j=1,\dots ,p\}$ singular values of $\widehat{\mathit{G}}$. This TIC modification of Thomas et al. (2020) is numerically examined in Section 4.1, for overparameterized linear regression; it rather estimates the number of nonzero eigenvalues in $\widehat{\mathit{G}}=n^{-1}{\mathit{X}}^{\top }\mathit{X}$, that is, the number of nonzero singular values of $\mathit{X}$, but not the generalization gap $\Delta (\alpha ;\mathit{X})$.

Instead of modifying TIC, we can employ RIC proposed by Shibata (1976). RIC replaces the Hessian matrix $\widehat{\mathit{G}}$ in TIC with the regularized inverse matrix ${(\widehat{\mathit{G}}+\alpha \mathit{I})}^{-1}$. Moody (1992) and Murata, Yoshizawa, and Amari (1994) further generalized RIC to arbitrary loss functions and demonstrated its application in shallow neural network models with a small number of hidden units. RIC has much in common with FV; specifically, RIC is also an asymptotically unbiased estimator of the Gibbs generalization gap $\Delta (\alpha ;\mathit{X})$ for overparameterized linear regression, and numerically, RIC behaves similarly to FV (and LFV) for the experiments discussed in Section 4.1. However, RIC still requires the inverse matrix of the p × p matrix $\widehat{\mathit{G}}+\alpha \mathit{I}$, and its computational cost is intensive in the overparameterized setting ($p\ge n$).

4 Numerical Experiments

This section presents the numerical evaluation of the Langevin FV (13) using both linear and nonlinear models.

4.1 LFV and Baselines in Linear Models

We first evaluate LFV (13) and compare it with some baselines in linear overparameterized models. The set-up of the numerical experiments is summarized as follows:

• Synthetic data generation: For $p=1.5n$, orthogonal matrices $\mathit{U}\in {\mathbb{R}}^{n\times n},\mathit{V}\in {\mathbb{R}}^{p\times n}$ are iid from the uniform distribution over the set of orthogonal matrices satisfying ${\mathit{U}}^{\top }\mathit{U}={\mathit{V}}^{\top }\mathit{V}={\mathit{I}}_n$, respectively. Entries of ${\mathit{\beta }}_0\in {\mathbb{R}}^p$ are iid from $N(0,1/p)$. With given singular values ${\left\{s_i\right\}}_{i=1}^n$, the design matrix $\mathit{X}:=\mathit{US}{\mathit{V}}^{\top }$ is computed. $\mathit{y}$ is generated 50 times from $N_n(\mathit{X}{\mathit{\beta }}_0,{\mathit{I}}_n)$, that is, ${\sigma }_0^2=1$.

• Singular values: We consider three different types of singular values: (i) $\mathit{X}$ has a fixed intrinsic dimension $d_{*}=10$, that is, $s_1=\cdots =s_{10}=n^{1/2}$ and $s_{11}=s_{12}=\cdots =s_n=0$, (ii) $s_i=n^{1/2}{i}^{-1}$, and (iii) $s_i=n^{1/2}{i}^{-1/2}$. Theorem 3 proves the convergence of FV in settings (i) and (ii) since $\underset{n\to \infty }{\lim }{n}^{-1}{\sum }_{i=1}^n{s}_i^2<\infty$ in (i) and (ii), whereas it does not prove the convergence in setting (iii) since $\underset{n\to \infty }{\lim }{n}^{-1}{\sum }_{i=1}^n{s}_i^2$ is not bounded in (iii).

• Evaluation: We evaluate the following baselines and LFV with 50 experiments. Throughout these experiments, we employ $\alpha =0.1$ for ridge regularization (2).

  1. TIC penalty (Takeuchi 1976) is extended to the overparameterization setting (Thomas et al. 2020) by replacing the inverse of the Hessian matrix $\widehat{\mathit{G}}$ in TIC with the κ-generalized inverse matrix ${\widehat{\mathit{G}}}_{\kappa }^{+}$. See Section 3.4 for the definition of $\text{TIC}(\kappa ):=\text{tr}\left\{\widehat{\mathit{F}}{\widehat{\mathit{G}}}_{\kappa }^{+}\right\}$.

  2. FV is empirically computed with $T=15n$ samples of $\mathit{\beta }$ generated from the quasi-posterior (1). Its expectation shown in Lemma 2 is also computed.

  3. LFV is computed with $\delta =1/(10n)$ and $T=15n$ Langevin samples $\mathit{\gamma }$ (see (10)).

The results are presented in the following Tables 1–3. Overall, FV and LFV were able to more effectively estimate the generalization gap $\Delta (\alpha ;\mathit{X})$ ($p=1.5n,n\to \infty$), as compared to TIC, and their biases were not drastically different. TIC was entirely inaccurate in settings (ii) and (iii); together with the result of (i), TIC was able to estimate the number of nonzero eigenvalues in the Hessian matrix $\widehat{\mathit{G}}=n^{-1}{\mathit{X}}^{\top }\mathit{X}$ (i.e., the number of nonzero singular values of $\mathit{X}$), which is generally different from the generalization gap $\Delta (\alpha ;\mathit{X})$.

Table 1 The generalization gap $\Delta (\alpha ;\mathit{X})$ and its estimates (TIC, FV, LFV) with the standard deviations in setting (i) where $s_i=n^{1/2}$ for $i=1,2,\dots ,10$ and 0 otherwise.

n	100	200	300	400
$\Delta (\alpha ;\mathit{X})$	9.091	9.091	9.091	9.091
$\text{TIC} (\kappa =0)$	8.752 ± 1.284	9.272 ± 0.941	9.633 ± 0.829	9.536 ± 0.769
$\text{TIC} (\kappa =0.1)$	8.752 ± 1.284	9.272 ± 0.941	9.633 ± 0.829	9.536 ± 0.769
$\text{FV}(\alpha ;\mathit{X})$	8.432 ± 1.180	8.658 ± 0.864	8.907 ± 0.761	8.781 ± 0.693
${\mathbb{E}}_{\mathit{y}}[\text{FV}(\alpha ;\mathit{X})]$	8.527	8.810	8.902	8.946
$\text{LFV}(\alpha ;\mathit{X})$	8.317 ± 1.144	8.815 ± 0.931	9.044 ± 0.8463	8.951 ± 0.743

NOTE: The generalization gap and LFV estimate are gray-colored.

Table 2 The generalization gap $\Delta (\alpha ;\mathit{X})$ and its estimates (TIC, FV, LFV) with the standard deviations in setting (ii), where $s_i=n^{1/2}{i}^{-1}$.

n	200	400	600	800
$\Delta (\alpha ;\mathit{X})$	4.417	4.442	4.451	4.455
$\text{TIC} (\kappa =0)$	192.5 ± 18.01	390.1 ± 24.97	589.2 ± 35.37	787.6 ± 39.41
$\text{TIC} (\kappa =0.1)$	135.6 ± 12.83	194.1 ± 12.41	239.8 ± 14.52	277.7 ± 13.92
$\text{FV}(\alpha ;\mathit{X})$	4.295 ± 0.415	4.302 ± 0.354	4.351 ± 0.306	4.369 ± 0.255
${\mathbb{E}}_{\mathit{y}}[\text{FV}(\alpha ;\mathit{X})]$	4.305	4.385	4.412	4.425
$\text{LFV}(\alpha ;\mathit{X})$	3.985 ± 0.491	4.193 ± 0.467	4.272 ± 0.357	4.341 ± 0.305

NOTE: The generalization gap and LFV estimate are gray-colored.

Table 3 The generalization gap $\Delta (\alpha ;\mathit{X})$ and its estimates (TIC, FV, LFV) with the standard deviations in setting (iii), where $s_i=n^{1/2}{i}^{-1/2}$.

n	200	400	600	800
$\Delta (\alpha ;\mathit{X})$	29.98	36.66	40.63	43.46
$\text{TIC} (\kappa =0)$	149.7 ± 14.68	335.3 ± 21.68	527.6 ± 33.40	721.0 ± 37.64
$\text{TIC} (\kappa =0.1)$	149.7 ± 14.68	335.3 ± 21.68	527.6 ± 33.40	721.0 ± 37.64
$\text{FV}(\alpha ;\mathit{X})$	24.68 ± 2.250	32.32 ± 2.026	37.06 ± 2.292	40.33 ± 2.046
${\mathbb{E}}_{\mathit{y}}[\text{FV}(\alpha ;\mathit{X})]$	24.78	32.60	37.26	40.53
$\text{LFV}(\alpha ;\mathit{X})$	22.32 ± 2.092	30.77 ± 2.100	35.83 ± 2.346	39.42 ± 2.052

NOTE: The generalization gap and LFV estimate are gray-colored.

Remark 5

(Computational time). Let us mention computational time of FV and LFV. We check the computation time by the linear synthetic dataset with $n=100,p=2000,\alpha =0.01,s_i=n^{1/2}{i}^{-1/2},T=500$ five times. Among five experiments, the processing time for LFV is 5.33 ± 0.09 sec, and that for FV is 18.87 ± 0.18 sec, which implies that LFV is in fact faster than FV. For these experiments, we used AMD Ryzen 7 5700X processors. The computation is not parallelized for fair comparison.

4.2 LFV for Nonlinear Models

This section presents the evaluation of the Langevin FV for nonlinear neural networks via synthetic dataset experiments. We use the procedure in Section 3.3 in applying LFV to nonlinear models. The set-up is summarized as follows:

• Synthetic data: Let $n\in \{50,500,1000\},{\sigma }^2=1$. For $i=1,2,\dots ,n$, we generate ${\mathit{z}}_i\stackrel{\text{iid}}{\sim }{N}_d(0,{\mathit{I}}_d)$ and ${\epsilon }_i\stackrel{\text{iid}}{\sim }N(0,{\sigma }^2)$, and set ${\mu }_i:=\mu ({\mathit{z}}_i)$ and $y_i:={\mu }_i+{\epsilon }_i$ with a function $\mu (\mathit{z}):=3\tanh (〈\mathit{z},\mathit{1}〉/2)$.

• Neural network (NN): We employ a fully-connected one-hidden-layer NN: $$g_{\mathit{\theta }}(\mathit{z}):=〈{\mathit{\theta }}^{(2)},\tanh ({\mathit{\theta }}^{(1)}\mathit{z}+{\mathit{\theta }}^{(0)})〉$$with $M\in \{50,100,150\}$ hidden units, where $\mathit{\theta }=({\mathit{\theta }}^{(0)},{\mathit{\theta }}^{(1)},{\mathit{\theta }}^{(2)})\in {\mathbb{R}}^M\times {\mathbb{R}}^{M\times d}\times {\mathbb{R}}^M$ is a parameter vector (and thus the number of parameters is $p=M(d+2)$) and the function $\tanh (·)$ applies the hyperbolic tangent function $\tanh (z):=(\hspace{0.17em}\exp \hspace{0.17em}(z)-\hspace{0.17em}\exp \hspace{0.17em}(-z))/(\hspace{0.17em}\exp \hspace{0.17em}(z)+\hspace{0.17em}\exp \hspace{0.17em}(-z))$ entry-wise. Let ${\tilde{\mathit{\theta }}}_0$ be a parameter satisfying $g_{{\tilde{\mathit{\theta }}}_0}(\mathit{z})=\mu (\mathit{z})$. For each experiment, we employ the true parameter ${\mathit{\theta }}_0$ that is iid drawn from the element-wise independent normal distribution with mean ${\tilde{\mathit{\theta }}}_0$ and the element-wise variance 0.01. Further, we initialize the parameter $\mathit{\theta }$ of the NN to be trained by the element-wise independent normal distribution whose mean is ${\mathit{\theta }}_0$ with the larger element-wise variance 1, and update $\mathit{\theta }$ by a full-batch gradient descent with a learning rate 0.1, ridge-regularization coefficient $\alpha ={10}^{-3}$ and, $0.3n$ iterations.

• Langevin FV: For each setting $(d,M)\in \{5,10,15\}\times \{50,100,150\}$, we take the average of LFV over 25 experiments. In each experiment, we randomly generate ${\{({\mathit{z}}_i,y_i)\}}_{i=1}^n\subset {\mathbb{R}}^d\times \mathbb{R}$, train the NN $g_{\mathit{\theta }}$, and compute T = 3000 iterations of the Langevin process with $\delta ={10}^{-5}$. We discard the first $0.1T$ iterations, and employ the remaining $0.9T$ iterations to compute LFV.

• Generalization gap $\tilde{\Delta }$: For each setting $(d,M)\in \{5,10,15\}\times \{50,100,150\}$, we take average of the following generalization gap $$\tilde{\Delta }:=\frac{1}{2{\sigma }_0^2}\{{\mathbb{E}}_{{\mathit{y}}^{*}}\left(\frac{1}{n}\sum _{i=1}^n{\left\{y_i^{*}-{\mathit{g}}_{\widehat{\mathit{\theta }}}({\mathit{z}}_i)\right\}}^2\right)$$ $$\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}-\frac{1}{n}\sum _{i=1}^n{\left\{y_i-{\mathit{g}}_{\widehat{\mathit{\theta }}}({\mathit{z}}_i)\right\}}^2\}$$over 200 times experiments. In each experiment, we randomly generate ${\{({\mathit{z}}_i,y_i)\}}_{i=1}^n\subset {\mathbb{R}}^d\times \mathbb{R}$, train the NN $g_{\mathit{\theta }}$, and evaluate $\tilde{\Delta }$ by leveraging the ground-truth ${\mu }_i=\mu ({\mathit{z}}_i)$ and ${\sigma }_0^2=1$.

LFV and generalization loss $\tilde{\Delta }$ for a nonlinear NN with M hidden units are shown in Tables 4–6. LFV values for overparameterized models (i.e., $p=M(d+2)>n$) are in gray. It was observed that the value of LFV provides a rough estimate of the generalization gap even for nonlinear overparameterized NN models, though our overparameterized theories on FV and LFV are proved for linear models. Table 4 (n = 50) and Table 6 (n = 1000) suggested that the biases of LFV reduces as the sample size grows.

Table 4 The generalization gap and LFV for the NN model with n = 50 and T = 3000.

	M = 50		M = 100		M = 150
	LFV	$\tilde{\Delta }$	LFV	$\tilde{\Delta }$	LFV	$\tilde{\Delta }$
d = 5	2.90 ± 0.84	5.78	3.50 ± 1.13	5.84	3.85 ± 0.90	5.88
d = 10	3.98 ± 1.22	11.66	4.25 ± 1.16	11.80	5.79 ± 2.03	11.95
d = 15	4.73 ± 1.70	15.17	6.23 ± 1.65	15.55	7.28 ± 1.80	15.94

NOTE: LFV values for the overparameterized regime (i.e., $p=M(d+2)>n$) are gray-colored.

Table 5 The generalization gap and LFV for the NN model with n = 500 and T = 3000.

	M = 50		M = 100		M = 150
	LFV	$\tilde{\Delta }$	LFV	$\tilde{\Delta }$	LFV	$\tilde{\Delta }$
d = 5	7.72 ± 1.31	8.85	8.42 ± 1.18	8.00	9.16 ± 1.51	8.48
d = 10	14.64 ± 1.99	17.79	15.79 ± 1.69	18.97	17.58 ± 2.11	19.56
d = 15	20.34 ± 1.79	29.85	24.68 ± 2.07	30.86	26.32 ± 2.82	31.24

NOTE: LFV values for the overparameterized regime (i.e., $p=M(d+2)>n$) are gray-colored.

Table 6 The generalization gap and LFV for the NN model with n = 1000 and T = 3000.

	M = 50		M = 100		M = 150
	LFV	$\tilde{\Delta }$	LFV	$\tilde{\Delta }$	LFV	$\tilde{\Delta }$
d = 5	8.86 ± 1.20	9.13	9.99 ± 1.32	9.80	10.78 ± 1.38	10.15
d = 10	17.16 ± 1.56	23.46	19.09 ± 1.55	23.87	21.50 ± 2.41	23.99
d = 15	25.42 ± 2.04	31.70	30.25 ± 2.21	31.82	32.81 ± 2.30	31.87

NOTE: LFV values for the overparameterized regime (i.e., $p=M(d+2)>n$) are gray-colored.

4.3 LFV for Nonlinear Models Using Real Datasets

In this section, we compare LFV to cross-validation statistic using real datasets with small sample sizes, where we note that the cross-validation becomes prohibited as the sample size becomes larger. We use the procedure in Section 3.3 in applying LFV to nonlinear models. The set-up is summarized as follows.

• Real data: We collected the following seven regression datasets from KEEL dataset repository (Alcalá-Fdez et al. 2011): “machineCPU” ($n=208,d=6$), “wankara” ($n=320,d=9$), “baseball” ($n=336,d=16$), “dee” ($n=364,d=6$), “autoMPG6” ($n=391,d=5$), “autoMPG8” ($n=391,d=7$), and “stock” ($n=949,d=9$). These datasets are selected so that our NN structure described below becomes overparameterized. For each dataset, we standardized (scaling and centering) the design matrix $\mathit{X}$ and the target variable $\mathit{y}$.

• Neural network (NN): We employ the same architecture as the nonlinear NN in Section 4.2, with the number of hidden layers M = 100. For the neural netweork training, we first initialize the parameter $\mathit{\theta }$ randomly by the normal distribution $N_p(\mathit{0},d^{-1/2}{\mathit{I}}_p)$; we employ 20 different random initializations for each dataset experiment. We train the NN by the fullbatch gradient descent (using the entire dataset) with the learning rate 0.01, ridge-regularization coefficient $\alpha =0.1$. Gradient descent algorithm is terminated if $|{\text{loss}}_t-{\text{loss}}_{t-1}|/\left|{\text{loss}}_{t-1}\right|<{10}^{-6}$, where ${\text{loss}}_t$ denotes the training loss at the iteration t.

• Cross-validation (CV): We divide each dataset into the training set (90%) and test set (10%) uniformly randomly. Starting from the NN parameters already trained with the entire dataset as shown above, we train NN by gradient descent with the learning rate 0.001 using the divided training set. Gradient descent algorithm for the jth CV instance is terminated if $|{\text{loss}}_{j,t}-{\text{loss}}_{j,t-1}|/\left|{\text{loss}}_{j,t-1}\right|<{10}^{-6}$, where ${\text{loss}}_{j,t}$ denotes the training loss at the iteration t, using the training set in the jth CV instance. After the training, we compute 10-fold cross-validation statistic: $$\text{CV}:=\frac{1}{2{\widehat{\sigma }}_0^2}\frac{1}{10}\sum _{j=1}^{10}\{\mathrm{}\frac{1}{n_{\text{test}(j)}}\sum _{i=1}^{n_{\text{test}(j)}}\{y_{\text{test}(j),i}-g_{{\widehat{\mathit{\theta }}}_{\text{train}(j)}}({\mathit{x}}_{\text{test}(j),i}){\}}^2,\mathrm{}\}$$where ${\{({\mathit{x}}_{\text{test}(j),i},y_{\text{test}(j),i})\}}_{i=1}^{n_{\text{test}}}$ denotes the jth test set and ${\{({\mathit{x}}_{\text{train}(j),i},y_{\text{train}(j),i})\}}_{i=1}^{n_{\text{train}}}$ denotes the jth training set. $g_{{\widehat{\mathit{\theta }}}_{\text{train}(j)}}$ is the NN trained with the jth training set. Since the training of the NN was unstable in some of the training sets, we ignored a few instances of 10 training sets in CV.

• Langevin FV: After training the NN $g_{\mathit{\theta }}$, we compute $T=10,000$ iterations of the Langevin process with $\delta =5\times {10}^{-6}/{\widehat{\sigma }}_0^2$, where ${\widehat{\sigma }}_0^2:={(n-1)}^{-1}{\sum }_{i=1}^n{\{y_i-g_{\widehat{\mathit{\theta }}}({\mathit{x}}_i)\}}^2$ denotes the estimator of ${\sigma }_0^2$ for each dataset. We discard first 2500 iterations and compute LFV using the remaining 7500 iterations; we further compute $$\underset{(\text{trainingloss})}{\underbrace{\frac{1}{2{\widehat{\sigma }}_0^2}\sum _{i=1}^n{\{y_i-g_{{\widehat{\mathit{\theta }}}_{\text{train}}}({\mathit{x}}_i)\}}^2}}+\text{LFV}(\alpha ;\mathit{X}),$$where it can be regarded as a Langevin variant of WAIC (Watanabe 2010), and compare it to the CV statistic approximating the generalization loss.

• Experimental Setting: we compute $\text{CV}$ statistic and training loss + LFV for each dataset. More precisely, we compute the average and standard deviation of these values over different 20 random initializations.

Results are shown in Table 7. Overall, training loss + LFV roughly approximated the CV statistic; they are expected to be compatible as both of these values approximate the generalization loss. Further, we computed the Spearman’s rank correlation between the averaged training loss + LFV and the averaged CV statistic (over 20 random initializations), for seven datasets. The correlation was 0.75, while that between the training loss (without LFV) and CV statistic was –0.75.

Table 7 CV statistic and training loss + LFV over 20 different initializations. n denotes the number of samples, and $p=M(d+2)$ denotes the number of parameters used in NN.

Dataset	n	p	$\text{CV}$	Training loss + LFV
machineCPU	208	800	0.637 ± 0.034	0.646 ± 0.123
wankara	320	1100	0.417 ± 0.019	0.560 ± 0.068
baseball	336	1800	0.671 ± 0.011	0.661 ± 0.103
dee	364	800	0.484 ± 0.025	0.622 ± 0.107
autoMPG6	391	700	0.523 ± 0.045	0.610 ± 0.101
autoMPG8	391	900	0.431 ± 0.055	0.559 ± 0.056
stock	949	1100	0.526 ± 0.016	0.580 ± 0.122

In addition to the computational burden, we note that CV statistic became unstable even in these experiments, as the multiple NNs trained with different CV instances fell in different local optima. So we ignored several CV instances for more stable computation; LFV was more stably computed as LFV was computed with only a single training of the NN.

Therefore, we expect that the proposed LFV can be used as a substitute of the cross-validation statistic for real datasets, even with larger sample sizes.

5 Conclusion

In this article, we considered a Gibbs generalization gap estimation method for overparameterized models. We proved that FV works as an asymptotically unbiased estimator, even in an overparameterization setting. We proposed a Langevin approximation of FV for efficient computation and applied it to overparameterized linear regression and nonlinear neural network models.

Supplementary Materials

Supplementary material contains the proofs of theorems, and descriptions of source codes to reproduce experimental results.

Acknowledgments

We thank the editor, the AE, and two anonymous reviewers for constructive comments and suggestions. We also thank Eiki Shimizu for suggesting several references, and Tetsuya Takabatake and Yukito Iba for helpful discussions.

Additional information

Funding

A. Okuno was supported by JST CREST (JPMJCR21N3) and JSPS KAKENHI (21K17718, 22H05106). K. Yano was supported by JST CREST (JPMJCR1763), JSPS KAKENHI (19K20222, 21H05205, 21K12067), and MEXT (JPJ010217).