The Mondrian Kernel
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 kernelwidthselection procedure as the random features can be reused 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 Balog^{†}^{†}thanks: Also affiliated with MaxPlanck 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 largescale kernel methods focus on the primal optimization problem where the input data are mapped to a finitedimensional feature space and the weights are learned using fast linear optimization techniques, e.g., stochastic gradient descent. Rahimi and Recht [2007] proposed to approximate shiftinvariant kernels by mapping the inputs to socalled 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 nonlinearities; 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 (lengthscale) is often not known a priori and is found by crossvalidation, 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 blackbox prediction tasks. The Mondrian kernel resembles Mondrian forests, a decisionforest 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, nonlinear 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 axisaligned 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 axisaligned 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 axisaligned box is a timeindexed stochastic process taking values in guillotinepartitions of . It starts at time with the trivial partition of (no cuts) and as time progresses, new axisaligned 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].
2.3 Projectivity
If a Mondrian process runs on , what distribution of random partitions does it induce on an axisaligned 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.
2.4 Mondrian Process on
The Mondrian process on is defined implicitly as a timeindexed stochastic process such that its restriction to any axisaligned box is a Mondrian process as defined in section 2.2. Fortunately, this infinitedimensional 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 axisaligned box to a larger axisaligned 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 largescale learning. Rahimi and Recht [2007] proposed a fast approximation through a randomized construction of a lowdimensional 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 lengthscale), which we call the lifetime parameter of the kernel. We use a nonstandard 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 :

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

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 nonzero 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 axisaligned 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 nonempty partiton cells and each data point has a single nonzero feature, equal to .
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 nonzero 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.
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 (MondrianLaplace 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 axisaligned box spanned by and .
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 finitedimensional 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 .
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 .
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 nonempty 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:

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

incorporate the new data point , which is now a simple rank1 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 nonlinear 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 (nonempty 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 nonzero 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 nonlinear 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.
7 Related Work
The idea of Rahimi and Recht [2007] to approximate shiftinvariant 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 bettersuited 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 blackbox 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 lineartime 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:
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 nonzero features per data point, which is relevant for solvers such as Pegasos SVM [ShalevShwartz et al., 2007], whose running time scales with the number of nonzero 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.)
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 () nonzero features per data point, random binning and Mondrian kernel are still able to outperform the Fourier features.
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 heldout 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.
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 nonzero features per datapoint. Figure 8 shows the performance of each scheme on a heldout 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
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.
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/20072013) 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 semisupervised 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. Fastfoodapproximating 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 largescale 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 9780262194754.
 Scornet [2015] E. Scornet. Random forests and kernel methods. arXiv preprint arXiv:1502.03836v2, 2015.
 Seeger [2003] M. Seeger. Bayesian Gaussian Process Models: PACBayesian 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.
 ShalevShwartz et al. [2007] S. ShalevShwartz, Y. Singer, and N. Srebro. Pegasos: Primal Estimated subGrAdient 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 KDtrees. 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 axisaligned 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 axisaligned 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 higherdimensional axisaligned boxes:
Definition 3.
Let be an axisaligned 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 axisaligned 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 nonempty 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 nonzero 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 nonzero 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