Kernel methods for vector output

From the regularization perspective, the problem is to learn ${\displaystyle f_{*}}$ belonging to a reproducing kernel Hilbert space of vector-valued functions (${\displaystyle {\mathcal {H}}}$). This is similar to the scalar case of Tikhonov regularization, with some extra care in the notation.

	Vector-valued case	Scalar case
Reproducing kernel	${\displaystyle \mathbf {K}$ :{\mathcal {X}}\times {\mathcal {X}}\rightarrow \mathbb {R} ^{D\times D}}	${\displaystyle k:{\mathcal {X}}\times {\mathcal {X}}\rightarrow \mathbb {R} }$
Learning problem	${\displaystyle f_{*}=\operatorname {argmin} \sum \limits _{j=1}^{D}{\frac {1}{N}}\sum \limits _{i=1}^{N}(f_{j}(\mathbf {x_{i}} )-y_{j,i})^{2}+\lambda \Vert \mathbf {f} \Vert _{\mathbf {K} }^{2}}$	${\displaystyle f_{*}=\operatorname {argmin} {\frac {1}{N}}\sum \limits _{i=1}^{N}(f(\mathbf {x_{i}} )-y_{i})^{2}+\lambda \Vert \mathbf {f} \Vert _{k}^{2}}$
Solution (derived via the representer theorem${\displaystyle ^{\dagger }}$)	${\displaystyle f_{*}(\mathbf {x} )=\sum \limits _{i=1}^{N}\mathbf {K} (\mathbf {x_{i}} ,\mathbf {x} )c_{i}}$ with ${\displaystyle {\bar {\mathbf {c} }}=(\mathbf {K} (\mathbf {X} ,\mathbf {X} )+\lambda N\mathbf {(} I))^{-1}{\bar {\mathbf {y} }}}$, where ${\displaystyle {\bar {\mathbf {c} }}{\text{ and }}{\bar {\mathbf {y} }}}$ are the coefficients and output vectors concatenated to form ${\displaystyle ND}$ vectors and ${\displaystyle \mathbf {K} (\mathbf {X} ,\mathbf {X} ){\text{ is an }}ND\times ND}$ matrix of ${\displaystyle N\times N}$ blocks: ${\displaystyle (\mathbf {K} (\mathbf {x_{i}} ,\mathbf {x_{j}} ))_{d,d'}}$	${\displaystyle f_{*}(\mathbf {x} )=\sum \limits _{i=1}^{N}k(\mathbf {x_{i}} ,\mathbf {x} )c_{i}=\mathbf {k} _{\mathbf {x} }^{\intercal }\mathbf {c} }$ Solve for ${\displaystyle \mathbf {c} }$ by taking the derivative of the learning problem, setting it equal to zero, and substituting in the above expression for ${\displaystyle f_{*}}$: ${\displaystyle \mathbf {c} =(\mathbf {K} +\lambda I)^{-1}\mathbf {y} }$ where ${\displaystyle \mathbf {K} _{ij}=k(\mathbf {x_{i}} ,\mathbf {x_{j}} )=i^{\text{th}}{\text{ element of }}\mathbf {k} _{\mathbf {x_{j}} }}$

${\displaystyle ^{\dagger }}$It is possible, though non-trivial, to show that a representer theorem also holds for Tikhonov regularization in the vector-valued setting.[8]

Note, the matrix-valued kernel ${\displaystyle \mathbf {K} }$ can also be defined by a scalar kernel ${\displaystyle R}$ on the space ${\displaystyle {\mathcal {X}}\times \{1,\ldots ,D\}}$. An isometry exists between the Hilbert spaces associated with these two kernels:

${\displaystyle (\mathbf {K} (x,x'))_{d,d'}=R((x,d),(x',d'))}$

The estimator of the vector-valued regularization framework can also be derived from a Bayesian viewpoint using Gaussian process methods in the case of a finite dimensional Reproducing kernel Hilbert space. The derivation is similar to the scalar-valued case Bayesian interpretation of regularization. The vector-valued function ${\displaystyle {\textbf {f}}}$, consisting of ${\displaystyle D}$ outputs ${\displaystyle \left\{f_{d}\right\}_{d=1}^{D}}$, is assumed to follow a Gaussian process:

${\displaystyle {\textbf {f}}\sim {\mathcal {GP}}({\textbf {m}},{\textbf {K}})}$

where ${\displaystyle {\textbf {m}}:{\mathcal {X}}\to {\textbf {R}}^{D}}$ is now a vector of the mean functions ${\displaystyle \left\{m_{d}({\textbf {x}})\right\}_{d=1}^{D}}$ for the outputs and ${\displaystyle {\textbf {K}}}$ is a positive definite matrix-valued function with entry ${\displaystyle ({\textbf {K}}({\textbf {x}},{\textbf {x}}'))_{d,d'}}$ corresponding to the covariance between the outputs ${\displaystyle f_{d}({\textbf {x}})}$ and ${\displaystyle f_{d'}({\textbf {x}}')}$.

For a set of inputs ${\displaystyle {\textbf {X}}}$, the prior distribution over the vector ${\displaystyle {\textbf {f}}({\textbf {X}})}$ is given by ${\displaystyle {\mathcal {N}}({\textbf {m}}({\textbf {X}}),{\textbf {K}}({\textbf {X}},{\textbf {X}}))}$, where ${\displaystyle {\textbf {m}}({\textbf {X}})}$ is a vector that concatenates the mean vectors associated to the outputs and ${\displaystyle {\textbf {K}}({\textbf {X}},{\textbf {X}})}$ is a block-partitioned matrix. The distribution of the outputs is taken to be Gaussian:

${\displaystyle p({\textbf {y}}\mid {\textbf {f}},{\textbf {x}},\Sigma )={\mathcal {N}}({\textbf {f}}({\textbf {x}}),\Sigma )}$

where ${\displaystyle \Sigma \in {\mathcal {\textbf {R}}}^{D\times D}}$ is a diagonal matrix with elements ${\displaystyle \left\{\sigma _{d}^{2}\right\}_{d=1}^{D}}$ specifying the noise for each output. Using this form for the likelihood, the predictive distribution for a new vector ${\displaystyle {\textbf {x}}_{*}}$ is:

${\displaystyle p({\textbf {f}}({\textbf {x}}_{*})\mid {\textbf {S}},{\textbf {f}},{\textbf {x}}_{*},\phi )={\mathcal {N}}({\textbf {f}}_{*}({\textbf {x}}_{*}),{\textbf {K}}_{*}({\textbf {x}}_{*},{\textbf {x}}_{*}))}$

where ${\displaystyle {\textbf {S}}}$ is the training data, and ${\displaystyle \phi }$ is a set of hyperparameters for ${\displaystyle {\textbf {K}}({\textbf {x}},{\textbf {x}}')}$ and ${\displaystyle \Sigma }$.

Equations for ${\displaystyle {\textbf {f}}_{*}}$ and ${\displaystyle {\textbf {K}}_{*}}$ can then be obtained:

${\displaystyle {\textbf {f}}_{*}({\textbf {x}}_{*})={\textbf {K}}_{{\textbf {x}}_{*}}^{T}({\textbf {K}}({\textbf {X}},{\textbf {X}})+{\boldsymbol {\Sigma }})^{-1}{\bar {\textbf {y}}}}$

${\displaystyle {\textbf {K}}_{*}({\textbf {x}}_{*},{\textbf {x}}_{*})={\textbf {K}}({\textbf {x}}_{*},{\textbf {x}}_{*})-{\textbf {K}}_{{\textbf {x}}_{*}}({\textbf {K}}({\textbf {X}},{\textbf {X}})+{\boldsymbol {\Sigma }})^{-1}{\textbf {K}}_{{\textbf {x}}_{*}}^{T}}$

where ${\displaystyle {\boldsymbol {\Sigma }}=\Sigma \otimes {\textbf {I}}_{N},{\textbf {K}}_{{\textbf {x}}_{*}}\in {\mathcal {\textbf {R}}}^{D\times ND}}$ has entries ${\displaystyle ({\textbf {K}}({\textbf {x}}_{*},{\textbf {x}}_{j}))_{d,d'}}$ for ${\displaystyle j=1,\cdots ,N}$ and ${\displaystyle d,d'=1,\cdots ,D}$. Note that the predictor ${\displaystyle {\textbf {f}}^{*}}$ is identical to the predictor derived in the regularization framework. For non-Gaussian likelihoods different methods such as Laplace approximation and variational methods are needed to approximate the estimators.

A simple, but broadly applicable, class of multi-output kernels can be separated into the product of a kernel on the input-space and a kernel representing the correlations among the outputs:[8]

${\displaystyle (\mathbf {K} (\mathbf {x} ,\mathbf {x'} ))_{d,d'}=k(\mathbf {x} ,\mathbf {x'} )k_{T}(d,d')}$

${\displaystyle k}$: scalar kernel on ${\displaystyle {\mathcal {X}}\times {\mathcal {X}}}$

${\displaystyle k_{T}}$: scalar kernel on ${\displaystyle \{1,\ldots ,D\}\times \{1,\ldots ,D\}}$

In matrix form: ${\displaystyle \mathbf {K} (\mathbf {x} ,\mathbf {x'} )=k(\mathbf {x} ,\mathbf {x'} )\mathbf {B} }$    where ${\displaystyle \mathbf {B} }$ is a ${\displaystyle D\times D}$ symmetric and positive semi-definite matrix. Note, setting ${\displaystyle \mathbf {B} }$ to the identity matrix treats the outputs as unrelated and is equivalent to solving the scalar-output problems separately.

For a slightly more general form, adding several of these kernels yields sum of separable kernels (SoS kernels).

One way of obtaining ${\displaystyle k_{T}}$ is to specify a regularizer which limits the complexity of ${\displaystyle f}$ in a desirable way, and then derive the corresponding kernel. For certain regularizers, this kernel will turn out to be separable.

Mixed-effect regularizer

${\displaystyle R(\mathbf {f} )=A_{\omega }(C_{\omega }\sum \limits _{l=1}^{D}\|f_{l}\|_{k}^{2}+\omega D\sum \limits _{l=1}^{D}\|f_{l}-{\bar {f}}\|_{k}^{2})}$

where:

• ${\displaystyle A_{\omega }={\frac {1}{2(1-\omega )(1-\omega +\omega D)}}}$
• ${\displaystyle C_{\omega }=(2-2\omega +\omega D)}$
• ${\displaystyle {\bar {f}}={\frac {1}{D}}\sum \limits _{q=1}^{D}f_{q}}$
• ${\displaystyle K_{\omega }(x,x')=k(x,x')(\omega \mathbf {1} +(1-\omega )\mathbf {I} _{D}}$

where ${\displaystyle \mathbf {1} {\text{ is a }}D\times D}$ matrix with all entries equal to 1.

This regularizer is a combination of limiting the complexity of each component of the estimator (${\displaystyle f_{l}}$) and forcing each component of the estimator to be close to the mean of all the components. Setting ${\displaystyle \omega =0}$ treats all the components as independent and is the same as solving the scalar problems separately. Setting ${\displaystyle \omega =1}$ assumes all the components are explained by the same function.

Cluster-based regularizer

${\displaystyle R(\mathbf {f} )=\varepsilon _{1}\sum _{c=1}^{r}\sum _{l\in I(c)}\|f_{l}-{\bar {f_{c}}}\|_{k}^{2}+\varepsilon _{2}\sum \limits _{c=1}^{r}m_{c}\|{\bar {f_{c}}}\|_{k}^{2}}$

where:

• ${\displaystyle I(c)}$ is the index set of components that belong to cluster ${\displaystyle c}$
• ${\displaystyle m_{c}}$ is the cardinality of cluster ${\displaystyle c}$
• ${\displaystyle {\bar {f_{c}}}={\frac {1}{m_{c}}}\sum \limits _{q\in I(c)}f_{q}}$
• ${\displaystyle \mathbf {M} _{l,q}={\frac {1}{m_{c}}}}$ if ${\displaystyle l}$ and ${\displaystyle q}$ both belong to cluster ${\displaystyle c}$  (${\displaystyle \mathbf {M} _{l,q}=0}$ otherwise
• ${\displaystyle K(x,x')=k(x,x')\mathbf {G} ^{\dagger }}$

where ${\displaystyle \mathbf {G} _{l,q}=\varepsilon _{1}\delta _{lq}+(\varepsilon _{2}-\varepsilon _{1})\mathbf {M} _{l,q}}$

This regularizer divides the components into ${\displaystyle r}$ clusters and forces the components in each cluster to be similar.

Graph regularizer

${\displaystyle R(\mathbf {f} )={\frac {1}{2}}\sum \limits _{l,q=1}^{D}\Vert f_{l}-f_{q}\Vert _{k}^{2}\mathbf {M} _{lq}+\sum \limits _{l=1}^{D}\Vert f_{l}\Vert _{k}^{2}\mathbf {M} _{l,l}}$

where ${\displaystyle \mathbf {M} {\text{ is a }}D\times D}$ matrix of weights encoding the similarities between the components

${\displaystyle K(x,x')=k(x,x')\mathbf {L} ^{\dagger }}$

where ${\displaystyle \mathbf {L} =\mathbf {D} -\mathbf {M} }$,   ${\displaystyle \mathbf {D} _{l,q}=\delta _{l,q}(\sum \limits _{h=1}^{D}\mathbf {M} _{l,h}+\mathbf {M} _{l,q})}$

Note, ${\displaystyle \mathbf {L} }$ is the graph laplacian. See also: graph kernel.

Several approaches to learning ${\displaystyle \mathbf {B} }$ from data have been proposed.[8] These include: performing a preliminary inference step to estimate ${\displaystyle \mathbf {B} }$ from the training data,[9] a proposal to learn ${\displaystyle \mathbf {B} }$ and ${\displaystyle \mathbf {f} }$ together based on the cluster regularizer,[15] and sparsity-based approaches which assume only a few of the features are needed.[16] [17]

In LMC, outputs are expressed as linear combinations of independent random functions such that the resulting covariance function (over all inputs and outputs) is a valid positive semidefinite function. Assuming ${\displaystyle D}$ outputs ${\displaystyle \left\{f_{d}({\textbf {x}})\right\}_{d=1}^{D}}$ with ${\displaystyle {\textbf {x}}\in {\mathcal {\textbf {R}}}^{p}}$, each ${\displaystyle f_{d}}$ is expressed as:

${\displaystyle f_{d}({\textbf {x}})=\sum _{q=1}^{Q}{a_{d,q}u_{q}({\textbf {x}})}}$

where ${\displaystyle a_{d,q}}$ are scalar coefficients and the independent functions ${\displaystyle u_{q}({\textbf {x}})}$ have zero mean and covariance cov${\displaystyle [u_{q}({\textbf {x}}),u_{q'}({\textbf {x}}')]=k_{q}({\textbf {x}},{\textbf {x}}')}$ if ${\displaystyle q=q'}$ and 0 otherwise. The cross covariance between any two functions ${\displaystyle f_{d}({\textbf {x}})}$ and ${\displaystyle f_{d'}({\textbf {x}})}$ can then be written as:

${\displaystyle \operatorname {cov} [f_{d}({\textbf {x}}),f_{d'}({\textbf {x}}')]=\sum _{q=1}^{Q}{\sum _{i=1}^{R_{q}}{a_{d,q}^{i}a_{d',q}^{i}k_{q}({\textbf {x}},{\textbf {x}}')}}=\sum _{q=1}^{Q}{b_{d,d'}^{q}k_{q}({\textbf {x}},{\textbf {x}}')}}$

where the functions ${\displaystyle u_{q}^{i}({\textbf {x}})}$, with ${\displaystyle q=1,\cdots ,Q}$ and ${\displaystyle i=1,\cdots ,R_{q}}$ have zero mean and covariance cov${\displaystyle [u_{q}^{i}({\textbf {x}}),u_{q'}^{i'}({\textbf {x}})']=k_{q}({\textbf {x}},{\textbf {x}}')}$ if ${\displaystyle i=i'}$ and ${\displaystyle q=q'}$. But ${\displaystyle \operatorname {cov} [f_{d}({\textbf {x}}),f_{d'}({\textbf {x}}')]}$ is given by ${\displaystyle ({\textbf {K}}({\textbf {x}},{\textbf {x}}'))_{d,d'}}$. Thus the kernel ${\displaystyle {\textbf {K}}({\textbf {x}},{\textbf {x}}')}$ can now be expressed as

${\displaystyle {\textbf {K}}({\textbf {x}},{\textbf {x}}')=\sum _{q=1}^{Q}{{\textbf {B}}_{q}k_{q}({\textbf {x}},{\textbf {x}}')}}$

where each ${\displaystyle {\textbf {B}}_{q}\in {\mathcal {\textbf {R}}}^{D\times D}}$ is known as a coregionalization matrix. Therefore, the kernel derived from LMC is a sum of the products of two covariance functions, one that models the dependence between the outputs, independently of the input vector ${\displaystyle {\textbf {x}}}$ (the coregionalization matrix ${\displaystyle {\textbf {B}}_{q}}$), and one that models the input dependence, independently of ${\displaystyle \left\{f_{d}({\textbf {x}})\right\}_{d=1}^{D}}$(the covariance function ${\displaystyle k_{q}({\textbf {x}},{\textbf {x}}')}$).

The ICM is a simplified version of the LMC, with ${\displaystyle Q=1}$. ICM assumes that the elements ${\displaystyle b_{d,d'}^{q}}$ of the coregionalization matrix ${\displaystyle \mathbf {B} _{q}}$ can be written as ${\displaystyle b_{d,d'}^{q}=v_{d,d'}b_{q}}$, for some suitable coefficients ${\displaystyle v_{d,d'}}$. With this form for ${\displaystyle b_{d,d'}^{q}}$:

${\displaystyle \operatorname {cov} \left[f_{d}(\mathbf {x} ),f_{d'}(\mathbf {x} ')\right]=\sum _{q=1}^{Q}{v_{d,d'}b_{q}k_{q}(\mathbf {x} ,\mathbf {x} ')}=v_{d,d'}\sum _{q=1}^{Q}{b_{q}k_{q}(\mathbf {x} ,\mathbf {x} ')}=v_{d,d'}k(\mathbf {x} ,\mathbf {x} ')}$

where

${\displaystyle k(\mathbf {x} ,\mathbf {x} ')=\sum _{q=1}^{Q}{b_{q}k_{q}(\mathbf {x} ,\mathbf {x} ')}.}$

In this case, the coefficients

${\displaystyle v_{d,d'}=\sum _{i=1}^{R_{1}}{a_{d,1}^{i}a_{d',1}^{i}}=b_{d,d'}^{1}}$

and the kernel matrix for multiple outputs becomes ${\displaystyle \mathbf {K} (\mathbf {x} ,\mathbf {x} ')=k(\mathbf {x} ,\mathbf {x} ')\mathbf {B} }$. ICM is much more restrictive than the LMC since it assumes that each basic covariance ${\displaystyle k_{q}(\mathbf {x} ,\mathbf {x} ')}$ contributes equally to the construction of the autocovariances and cross covariances for the outputs. However, the computations required for the inference are greatly simplified.

Another simplified version of the LMC is the semiparametric latent factor model (SLFM), which corresponds to setting ${\displaystyle R_{q}=1}$ (instead of ${\displaystyle Q=1}$ as in ICM). Thus each latent function ${\displaystyle u_{q}}$ has its own covariance.

While simple, the structure of separable kernels can be too limiting for some problems.

Notable examples of non-separable kernels in the regularization literature include:

• Matrix-valued exponentiated quadratic (EQ) kernels designed to estimate divergence-free or curl-free vector fields (or a convex combination of the two)[8][18]
• Kernels defined by transformations[8][19]

In the Bayesian perspective, LMC produces a separable kernel because the output functions evaluated at a point ${\displaystyle {\textbf {x}}}$ only depend on the values of the latent functions at ${\displaystyle {\textbf {x}}}$. A non-trivial way to mix the latent functions is by convolving a base process with a smoothing kernel. If the base process is a Gaussian process, the convolved process is Gaussian as well. We can therefore exploit convolutions to construct covariance functions.[20] This method of producing non-separable kernels is known as process convolution. Process convolutions were introduced for multiple outputs in the machine learning community as "dependent Gaussian processes".[21]

Approached from the regularization perspective, parameter tuning is similar to the scalar-valued case and can generally be accomplished with cross validation. Solving the required linear system is typically expensive in memory and time. If the kernel is separable, a coordinate transform can convert ${\displaystyle \mathbf {K} (\mathbf {X} ,\mathbf {X} )}$ to a block-diagonal matrix, greatly reducing the computational burden by solving D independent subproblems (plus the eigendecomposition of ${\displaystyle \mathbf {B} }$). In particular, for a least squares loss function (Tikhonov regularization), there exists a closed form solution for ${\displaystyle {\bar {\mathbf {c} }}}$:[8][14]

${\displaystyle {\bar {\mathbf {c} }}^{d}=\left(k(\mathbf {X} ,\mathbf {X} )+{\frac {\lambda _{N}}{\sigma _{d}}}\mathbf {I} \right)^{-1}{\frac {{\bar {\mathbf {y} }}^{d}}{\sigma _{d}}}}$

There are many works related to parameter estimation for Gaussian processes. Some methods such as maximization of the marginal likelihood (also known as evidence approximation, type II maximum likelihood, empirical Bayes), and least squares give point estimates of the parameter vector ${\displaystyle \phi }$. There are also works employing a full Bayesian inference by assigning priors to ${\displaystyle \phi }$ and computing the posterior distribution through a sampling procedure. For non-Gaussian likelihoods, there is no closed form solution for the posterior distribution or for the marginal likelihood. However, the marginal likelihood can be approximated under a Laplace, variational Bayes or expectation propagation (EP) approximation frameworks for multiple output classification and used to find estimates for the hyperparameters.

The main computational problem in the Bayesian viewpoint is the same as the one appearing in regularization theory of inverting the matrix

${\displaystyle {\overline {\mathbf {K} (\mathbf {X} ,\mathbf {X} )}}=\mathbf {K} (\mathbf {X} ,\mathbf {X} )+{\boldsymbol {\Sigma }}.}$

This step is necessary for computing the marginal likelihood and the predictive distribution. For most proposed approximation methods to reduce computation, the computational efficiency gained is independent of the particular method employed (e.g. LMC, process convolution) used to compute the multi-output covariance matrix. A summary of different methods for reducing computational complexity in multi-output Gaussian processes is presented in.[8]