Year: 2020

Abstract:
In singular models, the optimal set of parameters forms an analytic set with singularities and classical statistical inference cannot be applied to such models. This is significant for deep learning as neural networks are singular and thus “dividing” by the determinant of the Hessian or employing the Laplace approximation are not appropriate. Despite its potential for addressing fundamental issues in deep learning, singular learning theory appears to have made little inroads into the developing canon of deep learning theory. Via a mix of theory and experiment, we present an invitation to singular learning theory as a vehicle for understanding deep learning and suggest important future work to make singular learning theory directly applicable to how deep learning is performed in practice.

text: In singular models, the optimal set of parameters forms an analytic set with singularities and classical statistical inference cannot be applied to such models. This is significant for deep learning as neural networks are singular and thus “dividing” by the determinant of the Hessian or employing the Laplace approximation are not appropriate. Despite its potential for addressing fundamental issues in deep learning, singular learning theory appears to have made little inroads into the developing canon of deep learning theory. Via a mix of theory and experiment, we present an invitation to singular learning theory as a vehicle for understanding deep learning and suggest important future work to make singular learning theory directly applicable to how deep learning is performed in practice.

heading: 1 INTRODUCTION
text: It has been understood for close to twenty years that neural networks are singular statistical models (Amari et al., 2003; Watanabe, 2007). This means, in particular, that the set of network weights equivalent to the true model under the Kullback-Leibler divergence forms a real analytic variety which fails to be an analytic manifold due to the presence of singularities. It has been shown by Sumio Watanabe that the geometry of these singularities controls quantities of interest in statistical learning theory, e.g., the generalisation error. Singular learning theory (Watanabe, 2009) is the study of singular models and requires very different tools from the study of regular statistical models. The breadth of knowledge demanded by singular learning theory – Bayesian statistics, empirical processes and algebraic geometry – is rewarded with profound and surprising results which reveal that singular models are different from regular models in practically important ways. To illustrate the relevance of singular learning theory to deep learning, each section of this paper illustrates a key takeaway idea1.
The real log canonical threshold (RLCT) is the correct way to count the effective number of parameters in a deep neural network (DNN) (Section 4). To every (model, truth, prior) triplet is associated a birational invariant known as the real log canonical threshold. The RLCT can be understood in simple cases as half the number of normal directions to the set of true parameters. We will explain why this matters more than the curvature of those directions (as measured for example by eigenvalues of the Hessian) laying bare some of the confusion over “flat” minima.
For singular models, the Bayes predictive distribution is superior to MAP and MLE (Section 5). In regular statistical models, the 1) Bayes predictive distribution, 2) maximum a posteriori (MAP) estimator, and 3) maximum likelihood estimator (MLE) have asymptotically equivalent generalisation error (as measured by the Kullback-Leibler divergence). This is not so in singular models. We illustrate in our experiments that even “being Bayesian” in just the final layers improves generalisation over MAP. Our experiments further confirm that the Laplace approximation of the predictive distribution Smith & Le (2017); Zhang et al. (2018) is not only theoretically inappropriate but performs poorly.
Simpler true distribution means lower RLCT (Section 6). In singular models the RLCT depends on the (model, truth, prior) triplet whereas in regular models it depends only on the (model, prior) pair. The RLCT increases as the complexity of the true distribution relative to the supposed model increases. We verify this experimentally with a simple family of ReLU and SiLU networks.
1The code to reproduce all experiments in the paper will be released on Github. For now, see the zip file.

heading: 2 RELATED WORK
text: In classical learning theory, generalisation is explained by measures of capacity such as the l2 norm, Radamacher complexity, and VC dimension (Bousquet et al., 2003). It has become clear however that these measures cannot capture the empirical success of DNNs (Zhang et al., 2017). For instance, over-parameterised neural networks can easily fit random labels (Zhang et al., 2017; Du et al., 2018; Allen-Zhu et al., 2019b) indicating that complexity measures such as Rademacher complexity are very large. There is also a slate of work on generalisation bounds in deep learning. Uniform convergence bounds (Neyshabur et al., 2015; Bartlett et al., 2017; Neyshabur & Li, 2019; Arora et al., 2018) usually cannot provide non-vacuous bounds. Data-dependent bounds (Brutzkus et al., 2018; Li & Liang, 2018; Allen-Zhu et al., 2019a) consider the “classifiability” of the data distribution in generalisation analysis of neural networks. Algorithm-dependent bounds (Daniely, 2017; Arora et al., 2019; Yehudai & Shamir, 2019; Cao & Gu, 2019) consider the relation of Gaussian initialisation and the training dynamics of (stochastic) gradient descent to kernel methods (Jacot et al., 2018).
In contrast to many of the aforementioned works, we are interested in estimating the conditional distribution q(y|x). Specifically, we measure the generalisation error of some estimate q̂n(y|x) in terms of the Kullback-Leibler divergence between q and q̂n, see (8). The next section gives a crash course on singular learning theory. The rest of the paper illustrates the key ideas listed in the introduction. Since we cover much ground in this short note, we will review other relevant work along the way, in particular literature on “flatness”, the Laplace approximation in deep learning, etc.

heading: 3 SINGULAR LEARNING THEORY
text: To understand why classical measures of capacity fail to say anything meaningful about DNNs, it is important to distinguish between two different types of statistical models. Recall we are interested in estimating the true (and unknown) conditional distribution q(y|x) with a class of models {p(y|x,w) : w ∈ W} where W ⊂ Rd is the parameter space. We say the model is identifiable if the mapping w 7→ p(y|x,w) is one-to-one. Let q(x) be the distribution of x. The Fisher information matrix associated with the model {p(y|x,w) : w ∈W} is the matrix-valued function on W defined by
I(w)ij =
∫ ∫ ∂
∂wi [log p(y|x,w)] ∂ ∂wj [log p(y|x,w)]q(y|x)q(x)dxdy,
if this integral is finite. Following the conventions in Watanabe (2009), we have the following bifurcation of statistical models. A statistical model p(y|x,w) is called regular if it is 1) identifiable and 2) has positive-definite Fisher information matrix. A statistical model is called strictly singular if it is not regular.
Let ϕ(w) be a prior on the model parameters w. To every (model, truth, prior) triplet, we can associate the zeta function, ζ(z) = ∫ K(w)zϕ(w) dw, z ∈ C, where K(w) is the Kullback-Leibler (KL) divergence between the model p(y|x,w) and the true distribution q(y|x):
K(w) := ∫ ∫ q(y|x) log q(y|x)
p(y|x,w) q(x) dx dy. (1)
For a (model, truth, prior) triplet (p(y|x,w), q(y|x), ϕ), let −λ be the maximum pole of the corresponding zeta function. We call λ the real log canonical threshold (RLCT) (Watanabe, 2009) of the (model, truth, prior) triplet. The RLCT is the central quantity of singular learning theory.
By Watanabe (2009, Theorem 6.4) the RLCT is equal to d/2 in regular statistical models and bounded above by d/2 in strictly singular models if realisability holds: let
W0 = {w ∈W : p(y|x,w) = q(y|x)} be the set of true parameters, we say q(y|x) is realisable by the model class if W0 is non-empty. The condition of realisability is critical to standard results in singular learning theory. Modifications to the theory are needed in the case that q(y|x) is not realisable, see the condition called relatively finite variance in Watanabe (2018).
Neural networks in singular learning theory. Let W ⊆ Rd be the space of weights of a neural network of some fixed architecture, and let f(x,w) : RN ×W −→ RM be the associated function.
We shall focus on the regression task and study the model
p(y|x,w) = 1 (2π)M/2 exp ( − 12‖y − f(x,w)‖ 2 )
(2)
but singular learning theory can also apply to classification, for instance. It is routine to check (see Appendix A.1) that for feedforward ReLU networks not only is the model strictly singular but the matrix I(w) is degenerate for all nontrivial weight vectors and the Hessian of K(w) is degenerate at every point of W0.
RLCT plays an important role in model selection. One of the most accessible results in singular learning theory is the work related to the widely-applicable Bayesian information criterion (WBIC) Watanabe (2013), which we briefly review here for completeness. Let Dn = {(xi, yi)}ni=1 be a dataset of input-output pairs. Let Ln(w) be the negative log likelihood
Ln(w) = − 1
n n∑ i=1 log p(yi|xi, w) (3)
and p(Dn|w) = exp(−nLn(w)). The marginal likelihood of a model {p(y|x,w) : w ∈ W} is given by p(Dn) = ∫ W p(Dn|w)ϕ(w) dw and can be loosely interpreted as the evidence for the model. Between two models, we should prefer the one with higher model evidence. However, since the marginal likelihood is an intractable integral over the parameter space of the model, one needs to consider some approximation.
The well-known Bayesian Information Criterion (BIC) derives from an asymptotic approximation of − log p(Dn) using the Laplace approximation, leading to BIC = nLn(wMLE) + d2 log n. Since we want the marginal likelihood of the data for some given model to be high one should almost never adopt a DNN according to the BIC, since in such models d may be very large. However, this argument contains a serious mathematical error: the Laplace approximation used to derive BIC only applies to regular statistical models, and DNNs are not regular. The correct criterion for both regular and strictly singular models was shown in Watanabe (2013) to be nLn(w0)+λ log nwherew0 ∈W0 and λ is the RLCT. Since DNNs are highly singular λ may be much smaller than d/2 (Section 6) it is possible for DNNs to have high marginal likelihood – consistent with their empirical success.

heading: 4 VOLUME DIMENSION, EFFECTIVE DEGREES OF FREEDOM, AND FLATNESS
text: Volume codimension. The easiest way to understand the RLCT is as a volume codimension (Watanabe, 2009, Theorem 7.1). Suppose that W ⊆ Rd and W0 is nonempty, i.e., the true distribution is realisable. We consider a special case in which the KL divergence in a neighborhood of every point v0 ∈W0 has an expression in local coordinates of the form
K(w) = d′∑ i=1 ciw 2 i , (4)
where the coefficients c1, . . . , cd′ > 0 may depend on v0 and d′ may be strictly less than d. If the model is regular then this is true with d = d′ and if it holds for d′ < d then we say that the pair (p(y|x,w), q(y|x)) is minimally singular. It follows that the set W0 ⊆ W of true parameters is a regular submanifold of codimension d′ (that is, W0 is a manifold of dimension d− d′ where W has dimension d). Under this hypothesis there are, near each true parameter v0 ∈ W0, exactly d − d′ directions in which v0 can be varied without changing the model p(y|x,w) and d′ directions in which varying the parameters does change the model. In this sense, there are d′ effective parameters near v0.
This number of effective parameters can be computed by an integral. Consider the volume of the set of almost true parameters V (t, v0) = ∫ K(w)<t
ϕ(w)dw where the integral is restricted to a small closed ball around v0. As long as the prior ϕ(w) is non-zero on W0 it does not affect the relevant features of the volume, so we may assume ϕ is constant on the region of integration in the first d′ directions and normal in the remaining directions, so up to a constant depending only on d′ we have
V (t, v0) ∝ td ′/2
√ c1 · · · cd′
(5)
and we can extract the exponent of t in this volume in the limit
d′ = 2 lim t→0
log { V (at, v0)/V (t, v0) } log(a)
(6)
for any a > 0, a 6= 1. We refer to the right hand side of (6) as the volume codimension at v0. The function K(w) has the special form (4) locally with d′ = d if the statistical model is regular (and realisable) and with d′ < d in some singular models such as reduced rank regression (Appendix A.2). While such a local form does not exist for a singular model generally (in particular for neural networks) nonetheless under natural conditions (Watanabe, 2009, Theorem 7.1) we have V (t, v0) = ctλ + o(tλ) where c is a constant. We assume that in a sufficiently small neighborhood of v0 the point RLCT λ at v0 (Watanabe, 2009, Definition 2.7) is less than or equal to the RLCT at every point in the neighborhood so that the multiplicity m = 1, see Section 7.6 of (Watanabe, 2009) for relevant discussion. It follows that the limit on the right hand side of (6) exists and is equal to λ. In particular λ = d′/2 in the minimally singular case.
Note that for strictly singular models such as DNNs 2λ may not be an integer. This may be disconcerting but the connection between the RLCT, generalisation error and volume dimension strongly suggests that 2λ is nonetheless the only geometrically meaningful “count” of the effective number of parameters near v0.
RLCT and likelihood vs temperature. Again working with the model in (2), consider the expectation over the posterior at temperature T as defined in (17) of the negative log likelihood (3)
E(T ) = E1/Tw [ nLn(w) ] = E1/Tw [ 1 2 n∑ i=1 ‖yi − f(xi, w)‖2 ] + nM 2 log(2π) .
Note that when n is large Ln(v0) ≈ M2 log(2π) for any v0 ∈ W0 so for T ≈ 0 the posterior concentrates around the set W0 of true parameters and E(T ) ≈ nM2 log(2π). Consider the increase ∆E = E(T + ∆T )− E(T ) corresponding to an increase in temperature ∆T . It can be shown that ∆E ≈ λ∆T where the reader should see (Watanabe, 2013, Corollary 3) for a precise statement. As the temperature increases, samples taken from the tempered posterior are more distant from W0 and the error E will increase. If λ is smaller then for a given increase in temperature the quantity E increases less: this is one way to understand intuitively why a model with smaller RLCT generalises better from the dataset Dn to the true distribution.
Flatness. It is folklore in the deep learning community that flatness of minima is related to generalisation (Hinton & Van Camp, 1993; Hochreiter & Schmidhuber, 1997) and this claim has been revisited in recent years (Chaudhari et al., 2017; Smith & Le, 2017; Jastrzebski et al., 2017; Zhang et al., 2018). In regular models this can be justified using the lower order terms of the asymptotic expansion of the Bayes free energy (Balasubramanian, 1997, §3.1) but the argument breaks down in strictly singular models, since for example the Laplace approximation of Zhang et al. (2018) is invalid. The point can be understood via an analysis of the version of the idea in (Hochreiter & Schmidhuber, 1997). Their measure of entropy compares the volume of the set of parameters with tolerable error t0 (our almost true parameters) to a standard volume
− log [V (t0, v0)
t d/2 0
] = d− d′
2 log(t0) +
1 2 d∑ i=1 log ci . (7)
Hence in the case d = d′ the quantity − 12 ∑ i log(ci) is a measure of the entropy of the set of true parameters near w0, a point made for example in Zhang et al. (2018). However when d′ < d this conception of entropy is inappropriate because of the d − d′ directions in which K(w) is flat near v0, which introduce the t0 dependence in (7).

heading: 5 GENERALISATION
text: The generalisation puzzle (Poggio et al., 2018) is one of the central mysteries of deep learning. Theoretical investigations into the matter is an active area of research Neyshabur et al. (2017). Many of the recent proposals of capacity measures for neural networks are based on the eigenspectrum of the (degenerate) Hessian, e.g., Thomas et al. (2019); Maddox et al. (2020). But this is not appropriate for singular models, and hence for DNNs.
Since we are interested in learning the distribution, our notion of generalisation is slightly different, being measured by the KL divergence. Precise statements regarding the generalisation behavior in singular models can be made using singular learning theory. Let the network weights be denoted θ rather than w for reasons that will become clear. Recall in the Bayesian paradigm, prediction proceeds via the so-called Bayes predictive distribution, p(y|x,Dn) = ∫ p(y|x, θ)p(θ|Dn) dθ.More commonly encountered in deep learning practice are the MAP and MLE point estimators. While in a regular statistical model, the three estimators 1) Bayes predictive distribution, 2) MAP, and 3) MLE have the same leading term in their asymptotic generalisation behavior, the same is not true in singular models. More precisely, let q̂n(y|x) be some estimate of the true unknown conditional density q(y|x) based on the dataset Dn. The generalisation error of the predictor q̂n(y|x) is
G(n) := KL(q(y|x)||q̂n(y|x)) = ∫ ∫
q(y|x) log q(y|x) q̂n(y|x) q(x) dy dx. (8)
To account for sampling variability, we will work with the average generalisation error, EnG(n), where En denotes expectation over the dataset Dn. By Watanabe (2009, Theorem 1.2 and Theorem 7.2), we have
EnG(n) = λ/n+ o(1/n) if q̂n is the Bayes predictive distribution, (9) where λ is the RLCT corresponding to the triplet (p(y|x, θ), q(y|x), ϕ(θ)). In contrast, we should note that Zhang et al. (2018) and Smith & Le (2017) rely on the Laplace approximation to explain the generalisation of the Bayes predictive distribution though both works acknowledge the Laplace approximate is inappropriate. For completeness, a quick sketch of the derivation of (9) is provided in Appendix A.4. Now by (Watanabe, 2009, Theorem 6.4) we have
EnG(n) = C/n+ o(1/n) if q̂n is the MAP or MLE, (10) whereC (different for MAP and MLE) is the maximum of some Gaussian process. For regular models, the MAP, MLE, and the Bayes predictive distribution have the same leading term for EnG(n) since λ = C = d/2. However in singular models, C is generally greater than λ, meaning we should prefer the Bayes predictive distribution for singular models.
That the RLCT has such a simple relationship to the Bayesian generalisation error is remarkable. On the other hand, the practical implications of (19) are limited since the Bayes predictive distribution is intractable. While approximations to the Bayesian predictive distribution, say via variational inference, might inherit a similar relationship between generalisation and the (variational) RLCT, serious theoretical developments will be required to rigorously establish this. The challenge comes from the fact that for approximate Bayesian predictive distributions, the free energy and generalisation error may have different learning coefficients λ. This was well documented in the case of a neural network with one hidden layer (Nakajima & Watanabe, 2007).
We set out to investigate whether certain very simple approximations of the Bayes predictive distribution can already demonstrate superiority over point estimators. Suppose the input-target relationship is modeled as in (2) but we write θ instead of w. We set q(x) = N(0, I3). For now consider the realisable case, q(y|x) = p(y|x, θ0) where θ0 is drawn randomly according to the default initialisation in PyTorch when model (2) is instantiated. We calculate EnG(n) using multiple datasets Dn and a large testing set, see Appendix A.5 for more details.
Since f is a hierarchical model, let’s write it as fθ(·) = h(g(·; v);w) with the dimension of w being relatively small. Let θMAP = (vMAP, wMAP) be the MAP estimate for θ using batch gradient descent. The idea of our simple approximate Bayesian scheme is to freeze the network weights at the MAP estimate for early layers and perform approximate Bayesian inference for the final layers2. e.g., freeze the parameters of g at vMAP and perform MCMC over w. Throughout the experiments, g : R3 → R3 is a feedforward ReLU block with each hidden layer having 5 hidden units and h : R3 → R3 is either BAx or BReLU(Ax) where A ∈ R3×r, B ∈ Rr×3. We set r = 3. We shall consider 1 or 5 hidden layers for g.
To approximate the Bayes predictive distribution, we perform either the Laplace approximation or the NUTS variant of HMC (Hoffman & Gelman, 2014) in the last two layers, i.e., performing inference over A,B in h(g(·; vMAP);A,B). Note that MCMC is operating in a space of 18 di-
2This is similar in spirit to Kristiadi et al. (2020) who claim that even “being Bayesian a little bit” fixes overconfidence. They approach this via the Laplace approximation for the final layer of a ReLU network. It is also worth noting that Kristiadi et al. (2020) do not attempt to formalise what it means to ”fix overconfidence”; the precise statement should be in terms of G(n).
mensions in this case, which is small enough for us to expect MCMC to perform well. We also implemented the Laplace approximation and NUTS in the last layer only, i.e. performing inference overB in h2(h1(g(·; vMAP);AMAP);B). Further implementation details of these approximate Bayesian schemes are found in Appendix A.5.
From the outset, we expect the Laplace approximation over w = (A,B) to be invalid since the model is singular. We do however expect the last-layer-only Laplace approximation over B to be sound. Next, we expect the MCMC approximation in either the last layer or last two layers to be superior to the Laplace approximations and to the MAP. We further expect the last-two-layers MCMC to have better generalisation than the last-layer-only MCMC since the former is closer to the Bayes predictive distribution. In summary, we anticipate the following performance order for
these five approximate Bayesian schemes (from worst to best): last-two-layers Laplace, last-layeronly Laplace, MAP, last-layer-only MCMC, last-two-layers MCMC.
The results displayed in Figure 1 are in line with our stated expectations above, except for the surprise that the last-layer-only MCMC approximation is often superior to the last-two-layers MCMC approximation. This may arise from the fact that MCMC finds the singular setting in the last-twolayers more challenging. In Figure 1, we clarify the effect of the network architecture by varying the following factors: 1) either 1 or 5 layers in g, and 2) ReLU or identity activation in h. Table 1 is a companion to Figure 1 and tabulates for each approximation scheme the slope of 1/n versus EnG(n), also known as the learning coefficient. The R2 corresponding to the linear fit is also provided. In Appendix A.5, we also show the corresponding results when 1) the data-generating mechanism and the assumed model do not satisfy the condition of realisability and/or 2) the MAP estimate is obtained via minibatch stochastic gradient descent instead of batch gradient descent.

heading: 6 SIMPLE FUNCTIONS AND COMPLEX SINGULARITIES
text: In singular models the RLCT may vary with the true distribution (in contrast to regular models) and in this section we examine this phenomenon in a simple example. As the true distribution becomes more complicated relative to the supposed model, the singularities of the analytic variety of true parameters should become simpler and hence the RLCT should increase (Watanabe, 2009, §7.6). Our experiments are inspired by (Watanabe, 2009, §7.2) where tanh(x) networks are considered and the true distribution (associated to the zero network) is held fixed while the number of hidden nodes is increased.
Consider the model p(y|x,w) in (2) where f(x,w) = c+ ∑H i=1 qi ReLU(〈wi, x〉+bi) is a two-layer ReLU network with weight vector w = ({wi}Hi=1, {bi}Hi=1, {qi}Hi=1, c) ∈ R4H+1 and wi ∈ R2, bi ∈ R, qi ∈ R for 1 ≤ i ≤ H . We let W be some compact neighborhood of the origin.
Given an integer 3 ≤ m ≤ H we define a network sm ∈W and qm(y|x) := p(y|x, sm) as follows. Let g ∈ SO(2) stand for rotation by 2π/m, set w1 = √ g (1, 0)T . The components of sm are the vectors wi = gi−1w1 for 1 ≤ i ≤ m and wi = 0 for i > m, bi = − 13 and qi = 1 for 1 ≤ i ≤ m and bi = qi = 0 for i > m, and finally c = 0. The factor of 13 ensures the relevant parts of the decision boundaries lie within X = [−1, 1]2. We let q(x) be the uniform distribution on X and define qm(x, y) = qm(y|x)q(x). The functions f(x, sm) are graphed in Figure 2. It is intuitively clear that the complexity of these true distributions increases with m.
We let ϕ be a normal distribution N(0, 502) and estimate the RLCTs of the triples (p, qm, ϕ). We conducted the experiments with H = 5, n = 1000. For each m ∈ {3, 4, 5}, Table 2 shows the
estimated RLCT. Algorithm 1 in Appendix A.3 details the estimation procedure which we base on (Watanabe, 2013, Theorem 4). As predicted the RLCT increases with m verifying that in this case, the simpler true distributions give rise to more complex singularities.
Note that the dimension of W is d = 21 and so if the model were regular the RLCT would be 10.5. It can be shown that when m = H the set of true parameters W0 ⊆ W is a regular submanifold of dimension m. If such a model were minimally singular its RLCT would be 12 ((4m + 1) −m) = 1 2 (3m+ 1). In the case m = 5 we observe an RLCT more than an order of magnitude less than the value 8 predicted by this formula. So the functionK does not behave like a quadratic form nearW0.
Strictly speaking it is incorrect to speak of the RLCT of a ReLU network because the function K(w) is not necessarily analytic (Example A.4). However we observe empirically that the predicted linear relationship between Eβw[nLn(w)] and 1/β holds in our small ReLU networks (see the R 2 values in Table 2) and that the RLCT estimates are close to those for the two-layer SiLU network (Hendrycks & Gimpel, 2016) which is analytic (the SiLU or sigmoid weighted linear unit is σ(x) = x(1 + e−τx)−1 which approaches the ReLU as τ → ∞. We use τ = 100.0 in our experiments). The competitive performance of SiLU on standard benchmarks (Ramachandran et al., 2017) shows that the non-analyticity of ReLU is probably not fundamental.

heading: 7 FUTURE DIRECTIONS
text: Deep neural networks are singular models, and that’s good: the presence of singularities is necessary for neural networks with large numbers of parameters to have low generalisation error. Singular learning theory clarifies how classical tools such as the Laplace approximation are not just inappropriate in deep learning on narrow technical grounds: the failure of this approximation and the existence of interesting phenomena like the generalisation puzzle have a common cause, namely the existence of degenerate critical points of the KL function K(w). Singular learning theory is a promising foundation for a mathematical theory of deep learning. However, much remains to be done. The important open problems include:
SGD vs the posterior. A number of works (ŞimŠekli, 2017; Mandt et al., 2017; Smith et al., 2018) suggest that mini-batch SGD may be governed by SDEs that have the posterior distribution as its stationary distribution and this may go towards understanding why SGD works so well for DNNs.
RLCT estimation for large networks. Theoretical RLCTs have been cataloged for small neural networks, albeit at significant effort3 (Aoyagi & Watanabe, 2005b;a). We believe RLCT estimation in these small networks should be standard benchmarks for any method that purports to approximate the Bayesian posterior of a neural network. No theoretical RLCTs or estimation procedure are known for modern DNNs. Although MCMC provides the gold standard it does not scale to large networks. The intractability of RLCT estimation for DNNs is not necessarily an obstacle to reaping the insights offered by singular learning theory. For instance, used in the context of model selection, the exact value of the RLCT is not as important as model selection consistency. We also demonstrated the utility of singular learning results such as (9) and (10) which can be exploited even without knowledge of the exact value of the RLCT.
Real-world distributions are unrealisable. The existence of power laws in neural language model training (Hestness et al., 2017; Kaplan et al., 2020) is one of the most remarkable experimental results in deep learning. These power laws may be a sign of interesting new phenomena in singular learning theory when the true distribution is unrealisable.

heading: A APPENDIX

heading: A.1 NEURAL NETWORKS ARE STRICTLY SINGULAR
text: Many-layered neural networks are strictly singular (Watanabe, 2009, §7.2). The degeneracy of the Hessian in deep learning has certainly been acknowledged in e.g., Sagun et al. (2016) which recognises the eigenspectrum is concentrated around zero and in Pennington & Worah (2018) which deliberately studies the Fisher information matrix of a single-hidden-layer, rather than multilayer, neural network.
We first explain how to think about a neural network in the context of singular learning theory. A feedforward network of depth c parametrises a function f : RN −→ RM of the form
f = Ac ◦ σc−1 ◦Ac−1 · · ·σ1 ◦A1 where the Al : Rdl−1 −→ Rdl are affine functions and σl : Rdl −→ Rdl is coordinate-wise some fixed nonlinearity σ : R −→ R. Let W be a compact subspace of Rd containing the origin, where Rd is the space of sequences of affine functions (Al)cl=1 with coordinates denoted w1, . . . , wd so that f may be viewed as a function f : RN × W −→ RM . We define p(y|x,w) as in (2). We assume the true distribution is realisable, q(y|x) = p(y|x,w0) and that a distribution q(x) on RN is fixed with respect to which p(x, y) = p(y|x)q(x) and q(x, y) = q(y|x)q(x). Given some prior ϕ(w) on W we may apply singular learning theory to the triplet (p, q, ϕ).
By straightforward calculations we obtain
K(w) = 12 ∫ ‖f(x,w)− f(x,w0)‖2q(x)dx (11)
∂2
∂wi∂wj K(w) = ∫ 〈 ∂ ∂wi f(x,w), ∂∂wj f(x,w) 〉 q(x)dx
+ ∫ 〈 f(x,w)− f(x,w0), ∂ 2
∂wi∂wj f(x,w)
〉 q(x)dx (12)
I(w)ij = 1
2(M−3)/2π(M−2)/2
∫ 〈 ∂ ∂wi f(x,w), ∂∂wj f(x,w) 〉 q(x)dx (13)
where 〈−,−〉 is the dot product. We assume q(x) is such that these integrals exist. It will be convenient below to introduce another set of coordinates forW . Let wljk denote the weight from the kth neuron in the (l− 1)th layer to the jth neuron in the lth layer and let blj denote the bias of the jth neuron in the lth layer. Here 1 ≤ l ≤ c and the input is layer zero. Let ulj and alj denote the value of the jth neuron in the lth layer before and after activation, respectively. Let ul and al denote the vectors with values ulj and a l j , respectively. Let dl denote the number of neurons in the lth layer. Then
ulj = dl−1∑ k=1 wljka l−1 k + b l j , 1 ≤ l ≤ c, 1 ≤ j ≤ dl
alj = σ(u l j) 1 ≤ l < c, 1 ≤ j ≤ dl
with the convention that a0 = x is the input and uc = y is the output.
In the case where σ = ReLU the partial derivatives ∂∂wj f do not exist on all of R N . However given w ∈ W we let D(w) denote the complement in RN of the union over all hidden nodes of the associated decision boundary, that is
RN \ D(w) = ⋃
1≤l<c ⋃ 1≤j≤dl {x ∈ RN : ulj(x) = 0} .
The partial derivative ∂∂wj f exists on the open subset {(x,w) : x ∈ D(w)} of R N ×W .
Lemma A.1. Suppose σ = ReLU and there are c > 1 layers. For any hidden neuron 1 ≤ j ≤ dl in layer l with 1 ≤ l < c there is a differential equation{ dl−1∑
k=1
wljk ∂
∂wljk + blj
∂ ∂blj − dl+1∑ i=1 wl+1ij ∂ ∂wl+1ij } f = 0
which holds on D(w) for any fixed w ∈W .
Proof. Without loss of generality assume M = 1, to simplify the notation. Let ei ∈ Rdl+1 denote a unit vector and let H(x) = ddx ReLU(x). Writing ∂f ∂ul+1 for a gradient vector
∂f ∂wl+1ij = 〈 ∂f ∂ul+1 , ∂ul+1 ∂wl+1ij 〉 = 〈 ∂f ∂ul+1 , aljei 〉 = ∂f ∂ul+1i uljH(u l j)
∂f ∂wljk = 〈 ∂f ∂ul+1 , ∂ul+1 ∂wljk 〉 = 〈 ∂f ∂ul+1 , dl+1∑ i=1 wl+1ij a l−1 k H(u l j)ei 〉 = dl+1∑ i=1 ∂f ∂ul+1i wl+1ij a l−1 k H(u l j)
∂f ∂blj = 〈 ∂f ∂ul+1 , ∂ul+1 ∂blj 〉 = 〈 ∂f ∂ul+1 , dl+1∑ i=1 wl+1ij H(u l j)ei 〉 = dl+1∑ i=1 ∂f ∂ul+1i wl+1ij H(u l j).
The claim immediately follows.
Lemma A.2. Suppose σ = ReLU, c > 1 and that w ∈ W has at least one weight or bias at a hidden node nonzero. Then the matrix I(w) is degenerate and if w ∈ W0 then the Hessian of K at w is also degenerate.
Proof. Let w ∈ W be given, and choose a hidden node where at least one of the incident weights (or bias) is nonzero. Then Lemma A.1 gives a nontrivial linear dependence relation ∑ i λi ∂ ∂wi f = 0 as functions on D(w). The rows of I(w) satisfy the same linear dependence relation. At a true parameter the second summand in (12) vanishes so by the same argument the Hessian is degenerate.
Remark A.3. Lemma A.2 implies that every true parameter for a nontrivial ReLU network is a degenerate critical point of K. Hence in the study of nontrivial ReLU networks it is never appropriate to divide by the determinant of the Hessian of K at a true parameter, and in particular Laplace or saddle-point approximations at a true parameter are invalid.
The well-known positive scale invariance of ReLU networks (Phuong & Lampert, 2020) is responsible for the linear dependence of Lemma A.1, in the precise sense that the given differential operator is the infinitesimal generator (Boothby, 1986, §IV.3) of the scaling symmetry. However, this is only one source of degeneracy or singularity in ReLU networks. The degeneracy, as measured by the RLCT, is much lower than one would expect on the basis of this symmetry alone (see Section 6).
Example A.4. In general the KL functionK(w) for ReLU networks is not analytic. For the minimal counterexample, let q(x) be uniform on [−N,N ] and zero outside and consider
K(b) = ∫ q(x)(ReLU(x− b)− ReLU(x))2dx .
It is easy to check that up to a scalar factor
K(b) =
{ − 23b
3 + b2N 0 ≤ b ≤ N − 13b
3 + b2N −N ≤ b ≤ 0 so that K is C2 but not C3 let alone analytic.

heading: A.2 REDUCED RANK REGRESSION
text: For reduced rank regression, the model is
p(y|x,w) = 1 (2πσ2)N/2 exp
( − 1
2σ2 |y −BAx|2
) ,
where x ∈ RM , y ∈ RN , A an M ×H matrix and B an H×N matrix; the parameter w denotes the entries of A and B, i.e. w = (A,B), and σ > 0 is a parameter which for the moment is irrelevant.
If the true distribution is realisable then there is w0 = (A0, B0) such that q(y|x) = p(y|x,w0). Without loss of generality assume q(x) is the uniform density. In this case the KL divergence from p(y|x,w) to q(y|x) is
K(w) = ∫ q(y|x) log q(y|x)
p(y|x,w) dxdy = ‖BA−B0A0‖2 (1 + E(w))
where the errorE is smooth andE(w) = O(‖BA−B0A0‖2) in any region where ‖BA−B0A0‖ < C, so K(w) is equivalent to ‖BA − B0A0‖2. We write K(w) = ‖BA − B0A0‖2 for simplicity below.
Now assume that B0A0 is symmetric and that B0 is square, i.e. N = H . Then the zero locus of K(w) is explicitly given as follows W0 = {(A,B) : detB 6= 0 and A = B−1B0A0}. It follows that W0 is globally a graph over GL(H;R). Indeed, the set (B−1B0A0, B) with B ∈ GL(H;R) is exactly W0. Thus W0 is a smooth H2-dimensional submanifold of RH
2 × RH×M . To prove that W0 is minimally singular in the sense of Section 4 it suffices to show that rank(D2A,BK) ≥ HM whereD2A,BK denotes the Hessian, but as it is no more difficult to do so, we find explicit local coordinates (u, v) near an arbitrary point (A,B) ∈ W0 for which {v = 0} = W0 and K(u, v) = a(u, v)|u|2 in this neighborhood, where a is a C∞ function with a ≥ c > 0 for some c. Write
A(v) = (B + v)−1B0A0.
Then u, v 7→ (A(v) + u,B + v) gives local coordinates on RH2 × RH×M near (A,B), and K(u, v) = |(B + v)((B + v)−1B0A0 + u)−B0A0|
= |B0A0 + (B + v)u−B0A0|2
= |(B + v)u|2, so for v sufficiently small (and hence B + v invertible) we can take a(u, v) = |(B + v)u|2/|u|2.

heading: A.3 RLCT ESTIMATION
text: In this section we detail the estimation procedure for the RLCT used in Section 6. Let Ln(w) be the negative log likelihood as in (3). Define the data likelihood at inverse temperature β > 0 to be
Algorithm 1 RLCT via Theorem 4 in Watanabe (2013) Input: range of β’s, set of training sets T each of size n, approximate samples {w1, . . . , wR} from pβ(w|Dn) for each training set Dn and each β for training set Dn ∈ T do
for β in range of β’s do Approximate Eβw[nLn(w)] with 1R ∑R i=1 nLn(wr) where w1, . . . , wR are approximate sam-
ples from pβ(w|Dn) end for Perform generalised least squares to fit λ in (18), call result λ̂(Dn)
end for Output: 1|T | ∑ Dn∈T λ̂(Dn)
pβ(Dn|w) = Πni=1p(yi|xi, w)β which can also be written pβ(Dn|w) = exp(−βnLn(w)). (14) The posterior distribution, at inverse temperature β, is defined as
pβ(w|Dn) = Πni=1p(yi|xi, w)βϕ(w)∫ W Πni=1p(yi|xi, w)βϕ(w) = pβ(Dn|w)ϕ(w) pβ(Dn) (15)
where ϕ is the prior distribution on the network weights w and pβ(Dn) = ∫ W pβ(Dn|w)ϕ(w) dw (16) is the marginal likelihood of the data at inverse temperature β. Finally, denote the expectation of a random variable R(w) with respect to the tempered posterior pβ(w|Dn) as
Eβw[R(w)] = ∫ W R(w)pβ(w|Dn) dw (17)
In the main text, we drop the superscript in the quantities (14), (15), (16), (17) when β = 1, e.g., p(Dn) rather than p1(Dn). Assuming the conditions of Theorem 4 in Watanabe (2013) hold, we have
Eβw[nLn(w)] = nLn(w0) + λ
β + Un
√ λ
2β +Op(1) (18)
where β0 is a positive constant and Un is a sequence of random variables satisfying EnUn = 0. In Algorithm 1, we describe an estimation procedure for the RLCT based on the asymptotic result in (18).
For the estimates in Table 2 the a posteriori distribution was approximated using the NUTS variant of Hamiltonian Monte Carlo (Hoffman & Gelman, 2014) where the first 1000 steps were omitted and 20, 000 samples were collected. Each λ̂(Dn) estimate in Algorithm 1 was performed by linear regression on the pairs {(1/βi,Eβiw [nLn(w)])}5i=1 where the five inverse temperatures βi are centered on the inverse temperature 1/ log(20000).

heading: A.4 CONNECTION BETWEEN RLCT AND GENERALISATION
text: For completeness, we sketch the derivation of (9) which gives the asymptotic expansion of the average generalisation error EnG(n) of the Bayes prediction distribution in singular models. The exposition is an amalgamation of various works published by Sumio Watanabe, but is mostly based on the textbook (Watanabe, 2009).
To understand the connection between the RLCT and G(n), we first define the so-called Bayes free energy as
F (n) = − log p(Dn)
whose expectation admits the following asymptotic expansion (Watanabe, 2009): EnF (n) = EnnSn + λ log n+ o(log n)
where Sn = − 1n ∑n i=1 log q(yi|xi) is the entropy. The expected Bayesian generalisation error is related to the Bayes free energy as follows EnG(n) = EF (n+ 1)− EF (n) Then for the average generalisation error, we have EnG(n) = λ/n+ o(1/n). (19) Since models with more complex singularities have smaller RLCTs, this would suggest that the more singular a model is, the better its generalisation (assuming one uses the Bayesian predictive distribution for prediction). In this connection it is interesting to note that simpler (relative to the model) true distributions lead to more singular models (Section 6).

heading: A.5 DETAILS FOR GENERALISATION ERROR EXPERIMENTS
text: Simulated data The distribution of x ∈ R3 is set to q(x) = N(0, I3). In the realisable case, y ∈ R3 is drawn according to q(y|x) = p(y|x, θ0). In the nonrealisable setting, we set q(y|x) ∝ exp{−||y− hw0(x)||2/2}, where w0 = (A0, B0) is drawn according to the PyTorch model initialisation of h. MAP training The MAP estimator is found via gradient descent using the mean-squared-error loss with either the full data set or minibatch set to 32. Training was set to 5000 epochs. No form of early stopping was employed.
Calculating the generalisation error Using a held-out-test set Tn′ = {(x′i, y′i)}n ′
i=1, we calculate the average generalisation error as
1
n′ n′∑ i=1 log q(y′i|x′i)− En 1 n′ n′∑ i=1 log q̂n(y ′ i|x′i) (20)
Assume the held-out test set is large enough so that the difference between EnG(n) and (20) is negligible. We will refer to them interchangeably as the average generalisation error. In our experiments we use n′ = 10, 000 and 30 draws of the dataset Dn to estimate En. Last layer(s) inference Without loss of generality, we discuss performing inference in the w parameters of h while freezing the parameters of g at the MAP estimate. The steps easily extend to performing inference over the final layer only of f = h ◦ g. Let x̃i = gvMAP(xi). Define a new transformed dataset D̃n = {(x̃i, yi)}ni=1. We take the prior on w to be standard Gaussian. Define the posterior over w given D̃n as: p(w|D̃n) ∝ p(D̃n|w)ϕ(w) = Πni=1 exp{−||yi − hw(x̃i)||2/2}ϕ(w) (21) Define the following approximation to the Bayesian predictive distribution
p̃(y|x,Dn) = ∫ p(y|x, (vMAP, w))p(w|D̃n) dw.
Let w1, . . . , wR be some approximate samples from p(w|D̃n). Then we approximate p̃(y|x,Dn) with
1
R R∑ r=1 p(y|x, (vMAP, wr))
where R is a large number, set to 1000 in our experiments. We consider the Laplace approximation and the NUTS variant of HMC for drawing samples from p(w|D̃n):
• Laplace in the last layer(s) Recall θMAP = (vMAP, wMAP) is the MAP estimate for fθ trained with the data Dn. With the Laplace approximation, we draw w1, . . . wR from the Gaussian
N(wMAP,Σ)
where Σ = (−∇2 log p(w|D̃n)|wMAP)−1 is the inverse Hessian4 of the negative log posterior evaluated at the MAP estimate of the mode. • MCMC in the last layer(s) We used the NUTS variant of HMC to draw samples from
(21) with the first 1000 samples discarded.. Our implementation used the pyro package in PyTorch.
4Following Kristiadi et al. (2020), the code for the exact Hessian calculation is borrowed from https: //github.com/f-dangel/hbp

Review summaries:
This paper is more like a review of singular learning theory and its implication on deep learning. The authors point out that deep neural networks are singular models and ways to characterize generalization error for regular models cannot produce satisfactory results in this setting. Then the authors introduce the singular learning theory, which has been developed for decades. Then, a series of topics for deep learning, such as flatness and generalization, are studied within the framework singular learning theory, with a combination of theoretical analysis and numerical experiments. The paper is clearly written and well organized.