The Mondrian Kernel

The Mondrian Kernel

Matej Balog
Department of Engineering
University of Cambridge &Balaji Lakshminarayanan  
Gatsby Unit
University College London \AND  Zoubin Ghahramani   
Department of Engineering
University of Cambridge &Daniel M. Roy
Department of Statistical Sciences
University of Toronto &Yee Whye Teh
Department of Statistics
University of Oxford
 Also affiliated with Max-Planck Institute for Intelligent Systems, Tübingen, Germany.
Abstract

We introduce the Mondrian kernel, a fast random feature approximation to the Laplace kernel. It is suitable for both batch and online learning, and admits a fast kernel-width-selection procedure as the random features can be re-used efficiently for all kernel widths. The features are constructed by sampling trees via a Mondrian process [Roy and Teh, 2009], and we highlight the connection to Mondrian forests [Lakshminarayanan et al., 2014], where trees are also sampled via a Mondrian process, but fit independently. This link provides a new insight into the relationship between kernel methods and random forests.

 


The Mondrian Kernel


 


Matej Balogthanks:  Also affiliated with Max-Planck Institute for Intelligent Systems, Tübingen, Germany. Department of Engineering University of Cambridge                        Balaji Lakshminarayanan Gatsby Unit University College London


  Zoubin Ghahramani Department of Engineering University of Cambridge                        Daniel M. Roy Department of Statistical Sciences University of Toronto                        Yee Whye Teh Department of Statistics University of Oxford


1 Introduction

Kernel methods such as support vector machines and Gaussian processes are very popular in machine learning. While early work relied on dual optimization, recent large-scale kernel methods focus on the primal optimization problem where the input data are mapped to a finite-dimensional feature space and the weights are learned using fast linear optimization techniques, e.g., stochastic gradient descent. Rahimi and Recht [2007] proposed to approximate shift-invariant kernels by mapping the inputs to so-called random features, constructed so that the inner product of two mapped data points approximates the kernel evaluated at those two points (which is the inner product in the feature space corresponding to the kernel). Rahimi and Recht [2007] proposed two random feature construction schemes: random Fourier features, where data points are projected onto random vectors drawn from the Fourier transform of the kernel and then passed through suitable non-linearities; and random binning, where the input space is partitioned by a random regular grid into bins and data points are mapped to indicator vectors identifying which bins they end up in. Both of these approaches require specifying the kernel hyperparameters in advance, so that the appropriate distribution is used for sampling the random vectors or random grids, respectively. However, a suitable kernel width (length-scale) is often not known a priori and is found by cross-validation, or, where available, marginal likelihood optimization. In practice, this entails constructing a new feature space and training a linear learner from scratch for each kernel width, which is computationally expensive. Using a suitable kernel width is often more important than the choice of kernel type [Schölkopf and Smola, 2001], so a fast kernel width selection method is desirable.

We describe a connection between the Laplace kernel and the Mondrian process [Roy and Teh, 2009], and leverage it to develop a random feature approximation to the Laplace kernel that addresses the kernel width selection problem. This approximation, which we call the Mondrian kernel, involves random partitioning of data points using a Mondrian process, which can be efficiently reused for all kernel widths. The method preserves the nonparametric nature of kernel learning and is also suitable for online learning.

The Mondrian kernel reveals an interesting link between kernel methods and decision forests [Breiman, 2001, Criminisi et al., 2012], another popular class of nonparametric methods for black-box prediction tasks. The Mondrian kernel resembles Mondrian forests, a decision-forest variant introduced by Lakshminarayanan et al. [2014], where a Mondrian process is used as the randomization mechanism. The efficiently trainable Mondrian forests excel in the online setting, where their distribution is identical to the corresponding batch Mondrian forest, and have been successfully applied to both classification and regression [Lakshminarayanan et al., 2014, 2016]. Mondrian forests and the Mondrian kernel both lead to randomized, non-linear learning algorithms whose randomness stems from a Mondrian process. The former fits parameters corresponding to different Mondrian trees independently, while the latter fits them jointly. We compare these methods theoretically and thus establish a novel connection between the Laplace kernel and Mondrian forests via the Mondrian kernel.

The contributions of this paper are:

  • a review of the Mondrian process using the simple notion of competing exponential clocks (Section 2);

  • a novel connection between the Mondrian process and the Laplace kernel (Section 3), yielding a fast approximation to learning with the Laplace kernel;

  • an efficient procedure for learning the kernel width from data (Section 4); and

  • a comparison between Mondrian kernel and Mondrian forest that provides another connection between kernel learning and random forests (Section 6).

2 Mondrian Process

For completeness, we review the Mondrian process [Roy and Teh, 2009, Roy, 2011, Chapter 5]. Although simple and perhaps well known to experts, our exposition through competing exponential clocks has not explicitly appeared in this form in the literature. Readers familiar with the Mondrian process may skip this section on first reading.

2.1 Terminology

An axis-aligned box is a Cartesian product of bounded intervals . Their total length is the linear dimension of . A guillotine partition of is a hierarchical partitioning of using axis-aligned cuts. Such a partition can be naturally represented using a strictly binary tree.

An exponential clock with rate takes a random time to ring after being started, where is the exponential distribution with rate (inverse mean) . The notion of competing exponential clocks refers to independent exponential clocks with rates , started at the same time. It can be shown that (1) the time until some clock rings has distribution, (2) it is the -th clock with probability proportional to , and (3) once a clock rings, the remaining clocks continue to run independently with their original distributions.

2.2 Generative Process

The Mondrian process on an axis-aligned box is a time-indexed stochastic process taking values in guillotine-partitions of . It starts at time with the trivial partition of (no cuts) and as time progresses, new axis-aligned cuts randomly appear, hierarchically splitting into more and more refined partitions. The process can be stopped at a lifetime , which amounts to ignoring any cuts that would appear after time .

To describe the distribution of times and locations of new cuts as time progresses, we associate an independent exponential clock with rate to each dimension of . Let be the first time when a clock rings and let be the dimension of that clock. If then this process terminates. Otherwise, a point is chosen uniformly at random from and is split into and by a hyperplane in dimension that is perpendicular to at point . After making this first cut, the remaining clocks are discarded and the generative process restarts recursively and independently on and . However, those processes start at time rather than and thus have less time left until the lifetime is reached.

The specification of the generative process on is now complete. Due to the properties of competing exponential clocks, the time until the first cut appears in has exponential distribution with rate equal to the linear dimension of and the dimension in which the cut is made is chosen proportional to . This confirms equivalence of our generative process to the one proposed by Roy and Teh [2009]. Finally, we note that a.s. the Mondrian process does not explode, i.e., for every lifetime , the process generates finitely many cuts with probability [Roy, 2011].

\includegraphics

[width=]figures/mondrian-sample-1x1

(a)
\includegraphics

[width=]figures/mondrian-sample-1x1-tree

(b)
Figure 1: (a) Sample of a Mondrian process on the axis-aligned box with lifetime . Numbers on the cuts (shown in green) indicate the times when they appeared. The first cut appeared at time , in dimension , at location . (b) Representing the Mondrian sample as a strictly binary tree, with new nodes (shown as circles) appearing as time (y-axis) progresses. The two numbers below each node show the rates of the two exponential clocks competing to split that node, with the winning clock’s rate shown in green.

2.3 Projectivity

If a Mondrian process runs on , what distribution of random partitions does it induce on an axis-aligned subbox ? (See Figure 1(a) for an illustration in dimensions.) The Mondrian process was constructed so that the answer is the Mondrian process itself [Roy, 2011]. Here we explain this projectivity property using the notion of competing exponential clocks. To argue that the resulting process on is indeed a Mondrian process, we show that the process running on generates cuts in in the same way as a Mondrian process running directly on would.

Recall that each dimension of is associated with an exponential clock with rate and if it rings first, the cut location is chosen uniformly at random from . This procedure can be equivalently represented using two competing clocks for each dimension (rather than just one):

  • Clock with rate . If this clock rings first, the cut location is chosen uniformly at random from .

  • Clock with rate . If it rings first, the cut location is sampled uniformly from .

(See Figure 1(b).) Note that the clocks represent the same cut distribution as a Mondrian process running on would. If a clock rings first, a cut is made outside of and all of remains on one side of this cut. None of the clocks have rung in that case and would usually be discarded and replaced with fresh clocks of identical rates, but by property (3) of competing exponential clocks, we can equivalently reuse these clocks (let them run) on the side of the cut containing . (Figure 1(c) shows a cut in dimension that misses and the reused clocks ). Hence, cuts outside do not affect the distribution of the first cut crossing , and this distribution is the same as if a Mondrian process were running just on . When a cut is made within (see Figure 1(d)), the process continues on both sides recursively and our argument proceeds inductively, confirming that the Mondrian process on generates cuts in in the same way as a Mondrian process on would.

\includegraphics

[scale=0.74]figures/consistency-fig1

(a)
\includegraphics

[scale=0.74]figures/consistency-fig2

(b)
\includegraphics

[scale=0.74]figures/consistency-fig3

(c)
\includegraphics

[scale=0.74]figures/consistency-fig4

(d)
Figure 2: (a) A Mondrian process running on generates cuts (dashed lines), some of which intersect (green lines) and thus induces a random partition of . (b) Representing the first cut distribution using competing exponential clocks: in each dimension , clock corresponds to the region where making a cut splits (shown in green) and clock to the (disconnected) region where making a cut misses (shown in red). (c) Cut outside : reusing the clocks , on the side of the cut containing . (d) Cut inside (shown in black): the argument proceeds by induction on both sides.

2.4 Mondrian Process on

The Mondrian process on is defined implicitly as a time-indexed stochastic process such that its restriction to any axis-aligned box is a Mondrian process as defined in section 2.2. Fortunately, this infinite-dimensional object can be compactly represented by instantiating the Mondrian process only in regions where we have observed data. As we observe new data points, the Mondrian sample can be extended using the conditional Mondrian algorithm [Roy and Teh, 2009], a simple and fast sampling procedure for extending a Mondrian sample in an axis-aligned box to a larger axis-aligned box . The conditional Mondrian is useful for online learning and prediction, as it can be used to extend Mondrian samples to (yet) unobserved parts of the input space [Lakshminarayanan et al., 2014].

3 Mondrian Kernel

For concreteness, our running example will be regression: the problem of learning a function from a set of training examples . However, the Mondrian kernel applies equally well to classification, or any other learning task.

Learning with kernels involves choosing a kernel function to act as a similarity measure between input data points. Evaluating on all pairs of data points takes operations, with some models also requiring a operation on an kernel matrix. This generally makes exact kernel methods unsuitable for large-scale learning. Rahimi and Recht [2007] proposed a fast approximation through a randomized construction of a low-dimensional feature map such that

and then using a linear learning method in the feature space implied by . For example, linear regression , where is the feature matrix with -th row , is solvable exactly in time linear in . In general, the primal problem also lends itself naturally to stochastic gradient descent approaches for learning .

We use the Mondrian process to construct a randomized feature map for the (isotropic) Laplace kernel:

Here is the inverse kernel width (inverse length-scale), which we call the lifetime parameter of the kernel. We use a non-standard parametrization as this lifetime parameter will be linked to the Mondrian process lifetime.

3.1 Mondrian Kernel

Consider the following randomized construction of a feature map :

  1. Sample a partition of via a Mondrian process on with lifetime . Label the cells of the generated partition by in arbitrary order.

  2. To encode a data point , look up the label of the partion cell falls into and set to be the (column) indicator vector that has a single non-zero entry at position , equal to .

The Mondrian process on generates infinitely many partition cells and cannot be stored in memory, but projectivity comes to the rescue. As we only ever need to evaluate on finitely many data points, it suffices to run the Mondrian on the smallest axis-aligned box containing all these points. Also, we only label partition cells containing at least one data point, in effect removing features that would be for all our data points. Then, the dimensionality of equals the number of non-empty partiton cells and each data point has a single non-zero feature, equal to .

\includegraphics [scale=0.7]figures/mondrian-features [0 1 0] [0 1 0] [0 0 1] [1 0 0]
Figure 3: Feature expansions of data points in .

However, note that the set of points on which the feature map is evaluated need not be known in advance and can even grow in an online fashion. Indeed, the conditional Mondrian algorithm discussed in section 2.4 allows us to extend Mondrian samples to larger boxes as necessary, and we can increase the dimensionality of whenever a data point is added to a previously empty partition cell.

This feature map induces a kernel

(1)

which we call a Mondrian kernel of order .

Instead of using a single Mondrian sample (partition), we can use independent samples and construct a feature map by concatenating and normalizing the feature maps obtained from each individual sample as above:

(2)

This feature expansion is sparse: every data point has exactly non-zero features. The corresponding kernel, which we call a Mondrian kernel of order , is

This is the empirical frequency with which points and end up in the same partition cell of a Mondrian sample.

1:for  to  do
2:     construct feature map section 3.1
3:join and rescale into equation (2)
4:map data to feature representations using
5:use linear learning method on
Algorithm 1 Mondrian kernel

3.2 Mondrian–laplace Link

By independence of the Mondrian samples, a.s.

with convergence at the standard rate . We thus define the Mondrian kernel of order as

Proposition 1 (Mondrian-Laplace link).

The Mondrian kernel of order coincides with the Laplace kernel.

Proof.

As (defined in (1)) is a binary random variable, equals the probability that and fall into the same partition cell of a Mondrian sample, which is equivalent to the sample having no cut in the minimal axis-aligned box spanned by and .

\includegraphics

[width=0.75]figures/mondrian-laplace

By projectivity, this probability is the same as the probability of not observing any cuts in a Mondrian process with lifetime running on just this minimal box. Noting that the linear dimension of this box is , we obtain

Note that the lifetime (inverse width) of the Laplace kernel corresponds to the lifetime of the Mondrian process used in the construction of the Mondrian kernel.

This link allows us to approximate the Laplace kernel with a Mondrian kernel , which, unlike the Laplace kernel, admits a finite-dimensional feature expansion. The finite order trades off kernel approximation error and computational costs (indirectly through the complexity of ).

The following result confirms that the convergence of the Mondrian kernel approximation is exponentially fast in uniformly on any fixed bounded input domain .

Proposition 2.

For any bounded input domain and , as ,

Proof.

Given in Supplement A. ∎

4 Fast Kernel Width Learning

This section discusses the main advantage of our Mondrian approximation to the Laplace kernel: the efficient learning of kernel width from data. In particular, the approximation allows for efficient evaluation of all kernel lifetimes (inverse widths) , where the terminal lifetime need not be fixed a priori.

4.1 Feature Space Reusal

We make the following recollections from earlier sections:

  • the Mondrian process runs through time, starting at time and only refining the generated partition as time progresses (cuts are never removed)

  • the Mondrian process with lifetime is obtained by ignoring any cuts that would occur after time

  • the lifetime of the Mondrian process used in constructing an explicit feature map for a Mondrian kernel corresponds to the lifetime (inverse width) of the Laplace kernel that it approximates

Running the Mondrian process from time to some terminal lifetime thus sweeps through feature spaces approximating all Laplace kernels with lifetimes . More concretely, we start with and the feature map corresponding to trivial partitions, i.e., for any data point , the vector has length and all entries set to the normalizer . As we increase , at discrete time points new cuts appear in the Mondrian samples used in constructing . Suppose that at some time , the partition cell corresponding to the -th feature in is split into two by a new cut that first appeared at this time . We update the feature map by removing the -th feature and appending two new features, one for each partition cell created by the split. See Figure 4 for an example with .

\includegraphics [scale=0.7]figures/mondrian-features-update \includegraphics [scale=0.9]figures/mondrian-features-lifetime-update-table
Figure 4: A new cut (shown in thick blue) appeared, splitting cell (cf. Figure 3) into two new cells and . The table shows the update to , with the removed feature in gray italics and the two new features in bold blue.

This procedure allows us to approximate all Laplace kernels with lifetimes without having to resample new feature spaces for each lifetime. The total computational cost is the same (up to a multiplicative constant) as of constructing a single feature space just for the terminal lifetime . This is because a strictly binary tree with leaves (partition cells in the -th Mondrian sample at time ) contains at most internal nodes (features that had to be removed at some time point ).

4.2 Linear Learner Retraining

Evaluating suitability of a lifetime (inverse kernel width) requires training and evaluating a linear model in the feature space implied by . This can also be done more efficiently than retraining a new model from scratch every time a new cut is added and updated. We discuss the example of ridge regression with exact solutions, and a general case of models trainable using gradient descent methods.

4.2.1 Ridge regression

The MAP weights of the primal ridge regression problem are , where is the regularized feature covariance matrix and is the regularization hyperparameter. Instead of inverting , it is numerically more stable to work with its Cholesky factor [Seeger, 2003]. Phrasing the problem as Bayesian linear regression with, say, observation noise variance and prior weights variance , we can also obtain the log marginal likelihood of the form

where the dependence on is implicit through .

When a new cut appears in one of the Mondrian samples and is updated by deleting the -th feature and appending two new ones, the corresponding update to the regularized feature covariance matrix is to delete its -th row and -th column, and append two new rows and columns. Then both and can be appropriately updated in time, faster than recomputation from scratch. Updating the Cholesky factor when the -th row and column are removed is slightly involved but can be achieved by first permuting the rows and columns so that the ones to be removed are the last ones [Seeger, 2004], after which the Cholesky factor is updated by deleting its last row and column. If is the number of features at the terminal lifetime , this update is performed times, for a total computational cost . Note that performing the inversion or Cholesky factorization at just the terminal lifetime would have the same time complexity.

After updating or , the optimal weights can be updated in time and the determinant of required for the marginal likelihood can be obtained from as the squared product of its diagonal elements in time. Exploiting sparsity of , evaluating the model on data points takes time.

Finally, we note that computing the marginal likelihood for all and combining it with a prior supported on allows Bayesian inference over the kernel width . We refer to Supplement B for more details.

4.2.2 Models trainable using gradient descent

Consider a linear model trained using a gradient descent method. If (an approximation to) the optimal weight vector is available and then is updated by removing the -th feature and appending two new features, a natural way of reinitializing the weights for subsequent gradient descent iterations is to remove the -th entry of and append two new entries, both set to the removed value (as points in the split cell are partitioned into the two new cells, this preserves all model predictions). Note that we have the freedom of choosing the number of gradient descent iterations after each cut is added, and we can opt to only evaluate the model (on a validation set, say) at several values on the first pass through . One iteration of stochastic gradient descent takes time thanks to sparsity of .

This efficient kernel width selection procedure can be especially useful with models where hyperparameters cannot be tweaked by marginal likelihood optimization (e.g., SVM).

5 Online Learning

In this section, we describe how the Mondrian kernel can be used for online learning. When a new data point arrives, incorporating it into existing Mondrian samples (using the conditional Mondrian algorithm discussed in section 2.4) can create new non-empty partition cells, increasing the dimensionality of the feature map . We set the new features to for all previous data points .

In our running example of ridge regression, exact primal updates can again be carried out efficiently. The inverse or Cholesky factor of the regularized feature covariance matrix can be updated in two steps:

  1. extend or to incorporate the new features (set to for all existing data points) in

  2. incorporate the new data point , which is now a simple rank-1 update on , so or can again be updated efficiently in time

We refer to Supplement C for more details.

With gradient descent trainable models, we maintain (an approximation to) the optimal weights directly. When a new data point arrives, we expand the dimensionality of as described above. The previously optimal weights can be padded with ’s in any newly added dimensions, and then passed to the gradient descent method as initialization.

6 Link to Mondrian Forest

We contrast Mondrian kernel with Mondrian forest [Lakshminarayanan et al., 2014, 2016], another non-linear learning method based on the Mondrian process. They both start by sampling independent Mondrians on to provide independent partitions of the data. However, these partitions are then used differently in the two models:

  • In a Mondrian forest, parameters of predictive distributions in each tree are fitted independently of all other trees. The prediction of the forest is the average prediction among the trees.

  • With Mondrian kernel, the weights of all random features are fitted jointly by a linear learning method.

Let count the leaves (non-empty partition cells) in the -th Mondrian sample and let be the total number of leaves. Let be the indicator of the partition cell in the -th sample into which the -th data point falls (as in section 3.1). Also, as in equation (2), let be the normalized concatenated feature encoding of the -th data point. Recall that each vector contains exactly non-zero entries, all of which equal the normalizer .

For simplicity, we restrict our attention to ridge regression in this section and compare the learning objective functions of Mondrian kernel and Mondrian forest.

6.1 Mondrian Kernel Objective

The primal ridge regression problem in the feature space implied by is

Decomposing , so that each (rescaled) subvector corresponds to features from the -th Mondrian, denoting by the “contribution” of the -th Mondrian to the prediction at the -th data point, and writing , the Mondrian kernel objective function can be restated as

(3)

6.2 Mondrian Forest Objective

Assuming a factorizing Gaussian prior over the leaves in each Mondrian tree (i.e., without the hierarchical smoothing used by Lakshminarayanan et al. [2016]), the predictive mean parameters in the leaves of the -th Mondrian tree are fitted by minimizing

where is the ratio of noise and prior variance in the predictive model. The parameters are disjoint for different trees, so these independent optimization problems are equivalent to minimizing the average of the individual objectives. Writing for the -th tree’s prediction at the -th data point and concatenating the parameters , the Mondrian forest objective can be stated as

(4)

6.3 Discussion

Comparing (3) and (4), we see that subject to regularization parameters (priors) chosen compatibly, the two objectives only differ in the contribution of an individual data point to the total loss:

Mondrian kernel:
Mondrian forest:

Specifically, the difference is in the order in which the averaging over Mondrian samples/trees and the non-linear loss function are applied. In both models predictions are given by , so the Mondrian kernel objective is consistent with the aim of minimizing empirical loss on the training data, while the forest objective minimizes average loss across trees, not the loss of the actual prediction (when ) [Ren et al., 2015].

Ren et al. [2015] address this inconsistency between learning and prediction by proposing to extend random forests with a global refinement step that optimizes all tree parameters jointly, minimizing the empirical training loss. Our approximation of the Laplace kernel via the Mondrian kernel can be interpreted as implementing this joint parameter fitting step on top of Mondrian forest, revealing a new connection between random forests and kernel methods.

\includegraphics

[width=0.95]figures/forest-kernel

7 Related Work

The idea of Rahimi and Recht [2007] to approximate shift-invariant kernels by constructing random features has been further developed by Le et al. [2013] and Yang et al. [2015], providing a faster method of constructing the random features when the input dimension is high. The fast method of Dai et al. [2014] can adapt the number of random features, making it better-suited for streaming data. To the best of our knowledge, these methods require random features to be reconstructed from scratch for each new kernel width value; however, our solution allows us to efficiently learn this hyperparameter for the Laplace kernel.

Decision forests are popular for black-box classification and regression thanks to their competitive accuracy and computational efficiency. The most popular variants are Breiman’s Random Forest [Breiman, 2001] and Extremely Randomized Trees [Geurts et al., 2006]. Breiman [2000] established a link between the Laplace kernel and random forests with an infinite number of trees, but unlike our work, made two additional strong assumptions, namely infinite data and a uniform distribution of features. From a computational perspective, Shen et al. [2006] approximated evaluation of an isotropic kernel using -trees, reducing computational complexity as well as memory requirements. Davies and Ghahramani [2014] constructed ‘supervised’ kernels using random forests and demonstrated that this can lead to linear-time inference. We refer to [Scornet, 2015] for a recent discussion on the connection between decision forests and kernel methods.

A key difference between decision forests and kernel methods is whether parameters are fit independently or jointly. In decision forests, the leaf node parameters for each tree are fit independently, whereas the weights of random features are fit jointly. Scornet [2015] shows that random forests can be interpreted as adaptive kernel estimates and discusses the theoretical properties of fitting parameters jointly. Ren et al. [2015] propose to extend random forests with a global refinement step, optimizing all tree parameters jointly to minimize empirical training loss.

The proposed Mondrian kernel establishes a link between Mondrian trees and Laplace kernel for finite data, without any assumptions on the distribution of the features. Unlike prior work, we exploit this connection to construct an adaptive random feature approximation and efficiently learn the kernel width.

8 Experiments

We conducted three sets of experiments, with these goals:

  1. verify that Mondrian kernel approximates the Laplace kernel, and compare to other random feature generation schemes (Section 8.1);

  2. demonstrate usefulness of our efficient kernel width selection procedure, showing that it can quickly learn a suitable kernel width from data (Section 8.2); and

  3. empirically compare the Mondrian kernel and Mondrian forests, supporting the insight into their relationship from Section 6 (Section 8.3).

With the exception of two experiments on synthetic data, we carried out our evaluation on the CPU dataset from [Rahimi and Recht, 2007], containing training and test points with attributes. Note that the CPU dataset is an adversarial choice here, as Rahimi and Recht [2007] report that random Fourier features perform better than binning schemes on this task. In all experiments, the ridge regularization constant was set to , the value used by Rahimi and Recht [2007], and the primal optimization problems were solved using stochastic gradient descent.

8.1 Laplace Kernel Approximation

First we examined the absolute kernel approximation error directly. To this end, we sampled data points uniformly at random in the unit square and computed the maximum absolute error over all pairs of points. The Laplace kernel and Mondrian kernels had a common lifetime (inverse width) , so that several widths fit into the input domain . We repeated the experiment times for each value of , showing the results in Figure 5. We plot the maximum error against the number of non-zero features per data point, which is relevant for solvers such as Pegasos SVM [Shalev-Shwartz et al., 2007], whose running time scales with the number of non-zero features per data point. Under this metric, the Mondrian kernel and Random binning converged to the Laplace kernel faster than random Fourier features, showing that in some cases they can be a useful option. (The error of Random Fourier features would decrease faster when measured against the total number of features, as Mondrian kernel and Random binning generate sparse feature expansions.)

\includegraphics

[scale=0.5]figures/2016-02-24_12-58-44_random_features_Mmax_50_N_100.pdf

Figure 5: Maximum absolute kernel approximation error on all pairs of data points in .

Second, we examined the approximation error indirectly via test set error on the CPU dataset. We repeated the experiment times for each value of and show the results in Figure 6. Even though Fourier features are better suited to this task, for a fast approximation with few () non-zero features per data point, random binning and Mondrian kernel are still able to outperform the Fourier features.

\includegraphics

[scale=0.5]figures/2016-02-27_17-06-14_convergence_RMSE_Mmax_50_CPU.pdf

Figure 6: Test set error on the CPU dataset. The horizontal line at indicates the error achieved with an exact, but expensive computation using the Laplace kernel.

8.2 Fast Kernel Width Learning

First, using a synthetic regression dataset generated from a Laplace kernel with known ground truth lifetime , we verified that the lifetime could be recovered using our kernel width selection procedure from Section 4. To this end, we let the procedure run until a terminal lifetime and plotted the error on a held-out validation set as a function of the lifetime . The result in Figure 7 shows that the ground truth kernel lifetime was recovered within an order of magnitude by selecting the lifetime minimizing validation set error. Moreover, this value of led to excellent performance on an independent test set.

\includegraphics

[scale=0.45]figures/2016-03-01_10-40-18_width_recovery_D_2_M_50.pdf

Figure 7: Recovering the ground truth lifetime by selecting the value minimizing validation set error.

Second, we evaluated our kernel width selection procedure on the CPU dataset in order to demonstrate its practical usefulness. While the Mondrian kernel allows to efficiently sweep through lifetimes , Fourier features and random binning need to be reconstructed and retrained for each attempted lifetime value. We started the Fourier features and random binning at , and in each step, we either doubled the maximum lifetime or halved the minimum lifetime considered so far, based on which direction seemed more promising. Once a good performing lifetime was found, we further optimized using a binary search procedure. All schemes were set to generate non-zero features per datapoint. Figure 8 shows the performance of each scheme on a held-out validation set as a function of computation time. The result suggests that our kernel width learning procedure can be used to discover suitable lifetimes (inverse kernel widths) at least an order of magnitude faster than random Fourier features or random binning.

8.3 Mondrian Kernel vs Forest

\includegraphics

[width=]figures/2016-02-29_10-20-11_performance_vs_time_CPU_M_350.pdf

Figure 8: Validation set error as a function of computation time. Even though Fourier features are better suited to the CPU dataset [Rahimi and Recht, 2007] and eventually outperform the Mondrian kernel, the latter discovers suitable kernel widths at least an order of magnitude faster.

We compared the performance of Mondrian kernel and “Mondrian forest” (quotes due to omission of hierarchical smoothing) based on the same Mondrian samples, using the CPU dataset and varying the lifetime . Recall that higher values of lead to more refined Mondrian partitions, allowing more structure in the data to be modeled, but also increasing the risk of overfitting. Figure 9 shows that Mondrian kernel exploits the joint fitting of parameters corresponding to different trees and achieves a lower test error at lower lifetime values, thus producing a more compact solution based on simpler partitions. Figure 10 shows the parameter values learned by Mondrian kernel and Mondrian forest at the lifetime . The distribution of weights learned by Mondrian kernel is more peaked around , as the joint fiting allows achieving more extreme predictions by adding together several smaller weights.

\includegraphics

[scale=0.464]figures/2016-05-23_19-23-27_kernel_vs_forest_CPU_M_50.pdf

Figure 9: Comparison of Mondrian kernel and Mondrian forest models based on the same set of Mondrian samples.
\includegraphics

[width=]figures/2016-05-23_19-27-16_kernel_vs_forest_weights_CPU_M_50.pdf

Figure 10: Weights learned by Mondrian forest and Mondrian kernel at the lifetime in Figure 9.

9 Conclusion

We presented the Mondrian kernel, a fast approximation to the Laplace kernel that admits efficient kernel width selection. When a different kernel or a different approximation is used, our procedure can provide a fast and simple way of initializing the kernel width for further optimization. While a Gaussian kernel is often considered a default choice, in many situations it imposes an inappropriately strong smoothness assumption on the modelled function and the Laplace kernel may in fact be a preferable option.

Our approach revealed a novel link between the Mondrian process and the Laplace kernel. We leave the discovery of similar links involving other kernels for future work.

Acknowledgements

We would like to thank Nilesh Tripuraneni for useful discussions. Part of this research was carried out while MB was at the University of Oxford. BL gratefully acknowledges generous funding from the Gatsby Charitable Foundation. ZG acknowledges funding from the Alan Turing Institute, Google, Microsoft Research and EPSRC Grant EP/N014162/1. DMR is supported by an NSERC Discovery Grant. YWT’s research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) ERC grant agreement no. 617071.

References

References

  • Breiman [2000] L. Breiman. Some infinity theory for predictor ensembles. Technical report, University of California at Berkeley, 2000.
  • Breiman [2001] L. Breiman. Random forests. Mach. Learn., 45:5–32, 2001.
  • Criminisi et al. [2012] A. Criminisi, J. Shotton, and E. Konukoglu. Decision forests: A unified framework for classification, regression, density estimation, manifold learning and semi-supervised learning. Found. Trends Comput. Graphics and Vision, 2012.
  • Dai et al. [2014] B. Dai, B. Xie, N. He, Y. Liang, A. Raj, M.-F. F. Balcan, and L. Song. Scalable kernel methods via doubly stochastic gradients. In Adv. Neural Information Proc. Systems (NIPS), 2014.
  • Davies and Ghahramani [2014] A. Davies and Z. Ghahramani. The random forest kernel and other kernels for big data from random partitions. arXiv preprint arXiv:1402.4293v1, 2014.
  • Geurts et al. [2006] P. Geurts, D. Ernst, and L. Wehenkel. Extremely randomized trees. Mach. Learn., 63(1):3–42, 2006.
  • Lakshminarayanan et al. [2014] B. Lakshminarayanan, D. M. Roy, and Y. W. Teh. Mondrian forests: Efficient online random forests. In Adv. Neural Information Proc. Systems (NIPS), 2014.
  • Lakshminarayanan et al. [2016] B. Lakshminarayanan, D. M. Roy, and Y. W. Teh. Mondrian forests for large scale regression when uncertainty matters. In Int. Conf. Artificial Intelligence Stat. (AISTATS), 2016.
  • Le et al. [2013] Q. Le, T. Sarlós, and A. Smola. Fastfood-approximating kernel expansions in loglinear time. In Proc. Int. Conf. Mach. Learn. (ICML), 2013.
  • Rahimi and Recht [2007] A. Rahimi and B. Recht. Random features for large-scale kernel machines. In Adv. Neural Information Proc. Systems (NIPS), 2007.
  • Ren et al. [2015] S. Ren, X. Cao, Y. Wei, and J. Sun. Global refinement of random forest. In Proc. IEEE Conference on Computer Vision and Pattern Recognition, pages 723–730, 2015.
  • Roy [2011] D. M. Roy. Computability, inference and modeling in probabilistic programming. PhD thesis, Massachusetts Institute of Technology, 2011.
  • Roy and Teh [2009] D. M. Roy and Y. W. Teh. The Mondrian process. In Adv. Neural Information Proc. Systems (NIPS), 2009.
  • Schölkopf and Smola [2001] B. Schölkopf and A. J. Smola. Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond. MIT Press, Cambridge, MA, USA, 2001. ISBN 978-0-262-19475-4.
  • Scornet [2015] E. Scornet. Random forests and kernel methods. arXiv preprint arXiv:1502.03836v2, 2015.
  • Seeger [2003] M. Seeger. Bayesian Gaussian Process Models: PAC-Bayesian Generalisation Error Bounds and Sparse Approximations. PhD thesis, University of Edinburgh, 2003.
  • Seeger [2004] M. Seeger. Low rank updates for the Cholesky decomposition. Technical report, University of California at Berkeley, 2004.
  • Shalev-Shwartz et al. [2007] S. Shalev-Shwartz, Y. Singer, and N. Srebro. Pegasos: Primal Estimated sub-GrAdient SOlver for SVM. In Proc. Int. Conf. Mach. Learn. (ICML), 2007.
  • Shen et al. [2006] Y. Shen, A. Ng, and M. Seeger. Fast Gaussian process regression using KD-trees. In Adv. Neural Information Proc. Systems (NIPS), 2006.
  • Yang et al. [2015] Z. Yang, A. J. Smola, L. Song, and A. G. Wilson. A la Carte - Learning Fast Kernels. In Int. Conf. Artificial Intelligence Stat. (AISTATS), 2015.

The Mondrian Kernel

Supplementary material

Appendix A Proofs

Definition 1.

The linear dimension of an axis-aligned box is .

Our first result is a tail bound on the number of partition cells generated by a Mondrian process. We will use it as a Lemma in Proposition 4, but it also confirms that with probability 1, the Mondrian process does not explode (does not generate infinitely many partition cells in finite time).

Proposition 3.

Let be a Mondrian process on an axis-aligned box . For , let be the number of partition cells generated by until time . Then

In particular, the Mondrian process does not explode.

Proof.

At any time , by lack of memory of the exponential distribution, the residual time until a partition cell splits into two has distribution and is independent of all other cells by construction of the Mondrian process. As , this cell splitting process is dominated by a Yule process with birth rate . The number of individuals at time of a Yule process with birth rate has geometric distribution with mean and Markov’s inequality yields

as claimed. Hence for any . ∎

We define an -grid covering a (closed) interval as a set of points at most distance apart, including the boundary points, and with minimal possible cardinality:

Definition 2.

Let be an interval of length and let . Define . An -grid covering is a set of points in such that , and for all .

Note that such an -grid exists by our choice of , as we can take, e.g., for . The next lemma bounds the probability that two arrivals of a Poisson process running on a bounded interval occur between two consecutive points of an -grid covering that interval.

Lemma 1.

Consider a Poisson process with rate running on a bounded interval . Let be an -grid covering of . Then the probability that two or more arrivals of the process occur between any two consecutive points of is at most .

Proof.

As the distance between any two consecutive points of the -grid is at most by definition, the number of arrivals in a line segment between such two points is dominated by a Poisson random variable with mean . As there are such segments, the sought probability can be upper bounded using a union bound as

and using twice, we obtain as claimed

Definition 2 also set us up for defining the concept of an -grid on higher-dimensional axis-aligned boxes:

Definition 3.

Let be an axis-aligned box and let . An -grid covering is a cartesian product , where each is an -grid covering of in the sense of Definition 2.

Proposition 4.

For any bounded input domain and , as ,

Proof.

By extending if necessary, we may assume without loss of generality that is an axis-aligned box with linear dimension .

Recall that a Mondrian kernel of order corresponds to a random features obtained from independent Mondrians with lifetime . Let be an -grid covering , where will be specified later. The proof will upper bound the probability of the following three “bad” events:

  • := { any of the Mondrian samples contains more than partition cells }

  • := { the common refinement of the Mondrian partitions, disregarding any potential cuts after generating cells in one Mondrian, has a partition cell that does not contain an element of }

  • := { -approximation fails on , i.e., for some , , }

The constant will be specified (optimized) later. Note that implies that all partition cells in the common refinement of all Mondrian partitions contain a grid point from . Since is constant in each such cell, making small enough, smoothness of the Laplace kernel will ensure that if holds then -approximation holds throughout .

Proposition 3 and a union bound over the Mondrian samples give immediately that

Note that the -grid contains at most grid points. Hoeffding’s inequality and a union bound over all pairs of grid points gives for any :

To upper bound the probability of , note that at any time , in each partition cell generated so far by any of the Mondrian processes, an exponential clock is associated to each dimension of the cell, and if that clock rings, the cell is split at a random location by a hyperplane lying in dimension . Consider the point process obtained by projecting the cut points from all partition cells onto their respective coordinate axes. If each Mondrian process generates no more than than partition cells until its lifetime is exhausted, the cut points on the -th coordinate axis come from at most partition cells, each having width at most in dimension . Therefore this point process on the -th coordinate axis can be thought of as taking a suitable subset of points generated by a Poisson point process with intensity . Thus by Lemma 1, the probability that two cut points in dimension fall between two adjacent coordinates of the -grid is upper bounded by . Observe that if this does not happen in any of the dimensions then all partition cells in the common refinement must contain a grid point from . Hence, taking the union bound over all dimensions,

Thus the probability of a “bad” event occuring is at most

and minimizing over gives

If holds then each cell in the common refinement of the Mondrian partitions contains an element of the -grid , and the Laplace kernel of lifetime changes by at most when moving from any point in to the nearest grid point in its partition cell (in the common refinement). Therefore, as long as (i.e., ), the event implies that -approximation holds throughout . The upper bound on above is minimized for

which tends to as and so for large enough , we do have . For these large enough it then holds that

Proposition 5.

In a Mondrian regression forest with a factorizing Gaussian prior over leaf predictions, the learning objective function can be stated as

Proof.

The predictive mean parameters in the leaves of the -th tree are fitted by solving

where is the ratio of noise and prior variance in the predictive model. The parameters are disjoint for different trees, so these independent optimization problems are equivalent to minimizing the average

where is the -th tree’s prediction at data point . Rewriting in terms of the squared loss and the normalized concatenated weights , the learning objective function becomes

Appendix B Bayesian kernel width learning

Section 4.2.1 described how in a ridge regression setting, the marginal likelihood can be efficiently computed for all . With a prior over the lifetime (inverse kernel width) whose support is included in , the posterior distribution over is

with normalizing constant

where is the sequence of times when new cuts appeared in any of the Mondrian samples. The predictive distribution at a new test point is obtained by marginalizing out :

where the mixing coefficients

can be precomputed and cached for faster predictions. The integrals can be readily evaluated if we have access to the cumulative distribution function of our prior , which we assume.

Appendix C Online learning

Mirroring Section 4.2.1, we discuss the example of ridge regression where exact online updates can be carried out. Assume we have access to the regularized feature covariance matrix and its inverse or Cholesky decomposition before a new data point arrives, and we wish to update these efficiently.

If the dimensionality of increases by due to creating new non-empty partition cells, we first append rows and columns to , containing only, except on the main diagonal we put . Correspondingly, or are updated by appending rows and columns, with non-zero entries only on the main diagonal. (These entries would equal in and in ). This ensures the feature map now incorporates all necessary features.

Noting that the -entry of counts data points belonging to partition cells and at the same time (this can be non-zero only if , correspond to different Mondrian samples), normalized by , and that the -entry of the outer product is if the new data point falls into both cells and , and otherwise, we see that

is a rank- update. Therefore both and can be updated efficiently in time and the new MAP weights