Gradient Flow

This blog is the first in a series in which we consider machine learning from four different viewpoints. We either use gradient descent or a fully Bayesian approach, and for each, we can choose to focus on either the network weights or the output function (figure 1). To make these four viewpoints easy to understand, we first consider the linear regression model for which all four viewpoints yield tractable solutions.

This blog concerns gradient descent for which we can write closed-form expressions for the evolution of the parameters and the function itself. By analyzing these expressions, we can gain insights into the trainability and convergence speed of the model. In part II and part III, we replace the linear model with a neural network. This leads to the development of the neural tangent kernel (NTK), which provides insights into the trainability and convergence speed of neural networks.

In part IV and part V of this series, we consider the Bayesian approach to model fitting. We show that for the linear regression model, it is possible to derive a closed-form expression for the posterior probability over the model parameters and marginalize over this distribution to make predictions. Alternatively, we can consider the output function as a Gaussian process, and make predictions using the conditional relation for normal distributions. In parts VI and part VII, we again replace the linear regression model with a neural network and this leads to Bayesian neural networks (Bayesian approach with parameters), and neural network Gaussian processes (Bayesian approach with output function).

*The code related to this blog post can be found on github here.

Linear Regression model

A linear regression model $\text{f}[\mathbf{x},\mathit{\varphi }]$ computes a single scalar output:

$$\begin{array}{}\text{(1)}& \text{f}[\mathbf{x},\mathit{\varphi }]={\mathbf{x}}^T\mathit{\varphi },\end{array}$$

where $\mathbf{x}\in {\mathbb{R}}^D$ is the input data vector and $\mathit{\varphi }\in {\mathbb{R}}^D$ are the parameters of model. We can extend this to an affine model by setting $\mathbf{x}\leftarrow [1{\mathbf{x}}^T{]}^T$ (i.e., inserting the value 1 at the start of the input vector), where now the first value of $\mathit{\varphi }$ controls the y-offset (bias) of the model. For convenience, and according to the usual convention in machine learning, we will still refer to this as a linear model.

Such a model can be trained using $I$ pairs $\{{\mathbf{x}}_i,y_i\}$ of training data and the parameters $\hat{\mathit{\varphi }}$ are typically chosen such that they minimize the least squares loss between the model predictions $\text{f}[{\mathbf{x}}_i,\mathit{\varphi }]$ and the ground truth target values $y_i$:

$$\begin{array}{rcl}\hat{\mathit{\varphi }}=\underset{\mathit{\varphi }}{argmin}[\text{L}[\mathit{\varphi }]]& =& \underset{\mathit{\varphi }}{argmin}[\frac{1}{2}\sum _{i=1}^I(\text{f}[{\mathbf{x}}_i,\mathit{\varphi }]-y_i{)}^2]\\ \text{(2)}& & =& \underset{\mathit{\varphi }}{argmin}[\frac{1}{2}\sum _{i=1}^I({\mathbf{x}}_i^T\mathit{\varphi }-y_i{)}^2],\end{array}$$

where we have substituted in the definition of the linear model in the second line. The factor 1/2 does not change the position of the minimum but simplifies subsequent expressions.

For convenience, we will store all of the training data vectors $\{{\mathbf{x}}_i\}$ in the columns of a matrix $\mathbf{X}\in {\mathbb{R}}^{D\times I}$ and the targets $\{y_i\}$ in a column vector $\mathbf{y}\in {\mathbb{R}}^{I\times 1}$ and then equivalently write:

$$\begin{array}{rcl}\hat{\mathit{\varphi }}& =& \underset{\mathit{\varphi }}{argmin}[\frac{1}{2}(\mathbf{f}[\mathbf{X},\mathit{\varphi }]–\mathbf{y}{)}^T(\mathbf{f}[\mathbf{X},\mathit{\varphi }]–\mathbf{y})]\\ \text{(3)}& & =& \underset{\mathit{\varphi }}{argmin}[\frac{1}{2}({\mathbf{X}}^T\mathit{\varphi }–\mathbf{y}{)}^T({\mathbf{X}}^T\mathit{\varphi }–\mathbf{y})].\end{array}$$

Here, the model $\mathbf{f}[\mathbf{X},\mathit{\varphi }]={\mathbf{X}}^T\mathit{\varphi }$ is applied to all of the data simultaneously and produces an $I\times 1$ column vector whose $i^{th}$ entry is $\text{f}[{\mathbf{x}}_i,\mathit{\varphi }]={\mathbf{x}}_i^T\mathit{\varphi }$. Figure 2 illustrates this loss function for the case where we are fitting a line to $I=3$ data points.

Gradient descent and gradient flow

In gradient descent, we randomly initialize the parameters of the model to values ${\mathit{\varphi }}_0$. We then measure the gradient of the loss function with respect to the parameters at this point and take a step of size $\alpha$ in the downhill (negative gradient) direction. By iterating these two steps, we hope to reach the point on the loss function with minimum loss.

More formally, when we perform gradient descent, we apply the update rule:

$$\begin{array}{}\text{(4)}& {\mathit{\varphi }}_{t+1}={\mathit{\varphi }}_t–\alpha \cdot \frac{\partial L}{\partial \mathit{\varphi }}\end{array}$$

where ${\mathit{\varphi }}_t$ represents the parameters at time $t$ and $\alpha$ is termed the learning rate. For the case of the linear regression model, the loss function is convex and so we are guaranteed to reach the global minimum eventually if the learning rate is sufficiently small.

Gradient flow

The above description is how we typically think about gradient descent in machine learning. However, in this blog we consider what happens when we use an infinitesimally small learning rate $\alpha$. Consider rearranging equation 4:

$$\begin{array}{}\text{(5)}& \frac{{\mathit{\varphi }}_{t+1}-{\mathit{\varphi }}_t}{\alpha }=-\frac{\partial L}{\partial \mathit{\varphi }},\end{array}$$

and then letting the learning rate $\alpha$ tend zero to get:

$$\begin{array}{}\text{(6)}& \frac{d\mathit{\varphi }}{dt}=-\frac{\partial L}{\partial \mathit{\varphi }}.\end{array}$$

This ordinary differential equation (ODE) is known as gradient flow and tells us how the parameters change over time We will see that for the linear regression model and for known initial parameters, we can solve this ODE to derive a closed-form expression for the parameter values at any given time during training.

Gradient flow for least squares loss

For the least squares loss function we can expand the right hand side of equation 6 to get the expression:

$$\begin{array}{rcl}\frac{d\mathit{\varphi }}{dt}& =& -\frac{\partial L}{\partial \mathit{\varphi }}\\ & =& -\frac{\partial }{\partial \mathit{\varphi }}[\frac{1}{2}(\mathbf{f}[\mathbf{X},\mathit{\varphi }]–\mathbf{y}{)}^T(\mathbf{f}[\mathbf{X},\mathit{\varphi }]–\mathbf{y})],\\ \text{(7)}& & =& -{\frac{\partial \mathbf{f}[\mathbf{X},\mathit{\varphi }]}{\partial \mathit{\varphi }}}^T(\mathbf{f}[\mathbf{X},\mathit{\varphi }]-\mathbf{y}),\end{array}$$

where we have used the chain rule on the quadratic loss function between lines two and three.

Evolution of residual

We now consider the residual vector (i.e., the differences between the model predictions $\mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]$ and the ground truth values $\mathbf{y}$). Just as equation 7 described how the parameters evolve, it is possible to write an ODE that describes how the residual vector $\mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]-\mathbf{y}$ changes over time.

$$\begin{array}{rcl}\frac{d}{dt}[\mathbf{f}[\mathbf{X},\mathit{\varphi }]-\mathbf{y}]& =& \frac{d\mathbf{f}[\mathbf{X},\mathit{\varphi }]}{dt}\\ & =& \frac{\partial \mathbf{f}[\mathbf{X},\mathit{\varphi }]}{\partial \mathit{\varphi }}\frac{d\mathit{\varphi }}{dt}\\ \text{(8)}& & =& -\left(\frac{\partial \mathbf{f}[\mathbf{X},\mathit{\varphi }]}{\partial \mathit{\varphi }}{\frac{\partial \mathbf{f}[\mathbf{X},\mathit{\varphi }]}{\partial \mathit{\varphi }}}^T\right)(\mathbf{f}[\mathbf{X},\mathit{\varphi }]-\mathbf{y}),\end{array}$$

where the first equality sign follows because only the first term $\mathbf{f}[\mathbf{X},\mathit{\varphi }]$ in the residual vector $\mathbf{f}[\mathbf{X},\mathit{\varphi }]-\mathbf{y}$ changes over time. The second line is an application of the chain rule, and the final expression follows from substituting in the result from equation 7.

Closed-form solution for function evolution

The ODE in equation 8 describes an exponential decay of the residual $\mathbf{f}[\mathbf{X},\mathit{\varphi }]-\mathbf{y}$. If the product of the gradients in the exponential term is constant, it has a simple solution (see this article) that is expressed in terms of the matrix exponential:

$$\begin{array}{}\text{(9)}& \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]-\mathbf{y}=\exp [-\frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]}{\partial \mathit{\varphi }}{\frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]}{\partial \mathit{\varphi }}}^T\cdot t](\mathbf{f}[\mathbf{X},{\mathit{\varphi }}_0]-\mathbf{y}),\end{array}$$

where $\mathbf{f}[\mathbf{X},{\mathit{\varphi }}_0]-\mathbf{y}$ is the initial residual vector that contains to the signed errors with the initial parameters ${\mathit{\varphi }}_0$. When $t=0$, it can be verified that both the left and right sides correspond to this initial vector.

Rearranging, this equation, we can get a closed-form solution for the evolution of the function $\mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]$ at the training points $\mathbf{X}$:

$$\begin{array}{}\text{(10)}& \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]& =& \mathbf{y}+\exp [-\frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]}{\partial \mathit{\varphi }}{\frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]}{\partial \mathit{\varphi }}}^T\cdot t](\mathbf{f}[\mathbf{X},{\mathit{\varphi }}_0]-\mathbf{y}).\end{array}$$

This tells us how the function at the data points evolves as gradient descent proceeds. By analyzing these expressions we can get information about the trainability of the model (i.e., the ability of the model to fit the training data exactly) as well as the convergence speed of the training process (i.e., how fast this happens).

Trainability and convergence in the linear model

Let’s examine the implications of this expression for the linear model for which:

$$\begin{array}{rcl}\mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]& =& {\mathbf{X}}^T{\mathit{\varphi }}_t\\ \text{(11)}& \frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]}{\partial {\mathit{\varphi }}_t}& =& {\mathbf{X}}^T.\end{array}$$

In this case the derivative does not depend on the current parameters ${\mathit{\varphi }}_t$. Substituting these expressions into equation 10 yields:

$$\begin{array}{}\text{(12)}& \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]& =& \mathbf{y}+\exp [-{\mathbf{X}}^T\mathbf{X}\cdot t]({\mathbf{X}}^T{\mathit{\varphi }}_0-\mathbf{y}).\end{array}$$

The evolution of this function is depicted in figure 3.

Trainability: We can draw conclusions about the trainability of the model by examining the exponential term in equation 12.

This aligns with our expectations. The data matrix $\mathbf{X}$ has size $(D+1)\times I$ where $D$ is the input data dimension (recall that we appended a 1 to each data example). Hence, the matrix ${\mathbf{X}}^T\mathbf{X}$ will be full rank if $I\le D+1$ (i.e., if the number of data points is less than or equal to the number of parameters). For the case, of fitting a line there are $D+1=2$ parameters and we can fit $I=1$ or $I=2$ points exactly. However, a line cannot pass exactly through $I>2$ data points in general position.

Convergence speed: Equation 12 can also tell us about the speed of convergence. The time taken to converge depends on the (i) the eigenvalues of ${\mathbf{X}}^T\mathbf{X}$ (which tell us how fast the function changes) and (ii) the magnitude of the initial error ${\mathbf{X}}^T{\mathit{\varphi }}_0-\mathbf{y}$ when projected onto the corresponding eigenvectors (which tells us how far we need to change). Discarding the second of these quantities (which depends on the particular initialization), we can loosely say that the time to converge depends on the smallest eigenvalue of ${\mathbf{X}}^T\mathbf{X}$.

These results are perhaps not surprising for the the linear regression model, but in part II of this series, we’ll apply the same ideas to gain insight into the trainability and convergence of neural networks.

Evolution of parameters

Let’s assume that the model is close to linear. In this case, it is well approximated by retaining only the linear term from a Taylor expansion around the initial parameters ${\mathit{\varphi }}_0$:

$$\begin{array}{}\text{(13)}& \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]\approx \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_0]+\frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_0]}{\partial \mathit{\varphi }}({\mathit{\varphi }}_t-{\mathit{\varphi }}_0)\end{array}$$

We replace the left hand side of equation 10 with this expression to yield:

$$\begin{array}{rcl}\text{(14)}& \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_0]+\frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_0]}{\partial \mathit{\varphi }}({\mathit{\varphi }}_t-{\mathit{\varphi }}_0)& & & & =\mathbf{y}+\exp [-\frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]}{\partial \mathit{\varphi }}{\frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]}{\partial \mathit{\varphi }}}^T\cdot t](\mathbf{f}[\mathbf{X},{\mathit{\varphi }}_0]-\mathbf{y}).\end{array}$$

Moving the first term to the right-hand side, we get:

$$\begin{array}{rcl}\text{(15)}& \frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_0]}{\partial \mathit{\varphi }}({\mathit{\varphi }}_t-{\mathit{\varphi }}_0)& & & & =-(\mathbf{I}-\exp [-\frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]}{\partial \mathit{\varphi }}{\frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]}{\partial \mathit{\varphi }}}^T\cdot t])(\mathbf{f}[\mathbf{X},{\mathit{\varphi }}_0]-\mathbf{y}).\end{array}$$

Finally, we can solve for ${\mathit{\varphi }}_t$:

$$\begin{array}{}\text{(16)}& {\mathit{\varphi }}_t& =& {\mathit{\varphi }}_0–{\frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_0]}{\partial \mathit{\varphi }}}^{+}(\mathbf{I}-\exp [-\frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]}{\partial \mathit{\varphi }}{\frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]}{\partial \mathit{\varphi }}}^T\cdot t])(\mathbf{f}[\mathbf{X},{\mathit{\varphi }}_0]-\mathbf{y}),\end{array}$$

where ${\mathbf{A}}^{+}$ denotes the Moore-Penrose pseudoinverse of $\mathbf{A}$. When the the system is overdetermined ($\mathbf{A}$ is portrait), this is given by ${\mathbf{A}}^{+}=({\mathbf{A}}^T\mathbf{A}{)}^{-1}{\mathbf{A}}^T$. For the underdetermined case ($\mathbf{A}$ is landscape), this is given by ${\mathbf{A}}^{+}={\mathbf{A}}^T(\mathbf{A}{\mathbf{A}}^T{)}^{-1}$.

Parameter evolution for linear model

Let’s examine parameter evolution for the linear model, where:

$$\begin{array}{rcl}\mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]& =& {\mathbf{X}}^T{\mathit{\varphi }}_t\\ \text{(17)}& \frac{\partial \mathbf{f}[\mathbf{X},{\mathit{\varphi }}_t]}{\partial {\mathit{\varphi }}_t}& =& {\mathbf{X}}^T.\end{array}$$

Here, the linear approximation is certainly valid (since the second derivatives are zero). Substituting these expressions into equation 16 for the overdetermined case, we get:

$$\begin{array}{}\text{(18)}& {\mathit{\varphi }}_t& =& {\mathit{\varphi }}_0-(\mathbf{X}{\mathbf{X}}^T{)}^{-1}\mathbf{X}(\mathbf{I}-\exp [-{\mathbf{X}}^T\mathbf{X}\cdot t])({\mathbf{X}}^T{\mathit{\varphi }}_0-\mathbf{y}).\end{array}$$

Figure 4 plots the evolution of the parameters for the fitting a line to the three points from figure 2.

Comparison to closed-form solution

For the linear model, it is possible to compute the least squares solution for the parameters in closed form (see end of blog for proof) using the relation:

$$\begin{array}{}\text{(19)}& \mathit{\varphi }=(\mathbf{X}{\mathbf{X}}^T{)}^{-1}\mathbf{X}\mathbf{y}.\end{array}$$

If the closed-form solution for the parameter evolution in equation 18 is correct, then we should retrieve this expression when $t\to \mathrm{\infty }$. To see that this is true, we re-arrange equation 18 into two terms:

$$\begin{array}{rcl}{\mathit{\varphi }}_t& =& {\mathit{\varphi }}_0-(\mathbf{X}{\mathbf{X}}^T{)}^{-1}\mathbf{X}(\mathbf{I}-\exp [-{\mathbf{X}}^T\mathbf{X}\cdot t])({\mathbf{X}}^T{\mathit{\varphi }}_0-\mathbf{y})\\ & =& {\mathit{\varphi }}_0-(\mathbf{X}{\mathbf{X}}^T{)}^{-1}(\mathbf{X}{\mathbf{X}}^T{\mathit{\varphi }}_0-\mathbf{X}\mathbf{y})+(\mathbf{X}{\mathbf{X}}^T{)}^{-1}\mathbf{X}\exp [-{\mathbf{X}}^T\mathbf{X}\cdot t]({\mathbf{X}}^T{\mathit{\varphi }}_0-\mathbf{y})\\ & =& {\mathit{\varphi }}_0-({\mathit{\varphi }}_0-(\mathbf{X}{\mathbf{X}}^T{)}^{-1})\mathbf{X}\mathbf{y})+(\mathbf{X}{\mathbf{X}}^T{)}^{-1}\mathbf{X}\exp [-{\mathbf{X}}^T\mathbf{X}\cdot t]({\mathbf{X}}^T{\mathit{\varphi }}_0-\mathbf{y})\\ \text{(20)}& & =& (\mathbf{X}{\mathbf{X}}^T{)}^{-1}\mathbf{X}\mathbf{y}+(\mathbf{X}{\mathbf{X}}^T{)}^{-1}\mathbf{X}\exp [-{\mathbf{X}}^T\mathbf{X}\cdot t]({\mathbf{X}}^T{\mathit{\varphi }}_0-\mathbf{y}),\end{array}$$

where we have separated out the two terms from the inner bracket in line two and simplified the first of these terms in the subsequent lines.

The first term is the same as the closed-form solution, so if the second term becomes zero as $t\to \mathrm{\infty }$, the two expressions agree. However, it’s not obvious this is the case. When there are more training points than parameters (the matrix $\mathbf{X}$ is landscape), $\mathbf{X}{\mathbf{X}}^T$ is full rank and a closed-form expression exists. However, this means that the term ${\mathbf{X}}^T\mathbf{X}$ in the exponential is not full rank and only part of the subspace tends to zero as $t\to \mathrm{\infty }$. Fortunately, the part that does not tend to zero is the null space of $\mathbf{X}$, which is removed by the pre-multiplication by $\mathbf{X}$. It follows that the term $\mathbf{X}\exp [-{\mathbf{X}}^T\mathbf{X}\cdot t]$ becomes zero when $t$ becomes large, and hence:

$$\begin{array}{}\text{(21)}& {\mathit{\varphi }}_{\mathrm{\infty }}& =& (\mathbf{X}{\mathbf{X}}^T{)}^{-1}\mathbf{X}\mathbf{y},\end{array}$$

as required.

Evolution of model predictions

Equation 18 provides a closed-form solution for the evolution of the parameters ${\mathit{\varphi }}_t$ over time. It follows that we can also get a closed-form solution for the evolution of the model predictions $y^{\ast }=\text{f}[{\mathbf{x}}^{\ast },{\mathit{\varphi }}_t]$ for a new points ${\mathbf{x}}^{\ast }$ by passing the solution for the parameters into the model.

For the linear model, the prediction for a new point ${\mathbf{x}}^{\ast }$ is given by $y^{\ast }={\mathbf{x}}^{\ast T}\mathit{\varphi }$, which gives the solution:

$$\begin{array}{}\text{(22)}& y_t^{\ast }& =& {\mathbf{x}}^{\ast T}{\mathit{\varphi }}_0-{\mathbf{x}}^{\ast T}(\mathbf{X}{\mathbf{X}}^T{)}^{-1}\mathbf{X}(\mathbf{I}-\exp [-{\mathbf{X}}^T\mathbf{X}\cdot t])({\mathbf{X}}^T{\mathit{\varphi }}_0-\mathbf{y}),\end{array}$$

and when $t\to \mathrm{\infty }$, the final predictions are given by:

$$\begin{array}{}\text{(23)}& {\hat{y}}^{\ast }={\mathbf{x}}^{\ast T}(\mathbf{X}{\mathbf{X}}^T{)}^{-1}\mathbf{X}\mathbf{y}.\end{array}$$

The evolution of the function for the linear model is illustrated in figure 5.

Evolution of prediction covariance

Finally, we note that the evolution of the predictions $y_t^{\ast }$ for a new point ${\mathbf{x}}^{\ast }$ can be written as an affine transformation of the initial parameters ${\varphi }_{\mathbf{0}}$:

$$\begin{array}{rcl}{y}_t^{\ast }& =& {\mathbf{x}}^{\ast T}{\mathit{\varphi }}_0-{\mathbf{x}}^{\ast T}(\mathbf{X}{\mathbf{X}}^T{)}^{-1}\mathbf{X}(\mathbf{I}-\exp [-{\mathbf{X}}^T\mathbf{X}\cdot t])({\mathbf{X}}^T{\mathit{\varphi }}_0-\mathbf{y})\\ & =& {\mathbf{x}}^{\ast T}(\mathbf{X}{\mathbf{X}}^T{)}^{-1}\mathbf{X}(\mathbf{I}-\exp [-{\mathbf{X}}^T\mathbf{X}\cdot t])\mathbf{y}\\ & & +[1,-{\mathbf{x}}^{\ast T}(\mathbf{X}{\mathbf{X}}^T{)}^{-1}\mathbf{X}(\mathbf{I}-\exp [-{\mathbf{X}}^T\mathbf{X}\cdot t])]\left[\begin{array}{c}{\mathbf{x}}^{T\ast }\\ {\mathbf{X}}^T\end{array}\right]{\mathit{\varphi }}_0\\ \text{(24)}& & =& a+\mathbf{b}{\varphi }_0,\end{array}$$

where we have denoted the two components of the affine transform by the scalar $a$ and the row vector $\mathbf{b}$.

Now let’s define a normally distributed prior over the initial parameters ${\mathit{\varphi }}_0$:

$$\begin{array}{}\text{(25)}& Pr({\mathit{\varphi }}_0)=\text{Norm}[\mathbf{0},{\sigma }_p^2\mathbf{I}],\end{array}$$

with mean zero and spherical covariance ${\sigma }_p^2\mathbf{I}$. The normally distributed prior implies a distribution over the initial model outputs. Because the prior distribution is normal and the transformation of the prior is affine, we can use the relation (see proof at end of blog) that if $Pr(\mathbf{z})={\text{Norm}}_{\mathbf{z}}[\mathit{\mu },\mathbf{\Sigma }]$, then the affine transform $\mathbf{a}+\mathbf{B}\mathbf{z}$ is distributed as:

$$\begin{array}{}\text{(26)}& Pr(\mathbf{a}+\mathbf{B}\mathbf{z})=\text{Norm}[\mathbf{a}+\mathbf{B}\mathit{\mu },\mathbf{B}\mathbf{\Sigma }{\mathbf{B}}^T].\end{array}$$

and so can compute the mean and variance of the distribution $Pr(y_t^{\ast })$ over the predictions as training proceeds. In practice, it is typical to add a value ${\sigma }^2$ to the variance to account for noise in training data. This is illustrated in figure 6.

*The code related to this blog post can be found on GitHub here.

Discussion

This blog has considered gradient descent training using gradient flow, using linear regression as a working example. The main results are that:

In the next part of this series, we’ll develop the neural tangent kernel by applying these same ideas to a neural network model.

*The code related to this blog post can be found on github here.

Least squares solution for linear model

For a linear model, it is possible to find the final parameters in closed form. Writing out the loss function $L[\mathit{\varphi }]$ explicitly, we have:

$$\begin{array}{rcl}L[\mathit{\varphi }]& =& ({\mathbf{X}}^T\mathit{\varphi }–\mathbf{y}{)}^T({\mathbf{X}}^T\mathit{\varphi }–\mathbf{y})\\ & =& {\mathit{\varphi }}^T\mathbf{X}{\mathbf{X}}^T\mathit{\varphi }–{\mathit{\varphi }}^T\mathbf{X}\mathbf{y}-{\mathbf{y}}^T{\mathbf{X}}^T\mathit{\varphi }–{\mathbf{y}}^T\mathbf{y}\\ & =& {\mathit{\varphi }}^T\mathbf{X}{\mathbf{X}}^T\mathit{\varphi }–2{\mathit{\varphi }}^T\mathbf{X}\mathbf{y}–{\mathbf{y}}^T\mathbf{y},\end{array}$$

where we have expanded out terms in the second line. We note that the second and third terms are both scalars and so are equal to their own transposes and these are combined into a single term in the third line.

We can find a closed-form solution for the best parameters by computing the derivative $dL/d\mathit{\varphi }$ of the scalar loss $L$ with respect to the parameter vector $\mathit{\varphi }$, setting the result to zero:

$$2\mathbf{X}{\mathbf{X}}^T\mathit{\varphi }–2\mathbf{X}\mathbf{y}=0$$

which can be rearranged to find a solution for the best parameters $\mathit{\varphi }$.

$$\hat{\mathit{\varphi }}=(\mathbf{X}{\mathbf{X}}^T{)}^{-1}\mathbf{X}\mathbf{y}.$$

This result can only be computed if there are at least as many columns $I$ in $\mathbf{X}$ as there are unknown parameters $D+1$ in $\mathit{\varphi }$. Otherwise, the matrix ${\mathbf{A}}^T\mathbf{A}$ will be singular. In cases like this (e.g., trying to fit a line with $D+1=2$ parameters to $I=1$ data points) the solution is ambiguous.

Affine transformations of normal distributions

The multivariate normal distribution is defined as:

$$Pr(\mathbf{x})=\frac{1}{(2\pi {)}^D|\mathbf{\Sigma }{|}^{1/2}}\exp [-0.5(\mathbf{x}-\mathit{\mu }){\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mathit{\mu })].$$

To transform the distribution as $\mathbf{y}=\mathbf{Ax}+\mathbf{b}$, we substitute in $\mathbf{x}={\mathbf{A}}^{-1}(\mathbf{y}-\mathbf{b})$ and divide by the determinant of the Jacobian of the transformation, which for this case is just $|\mathbf{A}|$ giving:

$$Pr(\mathbf{x})=\frac{1}{(2\pi {)}^D|\mathbf{\Sigma }{|}^{1/2}|\mathbf{A}|}\exp [-0.5({\mathbf{A}}^{-1}(\mathbf{y}-\mathbf{b})-\mathit{\mu }){\mathbf{\Sigma }}^{-1}({\mathbf{A}}^{-1}(\mathbf{y}-\mathbf{b})-\mathit{\mu })].$$

We’ll first work with the constant $\kappa$ from outside the exponential to show that it is the constant for a distribution with variance $\mathbf{A}\mathbf{\Sigma }{\mathbf{A}}^T$. We have:

$$\begin{array}{rcl}\kappa & =& \frac{1}{(2\pi {)}^D|\mathbf{\Sigma }{|}^{1/2}|\mathbf{A}|}\\ & =& \frac{1}{(2\pi {)}^D|\mathbf{\Sigma }{|}^{1/2}|\mathbf{A}{|}^{1/2}|\mathbf{A}{|}^{1/2}}\\ & =& \frac{1}{(2\pi {)}^D|\mathbf{\Sigma }{|}^{1/2}|\mathbf{A}{|}^{1/2}|{\mathbf{A}}^T{|}^{1/2}}\\ & =& \frac{1}{(2\pi {)}^D|\mathbf{A}{|}^{1/2}|\mathbf{\Sigma }{|}^{1/2}|{\mathbf{A}}^T{|}^{1/2}}\\ & =& \frac{1}{(2\pi {)}^D|\mathbf{A}\mathbf{\Sigma }{\mathbf{A}}^T{|}^{1/2}},\end{array}$$

which is exactly the required constant.

Now we’ll work with the quadratic in the exponential term to show that it corresponds to a normal distribution in $\mathbf{y}$ with variance $\mathbf{A}\mathbf{\Sigma }{\mathbf{A}}^T$ and mean $\mathbf{A}\mathit{\mu }+\mathbf{b}$.

$$\begin{array}{rcl}& & ({\mathbf{A}}^{-1}(\mathbf{y}-\mathbf{b})-\mathit{\mu }){\mathbf{\Sigma }}^{-1}({\mathbf{A}}^{-1}(\mathbf{y}-\mathbf{b})-\mathit{\mu })\\ & =& {\mathbf{y}}^T{\mathbf{A}}^{-T}{\mathbf{\Sigma }}^{-1}{\mathbf{A}}^{-1}\mathbf{y}-2({\mathbf{b}}^T{\mathbf{A}}^{-T}+{\mathit{\mu }}^T){\mathbf{\Sigma }}^{-1}{\mathbf{A}}^{-1}\mathbf{y}+({\mathbf{A}}^{-1}\mathbf{b}+\mathit{\mu }{)}^T{\mathbf{\Sigma }}^{-1}({\mathbf{A}}^{-1}\mathbf{b}+\mathit{\mu })\\ & =& {\mathbf{y}}^T{\mathbf{A}}^{-T}{\mathbf{\Sigma }}^{-1}{\mathbf{A}}^{-1}\mathbf{y}-2({\mathbf{b}}^T+{\mathit{\mu }}^T{\mathbf{A}}^T){\mathbf{A}}^{-T}{\mathbf{\Sigma }}^{-1}{\mathbf{A}}^{-1}\mathbf{y}\\ & & +({\mathbf{A}}^{-1}\mathbf{b}+\mathbf{A}\mathit{\mu }{)}^T{\mathbf{A}}^{-T}{\mathbf{\Sigma }}^{-1}{\mathbf{A}}^{-1}(\mathbf{b}+\mathbf{A}\mathit{\mu })\\ & =& {(\mathbf{y}-(\mathbf{A}\mathit{\mu }+\mathbf{b}))}^T{\left(\mathbf{A}\mathbf{\Sigma }{\mathbf{A}}^T\right)}^{-1}(\mathbf{y}-(\mathbf{A}\mathit{\mu }+\mathbf{b})).\end{array}$$

This is clearly the quadratic term from a normal distribution in $\mathbf{y}$ with variance $\mathbf{A}\mathbf{\Sigma }{\mathbf{A}}^T$ and mean $\mathbf{A}\mathit{\mu }+\mathbf{b}$ as required.

Thanks to Lechao Xiao who explained the evolution of the parameters to me.