Kernel Bayes' rule

Introduction

Kernel methods have long provided powerful tools for generalizing linear statistical approaches to nonlinear settings, through an embedding of the sample to a high dimensional feature space, namely a reproducing kernel Hilbert space (RKHS) . Examples include support vector machines, kernel PCA, and kernel CCA, among others. In these cases, data are mapped via a canonical feature map to a reproducing kernel Hilbert space (of high or even infinite dimension), in which the linear operations that define the algorithms are implemented. The inner product between feature mappings need never be computed explicitly, but is given by a positive definite kernel function unique to the RKHS: this permits efficient computation without the need to deal explicitly with the feature representation.

The mappings of individual points to a feature space may be generalized to mappings of probability measures [e.g. 3, Chapter 4]. We call such mappings the kernel means of the underlying random variables. With an appropriate choice of positive definite kernel, the kernel mean on the RKHS uniquely determines the distribution of the variable , and statistical inference problems on distributions can be solved via operations on the kernel means. Applications of this approach include homogeneity testing , where the empirical means on the RKHS are compared directly, and independence testing , where the mean of the joint distribution on the feature space is compared with that of the product of the marginals. Representations of conditional dependence may also be defined in RKHS, and have been used in conditional independence tests .

In this paper, we propose a novel, nonparametric approach to Bayesian inference, making use of kernel means of probabilities. In applying Bayes’ rule, we compute the posterior probability of $x$ in $\mathcal{X}$ given observation $y$ in $\mathcal{Y}$;

where $\pi(x)$ and $p(y|x)$ are the density functions of the prior and the likelihood of $y$ given $x$, respectively, with respective base measures $\nu_{\mathcal{X}}$ and $\nu_{\mathcal{Y}}$, and the normalization factor $q_{\mathcal{Y}}$(y) is given by

Our main result is a nonparametric estimate of the kernel mean posterior, given kernel mean representations of the prior and likelihood.

A valuable property of the kernel Bayes’ rule is that the kernel posterior mean is estimated nonparametrically from data; specifically, the prior and the likelihood are represented in the form of samples from the prior and the joint probability that gives the likelihood, respectively. This confers an important benefit: we can still perform Bayesian inference by making sufficient observations on the system, even in the absence of a specific parametric model of the relation between variables. More generally, if we can sample from the model, we do not require explicit density functions for inference. Such situations are typically seen when the prior or likelihood is given by a random process: Approximate Bayesian Computation is widely applied in population genetics, where the likelihood is given by a branching process, and nonparametric Bayesian inference often uses a process prior with sampling methods. Alternatively, a parametric model may be known, however it might be of sufficient complexity to require Markov chain Monte Carlo or sequential Monte Carlo for inference. The present kernel approach provides an alternative strategy for Bayesian inference in these settings. We demonstrate rates of consistency for our posterior kernel mean estimate, and for the expectation of functions computed using this estimate.

An alternative to the kernel mean representation would be to use nonparametric density estimates for the posterior. Classical approaches include kernel density estimation (KDE) or distribution estimation on a finite partition of the domain. These methods are known to perform poorly on high dimensional data, however. By contrast, the proposed kernel mean representation is defined as an integral or moment of the distribution, taking the form of a function in an RKHS. Thus, it is more akin to the characteristic function approach (see e.g. ) to representing probabilities. A well conditioned empirical estimate of the characteristic function can be difficult to obtain, especially for conditional probabilities. By contrast, the kernel mean has a straightforward empirical estimate, and conditioning and marginalization can be implemented easily, at a reasonable computational cost.

The proposed method of realizing Bayes’ rule is an extension of the approach used in for state-space models. In this earlier work, a heuristic approximation was used, where the kernel mean of the new hidden state was estimated by adding kernel mean estimates from the previous hidden state and the observation. Another relevant work is the belief propagation approach in , which covers the simpler case of a uniform prior.

This paper is organized as follows. We begin in Section 2 with a review of RKHS terminology and of kernel mean embeddings. In Section 3, we derive an expression for Bayes’ rule in terms of kernel means, and provide consistency guarantees. We apply the kernel Bayes’ rule in Section 4 to various inference problems, with numerical results and comparisons with existing methods in Section 5. Our proofs are contained in Section 6 (including proofs of the consistency results of Section 3).

Preliminaries: positive definite kernel and probabilities

Throughout this paper, all Hilbert spaces are assumed to be separable. For an operator $A$ on a Hilbert space, the range is denoted by $\mathcal{R}(A)$. The linear hull of a subset $S$ in a vector space is denoted by ${\rm Span}S$.

A positive definite kernel on $\Omega$ is said to be bounded if there is $M>0$ such that $k(x,x)\leq M$ for any $x\in\Omega$.

Let $(\mathcal{X},\mathcal{B}_{\mathcal{X}})$ be a measurable space, $X$ be a random variable taking values in $\mathcal{X}$ with distribution $P_{X}$, and $k$ be a measurable positive definite kernel on $\mathcal{X}$ such that $E[\sqrt{k(X,X)}]<\infty$. The associated RKHS is denoted by $\mathcal{H}$. The kernel mean $m_{X}^{k}$ (also written $m_{P_{X}}^{k}$) of $X$ on the RKHS $\mathcal{H}$ is defined by the mean of the $\mathcal{H}$-valued random variable $k(\cdot,X)$. The existence of the kernel mean is guaranteed by $E[\|k(\cdot,X)\|]=E[\sqrt{k(X,X)}]<\infty$. We usually write $m_{X}$ for $m_{X}^{k}$ for simplicity, where there is no ambiguity. By the reproducing property, the kernel mean satisfies the relation

for any $f\in\mathcal{H}$. Plugging $f=k(\cdot,u)$ into this relation derives

which shows the explicit functional form. The kernel mean $m_{X}$ is also denoted by $m_{P_{X}}$, as it depends only on the distribution $P_{X}$ with $k$ fixed.

Let $(\mathcal{X},\mathcal{B}_{\mathcal{X}})$ and $(\mathcal{Y},\mathcal{B}_{\mathcal{Y}})$ be measurable spaces, $(X,Y)$ be a random variable on $\mathcal{X}\times\mathcal{Y}$ with distribution $P$, and $k_{\mathcal{X}}$ and $k_{\mathcal{Y}}$ be measurable positive definite kernels with respective RKHS ${\mathcal{H}_{\mathcal{X}}}$ and ${\mathcal{H}_{\mathcal{Y}}}$ such that $E[k_{\mathcal{X}}(X,X)]<\infty$ and $E[k_{\mathcal{Y}}(Y,Y)]<\infty$. The (uncentered) covariance operator $C_{YX}:{\mathcal{H}_{\mathcal{X}}}\to{\mathcal{H}_{\mathcal{Y}}}$ is defined as the linear operator that satisfies

for all $f\in{\mathcal{H}_{\mathcal{X}}},g\in{\mathcal{H}_{\mathcal{Y}}}$. This operator $C_{YX}$ can be identified with $m_{(YX)}$ in the product space ${\mathcal{H}_{\mathcal{Y}}}\otimes{\mathcal{H}_{\mathcal{X}}}$, which is given by the product kernel $k_{\mathcal{Y}}k_{\mathcal{X}}$ on $\mathcal{Y}\times\mathcal{X}$ , by the standard identification between the linear maps and the tensor product. We also define $C_{XX}$ for the operator on ${\mathcal{H}_{\mathcal{X}}}$ that satisfies $\langle f_{2},C_{XX}f_{1}\rangle=E[f_{2}(X)f_{1}(X)]$ for any $f_{1},f_{2}\in{\mathcal{H}_{\mathcal{X}}}$. Similarly to Eq. (4), the explicit integral expressions for $C_{YX}$ and $C_{XX}$ are given by

Throughout this paper, when positive definite kernels on a measurable space are discussed, the following assumption is made:

Positive definite kernels are bounded and measurable.

Under this assumption, the mean and covariance always exist with arbitrary probabilities.

Given i.i.d. sample $(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$ with law $P$, the empirical estimator of the kernel mean and covariance operator are given straightforwardly by

where $\widehat{C}^{(n)}_{YX}$ is written in tensor form. It is known that these estimators are $\sqrt{n}$-consistent in appropriate norms, and $\sqrt{n}(\widehat{m}^{(n)}_{X}-m_{X})$ converges to a Gaussian process on ${\mathcal{H}_{\mathcal{X}}}$ [3, Sec. 9.1]. While we may use non-i.i.d. samples for numerical examples in Section 5, in our theoretical analysis we always assume i.i.d. samples for simplicity.

Kernel expression of Bayes’ rule

Let $(\mathcal{X},\mathcal{B}_{\mathcal{X}})$ and $(\mathcal{Y},\mathcal{B}_{\mathcal{Y}})$ be measurable spaces, $(X,Y)$ be a random variable on $\mathcal{X}\times\mathcal{Y}$ with distribution $P$, and $k_{\mathcal{X}}$ and $k_{\mathcal{Y}}$ be positive definite kernels on $\mathcal{X}$ and $\mathcal{Y}$, respectively, with respective RKHS ${\mathcal{H}_{\mathcal{X}}}$ and ${\mathcal{H}_{\mathcal{Y}}}$. Let $\Pi$ be a probability on $(\mathcal{X},\mathcal{B}_{\mathcal{X}})$, which serves as a prior distribution. For each $x\in\mathcal{X}$, define a probability $P_{Y|x}$ on $(\mathcal{Y},\mathcal{B}_{\mathcal{Y}})$ by $P_{Y|x}(B)=E[I_{B}(Y)|X=x]$, where $I_{B}$ is the index function of a measurable set $B\in\mathcal{B}_{\mathcal{Y}}$. The prior $\Pi$ and the family $\{P_{Y|x}\mid x\in\mathcal{X}\}$ defines the joint distribution $Q$ on $\mathcal{X}\times\mathcal{Y}$ by

for any $A\in\mathcal{B}_{\mathcal{X}}$ and $B\in\mathcal{B}_{\mathcal{Y}}$, and its marginal distribution $Q_{\mathcal{Y}}$ by $Q_{\mathcal{Y}}(B)=Q(\mathcal{X}\times B)$. Throughout this paper, it is assumed that $P_{Y|x}$ and $Q$ are well-defined under some regularity conditions. Let $(Z,W)$ be a random variable on $\mathcal{X}\times\mathcal{Y}$ with distribution $Q$. It is also assumed that the sigma algebra generated by $W$ includes every point $\{y\}$ ($y\in\mathcal{Y}$). For $y\in\mathcal{Y}$, the posterior probability given $y$ is defined by the conditional probability

If the probability distributions have density functions with respect to measures $\nu_{\mathcal{X}}$ on $\mathcal{X}$ and $\nu_{\mathcal{Y}}$ on $\mathcal{Y}$, namely, if the p.d.f. of $P$ and $\Pi$ are given by $p(x,y)$ and $\pi(x)$, respectively, Eq. (6) is reduced to the well known form Eq. (1).

The goal of this subsection is to derive an estimator of the kernel mean of posterior $m_{Q_{\mathcal{X}}|y}$. The following theorem is fundamental to discuss conditional probabilities with positive definite kernels.

If $E[g(Y)|X=\cdot]\in{\mathcal{H}_{\mathcal{X}}}$ holds for $g\in{\mathcal{H}_{\mathcal{Y}}}$, then

If $C_{XX}$ is injective, i.e., if the function $f\in{\mathcal{H}_{\mathcal{X}}}$ with $C_{XX}f=C_{XY}g$ is unique, the above relation can be expressed as

Noting $\langle C_{XX}f,f\rangle=E[f(X)^{2}]$, it is easy to see that $C_{XX}$ is injective, if $\mathcal{X}$ is a topological space, $k_{\mathcal{X}}$ is a continuous kernel, and ${\rm Supp}(P_{X})=\mathcal{X}$, where ${\rm Supp}(P_{X})$ is the support of $P_{X}$.

From Theorem 3.1, we have the following result, which expresses the kernel mean of $Q_{\mathcal{Y}}$.

Let $m_{\Pi}$ and $m_{Q_{\mathcal{Y}}}$ be the kernel means of $\Pi$ in ${\mathcal{H}_{\mathcal{X}}}$ and $Q_{\mathcal{Y}}$ in ${\mathcal{H}_{\mathcal{Y}}}$, respectively. If $C_{XX}$ is injective, $m_{\Pi}\in\mathcal{R}(C_{XX})$, and $E[g(Y)|X=\cdot]\in{\mathcal{H}_{\mathcal{X}}}$ for any $g\in{\mathcal{H}_{\mathcal{Y}}}$, then

Take $f\in{\mathcal{H}_{\mathcal{X}}}$ such that $f=C_{XX}^{-1}m_{\Pi}$. For any $g\in{\mathcal{H}_{\mathcal{Y}}}$, $\langle C_{YX}f,g\rangle=\langle f,C_{XY}g\rangle=\langle f,C_{XX}E[g(Y)|X=\cdot]\rangle=\langle C_{XX}f,E[g(Y)|X=\cdot]\rangle=\langle m_{\Pi},E[g(Y)|X=\cdot]\rangle=\langle m_{Q_{\mathcal{Y}}},g\rangle$, which implies $C_{YX}f=m_{Q_{\mathcal{Y}}}$. ∎

As discussed in , the operator $C_{YX}C_{XX}^{-1}$ can be regarded as the kernel expression of the conditional probability $P_{Y|x}$ or $p(y|x)$.

Note, however, that the assumption $E[g(Y)|X=\cdot]\in{\mathcal{H}_{\mathcal{X}}}$ may not hold in general; we can easily give counterexamples in the case of Gaussian kernelsSuppose that ${\mathcal{H}_{\mathcal{X}}}$ and ${\mathcal{H}_{\mathcal{Y}}}$ are given by Gaussian kernel, and that $X$ and $Y$ are independent. Then, $E[g(Y)|X=x]$ is a constant function of $x$, which is known not to be included in a RKHS given by a Gaussian kernel [38, Corollary 4.44].. In the following, we nonetheless derive a population expression of Bayes’ rule under this strong assumption, use it as a prototype for defining an empirical estimator, and prove its consistency.

In deriving kernel realization of Bayes’ rule, we will use the following tensor representation of the joint probability $Q$, based on Theorem 3.2:

In the above equation, the covariance operator $C_{(YX)X}:{\mathcal{H}_{\mathcal{X}}}\to{\mathcal{H}_{\mathcal{Y}}}\otimes{\mathcal{H}_{\mathcal{X}}}$ is defined by the random variable $((Y,X),X)$ taking values on $(\mathcal{Y}\times\mathcal{X})\times\mathcal{X}$.

In many applications of Bayesian inference, the probability conditioned on a particular value should be computed. By plugging the point measure at $x$ into $\Pi$ in Eq. (8), we have a population expression

If we replace $P$ by $Q$ and $x$ by $y$ in Eq. (10), we obtain

This is exactly the kernel mean expression of the posterior, and the next step is to provide a way of deriving the covariance operators $C_{ZW}$ and $C_{WW}$. Recall that the kernel mean $m_{Q}=m_{(ZW)}\in{\mathcal{H}_{\mathcal{X}}}\otimes{\mathcal{H}_{\mathcal{Y}}}$ can be identified with the covariance operator $C_{ZW}:{\mathcal{H}_{\mathcal{Y}}}\to{\mathcal{H}_{\mathcal{X}}}$, and $m_{(WW)}$, which is the kernel mean on the product space ${\mathcal{H}_{\mathcal{Y}}}\otimes{\mathcal{H}_{\mathcal{Y}}}$, with $C_{WW}$. Then from Eq. (9) and the similar expression $m_{(WW)}=C_{(YY)X}C_{XX}^{-1}m_{\Pi}$, we are able to obtain the operators in Eq. (11), and thus the kernel mean of the posterior.

The above argument can be rigorously implemented, if empirical estimators are considered. Let $(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$ be an i.i.d. sample with law $P$. Since the kernel method needs to express the information of variables in terms of Gram matrices given by data points, we assume that the prior is also expressed in the form of an empirical estimate, and that we have a consistent estimator of $m_{\Pi}$ in the form

where $\varepsilon_{n}$ is the coefficient of the Tikhonov-type regularization for operator inversion, and $I$ is the identity operator. The empirical estimators $\widehat{C}_{ZW}$ and $\widehat{C}_{WW}$ for $C_{ZW}$ and $C_{WW}$ are identified with $\widehat{m}_{(ZW)}$ and $\widehat{m}_{(WW)}$, respectively. In the following, $G_{X}$ and $G_{Y}$ denote the Gram matrices $(k_{\mathcal{X}}(X_{i},X_{j}))$ and $(k_{\mathcal{Y}}(Y_{i},Y_{j}))$, respectively, and $I_{n}$ is the identity matrix of size $n$.

The Gram matrix expressions of $\widehat{C}_{ZW}$ and $\widehat{C}_{WW}$ are given by

The proof is similar to that of Proposition 3.4 below, and is omitted. The expressions in Proposition 3.3 imply that the probabilities $Q$ and $Q_{\mathcal{Y}}$ are estimated by the weighted samples $\{((X_{i},Y_{i}),\widehat{\mu}_{i})\}_{i=1}^{n}$ and $\{(Y_{i},\widehat{\mu}_{i})\}_{i=1}^{n}$, respectively, with common weights. Since the weight $\widehat{\mu}_{i}$ may be negative, in applying Eq. (11) the operator inversion in the form $(\widehat{C}_{WW}+\delta_{n}I)^{-1}$ may be impossible or unstable. We thus use another type of Tikhonov regularization, thus obtaining the estimator

For any $y\in\mathcal{Y}$, the Gram matrix expression of $\widehat{m}_{Q_{\mathcal{X}}|y}$ is given by

where $\Lambda={\rm diag}(\widehat{\mu})$ is a diagonal matrix with elements $\widehat{\mu}_{i}$ in Eq. (12), ${\bf k}_{X}=(k_{\mathcal{X}}(\cdot,X_{1}),\ldots,k_{\mathcal{X}}(\cdot,X_{n}))^{T}\in{\mathcal{H}_{\mathcal{X}}}^{n}$, and ${\bf k}_{Y}=(k_{\mathcal{Y}}(\cdot,Y_{1}),\ldots,k_{\mathcal{Y}}(\cdot,Y_{n}))^{T}\in{\mathcal{H}_{\mathcal{Y}}}^{n}$.

Let $h=(\widehat{C}_{WW}^{2}+\delta_{n}I)^{-1}\widehat{C}_{WW}k_{\mathcal{Y}}(\cdot,y)$, and decompose it as $h=\sum_{i=1}^{n}\alpha_{i}k_{\mathcal{Y}}(\cdot,Y_{i})+h_{\perp}=\alpha^{T}{\bf k}_{Y}+h_{\perp}$, where $h_{\perp}$ is orthogonal to ${\rm Span}\{k_{\mathcal{Y}}(\cdot,Y_{i})\}_{i=1}^{n}$. Expansion of $(\widehat{C}_{WW}^{2}+\delta_{n}I)h=\widehat{C}_{WW}k_{\mathcal{Y}}(\cdot,y)$ gives ${\bf k}_{Y}^{T}(\Lambda G_{Y})^{2}\alpha+\delta_{n}{\bf k}_{Y}^{T}\alpha+\delta_{n}h_{\perp}={\bf k}_{Y}^{T}\Lambda{\bf k}_{Y}(y)$. Taking the inner product with $k_{\mathcal{Y}}(\cdot,Y_{j})$, we have

The coefficient $\rho$ in $\widehat{m}_{Q_{\mathcal{X}}|y}=\widehat{C}_{ZW}h=\sum_{i=1}^{n}\rho_{i}k_{\mathcal{X}}(\cdot,X_{i})$ is given by $\rho=\Lambda G_{Y}\alpha$, and thus

We call Eqs.(13) and (14) the kernel Bayes’ rule (KBR). The required computations are summarized in Figure 1. The KBR uses a weighted sample to represent the posterior; it is similar in this respect to sampling methods such as importance sampling and sequential Monte Carlo (). The KBR method, however, does not generate samples of the posterior, but updates the weights of a sample by matrix computation. We will give some experimental comparisons between KBR and sampling methods in Section 5.1.

If our aim is to estimate the expectation of a function $f\in{\mathcal{H}_{\mathcal{X}}}$ with respect to the posterior, the reproducing property Eq. (3) gives an estimator

2 Consistency of the KBR estimator

where $\mathbf{f}_{X}^{T}R_{X|Y}\mathbf{k}_{Y}(y)$ is given by Eq. (15).

It is possible to extend the covariance operator $C_{WW}$ to one defined on $L^{2}(Q_{\mathcal{Y}})$ by

If we consider the convergence on average over $y$, we have a slightly better rate on the consistency of the KBR estimator in $L^{2}(Q_{\mathcal{Y}})$.

Bayesian inference with Kernel Bayes’ Rule

In Bayesian inference, we are usually interested in finding a point estimate such as the MAP solution, the expectation of a function under the posterior, or other properties of the distribution. Given that KBR provides a posterior estimate in the form of a kernel mean (which uniquely determines the distribution when a characteristic kernel is used), we now describe how our kernel approach applies to problems in Bayesian inference.

First, we have already seen that a consistent estimator for the expectation of $f\in{\mathcal{H}_{\mathcal{X}}}$ can be defined with respect to the posterior. On the other hand, unless $f\in{\mathcal{H}_{\mathcal{X}}}$ holds, there is no theoretical guarantee that it gives a good estimate. In Section 5.1, we discuss some experimental results in such situations.

To obtain a point estimate of the posterior on $x$, it is proposed in to use the preimage $\widehat{x}=\arg\min_{x}\|k_{\mathcal{X}}(\cdot,x)-{\bf k}_{X}^{T}R_{X|Y}{\bf k}_{Y}(y)\|^{2}_{{\mathcal{H}_{\mathcal{X}}}}$, which represents the posterior mean most effectively by one point. We use this approach in the present paper when point estimates are considered. In the case of the Gaussian kernel $\exp(-\|x-y\|^{2}/(2\sigma^{2}))$, the fixed point method

where $\rho=R_{X|Y}{\bf k}_{Y}(y)$, can be used to optimize $x$ sequentially . This method usually converges very fast, although no theoretical guarantee exists for the convergence to the globally optimal point, as is usual in non-convex optimization.

A notable property of KBR is that the prior and likelihood are represented in terms of samples. Thus, unlike many approaches to Bayesian inference, precise knowledge of the prior and likelihood distributions is not needed, once samples are obtained. The following are typical situations where the KBR approach is advantageous:

The probabilistic relation among variables is difficult to realize with a simple parametric model, while we can obtain samples of the variables easily. We will see such an example in Section 4.3.

The probability density function of the prior and/or likelihood is hard to obtain explicitly, but sampling is possible:

In the field of population genetics, Bayesian inference is used with a likelihood expressed by branching processes to model the split of species, for which the explicit density is hard to obtain. Approximate Bayesian Computation (ABC) is a popular method for approximately sampling from a posterior without knowing the functional form .

Another interesting application along these lines is nonparametric Bayesian inference ( and references therein), in which the prior is typically given in the form of a process without a density form. In this case, sampling methods are often applied ( among others). Alternatively, the posterior may be approximated using variational methods .

We will present an experimental comparison of KBR and ABC in Section 5.2.

Even if explicit forms for the likelihood and prior are available, and standard sampling methods such as MCMC or sequential MC are applicable, the computation of a posterior estimate given $y$ might still be computationally costly, making real-time applications unfeasible. Using KBR, however, the expectation of a function of the posterior given different $y$ is obtained simply by taking the inner product as in Eq. (15), once $\mathbf{f}_{X}^{T}R_{X|Y}$ has been computed.

2 Discussions concerning implementation

When implementing KBR, a number of factors should be borne in mind to ensure good performance. First, in common with many nonparametric approaches, KBR requires training data in the region of the new “test” points for results to be meaningful. In other words, if the point on which we condition appears in a region far from the sample used for the estimation, the posterior estimator will be unreliable.

Second, in computing the posterior in KBR, Gram matrix inversion is necessary, which would cost $O(n^{3})$ for sample size $n$ if attempted directly. Substantial cost reductions can be achieved if the Gram matrices are approximated by low rank matrix approximations. A popular choice is the incomplete Cholesky decomposition , which approximates a Gram matrix in the form of $\Gamma\Gamma^{T}$ with $n\times r$ matrix $\Gamma$ ($r\ll n$) at cost $O(nr^{2})$. Using this and the Woodbury identity, the KBR can be approximately computed at cost $O(nr^{2})$.

Third, kernel choice or model selection is key to the effectiveness of any kernel method. In the case of KBR, we have three model parameters: the kernel (or its parameter, e.g. the bandwidth), the regularization parameter $\varepsilon_{n}$, and $\delta_{n}$. The strategy for parameter selection depends on how the posterior is to be used in the inference problem. If it is to be applied in regression, we can use standard cross-validation. In the filtering experiments in Section 5, we use a validation method where we divide the training sample in two.

A more general model selection approach can also be formulated, by creating a new regression problem for the purpose. Suppose the prior $\Pi$ is given by the marginal $P_{X}$ of $P$. The posterior ${Q}_{\mathcal{X}|y}$ averaged with respect to $P_{Y}$ is then equal to the marginal $P_{X}$ itself. We are thus able to compare the discrepancy of the empirical kernel mean of $P_{X}$ and the average of the estimators $\widehat{m}_{Q_{\mathcal{X}|y=Y_{i}}}$ over $Y_{i}$. This leads to a $K$-fold cross validation approach: for a partition of $\{1,\ldots,n\}$ into $K$ disjoint subsets $\{T_{a}\}_{a=1}^{K}$, let $\widehat{m}_{Q_{\mathcal{X}|y}}^{[-a]}$ be the kernel mean of posterior computed using Gram matrices on data $\{(X_{i},Y_{i})\}_{i\notin T_{a}}$, and based on the prior mean $\widehat{m}_{X}^{[-a]}$ with data $\{X_{i}\}_{i\notin T_{a}}$. We can then cross validate by minimizing \sum_{a=1}^{K}\bigl{\|}\frac{1}{|T_{a}|}\sum_{j\in T_{a}}\widehat{m}_{Q_{\mathcal{X}|y=Y_{j}}}^{[-a]}-\widehat{m}_{X}^{[a]}\bigr{\|}^{2}_{{\mathcal{H}_{\mathcal{X}}}}, where $\widehat{m}_{X}^{[a]}=\frac{1}{|T_{a}|}\sum_{j\in T_{a}}k_{\mathcal{X}}(\cdot,X_{j})$.

3 Application to nonparametric state-space model

We next describe how KBR may be used in a particular application: namely, inference in a general time invariant state-space model,

where $Y_{t}$ is an observable variable, and $X_{t}$ is a hidden state variable. We begin with a brief review of alternative strategies for inference in state-space models with complex dynamics, for which linear models are not suitable. The extended Kalman filter (EKF) and unscented Kalman filter (UKF, ) are nonlinear extensions of the standard linear Kalman filter, and are well established in this setting. Alternatively, nonparametric estimates of conditional density functions can be employed, including kernel density estimation or distribution estimates on a partitioning of the space . The latter nonparametric approaches are effective only for low-dimensional cases, however. Most relevant to this paper are and , in which the kernel means and covariance operators are used to implement the nonparametric HMM.

In this paper, we apply the KBR for inference in the nonparametric state-space model. We do not assume the conditional probabilities $p(Y_{t}|X_{t})$ and $q(X_{t+1}|X_{t})$ to be known explicitly, nor do we estimate them with simple parametric models. Rather, we assume a sample $(X_{1},Y_{1}),\ldots,(X_{T+1},Y_{T+1})$ is given for both the observable and hidden variables in the training phase. The conditional probability for observation process $p(y|x)$ and the transition $q(x_{t+1}|x_{t})$ are represented by the empirical covariance operators as computed on the training sample,

While the sample is not i.i.d., we can use the empirical covariances, which are consistent by the mixing property of Markov models.

where the coefficients $\widehat{\mu}^{(t+1)}=(\widehat{\mu}^{(t+1)}_{i})_{i=1}^{T}$ are given by

In sequential filtering, a substantial reduction in computational cost can be achieved by low rank matrix approximations, as discussed above. Given an approximation of rank $r$ for the Gram matrices and transfer matrix, and employing the Woodbury identity, the computation costs just $O(Tr^{2})$ for each time step.

4 Bayesian computation without likelihood

We next address the setting where the likelihood is not known in analytic form, but sampling is possible. In this case, Approximate Bayesian Computation (ABC) is a popular method for Bayesian inference. The simplest form of ABC, which is called the rejection method, generates a sample from $q(Z|W=y)$ as follows: (i) generate a sample $X_{t}$ from the prior $\Pi$, (ii) generate a sample $Y_{t}$ from $P(Y|X_{t})$, (iii) if $D(y,Y_{t})<\tau$, accept $X_{t}$; otherwise reject, (iv) go to (i). In step (iii), $D$ is a distance measure of the space $\mathcal{X}$, and $\tau$ is tolerance to acceptance.

In the same setting as ABC, KBR gives the following sampling-based method for computing the kernel posterior mean:

Generate a sample $X_{1},\ldots,X_{n}$ from the prior $\Pi$.

Generate a sample $Y_{t}$ from $P(Y|X_{t})$ ($t=1,\ldots,n$).

Compute Gram matrices $G_{X}$ and $G_{Y}$ with $(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$, and $R_{X|Y}{\bf k}_{Y}(y)$.

Alternatively, since $(X_{t},Y_{t})$ is an sample from $Q$, it is possible to use Eq. (10) for the kernel mean of the conditional probability $q(x|y)$. As in , the estimator is given by

The distribution of a sample generated by ABC approaches to the true posterior if $\tau$ goes to zero, while empirical estimates via the kernel approaches converge to the true posterior mean in the limit of infinite sample size. The efficiency of ABC, however, can be arbitrarily poor for small $\tau$, since a sample $X_{t}$ is then rarely accepted in Step (iii).

The ABC method generates a sample, hence any statistics based on the posterior can be approximated. Given a posterior mean obtained by one of the kernel methods, however, we may only obtain expectations of functions in the RKHS, meaning that certain statistics (such as confidence intervals) are not straightforward to obtain. In Section 5.2, we present an experimental evaluation of the trade-off between computation time and accuracy for ABC and KBR.

Numerical Examples

2 Bayesian computation without likelihood

We compare ABC and the kernel methods, KBR and conditional mean, in terms of estimation accuracy and computational time, since they have an obvious tradeoff. To compute the estimation accuracy rigorously, the ground truth is needed: thus we use Gaussian distributions for the true prior and likelihood, which makes the posterior easy to compute in closed form. The samples are taken from the same model used in Section 5.1, and $\int xq(x|y)dx$ is evaluated at 10 different points of $y$. We performed 10 random runs with different random generation of the true distributions.

For ABC, we used only the rejection method; while there are more advanced sampling schemes , their implementation is dependent on the problem being solved. Various values for the acceptance region $\tau$ are used, and the accuracy and computational time are shown in Fig. 3 together with total sizes of the generated samples. For the kernel methods, the sample size $n$ is varied. The regularization parameters are given by $\varepsilon_{n}=0.01/n$ and $\delta_{n}=2\varepsilon_{n}$ for KBR, and $\varepsilon_{n}=0.01/\sqrt{n}$ for the conditional kernel mean. The kernels in the kernel methods are Gaussian kernels for which the bandwidth parameters are chosen by the median of the pairwise distances on the data (). The incomplete Cholesky decomposition is employed for the low-rank approximation. The results indicate that kernel methods achieve more accurate results than ABC at a given computational cost, and the conditional kernel mean shows better results.

3 Filtering problems

We next compare the KBR filtering method (proposed in Section 4.3) with EKF and UKF on synthetic data.

KBR has the regularization parameters $\varepsilon_{T},\delta_{T}$, and kernel parameters for $k_{\mathcal{X}}$ and $k_{\mathcal{Y}}$ (e.g., the bandwidth parameter for an RBF kernel). Under the assumption that a training sample is available, cross-validation can be performed on the training sample to select the parameters. By dividing the training sample into two, one half is used to estimate the covariance operators Eq. (17) with a candidate parameter set, and the other half to evaluate the estimation errors. To reduce the search space and attendant computational cost, we used a simpler procedure, setting $\delta_{T}=2\varepsilon_{T}$, and using the Gaussian kernel bandwidths $\beta\sigma_{\mathcal{X}}$ and $\beta\sigma_{\mathcal{Y}}$, where $\sigma_{\mathcal{X}}$ and $\sigma_{\mathcal{Y}}$ are the median of pairwise distances in the training samples (). This leaves only two parameters $\beta$ and $\varepsilon_{T}$ to be tuned.

where $\eta>0$ is an increment of the angle and $\zeta_{t}\sim N(0,\sigma_{h}^{2}I_{2})$ is independent process noise. Note that the dynamics of $(u_{t},v_{t})$ are nonlinear even for $b=0$. The observation $Y_{t}$ follows

where $\xi_{t}$ is independent noise. The two dynamics are defined as follows. (a) (rotation with noisy observation) $\eta=0.3$, $b=0$, $\sigma_{h}=\sigma_{o}=0.2$. (b) (oscillatory rotation with noisy observation) $\eta=0.4$, $b=0.4$, $M=8$, $\sigma_{h}=\sigma_{o}=0.2$. (See Fig.5).

We assume the correct dynamics are known to the EKF and UKF. The results are shown in Fig. 4. In all the cases, EKF and UKF show unrecognizably small difference. The dynamics in (a) are weakly nonlinear, and KBR has slightly worse MSE than EKF and UKF. For dataset (b), which has strong nonlinearity, KBR outperforms the nonlinear Kalman filter for $T\geq 200$.

In our second synthetic example, we applied the KBR filter to the camera rotation problem used in Song et al. . The angle of a camera, which is located at a fixed position, is a hidden variable, and movie frames recorded by the camera are observed. The data are generated virtually using a computer graphics environment. As in , we are given 3600 downsampled frames of $20\times 20$ RGB pixels ($Y_{t}\in^{1200}$), where the first 1800 frames are used for training, and the second half are used to test the filter. We make the data noisy by adding Gaussian noise $N(0,\sigma^{2})$ to $Y_{t}$.

Our experiments cover two settings. In the first, we assume we do not know that the hidden state $S_{t}$ is included in $SO(3)$, but only that it is a general $3\times 3$ matrix. In this case, we use the Kalman filter by estimating the relations under a linear assumption, and the KBR filter with Gaussian kernels for $S_{t}$ and $X_{t}$ as Euclidean vectors. In the second setting, we exploit the fact that $S_{t}\in SO(3)$: for the Kalman Filter, $S_{t}$ is represented by a quanternion, which is a standard vector representation of rotations; for the KBR filter the kernel $k(A,B)={\rm Tr}[AB^{T}]$ is used for $S_{t}$, and $S_{t}$ is estimated within $SO(3)$. Table 1 shows the Frobenius norms between the estimated matrix and the true one. The KBR filter significantly outperforms the EKF, since KBR has the advantage in extracting the complex nonlinear dependence between the observation and the hidden state.

Proofs

The proof idea for the consistency rates of the KBR estimators is similar to , in which the basic techniques are taken from the general theory of regularization .

The first preliminary result is a rate of convergence for the mean transition in Theorem 3.2. In the following $\mathcal{R}(C_{XX}^{0})$ means ${\mathcal{H}_{\mathcal{X}}}$.

Assume that $\pi/p_{X}\in\mathcal{R}(C_{XX}^{\beta})$ for some $\beta\geq 0$, where $\pi$ and $p_{X}$ are the p.d.f. of $\Pi$ and $P_{X}$, respectively. Let $\widehat{m}_{\Pi}^{(n)}$ be an estimator of $m_{\Pi}$ such that $\|\widehat{m}_{\Pi}^{(n)}-m_{\Pi}\|_{\mathcal{H}_{\mathcal{X}}}=O_{p}(n^{-\alpha})$ as $n\to\infty$ for some $0<\alpha\leq 1/2$. Then, with $\varepsilon_{n}=n^{-\max\{\frac{2}{3}\alpha,\frac{\alpha}{1+\beta}\}}$, we have

Take $\eta\in{\mathcal{H}_{\mathcal{X}}}$ such that $\pi/p_{X}=C_{XX}^{\beta}\eta$. Then, we have

First we show the rate of the estimation error:

as $n\to\infty$. By using $B^{-1}-A^{-1}=B^{-1}(A-B)A^{-1}$ for any invertible operators $A$ and $B$, the left hand side of Eq. (21) is upper bounded by

By the decomposition $\widehat{C}^{(n)}_{YX}=\widehat{C}_{YY}^{(n)1/2}\widehat{W}_{YX}^{(n)}\widehat{C}_{XX}^{(n)1/2}$ with $\|\widehat{W}_{YX}^{(n)}\|\leq 1$ , we have \|\widehat{C}^{(n)}_{YX}\bigl{(}\widehat{C}^{(n)}_{XX}+\varepsilon_{n}I\bigr{)}^{-1}\|=O_{p}(\varepsilon_{n}^{-1/2}), which implies the first term is of $O_{p}(n^{-\alpha}\varepsilon_{n}^{-1/2})$. From the $\sqrt{n}$ consistency of the covariance operators and $m_{\Pi}=C_{XX}^{\beta+1}\eta$, a similar argument to the first term proves that the second and third terms are of the order $O_{p}(n^{-1/2})$ and $O_{p}(n^{-1/2}\varepsilon_{n}^{-1/2})$, respectively, which means Eq. (21).

Next, we show the rate for the approximation error

Let $C_{YX}=C_{YY}^{1/2}W_{YX}C_{XX}^{1/2}$ be the decomposition with $\|W_{YX}\|\leq 1$. It follows from Eq. (20) and the relation

that the left hand side of Eq. (22) is upper bounded by

By the eigendecomposition $C_{XX}=\sum_{i}\lambda_{i}\phi_{i}\langle\phi_{i},\cdot\rangle$, where $\{\lambda_{i}\}$ are the positive eigenvalues and $\{\phi_{i}\}$ are the corresponding unit eigenvectors, the expansion

holds. If $0\leq\beta<1/2$, we have $\frac{\varepsilon_{n}\lambda_{i}^{(2\beta+1)/2}}{\lambda_{i}+\varepsilon_{n}}=\frac{\lambda_{i}^{(2\beta+1)/2}}{(\lambda_{i}+\varepsilon_{n})^{(2\beta+1)/2}}\frac{\varepsilon_{n}^{(1-2\beta)/2}}{(\lambda_{i}+\varepsilon_{n})^{(1-2\beta)/2}}\varepsilon_{n}^{(2\beta+1)/2}\leq\varepsilon_{n}^{(2\beta+1)/2}$. If $\beta\geq 1/2$, then $\frac{\varepsilon_{n}\lambda_{i}^{(2\beta+1)/2}}{\lambda_{i}+\varepsilon_{n}}\leq\|C_{XX}\|\varepsilon_{n}$. The dominated convergence theorem shows that the the above sum converges to zero of the order $O(\varepsilon_{n}^{\min\{2\beta+1,2\}})$ as $\varepsilon_{n}\to 0$.

From Eqs. (21) and (22), the optimal order of $\varepsilon_{n}$ and the optimal rate of consistency are given as claimed. ∎

The following theorem shows the consistency rate of the estimator used in the conditioning step Eq. (11).

Let $f$ be a function in ${\mathcal{H}_{\mathcal{X}}}$, and $(Z,W)$ be a random variable taking values in $\mathcal{X}\times\mathcal{Y}$. Assume that $E[f(Z)|W=\cdot]\in\mathcal{R}(C_{WW}^{\nu})$ for some $\nu\geq 0$, and $\widehat{C}^{(n)}_{WZ}:{\mathcal{H}_{\mathcal{X}}}\to{\mathcal{H}_{\mathcal{Y}}}$ and $\widehat{C}^{(n)}_{WW}:{\mathcal{H}_{\mathcal{Y}}}\to{\mathcal{H}_{\mathcal{Y}}}$ be compact operators, which may not be positive definite, such that $\|\widehat{C}^{(n)}_{WZ}-C_{WZ}\|=O_{p}(n^{-\gamma})$ and $\|\widehat{C}^{(n)}_{WW}-C_{WW}\|=O_{p}(n^{-\gamma})$ for some $\gamma>0$. Then, for a positive sequence $\delta_{n}=n^{-\max\{\frac{4}{9}\gamma,\frac{4}{2\nu+5}\gamma\}}$, we have as $n\to\infty$

Let $\eta\in{\mathcal{H}_{\mathcal{X}}}$ such that $E[f(Z)|W=\cdot]=C_{WW}^{\nu}\eta$. First we show

The left hand side of Eq. (23) is upper bounded by

Let $\widehat{C}^{(n)}_{WW}=\sum_{i}\lambda_{i}\phi_{i}\langle\phi_{i},\cdot\rangle$ be the eigendecomposition, where $\{\phi_{i}\}$ is the unit eigenvectors and $\{\lambda_{i}\}$ is the corresponding eigenvalues. From \bigl{|}\lambda_{i}/(\lambda_{i}^{2}+\delta_{n})\bigr{|}=1/|\lambda_{i}+\delta_{n}/\lambda_{i}|\leq 1/(2\sqrt{|\lambda_{i}|}\sqrt{\delta_{n}/|\lambda_{i}|})=1/(2\sqrt{\delta_{n}}), we have \|\widehat{C}^{(n)}_{WW}\bigl{(}(\widehat{C}^{(n)}_{WW})^{2}+\delta_{n}I\bigr{)}^{-1}\|\leq 1/(2\sqrt{\delta_{n}}), and thus the first term of the above bound is of $O_{p}(n^{-\gamma}\delta_{n}^{-1/2})$. A similar argument by the eigendecomposition of $C_{WW}$ combined with the decomposition $C_{WZ}=C_{WW}^{1/2}U_{WZ}C_{ZZ}^{1/2}$ with $\|U_{WZ}\|\leq 1$ shows that the second term is of $O_{p}(n^{-\gamma}\delta_{n}^{-3/4})$. From the fact $\|(\widehat{C}^{(n)}_{WW})^{2}-C_{WW}^{2}\|\leq\|\widehat{C}^{(n)}_{WW}(\widehat{C}^{(n)}_{WW}-C_{WW})\|+\|(\widehat{C}^{(n)}_{WW}-C_{WW})C_{WW}\|=O_{p}(n^{-\gamma})$, the third term is of $O_{p}(n^{-\gamma}\delta_{n}^{-5/4})$. This implies Eq. (23).

From $E[f(Z)|W=\cdot]=C_{WW}^{\nu}\eta$ and $C_{WZ}f=C_{WW}E[f(Z)|W=\cdot]=C_{WW}^{\nu+1}\eta$, the convergence rate

can be proved by the same way as Eq. (22).

Combination of Eqs.(23) and (24) proves the assertion. ∎

Note that for $f,g\in{\mathcal{H}_{\mathcal{X}}}$ we have $(f,g)_{L^{2}(Q_{\mathcal{Y}})}=E[f(W)g(W)]=\langle f,C_{WW}g\rangle_{\mathcal{H}_{\mathcal{X}}}$. It follows that the left hand side of the assertion is equal to

First, by the similar argument to the proof of Eq. (23), it is easy to show that the rate of the estimation error is given by

The consistency of KBR follows by combining the above theorems.

Let $f$ be a function in ${\mathcal{H}_{\mathcal{X}}}$, $(Z,W)$ be a random variable that has the distribution $Q$ with p.d.f. $p(y|x)\pi(x)$, and $\widehat{m}_{\Pi}^{(n)}$ be an estimator of $m_{\Pi}$ such that $\|\widehat{m}_{\Pi}^{(n)}-m_{\Pi}\|_{\mathcal{H}_{\mathcal{X}}}=O_{p}(n^{-\alpha})$ ($n\to\infty$) for some $0<\alpha\leq 1/2$. Assume that $\pi/p_{X}\in\mathcal{R}(C_{XX}^{\beta})$ with $\beta\geq 0$, and $E[f(Z)|W=\cdot]\in\mathcal{R}(C_{WW}^{\nu})$ for some $\nu\geq 0$. For the regularization constants $\varepsilon_{n}=n^{-\max\{\frac{2}{3}\alpha,\frac{1}{1+\beta}\alpha\}}$ and $\delta_{n}=n^{-\max\{\frac{4}{9}\gamma,\frac{4}{2\nu+5}\gamma\}}$, where $\gamma=\min\{\frac{2}{3}\alpha,\frac{2\beta+1}{2\beta+2}\alpha\}$, we have for any $y\in\mathcal{Y}$

where $\mathbf{f}_{X}^{T}R_{X|Y}\mathbf{k}_{Y}(y)$ is given by Eq. (14).

By applying Theorem 6.1 to $Y=(Y,X)$ and $Y=(Y,Y)$, we see that both of $\|\widehat{C}_{WZ}-C_{WZ}\|$ and $\|\widehat{C}_{WW}-C_{WW}\|$ are of $O_{p}(n^{-\gamma})$. Since

combination of Theorems 6.1 and 6.2 proves the theorem. ∎

The next theorem shows the rate on average w.r.t. $Q_{\mathcal{Y}}$. The proof is similar to the above theorem, and omitted.

We also have consistency of the estimator for the kernel mean of posterior $m_{Q_{\mathcal{X}|y}}$, if we make stronger assumptions. First, we formulate the expectation with the posterior in terms of operators. Let $(Z,W)$ be a random variable with distribution $Q$. Assume that for any $f\in{\mathcal{H}_{\mathcal{X}}}$ the conditional expectation $E[f(Z)|W=\cdot]$ is included in ${\mathcal{H}_{\mathcal{Y}}}$. We then have a linear operator $S$ defined by

If we further assume that $S$ is bounded, the adjoint operator $S^{*}:{\mathcal{H}_{\mathcal{Y}}}\to{\mathcal{H}_{\mathcal{X}}}$ satisfies

for any $y\in\mathcal{Y}$, and thus $S^{*}k_{\mathcal{Y}}(\cdot,y)$ is equal to the kernel mean of the conditional probability of $Z$ given $W=y$.

We make the following further assumptions: Assumption (S)

The covariance operator $C_{WW}$ is injective.

There exists $\nu>0$ such that for any $f\in{\mathcal{H}_{\mathcal{X}}}$ there is $\eta_{f}\in{\mathcal{H}_{\mathcal{X}}}$ with $Sf=C_{WW}^{\nu}\eta_{f}$, and the linear map

Let $(Z,W)$ be a random variable that has the distribution $Q$ with p.d.f. $p(y|x)\pi(x)$, and $\widehat{m}_{\Pi}^{(n)}$ be an estimator of $m_{\Pi}$ such that $\|\widehat{m}_{\Pi}^{(n)}-m_{\Pi}\|_{\mathcal{H}_{\mathcal{X}}}=O_{p}(n^{-\alpha})$ ($n\to\infty$) for some $0<\alpha\leq 1/2$. Assume (S) above, and $\pi/p_{X}\in\mathcal{R}(C_{XX}^{\beta})$ with some $\beta\geq 0$. For the regularization constants $\varepsilon_{n}=n^{-\max\{\frac{2}{3}\alpha,\frac{1}{1+\beta}\alpha\}}$ and $\delta_{n}=n^{-\max\{\frac{4}{9}\gamma,\frac{4}{2\nu+5}\gamma\}}$, where $\gamma=\min\{\frac{2}{3}\alpha,\frac{2\beta+1}{2\beta+2}\alpha\}$, we have for any $y\in\mathcal{Y}$

as $n\to\infty$, where $m_{Q_{\mathcal{X}}|y}$ is the kernel mean of the posterior given $y$.

First, in a similar manner to the proof of Eq. (23), we have

is proved. The left hand side of Eq. (25) is upper-bounded by

It follows from Theorem 3.1 that $C_{WZ}=C_{WW}S$, and thus $\|C_{WW}(C_{WW}^{2}+\delta_{n}I)^{-1}C_{WZ}-S\|=\|C_{WW}(C_{WW}^{2}+\delta_{n}I)^{-1}C_{WW}S-S\|\leq\delta_{n}\|(C_{WW}^{2}+\delta_{n}I)^{-1}C_{WW}^{\nu}\|\,\|C_{WW}^{-\nu}S\|$. The eigendecomposition of $C_{WW}$ together with the inequality $\frac{\delta_{n}\lambda^{\nu}}{\lambda^{2}+\delta_{n}}\leq\delta_{n}^{\min\{1,\nu/2\}}$ ($\lambda\geq 0$) completes the proof. ∎

We thank Arnaud Doucet, Lorenzo Rosasco, Yee Whye Teh and Shuhei Mano for their helpful comments.

References