Deep Convolutional Networks as shallow Gaussian Processes

In plain words

A neural network is usually a single fixed function that gets trained. If you start it from random weights and make every layer very wide, something simpler happens: the network's output across different inputs behaves like a Gaussian process. A Gaussian process is a well-understood probability model that gives a prediction plus a measure of how sure it is, and it needs no training by gradient descent.

This paper works out that correspondence for convolutional networks, the kind that look at images with small sliding filters, and for residual networks, which add shortcut connections. The result is a kernel, a formula for how similar two images look to the model. The kernel has almost no settings to tune, and it is cheap to compute because it only needs the variances of the activations, not every pairwise covariance.

Used as a classifier on handwritten digits, the resulting model is competitive with earlier Gaussian-process methods that do not use a learned neural network. It is not as accurate as a trained residual network.

TL;DR

• With random Gaussian filters and biases, a convolutional or residual network becomes a Gaussian process as the number of channels goes to infinity §2.
• The kernel is a recursion over layers that needs only per-location variances, so it costs about one forward pass per pair of images with one filter per layer §3.3.
• The 32-layer residual kernel gets 0.84% MNIST test error, author-reported, with 50,000 training images and no data augmentation §4 tab. 1.
• Use it when you want exact Bayesian predictions with two hyperparameters, and not when you need the best accuracy on images.

Problem

Bayesian inference is exact and gives calibrated uncertainty for Gaussian processes, but a convolutional network has too many parameters for exact Bayesian inference, and the usual approximations are rough §1. Earlier ways to combine convolutions with Gaussian processes either put a large neural-network kernel on top, which reintroduces overfitting risk, or built convolutional Gaussian processes that are expensive to evaluate §1 §5. Fully connected networks have a known Gaussian-process limit, but they do not exploit translation invariance, which is why convolutional and residual networks were the open case §1.

Method

The network takes an input image with C⁰ channels and H⁰·W⁰ locations. Each hidden layer applies a bias and a sum over input channels of filter responses, followed by an elementwise nonlinearity φ eq. 1 eq. 2. Weights and biases of each layer are independent zero-mean Gaussians eq. 3. The weight variance is σ_w² divided by the number of input channels, and the bias variance is σ_b² eq. 3 §2.1. The paper notes that the weight variance may also be divided by the filter size, which gives He initialisation; the main derivation does not depend on that choice §2.1.

Limit. The argument takes the channel count to infinity one layer at a time, starting from the first layer. At each step the feature maps of a layer become independent Gaussian processes over the inputs, using the multivariate central limit theorem §2.2 eq. 4 eq. 6. The appendix extends the argument to the case where all layers grow at the same time, reusing results from Matthews et al. (2018) without reproving them app. §7.2.

Kernel. For a pair of inputs, only the covariance at corresponding locations is needed, not the covariance between all pairs of locations. Under the prior, the downstream weights are independent with zero mean, so the covariance reduces to a sum over the filter patch of the expected product of activations §3.1 eq. 7 eq. 8 eq. 9. Writing the layer-one kernel and the recursion gives the following eq. 10 eq. 11:

K^{1}_\mu(x, x') = \sigma_b^2 + \frac{\sigma_w^2}{C^0} \sum_{c=1}^{C^0} \sum_{\nu \in \text{patch } \mu} x_{c,\nu}\, x'_{c,\nu} \qquad [\text{eq. 10}]

K^{\ell+1}_\mu(x, x') = \sigma_b^2 + \sigma_w^2 \sum_{\nu \in \text{patch } \mu} V^{\ell}_\nu(x, x') \qquad [\text{eq. 11}]

where V^ℓ is the expected product of post-nonlinearity activations at location ν eq. 12. Only the three entries K(x,x), K(x,x') and K(x',x') are needed at each location §3.2.

ReLU. For the nonlinearity used in the main experiments, V has a closed form in terms of the three kernel entries and the angle θ between the two inputs eq. 13. The text writes φ as max(0, x), but the constant in eq. 13 matches √2·max(0, x); see Pitfalls.

Erf. For φ = erf, V is given by a different closed form in app. §7.4, also in terms of the three kernel entries.

Algorithm. alg. 1 iterates the layer-one initialisation, the closed-form expectation and the recursion, and returns one number for the output.

# ConvNet GP kernel  [alg. 1] [eq. 10] [eq. 11] [eq. 13]
# Input:  images x, x2 of shape (C0, H0 * W0); number of layers L
# Output: scalar kernel value k(x, x2)

def convnet_kernel(x, x2, n_layers, sigma_w2, sigma_b2):
    # Layer-1 covariances at every patch location  [eq. 10]
    # Channel average of products, then summed over the filter patch
    k_xx = sigma_b2 + sigma_w2 * patch_sum(mean_channels(x * x))
    k_yy = sigma_b2 + sigma_w2 * patch_sum(mean_channels(x2 * x2))
    k_xy = sigma_b2 + sigma_w2 * patch_sum(mean_channels(x * x2))
    for layer in range(n_layers):
        # Expected product of activations, per location  [eq. 13]
        # relu_expectation uses the sqrt(2)*ReLU form (see Pitfalls)
        v_xy = relu_expectation(k_xy, k_xx, k_yy)
        # The diagonal case satisfies v_xx = k_xx for this nonlinearity
        # Propagate all three entries through the next filter  [eq. 11]
        k_xx = sigma_b2 + sigma_w2 * patch_sum(k_xx)
        k_yy = sigma_b2 + sigma_w2 * patch_sum(k_yy)
        k_xy = sigma_b2 + sigma_w2 * patch_sum(v_xy)
    # Read-out: the paper gives one output location [alg. 1]
    # The code averages over the final locations (see Paper vs code)
    return mean(k_xy)

Efficiency. The three kernel entries per location suffice, so the cost is about a forward pass per pair of inputs with one filter per layer §3.3. The paper also states that the variances at every location can be computed at once as a convolution §3.3. Asymptotically the kernel matrix costs O(N²·L·D) and inverting it O(N³) for N examples, L layers and input dimension D §4.

Residual networks. Adding a skip connection of span s changes the recursion to add the kernel from s layers back, plus the usual bias and branch terms eq. 14 eq. 15. This matches pre-activation shortcuts §3.4.

# Residual kernel step  [eq. 14] [eq. 15]
# k_skip is the kernel s layers back; v_xy is the branch term from [eq. 13]

def residual_step(k_skip, v_xy, sigma_w2, sigma_b2):
    # The skip adds the stored kernel, not the branch variance  [eq. 15]
    return k_skip + sigma_b2 + sigma_w2 * patch_sum(v_xy)

Pooling is not used. Average pooling needs off-diagonal covariances within a filter, which removes the efficiency gain, so the ResNet kernel ends with a dense layer in its place §4.

Classification is handled as regression. The targets are one-hot vectors with entries +1 and −1, and prediction is the largest output. The posterior is exact GP regression with noiseless targets, using the N×N kernel matrix §4.

Hyperparameters

Other settings the code uses, not stated in the paper: MNIST pixel values as produced by TensorFlow's loader save_kernels.py; ResNet depth 32 gives 3 blocks of 5 residual units with strides 1, 2, 2 [code: save_kernels_resnet.py#L36-L38] save_kernels_resnet.py.

Results

Author-reported. Setup: MNIST, 50,000 training and 10,000 validation images, 10,000 test images, exact GP regression on ±1 targets, no augmentation §4. Hyperparameters chosen on validation, no retraining afterwards §4.

Other results:

• The ResNet GP beats every comparable non-parametric method in the table, including the NNGP, which the paper calls state of the art for non-convolutional networks §4. It does not reach the parametric methods §4.
• Sampling 10,000 finite 32-layer ResNets with 3, 10, 30 and 100 first-layer channels: the output distribution approaches the limiting density around 100 channels, and the moments already match closely at 30 channels §4 fig. 2.
• Calibration: with a proper classification likelihood (RobustMax in GPflow) and 1000 inducing points, the paper's kernel gives 2.4% error against 3.4% for an RBF kernel, and both have calibration curves near the diagonal app. §7.3 fig. 3.
• Cost: computing the kernel matrix and the validation and test kernels took 3 h 40 min on two Tesla P100 GPUs for the ResNet. Inverting it and computing validation and test error took 43.25 ± 8.8 s on one P100 §4.

Paper vs code

Code: official repository convnets-as-gps/convnets-as-gps, cloned at commit 778c76c6a4adf91fa252ddc9bf1ee2aebbb27ef8. It is linked from the paper's abstract page footnote, so it is taken as the official code. Its LICENSE file is BSD 2-Clause, but setup.py says "Apache License 2.0" and gives a different repository URL setup.py. This card uses the LICENSE file.

• Nonlinearity constant. The text defines φ as max(0, x) §3.2. The ConvNet path uses ExReLU(multiply_by_sqrt2=True), which is √2·max(0, x) save_kernels.py exkern.py; its closed form divides by π, which matches eq. 13. The ResNet path uses plain ReLU with the 1/(2π) form save_kernels_resnet.py exkern.py. The paper does not say which constant its ResNet kernel uses.
• Weight-variance normalisation. The ConvNet kernel uses σ_w² times the filter sum over the patch, with channels averaged at the input, which matches eq. 10 and eq. 11 [code: deep_ckern/init.py#L86-L87] [code: deep_ckern/init.py#L139]. The ResNet convolutions divide the weight variance by the fan-in. The kernel tracks one channel, so the fan-in it uses is k², which is 9 for 3×3 filters. This is the He-style option the paper mentions but does not use in its equations resnet.py. Because residual sums combine branch and shortcut, this scale is not only a global factor; see Pitfalls.
• Read-out. The paper's output is one location with no normalisation stated alg. 1. The code averages over the final spatial locations for both the ConvNet and the ResNet [code: deep_ckern/init.py#L116-L118] [code: deep_ckern/resnet.py#L88-L89].
• Skip connections (ConvNet path). The code adds a stored variance before the linear step, every skip_freq layers [code: deep_ckern/init.py#L93-L99]. I did not verify that this equals eq. 15 term by term.
• Search ranges. The paper does not give them. The shipped sampler draws σ_w² from 0.5 to 8.5 and σ_b² from 0.2 to 8.2; layers 2 to 8 (no skip) or 4 to 15 (skip allowed); filters 3 to 7 or 3 to 6; stride 1 to 3; padding VALID or SAME; ReLU or erf save_kernels.py. When a skip is chosen, the code forces SAME padding and stride 1 save_kernels.py.
• Skip assignment. Every seed below 56 gets a positive skip frequency regardless of the coin flip, so the ConvNet GP rows with no skip cannot be identified from the script's sampling alone save_kernels.py. The paper does not say how the two groups were separated.
• Random search size. The paper reports 27 kernels per row tab. 1. The run script loops over seeds 1 to 60 plus 1234 and says 27 is the number that survived failures with VALID padding run_all_experiments.bash. Seed 1234 uses a fixed architecture with 5×5 and 2×2 filters that the paper does not describe save_kernels.py.
• Training size. The paper uses 50,000 training and 10,000 validation images §4. The run script sets 400 and 400 and says to raise them to 50,000 and 10,000 for the paper's experiments run_all_experiments.bash. The Python defaults are 50,000 and 10,000 save_kernels.py.
• Selection step. The paper picks the best validation run §4. The repository writes validation and test errors per seed, but I found no script that picks the best run classify_gp.py.
• Import. classify_gp.py calls scipy.linalg.solve at line 79, but I found no scipy import in the file or in the modules it imports. As shipped, this line would raise NameError, unless the environment injects scipy some other way.
• Regression. The solve uses no jitter by default classify_gp.py, matching the noiseless setting in the paper §4.
• Dataset path. The run uses TensorFlow's MNIST loader, which reads and downloads into a relative MNIST_data/ directory save_kernels.py.
• Hidden costs in the README. About 15 GB of disk for one run, and about 24 GB of RAM for classification at 64-bit precision, because the main kernel matrix is about 12 GB README.md.

The erf closed form, the channel-averaged input variance and the bias term of the ConvNet path match the paper's formulas.

Pitfalls

• The nonlinearity is ambiguous. Text says max(0, x); the displayed constant in eq. 13 corresponds to √2·max(0, x). Plugging plain ReLU into eq. 13 halves the expected products V, both on the diagonal and across inputs, so the propagated covariances change at every layer. The authors' ConvNet code uses the √2 form.
• The σ_w² note in Table 2 is not general. The table says σ_w² for the ResNet does not matter except for numerical stability, because the scale can be moved out of a plain ReLU stack tab. 2. With residual sums in eq. 15, the branch scale changes the balance against the shortcut, so the code's ResNet choices (ReLU constant and fan-in divisions) are not a pure rescaling of the paper's equations. Keep the code's convention if you want to reproduce the 0.84%.
• Output scale. The paper does not state how the dense output's variance is normalised. The code averages over locations, which changes the kernel by a constant factor. For a noiseless GP regression argmax, a constant factor on the whole kernel leaves predictions unchanged; a per-layer factor does not.
• Cost is front-loaded. Computing the kernel matrix dominated the authors' time (3 h 40 min) far more than solving it (43 s) §4. The README gives about 15 GB of disk per run and about 12 GB for the main MNIST kernel matrix README.md.
• Regression stands in for classification. The model is trained as squared-error regression on ±1 targets, not with a classification likelihood §4. Calibration uses a separate experiment with a proper likelihood app. §7.3. Compare error rates with that in mind.
• Limits. The main-text proof takes layer limits one at a time. The appendix covers simultaneous limits, but relies on Matthews et al. (2018) for several lemmas and states this is an extension app. §7.2. Read the appendix's assumptions before relying on the limit for a specific architecture.
• Random sampling. The 27 kernels are not a fixed architecture. Expect failures for some architectures, and the reported count depends on which ones failed run_all_experiments.bash.
• Appendix numbering. Section numbers in the appendix (§7.x) follow the LaTeX source, where \section{Appendix} is numbered by the article class. The compiled PDF may letter them differently. Check the heading names.
• Equation numbers. Equation numbers are counted from the LaTeX source in order, including the unlabelled align environments in §2.2. The compiled PDF should match, but check when citing.

Lineage

• Builds on: Neal (1996) for single-hidden-layer networks as Gaussian processes; Lee et al. (2017) for deep fully connected networks with infinite width, whose layer-wise argument this paper follows §2.2; Matthews et al. (2018), whose simultaneous-limit argument the appendix extends app. §7.2. Neither has a card, so they are named and not linked.
• Compared against: NNGP (Lee et al. 2017); convolutional GPs (van der Wilk et al. 2017); deep convolutional GPs (Kumar et al. 2018); deep kernel learning (Wilson et al. 2016; Bradshaw et al. 2017); a standard ResNet (He et al. 2016; Chen et al. 2018) tab. 1. None has a card, so none is linked.
• Superseded by: not stated in the paper.