Ridgeless kernel regression and kernel gradient descent

Kernel methods are machine learning algorithms which use only pairwise relations between input points. Kernel methods do not depend on the concrete values of the inputs; they only depend on the relations between the inputs and other inputs (such as the training set). These pairwise relations are fully captured by the kernel function: a symmetric, positive-semidefinite function of two inputs which represents some notion of similarity between the two inputs. A fully equivalent condition is that there exists some feature map ${\displaystyle {\mathbf {x}}\mapsto \psi ({\mathbf {x}})}$ such that the kernel function can be written as a dot product of the mapped inputs

$${\displaystyle K({\mathbf {x} },{\mathbf {x} }')=\psi ({\mathbf {x} })\cdot \psi ({\mathbf {x} }').}$$

The properties of a kernel method depend on the choice of kernel function. (Note that ${\displaystyle \psi ({\mathbf {x}})}$ may have higher dimension than ${\displaystyle \mathbf {x} }$.) As a relevant example, consider linear regression. This is the task of estimating ${\displaystyle {\mathbf {w}}^{*}}$ given ${\displaystyle N}$ samples ${\displaystyle ({\mathbf {x}}_{i},y_{i})}$ generated from ${\displaystyle y^{*}({\mathbf {x}})={\mathbf {w}}^{*}\cdot {\mathbf {x}}}$, where each ${\displaystyle \mathbf {x} _{i}}$ is drawn according to some input data distribution. In this setup, ${\displaystyle {\mathbf {w}}^{*}}$ is the weight vector which defines the true function ${\displaystyle y^{*}}$; we wish to use the training samples to develop a model ${\displaystyle \mathbf {\hat {w}} }$ which approximates ${\displaystyle {\mathbf {w}}^{*}}$. We do this by minimizing the mean-square error between our model and the training samples:

$${\displaystyle {\mathbf {\hat {w}} }=\arg \min _{\mathbf {w} }{\frac {1}{N}}\sum _{i=0}^{N}||y^{*}({\mathbf {x} }_{i})-{\mathbf {w} }\cdot {\mathbf {x} }_{i}||^{2}}$$

There exists an explicit solution for ${\displaystyle \mathbf {\hat {w}} }$ which minimizes the squared error: ${\displaystyle {\mathbf {\hat {w}}}=({\mathbf {X}}{\mathbf {X}}^{T})^{-1}{\mathbf {X}}{\mathbf {y}}}$, where ${\displaystyle {\mathbf {X}}}$ is the matrix whose columns are the training inputs, and ${\displaystyle {\mathbf {y}}}$ is the vector of training outputs. Then, the model can make predictions on new inputs: ${\displaystyle {\hat {y}}({\mathbf {x}})={\mathbf {\hat {w}}}\cdot {\mathbf {x}}}$. However, this result can be rewritten as: ${\displaystyle {\hat {y}}({\mathbf {x}})=({\mathbf {x}}^{T}{\mathbf {X}})({\mathbf {X}}^{T}{\mathbf {X}})^{-1}{\mathbf {y}}}$.^[7] Note that this dual solution is expressed solely in terms of the inner products between inputs. This motivates extending linear regression to settings in which, instead of directly taking inner products between inputs, we first transform the inputs according to a chosen feature map and then evaluate the inner products between the transformed inputs. As discussed above, this can be captured by a kernel function ${\displaystyle K({\mathbf {x}},{\mathbf {x}}')}$, since all kernel functions are inner products of feature-mapped inputs. This yields the ridgeless kernel regression estimator:

$${\displaystyle {\hat {y}}({\mathbf {x}})=K({\mathbf {x}},{\mathbf {X}})\;K({\mathbf {X}},{\mathbf {X}})^{-1}\;{\mathbf {y}}.}$$

If the kernel matrix ${\displaystyle K({\mathbf {X}},{\mathbf {X}})}$ is singular, one uses the Moore-Penrose pseudoinverse. The regression equations are called "ridgeless" because they lack a ridge regularization term.

In this view, linear regression is a special case of kernel regression with the identity feature map: ${\displaystyle \psi ({\mathbf {x}})={\mathbf {x}}}$. Equivalently, kernel regression is simply linear regression in the feature space (i.e. the range of the feature map defined by the chosen kernel). Note that kernel regression is typically a nonlinear regression in the input space, which is a major strength of the algorithm.

Just as it’s possible to perform linear regression using iterative optimization algorithms such as gradient descent, one can perform kernel regression using kernel gradient descent. This is equivalent to performing gradient descent in the feature space. It’s known that if the weight vector is initialized close to zero, least-squares gradient descent converges to the minimum-norm solution, i.e., the final weight vector has the minimum Euclidean norm of all the interpolating solutions. In the same way, kernel gradient descent yields the minimum-norm solution with respect to the RKHS norm. This is an example of the implicit regularization of gradient descent.

The NTK gives a rigorous connection between the inference performed by infinite-width ANNs and that performed by kernel methods: when the loss function is the least-squares loss, the inference performed by an ANN is in expectation equal to ridgeless kernel regression with respect to the NTK. This suggests that the performance of large ANNs in the NTK parametrization can be replicated by kernel methods for suitably chosen kernels.^[1][2]

Overparametrization, interpolation, and generalization

In overparametrized models, the number of tunable parameters exceeds the number of training samples. In this case, the model is able to memorize (perfectly fit) the training data. Therefore, overparametrized models interpolate the training data, achieving essentially zero training error.^[8]

Modern overparametrized models achieve low generalization error despite having the capacity to interpolate (memorize) the training set.^[8] This phenomenon can be understood by studying the generalization properties of high-dimensional kernel regression.

Kernel regression is typically viewed as a non-parametric learning algorithm, since there are no explicit parameters to tune once a kernel function has been chosen. An alternate view is to recall that kernel regression is simply linear regression in feature space, so the “effective” number of parameters is the dimension of the feature space. Therefore, studying kernels with high-dimensional feature maps can provide insights about strongly overparametrized models.

As an example, consider the problem of generalization. According to classical statistics, memorization should cause models to fit noisy signals in the training data, harming their performance on unseen data. To mitigate this, machine learning algorithms often introduce regularization to mitigate noise-fitting tendencies. Surprisingly, modern neural networks (which tend to be strongly overparametrized) seem to generalize well, even in the absence of explicit regularization.^[8][9] To study the generalization properties of overparametrized neural networks, one can exploit the infinite-width duality with ridgeless kernel regression. Recent works^[10][11][12] have derived equations describing the expected generalization error of high-dimensional kernel regression; these results immediately explain the generalization of sufficiently wide neural networks trained to convergence on least-squares.

Convergence to a global minimum

For a convex loss functional ${\displaystyle {\mathcal {C}}}$ with a global minimum, if the NTK remains positive-definite during training, the loss of the ANN ${\displaystyle {\mathcal {C}}\left(f\left(\cdot$ ;\theta \left(t\right)\right)\right)} converges to that minimum as ${\displaystyle t\to \infty }$. This positive-definiteness property has been shown in a number of cases, yielding the first proofs that large-width ANNs converge to global minima during training.^{[1][13][14][15][16][17]}

Case 1: Scalar output

An ANN with scalar output consists of a family of functions ${\displaystyle f\left(\cdot ,\theta \right):\mathbb {R} ^{n_{\mathrm {in} }}\to \mathbb {R} }$ parametrized by a vector of parameters ${\displaystyle \theta \in \mathbb {R} ^{P}}$.

The NTK is a kernel ${\displaystyle \Theta$ :\mathbb {R} ^{n_{\mathrm {in} }}\times \mathbb {R} ^{n_{\mathrm {in} }}\to \mathbb {R} } defined by

$${\displaystyle \Theta \left(x,y;\theta \right)=\sum _{p=1}^{P}\partial _{\theta _{p}}f\left(x;\theta \right)\partial _{\theta _{p}}f\left(y;\theta \right).}$$

In the language of kernel methods, the NTK ${\displaystyle \Theta }$ is the kernel associated with the feature map ${\displaystyle \left(x\mapsto \partial _{\theta _{p}}f\left(x;\theta \right)\right)_{p=1,\ldots ,P}}$. To see how this kernel drives the training dynamics of the ANN, consider a dataset ${\displaystyle \left(x_{i}\right)_{i=1,\ldots ,n}\subset \mathbb {R} ^{n_{\mathrm {in} }}}$ with scalar labels ${\displaystyle \left(z_{i}\right)_{i=1,\ldots ,n}\subset \mathbb {R} }$ and a loss function ${\displaystyle c:\mathbb {R} \times \mathbb {R} \to \mathbb {R} }$. Then the associated empirical loss, defined on functions ${\displaystyle f:\mathbb {R} ^{n_{\mathrm {in} }}\to \mathbb {R} }$, is given by

$${\displaystyle {\mathcal {C}}\left(f\right)=\sum _{i=1}^{n}c\left(f\left(x_{i}\right),z_{i}\right).}$$

When the ANN ${\displaystyle f\left(\cdot$ ;\theta \right):\mathbb {R} ^{n_{\mathrm {in} }}\to \mathbb {R} } is trained to fit the dataset (i.e. minimize ${\displaystyle {\mathcal {C}}}$) via continuous-time gradient descent, the parameters ${\displaystyle \left(\theta \left(t\right)\right)_{t\geq 0}}$ evolve through the ordinary differential equation:

$${\displaystyle \partial _{t}\theta \left(t\right)=-\nabla {\mathcal {C}}\left(f\left(\cdot$$ ;\theta \right)\right).}

During training the ANN output function follows an evolution differential equation given in terms of the NTK:

${\displaystyle \partial _{t}f\left(x;\theta \left(t\right)\right)=-\sum _{i=1}^{n}\Theta \left(x,x_{i};\theta \right)\partial _{w}c\left(w,z_{i}\right){\Big |}_{w=f\left(x_{i};\theta \left(t\right)\right)}.}$

This equation shows how the NTK drives the dynamics of ${\displaystyle f\left(\cdot$ ;\theta \left(t\right)\right)} in the space of functions ${\displaystyle \mathbb {R} ^{n_{\mathrm {in} }}\to \mathbb {R} }$ during training.

Case 2: Vector output

An ANN with vector output of size ${\displaystyle n_{\mathrm {out} }}$ consists in a family of functions ${\displaystyle f\left(\cdot$ ;\theta \right):\mathbb {R} ^{n_{\mathrm {in} }}\to \mathbb {R} ^{n_{\mathrm {out} }}} parametrized by a vector of parameters ${\displaystyle \theta \in \mathbb {R} ^{P}}$.

In this case, the NTK ${\displaystyle \Theta$ :\mathbb {R} ^{n_{\mathrm {in} }}\times \mathbb {R} ^{n_{\mathrm {in} }}\to {\mathcal {M}}_{n_{\mathrm {out} }}\left(\mathbb {R} \right)} is a matrix-valued kernel, with values in the space of ${\displaystyle n_{\mathrm {out} }\times n_{\mathrm {out} }}$ matrices, defined by

$${\displaystyle \Theta _{k,l}\left(x,y;\theta \right)=\sum _{p=1}^{P}\partial _{\theta _{p}}f_{k}\left(x;\theta \right)\partial _{\theta _{p}}f_{l}\left(y;\theta \right).}$$

Empirical risk minimization proceeds as in the scalar case, with the difference being that the loss function takes vector inputs ${\displaystyle c:\mathbb {R} ^{n_{\mathrm {out} }}\times \mathbb {R} ^{n_{\mathrm {out} }}\to \mathbb {R} }$. The training of ${\displaystyle f_{\theta \left(t\right)}}$ through continuous-time gradient descent yields the following evolution in function space driven by the NTK:

$${\displaystyle \partial _{t}f_{k}\left(x;\theta \left(t\right)\right)=-\sum _{i=1}^{n}\sum _{l=1}^{n_{\mathrm {out} }}\Theta _{k,l}\left(x,x_{i};\theta \right)\partial _{w_{l}}c\left(\left(w_{1},\ldots ,w_{n_{\mathrm {out} }}\right),z_{i}\right){\Big |}_{w=f\left(x_{i};\theta \left(t\right)\right)}.}$$

This generalizes the equation shown in case 1 for scalar outputs.

Interpretation

Each data point ${\displaystyle x_{i}}$ influences the evolution, of the output ${\displaystyle f\left(x;\theta \right)}$ for each input ${\displaystyle x}$, throughout the training. More concretely, with respect to example ${\displaystyle i}$, the NTK value ${\displaystyle \Theta \left(x,x_{i};\theta \right)}$ determines the influence of the loss gradient ${\displaystyle \partial _{w}c\left(w,z_{i}\right){\big |}_{w=f\left(x_{i};\theta \right)}}$ on the evolution of ANN output ${\displaystyle f\left(x;\theta \right)}$ through a gradient descent step. In the scalar case, this reads

$${\displaystyle f\left(x;\theta \left(t+\epsilon \right)\right)-f\left(x;\theta \left(t\right)\right)\approx \epsilon \sum _{i=1}^{n}\Theta \left(x,x_{i};\theta \left(t\right)\right)\partial _{w}c\left(w,z_{i}\right){\big |}_{w=f\left(x_{i};\theta \right)}.}$$

Wide fully-connected ANNs have a deterministic NTK, which remains constant throughout training

Consider an ANN with fully-connected layers ${\displaystyle \ell =0,\ldots ,L}$ of widths ${\displaystyle n_{0}=n_{\mathrm {in} },n_{1},\ldots ,n_{L}=n_{\mathrm {out} }}$, so that ${\displaystyle f\left(\cdot$ ;\theta \right)=R_{L-1}\circ \cdots \circ R_{0}} , where ${\displaystyle R_{\ell }=\sigma \circ A_{\ell }}$ is the composition of an affine transformation ${\displaystyle A_{i}}$ with the pointwise application of a nonlinearity ${\displaystyle \sigma$ :\mathbb {R} \to \mathbb {R} } , where ${\displaystyle \theta }$ parametrizes the maps ${\displaystyle A_{0},\ldots ,A_{L-1}}$. The parameters ${\displaystyle \theta \in \mathbb {R} ^{P}}$ are initialized randomly, in an independent, identically distributed way.

As the widths grow, the NTK's scale is affected by the exact parametrization of the ${\displaystyle A_{i}}$'s and by the parameter initialization. This motivates the so-called NTK parametrization ${\displaystyle A_{\ell }\left(x\right)={\frac {1}{\sqrt {n_{\ell }}}}W^{\left(\ell \right)}x+b^{\left(\ell \right)}}$. This parametrization ensures that if the parameters ${\displaystyle \theta \in \mathbb {R} ^{P}}$ are initialized as standard normal variables, the NTK has a finite nontrivial limit. In the large-width limit, the NTK converges to a deterministic (non-random) limit ${\displaystyle \Theta _{\infty }}$, which stays constant in time.

The NTK ${\displaystyle \Theta _{\infty }}$ is explicitly given by ${\displaystyle \Theta _{\infty }=\Theta ^{\left(L\right)}}$, where ${\displaystyle \Theta ^{\left(L\right)}}$ is determined by the set of recursive equations:

${\displaystyle {\begin{aligned}\Theta ^{\left(1\right)}\left(x,y\right)&=\Sigma ^{\left(1\right)}\left(x,y\right),\\\Sigma ^{\left(1\right)}\left(x,y\right)&={\frac {1}{n_{\mathrm {in} }}}x^{T}y+1,\\\Theta ^{\left(\ell +1\right)}\left(x,y\right)&=\Theta ^{\left(\ell \right)}\left(x,y\right){\dot {\Sigma }}^{\left(\ell +1\right)}\left(x,y\right)+\Sigma ^{\left(\ell +1\right)}\left(x,y\right),\\\Sigma ^{\left(\ell +1\right)}\left(x,y\right)&=L_{\Sigma ^{\left(\ell \right)}}^{\sigma }\left(x,y\right),\\{\dot {\Sigma }}^{\left(\ell +1\right)}\left(x,y\right)&=L_{\Sigma ^{\left(\ell \right)}}^{\dot {\sigma }},\end{aligned}}}$

where ${\displaystyle L_{K}^{f}}$ denotes the kernel defined in terms of the Gaussian expectation:

${\displaystyle L_{K}^{f}\left(x,y\right)=\mathbb {E} _{\left(X,Y\right)\sim {\mathcal {N}}\left(0,{\begin{pmatrix}K\left(x,x\right)&K\left(x,y\right)\\K\left(y,x\right)&K\left(y,y\right)\end{pmatrix}}\right)}\left[f\left(X\right)f\left(Y\right)\right].}$

In this formula the kernels ${\displaystyle \Sigma ^{\left(\ell \right)}}$ are the ANN's so-called activation kernels.^[27][28][5]

Wide fully connected networks are linear in their parameters throughout training

The NTK describes the evolution of neural networks under gradient descent in function space. Dual to this perspective is an understanding of how neural networks evolve in parameter space, since the NTK is defined in terms of the gradient of the ANN's outputs with respect to its parameters. In the infinite width limit, the connection between these two perspectives becomes especially interesting. The NTK remaining constant throughout training at large widths co-occurs with the ANN being well described throughout training by its first order Taylor expansion around its parameters at initialization:^[6]

${\displaystyle f\left(x;\theta (t)\right)=f\left(x;\theta (0)\right)+\nabla _{\theta }f\left(x;\theta (0)\right)\left(\theta (t)-\theta (0)\right)+{\mathcal {O}}\left(\min \left(n_{1}\dots n_{L-1}\right)^{-{\frac {1}{2}}}\right).}$