# Incorporating Prior Knowledge into Neural Networks through an Implicit
Composite Kernel

Ziyang Jiang  Tongshu Zheng  Yiling Liu  David Carlson

###### Abstract

It is challenging to guide neural network (NN) learning with prior knowledge.
In contrast, many known properties, such as spatial smoothness or seasonality,
are straightforward to model by choosing an appropriate kernel in a Gaussian
process (GP). Many deep learning applications could be enhanced by modeling
such known properties. For example, convolutional neural networks (CNNs) are
frequently used in remote sensing, which is subject to strong seasonal
effects. We propose to blend the strengths of deep learning and the clear
modeling capabilities of GPs by using a composite kernel that combines a
kernel implicitly defined by a neural network with a second kernel function
chosen to model known properties (e.g., seasonality). We implement this idea
by combining a deep network and an efficient mapping based on the Nyström
approximation, which we call Implicit Composite Kernel (ICK). We then adopt a
sample-then-optimize approach to approximate the full GP posterior
distribution. We demonstrate that ICK has superior performance and flexibility
on both synthetic and real-world data sets. We believe that ICK framework can
be used to include prior information into neural networks in many
applications.

Machine Learning, ICML

##  1 Introduction

In complex regression tasks, input data often contains _multiple sources of
information_. These sources can be presented in both high-dimensional (e.g.
images, audios, texts, etc.) and low-dimensional (e.g. timestamps, spatial
locations, etc.) forms. A common approach to learn from high-dimensional
information is to use neural networks (NNs) (Goodfellow et al., 2016; LeCun et
al., 2015), as NNs are powerful enough to capture the relationship between
complex high-dimensional data and target variables of interest. In many areas,
NNs are standard practice, such as the dominance of Convolutional Neural
Networks (CNNs) for image analysis (Jiang et al., 2022; Zheng et al., 2021,
2020). In contrast, for low-dimensional information, we usually have some
prior knowledge on how the information relates to the predictions. As a
concrete example, consider a remote sensing problem where we predict ground
measurements from satellite imagery with associated timestamps. A priori, we
expect the ground measurements to vary periodically with respect to time
between summer and winter due to seasonal effects. We would typically use a
CNN to capture the complex relationship between the imagery and the ground
measurements. In this case, we want to guide the learning of the CNN with our
prior knowledge about the seasonality. This is challenging because knowledge
represented in NNs pertains mainly to correlation between network units
instead of quantifiable statements (Marcus, 2018).

Conversely, Gaussian processes (GPs) have been used historically to
incorporate relevant prior beliefs by specifying the appropriate form of its
kernel (or covariance) function (Bishop & Nasrabadi, 2006; Williams &
Rasmussen, 2006). One approach to modeling multiple sources of information is
to assign a relevant kernel function to each source of information
respectively and combine them through addition or multiplication, resulting in
a composite kernel function (Duvenaud, 2014). This formulation means that
specifying a kernel to match prior beliefs on one source of information is
straightforward. Such composite kernel learning techniques are extensively
used in many application areas such as multi-media data (McFee et al., 2011),
neuroimaging (Zhang et al., 2011), spatial data analysis, and environmental
data analysis (Kim et al., 2005; Petelin et al., 2013). In view of the clear
modeling capabilities of GP, it is desirable to examine how a NN could be
imbued with the same modeling ease.

In recent years, researchers have come up with a variety of methods to
incorporate prior knowledge into NNs. These efforts can be broken into many
categories, such as those that add prior information through loss terms like
physics-informed NNs (Lagaris et al., 1998; Moseley et al., 2020). Here, we
focus on the major category of those methods that build integrated models of
NNs and GPs with various structures (Van der Wilk et al., 2017; Wilson et al.,
2016, 2011). Related to our proposed methodology, Pearce et al. (Pearce et
al., 2020) exploited the fact that a Bayesian neural network (BNN)
approximates a GP to construct additive and multiplicative kernels, but they
were limited to specific predefined kernels. Matsubara et al. (Matsubara et
al., 2020) then resolved this limitation by constructing priors of BNN
parameters based on the ridgelet transform and its dual, but they did not
explicitly show how their approach works for data with multiple sources of
information. To our knowledge, no existing approach allows a modeler to choose
any appropriate kernel over multiple sources.

We address this limitation by presenting a simple yet novel Implicit Composite
Kernel (ICK) framework, which processes high-dimensional information using a
kernel implicitly defined by a neural network and low-dimensional information
using a chosen kernel function. The low-dimensional kernels are mapped into
the neural network framework to create a straightforward and simple-to-learn
implementation. Our key results and contributions are:

  * •

We analytically show our ICK framework, under reasonable assumptions, is
approximately equivalent to sampling from a Gaussian process regression (GPR)
model with a composite kernel _a priori_.

  * •

We adopt a sample-then-optimize procedure to ICK to approximate the full
posterior distribution of a GP with a composite kernel.

  * •

We show that ICK yields better performance on prediction and forecasting
tasks, even with limited data.

  * •

We show that ICK can flexibly capture the patterns of the low-dimensional
information without bespoke pre-processing procedures or complex NN
structures.

Based on these contributions, we believe ICK is useful in learning from
complex _hybrid_ data with prior knowledge, especially in remote sensing and
spatial statistics.

##  2 Related Work

##### Equivalence between NNs and GPs

The equivalence between GPs and randomly initialized single-layer NNs with
infinite width was first shown by Neal (Neal, 1996). With the development of
modern deep learning, researchers further extended this relationship to deep
networks (Lee et al., 2017; Matthews et al., 2018) and convolutional neural
networks (CNNs) (Garriga-Alonso et al., 2018). This relationship is crucial
for showing the resemblance between GPR and our ICK framework, which is
discussed in Section 4.1.

##### NNs with prior knowledge

As mentioned before, one approach to equip NNs with prior knowledge is to
modify the loss function. For example, Lagaris et al. (Lagaris et al., 1998)
solved differential equations (DEs) using NNs by setting the loss to be a
function whose derivative satisfies the DE conditions. Another approach is to
build integrated models of NNs and kernel-based models. For example, Wilson et
al. (Wilson et al., 2011) implemented a regression network with GP priors over
latent variables and made inference by approximating the posterior using
Variational Bayes or sampling from the posterior using Gibbs sampling scheme.
Garnelo et al. (Garnelo et al., 2018) introduced a class of neural latent
variable models called Neural Processes (NPs) which are capable of learning
efficiently from the data and adapting rapidly to new observations. Our ICK
framework fuses prior knowledge into NNs by modulating the learnt features
using another set of features outputted from a kernel-based mapping, which can
also be viewed as an integrated model of NNs and kernel machines.

##### GP with composite kernels

Composite kernel GPs are widely used in both machine learning (Duvenaud, 2014;
Williams & Rasmussen, 2006) and geostatistical modeling (Datta et al., 2016;
Gelfand & Schliep, 2016). GPR in geostatistical modeling is also known as
_kriging_ (Journel & Huijbregts, 1976; Krige, 1951), which serves as a
surrogate model to replace expensive function evaluations. The inputs for a
composite GP are usually low-dimensional (e.g. spatial distance) as GPs do not
scale well with the number of samples for high-dimensional inputs (Bouhlel &
Martins, 2019; Bouhlel et al., 2016). To overcome this issue, Pearce et al.
(Pearce et al., 2020) and Matsubara et al. (Matsubara et al., 2020) developed
BNN analogues for composite GPs. Similar to these studies, our ICK framework
can also be viewed as a simulation for composite GPs.

##### Approximation methods for GP

For large data sets, approximation methods are needed as exact kernel learning
and inference scales $\mathcal{O}(N^{3})$. Nyström low-rank matrix
approximation (Drineas et al., 2005; Williams & Seeger, 2000) and Random
Fourier Features (Rahimi & Recht, 2007, 2008) are two commonly used methods. A
common technique is to choose inducing points as pseudo-inputs to efficiently
approximate the full kernel matrix (Snelson & Ghahramani, 2005; Hensman et
al., 2013). Our work is inspired by these approximation techniques and we use
them as _transformation functions_ to map the kernel matrix into latent space
representations in Section 4.2.

##  3 Background

Before elaborating on the ICK framework, we introduce our notation, briefly go
over the concepts of composite GPs, and describe the relationship between GPs
and NNs.

###  3.1 Problem Setup

To formalize the problem, we have a training data set which contains $N$ data
points
$\boldsymbol{X}=[\boldsymbol{x}_{i}]_{i=1}^{N}=[\boldsymbol{x}_{1},\boldsymbol{x}_{2},...,\boldsymbol{x}_{N}]^{T}$
and the corresponding labels of these data points are
$\boldsymbol{y}=[y_{i}]_{i=1}^{N}=[y_{1},y_{2},...,y_{N}]^{T}$ where
$y_{i}\in\mathbb{R}$. Each data point
$\boldsymbol{x}_{i}=\left\\{\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{i}^{(2)},...,\boldsymbol{x}_{i}^{(M)}\right\\}$
is composed of information from $M$ different sources where the $m^{th}$
source of information of the $i^{th}$ data point is denoted as
$\boldsymbol{x}_{i}^{(m)}\in\mathbb{R}^{D_{m}}$. Our goal is to learn a
function
$\hat{y}_{i}=f(\boldsymbol{x}_{i}):\left\\{\mathbb{R}^{D_{1}},\mathbb{R}^{D_{2}},...,\mathbb{R}^{D_{M}}\right\\}\rightarrow\mathbb{R}$
which takes in a data point $\boldsymbol{x}_{i}$ and outputs a predicted value
$\hat{y}_{i}$.

###  3.2 Composite GPs

A Gaussian process (GP) describes a distribution over functions (Williams &
Rasmussen, 2006). A key property of GP is that it is completely defined by a
mean function $\mu(\boldsymbol{x})$ and a kernel function
$K(\boldsymbol{x},\boldsymbol{x}^{\prime})$ where $\boldsymbol{x}$ and
$\boldsymbol{x}^{\prime}$ represent different samples from the training
dataset. The mean function $\mu(\boldsymbol{x})$ is often assumed to be zero
for simplicity. In that case, the outcome function is

$f(\boldsymbol{x})\sim\mathcal{GP}\left(0,K(\boldsymbol{x},\boldsymbol{x}^{\prime})\right).$
(1)

Any finite subset of random variables has a multivariate Gaussian distribution
with mean $\boldsymbol{0}$ and kernel matrix $\boldsymbol{K}$ whose entries
can be calculated as
$\boldsymbol{K}_{ij}=K(\boldsymbol{x}_{i},\boldsymbol{x}_{j})$ where $1\leq
i,j\leq N$. In many situations, the full kernel function is built by a
composite kernel by combining simple kernels through addition
$K^{\text{comp}}(\boldsymbol{x},\boldsymbol{x}^{\prime})=K_{1}(\boldsymbol{x},\boldsymbol{x}^{\prime})+K_{2}(\boldsymbol{x},\boldsymbol{x}^{\prime})$
or multiplication
$K^{\text{comp}}(\boldsymbol{x},\boldsymbol{x}^{\prime})=K_{1}(\boldsymbol{x},\boldsymbol{x}^{\prime})K_{2}(\boldsymbol{x},\boldsymbol{x}^{\prime})$
(Duvenaud, 2014). A useful property that ICK exploits is that $K_{1}$ and
$K_{2}$ can take different subparts of $\boldsymbol{x}$ as their inputs. For
example,
$K^{\text{comp}}(\boldsymbol{x},\boldsymbol{x}^{\prime})=K_{1}\left(\boldsymbol{x}^{(1)},{\boldsymbol{x}^{(1)}}^{\prime}\right)+K_{2}\left(\boldsymbol{x}^{(2)},{\boldsymbol{x}^{(2)}}^{\prime}\right)$
or
$K^{\text{comp}}(\boldsymbol{x},\boldsymbol{x}^{\prime})=K_{1}\left(\boldsymbol{x}^{(1)},{\boldsymbol{x}^{(1)}}^{\prime}\right)K_{2}\left(\boldsymbol{x}^{(2)},{\boldsymbol{x}^{(2)}}^{\prime}\right)$.

###  3.3 Correspondence between GPs and NNs

Neal (Neal, 1996) proved that a single-hidden layer network with infinite
width is _exactly equivalent_ to a GP over data indices $i=1,2,...,N$ under
the assumption that the weight and bias parameters of the hidden layer are
i.i.d. Gaussian with zero mean. This statement was then extended to deep NNs
(Lee et al., 2017; Matthews et al., 2018) and convolutional NNs (Garriga-
Alonso et al., 2018). Specifically, let
$\boldsymbol{z}=f_{\text{NN}}\left(\boldsymbol{x}^{(1)}\right):\mathbb{R}^{D_{1}}\rightarrow\mathbb{R}^{p}$
be the latent representation extracted from $\boldsymbol{x}^{(1)}$ where $p$
is the dimension of the extracted representation and $f_{\text{NN}}$ is a
neural network with zero-mean i.i.d. parameters and continuous activation
function $\phi$ which satisfies the linear envelope property

$|\phi(u)|\leq c+m|u|\quad\forall u\in\mathbb{R},$ (2)

if there exists $c,m\geq 0$. This includes many standard nonlinearities. The
$k^{th}$ entry of this representation will converge in distribution to a NN-
implied GP in the infinite width limit

$f_{\text{NN}}\left(\boldsymbol{x}^{(1)}\right)_{k}\xrightarrow[]{d}\mathcal{GP}\left(0,K_{\text{NN}}\left(\boldsymbol{x}^{(1)},{\boldsymbol{x}^{(1)}}^{\prime}\right)\right).$
(3)

Here $K_{\text{NN}}$ in Equation 3 denotes the covariance function of the
equivalent GP as the width of $f_{\text{NN}}$ goes to infinity and can be
computed numerically in a recursive manner (Lee et al., 2017). That is to say,
the $k^{th}$ component $z_{k}$ of the representation extracted by the network
has zero mean
$\mathbb{E}_{p\left(\boldsymbol{\theta}^{(1)}\right)}\left[z_{ik}^{(1)}\right]=0$
for all $i=1,2,...,N$ where $\theta$ represents the network parameters. The
covariance between $z_{ik}^{(1)}$ and $z_{jk}^{(1)}$ for _different_ data
indices $i,j=1,2,...,N$ can be approximated as
$\text{cov}\left(z_{ik}^{(1)},z_{jk}^{(1)}\right)=\mathbb{E}_{p\left(\boldsymbol{\theta}^{(1)}\right)}\left[z_{ik}^{(1)}z_{jk}^{(1)}\right]\approx
K_{\text{NN}}\left(\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{j}^{(1)}\right)$
where $\boldsymbol{x}_{i}^{(1)}$ and $\boldsymbol{x}_{j}^{(1)}$ are the
corresponding inputs in case the network width is finite.

Figure 1: Given data containing 2 sources of information
$\boldsymbol{x}^{(1)}$ and $\boldsymbol{x}^{(2)}$, we can process the data
using either (Left) a composite Gaussian process regression (GPR) model or
(Right) our ICK framework where $\boldsymbol{x}^{(1)}$ is processed with a
neural network $f_{\text{NN}}(\cdot)$ and $\boldsymbol{x}^{(2)}$ is processed
with $g(\cdot)$ where $g(\cdot)$ consists of a kernel function $K_{2}$ and
some transformation which maps the kernel matrix $\boldsymbol{K}_{2}$ into the
latent space.

##  4 Implicit Composite Kernel

We show the structure of a composite GPR model and our ICK framework in Figure
1. To make the illustration clear, we limit ourselves to data with information
from 2 different sources
$\boldsymbol{x}=\left\\{\boldsymbol{x}^{(1)},\boldsymbol{x}^{(2)}\right\\}$
where $\boldsymbol{x}^{(1)}$ is high-dimensional and $\boldsymbol{x}^{(2)}$ is
low-dimensional (i.e. $D_{1}\gg D_{2}$) with some known relationship with the
target $y$. We are inspired by composite GPR, which computes 2 different
kernel matrices $\boldsymbol{K}_{1}$ and $\boldsymbol{K}_{2}$ and then
combines them into a single composite kernel matrix
$\boldsymbol{K}^{\text{comp}}$. However, as discussed before, it is more
suitable to use a NN to learn from the high dimensional information
$\boldsymbol{x}^{(1)}$. In our ICK framework, we process
$\boldsymbol{x}^{(1)}$ with a NN
$f_{\text{NN}}(\cdot):\mathbb{R}^{D_{1}}\rightarrow\mathbb{R}^{p}$ and
$\boldsymbol{x}^{(2)}$ with a mapping
$g(\cdot):\mathbb{R}^{D_{2}}\rightarrow\mathbb{R}^{p}$, which consists of a
kernel function $K_{2}$ followed by a kernel-to-latent-space transformation
(described in Section 4.2), resulting in two latent representations
$\boldsymbol{z}^{(1)},\boldsymbol{z}^{(2)}\in\mathbb{R}^{p}$. Then, we make a
prediction $\hat{y}$ by doing an _inner product_ between these two
representations $\hat{y}=f_{\text{NN}}\left(\boldsymbol{x}^{(1)}\right)\cdot
g\left(\boldsymbol{x}^{(2)}\right)$. Finally, the parameters of both the NN
and the kernel function are learned via gradient-based optimization methods
(Bottou et al., 2018).

In the sections below, we first show the relationship between ICK and a
composite GPR model with a multiplicative kernel _a priori_ in Section 4.1,
which is used to motivate the model form. We then show how we implement the
kernel-to-latent-space transformation in Section 4.2. We then use the idea
from Lee et al. (Lee et al., 2017) that finite-width NNs trained with
stochastic gradient descent (SGD) mimic Bayesian inference on a GP. However,
without Bayesian formulation on the NN and kernel function parameters, the
prediction $\hat{y}$ can only be viewed as a _posterior sample_ corresponding
to an _approximate_ GP prior. To enable a full GP posterior _approximation_
for ICK, we adopt a sample-then-optimize procedure (Matthews et al., 2017) to
formulate an ensemble algorithm, as elaborated in Section 4.3.

###  4.1 Resemblance between Composite GPR and ICK

We will analytically prove the following theorem for data with information
from 2 different sources
$\boldsymbol{x}=\left\\{\boldsymbol{x}^{(1)},\boldsymbol{x}^{(2)}\right\\}$
for clarity, and we note this theorem can be straightforwardly extended to
$M>2$.

######  Theorem 4.1.

Let $f_{\text{NN}}:\mathbb{R}^{D_{1}}\rightarrow\mathbb{R}^{p}$ be a NN
function with random weights and
$g:\mathbb{R}^{D_{2}}\rightarrow\mathbb{R}^{p}$ be a mapping function, and
define an inner product between the representations
$\hat{y}=f_{\text{ICK}}\left(\boldsymbol{x}^{(1)},\boldsymbol{x}^{(2)}\right)=f_{\text{NN}}\left(\boldsymbol{x}^{(1)}\right)\cdot
g\left(\boldsymbol{x}^{(2)}\right)={\boldsymbol{z}^{(1)}}^{T}\boldsymbol{z}^{(2)}$.
Then $f_{\text{ICK}}$ will converge in distribution to a GP

$f_{\text{ICK}}\xrightarrow[]{d}\mathcal{GP}\left(0,K_{\text{NN}}^{(1)}\left(\boldsymbol{x}^{(1)},{\boldsymbol{x}^{(1)}}^{\prime}\right)K^{(2)}\left(\boldsymbol{x}^{(2)},{\boldsymbol{x}^{(2)}}^{\prime}\right)\right),$
(4)

if $f_{\text{NN}}$ is a neural network with zero-mean i.i.d. parameters and
continuous activation function $\phi$ which satisfies the linear envelope
property in Equation 2 and $g$ includes the following deterministic kernel-to-
latent-space transformation for all $1\leq i,j\leq N$

$K^{(2)}\left(\boldsymbol{x}_{i}^{(2)},\boldsymbol{x}_{j}^{(2)}\right)\approx{\boldsymbol{z}_{i}^{(2)}}^{T}\boldsymbol{z}_{j}^{(2)}=g\left(\boldsymbol{x}_{i}^{(2)}\right)^{T}g\left(\boldsymbol{x}_{j}^{(2)}\right),$
(5)

where $K_{\text{NN}}^{(1)}$ is a NN-implied kernel and $K^{(2)}$ is any valid
kernel of our choice.

To prove Theorem 1, we first make the assumption below.

######  Assumption 4.2.

For latent representations $\boldsymbol{z}_{i}^{(1)}$ and
$\boldsymbol{z}_{j}^{(1)}$ extracted from different data points
$\boldsymbol{x}_{i}$ and $\boldsymbol{x}_{j}$ where $i\neq j$, the
interactions between different entries of $\boldsymbol{z}_{i}^{(1)}$ and
$\boldsymbol{z}_{j}^{(1)}$ can be reasonably ignored. In other words, let
$\theta^{(1)}$ be the parameters of the neural network which takes in
$x^{(1)}$ and outputs $\boldsymbol{z}^{(1)}$, we have
$\mathbb{E}_{p\left(\theta^{(1)}\right)}\left[z_{ik}^{(1)}z_{jl}^{(1)}\right]=0$
for all $k\neq l$.

A detailed discussion of Assumption 4.2 and its reasonableness is provided in
Appendix A.2. With Assumption 4.2, let
$\boldsymbol{\Theta}=\left\\{\boldsymbol{\theta}^{(1)},\boldsymbol{\theta}^{(2)}\right\\}$
represent the parameters of ICK, we can calculate the covariance between
$\hat{y}_{i}$ and $\hat{y}_{j}$ for different data indices $i\neq j$:

$\displaystyle\text{cov}(\hat{y}_{i},\hat{y}_{j})$
$\displaystyle=\mathbb{E}_{p(\boldsymbol{\Theta})}[\hat{y}_{i}\hat{y}_{j}]-\mathbb{E}_{p(\boldsymbol{\Theta})}[\hat{y}_{i}]\mathbb{E}_{p(\boldsymbol{\Theta})}[\hat{y}_{j}]$
(6)
$\displaystyle=\textstyle\mathbb{E}_{p(\boldsymbol{\Theta})}\left[\left(\sum_{k=1}^{p}z_{ik}^{(1)}z_{ik}^{(2)}\right)\left(\sum_{k=1}^{p}z_{jk}^{(1)}z_{jk}^{(2)}\right)\right]$
(7)
$\displaystyle=\textstyle\mathbb{E}_{p(\boldsymbol{\Theta})}\left[\sum_{k=1}^{p}\sum_{l=1}^{p}z_{ik}^{(1)}z_{jl}^{(1)}z_{ik}^{(2)}z_{jl}^{(2)}\right]$
(8)
$\displaystyle=\textstyle\mathbb{E}_{p(\boldsymbol{\Theta})}\left[\sum_{k=1}^{p}z_{ik}^{(1)}z_{jk}^{(1)}z_{ik}^{(2)}z_{jk}^{(2)}\right]$
(9)
$\displaystyle=\textstyle\sum_{k=1}^{p}\mathbb{E}_{p\left(\boldsymbol{\theta}^{(1)}\right)}\left[z_{ik}^{(1)}z_{jk}^{(1)}\right]\mathbb{E}_{p\left(\boldsymbol{\theta}^{(2)}\right)}\left[z_{ik}^{(2)}z_{jk}^{(2)}\right]$
(10) $\displaystyle\approx\textstyle
K_{\text{NN}}^{(1)}\left(\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{j}^{(1)}\right)\sum_{k=1}^{p}\mathbb{E}_{p\left(\boldsymbol{\theta}^{(2)}\right)}\left[z_{ik}^{(2)}z_{jk}^{(2)}\right].$
(11)

Here, from Equation 6 to Equation 7, we use the statement
$\mathbb{E}_{p\left(\boldsymbol{\theta}^{(1)}\right)}\left[z_{ik}^{(1)}\right]=0$
in Section 3.3 and $\boldsymbol{\theta}^{(1)}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\boldsymbol{\theta}^{(2)}$, which leads to
$\mathbb{E}_{p(\boldsymbol{\Theta})}[\hat{y}_{i}]=\mathbb{E}_{p(\boldsymbol{\Theta})}[\hat{y}_{j}]=0$.
From Equation 8 to Equation 9, we get rid of all the cross terms under
Assumption 4.2 and $\boldsymbol{\theta}^{(1)}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\boldsymbol{\theta}^{(2)}$. Specifically, we
have
$\mathbb{E}_{p(\boldsymbol{\Theta})}\left[z_{ik}^{(1)}z_{jl}^{(1)}z_{ik}^{(2)}z_{jl}^{(2)}\right]=\mathbb{E}_{p\left(\boldsymbol{\theta}^{(1)}\right)}\left[z_{ik}^{(1)}z_{jl}^{(1)}\right]\mathbb{E}_{p\left(\boldsymbol{\theta}^{(2)}\right)}\left[z_{ik}^{(2)}z_{jl}^{(2)}\right]=0$
for all $k\neq l$. From Equation 9 to Equation 10, we again make use of
$\boldsymbol{\theta}^{(1)}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\boldsymbol{\theta}^{(2)}$. From Equation 10
to Equation 11, we use the statement
$\mathbb{E}_{p\left(\boldsymbol{\theta}^{(1)}\right)}\left[z_{ik}^{(1)}z_{jk}^{(1)}\right]\approx
K_{\text{NN}}\left(\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{j}^{(1)}\right)$
from Section 3.3. If the kernel-to-latent-space transformation in $g(\cdot)$
is _deterministic_ , we can remove the expectation sign from the summation
term in Equation 11 and the covariance can be further expressed as

$\displaystyle\begin{split}\text{cov}(\hat{y}_{i},\hat{y}_{j})&\approx
K_{\text{NN}}^{(1)}\left(\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{j}^{(1)}\right)\left({\boldsymbol{z}_{i}^{(2)}}^{T}\boldsymbol{z}_{j}^{(2)}\right)\\\
&=K_{\text{NN}}^{(1)}\left(\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{j}^{(1)}\right)K^{(2)}\left(\boldsymbol{x}_{i}^{(2)},\boldsymbol{x}_{j}^{(2)}\right),\end{split}$
(12)

which means that $\hat{y}$ approximately follows a GP with a multiplicative
composite kernel
$K^{\text{comp}}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})=K_{\text{NN}}^{(1)}\left(\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{j}^{(1)}\right)K^{(2)}\left(\boldsymbol{x}_{i}^{(2)},\boldsymbol{x}_{j}^{(2)}\right)$
a priori. This serves as a consise proof of Theorem 1. We also give a more
detailed proof Appendix A.

###  4.2 Kernel-to-latent-space Transformation

We now show how we can construct an appropriate mapping $g(\cdot)$ that
approximately satisfies the assumed form of (5) and is used in the derivation
of ICK from (11) to (12). Here we adopt two methods, Nyström approximation and
Random Fourier Features (RFF), to map the kernel matrix into the latent space.
Below, we give the formulations and results for the Nyström method, and give
the methods and results for RFF in Appendix B. According to Yang et al. (Yang
et al., 2012), the Nyström method will yield much better performance than RFF
if there exists a large gap in the eigen-spectrum of the kernel matrix. In our
applications, we also observe a large eigen-gap (see details in Appendix C)
and Nyström method does generalize much better than RFF. We name our framework
with N _y_ ström method and _r_ andom Fourier Features ICK _y_ and ICK _r_ ,
respectively.

Algorithm 1 Implicit Composite Kernel-Nyström (ICK _y_)

Input: data
$\boldsymbol{X}=\left\\{\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{i}^{(2)}\right\\}_{i=1}^{N}$,
targets $\boldsymbol{y}=[y_{i}]_{i=1}^{N}$,
$f_{\text{NN}},\boldsymbol{\theta}^{(1)},K^{(2)},\boldsymbol{\theta}^{(2)}$,
learning rate $\beta$

Sample a total of $N_{B}$ minibatches
$\\{\boldsymbol{X}_{B},\boldsymbol{y}_{B}\\}_{B=1}^{N_{B}}$

for $B$ from 1 to $N_{B}$ do

for $\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{i}^{(2)}$ in
$\boldsymbol{X}_{B},i=1,...,n_{B}$ do

$\boldsymbol{z}_{i}^{(1)}=f_{\text{NN}}\left(\boldsymbol{x}_{i}^{(1)}\right)$

Define inducing points
$\hat{\boldsymbol{x}}_{1}^{(2)},...,\hat{\boldsymbol{x}}_{p}^{(2)}$

Compute $\boldsymbol{K}_{p}$:
$\left(\boldsymbol{K}_{p}\right)_{jk}=K^{(2)}\left(\hat{\boldsymbol{x}}_{j}^{(2)},\hat{\boldsymbol{x}}_{k}^{(2)}\right)$

Do Cholesky decomposition
$\left(\boldsymbol{K}_{p}\right)^{-1}=\boldsymbol{U}^{T}\boldsymbol{U}$

Compute $\boldsymbol{K}_{np}$:
$\left(\boldsymbol{K}_{np}\right)_{jk}=K^{(2)}\left(\boldsymbol{x}_{j}^{(2)},\hat{\boldsymbol{x}}_{k}^{(2)}\right)$

$\boldsymbol{z}_{i}^{(2)}=\boldsymbol{U}\left(\boldsymbol{K}_{np}^{T}\right)_{:,i}$

$\hat{y}_{i}={\boldsymbol{z}_{i}^{(1)}}^{T}\boldsymbol{z}_{i}^{(2)}$

end for

$\hat{\boldsymbol{y}}_{B}=\texttt{concat}\left(\hat{y}_{1},...,\hat{y}_{n_{B}}\right)$

Compute loss
$\mathcal{L}=\mathcal{L}\left(\boldsymbol{y}_{B},\hat{\boldsymbol{y}}_{B}\right)$

$\boldsymbol{\theta}^{(1)}\leftarrow\boldsymbol{\theta}^{(1)}-\beta\nabla_{\boldsymbol{\theta}^{(1)}}\mathcal{L}$

$\boldsymbol{\theta}^{(2)}\leftarrow\boldsymbol{\theta}^{(2)}-\beta\nabla_{\boldsymbol{\theta}^{(2)}}\mathcal{L}$

end for

Return: Predictions $\hat{\boldsymbol{y}}$ and updated parameters
$\boldsymbol{\theta}^{(1)},\boldsymbol{\theta}^{(2)}$

####  4.2.1 Nyström Method

The main idea of Nyström method (Williams & Seeger, 2000) is to approximate
the kernel matrix $\boldsymbol{K}\in\mathbb{R}^{N\times N}$ with a much
smaller low-rank matrix $\boldsymbol{K}_{q}\in\mathbb{R}^{q\times q}$ where
$q\ll N$ so both the computational and space complexity of kernel learning can
be significantly reduced, yielding

$\boldsymbol{K}\approx\hat{\boldsymbol{K}}=\boldsymbol{K}_{nq}\boldsymbol{K}_{q}^{-1}\boldsymbol{K}_{nq}^{T}.$
(13)

The entries of $\boldsymbol{K}_{q}$ and $\boldsymbol{K}_{nq}$ can be
calculated as
$\left(\boldsymbol{K}_{q}\right)_{ij}=K(\hat{x}_{i},\hat{x}_{j}),i,j\in\\{1,2,...,q\\}$
and
$\left(\boldsymbol{K}_{nq}\right)_{ij}=K(x_{i},\hat{x}_{j}),i\in\\{1,2,...,N\\},j\in\\{1,2,...,q\\}$,
respectively. $x$ represents the original data points and $\hat{x}$ represents
pre-defined inducing points (or pseudo-inputs (Snelson & Ghahramani, 2005)).
In our study, these inducing points are chosen by defining an evenly spaced
vector over the range of original data points. By performing Cholesky
decomposition $\boldsymbol{K}_{q}^{-1}=\boldsymbol{U}^{T}\boldsymbol{U}$,
where $\boldsymbol{U}\in\mathbb{R}^{q\times q}$, $\hat{\boldsymbol{K}}$ is
expressed

$\displaystyle\begin{split}\textstyle\hat{\boldsymbol{K}}&=\boldsymbol{K}_{nq}\boldsymbol{K}_{q}^{-1}\boldsymbol{K}_{nq}^{T}=\boldsymbol{K}_{nq}\boldsymbol{U}^{T}\boldsymbol{U}\boldsymbol{K}_{nq}^{T}\\\
&=\left(\boldsymbol{U}\boldsymbol{K}_{nq}^{T}\right)^{T}\left(\boldsymbol{U}\boldsymbol{K}_{nq}^{T}\right).\end{split}$
(14)

Therefore, if we set $q=p$, then we can use
$\boldsymbol{z}_{i}\triangleq\boldsymbol{U}\left(\boldsymbol{K}_{np}^{T}\right)_{:,i}$
as a kernel-to-latent-space transformation because each element in
$\boldsymbol{K}$ approximately satisfies (5) as stated in Theorem 1:
$K\left(x_{i},x_{j}\right)=K_{ij}\approx\hat{K}_{ij}=\boldsymbol{z}_{i}^{T}\boldsymbol{z}_{j}$.
Conveniently, modern deep learning frameworks can propagate gradients through
the Cholesky operation, making it straightforward to update the kernel
parameters with gradient methods. Note that as we increase the number of
inducing points $p$, the approximation error between $\boldsymbol{K}$ and
$\hat{\boldsymbol{K}}$ decreases. However, it is not recommended to set $p$
very large as updating the Cholesky decomposition requires
$\mathcal{O}(p^{3})$ . The empirical impact of $p$ on computational time and
performance is shown in Appendix D. In our experiments, only mild values of
$p$ are necessary and the impact on computational is relatively small. The
full training procedure of ICK _y_ is shown in Algorithm 1.

###  4.3 Ensemble Formulation of ICK

We adopt a sample-then-optimize approach (Matthews et al., 2017) to construct
a deep ensemble posterior approximation for ICK. There are several ways to
enable GP posterior interpretation for a deep ensemble trained with SGD.
Specifically, let the final layer of each baselearner NN be $C$-dimensional
and denote the deep ensemble as $F=\\{f_{n_{e}}\\}_{n_{e}=1}^{N_{e}}$, where
$f_{n_{e}}$ is a baselearner NN and $N_{e}$ is the total number of
baselearners in the ensemble. After training, each $f_{n_{e}}$ can be viewed
as i.i.d. samples from a multi-output GP posterior with kernel function
$K_{F}$ in the infinite width limit and the ensemble $F$ thus represents
$N_{e}$ independent draws from the GP posterior, where $K_{F}$ can be either
the _Neural Network Gaussian Process_ (NNGP) kernel (Lee et al., 2017) if we
make only the last layer trainable or the _Neural Tangent Kernel_ (NTK) (He et
al., 2020; Jacot et al., 2018) if we add a randomized and untrainable function
to each baselearner.

This deep ensemble mechanism can be easily applied to ICK _y_. Specifically,
the final prediction of ICK _y_
$\hat{y}=f_{\text{ICK}y}\left(\boldsymbol{x}^{(1)},\boldsymbol{x}^{(2)}\right)$
can be viewed as a weighted sum of
$\hat{y}=\sum_{k=1}^{p}\alpha_{k}z_{k}^{(1)}$ where
$\alpha_{k}=z_{k}^{(2)}=g\left(\boldsymbol{x}^{(2)}\right)_{k}$. Hence, if we
construct an ensemble of ICK _y_ ,
$F_{\text{ICK}y}=\left\\{f_{1},...,f_{N_{e}}\right\\}$ (as shown in Algorithm
2), and $f_{\text{NN}}$ is appropriately initialized for all baselearners
$f_{s}\in F_{\text{ICK}y},s=1,...,N_{e}$, then all trained baselearners in
$F_{\text{ICK}y}$ can be _approximately_ viewed as i.i.d. posterior samples
from a single-output GP with a multiplicative kernel $K_{F}K^{(2)}$. $K_{F}$
will either be an NNGP or an NTK and $K^{(2)}$ comes from our defined kernel.
If the parameters of $f_{\text{NN}}$ (i.e. $\boldsymbol{\theta}^{(1)}$) are
independently drawn from a zero-mean Gaussian and are fixed except the last
layer (as specified by Lee et al. (Lee et al., 2019)) for all baselearners in
the ensemble, this would correspond to an NNGP. A more detailed proof of this
argument can be found in Appendix E.

To verify our argument, we train an ICK _y_ ensemble $F_{\text{ICK}y}$
containing 300 baselearners with SGD on the same 1D regression task as
provided by He et al. (He et al., 2020) and compare the predictive
distribution to its GP posterior counterpart based on the NNGP implementation
in the _neural_tangents_ (Novak et al., 2020) package. The GP posterior has
zero mean and kernel function $K_{\text{NNGP}}K^{(2)}$, where
$K_{\text{NNGP}}$ is the corresponding NNGP kernel of $f_{\text{NN}}$ and
$K^{(2)}$ is the chosen kernel function for the mapping $g$ in the ICK _y_
ensemble. Here, $K^{(2)}$ is set as an _exponential-sine-squared_ kernel with
period $T=2\pi$. As shown in Figure 2, the predictive distribution of ICK _y_
ensemble is very close to the analytic GP posterior, demonstrating the
_approximate_ equivalence between the two models (see Appendix L for more
evidence). Alternative strategies can also be used to approximate a GP
posterior with a multiplicative kernel between NTK and $K^{(2)}$.

Algorithm 2 ICK _y_ Ensemble

Input: data
$\boldsymbol{X}=\left\\{\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{i}^{(2)}\right\\}_{i=1}^{N}$,
targets $\boldsymbol{y}=[y_{i}]_{i=1}^{N}$, ensemble
$F_{\text{ICK}y}=\left\\{f_{1},...,f_{N_{e}}\right\\}$ where each function in
$F$ consists of $f_{\text{NN}},K^{(2)}$ with parameters
$\boldsymbol{\theta}^{(1)},\boldsymbol{\theta}^{(2)}$, respectively

for $s=1,...,N_{e}$ do

Apply proper initialization strategy to $f_{s}$

Perform Algorithm 1 on function $f_{s}$

end for

Return:
Predictive mean
$\hat{\boldsymbol{\mu}}=\frac{1}{N_{e}}\sum_{s=1}^{N_{e}}f_{s}(\boldsymbol{X})$
Predictive variance
$\hat{\boldsymbol{\sigma}}^{2}=\frac{1}{N_{e}}\sum_{s=1}^{N_{e}}\left(f_{s}(\boldsymbol{X})-\hat{\boldsymbol{\mu}}\right)^{2}$
Updated parameters $\boldsymbol{\theta}^{(1)},\boldsymbol{\theta}^{(2)}$ for
all functions $f$ in $F_{\text{ICK}y}$

##  5 Experimental Results

We evaluate ICK _y_ on 3 different data sets: a synthetic data set 111Code for
synthetic and worker productivity data:
https://anonymous.4open.science/r/ICK_NNGP-FF83/., a remote sensing data set
222Code for remote sensing and ICK _y_ ensemble:
https://anonymous.4open.science/r/ICK-8F87/., and a worker productivity data
set 111Code for synthetic and worker productivity data:
https://anonymous.4open.science/r/ICK_NNGP-FF83/. obtained from UCI Machine
Learning Repository (Dua & Graff, 2017). In all 3 experiments, our ICK _y_
framework only consists of 2 kernels (i.e. $M=2$), one NN-implied kernel and
one chosen kernel function with trainable parameters. To verify that ICK _y_
can work with more than 2 kernels, we create another synthetic data set with 3
kernels and show the results in Appendix F. The implementation details of all
the experiments in this section are provided in Appendix G and the data
accessibility and restrictions are provided in Appendix H. Besides regression
tasks, we also elaborate how ICK can adapt to classification tasks in Appendix
I.

Figure 2: Predictive distribution from ICK _y_ ensemble and its GP posterior
counterpart on a 1D regression task

###  5.1 Synthetic Data

To verify that ICK _y_ can simulate sampling from a GP with _multiplicative
kernel_ , we create a synthetic data set $y\sim\mathcal{GP}(0,K_{1}K_{2})$
containing 3000 data points where $\boldsymbol{x}^{(1)}\in[0,1]$ is the input
for the linear kernel $K_{1}$ and $\boldsymbol{x}^{(2)}\in[0,2]$ is the input
for the _spectral mixture_ kernel (Wilson & Adams, 2013) $K_{2}$ with 2
components. The data set is first randomly shuffled and then divided into
train and test set with a 50:50 ratio. With ICK _y_ , we process
$\boldsymbol{x}^{(1)}$ with a single-hidden-layer NN and
$\boldsymbol{x}^{(2)}$ with a spectral mixture kernel function.

We then compare ICK _y_ with two models: a multi-layer perceptron (MLP)
applied to the concatenated features and a novel multi-layer perceptron-random
forest (MLP-RF) joint model employed by Zheng et al. (2021), where MLP learns
from $\boldsymbol{x}^{(1)}$ and RF learns from $\boldsymbol{x}^{(2)}$. We
believe MLP-RF serves as a good benchmark model as it is a joint model with
similar architecture to our ICK _y_ framework. To see how ICK _y_ simulates
the spectral mixture kernel, we plot only $\boldsymbol{x}^{(2)}$ against the
predicted value of $y$ as shown in Figure 3. As can be seen from the figure,
plain MLP only captures the linear trend. MLP-RF only captures the mean of the
spectral mixture components. In contrast, our ICK _y_ framework captures both
the mean and the variance of the spectral mixture kernel. To test the
robustness of ICK _y_ , we also conduct experiments on another synthetic data
set in Appendix J to confirm that ICK _y_ can simulate an _additive kernel_.

1.22

Figure 3: Prediction of $y\sim\mathcal{GP}(0,K_{1}K_{2})$, where
$\boldsymbol{x}^{(1)}$ is input to a linear kernel $K_{1}$ and
$\boldsymbol{x}^{(2)}$ is input to a spectral mixture kernel $K_{2}$. We plot
$\boldsymbol{x}^{(2)}$ against the predicted $y$. We show results from a plain
MLP (left), MLP-RF (middle left), and ICK _y_ framework (middle right), and we
compare to the true values of $y$ (right).

###  5.2 Remote Sensing Data

We believe ICK _y_ will be particularly useful for remote sensing
applications. In this experiment, we collect remote sensing data from 51 air
quality monitoring (AQM) stations located in the National Capital Territory
(NCT) of Delhi and its satellite cities over the period from January 1, 2018
to June 30, 2020 (see Appendix H for notes on data availability). Each data
point $\boldsymbol{x}=\\{x,t\\}$ contains 2 sources of information: a three-
band natural color (red-blue-green) satellite image $x$ as the high-
dimensional information and the corresponding timestamp as the low-dimensional
information. Note that we convert the timestamps into numerical values $t$
(where the day 2018-01-01 corresponds to $t=0$) before feeding them into the
models. Our goal is to predict the ground-level PM2.5 concentration
$\hat{y}=f(x,t)$ using both sources of information.

We split the train, validation, and test data set based on $t$. Specifically,
we use all the data points with $t<365$ for training, $365\leq t<500$ for
validation, and $t\geq 500$ for testing. As PM2.5 varies with time on a yearly
basis, we use an _exponential-sine-squared_ kernel with a period of $T=365$
(days) to process the low-dimensional information $t$. The satellite images
are processed with a CNN. The results of ICK _y_ are then compared with 2
benchmarks as shown in Table 1: a Convolutional Neural Network-Random Forest
(CNN-RF) joint model (Zheng et al., 2021, 2020) (similar to the MLP-RF model
in Section 5.1, where RF learns the temporal variation of PM2.5 and CNN
captures the spatial variation of PM2.5 from satellite images) and a
_carefully designed_ CNN-RF model that maps $t$ into two new features,
$\sin(2\pi t/365)$ and $\cos(2\pi t/365)$, to explicitly model seasonality,
which we call Seasonal CNN-RF.

We note that the inner product operation in ICK is similar in mathematical
structure to attention-based mechanisms (Vaswani et al., 2017), which are
popular in many deep learning frameworks. Therefore, we introduce 2 variants
of ICK _y_ and compare them with another 4 benchmarks where the CNN is
replaced by an attention-based mechanism based off of a Vision Transformer
(ViT) (Dosovitskiy et al., 2020) and DeepVit (Zhou et al., 2021)
architectures. These models include ViT-RF, Seasonal ViT-RF, Seasonal DeepViT-
RF, and Seasonal MAE-ViT-RF, where ViT is pre-trained by a Masked Autoencoder
(He et al., 2022). “Seasonal” here denotes that we use the transformed
periodic representation defined above. We note that we are unaware of Vision
Transformers being used in this manner, and we want to evaluate whether it is
the model structure or the prior information that is improving the results.
Hence, these attention-based models, while not a primary contribution, are
novel and represent a good-faith effort to define models with similar forms.
As displayed in Table 1, both standard CNN-RF and ViT-RF models yields very
large errors on predicting PM2.5 values. After including seasonality, CNN-RF
performs significantly better and shares similar predictive performance with
CNN-ICK _y_. After replacing CNN with ViT, we observe slight improvement in
both RF-joint models (especially when pre-trained by MAE) and ICK _y_
variants. Among all the models we present here, ViT-ICK _y_ and DeepViT-ICK
_y_ achieve the highest correlation and the smallest prediction error,
respectively. To evaluate the uncertainty calibration capability of ICK _y_ ,
we construct ICK _y_ ensembles by following Algorithm 2 and compare them with
ensemble formulations of RF-joint model benchmarks. We use an evaluation
criterion called Mean Standardized Log Loss (MSLL) (Williams & Rasmussen,
2006),

$\textstyle\text{MSLL}=\frac{1}{2N}\sum_{i=1}^{N}\left(\log\left(2\pi\hat{\sigma}_{i}^{2}\right)+\frac{(y-\hat{\mu}_{i})^{2}}{\hat{\sigma}_{i}^{2}}\right),$
(15)

where $\hat{\mu}_{i}$ and $\hat{\sigma}_{i}^{2}$ represents the predictive
mean and variance for the $i^{th}$ data point, respectively. From Table 1, it
can be observed that CNN-ICK _y_ ensemble achieves the smallest MSLL. In
addition, we realize ViT-based ICK _y_ variants yield very large MSLL as they
make predictions which are far from the true PM2.5 labels with high confidence
(i.e. small variance). We provide details and possible explanations in
Appendix K.

Table 1: Correlation and error statistics of ICK _y_ and other joint deep models with both convolutional and attention-based architectures on the PM2.5 forecasting task. “S.” denotes seasonal variants. | R${}_{\text{Spear}}$ | RMSE | MAE | MSLL
---|---|---|---|---
CNN-RF | 0.00 | 194.63 | 185.83 | -
ViT-RF | 0.07 | 190.82 | 181.63 | -
| S. CNN-RF
---
0.62 | 53.36 | 39.38 | 96.77
| S. ViT-RF
---
0.66 | 56.45 | 41.73 | 14.69
| S. Deep-ViT-RF
---
0.65 | 56.36 | 42.46 | 17.63
| S. MAE-ViT-RF
---
0.67 | 53.87 | 40.78 | 31.09
CNN-ICK _y_ | 0.62 | 53.46 | 39.76 | 10.92
ViT-ICK _y_ | 0.68 | 56.56 | 41.41 | 12208
DeepViT-ICK _y_ | 0.66 | 52.41 | 35.93 | 38220

###  5.3 Worker Productivity Data

To see if our ICK _y_ framework generalizes to other domains, we acquire
another data set containing the normalized productivity and corresponding
features of garment workers from the UCI Machine Learning Repository. Imran et
al. (Al Imran et al., 2019) employed a dense MLP with 2 hidden layers to
predict the worker productivity with collected features such as date, team
number, targeted productivity, etc. To test our ICK _y_ framework, we separate
out the _date_ and use it as the low-dimensional information. The rest of the
features (excluding the temporal information) are then concatenated together
to serve as the high-dimensional information. Observing that the _daily
averaged_ worker productivity has an approximate _monthly trend_ , we again
use an _exponential-sine-squared_ kernel. The network architecture of ICK _y_
is the same as that of the two-hidden-layer MLP benchmark. To demonstrate the
strength of ICK _y_ compared to other methods that equip NNs with GPs, we also
add 2 additional benchmarks here: a Gaussian Neural Process (GNP) (Bruinsma et
al., 2021) and an Attentive Gaussian Neural Process (AGNP) (Kim et al., 2019).
Based on the results shown in Table 2, ICK _y_ has the best performance when
the period parameter of the kernel is set to be $T=30$ (days) and it
outperforms both MLP and NP benchmarks by almost one order of magnitude. When
we set $T=2$ or $T=7$, this improvement is less significant, which aligns with
our initial observation that the daily averaged productivity has a monthly
seasonal trend. It is also worth noting that the GNP benchmarks here yield
larger errors than MLPs. A possible explanation is that GNP does not allow
explicit assignment of a stationary kernel (as the kernel models a posterior
covariance) so it is hard for GNP to identify specific patterns in the data
such as seasonality without being given the pattern a priori.

Table 2: Prediction error of actual worker productivity on the test data set with ICK _y_ and other benchmark models (MLPs and NPs) | MSE $\downarrow(*10^{-3})$ | MAE $\downarrow(*10^{-2})$
---|---|---
MLP | 20.16 $\pm$ 1.26 | 9.93 $\pm$ 0.36
Cyclic MLP | 20.97 $\pm$ 1.98 | 10.16 $\pm$ 0.77
GNP | 57.25 $\pm$ 4.31 | 19.39 $\pm$ 0.94
AGNP | 43.11 $\pm$ 5.95 | 14.38 $\pm$ 0.88
ICK _y_ , $T=2$ | 3.43 $\pm$ 1.42 | 4.85 $\pm$ 1.00
ICK _y_ , $T=7$ | 0.44 $\pm$ 0.13 | 1.43 $\pm$ 0.15
ICK _y_ , $T=30$ | 0.31 $\pm$ 0.09 | 1.17 $\pm$ 0.14

##  6 Discussion

##### Efficiency and Flexibility of ICK

Compared to exact composite GP models which scale $\mathcal{O}(N^{3})$, the
training process of our ICK framework is more efficient as it leverages
standard backpropagation to learn both the paramters of NN and the kernel
function. Besides efficiency, our ICK framework is more flexible compared to
other joint models (i.e. BNNs and CNN-RF). To be specific, the BNNs
implemented by Pearce et al. (Pearce et al., 2020) cannot simulate complicated
kernels such as the spectral mixture kernel we use in Section 5.1. The CNN-RF
joint model implemented by Zheng et al. (Zheng et al., 2021) requires us to
carefully design the input pre-processing procedure. There is a potential
concern that ICKy may run into computational challenges when a large number of
inducing points are required. This was not a problem in our experiments, but
in large scale models this could be tackled by considering conjugate gradient
methods, which have been recently popular in GP inference (Gardner et al.,
2018).

##### Limitations

A limitation of ICK lies in our method of combining latent representations as
the nature of inner product (i.e. the effect of multiplying small numbers) may
cause _vanishing gradient_ problems when we have a large number of sources of
information (i.e. $M$ is large). Furthermore, in some cases, the properties of
predictive posterior given by ICK ensemble are dominated by the neural network
as shown in Appendix K. It can sometimes be challenging to choose an
appropriate NN architecture to make sure it does not interfere with other
specified kernels.

##### Broader Impacts

We believe our framework is extensively applicable to regression problems in
many fields of study involving high-dimensional data and multiple sources of
information with perceptible trends, such as remote sensing, spatial
statistics, or clinical diagnosis. Also, we are not aware of any negative
societal impacts of our work.

##  7 Conclusion

This paper presents a novel yet surprisingly simple Implicit Composite Kernel
(ICK) framework to learn from _hybrid_ data containing both high-dimensional
information and low-dimensional information with prior knowledge. We first
analytically show the resemblance between ICK and composite GPR models and
then conduct experiments using both synthetic and real-world data. It appears
that ICK outperforms various benchmark models in our experiments with lowest
prediction errors and highest correlations even with very limited data.
Overall, we show that our ICK framework is exceptionally powerful when
learning from _hybrid_ data with prior knowledge incorporated, and we hope our
work can inspire more future research on joint machine learning models,
enhancing their performance, efficiency, flexibility, and generalization
capability.

## References

  * Al Imran et al. (2019) Al Imran, A., Amin, M. N., Rifat, M. R. I., and Mehreen, S.  Deep neural network approach for predicting the productivity of garment employees.  In _2019 6th International Conference on Control, Decision and Information Technologies (CoDIT)_ , pp. 1402–1407. IEEE, 2019.
  * Bishop & Nasrabadi (2006) Bishop, C. M. and Nasrabadi, N. M.  _Pattern recognition and machine learning_ , volume 4.  Springer, 2006.
  * Bottou et al. (2018) Bottou, L., Curtis, F. E., and Nocedal, J.  Optimization methods for large-scale machine learning.  _Siam Review_ , 60(2):223–311, 2018.
  * Bouhlel & Martins (2019) Bouhlel, M. A. and Martins, J. R.  Gradient-enhanced kriging for high-dimensional problems.  _Engineering with Computers_ , 35(1):157–173, 2019.
  * Bouhlel et al. (2016) Bouhlel, M. A., Bartoli, N., Otsmane, A., and Morlier, J.  Improving kriging surrogates of high-dimensional design models by partial least squares dimension reduction.  _Structural and Multidisciplinary Optimization_ , 53(5):935–952, 2016.
  * Bruinsma et al. (2021) Bruinsma, W. P., Requeima, J., Foong, A. Y., Gordon, J., and Turner, R. E.  The gaussian neural process.  _arXiv preprint arXiv:2101.03606_ , 2021.
  * Choudhuri et al. (2007) Choudhuri, N., Ghosal, S., and Roy, A.  Nonparametric binary regression using a gaussian process prior.  _Statistical Methodology_ , 4(2):227–243, 2007\.
  * Datta et al. (2016) Datta, A., Banerjee, S., Finley, A. O., and Gelfand, A. E.  Hierarchical nearest-neighbor gaussian process models for large geostatistical datasets.  _Journal of the American Statistical Association_ , 111(514):800–812, 2016.
  * Dosovitskiy et al. (2020) Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  _arXiv preprint arXiv:2010.11929_ , 2020.
  * Drineas et al. (2005) Drineas, P., Mahoney, M. W., and Cristianini, N.  On the nyström method for approximating a gram matrix for improved kernel-based learning.  _journal of machine learning research_ , 6(12), 2005.
  * Dua & Graff (2017) Dua, D. and Graff, C.  UCI machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Duvenaud (2014) Duvenaud, D.  _Automatic model construction with Gaussian processes_.  PhD thesis, University of Cambridge, 2014.
  * Feydy (2020) Feydy, J.  _Geometric data analysis, beyond convolutions_.  PhD thesis, Université Paris-Saclay Gif-sur-Yvette, France, 2020.
  * Gardner et al. (2018) Gardner, J., Pleiss, G., Weinberger, K. Q., Bindel, D., and Wilson, A. G.  Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration.  _Advances in neural information processing systems_ , 31, 2018.
  * Garnelo et al. (2018) Garnelo, M., Schwarz, J., Rosenbaum, D., Viola, F., Rezende, D. J., Eslami, S., and Teh, Y. W.  Neural processes.  _arXiv preprint arXiv:1807.01622_ , 2018.
  * Garriga-Alonso et al. (2018) Garriga-Alonso, A., Rasmussen, C. E., and Aitchison, L.  Deep convolutional networks as shallow gaussian processes.  _arXiv preprint arXiv:1808.05587_ , 2018.
  * Gelfand & Schliep (2016) Gelfand, A. E. and Schliep, E. M.  Spatial statistics and gaussian processes: A beautiful marriage.  _Spatial Statistics_ , 18:86–104, 2016.
  * Girolami & Rogers (2006) Girolami, M. and Rogers, S.  Variational bayesian multinomial probit regression with gaussian process priors.  _Neural Computation_ , 18(8):1790–1817, 2006\.
  * Goodfellow et al. (2016) Goodfellow, I., Bengio, Y., and Courville, A.  _Deep learning_.  MIT press, 2016.
  * He et al. (2020) He, B., Lakshminarayanan, B., and Teh, Y. W.  Bayesian deep ensembles via the neural tangent kernel.  _Advances in neural information processing systems_ , 33:1010–1022, 2020.
  * He et al. (2022) He, K., Chen, X., Xie, S., Li, Y., Dollár, P., and Girshick, R.  Masked autoencoders are scalable vision learners.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 16000–16009, 2022.
  * Hensman et al. (2013) Hensman, J., Fusi, N., and Lawrence, N. D.  Gaussian processes for big data.  _arXiv preprint arXiv:1309.6835_ , 2013.
  * Jacot et al. (2018) Jacot, A., Gabriel, F., and Hongler, C.  Neural tangent kernel: Convergence and generalization in neural networks.  _Advances in neural information processing systems_ , 31, 2018.
  * Jaderberg et al. (2015) Jaderberg, M., Simonyan, K., Zisserman, A., et al.  Spatial transformer networks.  _Advances in neural information processing systems_ , 28, 2015.
  * Jiang et al. (2022) Jiang, Z., Zheng, T., Bergin, M., and Carlson, D.  Improving spatial variation of ground-level pm2. 5 prediction with contrastive learning from satellite imagery.  _Science of Remote Sensing_ , pp. 100052, 2022.
  * Journel & Huijbregts (1976) Journel, A. G. and Huijbregts, C. J.  _Mining geostatistics_.  The Blackburn Press, 1976.
  * Kim et al. (2019) Kim, H., Mnih, A., Schwarz, J., Garnelo, M., Eslami, A., Rosenbaum, D., Vinyals, O., and Teh, Y. W.  Attentive neural processes.  _arXiv preprint arXiv:1901.05761_ , 2019.
  * Kim et al. (2005) Kim, H.-M., Mallick, B. K., and Holmes, C. C.  Analyzing nonstationary spatial data using piecewise gaussian processes.  _Journal of the American Statistical Association_ , 100(470):653–668, 2005.
  * Kingma & Ba (2014) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Krige (1951) Krige, D. G.  A statistical approach to some basic mine valuation problems on the witwatersrand.  _Journal of the Southern African Institute of Mining and Metallurgy_ , 52(6):119–139, 1951.
  * Lagaris et al. (1998) Lagaris, I. E., Likas, A., and Fotiadis, D. I.  Artificial neural networks for solving ordinary and partial differential equations.  _IEEE transactions on neural networks_ , 9(5):987–1000, 1998.
  * LeCun et al. (2015) LeCun, Y., Bengio, Y., and Hinton, G.  Deep learning.  _nature_ , 521(7553):436–444, 2015.
  * Lee et al. (2017) Lee, J., Bahri, Y., Novak, R., Schoenholz, S. S., Pennington, J., and Sohl-Dickstein, J.  Deep neural networks as gaussian processes.  _arXiv preprint arXiv:1711.00165_ , 2017.
  * Lee et al. (2019) Lee, J., Xiao, L., Schoenholz, S., Bahri, Y., Novak, R., Sohl-Dickstein, J., and Pennington, J.  Wide neural networks of any depth evolve as linear models under gradient descent.  _Advances in neural information processing systems_ , 32, 2019.
  * Maddison et al. (2016) Maddison, C. J., Mnih, A., and Teh, Y. W.  The concrete distribution: A continuous relaxation of discrete random variables.  _arXiv preprint arXiv:1611.00712_ , 2016.
  * Marcus (2018) Marcus, G.  Deep learning: A critical appraisal.  _arXiv preprint arXiv:1801.00631_ , 2018.
  * Markou et al. (2021) Markou, S., Requeima, J., Bruinsma, W., and Turner, R.  Efficient gaussian neural processes for regression.  _arXiv preprint arXiv:2108.09676_ , 2021.
  * Matsubara et al. (2020) Matsubara, T., Oates, C. J., and Briol, F.-X.  The ridgelet prior: A covariance function approach to prior specification for bayesian neural networks.  _arXiv preprint arXiv:2010.08488_ , 2020.
  * Matthews et al. (2017) Matthews, A. G. d. G., Hron, J., Turner, R. E., and Ghahramani, Z.  Sample-then-optimize posterior sampling for bayesian linear models.  In _NeurIPS Workshop on Advances in Approximate Bayesian Inference_ , 2017.
  * Matthews et al. (2018) Matthews, A. G. d. G., Rowland, M., Hron, J., Turner, R. E., and Ghahramani, Z.  Gaussian process behaviour in wide deep neural networks.  _arXiv preprint arXiv:1804.11271_ , 2018.
  * McFee et al. (2011) McFee, B., Lanckriet, G., and Jebara, T.  Learning multi-modal similarity.  _Journal of machine learning research_ , 12(2), 2011.
  * Moseley et al. (2020) Moseley, B., Markham, A., and Nissen-Meyer, T.  Solving the wave equation with physics-informed deep learning.  _arXiv preprint arXiv:2006.11894_ , 2020.
  * Neal (1996) Neal, R. M.  Priors for infinite networks.  In _Bayesian Learning for Neural Networks_ , pp. 29–53. Springer, 1996.
  * Novak et al. (2020) Novak, R., Xiao, L., Hron, J., Lee, J., Alemi, A. A., Sohl-Dickstein, J., and Schoenholz, S. S.  Neural tangents: Fast and easy infinite neural networks in python.  In _International Conference on Learning Representations_ , 2020.  URL https://github.com/google/neural-tangents.
  * Pearce et al. (2020) Pearce, T., Tsuchida, R., Zaki, M., Brintrup, A., and Neely, A.  Expressive priors in bayesian neural networks: Kernel combinations and periodic functions.  In _Uncertainty in artificial intelligence_ , pp. 134–144. PMLR, 2020.
  * Petelin et al. (2013) Petelin, D., Grancharova, A., and Kocijan, J.  Evolving gaussian process models for prediction of ozone concentration in the air.  _Simulation modelling practice and theory_ , 33:68–80, 2013\.
  * Rahimi & Recht (2007) Rahimi, A. and Recht, B.  Random features for large-scale kernel machines.  _Advances in neural information processing systems_ , 20, 2007.
  * Rahimi & Recht (2008) Rahimi, A. and Recht, B.  Weighted sums of random kitchen sinks: Replacing minimization with randomization in learning.  _Advances in neural information processing systems_ , 21, 2008.
  * Shi & Horvath (2006) Shi, T. and Horvath, S.  Unsupervised learning with random forest predictors.  _Journal of Computational and Graphical Statistics_ , 15(1):118–138, 2006.
  * Snelson & Ghahramani (2005) Snelson, E. and Ghahramani, Z.  Sparse gaussian processes using pseudo-inputs.  _Advances in neural information processing systems_ , 18, 2005.
  * Van der Wilk et al. (2017) Van der Wilk, M., Rasmussen, C. E., and Hensman, J.  Convolutional gaussian processes.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Vaswani et al. (2017) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., and Polosukhin, I.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.
  * Wang et al. (2020) Wang, H., He, H., and Katabi, D.  Continuously indexed domain adaptation.  _International Conference on Machine Learning_ , 2020.
  * Williams & Seeger (2000) Williams, C. and Seeger, M.  Using the nyström method to speed up kernel machines.  _Advances in neural information processing systems_ , 13, 2000.
  * Williams & Barber (1998) Williams, C. K. and Barber, D.  Bayesian classification with gaussian processes.  _IEEE Transactions on pattern analysis and machine intelligence_ , 20(12):1342–1351, 1998.
  * Williams & Rasmussen (2006) Williams, C. K. and Rasmussen, C. E.  _Gaussian processes for machine learning_ , volume 2.  MIT press Cambridge, MA, 2006.
  * Wilson & Adams (2013) Wilson, A. and Adams, R.  Gaussian process kernels for pattern discovery and extrapolation.  In _International conference on machine learning_ , pp. 1067–1075. PMLR, 2013.
  * Wilson et al. (2011) Wilson, A. G., Knowles, D. A., and Ghahramani, Z.  Gaussian process regression networks.  _arXiv preprint arXiv:1110.4411_ , 2011.
  * Wilson et al. (2016) Wilson, A. G., Hu, Z., Salakhutdinov, R., and Xing, E. P.  Deep kernel learning.  In _Artificial intelligence and statistics_ , pp. 370–378. PMLR, 2016.
  * Yang et al. (2012) Yang, T., Li, Y.-F., Mahdavi, M., Jin, R., and Zhou, Z.-H.  Nyström method vs random fourier features: A theoretical and empirical comparison.  _Advances in neural information processing systems_ , 25, 2012.
  * Zhang et al. (2011) Zhang, D., Wang, Y., Zhou, L., Yuan, H., Shen, D., Initiative, A. D. N., et al.  Multimodal classification of alzheimer’s disease and mild cognitive impairment.  _Neuroimage_ , 55(3):856–867, 2011.
  * Zheng et al. (2020) Zheng, T., Bergin, M. H., Hu, S., Miller, J., and Carlson, D. E.  Estimating ground-level pm2. 5 using micro-satellite images by a convolutional neural network and random forest approach.  _Atmospheric Environment_ , 230:117451, 2020.
  * Zheng et al. (2021) Zheng, T., Bergin, M., Wang, G., and Carlson, D.  Local pm2. 5 hotspot detector at 300 m resolution: A random forest–convolutional neural network joint model jointly trained on satellite images and meteorology.  _Remote Sensing_ , 13(7):1356, 2021.
  * Zhou et al. (2021) Zhou, D., Kang, B., Jin, X., Yang, L., Lian, X., Jiang, Z., Hou, Q., and Feng, J.  Deepvit: Towards deeper vision transformer.  _arXiv preprint arXiv:2103.11886_ , 2021.

##  Appendix A Further Discussion on Theorem 1

###  A.1 A More Detailed Proof for Theorem 1

In this section, we provide a more detailed proof for Theorem 1. Again,
suppose we have data
$\boldsymbol{X}=\left[\boldsymbol{x}_{i}\right]_{i=1}^{N}$ with information
from 2 different sources
$\boldsymbol{x}=\left\\{\boldsymbol{x}^{(1)},\boldsymbol{x}^{(2)}\right\\}$
where $\boldsymbol{x}^{(1)}\in\mathbb{R}^{D_{1}}$ is high-dimensional and
$\boldsymbol{x}^{(2)}\in\mathbb{R}^{D_{2}}$ is low-dimensional with some known
relationship with the targets $\boldsymbol{Y}=[y_{i}]_{i=1}^{N}$. With ICK
formulation, we have
$\boldsymbol{z}^{(1)}=f_{\text{NN}}\left(\boldsymbol{x}^{(1)}\right)\in\mathbb{R}^{p}$
where $f_{\text{NN}}$ is a NN-implied function with parameters
$\boldsymbol{\theta}^{(1)}$ and
$\boldsymbol{z}^{(2)}=g\left(\boldsymbol{x}^{(2)}\right)\in\mathbb{R}^{p}$
where $g$ is a kernel-to-latent-space mapping specified by a chosen kernel
function $K^{(2)}$ with parameters $\boldsymbol{\theta}^{(2)}$. In Section
3.3, it is stated that the latent representation $\boldsymbol{z}^{(1)}$ from
$f_{\text{NN}}$ will converge in distribution to a multi-output GP in the
infinite width limit. In case when the NN has finite width,
$\boldsymbol{z}^{(1)}$ will _approximately_ follows a GP with empirical NNGP
kernel function $\hat{K}^{\text{NNGP}}$

$\boldsymbol{z}^{(1)}\sim\mathcal{GP}_{\textit{approx}}\left(0,\hat{K}^{\text{NNGP}}\right).$
(16)

Let $\boldsymbol{Z}^{(1)}=\left[\boldsymbol{z}_{i}^{(1)}\right]_{i=1}^{N}$. If
Assumption 4.2 holds, then $\boldsymbol{Z}^{(1)}$ can be approximately treated
as a sample from a multivariate Gaussian as shown below

$\boldsymbol{Z}^{(1)}\sim\mathcal{N}\left(\boldsymbol{0},\hat{\boldsymbol{K}}^{\text{NNGP}}\right),$
(17)

where $\boldsymbol{Z}^{(1)}\in\mathbb{R}^{N\times p}$ and
$\hat{\boldsymbol{K}}^{\text{NNGP}}\in\mathbb{R}^{Np\times Np}$ is the
corresponding kernel matrix

$\hat{\boldsymbol{K}}^{\text{NNGP}}=\begin{bmatrix}\hat{\boldsymbol{K}}^{\text{NNGP}}_{1}&\boldsymbol{0}&\cdots&\boldsymbol{0}\\\
\boldsymbol{0}&\hat{\boldsymbol{K}}^{\text{NNGP}}_{2}&\cdots&\boldsymbol{0}\\\
\vdots&\vdots&\ddots&\vdots\\\
\boldsymbol{0}&\boldsymbol{0}&\cdots&\hat{\boldsymbol{K}}^{\text{NNGP}}_{p}\end{bmatrix},$
(18)

and
$\hat{\boldsymbol{K}}^{\text{NNGP}}_{k}=\left[\hat{K}^{\text{NNGP}}_{k}\left(x^{(1)}_{i},x^{(1)}_{j}\right)\right]_{i,j=1,...,N}\in\mathbb{R}^{N\times
N}$ for all $k=1,2,...,p$. Also, since all entries in $\boldsymbol{z}^{(1)}$
are from the same NN architecture with parameters drawn from the same
distribution, as we increase the NN width, all
$\hat{K}^{\text{NNGP}}_{1},...,\hat{K}^{\text{NNGP}}_{p}$ will converge to a
deterministic limit. Therefore, in the finite-width case, it is reasonable to
say that
$\hat{K}^{\text{NNGP}}_{1}\approx\cdots\approx\hat{K}^{\text{NNGP}}_{p}\approx\hat{K}^{\text{NNGP}}$.
Based on the marginalization property of GPs (Williams & Rasmussen, 2006), any
$\boldsymbol{z}_{i}^{(1)}\in\mathbb{R}^{p}$ in $\boldsymbol{Z}^{(1)}$ should
_approximately_ satisfy

$\boldsymbol{z}_{i}^{(1)}\sim\mathcal{N}\left(\boldsymbol{0},\hat{K}^{\text{NNGP}}\left(\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{i}^{(1)}\right)\boldsymbol{I}_{p}\right),$
(19)

where $\boldsymbol{I}_{p}$ is a $p\times p$ identity matrix. In other words,
each entry in $\boldsymbol{z}_{i}^{(1)}$ has a univariate Gaussian
distribution

$z_{ik}^{(1)}\sim\mathcal{N}\left(0,\hat{K}^{\text{NNGP}}\left(\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{i}^{(1)}\right)\right),k=1,2,...,p.$
(20)

Let $\alpha_{ik}=z_{ik}^{(2)}=g\left(\boldsymbol{x}_{i}^{(2)}\right)_{k}$,
since the final output is
$\hat{y}=f_{\text{ICK}}\left(\boldsymbol{x}^{(1)},\boldsymbol{x}^{(2)}\right)={\boldsymbol{z}^{(1)}}^{T}\boldsymbol{z}^{(2)}$,
for the $i^{th}$ data point, $\hat{y}_{i}$ can be viewed as a weighted sum of
$p$ independent Gaussian random variables under Assumption 4.2

$\hat{y}_{i}\sim\sum_{k=1}^{p}\alpha_{ik}\mathcal{N}\left(0,\hat{K}^{\text{NNGP}}\left(\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{i}^{(1)}\right)\right)=\sum_{k=1}^{p}\mathcal{N}\left(0,\alpha_{ik}^{2}\hat{K}^{\text{NNGP}}\left(\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{i}^{(1)}\right)\right)=\mathcal{N}\left(0,\hat{K}^{\text{NNGP}}\left(\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{i}^{(1)}\right)\sum_{k=1}^{p}\alpha_{ik}^{2}\right).$
(21)

Therefore, the joint distribution of the final predictions of the whole
training set $\hat{\boldsymbol{Y}}=[\hat{y}_{i}]_{i=1}^{N}$ will be a
multivariate Gaussian as given below

$\displaystyle\hat{\boldsymbol{Y}}$
$\displaystyle\sim\mathcal{N}\left(\boldsymbol{0},\sum_{k=1}^{p}\begin{bmatrix}\alpha_{1k}^{2}\hat{K}^{\text{NNGP}}\left(\boldsymbol{x}_{1}^{(1)},\boldsymbol{x}_{1}^{(1)}\right)&\alpha_{1k}\alpha_{2k}\hat{K}^{\text{NNGP}}\left(\boldsymbol{x}_{1}^{(1)},\boldsymbol{x}_{2}^{(1)}\right)&\cdots&\alpha_{1k}\alpha_{Nk}\hat{K}^{\text{NNGP}}\left(\boldsymbol{x}_{1}^{(1)},\boldsymbol{x}_{N}^{(1)}\right)\\\
\alpha_{2k}\alpha_{1k}\hat{K}^{\text{NNGP}}\left(\boldsymbol{x}_{2}^{(1)},\boldsymbol{x}_{1}^{(1)}\right)&\alpha_{2k}^{2}\hat{K}^{\text{NNGP}}\left(\boldsymbol{x}_{2}^{(1)},\boldsymbol{x}_{2}^{(1)}\right)&\cdots&\alpha_{2k}\alpha_{Nk}\hat{K}^{\text{NNGP}}\left(\boldsymbol{x}_{2}^{(1)},\boldsymbol{x}_{N}^{(1)}\right)\\\
\vdots&\vdots&\ddots&\vdots\\\
\alpha_{Nk}\alpha_{1k}\hat{K}^{\text{NNGP}}\left(\boldsymbol{x}_{N}^{(1)},\boldsymbol{x}_{1}^{(1)}\right)&\alpha_{Nk}\alpha_{2k}\hat{K}^{\text{NNGP}}\left(\boldsymbol{x}_{N}^{(1)},\boldsymbol{x}_{2}^{(1)}\right)&\cdots&\alpha_{Nk}^{2}\hat{K}^{\text{NNGP}}\left(\boldsymbol{x}_{N}^{(1)},\boldsymbol{x}_{N}^{(1)}\right)\end{bmatrix}\right),$
(22)

which can be simplified as

$\hat{\boldsymbol{Y}}\sim\mathcal{N}\left(\boldsymbol{0},\sum_{k=1}^{p}\boldsymbol{\alpha}_{k}\boldsymbol{\alpha}_{k}^{T}*\hat{\boldsymbol{K}}^{\text{NNGP}}\right)=\mathcal{N}\left(\boldsymbol{0},\boldsymbol{\alpha}\boldsymbol{\alpha}^{T}\odot\hat{\boldsymbol{K}}^{\text{NNGP}}\right),$
(23)

where $\boldsymbol{\alpha}_{k}=[\alpha_{ik}]_{i=1}^{N}\in\mathbb{R}^{N}$,
$\boldsymbol{\alpha}=[\alpha_{ik}]_{\begin{subarray}{c}i=1,...,N\\\
k=1,...,p\end{subarray}}\in\mathbb{R}^{N\times p}$,
$\hat{\boldsymbol{K}}^{\text{NNGP}}=\left[\hat{K}^{\text{NNGP}}\left(\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{j}^{(1)}\right)\right]_{i,j=1,...,N}$,
and ”$\odot$” represents elementwise multiplication. Therefore, if $g$
includes the kernel-to-latent-space mapping
$K^{(2)}\left(\boldsymbol{x}_{i}^{(2)},\boldsymbol{x}_{j}^{(2)}\right)\approx{\boldsymbol{z}_{i}^{(2)}}^{T}\boldsymbol{z}_{j}^{(2)}=g\left(\boldsymbol{x}_{i}^{(2)}\right)^{T}g\left(\boldsymbol{x}_{j}^{(2)}\right),i,j=1,...,N$
as shown in Equation 5, then
$\boldsymbol{K}^{(2)}=\boldsymbol{\alpha}\boldsymbol{\alpha}^{T}$ and we
derive that

$\hat{\boldsymbol{Y}}\sim\mathcal{N}\left(\boldsymbol{0},\hat{\boldsymbol{K}}^{\text{NNGP}}\odot\boldsymbol{K}^{(2)}\right).$
(24)

Since this derivation reasonably holds for any finite set of $\boldsymbol{X}$
and $\boldsymbol{Y}$, we say the ICK framework is _approximately_ equivalent
to a GP with zero mean and a multiplicative kernel between the NNGP kernel
matrix $\hat{\boldsymbol{K}}^{\text{NNGP}}$ and the user-specified kernel
matrix $\boldsymbol{K}^{(2)}$.

###  A.2 A More Detailed Discussion on Assumption 4.2

As stated in Section 4.1 and the previous section, Theorem 4.1 relies heavily
on Assumption 4.2 which assumes that all entries in $\boldsymbol{z}^{(1)}$
(the representation extracted by the neural network) are independent with each
other. To show the validity of this assumption, we can explicitly write the
$k^{th}$ entry in $\boldsymbol{z}^{(1)}$ as follows:

$z_{k}^{(1)}(\cdot)=b_{k}^{(1)}+\sum_{r=1}^{N_{z^{(1)}}}W_{kr}^{(1)}\phi\left({z^{\prime}}^{(1)}_{r}(\cdot)\right),$
(25)

where $W_{kr}^{(1)}$ and $b_{k}^{(1)}$ are i.i.d. Gaussian weight and bias
parameters and ${z^{\prime}}^{(1)}_{r}$ is the $r^{th}$ entry in the previous
layer. We see that $z_{k}^{(1)}$ is a linear combination of i.i.d. Gaussian
terms, and we can show that $z_{k}^{(1)}(x_{i})$ and $z_{l}^{(1)}(x_{j})$ are
independent for two different inputs $x_{i}$ and $x_{j}$ and $k\neq l$ by
computing the covariance as shown below.

$\displaystyle\text{cov}\left(z_{k}^{(1)}(x_{i}),z_{l}^{(1)}(x_{j})\right)$
$\displaystyle=\mathbb{E}\left[\left(b_{k}^{(1)}+\sum_{r=1}^{N_{z^{(1)}}}W_{kr}^{(1)}\phi\left({z^{\prime}}^{(1)}_{r}(x_{i})\right)\right)\left(b_{l}^{(1)}+\sum_{s=1}^{N_{z^{(1)}}}W_{ls}^{(1)}\phi\left({z^{\prime}}^{(1)}_{s}(x_{j})\right)\right)\right]$
$\displaystyle=\mathbb{E}\left[b_{k}^{(1)}b_{l}^{(1)}\right]+\mathbb{E}\left[b_{k}^{(1)}\sum_{s=1}^{N_{z^{(1)}}}W_{ls}^{(1)}\phi\left({z^{\prime}}^{(1)}_{s}(x_{j})\right)\right]+\mathbb{E}\left[b_{l}^{(1)}\sum_{r=1}^{N_{z^{(1)}}}W_{kr}^{(1)}\phi\left({z^{\prime}}^{(1)}_{r}(x_{i})\right)\right]$
$\displaystyle+\mathbb{E}\left[\left(\sum_{r=1}^{N_{z^{(1)}}}W_{kr}^{(1)}\phi\left({z^{\prime}}^{(1)}_{r}(x_{i})\right)\right)\left(\sum_{s=1}^{N_{z^{(1)}}}W_{ls}^{(1)}\phi\left({z^{\prime}}^{(1)}_{s}(x_{j})\right)\right)\right]$

Since the weights and biases are all i.i.d. Gaussian with zero mean _a priori_
and they are independent of the outputs from the previous layer, we have:

$\displaystyle\mathbb{E}\left[b_{k}^{(1)}b_{l}^{(1)}\right]=\mathbb{E}\left[b_{k}^{(1)}\right]\mathbb{E}\left[b_{l}^{(1)}\right]=0$
$\displaystyle\mathbb{E}\left[b_{k}^{(1)}\sum_{s=1}^{N_{z^{(1)}}}W_{ls}^{(1)}\phi\left({z^{\prime}}^{(1)}_{s}(x_{j})\right)\right]=\mathbb{E}\left[b_{k}^{(1)}\right]\sum_{s=1}^{N_{z^{(1)}}}\mathbb{E}\left[W_{ls}^{(1)}\phi\left({z^{\prime}}^{(1)}_{s}(x_{j})\right)\right]=0$
$\displaystyle\mathbb{E}\left[b_{l}^{(1)}\sum_{r=1}^{N_{z^{(1)}}}W_{kr}^{(1)}\phi\left({z^{\prime}}^{(1)}_{r}(x_{i})\right)\right]=\mathbb{E}\left[b_{l}^{(1)}\right]\sum_{r=1}^{N_{z^{(1)}}}\mathbb{E}\left[W_{kr}^{(1)}\phi\left({z^{\prime}}^{(1)}_{r}(x_{i})\right)\right]=0$
$\displaystyle\mathbb{E}\left[\left(\sum_{r=1}^{N_{z^{(1)}}}W_{kr}^{(1)}\phi\left({z^{\prime}}^{(1)}_{r}(x_{i})\right)\right)\left(\sum_{s=1}^{N_{z^{(1)}}}W_{ls}^{(1)}\phi\left({z^{\prime}}^{(1)}_{s}(x_{j})\right)\right)\right]$
$\displaystyle=\sum_{r=1}^{N_{z^{(1)}}}\sum_{s=1}^{N_{z^{(1)}}}\mathbb{E}\left[W_{kr}^{(1)}W_{ls}^{(1)}\phi\left({z^{\prime}}^{(1)}_{r}(x_{i})\right)\phi\left({z^{\prime}}^{(1)}_{s}(x_{j})\right)\right]$
$\displaystyle=\sum_{r=1}^{N_{z^{(1)}}}\sum_{s=1}^{N_{z^{(1)}}}\mathbb{E}\left[W_{kr}^{(1)}\right]\mathbb{E}\left[W_{ls}^{(1)}\right]\mathbb{E}\left[\phi\left({z^{\prime}}^{(1)}_{r}(x_{i})\right)\phi\left({z^{\prime}}^{(1)}_{s}(x_{j})\right)\right]=0.$

Therefore, we have
$\text{cov}\left(z_{k}^{(1)}(x_{i}),z_{l}^{(1)}(x_{j})\right)=0$, meaning that
$z_{ik}^{(1)}$ and $z_{jl}^{(1)}$ are uncorrelated with each other _a priori_
, which shows the validity of Assumption 4.2. However, we do not claim the
validity of this assumption _a posteriori_.

##  Appendix B Random Fourier Features

###  B.1 Methodology

Random Fourier Features (RFF) is another popular approximation method used for
kernel learning (Rahimi & Recht, 2007). Unlike the Nyström method which
approximates the entire kernel matrix, RFF directly approximates the kernel
function $K$ using some randomized feature mapping
$\phi:\mathbb{R}^{D_{m}}\rightarrow\mathbb{R}^{2d_{m}}$ such that
$K\left(x_{i}^{(m)},x_{j}^{(m)}\right)\approx\phi\left(x_{i}^{(m)}\right)^{T}\phi\left(x_{j}^{(m)}\right)$.
To obtain the feature mapping $\phi$, based on Bochner’s theorem, we first
compute the Fourier transform $p$ of kernel $K$

$p(\omega)=\frac{1}{(2\pi)^{D_{m}}}\int_{-\infty}^{+\infty}e^{-j\omega^{T}\delta}K(\delta)d\delta,$
(26)

where $\delta=x_{i}^{(m)}-x_{j}^{(m)}$. Then we draw $d_{m}$ i.i.d. samples
$\omega_{1},\omega_{2},...,\omega_{d_{m}}$ from $p(\omega)$ and construct the
feature mapping $\phi$ as follows

$\phi\left(x^{(m)}\right)\equiv\textstyle
d_{m}^{-1/2}\left[\cos\left(\omega_{1}^{T}x^{(m)}\right),...,\cos\left(\omega_{d_{m}}^{T}x^{(m)}\right),\sin\left(\omega_{1}^{T}x^{(m)}\right),...,\sin\left(\omega_{d_{m}}^{T}x^{(m)}\right)\right].$
(27)

Since $\phi\left(x^{(m)}\right)\in\mathbb{R}^{2d_{m}}$, we need to set
$d_{m}=p/2$ when using RFF as a kernel-to-latent-space transformation. In
addition, since RFF involves sampling from a distribution, the kernel
parameters are thus not directly differentiable and we need to apply a
reparameterization trick (Maddison et al., 2016) to learn those parameters.

(a) $\hat{y}=f_{\text{NN}}\left(x^{(1)}\right)$

(b) $\hat{y}=f_{\text{ICK\emph{r}}}\left(x^{(1)},x^{(2)}\right)$

(c) $\hat{y}=f_{\text{ICK\emph{r}}}\left(x^{(1)},x^{(2)},x^{(3)}\right)$

Figure B1: Scatter plots of the true values of $y$ against the predicted
values of $y$ using our ICK _r_ framework with (a) one source of information,
(b) 2 sources of information, and (c) 3 sources of information. Note that here
we use RFF for kernel-to-latent-space transformation.

###  B.2 Experimental Results

####  B.2.1 Synthetic Data

We use the same toy data set where each data point
$\boldsymbol{x}=\left\\{x^{(1)},x^{(2)},x^{(3)}\right\\}$ contains 3 sources
of information as described in Appendix F. Also, we use the same types of
kernels as those in ICK _y_ as discussed in Appendix F. The only difference
here is that we use RFF instead of Nyström method to transform the kernel
matrix into the latent space in ICK _r_ framework.

The results are displayed in Figure B1. It can be observed that when we add in
only the side information $x^{(2)}$ along with the _exponential sine squared_
kernel, both the correlation and the predictive performance are improved
(though not as good as the results from ICK _y_ as shown in Figure F1).
However, after we further include $x^{(3)}$ with the _RBF_ kernel, we realize
that the parameters of ICK _r_ become very hard to optimize and it fails to
make valid predictions and starts to guess randomly around zero.

(a) ICK _y_

(b) ICK _r_

Figure B2: Density plots of the true PM2.5 concentrations against the
forecasted PM2.5 concentrations for $t\geq 500$ using our ICK framework with
(a) ICK _y_ and (b) ICK _r_

####  B.2.2 Remote Sensing Data

We also try ICK _r_ on the _forecasting_ task using the remote sensing data
(see Section 5.2) and compare the results with those from ICK _y_. Each data
point $\boldsymbol{x}=\\{x,t\\}$ contains a satellite image $x$ as the high-
dimensional information and its corresponding timestamp $t$ as the low-
dimensional information. The satellite images are processed with a two-layer
CNN and the timestamps are processed with an _exponential-sine-squared_ kernel
with a period of $T=365$ (days). As can be observed from Figure B2, ICK _r_
yields much higher error compared to ICK _y_.

##  Appendix C Estimated Kernel Matrix and its Eigen-spectrum

(a) $\boldsymbol{K}_{\text{true}}$

(b) $\boldsymbol{K}_{\text{ICK\emph{y}}}$

(c)
$\left|\boldsymbol{K}_{\text{true}}-\boldsymbol{K}_{\text{ICK\emph{y}}}\right|$

Figure C1: Visualization of (a) True matrix (b) estimated matrix by our ICK
_y_ framework, and (c) absolute difference between the true and estimated
matrix for the spectral mixture kernel Figure C2: Eigenvalues of the kernel
matrix computed from the first 4 batches of training data

We first examine whether ICK _y_ can retrieve the spectral mixture kernel in
the prediction task. After fitting the parameters of the spectral mixture
kernel in ICK _y_ , we compute the kernel matrix
$\boldsymbol{K}_{\text{ICK\emph{y}}}$ using these learned parameters and
compare it with the true kernel matrix $\boldsymbol{K}_{\text{true}}$ by
calculating the absolute difference between them as displayed in Figure C1. As
can be observed, $\boldsymbol{K}_{\text{ICK\emph{y}}}$ and
$\boldsymbol{K}_{\text{true}}$ are similar and their absolute difference is
relatively small, indicating that ICK _y_ can approximately retrieve the
spectral mixture kernel.

Yang et al. (Yang et al., 2012) studied the fundamental difference between
Nyström method and Random Fourier Features (RFF). They conclude that Nyström-
method-based approaches can yield much better generalization error bound than
RFF-based approaches if there exists a large gap in the eigen-spectrum of the
kernel matrix. This phenomenon is mainly caused by how these two methods
construct their basis functions. In particular, the basis functions used by
RFF are sampled from a Gaussian distribution that is independent from the
training examples, while the basis functions used by the Nyström method are
sampled from the training samples so they are data-dependent. In our synthetic
data experiments, we train our ICK framework using a batch size of 50. The
eigenvalues of the kernel matrices computed from the first 4 batches of the
synthetic data set are displayed in Figure C2. It can be observed that the
first few eigenvalues of the kernel matrix are much larger than the remaining
eigenvalues. Namely, there exists a large gap in the eigen-spectrum of the
kernel matrix, which helps explain why ICK _y_ has a much better performance
than ICK _r_.

##  Appendix D Number of Inducing Points

As discussed in Section 4.2.1, as we increase the number of inducing points
$p$, we expect the approximation error between the true kernel matrix
$\boldsymbol{K}$ and the approximated kernel matrix $\hat{\boldsymbol{K}}$ to
decrease. Here, we empirically show how the value of $p$ impacts our
predictions. In Figure 1(a), we plot the prediction error of
$\hat{y}=f_{\text{ICK\emph{y}}}\left(x^{(1)},x^{(2)},x^{(3)}\right)$ against
the number of inducing points using the synthetic data generated in Appendix
F. As can be observed, the prediction error drops sharply as we raise $p$ from
a small value (e.g. $p=2$). When $p$ is relatively large, increasing $p$
yields smaller improvement on the predictions. Additionally, in Figure 1(b),
we plot the total training time against $p$. The total training time is
dependent on how long a single iteration takes and the total number of epochs
required. We note that once $p>80$ the training time is relatively flat, which
is due to the fact that the total computation in the Cholesky is less than the
computation in the neural network. Interestingly, it appears that when $p$ is
very small, ICK _y_ takes longer to converge due to the need for many more
epochs. As we increase $p$, the training time goes down and then goes up again
due to the computational complexity, i.e. $\mathcal{O}(p^{3})$, of the
Cholesky decomposition. Based on these observations, we are not concerned
about the computational complexity for reasonable values of $p$.

(a) Prediction error

(b) Training time

Figure D1: Plots of (a) prediction error and (b) training time of
$\hat{y}=f_{\text{ICK\emph{y}}}\left(x^{(1)},x^{(2)},x^{(3)}\right)$ against
the number of inducing points $p$

##  Appendix E Applying Sample-then-optimize Procedure to ICK

As elaborated in Section 4.3, a deep ensemble with proper initialization
scheme will have a GP posterior interpretation in the infinite width limit
when trained by a squared-error loss (He et al., 2020), which is an example of
the ”sample-then-optimize” procedure of Matthews et al. (Matthews et al.,
2017) Algorithm 2 is an alternative to replace $f_{\text{NN}}$ with such a
deep ensemble $F=\left\\{f_{n_{e}}\right\\}_{n_{e}=1}^{N_{e}}$ where the
dimension of the final readout layer of each $f_{n_{e}}$ is $p$. In this case,
each baselearner $f_{n_{e}}$ can be viewed as an i.i.d. sample from a multi-
output GP in the infinite width limit. In the finite-width case, this
relationship becomes approximate:

$f_{n_{e}}\sim\mathcal{GP}_{\textit{approx}}\left(0,K_{F}\right),$ (28)

As stated by Lee et al. (Lee et al., 2019), if all the parameters in
$f_{n_{e}}$ are randomly drawn from a Gaussian distribution and are all fixed
except the last layer, then after training $F$ on a squared-error loss, we
will have $K_{F}\rightarrow K^{\text{NNGP}}$ where $K^{\text{NNGP}}$ is the
NNGP kernel. He et el. (He et al., 2020) also proposed to add a random and
untrainable function $\delta(\cdot)$ to the output of $f_{n_{e}}$ so we have
$K_{F}\rightarrow K^{\text{NTK}}$ where $K^{\text{NTK}}$ is the limiting NTK.

With ICK formulation, by following the proof given in Appendix A, we can
derive that, for each baselearner in the ICK ensemble, the joint distribution
of the final predictions of the whole training set
$\hat{\boldsymbol{Y}}_{n_{e}}$ will again be a multivariate Gaussian as shown
below if Assumption 4.2 holds _a posteriori_ :

$\hat{\boldsymbol{Y}}_{n_{e}}\sim\mathcal{N}\left(\boldsymbol{0},\boldsymbol{K}_{F}\odot\boldsymbol{K}^{(2)}\right),$
(29)

where
$\left(\boldsymbol{K}_{F}\right)_{ij}=K_{F}\left(\boldsymbol{x}_{i}^{(1)},\boldsymbol{x}_{j}^{(1)}\right)$.
Since Equation 29 holds for any finite input data set, we can conclude that
each baselearner $f_{s}$ in the ICK ensemble
$F_{\text{ICK}y}=\\{f_{s}\\}_{s=1}^{N_{e}}$ can be _approximately_ viewed as
an i.i.d. sample from a single-output GP

$f_{s}\sim\mathcal{GP}_{\textit{approx}}\left(0,K_{F}K^{(2)}\right)=\mathcal{GP}_{\textit{approx}}\left(0,K^{\text{comp}}\right).$
(30)

Building an ensemble $F_{\text{ICK}y}=\\{f_{s}\\}_{s=1}^{N_{e}}$ is thus
equivalent to performing a Monte Carlo approximation to a GP predictive
posterior distribution whose mean and covariance matrix for a test data set
$\boldsymbol{X}^{*}=\left[\boldsymbol{x}_{i}^{*}\right]_{i=1}^{N^{*}}$ are

$\displaystyle f_{s}(\boldsymbol{X}^{*}|\boldsymbol{X})$
$\displaystyle\sim\mathcal{N}\left(\boldsymbol{\mu}^{*},\boldsymbol{\Sigma}^{*}\right),$
(31) $\displaystyle\boldsymbol{\mu}^{*}$
$\displaystyle=\boldsymbol{K}^{\text{comp}}_{X^{*}X}\left(\boldsymbol{K}^{\text{comp}}_{XX}\right)^{-1}\boldsymbol{Y},$
(32) $\displaystyle\boldsymbol{\Sigma}^{*}$
$\displaystyle=\boldsymbol{K}^{\text{comp}}_{X^{*}X^{*}}-\boldsymbol{K}^{\text{comp}}_{X^{*}X}\left(\boldsymbol{K}^{\text{comp}}_{XX}\right)^{-1}\boldsymbol{K}^{\text{comp}}_{XX^{*}},$
(33)

where
$\left(\boldsymbol{K}^{\text{comp}}_{X^{*}X}\right)_{ij}=K^{\text{comp}}\left(\boldsymbol{x}_{i}^{*},\boldsymbol{x}_{j}\right)$,
$\left(\boldsymbol{K}^{\text{comp}}_{XX}\right)_{ij}=K^{\text{comp}}\left(\boldsymbol{x}_{i},\boldsymbol{x}_{j}\right)$,
$\left(\boldsymbol{K}^{\text{comp}}_{X^{*}X^{*}}\right)_{ij}=K^{\text{comp}}\left(\boldsymbol{x}_{i}^{*},\boldsymbol{x}_{j}^{*}\right)$,
and
$\boldsymbol{K}^{\text{comp}}_{XX^{*}}=\left(\boldsymbol{K}^{\text{comp}}_{X^{*}X}\right)^{T}$.
In other words, we can approximate the GP predictive posterior using the
predictive mean and variance generated by Algorithm 2 as follows:

$\displaystyle\hat{\boldsymbol{\mu}}^{*}$
$\displaystyle=\frac{1}{N_{e}}\sum_{n=1}^{N_{e}}f_{s}(\boldsymbol{X}^{*}|\boldsymbol{X})\approx\mathbb{E}\left[f_{s}(\boldsymbol{X}^{*}|\boldsymbol{X})\right]=\boldsymbol{\mu}^{*},$
(34) $\displaystyle{\hat{\boldsymbol{\sigma}}}^{*2}$
$\displaystyle=\frac{1}{N_{e}}\sum_{n=1}^{N_{e}}\left[f_{s}(\boldsymbol{X}^{*}|\boldsymbol{X})-\hat{\boldsymbol{\mu}}^{*}\right]^{2}\approx\mathbb{V}\left[f_{s}(\boldsymbol{X}^{*}|\boldsymbol{X})\right]=\text{diag}\left(\boldsymbol{\Sigma}^{*}\right).$
(35)

##  Appendix F ICK with More Than Two Kernels

(a) $\hat{y}=f_{\text{NN}}\left(x^{(1)}\right)$

(b) $\hat{y}=f_{\text{ICK\emph{y}}}\left(x^{(1)},x^{(2)}\right)$

(c) $\hat{y}=f_{\text{ICK\emph{y}}}\left(x^{(1)},x^{(2)},x^{(3)}\right)$

Figure F1: Scatter plots of the true values of $y$ against the predicted
values of $y$ using our ICK _y_ framework with (a) one source, (b) 2 sources,
and (c) 3 sources of information

Besides the visualization presented in Figure 1, we also show our ICK
framework for processing data
$\boldsymbol{x}=\left\\{x^{(1)},x^{(2)},...,x^{(M)}\right\\}$ with $M>2$
sources of information in Figure F2. Here $K^{(2)},...,K^{(M)}$ represent
different types of kernels with trainable parameters. The final prediction is
calculated by a chained inner product of all extracted representations
$\hat{y}=\sum_{k=1}^{p}\prod_{m=1}^{M}z_{k}^{(m)}$.

Figure F2: Given data containing $M$ sources of information
$\boldsymbol{x}=\left\\{x^{(1)},x^{(2)},...,x^{(M)}\right\\}$, we can process
the data using our ICK framework where high-dimensional information (e.g.
$x^{(1)}$ in the figure) is processed using a neural network and low-
dimensional information (e.g. $x^{(2)}$ in the figure) is processed using a
kernel function followed by Nyström or RFF transformation.

To confirm that ICK _y_ can work with more than 2 kernels, we construct
another synthetic data set containing 3000 data points in total. Each input
$\boldsymbol{x}=\left\\{x^{(1)},x^{(2)},x^{(3)}\right\\}$ has 3 sources of
information. The output $y$ is generated by
$y=x^{(3)}\tanh{\left(2x^{(1)}\cos^{2}\left(\pi
x^{(2)}/50\right)\right)}+\epsilon$ where $\epsilon$ is a Gaussian noise term.
We process $x^{(1)}$ with a small single-hidden-layer NN, $x^{(2)}$ with an
_exponential sine squared_ kernel, and $x^{(3)}$ with a _radial-basis
function_ (RBF) kernel. Figure F1 shows the prediction results as we
progressively add more sources of information into our ICK _y_ framework with
corresponding kernel functions. It can be observed that ICK _y_ yields both
smallest error and highest correlation with information from all 3 different
sources. Hence, ICK _y_ works well with the $M=3$ case and the regression
performance is improved as we add in more information related to the target.

##  Appendix G Experimental Details

###  G.1 Synthetic Data

We use the GPytorch package (Gardner et al., 2018) to generate the synthetic
data. Before feeding $x^{(1)}$ into MLP, we first map $x^{(1)}$ into higher
dimension using an unsupervised algorithm called Totally Random Trees
Embedding (Shi & Horvath, 2006). All the MLP structures in this experiment
(including those in MLP-RF and ICK _y_) contain one single fully connected
(FC) layer of width 1000, which serves as a simple benchmark since a one-
hidden-layer MLP can only capture linear relationship between the input and
output. For model training, we optimize a Mean Squared Error (MSE) objective
using Adam optimizer (Kingma & Ba, 2014) with a weight decay of 0.1.

###  G.2 Remote Sensing Data

The model architecture and training details are listed in Table 3. Here $p$
denotes the length of latent representations $\boldsymbol{z}$ as discussed in
Section 4 and $d_{\text{RT}}$ denotes the transformed dimension of timestamp
$t$ using the Random Trees Embedding as mentioned in Section G.1. Note that we
use stochastic gradient descent (SGD) optimizer with a momentum of 0.9 for ICK
_y_ as we realize that SGD helps ICK _y_ find a local minimum on the objective
more efficiently. We use MSE objective for ICK _y_ and all benchmark models in
this experiment.

Table 3: Model architecture and training details for remote sensing data experiment in Section 5.2 |  | Backbone architecture
---
details
| Output FC layers
---
dimension
Optimizer
CNN-RF |  | # Conv blocks = 2, # Channels = 16,
---
Kernel size = 3, Stride = 1
1000 + $d_{\text{RT}}$, 512, 512, 1 |  | Adam
---
$\beta_{1}=0.9$
$\beta_{2}=0.999$
ViT-RF |  | # Transformer blocks = 2,
---
# Attention heads = 8,
Dropout ratio = 0.1
1000 + $d_{\text{RT}}$, 512, 512, 1 |  | Adam
---
$\beta_{1}=0.9$
$\beta_{2}=0.999$
DeepViT-RF |  | # Transformer blocks = 2,
---
# Attention heads = 8,
Dropout ratio = 0.1
1000 + $d_{\text{RT}}$, 512, 512, 1 |  | Adam
---
$\beta_{1}=0.9$
$\beta_{2}=0.999$
MAE-ViT-RF |  | # Transformer blocks = 2,
---
# Attention heads = 8,
Dropout ratio = 0.1,
Masking ratio = 0.75
1000 + $d_{\text{RT}}$, 512, 512, 1 |  | Adam
---
$\beta_{1}=0.9$
$\beta_{2}=0.999$
CNN-ICKy |  | # Conv blocks = 2, # Channels = 16,
---
Kernel size = 3, Stride = 1
1000, 512, $p$ |  | SGD
---
momentum = 0.9
ViT-ICKy |  | # Transformer blocks = 2,
---
# Attention heads = 8,
Dropout ratio = 0.1
1000, 512, $p$ |  | SGD
---
momentum = 0.9
DeepViT-ICKy |  | # Transformer blocks = 2,
---
# Attention heads = 8,
Dropout ratio = 0.1
1000, 512, $p$ |  | SGD
---
momentum = 0.9

###  G.3 Worker Productivity Data

The MLPs (including the MLP part in ICK _y_) in this experiment share the same
structure as the one used in (Al Imran et al., 2019), which consist of 3
hidden layers of width 128, 32, and 32, respectively. For plain MLP, cyclic
MLP, and ICK _y_ , we use the mean absolute error (MAE) objective to put less
weight on the outliers and thus enhance the model performance. For GNP and
AGNP, we maximize a biased Monte Carlo estimate of the log-likelihood
objective as discussed in (Markou et al., 2021). All these objectives are
optimized by an Adam optimizer with $\beta_{1}=0.9$ and $\beta_{2}=0.999$.

##  Appendix H Generation, Accessibility, and Restrictions of the Data

All experiments are conducted on a computer cluster equipped with a GeForce
RTX 2080 Ti GPU. The synthetic data $y\sim\mathcal{GP}(0,K_{1}K_{2})$ in
Section 5.1 and $y\sim\mathcal{GP}(0,K_{1}+K_{2})$ in Appendix J are generated
using the GPyTorch package. The remote sensing data in Section 5.2 is
downloaded using PlanetScope API whose content is protected by copyright
and/or other intellectual property laws. To access the data on PlanetScope,
the purchase of an end-user license is required. When this manuscript is
accepted, we will provide the codes we used to acquire the data. The UCI
machine learning repository data we use in Section 5.3 has an open access
license, meaning that the data is freely available online.

##  Appendix I Adapting ICK for Classification

While regression tasks are the primary motivation for this paper, there are
many ways to adapt GPR for classification tasks. For example, a binary
classification model can be created by using a sigmoid (Williams & Barber,
1998) or probit link (Choudhuri et al., 2007) on the output of the GP.
Succinctly, given a function
$f(\boldsymbol{x})\sim\mathcal{GP}\left(0,K(\boldsymbol{x},\boldsymbol{x}^{\prime})\right)$,
the binary outcome probability is be given as $p(y=1|f(x))=\sigma(f(x))$.
Likewise, a multiple classification model can be constructed by using a multi-
output GP (or multiple GPs) and putting the outputs through a softmax function
(Williams & Barber, 1998) or multinomial probit link (Girolami & Rogers,
2006). This strategy can be summarized by calculating $C$ different functions
$f_{c}(\boldsymbol{x})\sim\mathcal{GP}\left(0,K(\boldsymbol{x},\boldsymbol{x}^{\prime})\right)$
for $c=1,...,C$, where $C$ is the number of classes, and then calculating the
class probabilities through a link function,
$p(y|\boldsymbol{x})=\text{softmax}([f_{1}(\boldsymbol{x}),f_{2}(\boldsymbol{x}),...,f_{C}(\boldsymbol{x})])$.

This same logic can be used to construct a multiple classification model from
ICK _y_. Succinctly, let
$r_{c}=f_{\text{NN},c}(\boldsymbol{x}^{(1)})\odot\boldsymbol{z}^{(2)}_{c}$,
where $f_{\text{NN},c}$ denotes a neural network specific to the $c^{th}$
class and $z^{(2)}_{c}$ represents the Nyström approximation specific to the
kernel for the $c^{th}$ class. We note that often in a multi-output case the
kernel parameters are shared, and so $\boldsymbol{z}^{(2)}_{c}$ would be an
identical vector for each class. Then, the output probabilities for a data
sample as $p(y|\boldsymbol{x})=\text{softmax}([r_{1},\dots,r_{C}])$. This
framework is learned with a cross-entropy loss.

To provide proof-of-concept of this multiple classification strategy, we
implemented this model on a version of Rotating MNIST. In this task, a dataset
was created by rotating each image in the dataset by a uniform random value
$\phi\in[0,2\pi)$, thus creating a dataset with 60,000 images each with an
associated rotation covariate $\phi$. We implemented the above multiple
classification model with a periodic kernel over the rotation angle. This
strategy yielded an accuracy of 92.3% on the validation data. This is lower
than methods such as spatial transformers (Jaderberg et al., 2015) that report
accuracy greater than 99%. However, those models explicitly use the fact that
the information is simply rotated, whereas ICK is modeling a smooth
transformation in the prediction function as a function of angle. This ICK
classification model is much closer in concept to the way Rotating MNIST is
used to evaluate unsupervised domain adaptation. While the evaluation strategy
is different than our random validation set, the state-of-the-art accuracy on
unsupervised domain adaption is 87.1% (Wang et al., 2020). Due to the lack of
complete and fair comparisons, we are not claiming that ICK _y_ is state-of-
the-art for classification, but ICK _y_ ’s classification model does seem
reasonable and viable based upon this result.

##  Appendix J Additional Results on Synthetic Data

(a) $y=\mathcal{GP}(0,K_{1}K_{2})$

(b) $y=\mathcal{GP}(0,K_{1}+K_{2})$

Figure J1: Prediction error of plain MLP, MLP-RF, and ICK _y_ with different
amount of training data generated by (a) $y=\mathcal{GP}(0,K_{1}K_{2})$ and
(b) $y=\mathcal{GP}(0,K_{1}+K_{2})$

To demonstrate the superior performance of ICK _y_ more clearly, we also test
plain MLP, MLP-RF, and ICK _y_ on the prediction task in Section 5.1 using
different number of training samples. As displayed in Figure 1(a), ICK _y_
yields the smallest error among all the 3 frameworks even with very limited
data.

While ICK is theoretically designed to simulate posterior sampling from GPs
with multiplicative kernels, to show its flexibility, we also evaluate how
well it could simulate additive kernels. We conduct experiments using another
synthetic data set generated by an additive kernel
$y\sim\mathcal{GP}(0,K_{1}+K_{2})$ with the same training settings. As shown
in Figure J2, ICK _y_ again outperforms plain MLP and MLP-RF. Moreover, we
test plain MLP, MLP-RF, and ICK _y_ on this task using different number of
training samples. As displayed in Figure 1(b), ICK _y_ yields the smallest
error among all the 3 frameworks. Also, the performance gap between ICK _y_
and the other 2 benchmark models shrink as we feed in more training data.
Therefore, we conclude ICK _y_ is robust enough to simulate posterior sampling
from GPs with both additive and multiplicative kernels.

1.22

Figure J2: Prediction of $y\sim\mathcal{GP}(0,K_{1}+K_{2})$, where
$\boldsymbol{x}^{(1)}$ is input to a linear kernel $K_{1}$ and
$\boldsymbol{x}^{(2)}$ is input to a spectral mixture kernel $K_{2}$. We plot
$\boldsymbol{x}^{(2)}$ against the predicted $y$. We show results from a plain
MLP (left), MLP-RF (middle left), and ICK _y_ framework (middle right), and we
compare to the true values of $y$ (right).

##  Appendix K Time Series Visualization for Remote Sensing Experiment

(a) CNN-ICK _y_ (1.92M parameters)

(b) DeepViT-ICK _y_ (21.80M parameters)

Figure K1: Time series visualization of predictive mean and uncertainty of
PM2.5 in Section 5.2 (remote sensing experiment) for (a) CNN-ICK _y_ and (b)
DeepViT-ICK _y_

To explore the problem of posterior uncertainty calibration for ViT variants
of ICK _y_ (i.e. extremely large MSLL as shown in Table 1), we visualize the
results as time series by first grouping the predictions by the timestamp $t$,
taking the _minimum_ of the predictive variance, and plotting them along with
the corresponding true values and the predictive mean as shown in Figure K1.
The shaded region represents the confidence interval
$[\mu-2\sigma,\mu+2\sigma]$ where $\mu$ and $\sigma$ are the predictive mean
and standard deviation of PM2.5, respectively. We realize that CNN-ICK _y_
ensemble tends to yield much higher variance than DeepViT-ICK _y_ ensemble
when the predictive mean deviates from the true values. A plausible
explanation is that the DeepViT structure contains much more parameters than
the CNN structure (21.80M vs 1.92M) used in our remote sensing experiment,
which makes DeepViT-ICK _y_ overparameterized (The DeepViT architecture is set
to be consistent with the CNN architecture as shown in Table 3). To alleviate
this problem, we try reducing the number of transformer blocks in ViT and
DeepViT and we do observe a significant drop in MSLL.

##  Appendix L Distance Between ICK Ensemble Posterior and its NNGP
Counterpart

Figure L1: Wasserstein-1 distance between the predictive posterior of ICK _y_
ensemble and its analytic NNGP counterpart

Besides Figure 2 in Section 4.3, we estimate the Wasserstein-1 distance
between the the predictive posterior of ICK _y_ ensemble and its analytic NNGP
counterpart using the de-biased Sinkhorn divergence (Feydy, 2020). As
displayed in Figure L1, the Wasserstein-1 distance shows a decreasing trend as
we increase the number of estimators $N_{e}$ in the ICK _y_ ensemble (except
at $N_{e}=50$ where there might be some noise or perturbations in the sampling
process from the posterior). This serves as a more convincing evidence,
verifying the predictive posterior of ICK _y_ ensemble will converge to its
analytic NNGP counterpart. However, since we are using finite-width NN for
each baselearner in the ICK _y_ ensemble, there will be a gap between these
two distributions and the Wasserstein-1 distance is not guaranteed to converge
to zero.