Learning of Gaussian Processes in Distributed and Communication Limited Systems

TL;DR

This paper investigates optimal and practical communication-efficient algorithms for learning Gaussian Processes in distributed systems, demonstrating near-optimal performance with minimal communication and outperforming existing methods.

Contribution

It introduces an information-theoretic optimal scheme for inner-product estimation and practical quantization methods, improving distributed Gaussian Process learning.

Findings

Optimal vector quantization scheme derived for inner-product estimation.

Practical per-symbol quantization performs close to the optimal scheme.

Proposed methods outperform existing distributed GP learning schemes with minimal communication.

Abstract

It is of fundamental importance to find algorithms obtaining optimal performance for learning of statistical models in distributed and communication limited systems. Aiming at characterizing the optimal strategies, we consider learning of Gaussian Processes (GPs) in distributed systems as a pivotal example. We first address a very basic problem: how many bits are required to estimate the inner-products of Gaussian vectors across distributed machines? Using information theoretic bounds, we obtain an optimal solution for the problem which is based on vector quantization. Two suboptimal and more practical schemes are also presented as substitute for the vector quantization scheme. In particular, it is shown that the performance of one of the practical schemes which is called per-symbol quantization is very close to the optimal one. Schemes provided for the inner-product calculations are incorporated into our proposed distributed learning methods for GPs. Experimental results show that with spending few bits per symbol in our communication scheme, our proposed methods outperform previous zero rate distributed GP learning schemes such as Bayesian Committee Model (BCM) and Product of experts (PoE).

Equations

$$y_i=f(x_i)+{ϵ}_i,{ϵ}_i\sim \mathcal{N}(0,{\sigma }_{ϵ}^2).$$

$$m(x)$$

$$cov(x,x^{\prime })$$

$$k(x_i,x_j)=a{x}_i^T{x}_j+b,$$

$$G=\left[\begin{array}{cc}{G}_{11}& G_{12}\\ G_{21}& G_{22}\end{array}\right],$$

$$D$$

$$D=\frac{1}{n}\sum _{i=1}^n(x_i-{\hat{x}}_i{)}^T{S}_y(x_i-{\hat{x}}_i),$$

$$\underset{(\varphi ,\psi )}{\min }\mu R+(1-\mu )D,\mu \in [0,1].$$

$$\mathbb{P}(X,Y)$$

$$=\prod _{i=1}^n\mathbb{P}(x_i)\prod _{i=1}^n\mathbb{P}(y_i).$$

$${\mathbb{E}}_{\mathbb{P}(y)}\{D\}$$

$$R$$

$${\mathbb{E}}_{P(x,\hat{x})}\{d(x,\hat{x})\}\le D,$$

$$d(x,\hat{x})$$

$$R_{lb}(D)=h(x)-\frac{1}{2}\log (2\pi e{)}^d\det (Q_y^{-1})-\frac{1}{2}\sum _{i=1}^d\log q_i,$$

$$q_{i}=\left\{\begin{array}[]{lr}\lambda&\quad if\quad\lambda<\Lambda_{ii}\\ \Lambda_{ii}&\quad if\quad\lambda\geq\Lambda_{ii}\end{array}\right.,$$

$$I(x;\hat{x})$$

$$=h(x)-h(x-\hat{x}\mid \hat{x})$$

$$\stackrel{\text{(a)}}{\ge }h(x)-h(x-\hat{x})$$

$$\stackrel{\text{(b)}}{\ge }h(x)-\frac{1}{2}\log (2\pi e{)}^d\det (Q)$$

$$\underset{Q}{\min }\log \det (Q)$$

$${\mathbb{E}}_{P(x,\hat{x})}\{d(x,\hat{x})\}=\mathrm{tr}(Q{Q}_y)\le D,$$

$$Q_x-Q⪰0.$$

$$\mathcal{L}=\log \det (Q)-\lambda \mathrm{tr}(Q{Q}_y)-\mathrm{tr}\left(A^T(Q-Q_x)\right),$$

$$\lambda \ge 0,A⪰0$$

$$\frac{\partial \mathcal{L}}{\partial Q}=Q^{-1}-\lambda Q_y-A$$

$$\mathrm{tr}\left(A^T(Q-Q_x)\right)$$

$$\mathrm{tr}(Q{Q}_y)$$

$$A^T(Q-Q_x)=0.$$

$$(Q^{-1}-\lambda Q_y)(Q-Q_x)$$

$$Q_y^{-1/2}\left(Q^{-1}-\lambda Q_y\right)Q_y^{-1/2}{Q}_y^{1/2}\left(Q-Q_x\right)Q_y^{1/2}=0,$$

$$\begin{gathered} \left({\left(Q_y^{1/2}Q{Q}_y^{1/2}\right)}^{-1}-\lambda I\right) \\ \left(Q_y^{1/2}Q{Q}_y^{1/2}-Q_y^{1/2}{Q}_x{Q}_y^{1/2}\right)=0. \end{gathered}$$

$$Q_y^{1/2}{Q}_x{Q}_y^{1/2}=U\Lambda U^T,$$

$$\begin{gathered} U^T\left({\left(Q_y^{1/2}Q{Q}_y^{1/2}\right)}^{-1}-\lambda I\right)U{U}^T \\ \left(Q_y^{1/2}Q{Q}_y^{1/2}-U\Lambda U^T\right)U=0. \end{gathered}$$

$$({\tilde{Q}}^{-1}-\lambda I)(\tilde{Q}-\Lambda )=0.$$

$$\tilde{Q}$$

Full text

1 Introduction

Complex systems are ubiquitous; biological, social, transportation networks are examples that we face them on a daily basis. Managing and controlling such systems require learning from data observed and measured at possibly distinct points of the systems. Extracted Data have several intrinsic features: they contain various forms of attributes which are noisy and incomplete; they are massive and usually generated distributively across the networks. Effective and efficient utilization of available resources plays an important role in learning, managing, and controlling of such complex systems. Through a distributed and joint communication-computation design, this paper aims at presenting important machine learning instances where improvement in performance compared to conventional methods can be analytically and experimentally justified.

Many frameworks are developed to handle complex and distributed datasets efficiently. Hadoop Distributed File Systems (HDFS) [1], Google File System (GFS) [2] are examples of distributed storage systems. MapReduce [3] and Spark [4] are well known frameworks for distributed data processing systems. These platforms can be used to implement machine learning algorithms naturally [5].

Distributed processing and storage systems are designed to be universal by separating computation and communication tasks. The aim of communication is to exchange information between distributed machines efficiently, whereas the aim of computation is to run a distributed algorithm with the given communication platform. However, communication can be adjusted and modified based on the problem at hand achieving superior results. A new trend in research is started to design distributed and joint communication-computation optimal systems [6].

In this paper, we present some important machine learning instances in distributed settings and determine analytically and experimentally how one can obtain the optimal trade-off between performance and communication cost. The problem that we focus on is learning of Gaussian Processes (GPs) which are fundamental models applicable to many machine learning tasks. We propose efficient methods for GP learning and analyze the performance and communication cost of the proposed methods. Although GPs can be used for both classification and regression, we only consider regression in presented applications and experiments. However, the proposed methods are also applicable to other learning models such as distributed manifold learning, KNN-based algorithms, etc.

There exist abundant studies related to distributed learning. Here, we address a few of them which we have found closest to the setting considered in this paper. From information theoretic point of view, distributed statistical inference is addressed by Ahlswede [7], Amari and Han [8, 9, 10]. In their settings, data are assumed to be distributed between two machines and through minimal communication, it is desired to test a hypothesis or estimate a parameter optimally. They used the method of types [11] to obtain some useful error and communication bounds on several problems. A survey by Han et. al. [8] summarizes all early studies in that direction where hypothesis testing, parameter estimation, and classification problems in distributed settings are addressed.

Recently, several researchers worked on distributed statistical learning problems and used tools from information theory to obtain bounds on the performance of the algorithms: [12, 13, 14, 15]. For instance, Duchi et. al. [14] investigated the minimum amount of communication required for achieving centralized minimax risk in certain estimators. Rahman et. al. [12] analyzed the problem of distributed hypothesis testing from information theoretic perspective and prove the optimality of binning method for conditionally independent sources.

There are papers addressing the problem of distributed statistical inference from machine learning point of view [16, 17, 18, 14, 19]. For instance, Liu in [16] showed that how the maximum likelihood estimation of parameters of exponential families can be obtained distributively. Meng et. al. [17] proposed a distributed algorithm for learning of inverse covariance matrix in Gaussian graphical models. Broderick et. al. [18] proposed a distributed algorithm for calculating the posterior distribution using variational inference over distributed datasets. Authors of [19] have addressed the kernel PCA in distributed systems on a large amount of data.

Distributed machine learning problems are also related to distributed optimization problems as optimization of an objective function during their training or test procedure is part of learning tasks. Therefore, several studies on distributed optimization algorithms [20] can directly or indirectly be applied to many machine learning models. Alternating Direction Method of Multipliers (ADMM) [20] is one of the most successful distributed algorithms for solving optimization problems over distributed datasets. There are many papers that are based on the ADMM algorithm [21, 22, 23]. For example, Forero et. al. [22] proposed a re-parameterization for SVM learning in distributed manner based on ADMM.

As mentioned before, in this paper we have considered GP learning for distributed systems. There are many works that propose methods for learning GPs distributively [24]. Deisenroth et. al. [24] presented a distributed model for GP learning which is based on models called product of experts (PoE). PoE-based models assume that the local dataset in each machine is independent of the other local datasets. Thus, the marginal likelihood factorizes into product of individual local marginal likelihoods. In these models, the overall gram matrix will be block diagonal where the blocks are the local gram matrices of the machines. In PoE models, after training the model hyper-parameters, predictions of GP experts are combined. There are several combining methods such as PoE [25], Bayesian Committee Machine (BCM) [26] and the recently proposed method robust BCM (rBCM) [24].

PoE models need no data transmission during training which makes them favorable when communication is very restricted. However, assuming independence of local experts is unfavorable in many applications. For example, in RBF kernels that are non-increasing function of Euclidean distances of inputs, such assumption may decrease the performance of PoE models. This is due to the fact that data is distributed randomly between machines and nearby inputs may be located at distinct machines. In this case, transmission of inputs is necessary to achieve the optimal performance of the learning task.

There are other notable studies on localized and distributed learning of sparse GPs such as [27, 28, 29, 30, 31]. The main idea behind these works is to select or extract $m\ll n$ inducing points from the training datasets so that the trained sparse GP is close to the full GP. Authors of [29] proposed a re-parameterization on [27] and [32] for applying sparse GPs on a network of multiple machines. This method for obtaining inducing points, in each iteration of the optimization procedure transmits all inducing and hyperparameters of the GPs to all machines. Due to the size of the inducing set and the number of iterations till the convergence, the communication cost of this method is usually very high. Moreover, the exact communication cost of this method is not calculable before running the training procedure.

This paper is organized as follows. Section 2 briefly describes background on GPs for the regression problem. In Section 3, we formally define the main problem of the paper. In Section 4, we propose three methods for solving the inner-product estimation problem in a distributed system. In section 5, we show how the proposed methods of Section 4 can be used for distributed learning of GPs. Finally, in Section 6, we evaluate our methods over synthetic and real datasets.

2 Gaussian Processes

Statistical machine learning models can be divided into parametric and non-parametric models. GPs are powerful non-parametric models that are used in many supervised learning tasks such as classification and regression [31]. Basically, a GP is a stochastic process that any finite collection of its points is jointly Gaussian. To learn a GP model, it is sufficient to estimate its mean and covariance functions. The covariance matrix of a GP is sometimes called gram matrix which is denoted by $G$ through out this paper.

In regression problems, we wish to estimate a real-valued function $f(x)$ where $x\in\mathbb{R}^{d}$. We can use a GP as a prior distribution over these real-valued functions. Suppose our training dataset is $\{(x_{i},y_{i})\}_{i=1}^{n}$, which we call $\{x_{i}\}$ input samples and $\{y_{i}\}$ target observations. The training examples are noisy realization of the target function $f(x)$ where each $y_{i}$ is obtained by adding a Gaussian noise to $f(x)$ at input $x_{i}$ as follows:

[TABLE]

The posterior process is a GP with the following mean and covariance functions:

[TABLE]

where $G_{xn}=[k(x,x_{1}),\cdots,k(x,x_{n})]$ is an $n$-dimensional row vector of kernel values between $x$ and the training inputs, $\mathbf{y}$ is the vector of all target observations $\mathbf{y}=[y_{1},\cdots,y_{n}]^{T}$, and $G_{nx}=G_{xn}^{T}$. Then, the predictive distribution for a new input $x_{*}$ is $\mathbb{P}(y_{*}|x_{*},\mathbf{y})=\mathcal{N}\left(m(x_{*}),cov(x_{*},x_{*})+\sigma_{\epsilon}^{2}\right)$. Note that $G$ is an $n\times n$ matrix where $G_{ij}$ is the correlation of $f(x_{i})$ and $f(x_{j})$.

In the training procedure, the hyperparameters of the model (e.g. $\sigma_{\epsilon}^{2}$, and parameters of the kernel function) need to be trained. Usually, the hyperparameters are obtained by maximizing the marginal likelihood $\mathbb{P}(\mathbf{y})$. Although there exist many kernel functions such as linear, polynomial, squared exponential, Matern, and etc [31], we select the following linear kernel function:

[TABLE]

where $a$ and $b$ are the hyperparameters of the kernel. Although our proposed methods for GP learning is based on the linear kernel, they can be readily applied on any Radial Basis Function (RBF) kernels as experimentally will showed in Section 6. It is worth mentioning that kernels that are only a function of $\|x_{i}-x_{j}\|$ are called RBFs.

3 Problem Definition

Suppose we have only two machines in our distributed system so that each machine stores some portion of the whole dataset. Then the gram matrix can be partitioned as follows:

[TABLE]

where $G_{11}$ (resp. $G_{22}$) is the gram matrix of data stored in the first (resp. second) machine and $G_{12}$ is the cross gram matrix between data of the first and second machines. Moreover, $G_{12}=G_{21}^{T}$. Note that $G_{11}$ and $G_{22}$ can be calculated locally without any data transfer between machines. However, $G_{12}$ and its transpose $G_{21}$ must be calculated cooperatively by the two machines which requires intensive communication.

Let the first machine be responsible for learning the GP model. This machine calculates $G_{11}$ locally, whereas $G_{12}$ and $G_{22}$ need to be calculated using messages received from the second machine. In Section 5, we will show that if $G_{12}$ is available, then $G_{22}$ can be efficiently estimated. Thus, our main goal is to find an efficient method to obtain $G_{12}$ in distributed manner with minimum distortion.

In fact, we seek for an algorithm that has low communication cost and at the same time has low distortion in reconstruction of the target matrix. Without any limitation on communication cost, we gather the whole dataset at the first machine which leads to the best performance with zero distortion. In Section 4, we propose novel methods respecting imposed communication constraints and achieving optimal/suboptimal performance in the learning procedure.

In our basic setting, the gram matrix of GP is first estimated by communicating information bits between machines. The estimated gram matrix is then used to learn the latent function. Therefore, minimum distortion in estimating of the gram matrix is desired. By considering a linear kernel function for the GP model, our goal is to calculate the inner product of two datasets that stored in different machines with minimum distortion.

To put it formally, consider two datasets $X=\{x_{1},\cdots,x_{n}\}$ and $Y=\{y_{1},\cdots,y_{n}\}$ stored in machines $\mathcal{M}_{x}$ and $\mathcal{M}_{y}$, respectively. Samples $x_{i}$ and $y_{i}$ are $d$-dimensional random vectors with zero means. We would like to calculate the inner product $\langle x_{i},y_{j}\rangle$ for any $x_{i}\in X$ and $y_{j}\in Y$ in machine $\mathcal{M}_{y}$. Exact calculation of inner products requires transmission of the whole dataset $X$ to machine $\mathcal{M}_{y}$ which is not efficient in many practical situations due to its communication cost. Lowering the communication cost, however, results in reduction in the quality of computed inner products. Therefore, finding the tradeoff between the communication cost and the quality of calculated inner products plays an important role in the design of practical systems.

To characterize the tradeoff, we limit the number of transmission bits per observation to $R$ bits; leading to the total of $nR$ transmission bits. These bits are used to design an encoding map $\phi$ from X to an index set $\mathcal{I}=\left\{1,2,\cdots,2^{nR}\right\}$. The encoded index is transmitted to machine $\mathcal{M}_{y}$ where a decoding map $\psi$ is used to map the index to another set $\hat{X}=\{\hat{x}_{1},\cdots,\hat{x}_{n}\}$ which is called the set of reconstructed points. In other words, each $x_{i}$ is quantized to $\hat{x}_{i}$.

We would like to design $\phi$ and $\psi$ such that the distortion in computing the inner products is minimized. The distortion is measured by

[TABLE]

that is,

[TABLE]

where $S_{y}$ is the sample covariance matrix of $Y$.

Definition 1**.**

*Achievable $(R,D)$

A rate-distortion pair $(R,D)$ is achievable if there exist encoding and decoding maps $\phi$ and $\psi$ such that the transmission rate is limited by $R$ and the distortion in reconstructing the inner products is limited by $D$.*

In this paper, we are interested in characterizing the set of all achievable rate-distortion pairs. The region of all achievable rate-distortion pairs is a convex region. To show this, one needs to prove that if $(R_{1},D_{1})$ and $(R_{2},D_{2})$ are two achievable pairs, then $(\lambda R_{1}+(1-\lambda)R_{2},\lambda D_{1}+(1-\lambda)D_{2})$ is an achievable pair for any $\lambda\in[0,1]$. A coding scheme which uses both schemes for different portions of data can achieve this. Let $\lambda$ portion of data be encoded by the first scheme and the rest with the second scheme. It is easy to show that the overall rate-distortion pair $(\lambda R_{1}+(1-\lambda)R_{2},\lambda D_{1}+(1-\lambda)D_{2})$ is indeed achievable.

Extreme points of the achievable region can be obtained by solving the following optimization problem

[TABLE]

4 Transmission Schemes

In this section, we present three different schemes for transmitting the dataset $X$ to machine $\mathcal{M}_{y}$ with low communication cost and distortion. The first scheme is based on an optimal compression method which characterizes rate-distortion region theoretically for normally distributed samples. Although the scheme is optimal, it is not practical due to computation complexity required for block coding. Thus, we propose another scheme that is an approximate version of the optimal one. This method uses per symbol coding instead of block coding. Experiments show its performance is near optimal and it is applicable for many distributed settings. The third scheme is even more simpler and does not require any nonlinear map. It is based on ideas from dimension reduction schemes which are tailored to our distributed setting. This method has more communication cost than the previous two schemes, but it is very simple and fairly efficient in many applications.

4.1 Compression

Let us assume the datasets $X$ and $Y$ are drawn from two independent distributions as

[TABLE]

Taking the expectation of distortion in (7) with respect to $y$, we obtain

[TABLE]

where $Q_{y}$ is the covariance matrix of $Y$ defined by $Q_{y}=\mathbb{E}\{y^{T}y\}.$ It is worth mentioning that under the preceding conditions, the covariance matrix of $y$ is the only parameter influencing the achievable rate-distortion region and therefore the knowledge of the exact distribution of $Y$ is not required.

We would like to characterize rate-distortion region for the above setting. This problem is an instance of the rate-distortion problem in information theory where the distortion is defined by (7). The optimal coding scheme achieving points on the boundary of the achievable region is obtained by solving the following optimization problem [11]:

[TABLE]

where

[TABLE]

Suppose $P^{*}(\hat{x}|x)$ is the solution of (11). The coding scheme with the rate $R$ and distortion $D$ is obtained as follows:

1. Generating Codebook: Calculate $P^{*}(\hat{x})=\int P(x)P^{*}(\hat{x}|x)dx$. Generate a codebook consisting of $2^{nR}$ sets of size $n$ that are drawn i.i.d. from $\prod_{i}P^{*}(\hat{x}_{i})$. We denote the generated random sets by $\hat{X}_{1},\cdots,\hat{X}_{2^{nR}}$. The codebook is available at both transmitter ($\mathcal{M}_{x}$) and receiver ($\mathcal{M}_{x}$) machines. Note that the codebook can be generated at both machines by an equal initial seed which leads to identical codebooks. Thus, there is no need to send the codebook from a machine to the other.
2. Encoding: We are given a set $X$. Index it by $j$ if there exists a set in the codebook such that the pair $(X,\hat{X}_{j})$ is strongly jointly typical (for definition of strong typicality refer to [11]). If there is more than one such $j$, send the first in lexicographic order. If there is no such $j$, send 1 to the receiver. It is clear that transmitting the index $j$ requires $nR$ bits. Since by each transmission we send $n$ data samples, the average rate will be $R$ bits per each sample.
3. Decoding: At the receiver (machine $\mathcal{M}_{y}$ in our case), select $\hat{X}_{j}$ as the reproduced set for $X$.

For the proof that this coding scheme asymptotically (when $n\to\infty$) has rate $R$ and distortion $D$ see [11]. Thus, by solving the optimization problem (11) and using the above coding scheme, we obtain the optimal method for transmitting the set $X$ with rate $R$ and distortion $D$. Unfortunately, in many cases, there is no closed form solution for (11), but we can obtain a lower bound for it (Theorem 1).

The following theorem presents a lower bound for rate-distortion function (11).

Theorem 1**.**

Let $\{x_{1},\cdots,x_{n}\}$ be $n$ $i.i.d.$ $d$-dimensional zero mean random vectors with covariance matrix $Q_{x}$, and let the distortion measure be $d(X,\hat{X})=\sum_{i=1}^{n}(x_{i}-\hat{x}_{i})^{T}Q_{y}(x_{i}-\hat{x}_{i})$, where $Q_{y}$ is a symmetric positive definite matrix. Then the rate distortion function (11) is lower bounded by

[TABLE]

where

[TABLE]

$h(x)$* is entropy of variable $x$, $\Lambda$ is the diagonal matrix of eigenvalues of $Q_{x}Q_{y}$, and $\lambda$ is chosen so that $\sum_{i=1}^{d}q_{i}=D$.*

Proof.

The lower bound is obtained as follows. Let $Q=\mathbb{E}\{(x-\hat{x})(x-\hat{x})^{T}\}$, we have

[TABLE]

where (a) follows from the fact that conditioning does not increase the entropy; (b) follows from the fact that Gaussian distribution maximizes the entropy conditioned on a covariance matrix.

One can find the best lower bound in (15) by minimizing $\log\det(Q)$ over all appropriate covariance matrices $Q$. This amounts to solving the following optimization problem:

[TABLE]

where the second constraint comes from the fact that independence of $u=x-\hat{x}$ and $\hat{x}$ results in an equality in (a). More precisely, $x=u+\hat{x}$ where $u$ and $\hat{x}$ are independent; this yields $Q_{x}=Q_{u}+Q$ and thus, $Q_{x}\succeq Q$.

It can be shown that the above problem is convex and therefore KKT conditions provide necessary and sufficient conditions for the optimal solution. To obtain the KKT conditions, we first obtain the Lagrangian as

[TABLE]

where $\lambda$ and $A$ are the Lagrange multipliers. Using $\mathcal{L}$, the KKT conditions amount to

[TABLE]

Since $A^{T}(Q-Q_{x})$ is positive semi-definite, (19) yields

[TABLE]

By substituting $A$ from (18) in (21), we obtain:

[TABLE]

Since $Q_{y}$ is positive definite, it can be decomposed as $Q_{y}=Q_{y}^{1/2}Q_{y}^{1/2}$ where $Q_{y}^{1/2}$ is the symmetric squared root of $Q_{y}$. Multiplying (22) by $Q^{-1/2}$ from left, $Q^{1/2}$ from right and $Q_{y}^{-1/2}Q_{y}^{1/2}$ from middle, yields

[TABLE]

that is

[TABLE]

We can decompose $Q_{y}^{1/2}Q_{x}Q_{y}^{1/2}$ as follows

[TABLE]

where $U$ is a unitary matrix. Since the eigenvalues of $Q_{y}^{1/2}Q_{x}Q_{y}^{1/2}$ are equal to eigenvalues of $Q_{x}Q_{y}$, $\Lambda$ is a diagonal matrix which includes eigenvalues of $Q_{x}Q_{y}$. By substituting (25) in (24) we can rewrite (24) as follows

[TABLE]

Here, (24) is multiplied by $U^{T}$ from left, $U$ from right, and $UU^{T}$ from the middle. Defining $\tilde{Q}=U^{T}Q_{y}^{1/2}QQ_{y}^{1/2}U$, from (26) obtains

[TABLE]

The above equation forces $\tilde{Q}$ to be diagonal. Each element on diagonal of $\tilde{Q}$ is whether $1/\lambda$ or $\Lambda_{ii}$. Clearly, since $Q=Q_{y}^{-1/2}U\tilde{Q}U^{T}Q_{y}^{-1/2}$, any solution for $Q$ can be mapped uniquely to $\tilde{Q}$ and vice versa. It is easy to show that $\tilde{Q}$ which is defined as follows satisfies all the KKT conditions and therefore, after mapping it to $Q$, it is an optimal solution for (16)

[TABLE]

where $\lambda$ is chosen so that:

[TABLE]

Finally, by substituting $Q=Q_{y}^{-1/2}U\tilde{Q}U^{T}Q_{y}^{-1/2}$ in (15) we obtain the lower bound (13). ∎

Theorem 1 gives a lower bound on the rate distortion function (11). In some cases such as the Gaussian distribution the lower bound is achievable (see Theorem 2), but it is not the case in general. Although this theorem does not propose a coding scheme for data compression, it gives some useful clues for designing low communication rate coding schemes. According to the lower bound $R_{lb}(D)$ in (13), eigenvalues of the product matrix $Q_{x}Q_{y}$ plays an important role in bit rate allocation. In other words, dimensions with larger eigenvalues need more bits compared to dimensions with small eigenvalues. This intuition lead us to propose an approximate coding scheme in the next section for Gaussian datasets.

The following theorem shows that the lower bound (13) is achievable when the dataset $X$ is drawn from a normal distribution. In this case, the achievable rate distortion region can be fully characterized as the following theorem presents.

Theorem 2**.**

If $\{x_{1},\cdots,x_{n}\}$ are drawn $i.i.d.$ from a normal distribution $\mathcal{N}(0,Q_{x})$, then the lower bound (13) is achievable if we choose $\{\hat{x}_{i}\}$ as follow:

[TABLE]

where $\hat{x}_{i}$ and $z_{i}$ are independent, and $Q$ is chosen so that $\operatorname{tr}(QQ_{y})=D$ and $Q\preceq Q_{x}$.

Proof.

To find the conditional distribution $\mathbb{P}(\hat{x}|x)$ that achieves the lower bound (13), it is usually more convenient to use conditional distribution $\mathbb{P}(x|\hat{x})$. If we choose the joint distribution of $x_{i}$ and $\hat{x}_{i}$ using (32), then

[TABLE]

By setting $Q=Q_{y}^{-1/2}U\tilde{Q}U^{T}Q_{y}^{-1/2}$, where $\tilde{Q}$ and $U$ is determined by (27) and (25), respectively; it is easy to show that $Q$ satisfies $\operatorname{tr}(QQ_{y})=D$ and $Q\preceq Q_{x}$. Then, we have

[TABLE]

which is equal to the lower bound (13). ∎

The theorem 2 shows that the optimal compression for Gaussian distribution must be performed using the product of covariance matrices $Q_{x}$ and $Q_{y}$. We will see that this result is similar to Theorem 3 where the optimal linear dimension reduction is performed by product of the sample covariance matrices $S_{x}$ and $S_{y}$. Theorem 2 is in fact the well known reverse-water filling problem on eigenvalues of matrix $Q_{x}Q_{y}$. In other words, if we decompose $Q_{y}^{1/2}Q_{x}Q_{y}^{1/2}$ as $U\Lambda U^{T}$ where $\Lambda$ is a diagonal matrix of eigenvalues of $Q_{x}Q_{y}$, then we should consume more bits for dimensions with larger corresponding eigenvalue $\Lambda_{ii}$. According to the theorem, for dimensions with eigenvalues less than a threshold $\lambda$ no bits are allocated (i.e. they are not encoded and transmitted), and for dimensions greater than the threshold we consume bits to reduce distortion for them. In fact, $q_{i}$ in Theorem 2 is the distortion value for $i$-th dimension after compression.

According to Theorem 2, the machine $\mathcal{M}_{x}$ needs $Q_{y}$ for compressing dataset $X$, thus, if $Q_{y}$ is not known for the machine $\mathcal{M}_{x}$, transmitting it has $\mathcal{O}$($d^{2}$) bits communication cost. After compression, transmitting the compressed dataset $\hat{X}=\{\hat{x}_{1},\cdots,\hat{x}_{n}\}$ requires $\mathcal{O}$($nR$) bits.

Although the compression method proposed by Theorem 2 is optimal and simulation results in section 6 (see Figure 2) shows that it significantly outperforms the dimension reduction method, it is unfortunately not practical in real world due to exponential computation complexity for encoding. To overcome this issue, in the next section we have proposed an applicable approximate method for compression which has near optimal performance.

4.2 Per-Symbol Compression Scheme

Although the proposed method for compression in Theorem 2 is optimal for normal data samples, it is not practically favorable to implement in many cases as it requires vector quantizers which are designed to work on blocks of data samples. Thus, a per-symbol quantizer which is easier to implement is favorable.

In this section, we propose a simple yet effective quantizer for normal data samples which operates near optimal in many practical situations. In the proposed method, we first transform samples of the dataset $X$ by $U^{T}Q_{y}^{1/2}$, where $Q_{y}^{1/2}$ is the squared root of $Q_{y}$, and $U$ is obtained by solving the following singular value decomposition problem:

[TABLE]

It is easy to see that the above transformation makes dimensions of samples $\{x_{i}\}$ uncorrelated. Hence, for Gaussian samples the dimensions are independent. We denote the transformed dataset by $X^{\prime}=\{x_{1}^{\prime},\cdots,x_{n}^{\prime}\}$. Thus, we can consider each $x_{i}^{\prime}$ as $d$ independent normal random variables. Finally, we quantize each dimension separately. Our coding scheme is as follows:

1. Encoder: Transform dataset $X=\{x_{1},\cdots,x_{n}\}$ to $X^{\prime}=\{x_{1}^{\prime},\cdots,x_{n}^{\prime}\}$ using $x_{i}^{\prime}=U^{T}Q_{y}^{1/2}x_{i}$. Quantize each dimension of $\{x_{i}^{\prime}\}$ samples using a scalar quantizer. Transmit quantized samples denoting by $\hat{X}^{\prime}=\{\hat{x}_{1}^{\prime},\cdots,\hat{x}_{n}^{\prime}\}$ to the receiver.
2. Decoder: Transform received data samples using

[TABLE]

to reconstruct $\hat{x}_{i}$.

In fact, we quantize the transformed version of samples in the encoder and apply the inverse transform on received samples at the decoder. The expected distortion of this method is calculated by taking expectation from (7) as follows:

[TABLE]

where $Q^{\prime}$ is the covariance matrix of $x^{\prime}-\hat{x}^{\prime}$. Since the dimensions of $x^{\prime}$ are independent for normal distribution and the quantization is done separately on each dimension, $Q^{\prime}$ is diagonal, hence by denoting $Q^{\prime}=diag(\sigma_{1},\cdots,\sigma_{d})$ the expected distortion is $\sum_{i=1}^{d}\sigma_{i}$. In fact, $\sigma_{i}$ is the distortion of the $i$-th dimension. If we want to quantize each sample by $R$ bits, we should allot more bits to the dimensions with larger distortions. More precisely, for finding the optimal bit allocation to each dimension, the following integer optimization problem should be solved:

[TABLE]

where $R_{i}$ is the number of bits allocated to $i$-th dimension. For solving the above minimization problem, $\{\sigma_{i}\}$ be calculated first.

To relate $\{\sigma_{i}\}$ and $\{R_{i}\}$, we use a simple scalar quantizer for each dimension. Assume we have a scalar random variable $u\sim\mathbb{P}(u)$ which we want to quantize it by $R$ bits. We create $2^{R}$ equally probable bins over the real axis and use bins centroids as the reproduction points at the decoder (receiver). To compress $u$, we send index of the bin that $u$ belongs to it by $R$ bits.

Assume we denote the $i$-th bin interval by $(a_{i},a_{i+1})$ and its centroid by $c_{i}$, then

[TABLE]

where $\{a_{i}\}$ are determined so that $\int_{a_{i}}^{a_{i+1}}P(u)du=2^{-R}$. If $u\sim\mathcal{N}(0,\sigma_{u}^{2})$, then

[TABLE]

If $\{\alpha_{i}\}$ is the set of bins boundaries for standard normal distribution, then the boundaries for $u$ is obtained by $a_{i}=\sigma_{u}\alpha_{i}$. Note that the $\{\alpha_{i}\}$ is only function of $R$. Hence, we can calculate $c_{i}$ using

[TABLE]

The above equation shows that the centroids of any normal distribution are simply obtained by multiplying the centroids of standard normal distribution by standard deviation $\sigma_{u}$.

The expected reconstruction error for this scalar quantizer is calculated as follows:

[TABLE]

where $\sigma_{c}^{2}$ is variance of the centroids.

In our problem setting, from (37) and (42), the distortion of $i$-th dimension is $\sigma_{i}(R_{i})=e(\Lambda_{ii},R_{i})$, where $R_{i}$ is the number of bits allocated to the $i$-th dimension and $\Lambda_{ii}$ is the variance of $i$-th dimension of the transformed dataset $X^{\prime}$. Now, we can solve the integer program in (38). Defining $\Delta\sigma_{i}=e(\Lambda_{ii},R_{i})-e(\Lambda_{ii},R_{i}+1)$ as the amount of decrease in distortion of $i$-th dimension by adding a new bit, then the greedy algorithm presented in Algorithm 1, obtains the optimal bit allocation. In this algorithm, we allocate one bit in each iteration to the dimension with largest $\Delta\sigma$ and update the distortion of that dimension. We iterate this procedure until all available bits are allocated.

To show that the above greedy bit allocation algorithm is optimal, we rely on the fact that $\Delta\sigma$ is a decreasing function of $R$. The plot of $\Delta\sigma$ against the bit rate $R$ shows that it is indeed decreasing exponentially (Figure 1). Note that the figure is plotted for standard normal distribution and behaves similarly for any normal distribution.

We denote the optimal bit allocation by $S_{opt}=(R_{1}^{*},\cdots,R_{d}^{*})$ and the greedy solution by $(R_{1}^{g},\cdots,R_{d}^{g})$. We define the distance between the greedy and optimal solutions as follows:

[TABLE]

If the problem has multiple optimal solutions, we select a solution with the minimum distance value $r_{min}$. If $r_{min}=0$, then the greedy solution is also optimal. Hence, we assume that $r_{min}>0$ and show that it leads to a contradiction.

If $r_{min}>0$, then there exist two bit rates $R_{i}^{g}$ and $R_{j}^{g}$ so that $R_{i}^{g}>R_{i}^{*}$ and $R_{j}^{g}<R_{j}^{*}$. Assume that in iteration $t$ of the greedy algorithm, the bit rate of dimension $i$ exceeds $R_{i}^{*}$, and the bit rate of dimension $j$ is $R_{j}^{t}$. Thus,

[TABLE]

On the other hand, it is clear that $R_{j}^{t}\leq R_{j}^{g}<R_{j}^{*}$. Hence,

[TABLE]

From (44) and (45), we have

[TABLE]

Thus, by removing a bit from $R_{j}^{*}$ and adding it to $R_{i}^{*}$ we obtain a new solution $S_{new}$ which its error value is less than or equal to error of $S_{opt}$. Since $S_{opt}$ is optimal, $S_{new}$ could not have smaller error value. But its distance with the greedy solution is $r_{min}-2$ which is a contradiction.

The proposed per-symbol coding scheme has $\mathcal{O}$($R$) time complexity for encoding (using binary search algorithm over the codebook) and $\mathcal{O}$($1$) for decoding. However, its space complexity is $\mathcal{O}$($2^{R}$) that is required for storing the bins intervals and centroids (codebook). This method requires to transmit matrices $Q_{x}$ and $Q_{y}$ between machine $\mathcal{M}_{x}$ and machine $\mathcal{M}_{y}$. Thus, the overall communication cost is $\mathcal{O}$($2d^{2}+Rn)$ bits.

Although in the proposed per-symbol method the underlying distribution of samples of $X$ is considered to be normal, we can use it for other distributions. However, the normal distribution is used to make the dimensions independent and to prove the optimality of the greedy bit loading. Experimental results in Section 6 reveal that the performance of this method for other distributions is desirable.

4.3 Dimension Reduction

The third approach proposed in this paper for reducing communication cost for transmitting the dataset $X$ is to represent them in a lower dimensional space. In this section we propose a linear dimension reduction method that efficiently minimizes the inner product distortion. Interestingly, PCA, the conventional dimension reduction method, falls from optimality as it is shown in Section 6.

Our goal is to construct a lower dimensional vector $\hat{x}_{i}$ for $x_{i}$ so that the error function (7) be minimized. In order to represent $x_{i}$ by $\hat{x}_{i}$ which lies in a lower dimensional space, $\hat{x}_{i}$ is defined as a weighted sum of some basis vectors as follows

[TABLE]

where $\{u_{k}\}$ are orthogonal $d$-dimensional basis vectors so that $\|u_{k}\|^{2}=1$, $z_{ik}$ is the weight of vector $u_{k}$ for $\hat{x}_{i}$, and $m$ is dimension of the new space. The equation (47) can be rewritten in vector-matrix form as

[TABLE]

where $z_{i}=[z_{i1},\cdots,z_{im}]^{T}$, and $U=[u_{1},\cdots,u_{m}]$ is a $d\times m$ matrix that $u_{k}$ forms the $k$-th column, thus, $U^{T}U=I$. In fact, $z_{i}\in\mathbb{R}^{m}$ is the corresponding representation of $x_{i}$ in the $m$-dimensional space. We would like to minimize the error function (7), with respect to the matrix $U$ and vectors $z_{i}$, $i=1,\cdots,n$. More precisely, we want to find the solution of the following optimization problem:

[TABLE]

Following theorem gives the optimal values of this problem.

Theorem 3**.**

The optimal basis vectors for minimizing $D$ in (7) are the first $m$ right eigenvectors of matrix $S_{x}S_{y}$ corresponding to $m$ largest eigenvalues, where $S_{x}$ and $S_{y}$ are sample covariance matrices of $X$ and $Y$, respectively. The vector $z_{i}\in\mathbb{R}^{m}$ is obtained by

[TABLE]

and the distortion value is the sum of $(d-m)$ smallest eigenvalues of $S_{x}S_{y}$.

Proof.

Substituting (48) in (7) yields

[TABLE]

We first optimize $D$ with respect to $z_{i}$ considering $U$ as a constant matrix. To this end, by taking the derivative with respect to $z_{i}$ and setting it to zero, we obtain

[TABLE]

where

[TABLE]

By substituting (52) in (51), we obtain the error function $D$ as a function of $U$ as

[TABLE]

The first term in (54) is constant and can be ignored. Next, we obtain $U$ by introducing the Lagrange multiplier $B$ as follows:

[TABLE]

where $B$ is an $m\times m$ symmetric matrix that comes from the constraint $U^{T}U=I$. Taking derivative from (55) with respect to $U$, yields

[TABLE]

By setting it to zero and multiplying by $U^{T}$ from the left, yields $B=0$. Hence, we have

[TABLE]

Multiplying by $S_{y}^{-1}$ and $A^{-1}$ from left and right, respectively, yields

[TABLE]

It is easy to show that replacing $S_{x}S_{y}U$ with $U\Lambda$ where $\Lambda$ is a diagonal matrix which includes $m$ eigenvalues of $S_{x}S_{y}$, satisfies (58). In other words, if columns of $U$ are the right eigenvectors of $S_{x}S_{y}$, then it satisfies (58). Thus, we would like to solve the following eigenvector problem:

[TABLE]

By substituting (59) in (54), we obtain

[TABLE]

Clearly for minimizing $D$, the trace of $\Lambda$ should be maximized, thus, the $m$ right eigenvectors of $S_{x}S_{y}$ that correspond to the $m$ largest eigenvalues must be selected as the columns of $U$. ∎

The above theorem shows that the optimal dimension reduction can be obtained from factorization of $S_{x}S_{y}$. Therefore, $S_{y}$ is first transmitted to machine $\mathcal{M}_{x}$. Transmitting $S_{y}$ from machine $\mathcal{M}_{y}$ requires communicating $\mathcal{O}$($d^{2}$) bits. After dimension reduction, transmitting $\{z_{i}\}$ and $U$ to the machine $\mathcal{M}_{y}$ requires $\mathcal{O}$($mn+md$) bits, thus, the overall communication cost is $\mathcal{O}$($d^{2}+mn+md$) bits. Theorem 3 shows that if the dimensions in $Y$ dataset are highly correlated, then the error value is significantly lower than PCA and conversely for loosely correlated dimensions the method is close to the PCA so that if $S_{y}=I$ then it will be completely equivalent to PCA.

One more interesting observation is that the solution to dimension reduction in Theorem 3 is very similar to that of Theorem 1 in the sense that both rely on $S_{x}S_{y}$ (or $Q_{x}Q_{y}$) for choosing their best strategy.

It is worth mentioning that if $S_{x}\simeq S_{y}$ (e.g. $X$ and $Y$ have identical distributions) there is no need to transmitting $S_{y}$, thus the communication cost reduces by $\mathcal{O}$($d^{2}$) bits. In this case, the machine $\mathcal{M}_{x}$ can perform the proposed dimension reduction by estimating $\hat{S}_{y}=S_{x}$. But in this case, our proposed method is equivalent to the PCA. Because, eigenvectors of $S_{x}^{2}$ are identical to that of $S_{x}$ and the eigenvalues are squared. Since $S_{x}$ is positive definite, order of eigenvalues for $S_{x}$ and $S_{x}^{2}$ are identical. This shows that our proposed method for dimension reduction outperforms the PCA when the local datasets on the machines differ from each other. This will be shown experimentally in Section 6.

5 Distributed Learning of GPs

In this section, we make use of the three methods proposed in Section 4 to learn GP models over distributed datasets. Two communication models are considered: single-center and broadcast models. In the first model, we assume that there is a centeral machine in the network which is responsible for learning the final GP model. In the second model, we assume that each machine in the network has a local GP model which is obtained from local dataset and the received data from all other machines. In this model, we also assume that each machine broadcasts its encoded data and all the machines receive the transmitted data.

Suppose the dataset $\mathcal{D}=\cup_{i=1}^{m}X_{i}$ is stored in $m$ machines $\mathcal{M}_{1}$, $\cdots$, $\mathcal{M}_{m}$ such that $X_{i}$ is stored in machine $\mathcal{M}_{i}$. For convenience and without loss of generality, we have assumed that each local dataset $X_{i}$ is of size $n$. Hence, $\mathcal{D}$ includes $N=mn$ observations in total. To learn a GP model in distributed systems, we first estimate the gram matrix which can be partitioned as follows:

[TABLE]

Clearly, $\mathcal{M}_{i}$ can calculate the gram matrix of its local dataset which is the sub-matrix $G_{ii}$ in (61). However, the other parts of $G$ should be calculated cooperatively among the machines. In the PoE-based models (BCM and rBCM) the off-diagonal sub-matrices are considered to be zero. In other words, they assume that the marginal likelihood $P(\mathbf{y}|X,\theta)$ can be factorized as a product of local marginal likelihoods. The assumption of independence is not valid in many practical situations. Therefore, we aim at learning the full $G$ with the help of data communication between machines.

We propose an algorithm based on Nyström method for gram matrix approximation. In this approximation, it is assumed that only $K<N$ rows (columns) of the gram matrix is given and the other rows (columns) must be approximated. Without loss of generality, we assume that the first $K$ rows of the gram matrix is available and the rest $N-K$ rows is approximated based on the first $K$ rows. Partitioning $G$ as

[TABLE]

and denoting the first $K$ rows of $G$ by $G_{KN}$, the Nyström approximation of $G$ is defined as

[TABLE]

where $G_{NK}$ is the transpose of $G_{KN}$ [31].

It can be shown that the two matrices $G$ and $\hat{G}$ agree on the first $K$ rows. In fact, the error in estimating $G$ comes from $\hat{G}_{(N-K)(N-K)}$. It can be shown that the error in $\hat{G}_{(N-K)(N-K)}$ is the Schur complement of $G_{KK}$ [31]. It worth mentioning that, Snelson and Ghahramani in [33] showed that the Nyström approximation can be improved by making it exact in the diagonal.

5.1 Single-Center Model

In the case of single-center model, we assume that one of the machines in the network, say $\mathcal{M}_{1}$, is the center node responsible for learning of the GP model. First, the center machine sends its local covariance matrix to other machines. Having known the covariance matrix of the center machine, each machine uses one of the coding schemes proposed in Section 4 to transmit its distorted data to the center. Finally, the center machine calculates the total gram matrix and learns hyper-parameters of the model.

For example, if $\mathcal{M}_{1}$ is the center node, it calculates its local gram matrix $G_{11}$ and local sample covariance matrix $S_{1}$. Then, it sends $S_{1}$ to other machines and each machine encodes its local dataset and transmits it to $\mathcal{M}_{1}$. Machine $\mathcal{M}_{1}$ uses the received distorted data to calculate the gram matrices $G_{12},G_{13},\cdots,G_{1m}$. Then, the full gram matrix $G$ is approximated by (63). Finally, using one of the gradient based optimization methods, the hyper-parameters of the model are learned, c.f. [31]. The prediction procedure is also performed entirely at the center machine.

5.2 Broadcast Model

In broadcast networks, each machine sends its local covariance matrix to other machines. Then, each machine broadcasts its local dataset to all other machines incorporating the methods proposed in Section 4. In particular, at machine $i$, sum of the local covariance matrices of other machines (i.e. $\sum_{j=1,j\neq i}^{m}S_{j}$) is used to obtain the best coding scheme for the proposed methods in Section 4.

Despite the single-center model, the gram matrix is approximated locally at each machine in the broadcasting model. Then, a predictive distribution can be obtained by combining these local GPs. In the broadcast model, the procedure that each machine performs is similar to the center machine in the single-center model. For example, consider the machine $\mathcal{M}_{1}$; this machine can calculate its local gram matrix $G_{11}$. The sub-matrices $G_{12},G_{13},\cdots,G_{1m}$ can be estimated using one of the methods presented in section 4. Thus, $\mathcal{M}_{1}$ can construct the first $n$ rows of the gram matrix $G$ and obtain the whole matrix by the Nyström approximation. The same procedure is done in other machines. After all, each machine has an approximate version of the full GP model.

In the prediction phase, a new test observation is sent to all machines. Each machine obtains a normal predictive distribution parameterized as $\mathcal{N}_{i}(\mu_{i},\Sigma_{i})$ for the target value of the input observation. Next, the local predictive distributions are aggregated to obtain the final predictive distribution at a fusion machine. The fusion machine can be one of the machines or a separate machine. It is also possible that all machines obtain the final predictive distribution if they access to all $\mu_{i}$’s and $\Sigma_{i}$’s. For merging local predictive distributions $\mathcal{N}_{i}(\mu_{i},\Sigma_{i})$ in the fusion machine we use the following objective function to obtain the final mean and covariance matrix:

[TABLE]

In [16], it is shown that the above objective function obtains the global maximum likelihood estimator for the situation in which $\mu_{i}$ and $\Sigma_{i}$ are maximum likelihood estimators of local dataset $X_{i}$. Although in our problem setting the $\mu_{i}$’s and $\Sigma_{i}$’s are not local maximum likelihood estimators, the objective function (64) intuitively is reasonable as a fusion criterion. However, by solving (64) we obtain the optimal values for $\mu$ and $\Sigma$ as follow

[TABLE]

Thus, each machine sends its calculated mean $\mu_{i}$ and covariance matrix $\Sigma_{i}$ to the fusion center, say machine $\mathcal{M}_{1}$, and it calculates the final predictive distribution using (65) and (66).

6 Experiments

In this section, we evaluate our proposed coding schemes presented in Section 4 over synthetic and real datasets. After comparing the distortion of the coding schemes in recovering the inner products, we evaluate their performance in training GP regression models over distributed datasets.

Figure 2 shows the distortion value in (7) against the communication rate for the proposed methods in Section 4. The dataset is generated synthetically by sampling from a $20$-dimensional Gaussian distribution by a random covariance matrix. The dataset is distributed between two machines. As can be seen, the distortion of the per-symbol compression scheme is close to the lower bound on the distortion which is obtained by Theorem 1; and is significantly lower than the proposed dimension reduction method. For the distortion curve of the dimension reduction method, we have assumed that each dimension in lower space is represented by $16$ bits without any quantization distortion. Thus, the distortion for dimension reduction is only due to loss of information by the lost dimensions.

From the experiment, it is evident that after expending $70$ bits per symbol or $3.5$ bits per dimension, the distortion is almost zero for the optimal distortion curve. This value becomes $5$ bits per dimension for the per-symbol compression scheme and $12.5$ bits for the dimension reduction scheme.

In Figure 3, we have compared two different dimension reduction methods: our proposed scheme and the one based on PCA. In the first experiment, Figure 3a we have generated two datasets $X$ and $Y$ from two multivariate Gaussian distributions with different covariance matrices. Whereas, in the second experiment, Figure 3b we have drawn both datasets from single Gaussian distribution. In the first case, since the covariance matrices are different, our method has lower distortion than PCA. But in the second case, the two methods have approximately the same distortion values. The reason is that in this situation both local covariance matrices $S_{x}$ and $S_{y}$ are approximately equal and close to the true covariance matrix (specially for large datasets). Thus, the basis vectors obtained by PCA and our method are identical. However, when the sample covariance matrices $S_{x}$ and $S_{y}$ are far from each other, our method outperforms PCA.

The comparison of the two schemes on real datasets is presented in Figure 3c and Figure 3d where the popular hand-written digits dataset known as MNIST is used. In the first experiment (Figure 3c), the images of digits $6$ and $7$ are stored in the first and second machines, respectively. In the second experiment (Figure 3d), the images are split uniformly between two machines. Figure 3c shows that under the first condition our method outperforms PCA. On the other hand, Figure 3d shows that under the second condition, the two methods have equal distortion.

A simple $1$-dimensional database is used for evaluation of GP training on quantized datasets. Figure 4 shows the GP trained on a synthetic $1$-dimensional dataset of size $N=200$. In this experiment, the whole dataset is quantized by rates $R=1$ to $R=8$ bits. The quantization method explained in Section 4.2 is used for this experiment. In the case of $R=1$, the true GP has a local maximum around $-5$, while the quantized GP has a local minimum at that region. This phenomenon is observed again around $5$ where true GP has local minimum while the quantized GP has local maximum. Thus, the divergence of posterior distribution of true and quantized GPs for $R=1$ is very high. Also as can be seen from Figure 4a the standard deviation of quantized GP is greater than true GP in the region of data. This experiment shows that GP training using very low quality train data may not be favorable.

In despite of the case $R=1$, the quantized GP for $R=2$ is more close to the true GP. Here we have no reverse peaks between the true and quantized GPs. However, it seems that the divergence between the two GPs is relatively high. For rate $R=3$, the mean function is fairly close to the true version, but the standard deviation in the region of training inputs is slightly larger than the true one. For rates $R\geq 6$, both mean and standard deviation of quantized GP is very close to the true GP. This experiment shows that by consuming a few bits for each data sample, we can reach to the result of full GP model, but if the bit rate is extremely low, say $1$ bit/sample in this experiment, the result may be unappealing.

In order to evaluate the trained GPs on quantized data we have used three popular real datasets: the SARCOS ($44,484$ training, $4,449$ test, 21 attributes), the KIN40K ($10,000$ training, $30,000$ test, 8 attributes), and the ABOLONE ($3,133$ training, $1,044$ test, 8 attributes). Learning several regression models on these datasets has been studied in previous works such as [27]. Since the whole datasets are too large to be trained, 1000 samples selected at random are used for the training GPs. To make a reference, a subset of unquantized data (SD) is used to train GPs at a center machine. We refer GPs trained using SD method as full GP. The training dataset is randomly distributed across $40$ machines. Inputs are normalized to have zero mean and unit variance on the training set. Outputs (target values) are also centered around zero. We compare our method with previous distributed GP models that are based on PoE [25] such as BCM [26] and rBCM [24].

Figure 5 shows the error of regression for SARCOS and ABOLONE datasets. We have used the linear kernel function in (4) for our GP model. The regression error is calculated using the standard mean squared error (SMSE) given by $\frac{1}{N}\frac{\|\mathbf{y}_{*}-\hat{\mathbf{y}}_{*}\|^{2}}{var(\mathbf{y}_{*})}$, where $\mathbf{y}_{*}$ and $\hat{\mathbf{y}}_{*}$ are the target and prediction vectors, respectively. As can be seen from Figure 5a for SARCOS dataset the broadcast model crosses the rBCM error line at $R=16$. Since the input dimension of SARCOS is $21$, by consuming about $0.8$ bit for each dimension, the broadcast model outperforms rBCM. The broadcast model is near the full GP around $40$ bits per sample or about 2 bits per dimension.

Figure 5b shows similar results on ABOLONE dataset. The broadcast model outperforms rBCM around $R=16$, which leads to $2$ bits for each dimension. The broadcast model converges to the full GP around $50$ bits (i.e. about $6$ bits per dimension).

The GP model with linear kernel function in (4) does not converge for KIN40K dataset even for the full GP, thus, we could not plot its results.

Figure 6 shows the result of a similar experiment on the three datasets for GP models with squared exponential kernel function:

[TABLE]

where $\sigma_{s}$ and $\ell$ are the hyper-parameters of the model. Figure 6a shows that our methods outperforms rBCM after spending about $25$ bits per sample on SARCOS dataset in the broadcast GP. As mentioned before, the input dimension of SARCOS dataset is $21$. Hence, each dimension is encode on the average by about $1$ bit. For single-center GP, the error curve crosses the rBCM error line when about $4$ bits are used per dimension.

Figure 6b shows the result of the same experiment for the KIN40K dataset. Here, for bit rate around $20$ we outperform both BCM and rBCM in the broadcast GP. The single-center GP has slightly more error values than the broadcast GP. The input space of the KIN40K dataset is $8$-dimensional, thus, our broadcast GP model requires on the average about $2.5$ bits per dimension to outperform both BCM and rBCM. The larger average bit rate shows that the attributes of KIN40K are more informative than that of SARCOS. The error curves of ABOLONE dataset are shown in Figure 6c. For the three datasets we see that the zero rate error for the broadcast model is smaller than the error at rate $R=5$. This shows that the highly distorted data could decrease the performance of the model.

The preceding experiments reveal that in case of highly limited communication capacity our scheme performs worse than the rBCM method where none of data points is transferred between the machines. This suggests that when the communication is highly limiting, it is not a good strategy to transmit large number of data points with low qualities. It is better to transmit lower number of data samples with an acceptable quality.

To overcome this issue, one idea is to use our quantization method on a subset of samples. To select effective representatives of samples, we have used the method of [27] in which inducing (variational) variables are found. We, then, quantize the inducing variables and transmit them to the center machine. In fact, each machine runs locally the method of Titsias [27] to find inducing variables, then applies our proposed quantization methods on them. Figure 7 shows the error of single-center GP model which adopts this method for distributed GP learning on KIN40K dataset. As can be seen from the figure, the results are significantly better than non-sparse models in Figure 6.

It is also possible to find a global inducing set by communication between machines using method of [29] which is based on an iterative algorithm. The communication cost of [29] depends on the size of inducing set as well as the rate of convergence and is relatively high.

7 Conclusion

In this paper, we have proposed three methods for computing inner-products between Gaussian vectors with minimum distortion in distributed systems. We have shown that the well-known reverse water-filling algorithm on the product of local covariance matrices is the optimal coding scheme for Gaussian variables and it is a lower bound for other distributions. We have also proposed a simple per-symbol coding scheme which operates near the optimal one. The third proposed method is based on a dimension reduction scheme that is similar to the PCA except its objective function is the inner-product distortion. It is also proved that the optimal basis vectors for dimension reduction is obtained from product of local covariance matrices. This is similar to the result of the reverse water-filling scheme where the product of local matrices play the key role in the procedure of encoding.

Through several experiments, we have observed that estimation error in GP regressions drops significantly by spending a few bits per each symbol resulting in near the full GP. We also experimentally show that by combining the sparse GP models and the proposed quantization methods, we can obtain a model outperforming PoE-based models even for very low bit rates. This suggests that optimizing inducing points in sparse GPs alongside their quantization will result in a better performance in distributed settings.