On Infinite-Width Hypernetworks
Etai Littwin Tomer Galanti Lior Wolf Greg Yang
cs.LG stat.ML



















 
    
       
            
       
       
       

      





    theoremTheoremlemmaLemmacorollaryCorollarydefinitionDefinitionremarkRemark[table]skip=10pt plus 0.01pt
        On Infinite-Width Hypernetworks
            
        Etai LittwinEqual Contribution
        
        School of Computer Science 
        
        Tel Aviv University 
        
        Tel Aviv, Israel 
        
        
        Tomer Galanti^*
        
        School of Computer Science 
        
        Tel Aviv University 
        
        Tel Aviv, Israel 
        
        
        Lior Wolf 
        
        School of Computer Science 
        
        Tel Aviv University 
        
        Tel Aviv, Israel 
        
        
        Greg Yang 
        
        Microsoft Research AI 
        
        
            December 30, 2023
        ===========================================================================================================================================================================================================================================================================================================================
        
         Hypernetworks are architectures that produce the weights of a task-specific  primary network. A notable application of hypernetworks in the recent literature involves learning to output functional representations. In these scenarios, the hypernetwork learns a representation corresponding to the weights of a shallow MLP, which typically encodes shape or image information. While such representations have seen considerable success in practice, they remain lacking in the theoretical guarantees in the wide regime of the standard architectures. In this work, we study wide over-parameterized hypernetworks. We show that unlike typical architectures, infinitely wide hypernetworks do not guarantee convergence to a global minima under gradient descent. We further show that convexity can be achieved by increasing the dimensionality of the hypernetwork's output, to represent wide MLPs. In the dually infinite-width regime, we identify the functional priors of these architectures by deriving their corresponding GP and NTK kernels, the latter of which we refer to as the  hyperkernel. As part of this study, we make a mathematical contribution by deriving tight bounds on high order Taylor expansion terms of standard fully connected ReLU networks.
        
        
        
        
        
        
        § INTRODUCTION
        r0.28
            < g r a p h i c s >
        
        In this work, we analyze the training dynamics of over-parameterized meta networks, which are networks that output the weights of other networks, often referred to as  hypernetworks.
        In the typical framework, a function h involves two networks, f and g. The hypernetworkf takes the input x (typically an image) and returns the weights of the primary network, g, which then takes the input z and returns the output of h.
        
        The literature of hypernetworks is roughly divided into two main categories. In the functional representation literature <cit.> the input to the hypernetwork f is typically an image. For shape reconstruction tasks, the  network g represents the shape via a signed distance field, where the input are coordinates in 3D space. In image completion tasks, the inputs to g are image coordinates, and the output is the corresponding pixel intensity.In these settings, f is typically a large network and g is typically a shallow fully connected network.
        In the second category <cit.>, hypernetworks are typically used for hyper-parameter search, where x is being treated as a hyperparameter descriptor and is optimized alongside with the network's weights. In this paper, we consider models corresponding to the first group of methods.
        
        
        
        Following a prominent thread in the recent literature, our study takes place in the regime of wide networks. <cit.> recently showed that, when the width of the network approaches infinity, the gradient-descent training dynamics of a fully connected network f can be characterized by a kernel, called the Neural Tangent Kernel (or NTK for short). In other words, as the width of each layer approaches infinity,
        provided with proper scaling and initialization of the weights, it holds that:
        
            ∂ f(x;w)/∂ w·∂^⊤ f(x';w)/∂ w→Θ^f(x,x')
        
        as the width, n of f tends to infinity. Here, w are the weights of the network f(x;w). As shown in <cit.>, as the width tends to infinity, when minimizing the squared loss using gradient descent, the evolution through time of the function computed by the network follows the dynamics of kernel gradient descent with kernel Θ^f. To prove this phenomenon, various papers <cit.> introduce a Taylor expansion of the network output around the point of initialization and consider its values.
        It is shown that the first-order term is deterministic during the SGD optimization and the higher-order terms converge to zero as the width n tends to infinity. 
        
        
        
        
        
        
        
        
        
        A natural question that arises when considering hypernetworks is whether a similar “wide” regime exists, where trained and untrained networks may be functionally approximated by kernels. If so, since this architecture involves two networks, the “wide” regime needs a more refined definition, taking into account both networks.
        
        Our contributions:
        
          *  We show that infinitely wide hypernetworks can induce highly non-convex training dynamics under gradient descent. The complexity of the optimization problem is highly dependent on the architecture of the primary network g, which may considerably impair the trainability of the architecture if not defined appropriately.
            
            
          *  However, when the widths of both the hypernetwork f and the primary network g tend to infinity, the optimization dynamics of the hypernetwork simplifies, and its neural tangent kernel (which we call the hyperkernel) has a well defined infinite-width limit governing the network evolution.
            
          *  We verify our theory empirically and also demonstrate the utility of this hyperkernel on several functional representation tasks. Consistent with prior observations on kernel methods, the hypernetwork induced kernels also outperforms a trained hypernetwork when training data is small. 
            
            
          *  We make a technical contribution by deriving asymptotically tight bounds on high order Taylor expansion terms in ReLU MLPs. Our result partially settles a conjecture posed in <cit.> regarding the asymptotic behavior of general correlation functions.
            
        
        
        
        
         §.§ Related Works
        
        
          
        Hypernetworks
        Hypernetworks were first introduced under this name in <cit.>, are networks that generate the weights of a second  primary network that computes the actual task. However, the idea of having one network predict the weights of another was proposed earlier and has reemerged multiple times <cit.>. 
        The tool can naturally be applied for image representations tasks. 
        In <cit.>, they applied hypernetworks for 3D shape reconstruction from a single image. In <cit.> hypernetworks were shown to be useful for learning shared image representations. Hypernetworks were also shown to be effective in non-image domains. For instance, hypernetworks achieve state of the art results on the task of decoding error correcting codes <cit.>.
        Several publications consider a different framework, in which, the inputs x of the hypernetwork are optimized alongside to the weights of the hypernetwork. In this setting, hypernetworks were recently used for continuous learning by <cit.>.  Hypernetworks can be efficiently used for neural architecture search, as was demonstrated by <cit.>, where a feedforward regression (with network f) replaces direct gradient-based learning of the weights of the primary network while its architecture is being explored. Lorraine et al. applied hypernetworks for hyperparameters selection <cit.>.
        
        Despite their success and increasing prominence, little theoretical work was done in order to better
        understand hypernetworks and their behavior. A recent paper <cit.> studies the role of multiplicative interaction within a unifying framework to describe a range of classical and modern neural network architectural motifs, such as gating, attention layers,
        hypernetworks, and dynamic convolutions amongst others. It is shown that standard neural networks are a strict subset of neural networks with multiplicative interactions. In <cit.> the authors theoretically study the modular properties of hypernetworks. In particular, they show that compared to standard embedding methods, hypernetworks are exponentially more expressive when the primary network is of small complexity. In this work, we provide a complementary perspective and show that a shallow primary network is a requirement for successful training. <cit.> showed that applying standard initializations on a hypernetwork produces sub-optimal initialization of the primary network. A principled technique for weight initialization in hypernetworks is then developed.
        
        
        
        
        
          
        Gaussian Processes and Neural Tangent Kernel
        The connection between infinitely wide neural networks, Gaussian processes and kernel methods, has been the focus of many recent papers <cit.>. Empirical support has demonstrated the power of CNTK (convolutional neural tangent kernel) on popular datasets, demonstrating new state of the art results for kernel methods <cit.>. <cit.> showed that ReLU ResNets <cit.> can have NTK convergence occur even when depth and width simultaneously tend to infinity, provided proper initialization. In this work, we extend the kernel analysis of networks to hypernetworks, and characterize the regime in which the kernels converge and training dynamics simplify.
        
        
        
        
        
        
        § SETUP
        
        
        In this section, we introduce the setting of the analysis considered in this paper. We begin by defining fully connected neural networks and hypernetworks in the context of the NTK framework.  Neural networks    In the NTK framework, a fully connected neural network, f(x;w) = y^L(x), is defined in the following manner: 
        
            y^l(x) = √(1/n_l-1) W^l q^l-1(x) 
            
            q^l(x) = √(2)·σ(y^l(x))
             and  q^0(x) = x  ,
        
        where σ : ℝ→ℝ is the activation function of f. Throughout the paper, we specifically take σ to be a piece-wise linear function with a finite number of pieces (e.g., the ReLU activation ReLU(x) := max(0,x) and the Leaky ReLU activation ReLU_α(x) = 
         
        x      if  x ≥ 0
        α x     if  x<0). The weight matrices W^l ∈ℝ^n_l × n_l-1 are trainable variables, initialized independently according to a standard normal distribution, W^l_i,j∼𝒩(0,1). The width of f is denoted by n := min(n_1,…,n_L-1). The parameters w are aggregated as a long vector w = (vec(W^1),…,vec(W^L)). The coefficients √(1/n_l-1) serve for normalizing the activations of each layer. This parametrization is nonstandard, and we will refer to it as the NTK parameterization. It has already been employed in several recent works <cit.>. For simplicity, in many cases, we will omit to specify the weights w associated with our model.
        
         Hypernetworks    Given the input tuple u = (x,z) ∈ℝ^n_0 + m_0, we consider models of the form: h(u;w) := g(z;f(x;w)), where f(x;w)  and g(z;v) are two neural network architectures with depth L and H respectively. The function f(x;w) referred to as  hypernetwork, takes the input x and computes the weights v = f(x;w) of a second neural network g(z;v), referred as the  primary network, which is assumed to output a scalar. As before, the variable w stands for a vector of trainable parameters (v is not trained directly and is given by f).
        
        We parameterize the primary network g(z;v) = g^H(z;v) as follows:
        
            g^l(z;v) = √(1/m_l-1) V^l · a^l-1(z;v) 
            
            a^l(z;v) = √(2)·ϕ(g^l(z;v))
             and  a^0(z) = z
        
        Here, the weights of the primary network V^l(x) ∈ℝ^m_l × m_l-1 are given in a concatenated vector form by the output of the hypernetwork f(x;w) =  v = (vec(V^1),…,vec(V^H)). The output dimension of the hypernetwork f is therefore n_L = ∑^H_i=1 m_i· m_i-1. We denote by f^d(x;w) := V^d(x;w) := V^d the d'th output matrix of f(x;w). The width of g is denoted by m := min(m_1,…,m_H-1). The function ϕ is an element-wise continuously differentiable function or a piece-wise linear function with a finite number of pieces.
        
         Optimization    Let S = {(u_i,y_i)}^N_i=1, where u_i = (x_i,z_i) be some dataset and let ℓ(a,b) := |a-b|^p be the ℓ^p-loss function. For a given hypernetwork h(u;w), we are interested in selecting the parameters w that minimize the empirical risk:
        
            c(w) 
            
            := ∑^N_i=1ℓ(h(u_i;w), y_i)
        
        
        For simplicity, oftentimes we will simply write ℓ_i(a) := ℓ(a,y_i) and h_i(w) := h(u_i) := h(u_i;w), depending on the context. In order to minimize the empirical error c(w), we consider the SGD method with learning rate μ > 0 and step of the form w_t+1← w_t - μ∇_w ℓ_j_t(h_j_t(w_t))for some index j_t ∼ U[N] that is selected uniformly at random for the t'th iteration. A continuous version of the GD method is the gradient flow method, in which ẇ = -μ∇_w c(w). In recent works <cit.>, the optimization dynamics of the gradient method for standard fully-connected neural networks was analyzed, as the network width tends to infinity. In our work, since hypernetworks consist of two interacting neural networks, there are multiple ways in which the size can tend to infinity. We consider two cases: (i) the width of f tends to infinity and that of g is fixed and (ii) the width of both f and g tend to infinity. 
        
        
        
        
        
        
        
        
        § DYNAMICS OF HYPERNETWORKS
        
        
          
        Infinitely wide f without infinitely wide g induces non-convex optimization
        In the NTK literature, it is common to adopt a functional view of the network evolution by analyzing the dynamics of the output of the network along with the cost, typically a convex function, as a function of that output. 
        
        
        
        
        
        
        In the hypernetwork case, this presents us with two possible viewpoints of the same optimization problem of h(u)=g(z;f(x)). On one hand, since only the hypernetwork f contains the trainable parameters, we can view the optimization of h under the loss ℓ as training of f under the loss ℓ∘ g.
        The classical NTK theory would imply that f evolves linearly when its width tends to infinity, but because ℓ∘ g is in general not convex anymore, even when ℓ originally is, an infinitely wide f without an infinitely wide g does not guarantee convergence to a global optimum.
        
        In what follows, we make this point precise by characterising how nonlinear the dynamics becomes in terms of the depth of g.
        
        
        
        
        
        
        After a single stochastic gradient descent step with learning rate μ, the hypernetwork output for example i is given by h_i(w - μ∇_w ℓ_j). When computing the Taylor approximation around w with respect to the function h at the point w' = w - μ∇_w ℓ_j, it holds that:
        
        
        
        
        
        
        
            h_i(w - μ∇_w ℓ_j) 
            = ∑_r=0^∞1/r!⟨∇^(r)_w h_i, (-μ∇_w ℓ_j)^r⟩  = ∑_r=0^∞1/r!(-μ∂ℓ_j/∂ h_j)^r ·𝒦^(r)_i,j
        
        where 𝒦^(r)_i,j := ⟨∇^(r)_w h_i, (∇_w h_j)^r⟩,
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        and ∇^(r)_w h_i is the r tensor that holds the r'th derivative of the output h_i. The terms ⟨∇^(r)_w h_i, (-μ∇_w ℓ_j)^r⟩ are the multivariate extensions of the Taylor expansion terms h^(r)_i(w) (w'-w)^r, and take the general form of correlation functions as introduced in Eq. 5 in the appendix. This equation holds for neural networks with smooth activation functions (including hypernetworks), and holds in approximation for piece-wise linear activation functions.
        
        
        Previous works have shown that, if h is a wide fully connected network, the first order term (r=1) converges to the NTK, while higher order terms (r>1) scale like 𝒪(1/√(n)) <cit.>.
        Hence, for large widths and small learning rates, these higher order terms vanish, and the loss surface appears deterministic and linear at initialization, and remains so during training.
        
        
        However, the situation is more complex for hypernetworks.
        As shown in the following theorem, for infinitely wide hypernetworks and finite primary network, the behaviour depends on the depth and width of the generated primary network. Specifically, when the primary network is deep and narrow, the higher order terms in Eq. <ref> may not vanish, and parameter dynamics can be highly non-convex.
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        [Higher order terms for hypernetworks]theoreminithypernet
        Let h(u) = g(z;f(x)) for a hypernetwork f and an primary network g. Then, we have:
        
        
            𝒦^(r)_i,j∼
                 n^H-r   if r > H
            
                 1    otherwise.
        
        
        
        
        
        
        
        Thm. <ref> illustrates the effect of the depth of the primary network g on the evolution of the output h. The larger H is, the more non-linear the evolution is, even when f is infinitely wide. Indeed, we observe empirically that when f is wide and kept fixed, a deeper g incurs slower training, and lower overall test performance as illustrated in Fig. <ref>. 
        
        
        As a special case of this theorem, when taking H=1, we can also derive the asymptotic behaviour of 𝒦^(r)_i,j∼ n^1-r for a neural network h. This provides a tighter bound than the previously conjectured 𝒪(1/n) upper bound <cit.>. The following remark is a consequence of this result and is validated in the supplementary material.
        
        
        
        
        The r'th order term of the Taylor expansion in Eq. <ref> is of order 𝒪(μ^r/r! · n^r-1) instead of the previously postulated 𝒪(μ^r/r! · n). Therefore, it is evident that for any choice μ = o(√(n)), all of the high order terms tend to zero as n →∞. This is opposed the previous bound, which guarantees that all of the high order terms tend to zero as n →∞ only when μ is constant. 
        
        
        
        
        
        
        
        
        
        
        § DUALLY INFINITE HYPERNETWORKS
        It has been shown by <cit.> that NNGPs and neural tangent kernels fully characterise the training dynamics of infinitely wide networks. As a result, in various publications <cit.>, these kernels are being treated as functional priors of neural networks. In the previous section, we have shown that the Taylor expansion of the hypernetwork is non-linear when the size of the primary network is finite. In this section, we consider the case when both hyper and primary networks are infinitely wide, with the intention of gaining insight into the functional prior of wide hypernetworks. For this purpose, we draw a formal correspondence between infinitely wide hypernetworks and GPs, and use this connection to derive the corresponding neural tangent kernel. 
        
        
        
         §.§ The NNGP kernel
        
        
        Previous work have shown the equivalence between popular architectures, and Gaussian processes, when the width of the architecture tends to infinity. This equivalence has sparked renewed interest in kernel methods, through the corresponding NNGP kernel, and the Neural Tangent Kernel (NTK) induced by the architecture, which fully characterise the training dynamics of infinitely wide networks. This equivalence has recently been unified to encompass most architectures which use a pre-defined set of generic computational blocks <cit.>. Hypernetworks represent a different class of neural networks where the parameters contain randomly initialized matrices except the last layer whose parameters are aggregated as a rank 3 tensor. All of the matrices/tensors dimensions tend to infinity.
        This means the results of <cit.> do not apply to hypernetworks.
        
        Nevertheless, by considering sequential limit taking, where we take the limit of the width of f ahead of the width of g, we show the output of f achieves a GP behaviour, essentially feeding g with Gaussian distributed weights with adaptive variances.
        A formal argument is presented in the following theorem.
        [Hypernetworks as GPs]theoremGP
        Let h(u) = g(z;f(x)) be a hypernetwork. For any pair of inputs u = (x,z) and u' = (x',z'), let  Σ^0(z,z') = z^⊤ z'/m_0, S^0(x,x') = x^⊤ x'/n_0. Then, it holds for any unit i in layer 0<l≤ H of the primary network:
        
            g^l_i(z;f(x)) d⟶𝒢_i^l(u)
        
        as m,n →∞ sequentially. Here, {𝒢_i^l(u)}_i=1^m_l are independent Gaussian processes, such that, (𝒢_i^l(u), 𝒢_i^l(u')) ∼𝒩(0,Λ^l(u,u')) defined by the following recursion:
        
            Λ^l+1(u,u') = 
            [   Σ^l(u,u)  Σ^l(u',u);  Σ^l(u,u') Σ^l(u',u') ][   S^L(x,x)  S^L(x',x);  S^L(x,x') S^L(x',x') ]
        
            Σ^l(u,u') =  2_(u,v)∼𝒩(0,Λ^l)[σ(u)·σ(v)]
        
        where S^L(x,x') is defined recursively:
        
            S^l(x,x') = 2_(u,v)∼𝒩(0,Γ^l)[σ(u)·σ(v)]   and  Γ^l(x,x') = [   S^l(x,x)  S^l(x',x);  S^l(x,x') S^l(x',x') ]
        
        
        In other words, the NNGP kernel, governing the behaviour of wide untrained hypernetworks, is given by the Hadamard product of the GP kernels of f and g (see Eq. <ref>). 
        
        As a consequence of the above theorem, we observe that the NNGP kernel of h at each layer, Λ^l(u,u'), is simply a function of Σ^0(z,z'), S^0(x,x').
        
        corollarycompGP
        Let h(u) = g(z;f(x)) be a hypernetwork.
        
        For any 0< l≤ H, there exists a function ℱ^l, such that, for all pairs of inputs u = (x,z) and u' = (x',z'), it holds that:
        
            Λ^H(u,u') = ℱ(Σ^0(z,z'),S^0(x,x'))
        
         
        The factorization of the NNGP kernel into a function of Σ^0(z,z') and S^0(x,x') provides a convenient way to explicitly encode useful invariances into the kernel.
        As an example, in the following remark, we investigate the behaviour of the NNGP kernel of h, when the inputs z are preprocessed random random Fourier features as suggested by <cit.>. 
        
        remarkperiodicGP
        Let p(z) = [cos(W^1_i z + b^1_i)]^k_i=1 be a Fourier features preprocessing, where W^1_i,j∼𝒩(0,1) and biases b_i∼ U[-π,π]. Let h(u) = g(p(z);f(x)) be a hypernetwork, with z preprocessed according to p. Let u = (x,z) and u'=(x',z') be two pairs of inputs. Then,  Λ^l(u,u') is a function of exp[-z-z'^2_2/2] and S^L(x,x'). 
        
        
        The above remark shows that for any given inputs x,x', the NNGP kernel depends on z,z' only through the distance between z and z', which has been shown to be especially useful in implicit neural representation <cit.>.
        
        
        
        
        
        
        
        
        
        
        We next derive the corresponding neural tangent kernel of hypernetworks, referred to as hyperkernels.
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
         §.§ The Hyperkernel
        
        
        Recall the definition of the NTK as the infinite width limit of the Jacobian inner product given by:
        
            𝒦^h(u,u') = ∂ h(u)/∂ w·∂^⊤ h(u')/∂ w = ∂ g(z;f(x))/∂ f(x)·𝒦^f(x,x') ·∂^⊤ g(z';f(x'))/∂ f(x')
        
        where 𝒦^f(x,x') := ∂ f(x)/∂ w·∂^⊤ f(x')/∂ w and 𝒦^g(u,u') := ∂ g(z;f(x))/∂ f(x)·∂^⊤ g(z';f(x'))/∂ f(x'). In the following theorem we show that 𝒦^h(u,u') converges in probability at initialization to a limiting kernel in the sequentially infinite width limit of f and g, denoted by Θ^h(u,u'). 
        
        Furthermore, we show that the hyperkernel is decomposed to the Hadamard product between the kernels corresponding to f and g. In addition, we show that the derivative of the hyperkernel with respect to time tends to zero at initialization.
        
        [Hyperkernel decomposition and convergence at initialization]theoremhyperkernel
        Let h(u;w) = g(z;f(x;w)) be a hypernetwork. Then,
        
        
            𝒦^h(u,u') p⟶Θ^h(u,u')
        
        where:
        
            Θ^h(u,u') = Θ^f(x,x') ·Θ^g(u,u',S^L(x,x'))
        
        such that:
        
            𝒦^f(x,x') p⟶Θ^f(x,x')· I  and 𝒦^g(u,u') p⟶Θ^g(u,u',S^L(x,x'))
        
        
        
        
        moreover, if w evolves throughout gradient flow, we have:
        
            ∂𝒦^h(u,u')/∂ t|_t=0p⟶ 0
        
        where the limits are taken with respect to m,n→∞ sequentially.
        
        
        
        
        As a consequence of Thm. <ref>, when applying a Fourier features preprocessing to z, one obtains that Θ^g(u,u') becomes shift invariant. 
        
        remarkperiodicNTK
        Let p(z) be as in Remark <ref>. Let h(u) = g(p(z);f(x)) be a hypernetwork, where z is preprocessed according to p. Let u=(x,z) and u'=(x',z') be two pairs of inputs. Then, Θ^g(u,u') is a function of exp[-z-z'^2_2/2] and S^0(x,x').
        
        Note that Θ^f is the standard limiting NTK of f and depends only on the inputs {x_i}_i=1^N. However from Eq. <ref>, the term Θ^g requires the computation of the NNGP kernel of f in advance in order to compute the terms {Σ^l,Σ̇^l}. This form provides an intuitive factorization of the hyperkernel into a term Θ^f which depends on the meta function and data, and Θ^g which can be though of as a conditional term.
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        § EXPERIMENTS
        
        Our experiments are divided into two main parts. In the first part, we validate the ideas presented in our theoretical analysis and study the effect of the width and depth of g on the optimization of a hypernetwork. In the second part, we evaluate the performance of the NNGP and NTK kernels on image representation tasks. For further implementation details on all experiments see Appendix A.
        
         §.§ Convergence of the Hyperkernel
        
        
        
        We verified our results of Thm. <ref> by constructing a simple hypernetwork, for which both f and g are four layered fully connected networks with ReLU activations. For the input of f, we used a fixed 2D vector x = (1,-1). The input z(θ) of g varied according to z(θ) = (sin(θ),cos(θ)), where θ∈ [-π/2,π/2]. We then compute the empirical hyperkernel as follows, while varying the width of both f and g:
        
            𝒦^h(u,u') = ∇_w h(x,z(θ))·∇_w^⊤ h(x,z(θ))
        
        Results are presented in Fig. <ref>. As can be seen, convergence to a fixed kernel is only observed in the dually wide regime, as stated in Thm. <ref>.
        
        
         §.§ Training Dynamics
        
        
        We consider a rotation prediction task. In this task, the hypernetwork f is provided with a randomly rotated image x and the primary network g is provided with a rotated version z of x with a random angle α. The setting is cast into a classification task, where the goal is to predict the closest value to α/360 within {α_i = 30i/360 | i = 0,…,11}. We experimented with the MNIST <cit.> and CIFAR10 <cit.> datasets. For each dataset we took 10000 training samples only.
        
        We investigate the effect of the depth and width of g on the training dynamics of a hypernetwork. We compared the performance of hypernetworks of various architectures to investigate the effect of the depth and width of g. The architectures of the hypernetwork and the primary network are as follows. The hypernetwork, f, is a fully-connected ReLU neural network of depth 4 and width 200. The inputs of f are flattened vectors of dimension c · h^2, where c specifies the number of channels and h the height/width of each image (c=1 and h=28 for MNIST and c=3 and h=32 for CIFAR10). The primary network g is a fully-connected ReLU neural network of depth ∈{3,6,8}. Since the MNIST rotations dataset is simpler, we varied the width of g in ∈{10,50,100} and for the the CIFAR10 variation we selected the width of g to be ∈{100,200,300}. The network outputs 12 values and is trained using the cross-entropy loss. 
        
        We trained the hypernetworks for 100 epochs, using the SGD method with batch size 100 and learning rate μ=0.01. For completeness, we conducted a sensitivity study on the learning rate for both datasets, to show that the reported behaviour is consistent for any chosen learning rate, see appendix. In Fig. <ref>(a-c) we compare the performance of the various architectures on the MNIST and CIFAR10 rotations prediction tasks. The performance is computed as an average and standard deviation (error bars) over 100 runs. As can be seen, we observe a clear improvement in test performance as the width of g increases, especially at the initialization. When comparing the three plots, we observe that when f is wide and kept fixed, a deeper g incurs slower training, and lower overall test performance. This is aligned with the conclusions of Thm. <ref>. 
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
         §.§ Image representation and Inpainting
        We compare the performance of a hypernetwork and kernel regression with the hyperkernel on two visual tasks: functional image representation and inpainting. In the MNIST image representation task, the goal of the hypernetwork f is to represent an input image via the network g, which receives image pixel coordinates and outputs pixel values. In the inpainting task, the goal is the same where only half of the image is observed by f.
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
         Problem Setup  
        We cast these problems as a meta-learning problem, where f receives an image, and the goal of the primary network g:[28]^2 → [0,1] is then to learn a conditional mapping from pixel coordinates to pixel values for all the pixels in the image, with the MSE as the metric. Our training dataset S = {(u_i,y_i)}_i=1^N then consists of samples u_i = (x_i,z_i), such that, x_i are images, and z_i are random pixel location (i.e., a tuple ∈ [28]^2), and y_i is a label specifying the pixel value at the specified location (normalized between 0 and 1). In both experiments, training was done on randomly sampled training data of varying size, specified by N.
        
         Evaluation  
        We evaluate the performance of both training a hypernetwork, and using kernel regression with the hyperkernel. For kernel regression, we use the following formula to infer the pixel value of a test point u:
        
            (Θ^h(u,u_1),..., Θ^h(u,u_N))·(Θ^h(U,U) + ϵ· I)^-1· Y
        
        where Θ^h(U,U) = (Θ^h(u_i,u_j))_i,j∈ [N] is the hyperkernel matrix evaluated on all of the training data and Y = (y_i)^N_i=1 is the vector of labels in the training dataset and ϵ = 0.001. 
        
        In Tab. <ref>, we compare the results of the hyperkernel and the NNGP kernel with the corresponding hypernetwork. The reported numbers are averages over 20 runs. As can be seen, in the case of a small dataset, the kernels outperforms the hypernetwork, and the NNGP outperforms the rest.
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        § CONCLUSIONS
        
        
        
        In this paper, we apply the well-established large width analysis to hypernetwork type models. For the class of models analyzed, we have shown that a wide hypernetwork must be coupled with a wide primary network in order achieve a simplified, convex training dynamics as in standard architectures. 
        
        The deeper g is, the more complicated the evolution is. In the dually infinite case, when the widths of both the hyper and primary networks tend to infinity, the optimization of the hypernetwork become convex and is governed by the proposed hyperkernel.
        The analysis presented in this paper is limited to a specific type of hypernetworks used in the literature, typically found in the functional neural representation literature, and we leave the extension of this work to additional types of hyper models to future work.
        
        Some of the tools developed in this study, also apply for regular NTKs. Specifically, <cit.> provide a conjecture, for which one of its consequences is that 
        𝒦^(r)_i,j = 𝒪(1/n). In Thm. <ref> we prove that this hypothesized upper bound is increasingly loose as r increases, and prove an asymptotic behaviour in the order of 𝒦^(r)_i,j∼ 1/n^r-1. 
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        § BROADER IMPACT
        
        This work improves our understanding and design of hypernetworks and hopefully will help us improve the transparency of machine learning involving them. Beyond that, this work falls under the category of basic research and does not seem to have particular societal or ethical implications.
        
        § ACKNOWLEDGEMENTS AND FUNDING DISCLOSURE
        
        
        This project has received funding from the European Research Council (ERC) under the European
        Union’s Horizon 2020 research and innovation programme (grant ERC CoG 725974). The contribution of Tomer Galanti is part of Ph.D. thesis research conducted at Tel Aviv University.
        
        
        plain
        
        § IMPLEMENTATION DETAILS
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
         §.§ Convergence of the Hyperkernel
        
        
        
        
        In Fig. 1(a) (main text) we plot the variance of the kernel values 𝒦^h(u,u') in log-scale, as a function of the width of both f and g. The variance was computed empirically over k=100 normally distributed samples w. As can be seen, the variance of the kernel tends to zero only when both widths increase. In Fig. 1(b) (main text) we plot the value of 𝒦^h(u,u') and its variance for a fixed hypernetwork f of width 500 and g of width 10 or 800. The x-axis specifies the value of θ∈ [-π/2,π/2] and the y-axis specifies the value of the kernel. As can be seen, the expected value of the empirical kernel, 𝒦^h(u,u'), is equal to the width-limit kernel (e.g., theoretical kernel) for both widths 10 and 800. In addition, convergence of the width-limit kernel is guaranteed only when the widths of both networks increase, highlighting the importance of wide architectures for both the hyper and implicit networks for stable training. 
        
        
        
         §.§ Image Completion and Impainting
         Architectures  
        In both tasks, we used fully connected architectures, where f contains two hidden layers, and g contains one hidden layer. The hyperkernel used corresponds to the infinite width limit of the same architecture. For the input of g, we used random Fourier features <cit.> of the pixel coordinates as inputs for both the hyperkernel and the hypernetwork. To ease on the computational burden of computing the full kernel matrix Θ^h(U,U) when evaluating the hyperkernel, we compute smaller kernel matrices on subsets of the data {U^s_i} = {x_i^s,z_i^s}_s∈ [10], where each subset contains 1k input images {x^s_i}, and 20 random image coordinates per input, producing a kernel matrix of size 20k×20k. The final output prediction is then given by:  
        
            1/10∑_s(Θ^h(u,u_1^s),..., Θ^h(u,u_N^s))·(Θ^h(U^s,U^s) + ϵ· I)^-1· Y^s
        
        where Y^s are the corresponding labels of the subset U^s.
        For the hypernetwork evaluation, we used the same inputs {U^s_i}_s∈[10] to train the hypernetwork using a batchsize of 20, and a learning rate of 0.01 which was found to produce the best results.  
        
        
        
        § ADDITIONAL EXPERIMENTS
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
         §.§ Sensitivity Study
        
        
        To further demonstrate the behavior reported in Fig. 2 (main text), 
        we verified that it is consistent regardless of the value of the learning rate. We used the same architectures as in the rotations prediction experiments, i.e., f is a fully-connected ReLU neural network of depth 4 and width 200 and g is of depth ∈{3,6,8} and width ∈{50,100,200}. We vary the learning rate: μ = 10^j-7, for j=0,…,7. For each value of the learning rate, we report the average performance (and standard deviation over 100 runs) of the various architectures after 40 epochs of training. 
        
        As can be seen in Fig. <ref>, when f is wide and kept fixed, there is a clear improvement in test performance as the width of g increases, for every learning rate in which the networks provide non-trivial performance. When f is wide and kept fixed, a deeper g incurs slower training and lower overall test performance. We note that it might seem that the width of g does not affect the performance when the learning rate is μ=0.01 in all settings except Figs. <ref>(c,f). Indeed, we can verify from Fig. 2 (main text) that the performance at epoch 40 is indeed similar for different widths. However, for earlier epochs, the performance improves for shallower and wider architectures.
        
        
        
        
         §.§ Training wide networks with a large learning rate
         Remark 1 (main text) states that one is able to train wide networks with a learning rate μ = o(n). To validate this remark, we trained shallow networks of varying width n ∈{10^2,10^3,10^4,2.5· 10^5} with learning rate μ = √(n) on MNIST. As can be seen in Fig. <ref>, training those networks is possible despite the very large learning rate. In fact, we observe that the accuracy rate and loss improve as we increase the width of the network.
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        § CORRELATION FUNCTIONS
        
        
        Correlation functions are products of general high order tensors representing high order derivatives of a networks output with respect to the weights. 
        In <cit.> a conjecture is posed on the order of magnitude of general correlation functions involving high order derivative tensors, which arise when analysing the dynamics of gradient descent. 
        Roughly speaking, given inputs {x_i}_i=1^r, the outputs of a neural network f(x_1;w),...,f(x_r;w) ∈ℝ with normally distributed parameters w∈ℝ^N, correlation functions takes the form:
        
            ∑_η_k_0,...,η_k_r∈ [N]∏_j=1^rΓ_η_k_j+1,...,η_k_j+1(x_j)
        
        where 
        
            Γ_η_1,...,η_k(x_j) := ∂^k f(x_j;w)/∂ w_η_1...∂ w_η_k
        
        For instance, the following are two examples of correlation functions,
        
            f(x_1;w) ·∂ f(x_2;w)/∂ w_μ_1,
            ∂^2 f(x_1;w)/∂ w_μ_1∂ w_μ_2·∂ f(x_2;w)/∂ w_μ_1
        
        Computing the expected value of these correlation functions involve keeping track of various moments of normally distributed weights along paths, as done in recent finite width correction works <cit.>.  
        <cit.> employ the Feynman diagram to efficiently compute the expected values (order of magnitude) of general correlation functions, albeit at the cost of only being provably accurate for deep linear, or shallow ReLU networks. 
        
        In this work, we analyze the asymptotic behaviour correlation functions of the form: 
        
            𝒯^r(x_0,...,x_r)    := ∑_η_k_0...η_k_r∈ [N]Γ_η_k_1,...,η_k_r(x_0)∏_j=1^rΓ_η_k_j(x_j) 
               = ⟨∇^(r)_w f(x_0), ⊗^r_j=1∇_w f(x_j) ⟩
        
        
        
        
        
        where ∇^(r)_w f(x_0) is a rank r tensor, representing the r'th derivative of the output, and ⊗^r_j=1∇_w f(x_j) denotes outer products of the gradients for different examples. 
        terms of the form in Eq. <ref> represent high order terms in the multivariate Taylor expansion of outputs, and are, therefore, relevant for the full understanding of training dynamics. As a consequence of Thm. <ref>, we prove that 𝒯^r(x_0,...,x_r) ∼ 1/n^max(r-1,0) for vanilla neural networks, where n is the width of the network.
        As we have shown in Sec. 3, terms of the form in Eq. <ref> represent high order terms in the multivariate Taylor expansion of outputs, and are, therefore, relevant for the full understanding of training dynamics. As a consequence of Thm. <ref>, we prove that 𝒯^r(x_0,...,x_r) ∼ 1/n^max(r-1,0) for vanilla neural networks, where n is the width of the network.
        
        This result is a partial solution to the open problem suggested by <cit.>. In their paper, they conjecture the asymptotic behaviour of general correlation functions, and predict an upper bound on the asymptotic behaviour of terms of the form in Eq. <ref> in the order of 𝒪(1/n). Our results therefore proves a stronger version of the conjecture, while giving the exact behaviour as a function of width.
        
        
        
        
        § PROOFS OF THE MAIN RESULTS
         Terminology and Notations  
        Throughout the appendix, we denote by A ⊗ B and A ⊙ B the outer and Hadamard products of the tensors A and B (resp.). When considering the outer products of a sequence of tensors {A_i}^k_i=1, we denote, ⊗^k_i=1 A_i = A_1 ⊗…⊗ A_k. We denote by (x) := x/| x| the sign function. The notation X_n ∼ a_n states that X_n/a_n converges in distribution to some non-zero random variable X. 
        A convenient property of this notation is that it satisfies: X_n · Y_n ∼ a_n · b_n when X_n∼ a_n and Y_n ∼ b_n.
        
        Throughout the paper, we will make use of sequential limits and denote n_k,…,n_1 →∞ to express that n_1 tend to infinity, then n_2, and so on. For a given sequence of random variable {X_n}^∞_n=1, we denote by X_nd⟶ X (X_np⟶ X), when X_n converges in distribution (probability) to a random variable X.
        
        
         §.§ Useful Lemmas
        
        Let X_nd⟶ X. Then, (X_n) d⟶(X).
        
        We have: 
        
            lim_n→∞ℙ[(X_n) = 1] = lim_n→∞ℙ[X_n ≥ 0] = ℙ[X ≥ 0] = ℙ[(X) = 1]
        
        Hence, (X_n) converges in distribution to (X).
        
        
         §.§ Main Technical Lemma
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        In this section, we prove Lem. <ref>, which is the main technical lemma that enables us proving Thm. <ref>. Let f(x;w) be a neural network with H outputs {f^d(x;w)}^H_d=1. We would like to estimate the order of magnitude of the following expression:
        
            𝒯^l,i,d_n,i,d :=
            ⟨∂^k f^d(x_i;w)/∂ W^l_1…∂ W^l_k , ⊗^k_t=1∂ f^d_1(x_i_t;w)/∂ W^l_t⟩
        
        where d = (d_1,…,d_k), i=(i_1,…,i_k) and l = (l_1,…,l_k). For simplicity, when, i_1=…=i_k=j, we denote: 𝒯^l,d_n,i,j,d := 𝒯^l,i,d_n,i,d and 𝒯^l_n,i,j,d := 𝒯^l,i,d_n,i,d when  d_1=…=d_k=d as well. 
        
        To estimate the order of magnitude of the expression in Eq. <ref>, we provide an explicit expression for ∂^k f^d(x_i;w)/∂ W^l_1…∂ W^l_k. First, we note that for any w, such that, f^d(x_i;w) is k times continuously differentiable at w, for any set l := {l_1,…,l_k}, we have:
        
            ∂^k f^d(x_i;w)/∂ W^l_1…∂ W^l_k = ∂^k f^d(x_i;w)/∂ W^l'_1…∂ W^l'_k
        
        where the set l' := {l'_1,…,l'_k} is an ordered version of l, i.e., the two sets consist of the same elements but l'_1 < … < l'_k. In addition, we notice that for any multi-set l, such that, l_i = l_j for some i≠ j, then, 
        
            ∂^k f^d(x_i;w)/∂ W^l_1…∂ W^l_k = 0
        
        since f^d(x_i;w) is a neural network with a piece-wise linear activation function. Therefore, with no loss of generality, we consider l = {l_1,…,l_k}, such that, l_1<…<l_k.  It holds that:
        
            ∂^k f^d(x_i;w)/∂ W^l_1…∂ W^l_k = 1/√(n_l_1-1)q^l_1-1_i,d⊗𝒜^l_1→ l_2_i,d
        
        where 𝒜^l_1→ l_2_i,d is a 2k-1 tensor, defined as follows:
        
            𝒜^l_j→ l_j+1_i,d = 
            1/√(n_l_j+1-1)C^l_j→ l_j+1_i,d⊗𝒜^l_j+1→ l_j+2_i,d    1<j<k-1
            1/√(n_l_k-1)C^l_k-1→ l_k_i,d⊗ C^l_k→ L_i,d    j=k-1
        
        where:
        
            C^l_j→ l_j+1_i,d = √(2)Z^l_j+1-1_i,dP_i,d^l_j→ l_j+1 - 1    l_j+1≠ L
            
            P_i,d^l_j→ L    else
        
        and:
        
            P_i^u→ v = ∏_l = u^v-1(√(2/n_l)W^l+1Z^l_i)  and  Z^l_i = diag (σ̇(y^l(x_i)))
        
        The individual gradients can be expressed using:
        
            ∂ f^d_j_w(x_i_j)/∂ W^l_j= q^l_j-1_i_j,d_j⊗ C^l_j→ L_i_j,d_j/√(n_l_j-1)
        
        Note that the following holds for any u<v<h≤ L:
        
            C^u→ h_i,d = C^v→ hW^v/√(n_v-1)C^u→ v_i,d and  C^u→ L_i,d = C^v-1→ L_i,d P^u→ v-1_i,d
        
        In the following, given the sets l = {l_1,…,l_k}, i={i_1,…,i_k} and d = {d_1,…,d_k}, we derive the limit of 𝒯^l,i,d_n,i,d using elementary tensor algebra. By Eqs. <ref> and <ref>, we see that:
        
            𝒯^l,i,d_n,i,d   = ⟨⊗^k_t=1∂ f^d_t(x_i_t;w)/∂ W^l_t,q^l_1-1_i,d/√(n_l_1-1)⊗C^l_1→ l_2_i,d/√(n_l_2-1)⊗ ...⊗C^l_r-1→ l_k_i,d/√(n_l_k-1)⊗ C^l_k→ L_i,d⟩
               =1/n_l_1-1⟨ q^l_1-1_i,d , q^l_1-1_i_1,d_1⟩·⟨ C_i_k,d_k^l_k→ L , C_i,d^l_k→ L⟩∏_j=1^k-1⟨C_i_j,d_j^l_j→ L⊗ q_i_j+1,d_j+1^l_j+1-1/n_l_j+1-1, C_i,d^l_j→ l_j+1⟩
        
        
        We recall the analysis of <cit.> showing that in the infinite width limit, with n = min(n_1…,n_L-1) →∞, every pre-activation y^l(x) of f(x;w) at hidden layer l∈ [L] has all its coordinates tending to i.i.d. centered Gaussian processes of covariance Σ^l(x,x') : ℝ^n_0×ℝ^n_0→ℝ defined recursively as follows:
        
            Σ^0(x,x')    = x^⊤ x', 
            Λ^l(x,x')    = 
            [   Σ^l-1(x,x)  Σ^l-1(x,x');  Σ^l-1(x',x) Σ^l-1(x',x') ]∈ℝ^2× 2,
            Σ^l(x,x')    = 𝔼_(u,v)∼𝒩(0,Λ^l-1(x,x'))[σ(u)σ(v)]
        
        In addition, we define the derivative covariance as follows:
        
            Σ̇^l(x,x') = 𝔼_(u,v)∼𝒩(0,Λ^l-1(x,x'))[σ̇(u)σ̇(v)]
        
        when considering x = x_i and x' = x_j from the training set, we simply write Σ^l_i,j := Σ^l(x_i,x_j) and Σ̇^l_i,j = Σ̇^l(x_i,x_j).
        
        
        The following holds: 
        
            
          *  For n_L-1,…,n_1→∞, we have: P_i,d_1^u → L(P_j,d_2^u → L)^⊤d⟶∏_l=u^L-1Σ̇_i,j^lδ_d_1=d_2.
            
          *  For n_v,…,n_1→∞, we have: (q_i^v)^⊤ q_j^v/n_vd⟶Σ_i,j^v.
        
        Here, δ_T is an indicator that returns 1 if T is true and 0 otherwise.
        
        See <cit.>.
        
        Let k ≥ 0 and sets l={l_1,…,l_k}, i={i_1,…,i_k} and d={d_1,…,d_k}. We have:
        
            n^max(k-1,0)·𝒯^l,i,d_n,i,dd⟶δ_d·∏_j=1^k-1𝒢_j     k > 1
            
            const     k = 1
        
        as n→∞. Here, 𝒢_1,...,𝒢_k-1 are centered Gaussian variables with finite, non-zero variances, and δ_d := δ(d_1 = ... = d_k = d).
        
        The case k=0 is trivial. Let k≥ 1. By Eq. <ref>, it holds that:
        
            n^k-1𝒯^l,i,d_n,i,d
             
            =    n^k-1⟨ q^l_1-1_i,d , q^l_1-1_i_1,d_1⟩⟨ C_i_k,d_k^l_k→ L , C_i,d^l_k→ L⟩/n·∏_j=1^k-1⟨C_i_j,d_j^l_j→ L⊗ q_i_j+1,d_j+1^l_j+1-1/n, C_i,d^l_j→ l_j+1⟩
            
            =   ⟨ q^l_1-1_i,d , q^l_1-1_i_1,d_1⟩⟨ C_i_k,d_k^l_k→ L , C_i,d^l_k→ L⟩/n·∏_j=1^k-1⟨ C_i_j,d_j^l_j→ L⊗ q_i_j+1,d_j+1^l_j+1-1, C_i,d^l_j→ l_j+1⟩
        
        Note that intermediate activations do not depend on the index d_j, and so we remove the dependency on d_j in the relevant terms. Next, by applying Lem. <ref>,
        
            ⟨ q^l_1-1_i , q^l_1-1_i_1⟩⟨ C_i_k,d_k^l_k→ L , C_i,d^l_k→ L⟩/nd⟶Σ_i,i_1^l_1-1(∏_j=l_k^LΣ̇_i,i_k^l_j)δ_d
        
        Expanding the second term using Eq. <ref>:
        
            ⟨ C_i_j,d_j^l_j→ L⊗ q_i_j+1^l_j+1-1, C_i,d^l_j→ i_j+1⟩
               = C_i_j,d_j^l_j→ LC_i^l_j→ i_j+1q_i_j+1^l_j+1-1
               = C_i_j,d_j^l_j+1-1→ LP_i_j^l_j→ l_j+1-1 (P_i^l_j→ l_j+1-1)^⊤√(2)· Z_i^l_j+1-1q_i_j+1^l_j+1-1
               =√(2)·⟨ C_i_j,d_j^l_j+1-1→ L⊗ (Z_i^l_j+1-1q_i_j+1^l_j+1-1), P_i_j^l_j→ l_j+1-1 (P_i^l_j→ l_j+1-1)^⊤⟩
               =√(2)· C_i_j,d_j^l_j+1-1→ LP_i_j^l_j→ l_j+1-1 (P_i^l_j→ l_j+1-1)^⊤ Z_i^l_j+1-1q_i_j+1^l_j+1-1
               = ξ_j.
         
        The above expression is fully implementable in a Tensor Program (see <cit.>), and approaches a GP as width tend to infinity. In other words:
        
        
        
        
        
        
        
        
        
        
        
        
        
            ξ_j d⟶𝒢̂_j
        
        and denoting ξ = [ξ_1,...,ξ_k-1], and 𝒢̂ = [𝒢̂_1,...,𝒢̂_k-1], it holds using the multivariate Central Limit theorem:
        
            ξd⟶𝒢̂
        
        Using the Mann-Wald theorem <cit.> (where we take the mapping as the product pooling of ξ), we have that:
        
            ∏_j=1^k-1ξ_j d⟶∏_j=1^k-1𝒢̂_j
        
        Finally, by Slutsky's theorem,
        
            n^k-1𝒯^l,i,d_n,i,dd⟶Σ_i,i_1^l_1-1(∏_j=l_k^LΣ̇_i,i_k^l_j) ·(∏^k-1_j=1𝒢̂_j ) ·δ_d.
        
        Assigning 𝒢_j:= Σ_i,i_1^l_1-1(∏_j=l_k^LΣ̇_i,i_k^l_j)𝒢̂_j completes the proof.
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
         §.§ Proof of Thm. <ref>
        
        
        Since we assume that g is a finite neural network, i.e., m_l < ∞ for all l ∈ [H], throughout the proofs with no loss of generality we assume that m_1=…=m_H=1. 
        
        
        Let h(u;w) = g(z;f(x;w)) be a hypernetwork. We have: 
        
            𝒦^(r)_i,j =   ∑_α_1 + … + α_H = r
            α_1,…,α_H≥ 0r!/α_1 ! ⋯α_H! · z_i ·[ ∏^H-1_j=1ϕ̇(g^j_i) ]  ·∏^H_d=1⟨∇^(α_d)_w f^d_i , (∇_w h_j )^α_d⟩
        
        By the higher order product rule and the fact that the second derivative of a piece-wise linear function is 0 everywhere: 
        
            ∇^(r)_w h_i = ∑_α_1 + … + α_H = r
            α_1,…,α_H≥ 0r!/α_1 ! ⋯α_H! · z_i ·∇^(α_H)_w f^H_i ⊗^H-1_d=1 D_H-d
        
        where 
        
            D_d := ϕ̇(g^d_i) ·∇^(α_d)_w f^d_i
        
        In addition, by elementary tensor algebra, we have:
        
            𝒦^(r)_i,j =   ⟨∇^(r)_w h_i ,(∇_w h_j )^r ⟩
            
            =   ∑_α_1 + … + α_H = r
            α_1,…,α_H≥ 0r!/α_1 ! ⋯α_H!   z_i ·⟨∇^(α_H)_w f^H_i ·⊗^H-1_d=1 D_H-d , (∇_w h_j )^r ⟩
            
            =   ∑_α_1 + … + α_H = r
            α_1,…,α_H≥ 0r!/α_1 ! ⋯α_H!  z_i·⟨∇^(α_H)_w f^H_i, (∇_w h_j )^α_H⟩
                           ·∏^H-1_d=1⟨ϕ̇(g^H-d_i) ·∇^(α_H-d)_w f^H-d_i , (∇_w h_j )^α_H-d⟩
            
            =   ∑_α_1 + … + α_H = r
            α_1,…,α_H≥ 0r!/α_1 ! ⋯α_H! · z_i ·[ ∏^H-1_d=1ϕ̇(g^d_i) ] ·∏^H_d=1⟨∇^(α_d)_w f^d_i , (∇_w h_j )^α_d⟩
        
        Let h(u;w) = g(z;f(x;w)) be a hypernetwork. In addition, let,
        
            ∀ d∈ [H]:  h^d_j := 
            a^d-1_j ∏^H-d_t=1 f^H-t+1_j ·ϕ̇(g^H-t_j )
        
        We have:
        
            ⟨∇^(α_d)f^d_i , (∇_w h_j)^α_d⟩ = ∑_l∈ [L]^α_d∑_d∈ [H]^α_d(∏^α_d_k=1 h^d_k_j) ·𝒯^l,d_n,i,j,d
        
        We have:
        
            ⟨∇^(α_d)f^d_i , (∇_w h_j )^α_d⟩ = ∑_l∈ [L]^α_d⟨∂^α_d f^d_i /∂ W^l_1…∂ W^l_α_d, ⊗^α_d_k=1∂ h_j /∂ W^l_k⟩
        
        By the product rule:
        
            ∂ h_j /∂ W^l_k 
            =   ∑^H_d=1[∏^H-d_t=1 f^H-t+1_j ·ϕ̇(g^H-t_j) ] ·∂ f^d_j /∂ W^l_k· a^d-1_j =∑^H_d=1 h^d_j  ·∂ f^d_j/∂ W^l_k
        
        
        
        Hence, 
        
            ⊗^α_d_k=1∂ h_j/∂ W^l_k
            = ∑_d∈ [H]^α_d(∏^α_d_k=1 h^d_k_j)  ⊗^α_d_k=1∂ f^d_k_j/∂ W^l_k
        
        In particular, 
        
            ⟨∇^(α_d)f^d_i , (∇_w h_j)^α_d⟩ = ∑_l∈ [L]^α_d∑_d∈ [H]^α_d(∏^α_d_k=1 h^d_k_j) ·𝒯^l,d_n,i,j,d
        
        
        
        
        
        
        
        
        
        *
        
        
        Throughout the proof, in order to derive certain limits of various sequences of random variables, we implicitly make use of the Mann-Wald theorem <cit.>. For simplicity, oftentimes, we will avoid explicitly stating when this theorem is applied. As a general note, the repeated argument is as follows: terms, such as, n^max(α_d-1,0)·𝒯^l,d_n,i,j,d, 𝒬^d_n,j, g^d_i, etc', (see below) can be expressed as continuous mappings of jointly convergent random variables. Hence, they jointly converge, and continuous mappings over them converge as well. 
        
        By Lems. <ref> and <ref>, we have: 
        
            𝒦^(r)_i,j 
            = ∑_α_1 + … + α_H = r
            α_1,…,α_H≥ 0r!/α_1 ! ⋯α_H! · z_i ·[∏^H-1_d=1ϕ̇(g^d_i) ]  ·∏^H_d=1∑_l∈ [H]^α_d∑_d∈ [H]^α_d𝒬^d_n,j·𝒯^l,d_n,i,j,d
        
        where 𝒬^d_n,j := (∏^α_d_k=1 h^d_k_j). By the Mann-Wald theorem <cit.>, g^d_i converges to some random variable 𝒰^d_i. If ϕ̇ is a continuous function, then ϕ̇(g^d_i) converges to ϕ̇(𝒰^d_i). If ϕ is the ReLU activation function, by Lem. <ref>, ϕ̇(g^d_i) = (g^d_i) converges to (𝒰^d_i) in distribution. We notice that 𝒬^d_n,j converges in distribution to some random variable 𝒬^d_j. 
        
        
        
        The proof is divided into two cases: H=1 and H>1.
        
        
        
          
        Case H=1: First, we note that for H=1 and d∈ [H] (i.e., d=1), we have: 
        
            h^d_j    = a^d-1_j ·∏^H-d_t=1 f^H-t+1_j ·σ̇(g^H-t_j) = a^0_j = z_j
        
        In addition, ∏^H-1_d=1σ̇(g^d_i) = 1 as it is an empty product. Therefore, we can rewrite:
        
            𝒦^(r)_i,j = z_i · z^r_j ∑_l∈ [H]^r∑_d∈ [H]^r𝒯^l,d_n,i,j,d
        
        By Lem. <ref>, for r=1, the above tends to a constant as n→∞. For r>1, n^r-1·𝒯^l,d_n,i,j,d converges in distribution to zero for all d≠ (d,…,d) and converges to a non-constant random variable 𝒯^l_i,j,d otherwise. Hence, by the Mann-Wald theorem <cit.>,
        
            n^r-1·𝒦^(r)_i,jd⟶ z_i · z^r_j ∑_l∈ [H]^r𝒯^l_i,j,d
         
        which is a non-zero random variable. 
        
          
        Case H>1: 
        
        By Lem. <ref>, n^α_d-1·𝒯^l,d_n,i,j,d converges in distribution to zero for all d≠ (d,…,d). Therefore, in these cases, by Slutsky's theorem, n^α_d-1·𝒬^d_n,j·𝒯^l,d_n,i,j,d converges to zero in distribution. On the other hand, for each l∈ [H]^α_d, d∈ [H] and d = (d,…,d), by Lem. <ref>, we have: 
        
            n^α_d-1·𝒬^d_n,j·𝒯^l_n,i,j,dd⟶𝒬^d_j ·𝒯^l_i,j,d
        
        In particular, 
        
            n^max(α_d-1,0)∑_l∈ [H]^α_d∑_d∈ [H]^α_d·𝒬^d_n,j·𝒯^l_n,i,j,dd⟶∑_l∈ [H]^α_d∑_d ∈ [H]𝒬^d_j ·𝒯^l_i,j,d
        
        Consider the case where r ≥ H. In this case, for any α_1,…,α_H, such that, there are t>1 indices i ∈ [H], such that, α_i=0. The following random variable converges in distribution:
        
            X_n := n^r-(H-t)·∏^H_d=1∑_l∈ [H]^α_d∑_d∈ [H]^α_d𝒬^d_n,j·𝒯^l,d_n,i,j,d
        
        Therefore, by Slutsky's theorem:
        
            n^r-H·∏^H_d=1∑_l∈ [H]^α_d∑_d∈ [H]^α_d𝒬^d_n,j·𝒯^l,d_n,i,j,d  = n^-t· X_n d⟶ 0
        
        We have:
        
        
            n^r-H·⟨∇^(r)_w h_i , (∇_w h_j )^r⟩
            
            =    n^r-H∑_α_1 + … + α_H = r
            α_1,…,α_H≥ 0r!/α_1 ! ⋯α_H! · z_i ·[∏^H-1_d=1σ̇(g^d_i) ]  ·∏^H_d=1∑_l∈ [H]^α_d∑_d∈ [H]^α_d(∏^α_d_k=1 h^d_k_j)  𝒯^l,d_n,i,j,d
            
            =   ∑_α_1 + … + α_H = r
            α_1,…,α_H≥ 0r!/α_1 ! ⋯α_H! · z_i ·[∏^H-1_d=1σ̇(g^d_i) ]  ·∏^H_d=1 n^α_d-1∑_l∈ [H]^α_d∑_d∈ [H]^α_d𝒬^d_n,j·𝒯^l,d_n,i,j,d
            d⟶   ∑_α_1 + … + α_H = r
            α_1,…,α_H≥ 1r!/α_1 ! ⋯α_H! · z_i ·[∏^H-1_d=1(𝒰^d_i) ]  ·∏^H_d=1∑_l∈ [H]^α_d𝒬^d_j ·𝒯^l_i,j,d
        
        which is a non-constant random variable.
        
        Next, we consider the case when r ≤ H. 
        By Lem. <ref>, for any α_d ≥ 2, the term 𝒯^l,d_n,i,j,d tends to zero as n→∞. In addition, 𝒬^d_n,j converges in distribution. Therefore, for any α_d ≥ 2, we have:
        
            ∑_l∈ [L]^α_d∑_d∈ [H]^α_d𝒬^d_n,j·𝒯^l,d_n,i,j,dd⟶ 0
        
        Hence, for any α_1,…,α_H ≥ 0, such that, there is at least one α_d ≥ 2, we have:
        
            ∏^H_d=1∑_l∈ [L]^α_d∑_d∈ [H]^α_d𝒬^d_n,j·𝒯^l,d_n,i,j,dd⟶ 0
        
        
        
        
        
        
        
        
        On the other hand, for any 0≤α_1,…,α_H ≤ 1, the terms {𝒯^l,i,d_n,i,d}, {g^d_i} and {𝒬^d_n,j} converge jointly in distribution to some random variables {𝒯^l,i,d_i,d}, {(𝒰^d_i)} and {𝒬^d_j} as n →∞. Hence, 
        
            ⟨∇^(r)_w h_i, (∇_w h_j)^r⟩d⟶∑_α_1 + … + α_H = r
             0≤α_1,…,α_H ≤ 1 r! ·[∏^H-1_d=1(𝒰^d_i)]·∏^H_d=1∑_l∈ [L]^α_d∑_d∈ [H]^α_d𝒬^d_j·𝒯^l,d_i,j,d
        
        which is a non-constant random variable.
        
        
         §.§ Proofs of the Results in Sec. 4
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        *
        
         By <cit.>, taking the width n=min(n_1,...,n_L-1) to infinity, the outputs V^d(x;w) := f^d(x;w) are governed by a centered Gaussian process, such that, the entries V^d_i,j(x;w), given some input x, are independent and identically distributed. Moreover, it holds that:
        
            (V^d_i,j(x;w),V^d_i,j(x';w)) ∼𝒩(0,S^L(x,x')).
        
        with S^L(x,x') as defined in Eq. <ref>.
        For the function h(u;w) = g(z;f(x;w)), it holds for the first layer:
        
            g^1(z;f(x;w)) = √(1/m_0)V^1(x;w) z
        
        After taking the limit n=min(n_1,...,n_L-1) to infinity, the implicit network g is fed with Gaussian distributed weights.
        
        
        
        
        And so g^1(z;f(x;w)) also converges to a Gaussian process, such that:
        
            (g^1(z;f(x;w))_i,g^1(z';f(x';w))_i) ∼𝒩(0,Λ^1)
        
        where:
        
            Λ^1 = 1/m_0[    S^L(x,x) z^⊤ z   S^L(x',x)z'^⊤ z;  S^L(x,x') z^⊤ z' S^L(x',x')z'^⊤ z' ]
        
        In a similar fashion to the standard feed forward case, the pre-activations g^l(z;f(x;w)) converge to Gaussian processes as we let m = min(m_1,...,m_H-1) tend to infinity, with a covariance defined recursively:
        
            Σ^l(u,u') =  √(2)_(u,v)∼𝒩(0,Λ^l)[σ(u)σ(v)]
        
        where,
        
            Λ^l = 
            [     S^L(x,x)·Σ^l-1(u,u)   S^L(x',x)·Σ^l-1(u',u);   S^L(x,x')·Σ^l-1(u,u') S^L(x',x')·Σ^l-1(u',u') ]
        
        and
        
            Σ^0(z,z') = 1/m_0z^⊤ z'
        
        proving the claim.
        
        
        
        
        
        
        
        *
        
        
        We prove that Λ^l(u,u') is a function of S^0(x,x') and Σ^0(u,u') by induction. First, we note that Λ^1(u,u') is a function of S^L(x,x') and Σ^0(u,u') by definition. By the recursive definition of S^L(x,x'), it is a function of S^0(x,x'). Therefore, Λ^1(u,u') can be simply represented as a function of S^0(x,x') and Σ^0(u,u'). We assume by induction that Λ^l(u,u') is a function of S^0(x,x') and Σ^0(u,u'). We would like to show that Λ^l+1(u,u') is a function of S^0(x,x') and Σ^0(u,u'). By definition, Λ^l+1(u,u') is a function of S^L(x,x') and Σ^l(u,u'). In addition, Σ^l(u,u') is a function of Λ^l(u,u'). Hence, by induction, Σ^l(u,u') is simply a function of S^0(x,x') and Σ^0(u,u'). Since S^L(x,x') is a function of S^0(x,x'), we conclude that one can represent Λ^l+1(u,u') as a function of S^0(x,x') and Σ^0(u,u').
        
        
        *
        
        
        We note that: 
        
            Σ^0(p(z),p(z')) = 1/k p(z)^⊤ p(z) = 1/k∑^k_i=1cos(W^1_i z + b^1_i) cos(W^1_i z' + b^1_i)
        
        By Thm. 1 in <cit.>, we have:
        
            lim_k→∞Σ^0(p(z),p(z')) = exp[-z-z'^2_2/2]/2
        
        which is a function of exp[z-z'^2_2] as desired.
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        We make use of the following lemma in the proof of Thm. <ref>.
        
        
        Recall the parametrization of the implicit network:
        
            g^l_i := g^l(z_i;v) = √(1/m_l-1) f^l(x_i;w)· a_i^l-1
            
            a^l_i := a^l(z_i;v) = √(2)·σ(g^l_i)
             and  a^0_i := z_i
        
        For any pair u_i = {u_i}, we denote:
        
            P_i^l_1→ l_2 = ∏_l = l_1^l_2-1(√(2/m_l)V^l+1(x_i;w)· Z^l(z_i) )  and  Z^l(z) = diag (σ̇(g^l(z)))
        
        It holds that:
        
          * P_i^l_1 → l_2(P_j^l_1 → l_2)^⊤d⟶∏_l=l_1^l_2-1Σ̇^l(u_i,u_j)I.
            
          * ∂ h(u_i,w)/∂ v·∂^⊤ h(u_j,w)/∂ vd⟶∑_l=0^H-1(Σ^l(u_i,u_j)∏_h=l+1^H-1Σ̇^l(u_i,u_j)).
        
        where the limits are taken with respect to m,n →∞ sequentially.
        
        We have:
        
            P_i^l_1 → l_2(P_j^l_1 → l_2)^⊤
            
            =    P_i^l_1 → l_2-12/m_l_2-1V^l_2(x_i;w) 
            · Z^l_2-1(z_i)Z^l_2-1(z_j)V^l_2(x_j;w)^⊤ (P_j^l_1 → l_2-1)^⊤
        
        Note that it holds that when m,n→∞ sequentially, we have:
        
            2/m_l_2-1V^l_2(x_i;w)· Z^l_2-1(z_i)Z^l_2-1(z_j)V^l_2(x_j;w)^⊤
            d⟶   √(2)_(u,v)∼𝒩(0,Λ^l_2)[σ̇(̇u̇)̇σ̇(̇v̇)̇]I = Σ̇^l_2(u_i,u_j)I
        
        Applying the above recursively proves the first claim.
        Using the first claim, along with the derivation of the neural tangent kernel (see <cit.>) proves the second claim.
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        *
        
        
        
        
        Recalling that v = vec(f(x)) = [vec(V^1),...,vec(V^H)], concatenated into a single vector of length ∑_l=0^H-1m_l · m_l+1. The components of the inner matrix 𝒦^f(x,x') are given by:
        
            𝒦^f(x,x')(i,j) = ∑_l=1^L ⟨∂ v_i(x)/∂ w^l,∂ v_j(x')/∂ w^l⟩
        
        and it holds that in the infinite width limit, 𝒦^f(x,x') is a diagonal matrix:
        
            𝒦^f(x,x') d⟶Θ^f(x,x')· I
        
        By letting the widths n and m tend to infinity consecutively, by Lem. <ref>, it follows that:
        
            ∂ h(u;w)/∂ v·∂^⊤ h(u';w)/∂ vd⟶Θ^g(u,u',S^L(x,x'))
        
        Since 𝒦^f(x,x') = ∂ f(x;w)/∂ w·∂^⊤ f(x';w) /w converges to the diagonal matrix Θ^f(x,x')· I, the limit of 𝒦^h(u,u') is given by:
        
            𝒦^h(u,u') =   ∂ g(z;f(x;w))/∂ f(x;w)·∂ f(x;w)/∂ w·∂^⊤ f(x';w) /w·∂^⊤ g(z';f(x';w))/∂ f(x';w)
             
            =   ∂ h(u;w)/∂ v·∂ f(x;w)/∂ w·∂^⊤ f(x';w) /w·∂^⊤ h(u';w)/∂ v
            d⟶   Θ^f(x,x')·Θ^g(u,u',S^L(x,x'))
        
        where we used the results of Lem. <ref>.
        
        Next, we would like to prove that ∂𝒦^h(u,u')/∂ t|_t=0 = 0. For this purpose, we write the derivative explicitly:
        
            ∂𝒦^h(u,u')/∂ t = ∂ h(u;w)/∂ w·∂/∂ t∂^⊤ h(u';w)/∂ w 
            + ∂/∂ t∂ h(u;w)/∂ w·∂ ^⊤ h(u';w)/∂ w
        
        We notice that the two terms are the same up to changing between the inputs u and u'. Therefore, with no loss of generality, we can simply prove the convergence of the second term. We have:
        
            ∂/∂ t∂ h(u;w)/∂ w·∂ ^⊤ h(u';w)/∂ w
            
            =   [ ∂/∂ t(∂ h(u;w)/∂ f(x;w)·∂ f(x;w)/∂ w) ] ·∂ ^⊤ h(u';w)/∂ w
            
            =   [ ∂ h(u;w)/∂ f(x;w) ∂ t·∂ f(x;w)/∂ w + ∂ h(u;w)/∂ f(x;w)·∂ f(x;w)/∂ w ∂ t] ·∂ ^⊤ h(u';w)/∂ w
            
            =   ∂ h(u;w)/∂ f(x;w) ∂ t·∂ f(x;w)/∂ w·∂ ^⊤ h(u';w)/∂ w +  ∂ h(u;w)/∂ f(x;w)·∂ f(x;w)/∂ w ∂ t·∂ ^⊤ h(u';w)/∂ w
        
        We analyze each term separately. 
        
        
        
          
        Analyzing the first term By substituting ∂/∂ t = -μ∇_w c(w)∂^⊤/∂ w = -μ∇_w c(w)∂^⊤ f/∂ w∂^⊤/∂ f, we have:
        
            ∂ h(u;w)/∂ f(x;w) ∂ t·∂ f(x;w)/∂ w·∂ ^⊤ h(u';w)/∂ w
               =-μ∇_w c(w)∂^⊤ f(x;w)/∂ w·∂^2 h(u;w)/∂ f(x;w)∂ f(x;w)·∂ f(x;w)/∂ w·∂^⊤ f(x';w)/∂ w·∂^⊤ h(u';w)/∂ f(x';w)
               =-μ∇_w c(w)∂^⊤ f(x;w)/∂ w∂^2 h(u;w)/∂ f(x;w)∂ f(x;w)𝒦^f(x,x')·∂^⊤ h(u';w)/∂ f(x';w)
               =-μ∑^N_i=1∂ℓ(h(u_i;w),y_i)/∂ h(u_i;w)·∂ h(u_i;w)/∂ f(x;w)·𝒦^f(x,x_i)·∂^2 h(u;w)/∂ f(x;w)∂ f(x;w)·𝒦^f(x,x')·∂^⊤ h(u';w)/∂ f(x';w)
        
        It then follows:
        
            lim_n →∞∂ h(u;w)/∂ f(x;w) ∂ t·∂ f(x;w)/∂ w·∂^⊤ h(u';w)/∂ w
            
            =   -μ∑^N_i=1ℓ_i ·Θ^f(x,x_i)·Θ^f(x,x')lim_n →∞∂ h(u_i;w)/∂ f(x_i;w)·∂^2 h(u;w)/∂ f(x;w)∂ f(x;w)·∂ h(u';w)/∂ f(x';w)
        
        We notice that:
        
            lim_n→∞∂^⊤ h(u_i;w)/∂ f(x_i;w)·∂^2 h(u;w)/∂ f(x;w)∂ f(x;w)·∂ h(u';w)/∂ f(x';w)
            
            =   ∑_l_1,l_2lim_n →∞⟨∂^2 h(u;w)/∂ f^l_1(x;w)∂ f^l_2(x;w),  ∂ h(u_i;w)/∂ f^l_1(x_i;w)⊗∂ h(u';w)/∂ f^l_2(x';w)⟩
            
            :=   ∑_l_1,l_2𝒯^l_1,l_2_m(u,u_i,u')
        
        We recall that f^l(x;w) converges to a GP (as a function of x) as n →∞ <cit.>. Therefore, 𝒯^l_1,l_2_m(u,u_i,u') are special cases of the terms 𝒯^l,i,d_n,i,d (see Eq. <ref>) with weights that are distributed according to a GP instead of a normal distribution. In this case, we have: k=2, d=d_1=…=d_k=1, the neural network f^1 is replaced with h, the weights W^l are translated into f^l(x;w). We recall that the proof of Lem. <ref> showing that 𝒯^l,i,d_n,i,d = 𝒪_p(1/n^k-1) is simply based on Lem. <ref>. Since Lem. <ref> extends Lem. <ref> to our case, the proof of Lem. <ref> can be applied to show that 𝒯^l_1,l_2_m(u,u_i,u') ∼ 1/m.
        
        
        
          
        Analyzing the second term We would like to show that for any m>0, we have:
        
            ∂ h(u;w)/∂ f(x;w)·∂ f(x;w)/∂ w ∂ t·∂ ^⊤ h(u';w)/∂ wd⟶ 0
        
        as n→∞. Since ∂ w/∂ t = -μ∇_w c(w), we have:
        
            ∂ h(u;w)/∂ f(x;w)·∂ f(x;w)/∂ w ∂ t·∂ ^⊤ h(u';w)/∂ w
            
            =    -μ·∂ h(u;w)/∂ f(x;w)·∇_w c(w) ·∂^2 f(x;w)/∂ w^2·∂ ^⊤ h(u';w)/∂ w
            
            =    -μ·∂ h(u;w)/∂ f(x;w)·∇_w c(w) ·∂^2 f(x;w)/∂ w^2·∂^⊤ f(x';w)/∂ w·∂ ^⊤ h(u';w)/∂ f(x;w)
        
        In addition, we have:
        
            ∇_w c(w) = ∑^N_i=1∂ℓ(h(u_i;w),y_i)/∂ h(u_i;w)·∂ h(u_i;w)/∂ w
        
        We note that ∂ℓ(h(u_i;w),y_i)/∂ h(u_i;w) converges in distribution as m,n→∞. 
        Therefore, we can simply analyze the convergence of:
        
            ∑^N_i=1∂ h(u;w)/∂ f(x;w)·∂ h(u_i;w)/∂ w·∂^2 f(x;w)/∂ w^2·∂^⊤ f(x';w)/∂ w·∂ ^⊤ h(u';w)/∂ f(x;w)
        
        Since N is a constant, it is enough to show that each term converges to zero. We have:
        
            ∂ h(u;w)/∂ f(x;w)·∂ h(u_i;w)/∂ w·∂^2 f(x;w)/∂ w^2·∂^⊤ f(x';w)/∂ w·∂ ^⊤ h(u';w)/∂ f(x;w)
            
            =   ∂ h(u;w)/∂ f(x;w)·∂ h(u_i;w)/∂ f(x_i;w)·∂ f(x_i;w)/∂ w·∂^2 f(x;w)/∂ w^2·∂^⊤ f(x';w)/∂ w·∂ ^⊤ h(u';w)/∂ f(x;w)
            
            =   ∑_l,j,k∂ h(u;w)/∂ f(x;w)_l·∂ h(u_i;w)/∂ f(x_i;w)_j·∂ f(x_i;w)_j/∂ w·∂^2 f(x;w)_l/∂ w^2·∂^⊤ f(x';w)_k/∂ w·∂ ^⊤ h(u';w)/∂ f(x;w)_k
        
        where f(x;w)_j is the j'th output of f over x. In addition, the summation is done over the indices of the corresponding tensors. We note that for any m>0, the number of indices l,j,k is finite. We would like to show that each term in the sum tends to zero as n→∞. We can write:
        
            ∂ f(x_i;w)_j/∂ w·∂^2 f(x;w)_l/∂ w^2·∂^⊤ f(x';w)_k/∂ w 
            = ⟨∂^2 f(x;w)_l/∂ w^2,∂ f(x_i;w)_j/∂ w⊗∂ f(x';w)_k/∂ w⟩
         
        By Lem. <ref>, the term in Eq. <ref> tends to zero as n →∞. In addition, it is easy to see that ∂ h(u;w)/∂ f(x;w)_l, ∂ h(u_i;w)/∂ f(x_i;w)_j and ∂ ^⊤ h(u';w)/∂ f(x;w)_k converge to some random variables. Therefore, for any fixed m>0, the above sum converges to zero as n→∞.
