Paper deep dive
On the Equivalence of Random Network Distillation, Deep Ensembles, and Bayesian Inference
Moritz A. Zanger, Yijun Wu, Pascal R. Van der Vaart, Wendelin Böhmer, Matthijs T. J. Spaan
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 95%
Last extracted: 7/20/2026, 3:26:40 PM
Summary
This paper establishes theoretical connections between Random Network Distillation (RND), Deep Ensembles, and Bayesian Inference within the Neural Tangent Kernel (NTK) framework for infinite-width neural networks. It proves that RND's squared self-predictive error is equivalent to the predictive variance of a deep ensemble. Furthermore, by constructing a specific RND target function, the authors show that RND error distributions can mirror the posterior predictive distribution of Bayesian inference, enabling a novel posterior sampling algorithm that generates i.i.d. samples from the exact Bayesian posterior.
Entities (7)
Relation Signals (5)
Neural Tangent Kernel → usedforanalysisof → Random Network Distillation
confidence 96% · analyzing RND within the neural tangent kernel framework in the limit of infinite network width.
Random Network Distillation → equivalentto → Deep Ensembles
confidence 95% · The uncertainty signal from RND -- its squared self-predictive error -- is equivalent to the predictive variance of a deep ensemble.
Random Network Distillation → mirrors → Bayesian Inference
confidence 94% · the RND error distribution can be made to mirror the centered posterior predictive distribution of Bayesian inference with wide neural networks.
Posterior Sampling Algorithm → generatessamplesfrom → Bayesian Inference
confidence 93% · generates i.i.d. samples from an exact Bayesian posterior predictive distribution using this modified Bayesian RND model.
Infinite Network Width → enablesequivalenceof → Random Network Distillation
confidence 92% · in the limit of infinite network width... the squared prediction errors of standard RND coincide exactly with the variance of a deep ensemble.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Uncertainty quantification is central to safe and efficient deployments of deep learning models, yet many computationally practical methods lack lacking rigorous theoretical motivation. Random network distillation (RND) is a lightweight technique that measures novelty via prediction errors against a fixed random target. While empirically effective, it has remained unclear what uncertainties RND measures and how its estimates relate to other approaches, e.g. Bayesian inference or deep ensembles. This paper establishes these missing theoretical connections by analyzing RND within the neural tangent kernel framework in the limit of infinite network width. Our analysis reveals two central findings in this limit: (1) The uncertainty signal from RND -- its squared self-predictive error -- is equivalent to the predictive variance of a deep ensemble. (2) By constructing a specific RND target function, we show that the RND error distribution can be made to mirror the centered posterior predictive distribution of Bayesian inference with wide neural networks. Based on this equivalence, we moreover devise a posterior sampling algorithm that generates i.i.d. samples from an exact Bayesian posterior predictive distribution using this modified \textit{Bayesian RND} model. Collectively, our findings provide a unified theoretical perspective that places RND within the principled frameworks of deep ensembles and Bayesian inference, and offer new avenues for efficient yet theoretically grounded uncertainty quantification methods.
Tags
Links
- Source: https://arxiv.org/abs/2602.19964v2
- Canonical: https://arxiv.org/abs/2602.19964v2
Trouble viewing inline? Open PDF directly →
Full Text
101,341 characters extracted from source content.
Expand or collapse full text
On the Equivalence of Random Network Distillation, Deep Ensembles, and Bayesian Inference Moritz A. Zanger Delft University of Technology Delft, 2628 XE, The Netherlands Yijun Wu Delft University of Technology Delft, 2628 XE, The Netherlands Pascal R. Van der Vaart Delft University of Technology Delft, 2628 XE, The Netherlands Wendelin Boehmer Delft University of Technology Delft, 2628 XE, The Netherlands Matthijs T. J. Spaan Delft University of Technology Delft, 2628 XE, The Netherlands Abstract Uncertainty quantification is central to safe and efficient deployments of deep learning models, yet many computationally practical methods lack lacking rigorous theoretical motivation. Random network distillation (RND) is a lightweight technique that measures novelty via prediction errors against a fixed random target. While empirically effective, it has remained unclear what uncertainties RND measures and how its estimates relate to other approaches, e.g., Bayesian inference or deep ensembles. We establish these missing theoretical connections by analyzing RND within the neural tangent kernel framework in the limit of infinite network width. Our analysis reveals two central findings in this limit: (1) The uncertainty signal from RND—its squared self-predictive error—is equivalent to the predictive variance of a deep ensemble. (2) By constructing a specific RND target function, we show that the RND error distribution can be made to mirror the centered posterior predictive distribution of Bayesian inference with wide neural networks. Based on this equivalence, we moreover devise a posterior sampling algorithm that generates i.i.d. samples from an exact Bayesian posterior predictive distribution using this modified Bayesian RND model. Collectively, our findings provide a unified theoretical perspective that places RND within the principled frameworks of deep ensembles and Bayesian inference, and offer new avenues for efficient yet theoretically grounded uncertainty quantification methods. 1 Introduction Quantifying predictive uncertainty remains a cornerstone of reliable machine learning and underpins applications from safe robotics to efficiently exploring agents and autonomous scientific discovery. Bayesian inference is widely regarded as a theoretical gold‐standard to this end [neal1996bayesian, goan2020bayesian] but its application to neural networks is typically intractable in practice, requiring approximations of simplified posteriors through variational inference [VI, kingmaAutoEncodingVariationalBayes2014, galDropoutBayesianApproximation2016a, blei2017variational] or complex sampling mechanisms through Markov chain Monte Carlo approaches [MCMC, chen2014stochastic, liu2016stein, garriga2021exact]. Deep ensembles [dietterich2000ensemble, lakshminarayananSimpleScalablePredictive2017] on the other hand maintain several independently initialized models to quantify predictive variance as uncertainty. Due to their simplicity and relative practical reliability, deep ensembles have become a widely established alternative to Bayesian approaches for uncertainty quantification in deep learning [abdarReviewUncertaintyQuantification2021]. However, both ensemble methods and approximate Bayesian methods typically incur substantial computational and memory costs, in particular for larger-scale models, motivating more efficient alternatives. RND [burda2018exploration] offers one such approach: by training a predictor network to mimic the outputs of a fixed, randomly initialized target network, RND produces a simple novelty or uncertainty signal via the squared prediction error. Random network distillation (RND) has seen empirical success in exploration, out-of-distribution detection, and continual learning [burda2018exploration, nikulin2023anti, matthews2024craftax], yet the theoretical understanding of the nature of its uncertainty estimates remains blurry. In particular, it is unclear how—or whether—the RND error relates to the principled uncertainties produced for example by Bayesian inference or deep ensembles. In this paper, we establish these missing theoretical connections by analyzing random network distillation in the idealized setting of infinite network width. In particular, we establish a Gaussian process (GP) interpretation of the self-predictive RND errors in the limit of infinitely wide neural networks, drawing on Neural Tangent Kernel (NTK) theory [jacotNeuralTangentKernel2020, leeWideNeuralNetworks2020]. Our three main contributions are: 1. Ensemble equivalence with Standard RND: We prove that, in the idealized infinite width limit, the squared prediction errors of standard RND coincide exactly with the variance of a deep ensemble. 2. Posterior equivalence with Bayesian RND: By engineering the RND target function, we design a Bayesian RND variant whose error distribution matches that of the exact Bayesian posterior predictive distribution of a neural network in the limit of infinite width. 3. Posterior sampling with Bayesian RND: Based on a multi-headed Bayesian RND model, we devise a posterior sampling algorithm that produces i.i.d. samples of the exact Bayesian posterior predictive distribution of neural networks in the limit of infinite width. This unifying perspective on the uncertainty estimates produced by RND, deep ensembles, and Bayesian inference provides a novel understanding and theoretical support for the empirical effectiveness of RND and suggests avenues for future research directions towards principled Bayesian inference with minimal computational overhead. 2 Preliminaries We begin by establishing notation, defining RND formally, and briefly introducing the theoretical framework used in our analysis. In our analysis, we consider fully connected neural networks f(x;θt)f(x; _t) of L layers of widths n1,…,nL=n_1,…,n_L=n, parametrized by θt _t at time t. The forward computation of such networks is defined recursively with zil(x;θt≤l)z_i^l(x;θ^≤ l_t) denoting the i-th output of layer l and zil(x,θt≤l)=σbbil+σwnl−1∑j=1nl−1wijlxjl(x)xjl(x)=ϕ(zjl−1(x;θt≤l−1)), splitz_i^l(x,θ^≤ l_t)&= _bb_i^l+ _w n_l-1 Σ _ j=1 n_l-1w_ij^lx_j^l(x)\\ x_j^l(x)&=φ(z_j^l-1(x;θ^≤ l-1_t))\,, split (1) where θt≤lθ^≤ l_t denotes the parameters w1,b1,…,wl,bl\w^1,b^1,…,w^l,b^l\ up to layer l, σb _b and σw _w denote scaling parameters of the forward computation, and ϕ:ℝ→ℝφ:R R is a Lipschitz-continuous nonlinearity. In Eq. (1), n0=din_0=d_in and x1(x)=x^1(x)=x. The output of a scalar-output neural network is then given by f(x;θt)=zL(x;θt≤L)f(x; _t)=z^L(x;θ^≤ L_t) . We furthermore assume that parameters are initialized i.i.d. from a normal distribution θ0∼(0,I) _0 (0,I)111Also known as NTK-parametrization, where variance scalings σb _b and σw _w affect both forward and gradient computations, yielding well-behaved gradients in the infinite-width limit.). For convenience, we will sometimes overload notation to concatenate function outputs, for example indicating a set =xi∈ℝdini=1NDX=\x_i ^d_in\_i=1^N_D and the corresponding function output as a column vector f(;θt)=(f(xi;θt))i=1NDf(X; _t)=(f(x_i; _t))_i=1^N_D, where f(;θt)∈ℝND×Kf(X; _t) ^N_D× K or matrix-valued identities Σ(,)=(Σ(xi,xj))i,j=1ND (X,X)=( (x_i,x_j))_i,j=1^N_D , where Σ(,)∈ℝND×ND (X,X) ^N_D× N_D. For conciseness our notation will furthermore use a shorthand for covariance and kernel matrices denoting Σ≡Σ(,) _XX≡ (X,X). In the following we briefly review methods pertinent to this work. Random network distillation. Random network distillation [burda2018exploration] is an uncertainty quantification technique that employs two neural networks of identical architecture: A fixed, randomly initialized target network g(x;ψ0):ℝdin→ℝKg(x; _0):R^d_in ^K , and a predictor network u(x;ϑt)u(x; _t), where parameters ϑt _t are subject to optimization via gradient descent. The predictor is trained to minimize the expected squared difference to the target network’s output on a set of data points =xi∈ℝdini=1NDX=\x_i ^d_in\_i=1^N_D ℒrnd(ϑt)=12‖u(;ϑt)−g(;ψ0)‖22. _rnd( _t)= 12\|u(X; _t)-g(X; _0)\|^2_2\,. (2) It is common to design RND with a multi headed architecture with output dimension K and individual output heads ui(x;ϑt)i=1K\u_i(x; _t)\_i=1^K, and gi(x;ψ0)i=1K\g_i(x; _0)\_i=1^K, where the sum of squared prediction errors ϵi(x;ϑt,ψ0)=ui(x;ϑt)−gi(x;ψ0) _i(x; _t, _0)=u_i(x; _t)-g_i(x; _0) at a test point x serves as an uncertainty signal ϵ2(x;ϑt,ψ0)=1K∑i=1K(ui(x;ϑt)−gi(x;ψ0))2. ε^2(x; _t, _0)= 1K Σ _ i=1 K (u_i(x; _t)-g_i(x; _0) )^2\,. (3) Gaussian processes. In our analysis, we will frequently use GPs to model distributions over random functions: A univariate GP [rasmussen2006gp] defines a distribution over functions f0∼(μ0,Σ0)f^0 (μ^0, ^0) characterized by a mean function μ0:ℝdin→ℝμ^0:R^d_in R and a covariance (kernel) function Σ0:ℝdin×ℝdin→ℝ ^0:R^d_in×R^d_in R such that f0(T)f_0(X_T) follows a multivariate Gaussian distribution f0(T)∼(μ0(T),Σ0(T,T))f_0(X_T) (μ^0(X_T), ^0(X_T,X_T)) for any finite set of evaluation points T=xiTesti=1NTX_T=\x_i^Test\_i=1^N_T. We can condition a prior GP (μ0(T),Σ0(T,T))N(μ^0(X_T), ^0(X_T,X_T)) on training data =xii=1NDX=\x_i\_i=1^N_D and labels =yii=1NDY=\y_i\_i=1^N_D to obtain a posterior GP whose posterior predictive distribution is Gaussian with mean and covariance given by μ(T)=μ0(T)+ΣT0(Σ0)−1(−μ0()),ΣTT=ΣTT0−ΣT0(Σ0)−1ΣT0. splitμ(X_T)&=μ^0(X_T)+ ^0_X_TX( ^0_XX)^-1 (Y-μ^0(X) ),\\ _X_TX_T&= ^0_X_TX_T- ^0_X_TX( ^0_XX)^-1 ^0_XX_T\,. split (4) Learning dynamics with infinite width. We turn to analytical tools to establish solutions to the learning dynamics of neural networks in the limit of infinite width n→∞n→∞. Within this setting, we consider the training dynamics under gradient flow, the continuous-time limit of gradient descent dtθt=−∇θℒ(θt) ddt _t=- _θL( _t). Under gradient flow with a square loss ℒ(θt)=12‖f(;θt)−‖22L( _t)= 12\|f(X; _t)-Y\|_2^2, the evolution of the N f is described by a differential equation in function space dtf(x;θt) ddtf(x; _t) =∇θf(x;θt)⊤dtθt = _θf(x; _t) ddt _t =−∇θf(x;θt)⊤∇θf(;θt)(f(;θt)−) =- _θf(x; _t) _θf(X; _t)(f(X; _t)-Y) ≡−Θt(x,)(f(;θt)−). ≡- _t(x,X)(f(X; _t)-Y)\,. (5) The above learning dynamics are governed by a gradient similarity function, called the neural tangent kernel [NTK, jacotNeuralTangentKernel2020], Θt(x,x′)=∇θf(x;θt)⊤∇θf(x′;θt) _t(x,x )= _θf(x; _t) _θf(x ; _t). While this inner product is dynamic and therefore intractable in general, the limit of infinite network width yields a remarkable simplification: 1.) due to large number effects, the inner product kernel Θ0(x,x′) _0(x,x ) at initialization is deterministic despite the random initialization of θ0 _0; 2.) Θt(x,x′) _t(x,x ) remains constant throughout t under gradient flow [jacotNeuralTangentKernel2020, leeWideNeuralNetworks2020]. In particular, this means limn→∞Θ0(x,x′)=limn→∞Θt(x,x′)≡Θ(x,x′) _n→∞ _0(x,x )= _n→∞ _t(x,x )≡ (x,x ) and converts Eq. 2 into a linear ordinary differential equation, which can be solved analytically. It can be shown that, under mild conditions, f(x;θt)f(x; _t) converges to the kernel regression solution [see jacotNeuralTangentKernel2020, and Appendix B.1] f(x;θ∞) f(x; _∞) =f(x;θ0)−ΘxΘ−1(−f(;θ0)), =f(x; _0)- _xX _XX^-1 (Y-f(X; _0) )\,, (6) Moreover, lee2017deep show that both f(x;θ0)f(x; _0) and f(x;θ∞)f(x; _∞) are indeed GPs described by the neural network Gaussian process [NNGP, lee2017deep] f(x;θ0)∼(0,κxx′)f(x; _0) (0, _x ) and the converged GP defined in Theorem 2.1. Theorem 2.1. [leeWideNeuralNetworks2020](Distribution of post-convergence neural network functions) Let f(T;θ∞)f(X_T; _∞) be a N as defined in Eq.(1), and let TX_T be testpoints. For random initializations θ0∼(0,I) _0 (0,I), and in the limit n→∞n→∞, f(T;θ∞)f(X_T; _∞) distributes as a Gaussian with mean and covariance given by [f(T,θ∞)] [f(X_T, _∞)] =ΘTΘ−1, = _X_TX _XX^-1Y\,, ΣTTf(θ∞) ^f_X_TX_T( _∞) =κTT+ΘTΘ−1κΘ−1ΘT = _X_TX_T+ _X_TX _XX^-1 _XX _XX^-1 _XX_T −(ΘTΘ−1κT+h.c.), - ( _X_TX _XX^-1 _XX_T+h.c. )\,, where h.c. is the Hermitian conjugate of the preceding term. See also Appendix B.1.2 or leeWideNeuralNetworks2020. Note that the GP described in Theorem 2.1 represents the law by which an infinite ensemble of infinitely wide neural networks from i.i.d. initializations distributes after training on (,)(X,Y), but—as is—permits no Bayesian posterior interpretation, which is of the canonical form described in Eq. 4. 3 Equivalence of Random Network Distillation & Deep Ensembles We proceed to characterize formally the relationship between the error signals as measured by random network distillation and the predictive variance of deep neural network ensembles. Before treating multivariate output dimensions in section 3.1, we first consider scalar function outputs for simplicity, i.e. f,u,g:ℝdin→ℝf,u,g:R^d_in R with K=1K=1. This setup involves training a predictor u(x;ϑt)u(x; _t) to match a fixed random target function g(x;ψ0)g(x; _0). Intuitively, the expected errors ought to vanish for training points in X and remain non-zero elsewhere, inheriting the randomness and generalization behaviors of the functions u and g. Owing to the linear training dynamics in the NTK regime, the dynamics of the error evolution dtϵ(x;ϑt,ψ0) ddtε(x; _t, _0) become akin to those outlined in Eq. (2) as dtϵ(x;ϑt,ψ0) ddtε(x; _t, _0) =∇θu(x;ϑt)⊤dtϑt = _θu(x; _t) ddt _t =−∇ϑu(x;ϑt)⊤∇ϑℒrnd(ϑt) =- _ u(x; _t) _ L_rnd( _t) (7) =−Θt(x,)ϵ(x;ϑt,ψ0). =- _t(x,X)ε(x; _t, _0)\,. We then draw on the results of Theorem 2.1 to provide a probabilistic description of the self-predictive errors ϵ(x;ϑ∞,ψ0)ε(x; _∞, _0) of a converged RND model in the limit of infinite network width. Theorem 3.1. (Distribution of post-convergence RND errors) Under NTK parametrization, let u(x;ϑ∞)u(x; _∞) be a converged prediction network in t→∞t→∞, with data X and fixed target network g(;ψ0)g(X; _0). Let parameters ϑ0,ψ0 _0, _0 be drawn i.i.d. ϑ0,ψ0∼(0,I) _0, _0 (0,I), with the resulting NNGP u(x;ϑ0)∼(0,κu(x,x′))u(x; _0) (0,κ^u(x,x )) and g(x;ψ0)∼(0,κg(x,x′))g(x; _0) (0,κ^g(x,x )). The post-convergence RND error ϵ(T;ϑ∞,ψ0)ε(X_T; _∞, _0) is Gaussian with zero mean and covariance [ϵ(T,ϑ∞,ψ0)]=0, [ε(X_T, _∞, _0)]=0\,, ΣTTϵ(ϑ∞,ψ0)=κTTϵ+ΘTΘ−1κϵΘ−1ΘT ^ε_X_TX_T( _∞, _0)\!=\!κ^ε_X_TX_T\!+ _X_TX _XX^-1κ^ε_XX _XX^-1 _XX_T −(ΘTΘ−1κTϵ+h.c.), - ( _X_TX _XX^-1κ^ε_XX_T\!+h.c. )\,, where κxx′ϵ=κxx′u+κxx′gκ^ε_x =κ^u_x +κ^g_x is the covariance kernel of initialization errors ϵ(x;ϑ0,ψ0)=u(x;ϑ0)−g(x;ψ0)ε(x; _0, _0)=u(x; _0)-g(x; _0). Proof sketch. The error function u(x;ϑ∞)−g(x;ψ0)u(x; _∞)-g(x; _0) is a sum of the random post-convergence function u(x;ϑ∞)u(x; _∞) and the fixed random target function g(x;ψ0)g(x; _0). The latter g(x;ψ0)g(x; _0) is known to follow the NNGP. By the linearity of NTK learning dynamics, the online function u(x;ϑ∞)u(x; _∞) is an affine transformation of its initialization u(x;ϑ0)u(x; _0), which itself follows the NNGP. Moreover, this affine transformation is independent of g or ψ0 _0, such that the error ϵ(x;ϑ∞,ψ0)ε(x; _∞, _0) is a sum of two independent GPs and therefore a GP itself. The resulting GP has zero-mean and covariance with an altered prior NNGP kernel κϵ(x,x′)κ^ε(x,x ) composed of the online prior kernel κxx′uκ^u_x and the target prior kernel κxx′gκ^g_x . See also Appendix B.1.3. Corollary 3.2. (Equivalence in expectation between RND errors and ensemble variance) Under the conditions of Theorem 3.1, let ϵ(x;ϑ∞,ψ0)ε(x; _∞, _0) be the error function of a converged RND network with data X. Moreover, for a regression problem on X for some labels Y, let [f(x;θ∞)]V[f(x; _∞)] denote the variance of converged N functions random initializations. Furthermore, suppose an architectural equivalence between f, u, and g and i.i.d. parameter initialization θ0,ϑ0,ψ0∼(0,I) _0, _0, _0 (0,I). The expected norm of the RND error ϵ2(x;ϑ∞,ψ0)ε^2(x; _∞, _0) then coincides with the ensemble variance ϑ0,ψ0[ϵ2(x;ϑ∞,ψ0)]=θ0[f(x;θ∞)] _ _0, _0 [ε^2(x; _∞, _0) ]=V_ _0[f(x; _∞)] (8) Proof sketch. Corollary 3.2 follows straighforwardly from Theorem 3.1 by using κu(x,x′)=κg(x,x′)κ^u(x,x )=κ^g(x,x ). Taking the trace of the covariance matrix and dividing by 22, we recover the predictive ensemble variance θ0[f(x;θ∞)]V_ _0[f(x; _∞)]. Theorem 3.1 and Corollary 3.2 formally show that, for an architectural equivalence between ensemble, predictor and target network, the expected RND errors directly quantify the predictive variance of the corresponding infinite ensemble model described by Theorem 2.1. To the best of our knowledge, it is the first formal analysis of random network distillation in the NTK regime and reveals a first theoretical motivation for the popular algorithm: in the idealized infinite-width setting, expected RND errors exactly quantify the variance of deep ensembles for any input x. 3.1 Multi-Headed Random Network Distillation The analysis thus far has considered the average behavior of scalar network outputs for simplicity. While insightful in its own right, this setting does not reflect most common practical implementations of random network distillation and instead, if taken literally, would imply an ensemble of random network distillation models. To connect with common practical implementations that typically use multi-headed architectures for enhanced reliability and efficiency, we now seek to incorporate the probabilistic relation between different function outputs fi(x;θt)f_i(x; _t) and fj(x′;θt)f_j(x ; _t) of a N with shared hidden layers in the infinite-width limit. The result below identifies this relationship simply as a statistical independence between the different random network outputs fi(x;θt)f_i(x; _t) and fj(x′;θt)f_j(x ; _t) for any time t during gradient flow optimization. Proposition 3.3. (Independence of N functions) Under NTK parametrization and in the limit n→∞n→∞, the random functions fi(x;θt)f_i(x; _t) of a N with K output dimensions and shared hidden layers are mutually independent with covariance Σxx′ij(θt) _x ^ij( _t) =[fi(x;θt)fj(x′;θt)]=Σxx′f(θt)i=j,0i≠j, =E[f_i(x; _t)f_j(x ; _t)]= cases ^f_x ( _t)&\,i=j\,,\\ 0&\,i≠ j\,, cases on the interval t∈[0,∞)t∈[0,∞). Proof sketch. The property follows from known results that state the independence between output dimensions of the NNGP kernel κ and the NTK Θ [arora2019exact, lee2017deep, jacotNeuralTangentKernel2020]. For both kernels, the proof proceeds by induction, where the independence property between output dimensions is propagated layer-wise. The induction start is equal for both kernels, where first layer outputs, as well as gradients are linear transformations of the Gaussian first-layer weights. Both the NNGP and NTK permit a recursive formulation, through which the independence property can be propagated layer-wise, constituting the induction step. Combined with the learning dynamics of wide NNs, we can conclude that the individual function outputs of a multi-headed N, too, are statistically independent for any time t on the interval [0,∞)[0,∞). See Appendix B.1.4 or lee2017deep and jacotNeuralTangentKernel2020. Notably, this decoupling holds despite the shared hidden layers and is an artifact of the learning dynamics exhibited in the infinite width limit and the NTK regime. In the absence of feature learning, output functions become statistically independent despite sharing a network body. By virtue of this independence property, a translation of the earlier obtained single-function results on RND error distributions (Theorem 3.1 and Corollary 3.2) to the multi-headed setting is straightforward. Our next result thus establishes an equivalence between the errors of the multi-headed RND algorithm, a widely used architecture in practice, and the variance of a finite-sized deep ensemble. Theorem 3.4. (Distributional equivalence between multi-headed RND and finite deep ensembles) Under the conditions of Theorem 3.1, let ui(x;ϑ∞),gi(x;ψ0)u_i(x; _∞),g_i(x; _0) be the i-th output of predictor and target networks respectively with K output dimensions. Denote their sample mean RND error ϵ¯2(x;ϑ∞,ψ0)=1K∑i=1Kϵi2(x;ϑ∞,ψ0) ε^2(x; _∞, _0)= 1K _i=1^K _i^2(x; _∞, _0). Moreover, let f(x;θ∞i)i=1K+1\f(x;θ^i_∞)\_i=1^K+1 be an ensemble of K+1K+1 NNs from i.i.d. initial draws θ0 _0. Denote its sample variance σ¯f2(x;θ∞i…K+1)=1K∑i=1K+1(f(x;θ∞i)−1K+1∑j=1K+1f(x;θ∞j))2 σ_f^2(x; _∞^i… K+1)= 1K _i=1^K+1(f(x;θ^i_∞)- 1K+1 _j=1^K+1f(x;θ^j_∞))^2. The sample mean RND error and sample ensemble variance distribute to the same law 12ϵ¯2(x;ϑ∞,ψ0)=σ¯f2(x;θ∞i…K+1), 12 ε^2(x; _∞, _0) D= σ_f^2(x; _∞^i… K+1)\,, (9) where = D= indicates an equality in distribution, namely by a scaled Chi-squared distribution σ¯f2(x;θ∞i…K+1)∼Σxxf(θ∞)Kχ2(K) σ_f^2(x; _∞^i… K+1) ^f_x( _∞)Kχ^2(K) with scale Σxxf(θ∞) ^f_x( _∞) given by the analytical variance as given in Theorem 2.1. Proof sketch. By Proposition 3.3, the function heads ui(x;ϑ∞)i=1K\u_i(x; _∞)\_i=1^K are K independent predictors, each trained to match their independent targets gi(x;ψ0)g_i(x; _0). Thus, the errors ϵi(x;ϑ∞,ψ0)i=1K\ _i(x; _∞, _0)\_i=1^K are i.i.d. samples from the error distribution outlined in Proposition 3.2. In particular, ϵ¯2 ε^2 is the empirical mean of i.i.d. samples from a Gaussian which is known to be Chi-squared distributed. Similarly, we have that the ensemble f(x;θ∞i)i=1K+1\f(x;θ^i_∞)\_i=1^K+1 are K+1K+1 i.i.d. samples from the GP defined in Theorem 2.1, again yielding the known Chi-squared distribution for its sample variance σ¯f2(x;θ∞i…K+1) σ_f^2(x; _∞^i… K+1). See Appendix B.1.5. Theorem 3.4 establishes a distributional equality between the empirical error of a multi-headed RND architecture and the empirical variance of a finite ensemble of neural networks in the limit of infinite width, providing a theoretical motivation for the use of RND and its common multi-headed architecture as an uncertainty quantification technique. In a broader sense, we believe this analysis is insightful to many practitioners using random network distillation by establishing an intuitive link between theory and practice. Still, the NTK-based perspective applies to an inherently idealized regime and naturally opens up new avenues for investigation. Understanding the relationship between RND networks and deep ensembles at finite width, where feature learning impacts behavior, remains a critical open question beyond the scope of our current framework. Yet, intriguing possibilities also arise within the infinite-width setting itself: Could the properties of the RND target network be deliberately chosen or modified? Exploring different target initializations offers a computationally inexpensive lever to shape the uncertainty signal captured by RND. Indeed, pursuing this very direction, the next section investigates how a specific adaptation of the RND target function allows us to establish a direct correspondence not just with ensemble variance, but with the principled uncertainty quantification provided by Bayesian posterior inference. 4 Equivalence of Random Network Distillation & Bayesian Posteriors Having formulated an equivalence between standard random network distillation and deep ensemble variance, we now proceed to investigate how theoretical connections to the Bayesian inference framework can be established by invoking deliberate changes to the standard random network distillation algorithm, namely by modifying the fixed target function g. Our goal is to show that the RND error signal itself can, under specific conditions, be interpreted as a draw from a centered Bayesian posterior predictive distribution. To this end, we briefly recall Bayesian inference with the classical Gaussian linear model. We define a regression model as f(x;θ)=ϕ(x)⊤θf(x;θ)=φ(x) θ with a feature mapping ϕ:ℝdin→ℝdPφ:R^d_in R^d_P, and a prior distribution over the parameters p(θ)∼(0,Σ0)p(θ) (0, ^0). The prior distribution p(θ)p(θ) implicitly defines a GP prior f0(x;θ)∼(0,ϕ(x)⊤Σ0ϕ(x′))f^0(x;θ) (0,φ(x) ^0φ(x )), with the prior kernel Kxx′=ϕ(x)⊤Σ0ϕ(x′)K_x =φ(x) ^0φ(x ). Within this linear model222We use a noise-free regression model for ease of notation here, but extensions to the noisy case by including an observation noise term σn2I _n^2I in the kernel matrix inversions (cf. Eq. (10)-(11)) are straightforward., we look to infer a posterior distribution over functions given observations =xi∈ℝdini=1NDX=\x_i ^d_in\_i=1^N_D and labels =yi∈ℝi=1NDY=\y_i \_i=1^N_D. Owing to our prior choice, the corresponding posterior predictive distribution conditioned on ,X,Y is a GP with p(f|x,,)∼(KxK−1,Kxx−KxK−1Kx). p(f|x,X,Y) (K_xXK^-1_XXY,\,K_x-K_xXK^-1_XXK_Xx)\,. (10) When contrasting this identity with the GP governing the distribution of converged N functions of Theorem 3.1, one observes a disparity in the structure of the covariance functions. While Theorem 2.1 and Theorem 3.1, too, specify GPs, they do not permit an interpretation as a Bayesian posterior predictive distribution [leeWideNeuralNetworks2020] due to the presence of two (in general) distinct kernel functions, namely the NNGP kernel κ and the NTK Θ . However, inspection of Theorem (3.1) and Eq. (10) suggests a path: if the prior kernel components within Σxx′ϵ(ϑ∞,ψ0) ^ε_x ( _∞, _0), namely κxx′ϵ _x ^ε, are be aligned with the dynamics kernel Θxx′ _x (i.e., if κϵ∝Θκ^ε ), then the resulting covariance structure simplifies to the desired Bayesian posterior form of f(x;θ∞)∼(ΘxΘ−1,Θxx−ΘxΘ−1Θx). f(x; _∞) ( _xX _XX^-1Y,\, _x- _xX ^-1_XX _Xx )\,. (11) An important insight here is that Eq. 11 now is the exact Bayesian posterior predictive distribution of a neural network in the infinite width limit, which corresponds to a kernel regression model with the NTK as a GP prior (0,Θxx′)GP(0, _x ) and conditioned on the data (,)(X,Y). The idea of aligning the prior and dynamic kernels has been previously explored by heBayesianDeepEnsembles2020a to construct Bayesian ensembles where the predictive distribution of the ensemble matches the posterior predictive distribution of the NTK-GP. We propose that a similar alignment can be achieved in the RND framework by constructing the target function g(x;ψ0)g(x; _0) to assume a specific form. The idea is to design a target g~(x;ϑ0,ψ0) g(x; _0, _0) such that when a predictor u(x;ϑ0)u(x; _0) is trained to match it, the resulting “Bayesian” error distribution ϵb(x;ϑ∞,ϑ0,ψ0)=u(x;ϑ∞)−g~(x;ϑ0,ψ0)ε^b(x; _∞, _0, _0)=u(x; _∞)- g(x; _0, _0) behaves like a draw from the posterior of a Bayesian model whose prior kernel is the NTK Θxx′ _x itself 333The newly constructed target function g~(x;ϑ0,ψ0) g(x; _0, _0) uses both ϑ0 _0 and ψ0 _0 for reasons that will become clear in the remainder of section. . In the random network distillation algorithm, the prior kernel κxx′ϵbκ^ε^b_x of initialization errors ϵb(x;ϑ0,ϑ0,ψ0)=u(x;ϑ0)−g~(x;ϑ0,ψ0)ε^b(x; _0, _0, _0)=u(x; _0)- g(x; _0, _0) is given by the sum of the online prior kernel and the target prior kernel κxx′ϵb=κxx′u+κxx′g~κ^ε^b_x =κ^u_x +κ g_x (cf. Theorem 3.1), provided that u and g~ g follow independent GPs. To obtain an error prior kernel that aligns with the NTK such that κxx′ϵb=Θxx′κ^ε^b_x = _x , one may thus construct the target prior such that it satisfies κxx′g~=Θxx′−κxx′uκ g_x = _x -κ^u_x . To this end, a closer inspection of the relation between the NNGP kernel κxx′uκ^u_x and the NTK Θxx′ _x is instructive. For this purpose, we will view the online network u(x;ϑ0)u(x; _0) as a random feature model with its forward computation path as described in Eq. 1. Let in this scenario xL(x)x^L(x) denote the output vector, or the post-activations, before the final linear layer and denote the last-layer parameters at initialization t=0t=0 as (wL,bL)(w^L,b^L). We can write the N output at initialization u(x;ϑ0)u(x; _0) as u(x;ϑ0) u(x; _0) =σbbL+σwnL−1∑i=1nL−1wiLxiL(x), = _bb^L+ _w n_L-1 Σ _ i=1 n_L-1w_i^Lx_i^L(x)\,, (12) that is, as a simple linear model of the random final post-activations xL(x)x^L(x). Viewing the function in Eq. (12) as a random feature model leads to a central insight: since the last-layer weights and biases (wL,bL)(w^L,b^L) are assumed to be initialized i.i.d. from a standard normal (wL,bL)∼(0,I)(w^L,b^L) (0,I), Eq. (12) describes a (random) affine transformation of a Gaussian vector 444To see the correspondence in Eq. 13, first notice that due to the i.i.d. initialization of (wL,bL)(w^L,b^L), any cross-products (e.g., involving elements indexed with i≠ji≠ j) vanish in the expectation [u(x;ϑ0)u(x′;ϑ0)]E[u(x; _0)u(x ; _0)]. The expectation thus becomes [u(x;ϑ0)u(x′;ϑ0)]=w≤L,b≤L[σb2+σw2nL−1∑i=1nL−1xiL(x)xiL(x′)]E[u(x; _0)u(x ; _0)]=E_w^≤ L,b^≤ L[ _b^2+ _w^2n_L-1 _i=1^n_L-1x^L_i(x)x^L_i(x )]. By linearity, the expectation on the r.h.s. can be pulled inside the sum and by symmetry we have that w≤L,b≤L[xiL(x)xiL(x′)]E_w^≤ L,b^≤ L[x^L_i(x)x^L_i(x )] is independent of i, s.t. w≤L,b≤L[σw2nL−1∑i=1nL−1xiL(x)xiL(x′)]=σw2[xiL(x)xiL(x′)]E_w^≤ L,b^≤ L[ _w^2n_L-1 _i=1^n_L-1x^L_i(x)x^L_i(x )]= _w^2E[x^L_i(x)x^L_i(x )]. whose covariance in the limit n→∞n→∞ is quantified by the NNGP kernel κxx′uκ^u_x given by κxx′u=[u(x;ϑ0)u(x′;ϑ0)]=σb2+σw2[xiL(x)xiL(x′)]. κ^u_x =E[u(x; _0)u(x ; _0)]= _b^2+ _w^2E[x^L_i(x)x^L_i(x )]\,. (13) Let us now compare this expression for the the prior kernel κxx′uκ^u_x of the online network with its dynamics kernel Θxx′ _x . In particular, we will split the dynamics kernel Θxx′ _x into a last-layer component Θxx′L=∇wL,bLu(x;ϑ0)⊤∇wL,bLu(x′;ϑ0) _x ^L= _\w^L,b^L\u(x; _0) _\w^L,b^L\u(x ; _0) and a component summarizing all preceding parameters Θxx′≤L−1=∇ϑ≤L−1u(x;ϑ0)⊤∇ϑ≤L−1u(x′;ϑ0) _x ^≤ L-1= _ ^≤ L-1u(x; _0) _ ^≤ L-1u(x ; _0) such that Θxx′=Θxx′L+Θxx′≤L−1 _x = _x ^L+ _x ^≤ L-1. Since u(x;ϑ0)u(x; _0) is linear in the last-layer parameters wL,bL\w^L,b^L\ (cf. Eq. 12), we make the crucial observation that the last-layer NTK component Θxx′L ^L_x equals the NNGP prior kernel Θxx′L=κxx′u _x ^L=κ^u_x 555To see this correspondence, notice that the last-layer gradient inner product ∇wL,bLu(x;ϑ0)⊤∇wL,bLu(x′;ϑ0) _\w^L,b^L\u(x; _0) _\w^L,b^L\u(x ; _0) reduces to the sum σb2+σw2nL−1∑i=1nL−1xiL(x)xiL(x′) _b^2+ _w^2n_L-1 _i=1^n_L-1x^L_i(x)x^L_i(x ), where the r.h.s. sum tends to its expectation in the limit nL−1→∞n_L-1→∞ given that summands are identically distributed (as before by symmetry) and independent (which is shown more rigorously for example in Sec. B.1.4).. This property gives a clear instruction for engineering the prior kernel of the target network: by constructing κxx′g~κ g_x such that κxx′g~=Θxx′≤L−1κ g_x = _x ^≤ L-1 and independently from κxx′uκ^u_x , we obtain an error prior as κxx′ϵb=κxx′g~+κxx′u=Θxx′L+Θxx′≤L−1=Θxx′. κ^ε^b_x =κ g_x +κ^u_x = _x ^L+ _x ^≤ L-1= _x \,. (14) In the following, we will thus aim to construct a target function g~(x;ϑ0,ψ0) g(x; _0, _0) with the desired property κxx′g~=Θxx′≤L−1κ g_x = _x ^≤ L-1, in particular by modeling g~ g as a linear function in the feature space corresponding to gradients in earlier layers. This approach has also previously been explored by heBayesianDeepEnsembles2020a to obtain Bayesian ensembles. Proposition 4.1. (Bayesian RND target function) Under the conditions of Theorem 3.1, let u(x;ϑ0)u(x; _0) and g(x;ψ0)g(x; _0) be neural networks of L layers with parameters ϑ0,ψ0∼(0,I) _0, _0 (0,I) i.i.d. Moreover, let ψ0L=wL,bL _0^L=\w^L,b^L\ denote the last-layer parameters of ψ0 _0 and ψ0≤L−1 _0^≤ L-1 the parameters of all preceding layers. Suppose the target function g~(x;ϑ0,ψ0) g(x; _0, _0) is given by g~(x;ϑ0,ψ0)=∇ϑ0u(x;ϑ0)⊤ψ0∗, g(x; _0, _0)= _ _0u(x; _0) _0^*\,, where ψ0∗=ψ0≤L−1,0dim(ψ0L) _0^*=\ _0^≤ L-1,0_dim( _0^L)\ is a copy of ψ0 _0 with its last-layer weights set to 0. In the infinite width limit n→∞n→∞, g~(x;ϑ0,ψ0) g(x; _0, _0) distributes by construction as g~(x;ϑ0,ψ0)∼(0,κxx′g~) g(x; _0, _0) (0, _x g) where κxx′g~=Θxx′≤L−1 _x g= _x ^≤ L-1. Proof sketch. The function g~(x;ϑ0,ψ0) g(x; _0, _0) is by construction equivalent to a linear function with the (random) feature map ∇ϑ0≤L−1u(x;ϑ0) _ _0^≤ L-1u(x; _0) given by the gradient of parameters in the pre-final layers and with a parameter vector ψ0≤L−1 _0^≤ L-1. Conditioned on ϑ0 _0, the random function g~(x;ϑ0,ψ0) g(x; _0, _0) is thus an affine transformation of the Gaussian vector ψ0≤L−1 _0^≤ L-1 and thus a GP itself, at any width n. Using the central results by jacotNeuralTangentKernel2020 that Θ0,xx′→Θxx′ _0,x → _x as n→∞n→∞ and appealing to the bounded convergence theorem, the limiting distribution of the unconditioned random function g~(x;ϑ0,ψ0) g(x; _0, _0), too, becomes Gaussian with the deterministic covariance Θxx′≤L−1 ^≤ L-1_x . While the specific form of the kernel Θxx′≤L−1=Θxx′−Θxx′L _x ^≤ L-1= _x - _x ^L seems unusual as a standalone prior, it is crucially important in shaping the final error distribution. This is because with the altered “Bayesian” target function g~(x;ϑ0,ψ0) g(x; _0, _0) we can shape the covariance structure of errors at initialization by satisfying Eq. 14, appealing to Theorem (3.1). With the engineered target function g~(x;ϑ0,ψ0) g(x; _0, _0), the learning dynamics of an RND model where the predictor network u(x;ϑt)u(x; _t) learns to mimic g~(;ϑ0,ψ0) g(X; _0, _0) can be shaped in the desired way. Our central statement is that the distribution of the error between the converged predictor u(x;ϑ∞)u(x; _∞) and the target function g~(x;ϑ0,ψ0) g(x; _0, _0) will then no longer reflect the variance of deep ensembles trained with gradient descent, but will instead directly exhibit the statistics of a Bayesian posterior predictive distribution derived from the NTK-GP prior. Theorem 4.2 formalizes this result. Theorem 4.2. (Distribution of Bayesian RND errors) Under the conditions of Theorem 3.1, let u(x;ϑ∞)u(x; _∞) be a converged predictor network trained on data X with labels from the fixed target function g~(;ϑ0,ψ0) g(X; _0, _0) as defined in Proposition 4.1. Let parameters ϑ0,ψ0 _0, _0 be drawn i.i.d. ϑ0,ψ0∼(0,I) _0, _0 (0,I). The convergenced Bayesian RND error ϵb(T;ϑ∞,ϑ0,ψ0)=u(T;ϑ∞)−g~(T;ϑ0,ψ0)ε^b(X_T; _∞, _0, _0)=u(X_T; _∞)- g(X_T; _0, _0) on a test set TX_T is Gaussian with zero mean and covariance ΣTTϵb(ϑ∞,ϑ0,ψ0) ^ε^b_X_TX_T( _∞, _0, _0) =ΘTT−ΘTΘ−1ΘT, = _X_TX_T- _X_TX _XX^-1 _XX_T\,, and thus recovers the covariance of the exact Bayesian posterior predictive distribution of an infinitely wide neural network with the corresponding NTK Θxx′ _x . Proof sketch. The result follows by combining Theorem 3.1 and Proposition 4.1, provided that the GP governing the predictor initialization κxx′uκ^u_x and the target function κxx′g~κ g_x are independent. Owing to the fact that the parameters ϑ0 _0 and ψ0 _0 are drawn independently, the independence between u(x;ϑ0)u(x; _0) and g~(x;ϑ0,ψ0) g(x; _0, _0) is apparent by rewriting the covariance [u(x;ϑ0)g~(x;ϑ0,ψ0)]E[u(x; _0) g(x; _0, _0)] in terms of conditional expectations on ϑ0 _0 by the law of total expectation. Furthermore, since Θxx′=Θxx′L+Θxx′≤L−1 _x = ^L_x + ^≤ L-1_x and κxx′g~=Θxx′≤L−1κ g_x = ^≤ L-1_x , κxx′u=Θxx′Lκ^u_x = ^L_x , we have that κxx′ϵb=Θxx′κ^ε^b_x = _x . In other words, the GP kernel of initial errors aligns with the NTK of the online predictor, such that the distribution of post-convergence errors in Theorem 3.1 simplifies significantly. This same covariance function indeed also defines the posterior predictive distribution of infinitely wide neural networks as described by the GP with prior (0,Θxx′)GP(0, _x ) and conditioned on (,)(X,Y). Theorem 4.2 shows that with a specifically engineered target function, the RND error signal ϵb(x;ϑ∞,ϑ0,ψ0)=u(x;ϑ∞)−g~(x;ϑ0,ψ0)ε^b(x; _∞, _0, _0)=u(x; _∞)- g(x; _0, _0) is no longer just related to ensemble variance, but rather becomes a direct sample from the centered posterior predictive distribution of a Bayesian model whose prior kernel is the NTK itself. This novel result provides a direct bridge between RND and Bayesian inference in the limit of infinite network width, providing a useful insight: the error signal generated by this modified RND procedure is not merely a heuristic measure of distance, but is itself a random draw from the (centered) Bayesian posterior predictive distribution of an NTK-based GP. This direct distributional equivalence has immediate practical implications, for example prescribing rather straightforwardly how this Bayesian form of RND can be used for exact posterior sampling. By applying Proposition 3.3 to the multi-headed Bayesian RND architecture666In a multi-headed architecture, the Bayesian target function described in Proposition 4.2 becomes a JVP. Several common machine learning libraries (e.g., JAX [jax2018github] offer dedicated algorithms to compute such JVPs efficiently., in contrast to obtaining samples from deep ensembles as done in Theorem 3.4, we now obtain several independent samples from the centered posterior predictive distribution through ϵib(x;ϑ∞,ϑ0,ψ0)=ui(x;ϑ∞)−g~i(x;ϑ0,ψ0)ε^b_i(x; _∞, _0, _0)=u_i(x; _∞)- g_i(x; _0, _0). The below corollary details how this can be leveraged to conduct a posterior sampling procedure, requiring access only to a mean estimate and a single Bayesian RND model. Corollary 4.3 (Posterior Sampling via Bayesian RND). Let (μb(x),Σxx′b)N (μ^b(x)\,,\,\, ^b_x ) be the posterior predictive distribution of an infinitely wide neural network conditioned on x with mean μb(x)=ΘxΘ−1μ^b(x)= _xX ^-1_XXY and covariance Σxx′b=Θxx′−ΘxΘ−1Θx′ ^b_x = _x - _xX ^-1_XX _Xx . Suppose μ~(x;θ∞)≈μb(x) μ(x; _∞)≈μ^b(x) is an estimate of the mean function and let ϵib(x;ϑ∞,ϑ0,ψ0)i=1K\ε^b_i(x; _∞, _0, _0)\_i=1^K be error functions of a K-head Bayesian RND model as defined in Theorem 4.2. The following procedure generates (at most K) independent samples from the conditional posterior predictive distribution (μb(x),Σxx′b)N (μ^b(x)\,,\,\, ^b_x ): 1. sample i∼[1,K]i [1,K] 2. compute μ~i(x)=μ~(x;θ∞)+ϵib(x;ϑ∞,ϑ0,ψ0) μ_i(x)= μ(x; _∞)+ε^b_i(x; _∞, _0, _0) 3. μ~i(x) μ_i(x) is an i.i.d. sample from the conditional posterior predictive (μb(x),Σxx′b)N (μ^b(x)\,,\,\, ^b_x ) Proof sketch. The result follows directly from Theorem (4.2) and application of the independence argument of Proposition (3.3) to the multi-headed setting. Corollary 4.3 shows that, given an estimator of the posterior predictive mean, a modified Bayesian RND setup can be used to perform direct Bayesian posterior sampling in the NTK limit. By extension, this offers a pathway to performing exact Bayesian inference through the lens of network distillation, provided that the target and predictor networks initializations are handled deliberately. This completes our theoretical development, first showing an equivalence of RND in the NTK regime to ensemble variance and now, through specific modifications to its target function, to the generation of independent samples from exact Bayesian posterior predictive distributions. 5 Numerical Analysis We proceed with a numerical analysis to validate the thus far presented results. In the following, we study how predictive RND errors relate to predictive variances of deep ensembles in practice, both in the standard and Bayesian settings. To this end, we train two-layer connected neural networks with SiLU activations [elfwing2018sigmoid] on a synthetic dataset with N=10N=10 train and N~=5000 N=5000 test samples from an isotropic Gaussian xi∼(0,I3)x_i (0,I_3). Ensemble models are fit to a toy target function, and multiheaded RND models optimized as described above. The variance of the true underlying GP is approximated with Monte-Carlo estimates of 512 independent models and a single Bayesian RND model with 512 heads, such that a small residual amount of discrepancy is to be expected. Fig. 1 shows a stark decrease in average squared discrepancy between test evaluations of predictive ensemble variances and RND errors as model width increases, a trend in line with our theoretical derivations and present even at practical network widths. Further evaluations and details of this experiment are reported in Appendix C. Figure 1: Test-set errors between predictive variances of (Bayesian) ensembles and self-predictive errors of (Bayesian) RND vanish with large layer widths. 6 Related Work A substantial body of research studies the analytical learning dynamics of deep learning, particularly in the infinite-width limit. Central to our analysis are seminal works characterizing the NNGP [lee2017deep] at initialization, the dynamics-governing NTK [jacotNeuralTangentKernel2020], and the evolution of wide networks as linear models [leeWideNeuralNetworks2020, arora2019exact, chizat2018global]. This provides a theoretical framework for analytical descriptions of deep ensembles [lakshminarayananSimpleScalablePredictive2017, dietterich2000ensemble], with subsequent studies using NTK theory to precisely characterize ensemble variances under various conditions, including observation noise [yang2019scaling, kobayashi2022disentangling, calvo2024epistemic]. A central line of work for our paper is the connection between deep ensembles and Bayesian inference in infinite-width NTK regime. Notably, heBayesianDeepEnsembles2020a demonstrate how to construct “Bayesian ensembles”, an approach we adopt to construct “Bayesian RND” algorithms. The broader link between deep ensembles and approximations of Bayesian posteriors has been studied extensively [khan2019approximate, osawa2019practical, dangelo2021repulsive, osbandDeepExplorationRandomized2019, izmailov2021bayesian]. More recently, NTK-based approaches have been used for single-model uncertainty estimation [zanger2026contextual] or ad-hoc uncertainty quantification [wilson2025uncertainty]. Our work provides a theoretical basis for RND [burda2018exploration], which belongs to a class of computationally cheaper, single-model methods [pathakCuriosityDrivenExplorationSelfSupervised2017, lahlou2021deup, guo2022byol, sensoyEvidentialDeepLearning2018, van2020uncertainty, rudner2022tractable, laurent2022packed, tagasovska2019single]. Moreover, Uncertainty quantification from the lens of learning dynamics is moreover widespread in reinforcement learning (RL)[xiao2021understanding, cai2019neural, wai2020provably, lyle2022learning, yang2020provably], the original application domain of RND. Notably, zanger2026universal derive an RND-like estimator for value function uncertainty using NTK theory. More broadly, deep ensembles and Bayesian methods are widely used in RL, driving exploration [osbandDeepExplorationBootstrapped2016b, chen2017ucb, osbandDeepExplorationRandomized2019, nikolovInformationDirectedExplorationDeep2019, ishfaq2021randomized, zanger2024diverse]. 7 Conclusions In this work, we have established a novel theoretical understanding of random network distillation (RND) by connecting it to the principled uncertainty frameworks of deep ensembles and Bayesian inference. By analyzing these techniques within the unifying setting of infinitely wide neural networks, we provide a clear analytical interpretation for the empirically successful RND algorithm. Our analysis yields a twofold equivalence: first, we prove that the squared error of standard RND exactly recovers the predictive variance of deep ensembles in the NTK regime. Second, we demonstrate that the RND framework is more versatile; by deliberately designing the RND target function, the resulting error signal can be made to directly mirror the centered posterior predictive distribution of an NTK-governed GP, that is, the exact posterior predictive distribution of neural networks in the infinite width limit. This “Bayesian RND” variant furthermore allows for posterior sampling procedures that produce i.i.d. samples from this posterior. Our work thereby unifies RND, ensembles, and Bayesian inference under the same theoretical lens from an infinite width perspective. Crucially, our findings hold under the assumptions infinite-width and the NTK regime, a setting where networks effectively linearize and operate as kernel machines with a fixed kernel. This “lazy” training regime, while analytically tractable and predictive for very wide networks, does not capture the phenomenon of feature learning. The degree to which our established equivalences translate to practical, finite-width networks that learn features remains a significant open question. Conversely, this also suggest avenues for future research: deviations between RND, ensembles, and Bayesian posteriors in practice must arise from departures from the NTK regime. Characterizing specifically these deviations could lead to novel techniques and a deeper understanding of computationally efficient approaches that approximate Bayesian inference, operating well outside the kernelized infinite-width setting. Another exciting direction is the concept of target engineering as cheap way of studying priors for Bayesian deep learning, an actively studied field that garners widespread interested from the uncertainty quantification and Bayesian deep learning community. References On the Equivalence of Random Network Distillation, Deep Ensembles, and Bayesian Inference (Supplementary Material) Appendix A Limitations and Assumptions We provide an overview of the primary assumptions underpinning our analysis and discuss their relation to practical settings. The foremost assumption is that our analysis operates within the NTK regime. This framework presupposes the asymptotic limit of infinitely wide neural networks and a so-called NTK-parametrization of forward computations that ensures network dynamics linearize around their initialization, leading to “lazy” learning with kernel regression behavior. This idealized setting naturally deviates from practical implementations involving finite-width networks. Nonetheless, a significant body of work has demonstrated that predictions from NTK theory can remain remarkably accurate for sufficiently wide, modern architectures, providing a reasonable approximation of their behavior [e.g., lee2020finite, seleznova2022analyzing, samarin2020empirical]. Furthermore, our derivations assume training via full-batch gradient flow, which corresponds to gradient descent with an infinitesimal step size. This abstains from the use of stochastic minibatch optimizers, which are standard in practice. While beyond our current scope, extensions of NTK analysis to incorporate the effects of stochastic gradient noise do exist [e.g., yang2019scaling, cao2019generalization, nitanda2021optimal]. Finally, our analysis considers a fixed training dataset X. This contrasts with prominent applications of RND, particularly in online reinforcement learning, where the agent interacts with an environment and learns from an inherently non-stationary data stream. Characterizing how these equivalences with ensembles and Bayesian posteriors evolve under such distribution shifts remains an important open question. Appendix B Proofs This section provides extended proofs for our analysis of RND. B.1 Ensemble Equivalence Our first result states the equivalence of self-predictive errors of RND and predictive variance of deep ensembles in the infinite-width NTK regime. For completeness, we also include proofs or simplified proof sketches for known results that support our analysis. Theorem B.1. [jacotNeuralTangentKernel2020](Post-convergence neural network function) In the limit of infinite layer widths n→∞n ∞ and infinite time t→∞t ∞, the output function of a neural network f(x;θ∞)f(x; _∞) with NTK parametrization according to Eq. 1 is given by f(x;θ∞) f(x; _∞) =f(x;θ0)−ΘxΘ−1(−f(;θ0)), =f(x; _0)- _xX _XX^-1 (Y-f(X; _0) )\,, where we used the shorthand Θxx′≡Θ(x,x′) _x ≡ (x,x ). Proof sketch. By taking the infinite width limit n→∞n→∞, we obtain a linear ODE from Eq. (2). Through an exponential ansatz, its explicit solution with initial condition f(x;θ0)f(x; _0) is given by f(x;θt)=f(x;θ0)+ΘxΘ−1(I−e−tΘ)(−f(;θ0)).f(x; _t)=f(x; _0)+ _xX _XX^-1(I-e^-t _XX)(Y-f(X; _0)). Assuming the training Gram matrix Θ _XX is positive definite (and thus invertible), the exponential term decays to zero as t→∞t→∞, yielding the kernel regression formula in Proposition (B.1). See jacotNeuralTangentKernel2020 and Appendix B.1.1. B.1.1 Proof of Theorem B.1 Proof. The proof is centered around the learning dynamics of a neural network under gradient descent, whereby we assume the limit of infinitesimal step size for simplicity. This setting is also referred to as “gradient flow”. The driving force behind the learning dynamics of parameters θt _t is gradient flow optimization on the loss ℒ(θt) ( _t) =12‖f(,θt)−‖22, = 12\|\,f(X, _t)-Y\,\|^2_2, (15) with the subsequent evolution of parameters by dtθt ddt _t =−α∇θℒ(θt), =-α _θL( _t)\,, (16) where α is a learning rate. From this, we can obtain the parameter space differential equation dtθt ddt _t =−α∇θf(,θt)(f(,θt)−). =-α _θf(X, _t) (f(X, _t)-Y\ )\,. (17) In order to translate this expression to a function-space view through a first-order Taylor expansion of f around its initialization parameters θ0 _0: flin(x,θt)=f(x,θ0)+∇θf(x,θ0)⊤(θt−θ0). f_lin(x, _t)=f(x, _0)+ _θf(x, _0) ( _t- _0)\,. (18) The use of a linearized neural network function simplifies the analysis in two aspects: 1.) the linearization offers a simple translation of the parameter space evolution dtθt ddt _t to a function-space evolution and 2.) the linearized neural network function flin(x,θt)f_lin(x, _t) results in linear dynamics, simplifying the earlier derived differential equation to a linear ODE. The evolution of flinf_lin is then obtained by taking the time-derivative of Eq. (18) and plugging in the parameter evolution for a linearized function from Eq. (17) such that dtflin(x,θt) ddtf_lin(x, _t) =−α∇θf(x,θ0)⊤∇θf(,θ0)(flin(,θt)−). =-α _θf(x, _0) _θf(X, _0) (f_lin(X, _t)-Y )\,. (19) Let us denote the training error of flinf_lin at time t with δt=flin(,θt)− _t=f_lin(X, _t)-Y and accordingly write dtδt ddt _t =−αΘ0δt, =-α ^0_XX _t\,, (20) where Θ0 ^0_XX denotes the empirical tangent kernel Θ0=∇θf(,θ0)⊤∇θf(,θ0) ^0_XX= _θf(X, _0) _θf(X, _0) at initialization. The differential equation (20) is a linear ODE system to which an exponential ansatz provides the explicit solution δt=e−αtΘ0δ0, _t=e^-α t ^0_XX _0\,, (21) where eΘ=∑k=0∞1k!(Θ)Ke _XX= _k=0^∞ 1k!( _XX)^K is the matrix exponential. We plug this result back in the linearized function space differential equation 19 to obtain dtflin(x,θt) ddtf_lin(x, _t) =−αΘx0e−αtΘ0(f(,θ0)−). =-α _xX^0e^-α t ^0_XX (f(X, _0)-Y )\,. (22) In this form, we can solve for flin(x,θt)f_lin(x, _t) directly by integration flin(x,θt) f_lin(x, _t) =f(x,θ0)+∫0tddt′flin(x,θt′)dt′ =f(x, _0)+ _0^t ddt f_lin(x, _t ) \!dt (23) =f(x,θ0)+Θx0(Θ0)−1(e−αtΘ0−I)(f(,θ0)−). =f(x, _0)+ ^0_xX( ^0_XX)^-1 (e^-α t ^0_XX-I ) (f(X, _0)-Y )\,. (24) Remarkably, the linearized and true learning dynamics become increasingly aligned with increasing neural network width. jacotNeuralTangentKernel2020 and leeWideNeuralNetworks2020 show that as network width increases, the required individual movement of parameters θt−θ0 _t- _0 to effect sufficient movement in the output function f(x,θt)f(x, _t) decreases. In the limit of infinite width n→∞n→∞, the linearization of f then becomes exact limn→∞flin(x,θt)=f(x,θt) _n→∞f_lin(x, _t)=f(x, _t). Under the outlined training dynamics, the same limit furthermore causes the NTK to become deterministic (despite random weight initializations) and stationary limn→∞Θxx′0=Θxx′t=Θxx′ _n→∞ _x ^0= _x ^t= _x . Thus, the convergenced function at time t→∞t→∞ is described by f(x,θ∞) f(x, _∞) =f(x,θ0)−ΘxΘ−1(f(,θ0)−). =f(x, _0)- _xX _XX^-1 (f(X, _0)-Y )\,. (25) ∎ B.1.2 Proof of Theorem 2.1 We restate Theorem 2.1 for convenience. See 2.1 Proof sketch. We use the fact that f(x;θ∞)f(x; _∞) can be written as a linear combination of the test initialization f(x;θ0)f(x; _0) and the training initialization f(;θ0)f(X; _0). Both these identities are described probabilistically by the NNGP f(x;θ0)∼(0,κxx′)f(x; _0) (0, _x ), and f(;θ0)∼(0,κ)f(X; _0) (0, _XX). Applying a linear transformation to a GP yields another GP [rasmussen2006gp], meaning f(x;θ∞)f(x; _∞) also follows a GP. Propagating the prior covariance κ through the linear transformation described by Proposition B.1 reveals the expression for the post-convergence covariance function ΣTTf(θ∞) ^f_X_TX_T( _∞) given in Theorem 2.1. Proof. The proof builds on the previous result of Proposition B.1 providing a closed-form expression for the post-convergence function as a deterministic function of its initialization, here evaluated for a set of test points TX_T f(T,θ∞) f(X_T, _∞) =f(T,θ0)−ΘTΘ−1(f(,θ0)−). =f(X_T, _0)- _X_TX _XX^-1 (f(X, _0)-Y )\,. (26) To be precise, the post-convergence predictions f(T,θ∞)f(X_T, _∞) can be written as an affine transformation of the vector (f(T,θ0),f(,θ0)⊤)⊤(f(X_T, _0),f(X, _0) ) . This yields the block matrix equation (f(T,θ∞)f(,θ∞))= pmatrixf(X_T, _∞)\\ f(X, _∞) pmatrix= (I−Θ(T,)Θ(,)−100)(f(T,θ0)f(,θ0))+(Θ(T,)Θ(,)−1). pmatrixI&- (X_T,X) (X,X)^-1\\ 0&0 pmatrix pmatrixf(X_T, _0)\\ f(X, _0) pmatrix+ pmatrix (X_T,X) (X,X)^-1Y\\ Y pmatrix\,. (27) We recall that, at initialization, neural networks in the infinite width limit distribute to a GP called NNGP [lee2017deep] as f(T,θ0)∼(0,κTT)whereκTT=θ0[f(T,θ0)f(T,θ0)⊤]. f(X_T, _0) (0, _X_TX_T)\, _X_TX_T=E_ _0[f(X_T, _0)f(X_T, _0) ]\,. (28) The block eq. (B.1.2) thus describes an affine transformation of a GP itself. We have that affine transformations of multivariate Gaussian random variables X∼(μX,ΣX)X ( _X, _X) with Y=a+BXY=a+BX distribute Gaussian themselves with Y∼(a+BμX,BΣXB⊤)Y (a+B _X,\,\,B _XB ). Application to Eq. B.1.2 and rearrangement then yields the post-convergence GP with mean and covariance [f(T,θ∞)]=ΘTΘ−1, [f(X_T, _∞)]= _X_TX _XX^-1Y\,, (29) ΣTTf(θ∞)= ^f_X_TX_T( _∞)= κTT+ΘTΘ−1κΘ−1ΘT−(ΘTΘ−1κT+h.c.), _X_TX_T+ _X_TX _XX^-1 _XX _XX^-1 _XX_T- ( _X_TX _XX^-1 _XX_T+h.c. )\,, (30) where h.c. refers to the Hermitian conjugate of the preceding term. This completes the proof. ∎ B.1.3 Proof of Theorem 3.1 We restate Theorem 3.1 for convenience. See 3.1 Proof. This proposition considers the post-convergence distribution of self-predictive errors as produced by RND. The online predictor u(x;ϑt)u(x; _t) undergoes learning dynamics under the same conditions as outlined in the derivation of Proposition B.1, albeit with the self-predictive loss ℒ(ϑt) ( _t) =12‖u(,ϑt)−g(,ψ0)‖22. = 12\|\,u(X, _t)-g(X, _0)\,\|^2_2\,. (31) This, by analogy to Theorem B.1, implies that the online predictor u(x;ϑt)u(x; _t) converges as t→∞t→∞ to the function u(x,ϑ∞) u(x, _∞) =u(x,ϑ0)−ΘxΘ−1(u(,ϑ0)−g(,ψ0)). =u(x, _0)- _xX _XX^-1 (u(X, _0)-g(X, _0) )\,. (32) For a set of test points TX_T, the error ϵ(T;ϑ∞,ψ0)=u(T;ϑ∞)−g(T;ψ0)ε(X_T; _∞, _0)=u(X_T; _∞)-g(X_T; _0) at convergence can thus be written as the affine transformation ϵ(T;ϑ∞,ψ0) ε(X_T; _∞, _0) =ϵ(T;ϑ0,ψ0)−ΘTΘ−1ϵ(;ϑ0,ψ0). =ε(X_T; _0, _0)- _X_TX _XX^-1ε(X; _0, _0)\,. (33) and the corresponding block matrix equation (ϵ(T;ϑ∞,ψ0)ϵ(;ϑ∞,ψ0)) pmatrixε(X_T; _∞, _0)\\ ε(X; _∞, _0) pmatrix =(I−ΘTΘ−100)(ϵ(T;ϑ0,ψ0)ϵ(;ϑ0,ψ0)). = pmatrixI&- _X_TX _XX^-1\\ 0&0 pmatrix pmatrixε(X_T; _0, _0)\\ ε(X; _0, _0) pmatrix\,. (34) The errors accordingly are themselves Gaussian with ϵ(T;ϑ∞,ψ0)∼(0,κTTϵ)ε(X_T; _∞, _0) (0,κ^ε_X_TX_T) where κTTϵ=ϑ0,ψ0[ϵ(T;ϑ0,ψ0)ϵ(T;ϑ0,ψ0)⊤]κ^ε_X_TX_T=E_ _0, _0[ε(X_T; _0, _0)ε(X_T; _0, _0) ]. The latter term describes the distribution of self-predictive errors at initialization, which is a simple sum of two independent NNGP ϵ(T;ϑ0,ψ0)=u(T;ϑ0)−g(T;ψ0)ε(X_T; _0, _0)=u(X_T; _0)-g(X_T; _0) such that κTTϵ=κTTu+κTTgκ^ε_X_TX_T=κ^u_X_TX_T+κ^g_X_TX_T, completing the proof. ∎ B.1.4 Proof of Proposition 3.3 Before treating Proposition 3.3 we first derive two known results concerning the independence and recursive character of the NNGP kernel and the NTK. We assume forward computations of f(x;θt)f(x; _t) are defined according to Eq. 1. To avoid confusion with indices i,ji,j we will in this section use the notation κ(x,x′)κ(x,x ) rather than κxx′ _x to denote the function inputs x,x′x,x (and similarly for Θ(x,x′) (x,x )). Proposition B.2. [lee2017deep] (Recursive NNGP formulation) At initialization t=0t=0 and in the limit n→∞n→∞, the i-th output at layer l, zil(x;θ0≤l)z_i^l(x; _0^≤ l), converges to a GP with zero mean and covariance function κiil(x,x′) _i^l(x,x ) given by κii1(x,x′) _i^1(x,x ) =σw2n0x⊤x′+σb2,andkij1(x,x′)=0,if i≠j, = _w^2n_0x x + _b^2, k_ij^1(x,x )=0, i≠ j\,, (35) κiil(x,x′) _i^l(x,x ) =σb2+σw2zil−1∼(0,κiil−1)[ϕ(zil−1(x;θ0≤l−1))ϕ(zil−1(x′;θ0≤l−1))], = _b^2+ _w^2E_z_i^l-1 (0, _i^l-1)[φ(z_i^l-1(x; _0^≤ l-1))φ(z_i^l-1(x ; _0^≤ l-1))]\,, (36) andκijl(x,x′)=0,if i≠j, \,\,\,and _ij^l(x,x )=0, i≠ j\,, (37) and we have κiil(x,x′)=κl(x,x′),∀i _i^l(x,x )=κ^l(x,x )\,, ∀ i. Proof. We prove the proposition by induction. The induction assumption is that if outputs at layer l−1l-1 satisfy a GP structure zil−1∼(0,κl−1), z_i^l-1 (0,κ^l-1), (38) with the covariance function defined as κijl−1(x,x′) _ij^l-1(x,x ) =[zil−1(x;θ0≤l−1)zjl−1(x′;θ0≤l−1)]=kl−1(x,x′)if i=j,0if i≠j, =E[z_i^l-1(x; _0^≤ l-1)z_j^l-1(x ; _0^≤ l-1)]= casesk^l-1(x,x )& i=j\,,\\ 0& i≠ j\,,\\ cases (39) then, outputs at layer l follow zil(x)∼(0,κl), z_i^l(x) (0,κ^l), (40) where the NNGP kernel at layer l is given by: κiil(x,x′) _i^l(x,x ) =[zil(x;θ0≤l)zil(x′;θ0≤l)]=κl(x,x′),∀i, =E[z_i^l(x; _0^≤ l)z_i^l(x ; _0^≤ l)]=κ^l(x,x ), ∀ i, (41) κijl(x,x′) _ij^l(x,x ) =[zil(x;θ0≤l)zjl(x′;θ0≤l)]=0,if i≠j. =E[z_i^l(x; _0^≤ l)z_j^l(x ; _0^≤ l)]=0, i≠ j. (42) with the recursive definition κl(x,x′)=σb2+σw2zil−1∼(0,kl−1)[ϕ(zil−1(x;θ0≤l−1))ϕ(zil−1(x′;θ0≤l−1))]. κ^l(x,x )= _b^2+ _w^2E_z_i^l-1 (0,k^l-1)[φ(z_i^l-1(x; _0^≤ l-1))φ(z_i^l-1(x ; _0^≤ l-1))]. (43) Base case (l=1)(l=1). At layer l=1l=1 we have: zi1(x;θ0≤1)=σwn0∑j=1n0wij1xj+σbbi1. z_i^1(x; _0^≤ 1)= _w n_0 _j=1^n_0w_ij^1x_j+ _bb_i^1\,. (44) This is an affine transform of Gaussian random variables; thus, zi1(x;θ0≤1)z_i^1(x; _0^≤ 1) distributes Gaussian with zi1(x)∼(0,κ1), z_i^1(x) (0,κ^1), (45) with kernel κ1(x,x′)=σw2n0x⊤x′+σb2=κii1(x,x′),andκij1=0,if i≠j, κ^1(x,x )= _w^2n_0x x + _b^2=κ^1_i(x,x )\,, _ij^1=0, \,i≠ j\,, (46) where the independence follows from the fact that zi1(x;θ0≤1)z_i^1(x; _0^≤ 1) is computed from separate, independent rows of weights and biases. Induction step l>1l>1. For layers l>1l>1 we have zil(x;θ0≤l)=σbbil+σwnl−1∑j=1nl−1wijlxjl(x),xjl(x)=ϕ(zjl−1(x;θ0≤l−1)). z_i^l(x; _0^≤ l)= _bb_i^l+ _w n_l-1 _j=1^n_l-1w_ij^lx_j^l(x), x_j^l(x)=φ(z_j^l-1(x; _0^≤ l-1))\,. (47) By the induction assumption, zjl−1(x;θ0≤l−1)z_j^l-1(x; _0^≤ l-1) are generated by independent GP. Hence, xil(x)x_i^l(x) and xjl(x)x_j^l(x) are independent for i≠ji≠ j. Consequently, zil(x;θ0≤l)z_i^l(x; _0^≤ l) is a sum of independent random variables. By the CLT (as n1,…,nL→∞n_1,…,n_L→∞) the tuple zil(x;θ0≤l),zil(x′;θ0≤l)\z_i^l(x; _0^≤ l),z_i^l(x ; _0^≤ l)\ tends to be jointly Gaussian, with covariance given by: [zil(x;θ0≤l)zil(x′;θ0≤l)]= [z_i^l(x; _0^≤ l)z_i^l(x ; _0^≤ l)]= σb2+σw2zil−1∼(0,κl−1)[ϕ(zil−1(x;θ0≤l−1))ϕ(zil−1(x′;θ0≤l−1))]. _b^2+ _w^2E_z_i^l-1 (0,κ^l-1)[φ(z_i^l-1(x; _0^≤ l-1))φ(z_i^l-1(x ; _0^≤ l-1))]\,. (48) Moreover, as zilz_i^l and zjlz_j^l for i≠ji≠ j are defined through independent rows of the parameters wl,blw^l,b^l and independent pre-activations xl(x)x^l(x), we have κijl=[zil(x)zjl(x′)]=0,if i≠j, _ij^l=E[z_i^l(x)z_j^l(x )]=0, \,i≠ j, (49) and thus completing the proof. ∎ Proposition B.3. [jacotNeuralTangentKernel2020] (Recursive NTK formulation) In the limit n→∞n→∞, the neural tangent kernel Θiil(x,x′) ^l_i(x,x ) of the i-th output zil(x;θ0≤l)z_i^l(x; _0^≤ l) at layer l, defined as the gradient inner product Θiil(x,x′)=∇θlzil(x;θ0≤l)⊤∇θlzil(x′;θ0≤l), ^l_i(x,x )= _θ^lz_i^l(x; _0^≤ l) _θ^lz_i^l(x ; _0^≤ l)\,, (50) is given recursively by Θii1(x,x′) _i^1(x,x ) =κii1(x,x′)=σw2n0x⊤x′+σb2,andΘij1(x,x′)=0,if i≠j, = _i^1(x,x )= _w^2n_0x x + _b^2, _ij^1(x,x )=0, \,i≠ j\,, (51) Θiil(x,x′) _i^l(x,x ) =Θiil−1(x,x′)κ˙iil−1(x,x′)+κiil(x,x′), = _i^l-1(x,x ) κ_i^l-1(x,x )+ _i^l(x,x ), (52) where κ˙iil(x,x′) κ_i^l(x,x ) =σw2zil−1∼(0,κiil−1)[ϕ˙(zil−1(x;θ0≤l−1))ϕ˙(zil−1(x′;θ0≤l−1))], = _w^2E_z_i^l-1 (0, _i^l-1)[ φ(z_i^l-1(x; _0^≤ l-1)) φ(z_i^l-1(x ; _0^≤ l-1))]\,, (54) and Θijl(x,x′) _ij^l(x,x ) =∇θlzil(x;θ0≤l)⊤∇θlzjl(x′;θ0≤l)=0if i≠j. = _θ^lz_i^l(x; _0^≤ l) _θ^lz_j^l(x ; _0^≤ l)=0 \,i≠ j. (55) Proof. The proof is by induction. The induction assumption is that if gradients satisfy at layer l−1l-1 Θijl−1(x,x′)= _ij^l-1(x,x )= ∇θl−1zil−1(x;θ0≤l−1)⊤∇θl−1zjl−1(x′;θ0≤l−1)=Θl−1(x,x′)if i=j,0if i≠j, _θ^l-1z_i^l-1(x; _0^≤ l-1) _θ^l-1z_j^l-1(x ; _0^≤ l-1)= cases ^l-1(x,x )& i=j,\\ 0& i≠ j, cases (56) then at layer l we have Θijl(x,x′) _ij^l(x,x ) =Θiil−1(x,x′)κ˙iil(x,x′)+κiil(x,x′)if i=j,0if i≠j. = cases _i^l-1(x,x ) κ_i^l(x,x )+ _i^l(x,x )& i=j\,,\\ 0& i≠ j\,.\\ cases (57) Base case (l=1l=1). At layer l=1l=1, we have zi1(x;θ0≤1) z_i^1(x; _0^≤ 1) =σbbi1+σwn0∑jn0wij1xj, = _bb_i^1+ _w n_0 _j^n_0w_ij^1x_j, (58) and the gradient inner product is given by: ∇θ1zi1(x;θ0≤1)⊤∇θ1zi1(x′;θ0≤1)=σw2n0x⊤x′+σb2=κii1(x,x′). _θ^1z_i^1(x; _0^≤ 1) _θ^1z_i^1(x ; _0^≤ 1)= _w^2n_0x x + _b^2= _i^1(x,x ). (59) Inductive step (l>1l>1). For layers l>1l>1, we split parameters θl=θl−1∪wl,blθ^l=θ^l-1∪\w^l,b^l\ and split the inner product by Θiil(x,x′) _i^l(x,x ) =∇θl−1zil(x;θ0≤l)⊤∇θl−1zil(x′;θ0≤l)⏟l.h.s+∇wl,blzil(x;θ0≤l)⊤∇wl,blzil(x;θ0≤l)⏟r.h.s. = _θ^l-1z_i^l(x; _0^≤ l) _θ^l-1z_i^l(x ; _0^≤ l)_l.h.s+ _\w^l,b^l\z_i^l(x; _0^≤ l) _\w^l,b^l\z_i^l(x; _0^≤ l)_r.h.s. (60) Note that the above r.h.sr.h.s involves gradients w.r.t. last-layer parameters, i.e. the post-activation outputs of the previous layer, and by the same arguments as in the NNGP derivation of Proposition B.2, this is a sum of independent post activations s.t. in the limit nl−1→∞n_l-1 ∞ ∇wl,blzil(x;θ0≤l)⊤∇wl,blzjl(x′;θ0≤l) _\w^l,b^l\z_i^l(x; _0^≤ l) _\w^l,b^l\z_j^l(x ; _0^≤ l) =kiil(x,x′),i=j,0,i≠j. = casesk_i^l(x,x ),& i=j,\\ 0,& i≠ j. cases (61) For the l.h.s.l.h.s., we first apply chain rule to obtain ∇θl−1zil(x;θ0≤l)=σwnl−1∑jnl−1wijlϕ˙(zjl−1(x;θ0≤l−1))∇θl−1zjl−1(x;θ0≤l−1). _θ^l-1z_i^l(x; _0^≤ l)= _w n_l-1 _j^n_l-1w_ij^l φ(z_j^l-1(x; _0^≤ l-1)) _θ^l-1z_j^l-1(x; _0^≤ l-1)\,. (62) The gradient inner product of outputs i and j thus reduces to ∇θl−1zil(x;θ0≤l)⊤∇θl−1zjl(x′;θ0≤l)= _θ^l-1z_i^l(x; _0^≤ l) _θ^l-1z_j^l(x ; _0^≤ l)= σw2nl−1∑knl−1wiklwjklϕ˙(zkl−1(x;θ0≤l−1))ϕ˙(zkl−1(x′;θ0≤l−1))Θkkl−1(x,x′). _w^2n_l-1 _k^n_l-1w_ik^lw_jk^l φ(z_k^l-1(x; _0^≤ l-1)) φ(z_k^l-1(x ; _0^≤ l-1)) _k^l-1(x,x )\,. (63) By the induction assumption Θkkl−1(x,x′)=Θl−1(x,x′) _k^l-1(x,x )= ^l-1(x,x ) and again by the independence of the rows wilw^l_i and wjlw_j^l for i≠ji≠ j, the above expression converges in the limit nl−1→∞n_l-1 ∞ to an expectation with Θijl(x,x′)=Θl−1(x,x′)κ˙iil(x,x′)+κiil(x,x′)i=j,0i≠j, _ij^l(x,x )= cases ^l-1(x,x ) κ_i^l(x,x )+ _i^l(x,x )& i=j,\\ 0& i≠ j\,, cases (64) thereby completing the proof. ∎ We now restate Proposition 3.3 for convenience. See 3.3 Proof. We begin by deriving the training dynamics for the output fi(x;θt)f_i(x; _t) analogously to the proof of Proposition B.1. We denote by iY_i the labels used to train the function fi(x;θt)f_i(x; _t). By Proposition B.3, the training dynamics of fi(x;θt)f_i(x; _t) and fj(x;θt)f_j(x; _t) are decoupled for i≠ji≠ j and we can thus derive Eq. 23 analogously for individual output heads i. Taking the infinite width limit, we obtain at time t fi(x;θt)=fi(x;θ0)+Θii(x,)Θii(,)−1(e−αtΘii(,)−I)(fi(;θ0)−i). f_i(x; _t)=f_i(x; _0)+ _i(x,X) _i(X,X)^-1 (e^-α t _i(X,X)-I )(f_i(X; _0)-Y_i)\,. (65) Thus, the output head fi(x;θt)f_i(x; _t) at time t is a deterministic function of its own initialization only, which itself is characterized by a GP fi(x;θ0)∼(0,κii(x,x′))f_i(x; _0) (0, _i(x,x )) that is independent of output heads j≠ij≠ i by Proposition B.2. And thus, since fi(x;θt)f_i(x; _t) is an affine transform of its own independent initialization terms fi(x;θ0)f_i(x; _0) and fi(;θ0)f_i(X; _0), it too must follow an independent GP with θ0[fi(x;θt)fi(x′;θt)]=Σ(x,x′;θt)E_ _0[f_i(x; _t)f_i(x ; _t)]= (x,x ; _t) and in particular θ0[fi(x;θt)fj(x′;θt)]=0E_ _0[f_i(x; _t)f_j(x ; _t)]=0 if i≠ji≠ j. ∎ B.1.5 Proof of Theorem 3.4 We restate Theorem 3.4 for convenience. See 3.4 Proof. The proof follows by combining the results of Propositions (3.1) and (3.3). We define a multiheaded RND predictor with K output heads ui(x,ϑt)i=1K\u_i(x, _t)\_i=1^K and a fixed multiheaded target network gi(xt;ψ0)i=1K\g_i(x_t; _0)\_i=1^K of equivalent architecture as uiu_i (i.e., both corresponding to the same NTK Θ ) with the corresponding prediction errors ϵi(x;ϑt,ψ0)i=1K\ _i(x; _t, _0)\_i=1^K accordingly. Let ui(x,ϑt)u_i(x, _t) be trained such that each head i is trained to match the i-th target output gi(x;ψ0)g_i(x; _0). By Proposition 3.3, the predictions of online predictor heads ui(x,ϑt)i=1K\u_i(x, _t)\_i=1^K at time t and fixed target networks gi(xt;ψ0)i=1K\g_i(x_t; _0)\_i=1^K are each mutually independent with ϑ0[ui(x;ϑt)uj(x;ϑt)] _ _0[u_i(x; _t)u_j(x; _t)] =0,if i≠j, =0\,, \,i≠ j\,, (66) and ψ0[gi(x;ψ0)gj(x;ψ0)] _ _0[g_i(x; _0)g_j(x; _0)] =0,if i≠j. =0\,, \,i≠ j\,. (67) As a consequence, we also have that ϑ0,ψ0[ϵi(x;ϑt,ψ0)ϵj(x;ϑt,ψ0)]=0,if i≠j. _ _0, _0[ _i(x; _t, _0) _j(x; _t, _0)]=0\,, \,i≠ j\,. (68) As previously established in the proof of Proposition 3.3, the multi-headed functions ϵi(x;ϑt,ψ0)i=1K\ _i(x; _t, _0)\_i=1^K follow equivalent learning dynamics as their scalar-output counterparts. The post-convergence distribution of individual heads ϵi(x;ϑ∞,ψ0) _i(x; _∞, _0) must therefore equal the scalar-output post-convergence distribution established in Theorem 3.1. Consequently, the errors ϵi(x;ϑt,ψ0)i=1K\ _i(x; _t, _0)\_i=1^K are independent and identically distributed draws from a Gaussian with mean and covariance [ϵ(x,ϑ∞,ψ0)] [ε(x, _∞, _0)] =0, =0\,, Σxx′ϵ(ϑ∞,ψ0) ^ε_x ( _∞, _0) =κxx′ϵ+ΘxΘ−1κϵΘ−1Θx′−(ΘxΘ−1κx′ϵ+h.c.), =κ^ε_x + _xX _XX^-1κ^ε_XX _XX^-1 _Xx - ( _xX _XX^-1κ^ε_Xx +h.c. )\,, where κxx′ϵ=κxx′u+κxx′gκ^ε_x =κ^u_x +κ^g_x . The sample mean square 12ϵ¯2(x;ϑ∞,ψ0)=12K∑i=1Kϵi2(x;ϑ∞,ψ0) 12 ε^2(x; _∞, _0)= 12K _i=1^K _i^2(x; _∞, _0) is then known to follow a scaled Chi-squared distribution with K degrees of freedom 12ϵ¯2(x;ϑ∞,ψ0)∼12Σxxϵ(ϑ∞,ψ0)Kχ2(K) 12 ε^2(x; _∞, _0) 12 ^ε_x( _∞, _0)Kχ^2(K) (69) where Σxxϵ(ϑ∞,ψ0) ^ε_x( _∞, _0) is the variance of the GP described in Theorem 3.1. Conversely, a set of K+1K+1 independent neural networks arranged to a deep ensemble f(x;θ∞i)i=1K+1\f(x; _∞^i)\_i=1^K+1 in the infinite width limit n→∞n→∞ and at convergence t→∞t→∞ are by definition i.i.d. samples from the GP described in Theorem 2.1. As before, the empirical variance defined as σ¯f2(x;θ∞i…K+1)=1K∑i=1K+1(f(x;θ∞i)−1K+1∑j=1K+1f(x;θ∞j))2 σ_f^2(x; _∞^i… K+1)= 1K _i=1^K+1 (f(x;θ^i_∞)- 1K+1 _j=1^K+1f(x;θ^j_∞) )^2 distributes as a scaled Chi-squared distribution with K degrees of freedom σ¯f2(x;θ∞i…K+1)∼Σxxf(θ∞)Kχ2(K), σ_f^2(x; _∞^i… K+1) ^f_x( _∞)Kχ^2(K)\,, (70) where Σxxf(θ∞) ^f_x( _∞) is the variance of the GP described in Theorem 2.1. Finally, as we assume equal architecture and i.i.d. initialization of u, g, and f, we have that κxx′ϵ=κxx′u+κxx′g=2κxx′u=2κxx′κ^ε_x =κ^u_x +κ^g_x =2κ^u_x =2 _x and accordingly 12Σxxϵ(ϑ∞,ψ0)=Σxxf(θ∞) 12 ^ε_x( _∞, _0)= ^f_x( _∞), completing the proof. ∎ B.2 Posterior Equivalence This section contains proofs for results pertaining to the equivalence of self-predictive errors of “Bayesian RND” and the variance of Bayesian posterior predictive distributions of neural networks in the infinite width limit. B.2.1 Proof of Proposition 4.1 We restate Proposition 4.1 for convenience. See 4.1 Proof. The proof will show that in the limit n→∞n→∞ the function g~(x;ϑ0,ψ0) g(x; _0, _0) converges to a GP g~(x;ϑ0,ψ0)∼(0,Θxx′≤L−1) g(x; _0, _0) (0, ^≤ L-1_x ) by Lévy’s continuity theorem, which we recall informally below. Theorem B.4. (Lévy’s continuity theorem) Let Znn=1∞\Z_n\_n=1^∞ be a sequence of ℝnR^n-valued random variables. Their characteristic functions φZn(t) _Z_n(t) for some t∈ℝnt ^n are given by φZn(t)=[eit⊤Zn], _Z_n(t)=E[e^it Z_n]\,, (71) where i is the imaginary unit. If in the limit n→∞n→∞ the sequence of characteristic functions converges pointwise to a function φZn(t)→φ(t)∀t∈ℝn, _Z_n(t)→ (t) ∀ t ^n\,, (72) then ZnZ_n converges in distribution to a random variable Z Zn→Z, Z_n D→Z\,, (73) whose characteristic function is φZ(t)=φ(t) _Z(t)= (t) Rigorous proof can be found for example in durrett2019probability. We begin by rewriting the function g~(x;ϑ0,ψ0) g(x; _0, _0) as a linear model with g~(x;ϑ0,ψ0) g(x; _0, _0) =∇ϑu(x;ϑ0)⊤ψ0∗ = _ u(x; _0) _0^* (74) =∇ϑ≤L−1u(x;ϑ0)⊤ψ0≤L−1. = _ ^≤ L-1u(x; _0) _0^≤ L-1\,. (75) Since ψ0≤L−1 _0^≤ L-1 is an independent draw from ϑ0 _0 by assumption, g~(x;ϑ0,ψ0) g(x; _0, _0) is a random affine transform of the Gaussian vector ψ0≤L−1 _0^≤ L-1. For more precise treatment of the distribution of g~(x;ϑ0,ψ0) g(x; _0, _0), we write G~(T) G(X_T) to denote the random variable corresponding to the function evaluations of g~ g on a test set TX_T. Conditioned on ϑ0 _0 (i.e., fixing the affine transform), we thus have that G~(T)|ϑ0∼(0,Θ0,TT≤L−1) G(X_T)| _0 (0, ^≤ L-1_0,X_TX_T), where Θ0,TT≤L−1=∇ϑ≤L−1u(T;ϑ0)⊤∇ϑ≤L−1u(T;ϑ0) ^≤ L-1_0,X_TX_T= _ ^≤ L-1u(X_T; _0) _ ^≤ L-1u(X_T; _0) is the empirical NTK matrix of u. Note that this statement holds irrespective of the network width n. Next, we show that the unconditional law of G~(T) G(X_T), too, tends to a GP in the limit n→∞n→∞. To this end, we examine the distribution of the unconditioned random vector G~(T) G(X_T) through its characteristic function φG~(T)(t)=[eit⊤G~(T)]. _ G(X_T)(t)=E[e^it G(X_T)]\,. (76) This characteristic function φG~(T)(t) _ G(X_T)(t) uniquely defines the distribution of G~(T) G(X_T) [durrett2019probability]. By the law of total expectation, the characteristic function of the unconditional variable G~(T) G(X_T) can then be written as φG~(T)(t)=ϑ0[[eit⊤G~(T)|ϑ0]]. _ G(X_T)(t)=E_ _0 [E[e^it G(X_T)| _0] ]\,. (77) As stated above, the conditional distribution of G~(T)|ϑ0 G(X_T)| _0 is a zero-mean Gaussian with the empirical covariance Θ0,TT≤L−1 ^≤ L-1_0,X_TX_T, to which we can show the conditional characteristic function is given by [durrett2019probability] [eit⊤G~(T)|ϑ0]=e−12t⊤Θ0,TT≤L−1t. [e^it G(X_T)| _0]=e^- 12t ^≤ L-1_0,X_TX_Tt\,. (78) Plugging this back into Eq. 77 gives φG~(T)(t)=ϑ0[e−12t⊤Θ0,TT≤L−1t]. _ G(X_T)(t)=E_ _0[e^- 12t ^≤ L-1_0,X_TX_Tt]\,. (79) We now use the known result by jacotNeuralTangentKernel2020 that, as n→∞n→∞ we have that Θ0,TT→ΘTT _0,X_TX_T→ _X_TX_T in probability and accordingly Θ0,TT≤L−1→ΘTT≤L−1 ^≤ L-1_0,X_TX_T→ ^≤ L-1_X_TX_T converges to a deterministic kernel matrix. Moreover, since the Gram matrix Θ0,TT≤L−1 ^≤ L-1_0,X_TX_T is positive semidefinite in general, the term e−12t⊤Θ0,TT≤L−1te^- 12t ^≤ L-1_0,X_TX_Tt is bounded and continuous. By bounded convergence [durrett2019probability], we can then conclude that we also have convergence of the characteristic function through limn→∞φG~(T)(t) _n→∞ _ G(X_T)(t) =limn→∞ϑ0[e−12t⊤Θ0,TT≤L−1t] = _n→∞E_ _0[e^- 12t ^≤ L-1_0,X_TX_Tt] (80) =e−12t⊤ΘTT≤L−1t. =e^- 12t ^≤ L-1_X_TX_Tt\,. (81) As stated earlier, for a Gaussian random vector Z with Z∼(0,ΘTT≤L−1)Z (0, _X_TX_T^≤ L-1) its characteristic function is given by e−12t⊤ΘTT≤L−1te^- 12t ^≤ L-1_X_TX_Tt. Invoking Lévy’s continuity theorem, the pointwise convergence of φG~(T)(t) _ G(X_T)(t) to this exact limit φG~(T)(t)→e−12t⊤ΘTT≤L−1t _ G(X_T)(t)→ e^- 12t ^≤ L-1_X_TX_Tt then implies convergence in distribution of G~(T)→Z G(X_T) D→Z and we can thus conclude g~(x;ϑ0,ψ0)∼(0,Θxx′≤L−1) g(x; _0, _0) (0, ^≤ L-1_x ). ∎ B.2.2 Proof of Theorem 4.2 We restate Proposition 4.1 for convenience. See 4.2 Proof. The result follows from the independence of the two GP of interest in the limit n→∞n→∞. First, this is g~(x;ϑ0,ψ0)∼(0,Θxx′≤L−1) g(x; _0, _0) (0, _x ^≤ L-1) and second, u(x;ϑ0)∼(0,Θxx′L)u(x; _0) (0, _x ^L). In the following, we will show that the two GPs are in the limit n→∞n→∞ independent processes such that Eq. 14 applies. We first write for any two points x,x′x,x the covariance Cov[g~(x;ϑ0,ψ0),u(x′;ϑ0)]=[g~(x;ϑ0,ψ0)u(x′;ϑ0)]. [ g(x; _0, _0),u(x ; _0)]=E[ g(x; _0, _0)u(x ; _0)]\,. (82) As ψ0 _0 is drawn independently of ϑ0 _0, the conditional expectation can be written as [g~(x;ϑ0,ψ0)u(x′;ϑ0)|ϑ0] [ g(x; _0, _0)u(x ; _0)| _0] =u(x′;ϑ0)[g~(x;ϑ0,ψ0)|ϑ0] =u(x ; _0)E[ g(x; _0, _0)| _0] (83) =u(x′;ϑ0)[∇ϑ≤L−1u(x;ϑ0)⊤ψ0≤L−1|ϑ0] =u(x ; _0)E[ _ ^≤ L-1u(x; _0) _0^≤ L-1| _0] (84) =u(x′;ϑ0)⋅0, =u(x ; _0)· 0\,, (85) and by the law of total expectation [g~(x;ϑ0,ψ0)u(x′;ϑ0)] [ g(x; _0, _0)u(x ; _0)] =ϑ0[[g~(x;ϑ0,ψ0)u(x′;ϑ0)|ϑ0]] =E_ _0 [E[ g(x; _0, _0)u(x ; _0)| _0] ] (86) =0. =0\,. (87) We conclude that the two GP g~(x;ϑ0,ψ0)∼(0,Θxx′≤L−1) g(x; _0, _0) (0, _x ^≤ L-1) and u(x;ϑ0)∼(0,Θxx′L)u(x; _0) (0, _x ^L) are mutually independent such that the initialization kernel κxx′ϵbκ^ε^b_x is given as κxx′ϵb=Θxx′. κ^ε^b_x = _x \,. (88) This is because Θxx′=Θxx′L+Θxx′≤L−1 _x = ^L_x + ^≤ L-1_x and κxx′g~=Θxx′≤L−1κ g_x = ^≤ L-1_x , κxx′u=Θxx′Lκ^u_x = ^L_x are mutually independent. ∎ Appendix C Additional Experimental Details We report additional experimental details and evaluations. As outlined in the main text, we use two-layer fully connected neural networks with SiLU activations and NTK parametrization. All weights and biases are initialized as θ∼(0,I)θ (0,I). We use an ensemble of 512 models and a single multiheaded RND network with 512 heads. A synthetic dataset is generated with N=10N=10 train and N~=5000 N=5000 test samples from an isotropic Gaussian x∼(0,I3)x (0,I_3). We label training samples with a synthetic target function y(x)=x0+x1+x2−2∏i=13xi,y(x)=x^0+x^1+x^2-2 _i=1^3x^i, (89) where xix^i denotes the i-th component of vector x. All models are trained according to the algorithms outlined in the main text. For this, we use full-batch gradient descent with a learning rate of 0.10.1 for all models. Fig. 2 shows additional results of the same experiment, in which we plot individual test-set ensemble variances against RND errors. As the network width increases, ensemble variances and self-predictive RND errors become more correlated and well-calibrated in scale. Code for full reproduction will be released upon publication. (a) (b) Figure 2: (a) Scatter plot of test-set errors between predictive variances of ensembles and self-predictive errors of RND. As width increases, errors become more correlated and correctly calibrated in scale. (b) Likewise, for Bayesian ensembles and Bayesian RND.