Neural Tangent Kernel: Convergence and Generalization in Neural Networks

As a network gets wider, gradient descent moves its outputs through a kernel that stays fixed, so an infinitely wide network trains like kernel regression.

The paper defines the neural tangent kernel, the dot product of two inputs' parameter gradients, which sets how much a gradient step on one example changes the prediction on another. In the infinite-width limit this kernel is the same for every random initialization and does not change during training, so under a squared-error loss the network's outputs follow a linear equation that can be solved exactly.

A 2018 theory paper from EPFL and Imperial College: two theorems about networks whose hidden layers grow without bound, checked on points of a circle and on 512 MNIST digits.

Arthur Jacot, Franck Gabriel and Clément Hongler study fully connected networks trained by gradient descent and ask what the trained network computes, which the non-convex loss surface had kept out of reach. Their answer holds in a limit, every hidden layer infinitely wide, and in that limit it is exact. The page covers, in order:

1. the tangent kernel: how one gradient step on one example moves predictions on every other input;
2. kernel gradient descent, the paper's way of describing training in function space;
3. the 1/√n parametrization and the two theorems: the kernel becomes deterministic, then stops moving;
4. the least-squares solution, which fits the kernel's top principal components first;
5. what a trained infinitely wide network converges to, and what the result leaves unexplained.

What happens to the function during training

A network with parameter vector $\theta$ computes a function $f_\theta$. Gradient descent moves $\theta$, but the quantity anyone cares about is the function: what the network outputs on a training example, and on inputs it has never seen. The map from $\theta$ to $f_\theta$ is highly nonlinear, so the loss as a function of $\theta$ has saddle points and many minima, and nothing guarantees that gradient descent finds a good one.

The same loss written as a function of $f$ is simple. Least squares on $N$ examples is a convex bowl in the $N$ numbers $f(x_1), \dots, f(x_N)$. The difficulty is that gradient descent cannot move those numbers directly; it can only move $\theta$, and each parameter step drags the whole function along with it, on training inputs and test inputs alike, by amounts the neural tangent kernel measures.

Earlier work describes the network before training. Neal (1996) showed that a one-hidden-layer network with random weights becomes a Gaussian process as its width goes to infinity: its output at any finite set of inputs is jointly Gaussian, with a covariance (a kernel) determined by the architecture. Lee et al. and Matthews et al. extended this to deep networks in 2018. This paper asks what happens to that random function once gradient descent starts moving it.

How the tangent kernel works

Take one gradient step on one training example $(x_j, y_j)$ with loss $\tfrac12 (f(x_j) - y_j)^2$. Write $r = f(x_j) - y_j$ for the residual. The step changes each parameter by $\Delta\theta_p = -\eta\, r\, \partial f(x_j)/\partial\theta_p$, where $\eta$ is the learning rate. To first order in $\eta$, the output at any other input $x$ changes by the sum over parameters of (how much $f(x)$ depends on $\theta_p$) times (how much $\theta_p$ moved):

where

is the neural tangent kernel: the dot product of the two inputs' gradient vectors, each of length $P$, the number of parameters. (For networks with several outputs it is an $n_L \times n_L$ matrix per pair of inputs; the paper writes it as a sum of tensor products $\partial_{\theta_p} F \otimes \partial_{\theta_p} F$. This page uses one output, as the paper's experiments do.) Read $\Theta(x, x_j)$ as how far the output at $x$ moves, per unit of learning rate and residual, when a step pushes the output at $x_j$. A large positive value means the two inputs share parameters that both outputs depend on in the same direction; a negative value means a push at $x_j$ pulls $x$ down.

A linear model makes this exact. Take $f(x) = \theta_1 + \theta_2 x$. The gradient with respect to the parameters is $(1, x)$, so $\Theta(x, x') = 1 + x x'$, the same for every $\theta$. Start at $\theta = (0, 0)$, train on the single point $x_j = 2$, $y_j = 3$, with $\eta = 0.1$. The residual is $r = -3$, the step is $\Delta\theta = 0.1 \cdot 3 \cdot (1, 2) = (0.3, 0.6)$, and the new function is $0.3 + 0.6x$. At $x = 2$ it rose by 1.5, which is $0.1 \cdot 3 \cdot \Theta(2, 2)$ with $\Theta(2,2) = 5$. At $x = -1$ it fell by 0.3, which is $0.1 \cdot 3 \cdot \Theta(-1, 2)$ with $\Theta(-1, 2) = -1$.

For a network, (1) holds only to first order, and the kernel itself depends on $\theta$: the gradients in (2) are evaluated at the current parameters, so every step changes them. Figure 1 runs this on a real network. Its inputs are points on the unit circle, $x = (\cos\gamma, \sin\gamma)$, so a function of $x$ can be drawn as a curve over the angle $\gamma$. The network has two ReLU hidden layers of width $n$ and one output; at $n = 100$ it has 10,501 parameters, at $n = 1000$, 1,005,001. Watch the bottom panel, where the actual change in the function after one gradient step is drawn over the prediction of (1) computed from the network's own kernel.

At the default setting (width 100, training input at 46°, target 0.5 above the current output, step size chosen so that (1) predicts landing exactly on the target), the network's own kernel gives $\Theta(x_j, x_j) = 0.401$, the learning rate is $\eta = 1/0.401 = 2.49$, and the output at $x_j$ moves from 0.479 to 0.971 against a predicted 0.979. Across the circle the largest gap between the actual and predicted change is 5.3% of the largest predicted change. At width 10 the gap is 22%; at width 1000 it is 1.1%. Two things are visible. The change has the shape of the kernel column $\Theta(\cdot, x_j)$, a bump centred on the training input that decays with angular distance, so one example moves the network's output on nearby inputs most and on distant ones much less. And the first-order picture gets more accurate as the width grows, at a fixed step size. Pushing the step slider to 2 doubles the step and shows the second-order error growing with it.

Kernel gradient descent

Summing (1) over all training examples and letting the learning rate go to zero turns gradient descent into a differential equation in training time $t$ (with learning rate $\eta$, step $k$ corresponds to $t = \eta k$; the paper's experiments use $\eta = 1$, so $t$ counts steps). For a cost $C$ that depends on the function only through its values at the training inputs, write $d(x_j)$ for $N$ times the derivative of the cost with respect to the output $f(x_j)$ (for a cost that averages over the $N$ examples, this cancels the average's $1/N$). The paper shows the network function then evolves as

and calls the right-hand side, with any kernel $K$ in place of $\Theta$, the kernel gradient $-\nabla_K C$. For least squares, $C(f) = \frac{1}{2N}\sum_j (f(x_j) - y_j)^2$ and $d_t(x_j) = f_t(x_j) - y_j$, the residual. Equation (3) says the training signal lives only on the $N$ training inputs, and the kernel spreads it to every other input. Predictions on test inputs are therefore whatever the kernel carries over from the training residuals.

Along this flow the cost falls at a rate set by the kernel:

If the $N \times N$ Gram matrix $\Theta(x_i, x_j)$ is positive definite, the right side is strictly negative until the cost gradient on the data is zero. A convex cost bounded below, such as least squares, then converges to its global minimum. The non-convexity in parameter space does not enter, provided the kernel stays positive definite while training runs. For a finite network the proviso is not automatic, because $\Theta$ changes with $\theta$ as training runs.

Section 3.1 of the paper gives a model where the kernel does not change. Draw $P$ random functions $f^{(1)}, \dots, f^{(P)}$ whose covariance is a target kernel $K$, and fit only their linear combination

(the random-features model of Rahimi and Recht, 2007). The gradient with respect to $\theta_p$ is $f^{(p)}(x)/\sqrt{P}$, independent of $\theta$, so the tangent kernel is $\tilde K(x, x') = \frac1P\sum_p f^{(p)}(x) f^{(p)}(x')$ for the whole of training. It is random, because the features are, but it is an average of $P$ independent terms, and by the law of large numbers it converges to $K$ as $P \to \infty$. Gradient descent on this model is kernel gradient descent with a kernel close to $K$. The paper's main theorems say a deep network behaves like this model as its width grows, although its features are not fixed and its parameters enter nonlinearly.

The 1/√n parametrization

The paper uses layers numbered 0 (input) to $L$ (output) with widths $n_0, \dots, n_L$, a nonlinearity $\sigma$ applied entrywise, and writes each layer as

with $\alpha^{(0)}(x) = x$ and output $f_\theta(x) = \tilde\alpha^{(L)}(x)$. Every weight and bias is initialized from $\mathcal N(0, 1)$. The usual parametrization (LeCun initialization) instead draws weights with variance $1/n_\ell$ and omits the $1/\sqrt{n_\ell}$ factor. The two give the same function at initialization and the same set of representable functions; they differ in the gradients. In (6), the derivative of the output with respect to a weight carries the factor $1/\sqrt{n_\ell}$, so a hidden layer with $n_\ell n_{\ell+1}$ weights contributes a sum of $n_\ell n_{\ell+1}$ terms each of size about $1/(n_\ell n_{\ell+1})$ to the kernel (2), which stays finite as the width grows. Under LeCun initialization the same layer's contribution grows in proportion to the width, so the kernel has no finite limit. The $1/\sqrt{P}$ in the random-features model (5) has the same role.

The paper notes a side effect: the factors put $1/\sqrt{n}$ into every weight's gradient but not into the biases', which reduces the weights' influence on training relative to the biases. The scalar $\beta$ scales the biases down to balance the two. The experiments use $\beta = 0.1$ and learning rate 1.0, which the paper equates to a width-100 network in the usual parametrization trained at learning rate 0.01. At width 100 a weight $W$ enters the layer as $W/\sqrt{n}$. Its gradient carries one factor $1/\sqrt{n}$, and dividing the updated $W$ by $\sqrt{n}$ applies a second, so the effective weight $W/\sqrt{n}$ moves by $1/n = 0.01$ times the step a usual-parametrization weight takes at the same learning rate.

Why the NTK stops changing as width grows

The paper's two theorems are statements about the kernel (2) of the network (6) as the hidden widths $n_1, \dots, n_{L-1}$ go to infinity. The first concerns initialization. Proposition 1 restates the Gaussian-process result in this parametrization: each output tends to a centred Gaussian process with covariance $\Sigma^{(L)}$, built layer by layer,

The expectation is over a pair of jointly Gaussian numbers $(f(x), f(x'))$ with the previous layer's covariance: wide layers sum many independent terms, so each pre-activation is Gaussian by the central limit theorem, and the next layer's covariance is the average of $\sigma$ applied to those Gaussians. Theorem 1 says the NTK converges in probability to a deterministic kernel $\Theta^{(L)}_\infty$ (times the identity across outputs), given by

with $\dot\sigma$ the derivative of the nonlinearity. Read (9) as an accounting of parameters. The term $\Sigma^{(\ell+1)}$ is the contribution of layer $\ell$'s own weights and biases: the gradient of $\tilde\alpha^{(\ell+1)}$ with respect to them is the previous activation, and the dot product of activations averages to the covariance. The term $\Theta^{(\ell)}_\infty \dot\Sigma^{(\ell+1)}$ carries the contributions of all earlier layers through one more nonlinearity, scaled by how often both inputs pass through it together ($\dot\Sigma$). The limit depends only on the depth, the nonlinearity, $\beta$ and the initialization variance.

For ReLU both expectations have closed forms (the arc-cosine kernels of Cho and Saul, 2009). For a centred Gaussian pair with variances $a, b$ and correlation $\cos\phi$,

Work it through for the paper's setting: $L = 4$, inputs on the unit circle ($n_0 = 2$), $\beta = 0.1$, and $x = x' = (1, 0)$, so $\phi = 0$. Layer 1: $\Sigma^{(1)} = 1/2 + 0.01 = 0.51$ and $\Theta^{(1)}_\infty = 0.51$. Each ReLU layer halves the variance ($\mathbb{E}[\sigma(u)^2] = a/2$) and has $\dot\Sigma = 1/2$, so $\Sigma^{(2)} = 0.265$ and $\Theta^{(2)}_\infty = 0.51 \cdot 0.5 + 0.265 = 0.52$; then $\Sigma^{(3)} = 0.1425$, $\Theta^{(3)}_\infty = 0.4025$; then $\Sigma^{(4)} = 0.08125$, $\Theta^{(4)}_\infty = 0.2825$. For an input at 90° the limit is 0.0927, and at 180°, 0.0644: pushing the output at $(1, 0)$ up moves the output at the opposite point of the circle up by less than a quarter as much.

Theorem 2 is about training. The parameters move along a training direction $g_t$ in function space (for gradient descent, $g_t = -d_t$, minus the cost derivative from (3)); as long as $\int_0^T \lVert g_t\rVert\,dt$ stays bounded as the width grows, which least squares guarantees because the residual only shrinks, the NTK stays equal to $\Theta^{(L)}_\infty$ uniformly on $[0, T]$. The mechanism, in Remark 4: each individual activation in a hidden layer changes by an amount that shrinks with width, so the gradients in (2) barely change, yet the many small changes add up to an order-one change in the output, which is what lets the lower layers contribute to learning at all. With a fixed kernel, (3) becomes a linear equation in $f$, and the network trains exactly as kernel gradient descent with $\Theta^{(L)}_\infty$.

Figure 2 re-runs the paper's first experiment: a depth-4 ReLU network, $\Theta(x_0, x)$ for $x_0 = (1, 0)$ and $x$ around the unit circle, ten random initializations per width, before and after 200 gradient steps at learning rate 1.0 on the target $f^*(x) = x_1 x_2$. Look at two things as the width grows: how far the ten t = 0 curves spread around the dashed limit, and how far the trained curves move away from them.

At initialization the ten kernels at angle 0 have a standard deviation of 0.170 at width 20, 0.095 at 100, 0.040 at 500 and 0.020 at 2000: every fourfold increase in width halves the spread, the $1/\sqrt{n}$ rate of an average of independent terms. The mean is close to the limit 0.2825 at every width from 100 up. Training moves the kernel in the direction the paper calls inflation: the mean value at angle 0 rises from 0.21 to 1.21 at width 20, by 86% at width 100, 20% at 500 and 4.4% at 2000. Over the same 200 steps the full parameter vector moves by 14.9%, 3.2%, 0.69% and 0.17% of its length, roughly in proportion to $1/n$, because the distance travelled stays about the same while the length of $\theta$ grows with the number of parameters. All four widths reach a similar training loss after 200 steps, 0.016 to 0.022.

Proposition 2 handles the positive definiteness that (4) needs. For a non-polynomial Lipschitz nonlinearity (one with bounded slope, such as ReLU or tanh) and depth $L \ge 2$, the limiting NTK restricted to the unit sphere is positive definite, so for distinct training inputs on the sphere the Gram matrix is invertible and least-squares training converges to zero training error. The paper motivates the sphere by high-dimensional data, where centred points tend to have nearly equal norms.

Least squares and early stopping

With the kernel fixed and the least-squares cost, (3) is a linear ordinary differential equation for the function. Define the operator $\Pi$ by $\Pi(g)(x) = \frac1N \sum_i \Theta(x, x_i)\, g(x_i)$, so (3) reads $\partial_t f_t = -\Pi(f_t - f^*)$. Its solution is a matrix exponential:

On the training inputs, $\Pi$ is the Gram matrix divided by $N$. Its eigenvectors are the kernel principal components of the data (kernel PCA, Schölkopf et al. 1998), and its eigenvalues $\lambda_i$ are the variance each component captures. Split the initial error $f_0 - f^*$ into a part $\Delta^0$ that the kernel cannot see (zero on every training input) and one part $\Delta^i$ along each component. Each component then decays at its own rate:

A component with eigenvalue $\lambda_i$ loses a factor $e$ of its error every $1/\lambda_i$ time units. Large-eigenvalue components are fit first, and the paper reads this as a reason for early stopping: stopping at a moderate $t$ fits the top components and leaves the small-eigenvalue ones, which for smooth kernels are the high-frequency, noise-like directions, mostly unfit.

Evenly spaced points on a circle make the components easy to name. The limiting NTK depends only on the angle between two inputs, so with $N = 16$ evenly spaced points the Gram matrix is circulant (each row is the previous row shifted by one), and its eigenvectors are the sines and cosines of frequency $k = 0, \dots, 8$ sampled at the 16 points. For the depth-4 limit with $\beta = 0.1$, the eigenvalues of $\Pi$ are 0.124 for the constant ($k = 0$), 0.0440 for $k = 1$, 0.0183 for $k = 2$, 0.00653 for $k = 3$, falling to 0.00158 at $k = 8$. At learning rate 1, an error at frequency 1 loses a factor $e$ every 23 steps; at frequency 8, every 633 steps.

Figure 3 runs (12) exactly for that kernel on a target $f^*(\gamma) = 0.6\cos\gamma + 0.3\sin 2\gamma$ with Gaussian label noise. It shows the mean network, starting from $f_0 = 0$; the random part of $f_0$ decays by the same formula and returns in the next section. Press play and watch the bars: frequencies 1 and 2 hold the signal, and by $t = 100$ their residuals are down to 1% and 16% of their starting size, while bars 5 to 8, which hold only noise, still have 79% to 85% of theirs.

With the default noise level of 0.2, the test error against the clean target bottoms out at 0.0159 near $t = 74$, when training error is still 0.029. After that the slow bars, which hold mostly noise, get fit too: training error goes to zero and test error climbs to 0.0388 by $t = 31{,}623$, 2.4 times the best value. With no noise the test error keeps falling to 0.0001 and stopping early only costs accuracy. Raising the noise to 0.5 moves the best stopping time earlier, to $t \approx 34$ for the first noise draw.

The paper tests (13) on real data. It takes 512 MNIST digits, approximates the limiting kernel with a width-10,000 network, and finds the top three eigenvalues 0.0457, 0.00108 and 0.00078. The first component is nearly the constant function, because the kernel PCA is not centred, which explains the factor of 42 between the first two eigenvalues; with $\beta = 1$ instead of 0.1 that gap is about ten times larger, which the paper reports makes training harder. It then sets the target to the initial network plus half the second component, $f^* = f_{\theta(0)} + 0.5 f^{(2)}$. In the limit the error decays as $0.5\,e^{-\lambda_2 t}$, halving every $\ln 2/0.00108 \approx 642$ steps, along a straight line in function space. Networks of width 100, 1000 and 10,000 deviate from that line less as the width grows, and the narrower ones converge faster than the wider ones, which the paper attributes to the kernel inflation of Figure 2: multiplying the kernel by a factor has the same effect as multiplying the learning rate.

What a trained wide network converges to

Letting $t \to \infty$ in (12) gives the trained network in closed form. Write $\Theta_{XX}$ for the $N \times N$ Gram matrix on the training inputs, $\Theta_{xX}$ for the row of kernel values between a test input $x$ and the training inputs, $y^*$ for the targets and $y_0$ for the initial network's outputs on the training inputs. With a positive-definite Gram matrix,

The first term is kernel regression with the NTK: it finds weights $c = \Theta_{XX}^{-1} y^*$ for which the sum of kernel bumps $\sum_j c_j \Theta(x, x_j)$ hits every training target exactly, and predicts with that sum. It equals kernel ridge regression with the ridge penalty sent to zero, and the posterior mean of a Gaussian process with prior covariance $\Theta$ conditioned on noise-free observations. The second term is the part of the random initial function that training never removes. It equals zero at every training input, because there $\Theta_{xX}\Theta_{XX}^{-1}$ picks out that input's own value and the bracket cancels, and it is nonzero in between. Since $f_0$ is a Gaussian process at infinite width and (14) is linear in it, the trained network is Gaussian too, with mean given by the first term and a variance that vanishes on the training data.

Figure 4 re-runs the paper's second experiment: depth 4, four training points on the circle, 1000 steps at learning rate 1.0, ten networks each at widths 50 and 1000. The shaded band is the 10th to 90th percentile of the exact infinite-width Gaussian; its variance follows Lee et al. (2019, eq. 16), with $\Sigma^{(4)}$ as the covariance of $f_0$.

Before training, the networks are random functions with a standard deviation of 0.285 at every input, the square root of $\Sigma^{(4)}(x, x) = 0.08125$. After training every network passes through all four points, where the band has zero width, and between the points the networks spread as the band predicts, at width 50 as well as at width 1000. The paper reports the same: even at width 50 the limit gives a good description of the distribution of trained networks. In our re-run the average of the ten networks differs from the limiting mean by 0.021 to 0.025 (root mean square over the circle), within the noise of averaging ten draws whose standard deviation is about 0.07.

The paper calls the first term of (14) the MAP (maximum a posteriori) estimate under a Gaussian prior with covariance $\Theta$. That describes the mean only. The random network the training started from has covariance $\Sigma$, not $\Theta$, and Lee et al. (2019) show that the distribution of the trained network is the posterior of no probabilistic model, unless $\Sigma$ and $\Theta$ happen to be proportional. Training only the last layer would give the posterior of the Gaussian process with kernel $\Sigma$; training all layers adds the lower-layer terms of (9) to the kernel without changing the prior.

What the NTK does not explain

Theorem 2 holds as the widths go to infinity and for a fixed training horizon $T$. The proof in the appendix takes the widths to infinity one layer at a time, first $n_1$, then $n_2$, and assumes a twice-differentiable nonlinearity with bounded second derivative, while the experiments use ReLU (the paper notes the assumptions "do not seem to be strictly needed"). Yang (2019) proved convergence with the widths growing at equal rates, and Arora, Du, Hu, Li, Salakhutdinov and Wang (2019) gave a bound for ReLU networks that only needs every width to be large and covers infinite training time.

A constant kernel means the network's internal representations effectively do not change: the hidden activations move by amounts that vanish with width. Figure 2 shows the consequence at finite width, where the narrower networks' kernels moved by 20% to several hundred percent in 200 steps. Chizat, Oyallon and Bach (2019) named the frozen-kernel regime lazy training and showed that it comes from the choice of scaling: any model whose output is scaled up enough relative to its parameter changes trains this way. They also found that convolutional networks pushed into the lazy regime perform worse. The best convolutional NTK of Arora et al. (2019), an 11-layer architecture with global average pooling, reaches 77.43% on CIFAR-10, while the ResNet-110 of 2015 reached 93.57%. The μP explainer covers a parametrization chosen so that features keep changing as width grows, the alternative limit to this one.

Within its scope the NTK turns several observations about wide networks into theorems. Gradient descent on an over-parameterized network reaches zero training error despite the non-convex loss, because the kernel stays positive definite. Wide networks can fit any labels, including random ones as in the Rethinking Generalization experiments, because kernel regression with an invertible Gram matrix interpolates any targets; the paper's introduction makes the same comparison with kernel methods. And their predictions on new inputs are set by the kernel, so comparing limiting NTKs compares how architectures, depths and initializations generalize in this regime. The paper also notes that Theorem 2 allows training directions that come from another network, as in a GAN's generator update.

Computing the NTK

The limit (7)-(11) for ReLU takes a few lines. The listing computes $\Sigma^{(L)}$ and $\Theta^{(L)}_\infty$ for one pair of inputs and reproduces the worked example above; the comment at the end gives the empirical kernel of a finite network, which Figure 1 computes layer by layer and an autograd library computes as a product of Jacobians.

import numpy as np

def relu_moments(a, b, c):
    """E[relu(u)relu(v)] and E[step(u)step(v)] for var a, b and cov c."""
    s = np.sqrt(a * b)
    th = np.arccos(np.clip(c / s, -1.0, 1.0))
    e = s * (np.sin(th) + (np.pi - th) * np.cos(th)) / (2 * np.pi)
    return e, (np.pi - th) / (2 * np.pi)

def ntk_limit(x, y, L, beta=0.1):
    """Sigma^(L) and Theta^(L)_inf for one pair of inputs (Prop. 1, Thm 1)."""
    n0 = len(x)
    S = x @ y / n0 + beta**2           # Sigma^(1)(x, y)
    Sx = x @ x / n0 + beta**2          # Sigma^(1)(x, x)
    Sy = y @ y / n0 + beta**2          # Sigma^(1)(y, y)
    T = S                              # Theta^(1) = Sigma^(1)
    for _ in range(L - 1):
        e, edot = relu_moments(Sx, Sy, S)
        Sx = relu_moments(Sx, Sx, Sx)[0] + beta**2
        Sy = relu_moments(Sy, Sy, Sy)[0] + beta**2
        S = e + beta**2                # Sigma^(l+1)
        T = T * edot + S               # Theta^(l+1) = Theta^(l) Sdot + Sigma
    return S, T

x0 = np.array([1.0, 0.0])
print(ntk_limit(x0, x0, L=4))          # (0.08125, 0.2825)

# Empirical NTK of a finite net: sum over layers of
#   (a_l(x) . a_l(y) / n_l + beta^2) * (delta_{l+1}(x) . delta_{l+1}(y))
# where a_l are activations and delta_{l+1} = df/dz_{l+1} from backprop.
# With autograd: theta = J(x) @ J(y).T, J = df/dparams (one row per input).

For $N$ training inputs the limiting Gram matrix costs $N^2$ evaluations of this recursion, and (14) needs one $N \times N$ solve, the usual cost of kernel regression: cheap for the 4 to 512 points of the paper's experiments. The solve takes on the order of $N^3$ operations, about $10^{14}$ for the 50,000 CIFAR-10 training images that Arora et al. used. The empirical NTK of a finite network costs a Jacobian of size $N \times P$, which for a 1,005,001-parameter network and 1,000 inputs is a billion numbers.