Large scale multi-output multi-class classification using Gaussian processes

Multi-output Gaussian processes (MOGPs) can help to improve predictive performance for some output variables, by leveraging the correlation with other output variables. In this paper, our main motivation is to use multiple-output Gaussian processes to exploit correlations between outputs where each output is a multi-class classification problem. MOGPs have been mostly used for multi-output regression. There are some existing works that use MOGPs for other types of outputs, e.g., multi-output binary classification. However, MOGPs for multi-class classification has been less studied. The reason is twofold: 1) when using a softmax function, it is not clear how to scale it beyond the case of a few outputs; 2) most common type of data in multi-class classification problems consists of image data, and MOGPs are not specifically designed to image data. We thus propose a new MOGPs model called Multi-output Gaussian Processes with Augment & Reduce (MOGPs-AR) that can deal with large scale classification and downsized image input data. Large scale classification is achieved by subsampling both training data sets and classes in each output whereas downsized image input data is handled by incorporating a convolutional kernel into the new model. We show empirically that our proposed model outperforms single-output Gaussian processes in terms of different performance metrics and multi-output Gaussian processes in terms of scalability, both in synthetic and in real classification problems. We include an example with the Ommiglot dataset where we showcase the properties of our model.

Our main purpose in this paper is to use MOGPs to study the problem of multiple outputs where each output is a multi-class classification problem. The setting considered here goes beyond multi-label classification since we allow each output to potentially have its own inputs moving into the multi-task setting.
MOGPs have mainly been used for multi-output regression to predict continuous variables (Bonilla et al., 2008;Álvarez et al., 2012;Dai et al., 2017). In this setting, the assumption is that each output follows a Gaussian likelihood and the mean of the Gaussian likelihood is given by one output of the MOGP. Due to the properties of the Gaussian distribution, Bayesian inference is tractable in this case.
Beyond the muti-output regression problem, there is some research on other types of outputs in MOGPs. For example, Skolidis and Sanguinetti (2011) use MOGPs to model a setting where each output corresponds to a binary classification problem. Each binary outcome is modelled using a probit likelihood. The MOGP corresponds to the so called intrinsic coregionalisation model (ICM) (Bonilla et al., 2008). Since Bayesian inference is intractable in this model, the authors approximate posterior distributions using expectation-propagation and variational Bayes.
Several research works have addressed the case of multi-class classification using GPs. Previous works have used the softmax likelihood (Williams & Rasmussen, 2006;Kim & Ghahramani, 2006;Galy-Fajou et al., 2020), the multinomial probit likelihood function (Girolami & Rogers, 2006), the step function (Hernández-Lobato et al., 2011). Recently, Liu et al. (2019) can use all the above likelihoods through additive noise terms. The parameters in these likelihood functions are assumed to follow independent Gaussian processes. Another strand of works generalise this setting by allowing correlated Gaussian processes for the latent parameters of the likelihood functions, typically using MOGPs. Both Dezfouli and Bonilla (2015) and Chai (2012) use an ICM for a single-output multi-class classification problem modelled through a multinomial logistic likelihood, i.e. the softmax likelihood. In terms of Bayesian inference, Chai (2012) proposes a variational sparse approximation for the posterior distribution, and based on scalable automated variational inference, Dezfouli and Bonilla (2015) approximates the posterior distribution by a mixture of Gaussians. Moreno-Muñoz et al. (2018) build a heterogeneous multi-output Gaussian process, where each output has its own likelihood, through a linear model of coregionalisation (LMC) (Álvarez et al., 2012). Moreno-Muñoz et al. (2018) use an stochastic variational approach for Bayesian inference.
The approaches for single-output multi-class classification described above are restricted to the case where the number of classes is small. They scale poorly when the number of classes go beyond a few tens. Scalability is also poorly handled by the more general model of Moreno-Muñoz et al. (2018) for the multi-output multi-class classification case, where once again, problems that go beyond a few tens of classes are out of reach.
Our main contribution in this paper is that we introduce a new extension of multi-output GPs able to handle large scale multi-output multi-class classification problems, typically in the range of hundreds and even thousands of classes. We achieve scalability by subsampling both training input data and classes in each output, by using stochastic variational inference (Hensman et al., 2013;Moreno-Muñoz et al., 2018), and by choosing a softmax likelihood function via Gumbel noise error for all outputs. We refer to this model as Multioutput Gaussian Processes with Augment & Reduce (MOGPs-AR).
We also enable our MOGPs-AR to allow downsized images as input data. To efficiently deal with downsized images, we employ convolutional kernels (Van der Wilk et al., 2017), computing the entries of the kernel matrices using kernels over patches of the images and integrating these kernels within a MOGP. Since our model is able to capture both intraand inter-output dependencies, it also provides a means to perform transfer learning in the multi-task setting. We show an example of this multi-task learning ability of our model in the Ommiglot dataset. To the best of our knowledge, this is the first time that a multi-task multi-class Gaussian process model is used over such dataset.

Related work
As we mentioned early, the multi-class classification problem has been mainly studied using single-output GPs (Williams & Rasmussen, 2006;Kim & Ghahramani, 2006;Hernández-Lobato et al., 2011;Girolami & Rogers, 2006;Liu et al., 2019). The model introduced in this paper, MOGPs-AR, uses the softmax likelihood through additive noise errors, which is the same as Liu et al. (2019). However, MOGPs-AR solves multiple outputs problems together while the model in Liu et al. (2019), like all single-output GPs, only solves single output problems. Regarding single output problems, MOGPs-AR can also improve prediction using a correlation between all latent parameter functions whereas single-output GPs cannot capture the correlation.
The works more relevant to ours are Chai (2012); Dezfouli and Bonilla (2015); Skolidis and Sanguinetti (2011);Moreno-Muñoz et al. (2018). Both Chai (2012) and Dezfouli and Bonilla (2015) can only handle a single output multi-class classification problem even if they use MOGPs. Nevertheless, our model can tackle multiple outputs where each output is a multi-class classification problem. Skolidis and Sanguinetti (2011) only solve multioutput binary classification problems, which is different to ours. Compared with Skolidis and Sanguinetti (2011), our inference method is also suited to large scale data sets. Moreno-Muñoz et al. (2018) can tackle multi-output multi-class classification problems and develop a similar stochastic variational inference method as us. However, we are different to Moreno-Muñoz et al. (2018) since we can cope with a large number of classes by subsampling classes and also can deal with downsized images through convolutional kernels ( Van der Wilk et al., 2017).
The work by Panos et al. (2021) is much related to us since we use a similar subsampling method. Panos et al. (2021) extend a semiparametric latent model, a special case of LMC, to address the multi-label problem by using sigmoidal/Bernoulli likelihood for each latent parameter function. Panos et al. (2021) can doubly subsample data points and classes to reduce computational complexity based on stochastic variational inference, which is analogous to us. However, we are different in other aspects. First, we solve multi-class classification problems using the softmax likelihood instead of multi-label problems using sigmoidal/Bernoulli likelihood. Further, we can apply a convolutional kernel to handle downsized image data. Finally, our model can deal with multi-output problems instead of only tackling single output problems.

3 3 Methodology
In this section, we will derive the MOGPs-AR model. We first develop the LMC model with a convolutional kernel. We then define the softmax likelihood through augmenting noise data. We finally describe stochastic variational inference and the approximated predictive distribution for our model. We assume there are D different outputs ( Table 1 shows the description of our notation). The vector y ∈ R D groups all the D different outputs: , ..., D}) is a categorical variable and C d is the number of classes in the d-th output. Like Moreno-Muñoz et al. (2018), we also assume that those outputs are conditionally independent given parameters is defined by latent parameter functions: is c-th latent parameter function in the d-th output evaluated at x . We then obtain: (2) Gaussian process covariance function of u q (x) is a group of latent parameter functions defining the parameters in d (x).

Combining with convolutional kernel
We use the linear model of coregionalisation (LMC) and combine it with the convolutional kernel. The LMC is a popular model in MOGPs, where each output is expressed as a linear combination of a collection of Gaussian processes (Álvarez et al., 2012). The convolutional kernel (Van der Wilk et al., 2017) can effectively exploit features in images. We construct a convolutional structure for mutually independent latent functions U = u q (x) Q q=1 where u q follows a Gaussian process, Q is the number of the latent functions and each latent parameter function f c d (x) is a linear combination of the latent functions U . Here, we assume x ∈ R W×H is an image data point that has a v = W × H size where W and H are the width and height of the image separately. We also assume x [p] is the p th patch of x with patches of size E = w × h where w and h are the width and height of each patch, respectively.
After dividing an image into patches, we get a total of P = (W − w + 1) × (H − h + 1) patches. We begin with a patch response function u q x [p] ∶ R w×h → R , which maps a patch of size E = w × h to a real number in R . Then we add a weight for each patch response function and get a latent function u q (x) ∶ R W×H → R , where u q (x) is the sum of all patch responses with weights: u q (x) = ∑ p w p u q � x [p] � . Each function u q is drawn from an independent GP prior: u q (⋅) ∼ GP 0, k q (⋅, ⋅) , where k q (⋅, ⋅) can be any kernel function. In this paper, we use the radial basis function kernel with automatic relevance determination (RBF-ARD) (Williams & Rasmussen, 2006): ard is a variance parameter and l j is the length scale for the j-th input dimension. Therefore k q (⋅, ⋅) is RBF-ARD. When all length scales are the same, the kernel is called radial basis function kernel (RBF) (Lawrence & Hyvärinen, 2005). Hence, each f c d (x) is defined as where a i d,c,q ∈ R can be considered as a weight on U and we assume { q } Q q=1 are the hyperparameters for {k q (⋅, ⋅)} Q q=1 , with q being the hyperparameters for the kernel k q (⋅, ⋅) . R q represents the number of latent functions u i q (x) that are sampled independently and identically from the Gaussian processes u q (⋅) ∼ GP 0, k q (⋅, ⋅) . The difference between the convolutional kernel model and a more classic kernel, e.g., RBF, is that we use the convolutional structure term ∑ P p=1 w p u i q (x [p] ) instead of solely u i q (x) , where Fig. 1 shows an example of two images and how they are handled through the convolutional kernel. With q = 1, ..., Q

3
and i = 1, ..., R q , the function u i q (x) have a zero mean and covariance cov u i Let the mean function of f c d (x) be zero and the cross-covariance function of f c d (x) be Because u i q (⋅) is independently and identically drawn from u q (⋅) and U(⋅) are mutually independent For simplicity in the presentation, we assume that all outputs y d (x) have a collection of the same input vectors X = x n N n=1 ∈ R N×v . Our model also works for each output with a different input data set. For notation simplicity, we define . We consider two characters as two classes. The two images are one data point for each class separately. Left: The whole image is considered as an input data point x and the blue grid represents the p-th patch x [p] . Right: The whole image is considered as an input data point x ′ and the blue grid represents the p ′ -th patch x �[p � ] (Color figure online)

3
The prior distribution of f is given by known as a coregionalisation matrix and it controls the correlation between each latent parameter function.

Augmenting model by noise data
In this section, we generalise the model in the last subsection to cope with the multi-output multi-class classification problem using the softmax likelihood. We derive a softmax likelihood function through Gumbel noise error for a generic output y d .
We take the d-th output y d (x) with the latent parameter function We first add a Gumbel noise error to each of latent parameter functions f for each of the classes in the d-th output. We thus obtain: We then employ the internal step likelihood (Liu et al., 2019): where G and Φ G are the probability density function and the cumulative distribution function of the Gumbel distribution, respectively. (NB: We drop out the c in c d,i for convenience since all the c d,i are from the same Gumbel error distribution). Now, we assume the Gumbel error d, (20), we recover a softmax likelihood (Liu et al., 2019;Ruiz et al., 2018): The softmax likelihood is a common likelihood used in multi-class classification with Gaussian processes (Williams & Rasmussen, 2006). As we mentioned in expression (5), all outputs are conditional independent given the corresponding latent parameter functions so each output has its own likelihood expression (20).

Scalable variational inference
We have derived the LMC model with a convolutional kernel and used the softmax likelihood. However, there exists a computational challenge if there are a very large number of classes and training instances in each output. We thus use scalable variational inference to reduce the computational complexity by the techniques of inducing patches and subsampling, where we refer our model to Multi-output Gaussian Processes with Augment & Reduce (MOGPs-AR). Inducing patches can ease the computational complexity of the inversion of a kernel matrix from O N 3 to O NM 2 , where N is the number of data points per output and M is the number of inducing patches ( M ≪ N ). Subsampling reduces the computational complexity of our model using a subset of both training data and classes for each output during hyperparameters and parameters optimisation.

Inducing patches for MOGPs-AR
We assume we use image data sets in this section. We define the inducing patches ( Van der Wilk et al., 2017) at the latent functions U . If our input data sets are not image data sets, we use inducing points (Hensman et al., 2013). The difference between the inducing points and the inducing patches is the dimension size. The dimensions of the inducing points are the same as the input data, whereas the dimensions of the inducing patches match the patch of the images.
We first define a group of M inducing patches (Van der Wilk et al., 2017) We then obtain The latent parameter functions f c d are conditionally independent given u . We therefore obtain the conditional distribution p(f|u): where K uu ∈ R QM×QM is a block-diagonal matrix based on K q as each block.
Based on Moreno-Muñoz et al. (2018); Liu et al. (2019), we obtain the lower bound L for log p(y): where where q u q = N u q | u q , S u q and we refer u q , S u q Q q=1 as the variational hyperparameters that need to be optimised. Further, we get (see the Appendix A for detail): and S u is a block diagonal matrix with blocks given as S u q .
After calculation (NB: detail is in the Appendix A), we obtain: where

Reducing computational complexity by subsampling
To reduce the computational complexity of our model, we use only a subset of the data observations and a subset of the classes to estimate the parameters and hyperparameters. The optimal parameters and hyperparameters are chosen by maximising an unbiased estimator of L (37), where the unbiased estimator is obtained through a subset of both training data points and classes in each output. In our model, the hyperparameters are Z , a c,q , {w p } P p=1 and the variational parameters are u q , S u q Q q=1 for {q u q } q=Q q=1 . We refer all those parameters as Θ . To obtain the optimal values of the parameters Θ , we use gradient descent to maximise L with respect to Θ: where, for notation simplicity, we define We then estimate ∇ Θ L by randomly sampling a subset of data points ∇ ΘL is an unbiased estimator for ∇ Θ L where the computational complexity of MOGPs-AR is dominated by optimising the parameters through maximising L.
Our sampling strategy is in Algorithm 1. The computational complexity of MOGPs-AR mainly depends on the inversion of K u u with complexity O(QM 3 ) and products like Kf u with complexity O D|S||B|QM 2 ; If we do not use the subsampling of classes, we have to calculate products like K fu with a cost of O C|B|QM 2 , where the notations are defined as below: MOGPs-AR alleviates the computational complexity of the product Kf u from O C|B|QM 2 to O D|S||B|QM 2 .

Prediction
In this subsection, we derive the predictive distribution of MOGPs-AR. Considering a new test input x * in the d-th output, we assume p(u|y) ≈ q(u) and approximate the pre- parameter functions q f d x * are mutually independent, so we obtain We can use Monte Carlo to approximate the integral in the same way as Liu et al. (2019).

Experiments
In this section, we evaluate MOGPs-AR in various data sets. We apply MOGPs-AR in a synthetic data set to show its scalability in the number of classes compared to multi-output Gaussian processes. We also compare MOGPs-AR to other models in different real data sets. Further, to test MOGPs-AR the capacity in dealing with an image data set, we compare the performance of a convolutional kernel and RBF-ARD in our model. Baselines. We compare the MOGPs-AR with the following two single-output and one multi-output Gaussian process models: 1) A Gaussian process for multi-class classification model (G-M), an independent Gaussian process using the softmax likelihood. 2) A Gaussian process multi-class classification with additive noise model (G-A), an independent Gaussian process using the softmax likelihood via Gumbel noise. 3) A multi-output Gaussian process model for multi-class classification problems (MG-M), a standard linear model of coregionalisation for MOGPs using the softmax likelihood. For all the different models in this paper, we use RBF-ARD, unless otherwise stated. For all models, we use traditional inducing points (Hensman et al., 2013) unless mentioned otherwise. All models are trained using the Adam optimiser with 0.01 learning rate and trained by 4000 iterations (Kingma & Ba, 2014). We use the same 80% training and 20% validation data set to choose the optimal number Q of latent functions U , where we optimise again all hyperparameters during cross-validation, for MOGPs-AR and MG-M.
Evaluation Metrics. There are three different evaluation metrics in this paper: where true and prediction are sets of true and predicted pairs (input data point, class) separately (e.g., true,n = x n , y true,n where x n is the n-th input data point and y true,n is the F 1 ( l prediction , l true ) = 2 P( l prediction , l true ) × R( l prediction , l true ) P( l prediction , l true ) + R( l prediction , l true ) , corresponding class for x n . The l true and l prediction are subsets of true and prediction separately (e.g., l true = { x n , y true,n ∈ true | y true,n = l, n ∈ ℕ} ). The and ℕ are the sets of classes and input data points, respectively. The formulas use P( l prediction , l true ) = 0 or R( l prediction , l true ) = 0 if l prediction = � or l true = �. The synthetic data experiment was performed on a Dell PowerEdge C6320 with an Intel Xeon E5-2630 v3 at 2.40 GHz and 64GB of RAM. All real data experiments were performed on a PowerEdge R740XD Server with NVIDIA Tesla v100 32GB GDDR. 1

Synthetic data
In this subsection, we compare the performance of MOGPs-AR with MG-M on synthetic data where we generate a single output classification synthetic data set. 2 We create a 20 class data set by assigning a cluster of 100 points normally distributed, where each data point has five features, to each class. In total, there are 2000 samples. Since the synthetic data has 20 classes, we refer to it as S-20. We use 20 classes to compare MOGPs-AR with MG-M in terms of scalability. MOGPs-AR and MG-M use the same parameter setting (see Table 3) exclude that MOGPs-AR used a different number of subset classes.
We compare MOGPs-AR with MG-M in terms of training time and Recall-Weighted performance. Figure 2 shows the mean training time for MOGPs-AR is less than MG-M in five folds cross-validation. This is because the computational complexity of MOGPs-AR is less than MG-M. As we mentioned in 3.3.2, compared to MG-M, MOGPs-AR reduce the computational complexity of the product Kf u from O C|B|QM 2 to O D|S||B|QM 2 where D|S| ≪ C . Figure 2 empirically shows the mean training time of MOGPs-AR (1) with 596s is nearly onesixth of MG-M with 3641s. The mean training time in MOGPs-AR increases as the number of |S d | increases but it is still less than MG-M. While MOGPs-AR has less training time than MG-M, it has a similar performance in Recall-Weighted with MG-M for S-20. Even if we use a small subset of classes, e.g., five classes, MOGPs-AR also has a close performance to MG-M (see Fig. 2 right panel). The Recall-Weighted of MOGPs-AR slightly increases as the number of samples increase. Further, we notice that MOGPs-AR (17) has a better performance than MG-M. In theory, MOGPs-AR should have the same performance as MG-M. However, we can not perform convex optimisation for both MG-M and MOGPs-AR, so MOGPs-AR may outperform MG-M in various performance metrics in practice.

Single-output GP classification: four real data sets
We will use the following four real data sets to test the performance of the different GP classifiers: 1) Balance (Dua & Graff, 2017) is a data set for the results of psychology experiments. There are 625 data points with four discrete variables: Left-Weight, Left-Distance, Right-Weight and Right-Distance. The value of all four discrete variables ranges from one to five. The data set consists of three classes: the balance scale tipped to the right (R), tipped to the left (L) or be balanced (B). 2) CANE9 (Dua & Graff, 2017) contains 1080 documents of free text business descriptions of Brazilian companies. Those documents are divided into nine different categories. Each document has 856 integer variables (word frequency). 3) Mediamill (Snoek et al., 2006) is a multi-label data set for generic video indexing. To apply multi-classification, we only maintained one label, which is the first label to appear, for each data point. Further, we only use part of this data set since the original data set is highly imbalanced. We then obtain the number of data points for each class ranged from 31 to 545. In total, we have 6689 data points with 120 numeric features and 35 classes. 4) Bibtex data set (Katakis et al., 2008) is also a multi-label data that contains 7395 Bibtex entries with 1836 variables. Similarly, we only maintained one label, which is the first label to appear, obtaining 148 classes.
In all three performance measures, MOGPs-AR outperforms G-A and G-M on all four data sets (see Figure. 3). This is because MOGPs-AR can use each of latent parameter functions f , which is a linear combination of latent functions U , to predict each class. The underlying function of the latent functions U and B q can transfer knowledge between each class in the same output. However, G-A and G-M only have independent Gaussian processes that can not capture the similarity between each class. Further, Fig. 3 indicates that using a small subset of classes (e.g., MOGPs-AR(1) or MOGPs-AR(5)), MOGPs-AR obtains a similar result as MG-M for Balance, CANE9, Mediamill data sets as we have discussed as in 4.1.
Compared with single output Gaussian processes, MOGPs-AR can achieve around 10% improvement in terms of three performance metrics on Balance and CANE9 data set ( Fig. 3 upper panel). The optimal number (Q) of latent functions U is two and nine for the Balance and CANE9 data sets separately. Those latent functions share the knowledge between each class and help to improve the performance. There is also a connection between single output and multi-output Gaussian processes. Considering an extreme case, we assume there is only one class, Q=1 and B q = 1, MOGPs-AR and MG-M have the same structure as G-A and G-M in theory, respectively.
Regarding both Mediallmill (35 classes) and Bibtex (148 classes), MOGPs-AR has excellent performance compared to the single-output Gaussian processes and MG-M. For the Mediamill dataset, based on capturing dependency between each class, MOGPs-AR is nearly 6 times better than G-A and 4 times better than G-M in terms of F1-Weighted, where the mean of F1-Weighted is 0.04 for G-A, 0.08 for G-M and 0.25 for MOGPs-AR. Further, we cannot apply MG-M in the Bibtex data set since it is not able to compute K fu (out of memory). However, MOGPs-AR scales well since it only uses a subset of classes (MOGPs-AR (20)) for prediction.

Multi-output GPs classifications: UJIIndoorLoc
To compare the performance of MOGPs-AR in multi-output multi-class classification problems, we apply MOGPs-AR to UJIIndoorLoc Data Set (Torres-Sospedra et al., 2014). There are 21048 instances that rely on WIFI fingerprint for three buildings of Universitat Jaume I where Building I and Building II have four floors respectively and Building III has five floors. Each instance has 520 features based on signal strength intensity. We randomly sample 200 data points from each floor so there are 800 data points for Building I and Building II respectively and 1000 data points for Building III. Further, we standardise the data set for each Building. We make predictions for each floor depending on the 520 features. Since there are three buildings of Universitat Jaume I, we assume there is a strong correlation between each building. We regard each building as each output and different floors as different classes in our model. The UJIIndoorLoc is considered as a multi-output multi-class classification problem. In this experiment and following, we do not apply the MG-M model due to its computational complexity. MOGPs-AR can be an alternative model for MG-M so we only consider MOGPs-AR and two single output GP models. Figure 4 shows that MOGPs-AR outperforms single-output Gaussian processes in both Building I, II and III in all three performance measures. For example, MOGPs-AR can achieve around 50% improvement in terms of Recall-Weighted on Building I compared with single output Gaussian processes. The reason is that MOGPs-AR can capture dependencies of intraand inter-three buildings. The dependencies can help improve the prediction for all buildings.

3
The single-output Gaussian processes cannot use the dependency so the single output Gaussian process does not perform well in UJIIndoorLoc.
To investigate the correlation between intra-and inter output, we create a global absolute coregionalisation matrix. First, we create absolute coregionalisation matrices B abs q Q q=1 , where B abs q ∈ R C×C , by taking the absolute value of each entry in B q . Second, we obtain the mean of those absolute coregionalisation matrices: B = 1 Q ∑ Q q=1 B abs q and B ∈ R C×C . Since we are performing K-fold cross-validation, we have K different mean absolute coregionalisation matrices: , where B i ∈ R C×C refers to the mean absolute coregionalisation matrices during the i-th fold cross-validation. Further, we calculate the mean of B where B i,j ∈ R C i ×C j indicates the correlations for all latent parameter functions between i-th output and j-th output. At the end, in order to find the correlation for outputs independently, we calculate a scalar B i,j = 1 , which represents dependence between i-th output and j output. We therefore define a global absolute coregionalisation matrix (GACM ∈ R D×D ) as the following: Fig. 4 Performance in cross-validation (mean ± standard deviation) Figure 5a shows the correlation between each building captured by our model. We can notice there is a strong correlation between the different buildings. Building I and Building II have a relatively strong correlation compared to Building I and Building III, Building II and Building III. Building II has the strongest intra-output correlation while Building III has the smallest intra-output correlation among those three buildings.

Multi-output GPs classifications: Omniglot-dataset
We apply MOGPs-AR to Omniglot image data set (Lake et al., 2015). The Omniglot data set includes 1623 various handwritten characters from 50 distinct alphabets. Each of the 1623 characters was drawn by 20 different people (the total number of images is 32460). Although traditional MOGPs are not specifically designed to deal with image data, MOGPs-AR can handle image data by incorporating a convolutional kernel ( Van der Wilk et al., 2017). The size of each image is 105 × 105 pixels. To help speeding up the computation and reducing the computational complexity in the covolutional kernel, we resize the images from 105 × 105 to 20 × 20 as Santoro et al. (2016) did. We regard each alphabet as an output in our model. Each alphabet has different characters which are considered as different classes. Therefore, we consider the Omniglot data set as multi-output multi-class classification problems.

Ojibwe and blackfoot alphabets
To compare the performance of MOGPs-AR in multi-output multi-class classification problems and image input data, we first consider Ojibwe and Blackfoot alphabets as two different multi-class classification problems (see Fig. 6). Since the two alphabets are from Canadian Aboriginal Syllabics, we assume there is a strong correlation between them. Our model can capture the correlation through joint modelling of the two alphabets to improve predictive performance for each multi-class classification problem. There are 14 different characters in each output so there are 14 classes, and each class has 20 data points. We compare both the RBF-ARD kernel and the convolutional kernel. Table 4 shows the parameter setting in Omniglot data set.
In Fig. 7 we show that MOGPs-AR outperforms single-output Gaussian processes in both alphabets in terms of the convolutional kernel or RBF-ARD. The reason is that MOGPs-AR can capture the dependency between the two alphabets. The dependency can help improve the prediction for both alphabets. The single output Gaussian processes cannot use the dependency so the single output Gaussian process with either the convolutional kernel or RBF-ARD does not perform well in both Ojibwe and Blackfoot. The size of the mini-batch is too small that has also a negative influence on the single output Gaussian processes (Fig. 8). Especially, the values of the three performance metrics are closed to 0.05 for G-A with the convolutional kernel on Ojibwe.
Since both alphabets are from Canadian Aboriginal Syllabic we expect they have a strong correlation. Figure 5b indeed shows there is a similar global correlation between intra-and inter-output for both alphabets, which indicates that our model has the capacity of capturing the underline correlation among those correlated data sets. MOGPs-AR with the convolutional kernel outperforms MOGPs-AR with RBF-ARD in both alphabets in terms of three performance metrics (see Fig. 7). For example, MOGPs-AR improves the Recall-Weighted from 0.468 to 0.714 by changing RBF-ARD to the convolutional kernel on the Blackfoot alphabet. Moreover, we also combine G-M and G-A with the convolutional kernel and they also have stronger performance compared with RBF-ARD. In particular, G-M with the convolutional kernel obtains 0.5858 compared with 0.0857 using RBF-ARD in terms of Recall-Weighted on the Blackfoot alphabet. The performance of G-M with the convolutional kernel (0.5858) is better than MOGPs-AR with RBF-ARD (0.468) on the Blackfoot alphabet. The reason is that the convolutional kernel is more effectively capturing image-level features than the RBF-ARD kernel.
To investigate the effects of mini-batch size, we set up another experiment. We train again the exact same models with the parameters initialised in the same way as the experiment above but using different mini-batch sizes (e.g., 50, 70, 90). Since the convolution kernel provided better results in the previous experiments, we only show the results using the convolution kernel and the Recall-Weighted performance measure for both alphabets. Figure 8 shows that the size of the mini-batch has more influence on single output Gaussian processes than MOGPs-AR. A small size number for the mini-batch, e.g., 50, has a negative impact on G-M and G-A. However, MOGPs-AR has a slight increase in performance or keeps a similar result with the mini-batch size increasing. G-A and G-M improve the performance as the mini-batch size grows up from 50 to 90. When the size of the minibatch is 90, G-M has a similar performance with MOGPs-AR. However, when we consider the mini-batch of size 50, MOGPs-AR still can get good performance compared to single GPs.

All alphabets
In our final experiment, we apply MOGPs-AR in 50 alphabets in the original data set. There are 50 outputs with different number of classes in each output (for more details on the number of classes in each output see Table 2). The total number of classes in the 50 outputs are 1623. We follow Lake et al. (2015) and split the 50 alphabets into two sets: a background set and an evaluation set, where the background set has 30 alphabets (with a total of 964 classes) and the evaluation set has 20 alphabets (with a total of 659 classes). In order to apply MOGPs-AR in all 50 alphabets, we use a mini-batch size of nine data points for each output to train our model. The small mini-batch size has a negative impact on G-M and G-A so we only apply MOGPs-AR in this experiments. We apply MOGPs-AR for three different sets of alphabets: all alphabets, the Background alphabets and the Evaluation alphabets.
In Fig. 9, we empirically (50 different outputs and a total of 1623 classes of image data) show that MOGPs-AR has better scalability than traditional multi-output Gaussian processes. MOGPs-AR reduces the computational complexity by subsampling both training data points and classes in each output. Figure 9 also indicates that MOGPs-AR obtains good performance even if we choose a small size of mini-batch (nine) and only a small number of classes (one) in each output since it captures both intra-and inter-output correlation.
In most predictions, our model trained with the data of all alphabets could outperform one trained with the data of part of the alphabets. For example, our model trained using all alphabets improves the Recall-Weighted from 0.6096 to 0.6692 for the Aurek alphabet, compared with one using evaluation alphabets for training. The extra alphabets can help our model improve its performance.
However, there are exceptions to the scenario in the last paragraph. For example, for the Syriac (Estrangelo) alphabet, the values of the Recall-Weighted 0.5174 is less than 0.5283 where only use background alphabets for training our model. One likely reason is that our model assumes a correlation with all alphabets. However, the correlation with those alphabets may not exist or the correlation may hinder the predictive performance.

Conclusion
In this paper, we have introduced MOGPs-AR, a novel framework that allows the use of multioutput Gaussian processes for multi-output multi-class classification. MOGPs-AR can tackle large scale data sets and a large number of classes in each output. Further, when combined with the convolutional kernel, it is suited for downsized image data.
We experimentally show that MOGPs-AR has a similar result to MG-M that is a linear model of coregionalization and uses a similar stochastic variational inference method as us. However, the training time of MOGPs-AR is less than MG-M. Experimental results in various data sets also indicate that MOGPs-AR significantly improves the performance compared to single output Gaussian processes.
MOGPs-AR has good performance in extreme classification using a softmax function which is only suited to each instance associated with a single class. Because of the softmax function, MOGPs-AR can not deal with a multi-label problem where each data point belongs to multiple classes. It would be an interesting work for future research if we can generalise MOGPs-AR for the multi-label problem that has a strong correlation with extreme classification problems.
A practical application of Gaussian process models to realistic image recognition tasks is still an open research problem. For example, in terms of accuracy performance in a realistic RGB set CIFAR-10 (Krizhevsky & Hinton, 2009), the accuracy performance of Gaussian processes ( Van der Wilk et al., 2017;Blomqvist et al., 2019) is not as high as the state-of-the-art like deep learning. A potential extension would be to consider integrating the structural properties of deep learning architectures into our model by using deep kernel learning (Wilson et al., 2016).

Appendix A Complete derivation of the lower bound L
To compute the derivation of the lower bound L , we begin with the following: where q(f, u, ) = p(f|u)q(u)q( |f) . We assume q(f) ≈ p(f | y) so we obtain The above function means the latent parameter functions are mutually independent in q(f) . Then, we obtain: To get a tight bound, we derivative L with respect to d,i , We thus obtain the optimal value * d,i = P d,i + 1. After substitution of d,i by * d,i , there is Table 2 shows the number of data points and classes for each alphabet in the Omniglot data set. The columns of Background set and Evaluation set have shown 30 and 20 alphabets separately.

Appendix B Omniglot data
Mauricio Álvarez supervised the development of the research and provided feedback at all the stages of the process including editing the final manuscript. Chunchao Ma and Mauricio Álvarez contributed equally to the conception of the main research ideas developed in the manuscript.

Conflict of interest
The authors declare that they have no conflict of interest.

Ethics approval and Consent to participate
The authors declare that this research did not require Ethics approval or Consent to participate since it does not concern human participants or human or animal datasets.

Consent for publication
The authors of this manuscript consent to its publication.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.