---
abstract: 'We propose Bayesian extensions of two nonparametric regression methods which are kernel and mutual $k$-nearest neighbor regression methods. Derived based on Gaussian process models for regression, the extensions provide distributions for target value estimates and the framework to select the hyperparameters. It is shown that both the proposed methods asymptotically converge to kernel and mutual $k$-nearest neighbor regression method, respectively. The simulation results show that the proposed methods can select proper hyperparameters and are better than or comparable to the former methods for an artificial data set and a real world data set.'
author:
- 'Hyun-Chul Kim[^1]'
bibliography:
- 'paper.bib'
title: 'Bayesian Kernel and Mutual $k$-Nearest Neighbor Regression'
---

kernel regression, Bayesian kernel regression, bandwidth selection, mutual $k$-NN regression, Bayesian mutual $k$-NN regression, Gaussian processes, Bayesian model selection

Introduction
============

In regression analysis, it is analyzed how the response variable $y \in \mathbf{R}$ depends on the value of the observation vector ${\ensuremath{{\rm\bf{x}}}} \in \mathbf{R}^d$ [@opac-b1123996]. The classical approach for regression is a parametric regression approach, where a certain type of structure of the regression function is assumed to be known and can be described by finitely many parameters. A big drawback of the parametric regression is that it cannot approximate the function better than the best one in the assumed parametric structure. This drawback can be avoided by the nonparametric regression approach, where does not assume that the regression function can be determined by finitely many parameters.

One of well-known nonparametric regression methods is kernel regression [@opac-b1123996]. Kernel regression estimator has kernel function with an essential parameter called bandwidth. The performance of kernel regression is known to be largely dependent on the selection of the parameter called bandwidth [@herrmann2000variance]. Methods based on cross-validation for bandwidth selection in kernel regression has been proposed [@hardle1985optimal; @hardle1997multivariate]. It was shown to be asymptotically optimal. A Bayesian approach of averaging over bandwidth rather than selecting bandwidth has been also proposed [@zhang2009bayesian]. They set a leave-one-out likelihood and a prior for bandwidth, and derived the posterior estimate of bandwidth. Rather than selecting a single bandwidth, they averaged the estimator over bandwidth. Metropolis-Hastings algorithm was used for sampling bandwidth to approximate the posterior.

Kernel regression has been widely applied to various kinds of areas including the empirical finance. It was applied to estimate state-densities implicit in financial asset prices [@ait1998nonparametric; @ait2000nonparametric]. It was also applied to estimate call prices and exercise decisions [@broadie2000american], option prices [@aït2001goodness], and the probability of a crash [@fernandes2006financial]. Besides the empirical finance area, kernel regression and its data-adapted extensions have been applied to image processing, restoration, and enhancement [@Takeda07kernelregression]. Bilateral fitlering widely used in image processing is a data-adapted extension of kernel regression which uses image pixel values as inputs of kernel function [@Takeda07kernelregression; @tomasi1998bilateral; @elad2002origin]. Nonlocal-means filtering for super-resolution reconstruction is a data-adapted extension of kernel regression which uses pixel locations and extracted image patches as inputs of kernel function [@protter2009generalizing]. Another data-dependent extension of kernel regression is steering kernel regression which uses pixel locations and local gradients in the image as inputs of kernel function [@Takeda07kernelregression].

Another of well-known nonparametric regression methods is $k$-nearest neighbor ($k$-NN) regression [@opac-b1123996]. Mutual $k$-NN (M$k$NN) regression is a variate of $k$-NN regression based on mutual neighborship rather than one-sided neighborship [@JMLR:v14:guyader13a]. Even though it was not applied to regression task, M$k$NN concept existed even in 1970s. [@gowda1978agglomerative; @gowda1979condensed] has used M$k$NN methods for clustering. More recently, M$k$NN methods have been applied to classification [@liu2010new], outlier detection [@hautamaki2004outlier], object retrieval [@jegou2010accurate], and clustering of interval-valued symbolic patterns [@guru2006clustering]. [@ozaki2011using] used M$k$NN concept to semi-supervised classification of natural language data and showed that the case of using M$k$NN concept consistently outperform the case of using $k$-NN concept. For the regression case, [@JMLR:v14:guyader13a] has argued that M$k$NN regression might be less influenced by hubs which are data points appearing in the nearest neighbor list of many data points.

In this paper we propose a novel method for Bayesian kernel regression and a Bayesian method for mutual $k$-NN regression[[^2]]{}. It provides not only the distribution of target value, but also Bayesian model selection framework for bandwidth in kernel regression and $k$ in mutual $k$-NN regression. We take Laplacian-based covariance matrix and use Gaussian process model. We show that the mean of target values in the proposed methods is asymptotically the estimates of kernel regression and mutual $k$-NN regression, respectively

The paper is organized as follows. In section \[sec:kernel\_MkNN\_reg\] we describe kernel regression and mutual $k$-NN regression. In section \[sec:Bayes\_kernel\_MkNN\] we briefly explain Gaussian process regression with Laplacian-based covariance matrix, and based on that model propose Bayesian kernel regression and Bayesian mutual $k$-NN regression. In the section we also propose the Bayesian model selection methods for both cases. In section \[sec:sim\_res\] we show simulation results for an artificial data set and a real-word data set. In section \[sec:related\_works\] we discuss the works related to the proposed models. Finally a conclusion is drawn.

Kernel and Mutual $k$-Nearest Neighbor Regression {#sec:kernel_MkNN_reg}
=================================================

Kerenel Regression
------------------

Let us assume that we have the set of pairs of the observation vector and response variable $\mathcal{D}_n=\{ ({\ensuremath{{\rm\bf{x}}}}_1, y_1), \ldots, ({\ensuremath{{\rm\bf{x}}}}_n, y_n) \}$, where ${\ensuremath{{\rm\bf{x}}}}_i \in \mathbf{R}^d$ and $y_i \in \mathbf{R}$. Given a new observation vector ${\ensuremath{{\rm\bf{x}}}}$, the kernel estimate is defined as follows. $$\begin{aligned}
m^{\mathrm{ker}}_n({\ensuremath{{\rm\bf{x}}}}) = 
\left\{ \begin{array}{ll}
         \frac{\sum_{i=1}^n y_i k(\frac{{\ensuremath{{\rm\bf{x}}}} - {\ensuremath{{\rm\bf{x}}}}_i}{h_n})}
                {\sum_{i=1}^n k(\frac{{\ensuremath{{\rm\bf{x}}}} - {\ensuremath{{\rm\bf{x}}}}_i}{h_n})}
 & \mbox{if $\sum_{i=1}^n k(\frac{{\ensuremath{{\rm\bf{x}}}} - {\ensuremath{{\rm\bf{x}}}}_i}{h_n}) \neq 0$};\\
        0 & \mbox{otherwise,} \end{array} \right.
\label{eqn:ker_reg}\end{aligned}$$ where $h_n$ is a bandwidth only depending on the sample size $n$ and $k(\cdot)$ is a kernel function which maps $\mathcal{R}^d$ to $ \mathcal{R}^{+}\cup \{ 0 \}$ [@opac-b1123996]. In this paper we use the following kernel function. $$\begin{aligned}
    k({\ensuremath{{\rm\bf{z}}}}) = \exp ( - || {\ensuremath{{\rm\bf{z}}}}  ||^2 ).\end{aligned}$$

In the kernel function $k(\frac{{\ensuremath{{\rm\bf{x}}}}_i - {\ensuremath{{\rm\bf{x}}}}_j}{h_n})$, a single identical bandwidth for all the dimensions is used. In addition to the kernel function with a single bandwidth we use another kernel function, which uses an individual bandwidth for each dimension, $k(({\ensuremath{{\rm\bf{x}}}}_i - {\ensuremath{{\rm\bf{x}}}}_j) {{\ensuremath{{\rm\bf{H}}}}}^{-1})$ where ${\ensuremath{{\rm\bf{H}}}}$ is a $d\times d$ diagonal matrix with the diagonal elements $h_{n,1}, h_{n,2}, \ldots, h_{n,d}$. The former one is called the single bandwidth scheme, and the latter one is called the multiple bandwidth scheme. For the convenience, we also denote $k(\frac{{\ensuremath{{\rm\bf{x}}}}_i - {\ensuremath{{\rm\bf{x}}}}_j}{h_n})$, or $k(({\ensuremath{{\rm\bf{x}}}}_i - {\ensuremath{{\rm\bf{x}}}}_j) {{\ensuremath{{\rm\bf{H}}}}}^{-1})$ by $k({\ensuremath{{\rm\bf{x}}}}_i, {\ensuremath{{\rm\bf{x}}}}_j)$, and we mention whether it is under the single bandwidth scheme or the multiple bandwidth scheme.

Mutual $k$-Nearest Neighbor Regression
--------------------------------------

First, we describe $k$-nearest neighbor regression method [@opac-b1123996]. We denote a reordering of the elements of $\mathcal{D}_n$ by $({\ensuremath{{\rm\bf{x}}}}_{(1,n)}({\ensuremath{{\rm\bf{x}}}}),y_{(1,n)}({\ensuremath{{\rm\bf{x}}}})),\ldots,({\ensuremath{{\rm\bf{x}}}}_{(n,n)}({\ensuremath{{\rm\bf{x}}}}),y_{(n,n)}({\ensuremath{{\rm\bf{x}}}}) )$ according to increasing values of $||{\ensuremath{{\rm\bf{x}}}}_i - {\ensuremath{{\rm\bf{x}}}}||$, Euclidean distance between ${\ensuremath{{\rm\bf{x}}}}_i$ and ${\ensuremath{{\rm\bf{x}}}}$. Then, given ${\ensuremath{{\rm\bf{x}}}} \in \mathbf{R}^d$, the $k$-Nearest Neighbor ($k$-NN) estimate is defined by $$\begin{aligned}
   m^{k\mbox{NNR}}_n({\ensuremath{{\rm\bf{x}}}}) = \frac{1}{k}\sum_{i=1}^k y_{(i,n)}({\ensuremath{{\rm\bf{x}}}}).\end{aligned}$$

Mutual $k$-nearest neighbor regression method [@JMLR:v14:guyader13a] is a method to consider $k$ mutually nearest neighbor data points, rather than just $k$ nearest neighbor data points. Let $\mathcal{N}_k({\ensuremath{{\rm\bf{x}}}})$ be the set of the $k$ nearest neighbors of ${\ensuremath{{\rm\bf{x}}}}$ in $\mathcal{D}_n$, $\mathcal{N}'_k({\ensuremath{{\rm\bf{x}}}}_i)$ the set of $k$ nearest neighbors of ${\ensuremath{{\rm\bf{x}}}}_i$ in $(\mathcal{D}_n \backslash \{{\ensuremath{{\rm\bf{x}}}}_i \}) \cup \{ {\ensuremath{{\rm\bf{x}}}} \}.$ The set of Mutual $k$-Nearest Neighbors (M$k$NNs) of ${\ensuremath{{\rm\bf{x}}}}$ is defined as $$\begin{aligned}
    \mathcal{M}_k({\ensuremath{{\rm\bf{x}}}}) = \{ {\ensuremath{{\rm\bf{x}}}}_i \in \mathcal{N}_k({\ensuremath{{\rm\bf{x}}}}) : {\ensuremath{{\rm\bf{x}}}} \in \mathcal{N}'_k ({\ensuremath{{\rm\bf{x}}}}_i) \}.\end{aligned}$$ Then, the mutual $k$-nearest neighbor regression estimate is defined as $$\begin{aligned}
    m^{\mbox{M$k$NNR}}_n({\ensuremath{{\rm\bf{x}}}}) = \left\{ \begin{array}{ll}
         \frac{1}{M_k({\ensuremath{{\rm\bf{x}}}})}\sum_{i:{\ensuremath{{\rm\bf{x}}}}_i \in \mathcal{M}_k({\ensuremath{{\rm\bf{x}}}})}^k y_i & \mbox{if $M_k({\ensuremath{{\rm\bf{x}}}}) \neq 0$};\\
        0 & \mbox{if $M_k({\ensuremath{{\rm\bf{x}}}}) = 0$}.\end{array} \right.
\label{eqn:mut_nn_reg_est}\end{aligned}$$ where $M_k({\ensuremath{{\rm\bf{x}}}}) = | \mathcal{M}_k({\ensuremath{{\rm\bf{x}}}})| $.

Bayesian Kernel and Mutual $k$-NN Regression via Gaussian Processes {#sec:Bayes_kernel_MkNN}
===================================================================

Gaussian Process Regression
---------------------------

Assume that we have a data set $D$ of data points ${\ensuremath{{\rm\bf{x}}}}_i$ with continuous target values $y_i$: $D = \{({\ensuremath{{\rm\bf{x}}}}_i, y_i)|i = 1, 2,
\ldots, n\}$, $X = \{{\ensuremath{{\rm\bf{x}}}}_i|i = 1, 2, \ldots, n \}$, ${\ensuremath{{\rm\bf{y}}}}=[y_1,
y_2, \ldots, y_n]^T$. We assume that the observations of target values are nosiy, and set $y_i=f({\ensuremath{{\rm\bf{x}}}}_i)+\epsilon_i$, where $f(\cdot)$ is a target function to be estimated and $\epsilon_i \sim \mathcal{N}(0,v_1)$. A function $f(\cdot)$ to be estimated given $D$ is assumed to have Gaussian process prior, which means that any collection of functional values are assumed to be multivariate Gaussian [@Williams95gaussian; @Rasmussen:2005:GPM:1162254].

The prior for the function values ${\ensuremath{{\rm\bf{f}}}}$ $(=[f({\ensuremath{{\rm\bf{x}}}}_1), f({\ensuremath{{\rm\bf{x}}}}_2),  $ $\ldots, f({\ensuremath{{\rm\bf{x}}}}_n) ]^T)$ is assumed to be Gaussian: $$\begin{aligned}
p({\ensuremath{{\rm\bf{f}}}}|X, \Theta_f) = \mathcal{N}({\ensuremath{{\rm\bf{0}}}}, {\ensuremath{{\rm\bf{C}}}}_f).\end{aligned}$$ Then the density function for the target values can be described as follows. $$\begin{aligned}
    p({\ensuremath{{\rm\bf{y}}}}| X, \Theta) &= \mathcal{N}({\ensuremath{{\rm\bf{0}}}}, {\ensuremath{{\rm\bf{C}}}}_f + v_1 {\ensuremath{{\rm\bf{I}}}}) \\
                    &= \mathcal{N}({\ensuremath{{\rm\bf{0}}}}, {\ensuremath{{\rm\bf{C}}}}),\end{aligned}$$ where ${\ensuremath{{\rm\bf{C}}}}$ is a matrix whose elements $C_{ij}$ is a covariance function value $c({\ensuremath{{\rm\bf{x}}}}_i,{\ensuremath{{\rm\bf{x}}}}_j)$ of ${\ensuremath{{\rm\bf{x}}}}_i$, ${\ensuremath{{\rm\bf{x}}}}_j$ and $\Theta$ is the set of hyperparameters in the covariance function.

One of the widely used covariance functions is as follows: $$\begin{aligned}
c({\ensuremath{{\rm\bf{x}}}}_i, {\ensuremath{{\rm\bf{x}}}}_j) =& v_0 \exp \{ -\frac{1}{2}\sum_{m=1}^d l_m
(x_i^m-x_j^m)^2 \} + v_1\delta(i,j), \label{eqn:cov_ftn}\end{aligned}$$ where $x_i^m$ is the $m$th element of ${\ensuremath{{\rm\bf{x}}}}_i$. The hyperparameter $v_0$ specifies the overall vertical scale of variation of the target values, $v_1$ the noise variance of the target values, and $l_m$ the (inverse) lengthscale for feature dimension $m$. The covariance function described in Eq (\[eqn:cov\_ftn\]) enables GPR to estimate a nonlinear function. $v_1$ makes GPR robust to noise, and the optimized $l_1, \ldots, l_d$ can be used to determine relevant features.

It can be shown that GPR provides the following distribution of target value $f_{\mathrm{new}}(=f({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}))$ given a test data ${\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}$: $$\begin{aligned}
 p(f_{\mathrm{new}}|{\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}, D,\Theta)=\mathcal{N}({\ensuremath{{\rm\bf{k}}}}^T {\ensuremath{{\rm\bf{C}}}}^{-1} {\ensuremath{{\rm\bf{f}}}}, \kappa - {\ensuremath{{\rm\bf{k}}}}^T {\ensuremath{{\rm\bf{C}}}}^{-1} {\ensuremath{{\rm\bf{k}}}}),
 \label{eqn:result}\end{aligned}$$ where ${\ensuremath{{\rm\bf{k}}}} = [c({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}},{\ensuremath{{\rm\bf{x}}}}_1) \ldots c({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}},{\ensuremath{{\rm\bf{x}}}}_n)]^T$, $\kappa=c({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}},{\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}})$. The variance of the target value $f_{\mathrm{new}}$ is related to the degree of its uncertainty.

Since the set of hyperparameters $\Theta(=\{v_0,v_1,l_1,\ldots,l_d\})$ controls the model complexity, it is very important to select the most proper $\Theta$. The proper $\Theta$ can be obtained by maximizing the marginal likelihood $p({\ensuremath{{\rm\bf{y}}}}|X,\Theta)$ [@Williams95gaussian; @Gibbs97efficientimplementation; @Rasmussen:2005:GPM:1162254]. $$\begin{aligned}
    \log p({\ensuremath{{\rm\bf{y}}}} | {\ensuremath{{\rm\bf{X}}}}, \Theta) = -\frac{1}{2}{\ensuremath{{\rm\bf{y}}}}^\top {\ensuremath{{\rm\bf{C}}}}^{-1} {\ensuremath{{\rm\bf{y}}}} - \frac{1}{2} \log | {\ensuremath{{\rm\bf{C}}}} | - \frac{N  T}{2}  \log 2 \pi.
\label{eqn:loglikeli}\end{aligned}$$ Rather than choosing a single set of hyperparameters by optimization, we can average over the hyperparameters with MCMC methods [@Williams95gaussian; @Neal1997RC_GP].

Laplacian-based Covariance Matrix
---------------------------------

The combinatorial Laplacian ${\ensuremath{{\rm\bf{L}}}}$ is defined as follows. $$\begin{aligned}
{\ensuremath{{\rm\bf{L}}}} = {\ensuremath{{\rm\bf{D}}}} - {\ensuremath{{\rm\bf{W}}}},\label{eq:Delta}\end{aligned}$$where ${\ensuremath{{\rm\bf{W}}}}$ is an $N\times N$ edge-weight matrix with the edge weight between two points ${\ensuremath{{\rm\bf{x}}}}_i$,${\ensuremath{{\rm\bf{x}}}}_j$ given as $w_{ij}(=w({\ensuremath{{\rm\bf{x}}}}_i,{\ensuremath{{\rm\bf{x}}}}_j))$ and ${\ensuremath{{\rm\bf{D}}}} =
\mathrm{diag}(d_1,...,d_N)$ is a diagonal matrix with diagonal entries $d_i=\sum_{j} w_{ij}$.

Similarly to [@zhu2003semi], to avoid the singularity we use Laplacian-based covariance matrix as $$\begin{aligned}
   {\ensuremath{{\rm\bf{C}}}} =&  ({\ensuremath{{\rm\bf{L}}}} + \sigma^2 {\ensuremath{{\rm\bf{I}}}})^{-1} = \tilde{{\ensuremath{{\rm\bf{C}}}}}^{-1}.\end{aligned}$$ Then, we have Gaussian process prior as follows. $$\begin{aligned}
    p({\ensuremath{{\rm\bf{y}}}}| X, \Theta) &= \mathcal{N} ({\ensuremath{{\rm\bf{0}}}}, {\ensuremath{{\rm\bf{C}}}}),\end{aligned}$$ The predictive distiribution for $y_{\mathrm{new}}$ is as follows. $$\begin{aligned}
    p(y_{\mathrm{new}}| {\ensuremath{{\rm\bf{y}}}}, X, {\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}, \Theta) &= \frac{p({\ensuremath{{\rm\bf{y}}}}_{\mathrm{new}}| X, {\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}, \Theta)}{p({\ensuremath{{\rm\bf{y}}}}| X, {\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}, \Theta)}
    \nonumber \\
    &= \mathcal{N}_{{\ensuremath{{\rm\bf{y}}}}_{\mathrm{new}}} ({\ensuremath{{\rm\bf{0}}}}, {\ensuremath{{\rm\bf{C}}}}_{\mathrm{new}}) / \mathcal{N}_{{\ensuremath{{\rm\bf{y}}}}} ({\ensuremath{{\rm\bf{0}}}}, {\ensuremath{{\rm\bf{C}}}}) \nonumber \\
                &= \mathcal{N}_{{\ensuremath{{\rm\bf{y}}}}_{\mathrm{new}}} ({\ensuremath{{\rm\bf{0}}}}, \tilde{{\ensuremath{{\rm\bf{C}}}}}_{\mathrm{new}}^{-1}) /  \mathcal{N}_{{\ensuremath{{\rm\bf{y}}}}} ({\ensuremath{{\rm\bf{0}}}}, \tilde{ {\ensuremath{{\rm\bf{C}}}}}^{-1}) \nonumber \\
                &\propto \frac{\exp(- \frac{1}{2} [{\ensuremath{{\rm\bf{y}}}}^T y_{\mathrm{new}}] 
                \tilde{{\ensuremath{{\rm\bf{C}}}}}_{\mathrm{new}} [{\ensuremath{{\rm\bf{y}}}}^T y_{\mathrm{new}}]^T )}
                {\exp(- \frac{1}{2} {\ensuremath{{\rm\bf{y}}}}^T
                \tilde{{\ensuremath{{\rm\bf{C}}}}} {\ensuremath{{\rm\bf{y}}}} )} \nonumber \\
                &\propto \exp(- \frac{1}{2}\tilde{\kappa} y_{\mathrm{new}}^2 - \tilde{{\ensuremath{{\rm\bf{k}}}}}^T {\ensuremath{{\rm\bf{y}}}} y_{\mathrm{new}}) \nonumber \\
                &\propto  \mathcal{N}(-\frac{1}{\tilde{\kappa}} \tilde{{\ensuremath{{\rm\bf{k}}}}}^T {\ensuremath{{\rm\bf{y}}}},\frac{1}{\tilde{\kappa}}),\end{aligned}$$ where $$\begin{aligned}
\tilde{{\ensuremath{{\rm\bf{C}}}}}_{\mathrm{new}}&=\left[\begin{array}{ll}
\mathbf{\tilde{C}} & \tilde{{\ensuremath{{\rm\bf{k}}}}} \\
   \tilde{{\ensuremath{{\rm\bf{k}}}}}^T & \tilde{\kappa}
   \end{array} \right]
=\mathbf{C}_{\mathrm{new}}^{-1},  \\
\label{aaa}
   \tilde{\kappa} &= \sum_{i=1}^N w({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}},{\ensuremath{{\rm\bf{x}}}}_i) + \sigma^2, \\
   \tilde{ {\ensuremath{{\rm\bf{k}}}}}^T &= - [ w({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}, {\ensuremath{{\rm\bf{x}}}}_1), w({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}, {\ensuremath{{\rm\bf{x}}}}_2), \ldots, w({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}, {\ensuremath{{\rm\bf{x}}}}_N) ].\end{aligned}$$ The mean and variance of $y_{\mathrm{new}}$ is represented as $$\begin{aligned}
   \mu_{ y_{\mathrm{new}}}&= -\frac{1}{\tilde{\kappa}} \tilde{{\ensuremath{{\rm\bf{k}}}}}^T {\ensuremath{{\rm\bf{y}}}}_L 
                            =\frac{\sum_{i=1}^N w ({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}, {\ensuremath{{\rm\bf{x}}}}_i) y_i }{\sum_{i=1}^N w ({\ensuremath{{\rm\bf{x}}}}_\mathrm{new}, {\ensuremath{{\rm\bf{x}}}}_i)+\sigma^2}, \label{eq:transd_lap_mean} \\
   \sigma^2_{ y_{\mathrm{new}}}&= \frac{1}{\tilde{\kappa}}
                            =\frac{1}{\sum_{i=1}^N w({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}},{\ensuremath{{\rm\bf{x}}}}_i) + \sigma^2}. \label{eq:transd_lap_std}\end{aligned}$$ Eq (\[eq:transd\_lap\_mean\]) links Gaussian process regression with Laplacian-based covariance matrix to kernel and mutual $k$-NN regression, which will be described below.

Bayesian Kernel Regression
--------------------------

If we set $w_{ij}=w_{\mathrm{ker}} ({\ensuremath{{\rm\bf{x}}}}_i, {\ensuremath{{\rm\bf{x}}}}_j)=\sigma_0 k({\ensuremath{{\rm\bf{x}}}}_i, {\ensuremath{{\rm\bf{x}}}}_j)$, where $k(\cdot)$ is a kernel function in Eq (\[eqn:ker\_reg\]) and $\sigma_0>0$, we get the following theorem for the validity of the covariance matrix for Gaussian processes.

Covairance matrix $\tilde{{\ensuremath{{\rm\bf{C}}}}}$ with $w_{ij}=w_{\mathrm{ker}} ({\ensuremath{{\rm\bf{x}}}}_i, {\ensuremath{{\rm\bf{x}}}}_j)$ is valid for Gaussian processes if $\sigma^2>0$. \[thm:ker\_cov\_valid\]

\(1) Since Laplacian matrix ${\ensuremath{{\rm\bf{L}}}}(={\ensuremath{{\rm\bf{D}}}} - {\ensuremath{{\rm\bf{W}}}})$ is positive semidefinite [@merris1994laplacian], for $\sigma^2>0$ $\tilde{{\ensuremath{{\rm\bf{C}}}}}(={\ensuremath{{\rm\bf{L}}}} + \sigma^2 {\ensuremath{{\rm\bf{I}}}})$ is positive definite. So $\tilde{{\ensuremath{{\rm\bf{C}}}}}$ is positive definite.\
(2) Since $\tilde{{\ensuremath{{\rm\bf{C}}}}}^T = ({\ensuremath{{\rm\bf{D}}}} - {\ensuremath{{\rm\bf{W}}}} + \sigma^2 {\ensuremath{{\rm\bf{I}}}})^T
        = {\ensuremath{{\rm\bf{D}}}}^T - {\ensuremath{{\rm\bf{W}}}}^T + \sigma^2 {\ensuremath{{\rm\bf{I}}}}^T = {\ensuremath{{\rm\bf{D}}}} - {\ensuremath{{\rm\bf{W}}}} + \sigma^2 {\ensuremath{{\rm\bf{I}}}}
            = \tilde{{\ensuremath{{\rm\bf{C}}}}}$, $\tilde{{\ensuremath{{\rm\bf{C}}}}}$ is symmetric.\
From (1) & (2), by Theorem 7.5 in [@stefanica2014linear] $\tilde{{\ensuremath{{\rm\bf{C}}}}}$ is a valid covariance matrix. QED.

By applying Eq (\[eq:transd\_lap\_mean\]), Bayesian kernel regression estimate for ${\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}$ is defined as follows: $$\begin{aligned}
m^{\mathrm{Bker}}_n({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}) = 
\mu_{y_{\mathrm{new}},\mathrm{ker}},\end{aligned}$$ where $$\begin{aligned}
    \mu_{y_{\mathrm{new}},\mathrm{ker}} &=  \frac{\sum_{i=1}^N w_\mathrm{ker} ({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}, {\ensuremath{{\rm\bf{x}}}}_i) y_i }{\sum_{i=1}^N w_{\mathrm{ker}} ({\ensuremath{{\rm\bf{x}}}}_\mathrm{new}, {\ensuremath{{\rm\bf{x}}}}_i)+\sigma^2} \nonumber \\
        &= \frac{\sum_{i=1}^N k ({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}, {\ensuremath{{\rm\bf{x}}}}_i) y_i }{\sum_{i=1}^N  k ({\ensuremath{{\rm\bf{x}}}}_\mathrm{new}, {\ensuremath{{\rm\bf{x}}}}_i)+\sigma^2/\sigma_0}.
\label{eqn:Bayes_ker_reg}\end{aligned}$$

The following theorem shows that Bayesian kernel regression introduced above asymptotically converges to the traditional kernel regression.

$\mu_{y_{\mathrm{new}},{\mathrm{ker}}} (=-\frac{1}{\tilde{\kappa}} \tilde{{\ensuremath{{\rm\bf{k}}}}}^T {\ensuremath{{\rm\bf{y}}}})$ converges to kernel regression as $\sigma^2/\sigma_0$ approaches $0$.

In case $\sum_{i=1}^N  k ({\ensuremath{{\rm\bf{x}}}}_\mathrm{new}, {\ensuremath{{\rm\bf{x}}}}_i) = 0$, $ k ({\ensuremath{{\rm\bf{x}}}}_\mathrm{new}, {\ensuremath{{\rm\bf{x}}}}_i) = 0$ for all $i$. So it is trivial by Eq (\[eqn:ker\_reg\]) and (\[eqn:Bayes\_ker\_reg\]).

Otherwise, take a small positive $\epsilon< m^{\mathrm{ker}}_n({\ensuremath{{\rm\bf{x}}}})$. Set $\delta = \{ \sum_{i=1}^N k ({\ensuremath{{\rm\bf{x}}}}, {\ensuremath{{\rm\bf{x}}}}_i) \} / \{ \frac{m^{\mathrm{ker}}_n({\ensuremath{{\rm\bf{x}}}})}{\epsilon} -1 \}$. Then, if $||\sigma^2/\sigma_0||<\delta $, $ || \mu_{{\ensuremath{{\rm\bf{f}}}}_U,\mathrm{M}{\it k}\mathrm{NN}} - m^{\mathrm{ker}}_n({\ensuremath{{\rm\bf{x}}}})  || < \epsilon $. By the $(\epsilon,\delta)$ definition of the limit of a function, we get the statement in the theorem. QED.

Bayesian Mutual $k$-NN Regression
---------------------------------

Gaussian processes with Laplacian-based covariance matrix can be associated with mutual $k$-NN regression, by replacing $w_{ij}(=w({\ensuremath{{\rm\bf{x}}}}_i,{\ensuremath{{\rm\bf{x}}}}_j))$ with the function $$\begin{aligned}
  w_{\mathrm{M}\it{k}\mathrm{NN}}({\ensuremath{{\rm\bf{x}}}}_i,{\ensuremath{{\rm\bf{x}}}}_j) =  \sigma_0 \delta_{{\ensuremath{{\rm\bf{x}}}}_j \sim_k {\ensuremath{{\rm\bf{x}}}}_i} \cdot \delta_{{\ensuremath{{\rm\bf{x}}}}_i \sim_k {\ensuremath{{\rm\bf{x}}}}_j},\end{aligned}$$ where the relation $\sim_k$ is defined as $$\begin{aligned}
    {\ensuremath{{\rm\bf{x}}}}_i \sim_k {\ensuremath{{\rm\bf{x}}}}_j = \left\{ \begin{array}{ll}
 T  & \mbox{if $j \neq i$ and ${\ensuremath{{\rm\bf{x}}}}_j$ is a $k$-nearest neighbor} \\
    & \mbox{of ${\ensuremath{{\rm\bf{x}}}}_i$}; \\
   F   & \mbox{otherwise,} \end{array} \right.  \end{aligned}$$

Similarly to Bayesian kernel regression, we get the following theorem related to the validity of Laplacian-based covariance matrix for Gaussian processes.

Covairance matrix $\tilde{{\ensuremath{{\rm\bf{C}}}}}$ with $w_{ij}(=w_{\mathrm{M}\it{k}\mathrm{NN}}({\ensuremath{{\rm\bf{x}}}}_i,{\ensuremath{{\rm\bf{x}}}}_j))$ is valid for Gaussian processes if $\sigma^2>0$.

Similarly to the proof of Theorem \[thm:ker\_cov\_valid\], it can be shown that $\tilde{{\ensuremath{{\rm\bf{C}}}}}$ is positive definite and symmetric.

By applying Eq (\[eq:transd\_lap\_mean\]) like in Bayesian kernel regression, Bayesian mutual $k$-NN regression estimate for a given data ${\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}$ is defined as follows. $$\begin{aligned}
   m^{\mbox{BM}k\mbox{NNR}}_n({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}) =      \mu_{f_{\mathrm{new}},\mathrm{M}\it{k}\mathrm{NN}},\end{aligned}$$ where $$\begin{aligned}
   \mu_{f_{\mathrm{new}},\mathrm{M}\it{k}\mathrm{NN}}
=&\frac{\sum_{i=1}^N w_{\mathrm{M}\it{k}\mathrm{NN}} ({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}, {\ensuremath{{\rm\bf{x}}}}_i) y_i }
{\sum_{i=1}^N w_{\mathrm{M}\it{k}\mathrm{NN}} ({\ensuremath{{\rm\bf{x}}}}_{\mathrm{new}}, {\ensuremath{{\rm\bf{x}}}}_i)+\sigma^2} \nonumber \\
 =&\frac{\sum_{i=1}^N  \delta_{{\ensuremath{{\rm\bf{x}}}}_j \sim_k {\ensuremath{{\rm\bf{x}}}}_i} \cdot \delta_{{\ensuremath{{\rm\bf{x}}}}_i \sim_k {\ensuremath{{\rm\bf{x}}}}_j} y_i }
{\sum_{i=1}^N \delta_{{\ensuremath{{\rm\bf{x}}}}_j \sim_k {\ensuremath{{\rm\bf{x}}}}_i} \cdot \delta_{{\ensuremath{{\rm\bf{x}}}}_i \sim_k {\ensuremath{{\rm\bf{x}}}}_j} +\sigma^2/\sigma_0 }.
\label{eqn:Bayes_mut_knn_reg}\end{aligned}$$

The following theorem shows that Bayesian mutual $k$-nearest neighbor regression introduced above asymptotically converges to the traditional mutual $k$-NN regression.

$\mu_{f_{\mathrm{new}},{\mathrm{M}{\it k}\mathrm{NN}}} (=-\frac{1}{\tilde{\kappa}} \tilde{{\ensuremath{{\rm\bf{k}}}}}^T {\ensuremath{{\rm\bf{y}}}})$ converges to mutual $k$-$\mathrm{NN}$ regression as $\sigma^2/\sigma_0$ approaches $0$.

In case $\sum_{i=1}^N \delta_{{\ensuremath{{\rm\bf{x}}}}_j \sim_k {\ensuremath{{\rm\bf{x}}}}_i} \cdot \delta_{{\ensuremath{{\rm\bf{x}}}}_i \sim_k {\ensuremath{{\rm\bf{x}}}}_j} = 0$, $\delta_{{\ensuremath{{\rm\bf{x}}}}_j \sim_k {\ensuremath{{\rm\bf{x}}}}_i} \cdot \delta_{{\ensuremath{{\rm\bf{x}}}}_i \sim_k {\ensuremath{{\rm\bf{x}}}}_j} = 0$ for all $i$. So it is trivial by Eq (\[eqn:mut\_nn\_reg\_est\]) and (\[eqn:Bayes\_mut\_knn\_reg\]).

Otherwise, take a small positive $\epsilon< m^{\mbox{M}k\mbox{NNR}}_n({\ensuremath{{\rm\bf{x}}}})$.\
Set $\delta = \{ \sum_{i=1}^N  \delta_{{\ensuremath{{\rm\bf{x}}}}_j \sim_k {\ensuremath{{\rm\bf{x}}}}_i} \cdot \delta_{{\ensuremath{{\rm\bf{x}}}}_i \sim_k {\ensuremath{{\rm\bf{x}}}}_j}  \} / \{ \frac{m^{\mbox{M}k\mbox{NNR}}_n({\ensuremath{{\rm\bf{x}}}})}{\epsilon} -1 \}$. Then, if $||\sigma^2/\sigma_0||<\delta $, $ || \mu_{{\ensuremath{{\rm\bf{f}}}}_U,\mathrm{M}{\it k}\mathrm{NN}} - m^{\mbox{M}k\mbox{NNR}}_n({\ensuremath{{\rm\bf{x}}}})  || < \epsilon $. By the $(\epsilon,\delta)$ definition of the limit of a function, we get the statement in the theorem. QED.

Hyperparameter Selection
------------------------

There is a set of hyperparameters that should be selected in both the proposed methods. For Bayesian kernel regression, the set of hyperparameters is $\Theta=\{ h_n, \sigma_0, \sigma  \}$ for the single bandwidth scheme, or $\Theta=\{ h_{n,1}, h_{n,2}, \ldots, h_{n,d}, \sigma_0, \sigma  \}$ for the multiple bandwidth scheme. For Bayesian $k$-NN regression, the set of hyperparameters is $\Theta=\{ k, \sigma_0, \sigma  \}$, where $k$ is a interger greater than $0$. These sets of hyperparameters can be selected through the Bayesian evidence framework by maximizing the log of the marginal likelihood [@Gibbs97efficientimplementation] as follows. $$\begin{aligned}
    \Theta^*
        =& \mbox{argmax}_{\Theta} \mathcal{L}(\Theta),\end{aligned}$$ where $$\begin{aligned}
\mathcal{L}(\Theta) =& \log p({\ensuremath{{\rm\bf{y}}}}|\Theta) \\
            =& \log \{ |2\pi {\ensuremath{{\rm\bf{C}}}}|^{-\frac{1}{2}} \exp (-\frac{1}{2}{\ensuremath{{\rm\bf{y}}}}^T {\ensuremath{{\rm\bf{C}}}}^{-1} {\ensuremath{{\rm\bf{y}}}}) \} \\
            =& \frac{1}{2} \log |\tilde{{\ensuremath{{\rm\bf{C}}}}}| - \frac{1}{2} {\ensuremath{{\rm\bf{y}}}}^T \tilde{{\ensuremath{{\rm\bf{C}}}}} {\ensuremath{{\rm\bf{y}}}} - \frac{N}{2} \log 2\pi, \label{eqn:bkr_log_ev}\end{aligned}$$ where $\tilde{{\ensuremath{{\rm\bf{C}}}}}={\ensuremath{{\rm\bf{L}}}} + \sigma^2 {\ensuremath{{\rm\bf{I}}}}$.

The discrete hyperparameters (e.g., $k$) can be selected based on the value of $\mathcal{L}$ as $$\begin{aligned}
      K^*   &= \mbox{argmax}_k  \mathcal{L}(\{k,\sigma,\sigma_0\}). \end{aligned}$$ For the continuous hyperparameters (e.g., $\sigma$, $\sigma_0$, $h_n$), to optimize $\mathcal{L}$ with respect to $\Theta$ we can use the derivative $$\begin{aligned}
\frac{\partial \mathcal{L}}{\partial \theta} = \frac{1}{2}
\mathrm{trace}(\tilde{{\ensuremath{{\rm\bf{C}}}}}^{-1} \frac{\partial \tilde{{\ensuremath{{\rm\bf{C}}}}}}{\partial
\theta}) - \frac{1}{2}{\ensuremath{{\rm\bf{y}}}}^T \frac{\partial \tilde{{\ensuremath{{\rm\bf{C}}}}}}{\partial \theta} {\ensuremath{{\rm\bf{y}}}},\end{aligned}$$ where $\tilde{{\ensuremath{{\rm\bf{C}}}}}={\ensuremath{{\rm\bf{L}}}} + \sigma^2 {\ensuremath{{\rm\bf{I}}}}$.

On the other hand, the posterior distributions of the hyperparameters given the data can be inferred by the Bayesian method via Markov Chain Monte Carlo methods similarly to [@Neal1997RC_GP; @Williams95gaussian]. And the regression estimate can be averaged over the hyperparameters rather than obtained by one fixed set of hyperparameters. This would produce better results but cost more computational power. This approach has not been taken in this paper

Simulation Results {#sec:sim_res}
==================

We did the simulations for the proposed methods to observe how the methods work and show their usefulness. The simulations for both methods were performed for both an artificial data set and a real world data set.

Simulation Results for Bayesian Kernel Regression
-------------------------------------------------

First, we did the simulations for Bayesian kernel regression. To generate an artificial data set, we used the equation $\mathrm{sinc}(x) = \frac{\sin(\pi x)}{\pi x}$ for the sinc function. We took the points equally spaced with the interval 0.2 between -5 and 5. We made up the training set with those points as inputs and with the sinc function values at those points as target values. And we took the points equally spaced with the interval 0.1 between -5.01 and 5. We made up the test set with those points as inputs and with the sinc function values at those points as target values. We call this training and test data set as the sinc data set I. We made up another sinc data set (called the sinc data set II), which has the points equally spaced with the interval 0.5 (rather than 0.2) between -5 and 5 as input points of the training set and the same test set as in the sinc data set I.

For the sinc data set I and II, we applied the proposed Bayesian kernel regression. For the comparison we also applied the traditional kernel regression with the bandwidth selected by leave-one-out cross-validation[^3] [@hardle1985optimal]. (We call this bandwidth as CV bandwidth.) In addition we applied kernel regression with the bandwidth chosen in the proposed Bayesian kernel regression. (We call this bandwidth as B bandwidth.) For both data sets we tried the simulations repeatedly with different initial values for $\sigma_0,\sigma$, and found that one of the lowest marginal likelihoods is reached with the initial value 100, 1.

Figure \[fig.1\] shows the results of the three methods with the training data points for the sinc data set I. Apparently the proposed Bayesian kernel regression methods perform best. In Figure \[fig.2\] we can see how the hyperparameter can be selected in Bayesian kernel regression for the sinc data set I. They were plotted in two levels of scales to show that the selection of a specific value as well as averaging is reasonable for the hyperparameter because the evidence is highly peaked near the optimum as seen in Figure 2(b).

![The results of kernel regression with CV bandwidth, Bayesian kernel regression, and kernel regression with B bandwidth for the sinc data set I[]{data-label="fig.1"}](sinc_result_Bkr.eps){width="8cm"}

  ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
   ![ The use of evidence for hyperparameter selection in Bayesian kernel regression for the sinc data set I. Points x mean log evidence or evidence at some points. Circle points represent the optimum: (a) bandwidth vs. log evidence (b) bandwidth vs. evidence[]{data-label="fig.2"}](sinc_result_Bkr_logev.eps "fig:"){width="8cm"}
                                                                                                                                                                    (a)
    ![ The use of evidence for hyperparameter selection in Bayesian kernel regression for the sinc data set I. Points x mean log evidence or evidence at some points. Circle points represent the optimum: (a) bandwidth vs. log evidence (b) bandwidth vs. evidence[]{data-label="fig.2"}](sinc_result_Bkr_ev.eps "fig:"){width="8cm"}
                                                                                                                                                                    (b)
  ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Table \[table:Bkr\_sinc\] shows the performance of the methods applied to the sinc data set I and II. The performance was measured in terms of mean squared error between true target values and regression estimates. Bayesian kernel regression was better than kernel regression with CV bandwidth or B bandwidth.

  ----------------------------------------------- -------------------- -------------------
                      Methods                       sinc data set I     sinc data set II
   \[0.5ex\] Kernel regression with CV bandwidth       9.1840e-04           0.0060448
            Bayesian kernel regression             [**3.5371e-05**]{}   [**0.0012617**]{}
        Kernel regression with B bandwidth             0.0018597            0.018771
                      \[1ex\]                                          
  ----------------------------------------------- -------------------- -------------------

  :  Mean squared errors in kernel regression with CV bandwidth, Bayesian kernel regression, and kernel regression with B bandwidth for the sinc data set I and II

\[table:Bkr\_sinc\]

Next, we applied the three methods including the proposed method to a real world data set. We used the yacht hydrodynamics data set [@Lichman:2013] for the evaluation. The data set is related to prediction of residuary resistance of sailing yachts at the initial design stage, which are very important for evaluating the performance of the ship and for estimating the required propulsive power. It has 308 data points, each of which contains 6 inputs related to hull geometry coefficients and the Froude number, and a target value of residuary resistance per unit weight of displacement. Input data were normalized before they were used for the evaluation of all three methods.

Since the data set is 6-dimensional, we have two choices for bandwidth configuration. A single identical bandwidth can be used for all the dimensions (i.e. single bandwidth), or an individual bandwidth can be used for each dimension (i.e. multiple bandwidths). Bandwidth selections by both single and multiple bandwidth schemes were tried. The leave-one-out cross-validation for multiple bandwidths in the traditional kernel regression[^4] is expensive, but it was tried for comparison.

Table \[table:Bkr\_Yacht\] shows the performance of the three methods with the single bandwidth scheme. The 10 fold cross-validation was done for the performance evaluation. In the table means and standard deviations for the 10 fold cross-validation are shown. To speed up the optimisation and avoid unreasonable local minima we use the heuristic scheme to fix $\sigma_0$ to 1 and $\sigma$ to 0.0000001 and to update the bandwidth only. Bayesian kernel regression and kernel regression with B bandwidth performed better than kernel regression with CV bandwidth.

  ------------------------------------------------ -----------------------------
                      Methods                         MSE ($\mu \pm \sigma$)
   \[0.5ex\] Kernel regressison with CV bandwidth       33.540 $\pm$ 23.507
             Bayesian kernel regression             [**33.529 $\pm$ 23.504**]{}
         Kernel regression with B bandwidth         [**33.529 $\pm$ 23.504**]{}
                      \[1ex\]                      
  ------------------------------------------------ -----------------------------

  :  Means and standard deviations of mean squared errors in kernel regression with CV bandwidth, Bayesian kernel regression, and kernel regression with B bandwidth, in the single bandwidth scheme for the yacht hydrodynamics data set

\[table:Bkr\_Yacht\]

Table \[table:Bkr\_m\_Yacht\] shows the performance of the three methods with the multiple bandwidth scheme. As in the single bandwidth scheme, the 10 fold cross-validation was done for the performance evaluation. In the table means and standard deviations for the cross-validation are shown. We fix $\sigma_0$, $\sigma$ to the same values as in the single bandwidth scheme and update the bandwidths only. Comparing with the results in Table \[table:Bkr\_Yacht\] it is clear that the multiple bandwidth scheme is far better than the single bandwidth scheme. Bayesian kernel regression and kernel regression with B bandwidth performed better than kernel regression with CV bandwidth.

  ------------------------------------------------ -----------------------------
                      Methods                         MSE ($\mu \pm \sigma$)
   \[0.5ex\] Kernel regressison with CV bandwidth       1.3266 $\pm$ 1.1723
             Bayesian kernel regression             [**1.1436 $\pm$ 0.7894**]{}
         Kernel regression with B bandwidth         [**1.1436 $\pm$ 0.7894**]{}
                      \[1ex\]                      
  ------------------------------------------------ -----------------------------

  :  Means and standard deviations of mean squared errors in kernel regression with CV bandwidth, Bayesian kernel regression, and kernel regression with B bandwidth in the multiple bandwidth scheme for the yacht hydrodynamics data set

\[table:Bkr\_m\_Yacht\]

Simulation Results for Bayesian Mutual $k$-NN Regression
--------------------------------------------------------

Second, we did the simulations for Bayesian mutual $k$-NN regression. We used the sinc data set I and II like in the simulation for Bayesian kernel regression. For both the data sets, we applied the proposed Bayesian mutual $k$-NN regression. For the comparison we also applied the traditional $k$-NN regression and mutual $k$-NN regression where both $k$’s were selected by leave-one-out cross-validation. (We call this $k$ as CV $k$.) In addition we applied mutual $k$-NN regression with $k$ chosen in the proposed Bayesian mutual $k$-NN regression. (We call this $k$ as B $k$.) For both data sets we tried the simulation repeatedly with different initial values for $\sigma_0,\sigma$, and found that one of the lowest marginal likelihoods is reached with the initial value 300, 3.

Figure \[fig.3\] shows the results of Bayesian mutual $k$-NN regression and mutual $k$-NN regression with CV $k$ with the training data points for the sinc data set I. Figure \[fig.4\] shows how the leave-one-out cross-validation method works for mutual $k$-NN regression for the sinc data set I. In Figure \[fig.5\] we can see how the hyperparameter $k$ can be selected in Bayesian mutual $k$-NN regression for the sinc data set I. Marginal likelihoods were plotted in two kinds of scales to show that the selection of a specific value as well as averaging is reasonable for the hyperparameter $k$ because the evidence is highly peaked near the optimum as seen in Figure 5(b).

![The results of mutual $k$-NN with CV $k$, Bayesian mutual $k$-NN regression, mutual $k$-NN with B $k$ for sinc data set I[]{data-label="fig.3"}](BkNN_sinc.eps){width="8cm"}

![The leave-one-out cros-validation to select $k$ for mutual $k$-NN regression for sinc data set I[]{data-label="fig.4"}](knn_loocv_err_sinc.eps){width="8cm"}

  -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
   ![The use of evidence for the selection of $k$ in Bayesian mutual $k$-NN regression for the sinc data set I. Points x mean log evidence or evidence at some $k$. Circle points represent the optimum: (a) $k$ vs. log evidence, (b) $k$ vs. evidence []{data-label="fig.5"}](BkNN_logev_sinc.eps "fig:"){width="8cm"}
                                                                                                                                                            (a)
    ![The use of evidence for the selection of $k$ in Bayesian mutual $k$-NN regression for the sinc data set I. Points x mean log evidence or evidence at some $k$. Circle points represent the optimum: (a) $k$ vs. log evidence, (b) $k$ vs. evidence []{data-label="fig.5"}](BkNN_ev_sinc.eps "fig:"){width="8cm"}
                                                                                                                                                            (b)
  -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Table \[table:BMkNN\_sinc\] shows the performance of the methods applied to the sinc data set I and II. The performance was measured in terms of mean squared error between true target values and regression estimates. For the sinc data set I, Bayesian $k$-NN regression was a little bit worse than mutual $k$-NN regression with CV $k$, and mutual $k$-NN regression with B $k$ was as good as mutual $k$-NN regression with CV $k$. For the sinc data set II, Bayesian $k$-NN regression was the best.

  ----------------------------------------- ------------------- -------------------
                   Methods                    sinc data set I    sinc data set II
   \[0.5ex\] $k$-NN regression with CV $k$       0.0012619           0.0069711
    Mutual $k$-NN regression with CV $k$     [**0.0012573**]{}       0.0069573
         Bayesian $k$-NN regression              0.0012587       [**0.0059893**]{}
         $k$-NN regression with B $k$        [**0.0012573**]{}       0.0060497
                   \[1ex\]                                      
  ----------------------------------------- ------------------- -------------------

  :  Mean squared errors in $k$-NN with CV $k$, mutual $k$-NN with CV $k$, Bayesian mutual $k$-NN regression, mutual $k$-NN with B $k$ for the sinc data set I and the sinc data set II

\[table:BMkNN\_sinc\]

Next, we applied the four methods including the proposed method to a real world data set. Similarly to Bayesian kernel regression, we used the yacht hydrodynamics data set for the evaluation. Table \[table:BMkNN\_Yacht\] shows the performance of the three methods. The 10 fold cross-validation was done for the performance evaluation. In the table means and standard deviations for the cross-validation are shown. The hyperparamters $\sigma_0,\sigma$ was fixed to 0.1, 0.00001, and only $k$ was selected. All the three methods showed the same performance with the same $k(=2)$ selected.

  ----------------------------------------- ------------------------------
                   Methods                             MSE avg.
   \[0.5ex\] $k$-NN regression with CV $k$   [**39.127 $\pm$ 23.004** ]{}
    Mutual $k$-NN regression with CV $k$     [**39.127 $\pm$ 23.004** ]{}
      Bayesian mutual $k$-NN regression      [**39.127 $\pm$ 23.004** ]{}
     Mutual $k$-NN regression with B $k$     [**39.127 $\pm$ 23.004** ]{}
                   \[1ex\]                  
  ----------------------------------------- ------------------------------

  :  Means and standard deviations of mean squared errors in $k$-NN regression with CV $k$,mutual $k$-NN regression with CV $k$, Bayesian mutual $k$-NN regression, and mutual $k$-NN regression with B $k$ for the yacht hydrodynamics data set

\[table:BMkNN\_Yacht\]

Related Works {#sec:related_works}
=============

We discuss the former works related to the proposed methods in this paper. [@kim2013transductive] has proposed a transductive regression method with Gaussian processes and applied it to object pose estimation. They has defined Laplacian kernel similarly to Laplaican-based covariance matrix in this paper and proposed a Bayesian method to select the hyperparameters.

[@verbeek2006gaussian] used Gaussian field and Laplacians with a kernel related to $k$-NN and locally linear embedding. They applied their methods to facial pose estimation and object correspondence learning. They proposed active learning method based on entropy minimization, and a model selection scheme based on maximum likelihood.

[@zhu2003semi] extended Gaussian field to Gaussian processes for semi-supervised classification. They used graph Laplacians with various kinds of covariance functions including squared exponential and mutual $k$-NN type functions. The proposed a Bayesian method how to select the hyperparameters. They applied their methods to various binary classification problems.

None of the above works has shown the relationship between the Laplacian-based method and the traditional kernel and mutual $k$-NN regression. In this paper we have shown that the means of predictive target values in our proposed methods asymptotically converge to the kernel and mutual $k$-NN regression estimate. We have also proposed Bayesian model selection methods for key hyperparameters in both estimates. In [@kim2013transductive] and [@verbeek2006gaussian] symmetry as a covariance matrix related to $k$-NN was not checked, even though [@verbeek2006gaussian] used only symmetric matrixes for the simulation. We have presented the theorems to show the validity of the covariance matrixes used in the proposed methods.

Conclusion
==========

We have proposed Bayesian kernel and mutual $k$-NN regression methods. Those two regression methods work in the framework of Gaussian process regression. Comparing to the traditional kernel regression and mutual $k$-NN regression, it has advantages to provide not only a distribution for the target value but also the principled way to select hyperparameters. Even though the leave-one-out cross-validation can be done for the traditional methods, the performances of the proposed methods were better or comparble in terms of estimation accuracy and computational complexity. Especially, the multiple bandwidth scheme can be applied much more efficiently in Bayesian kernel regression than in traditional kernel regression. The simulation results for the artificial and real world data set show the superiority and efficacy of the proposed methods comparing with the traditional methods.

It is valuable to compare Gaussian process regression (GPR) with squared exponential covariance function (Eq (\[eqn:cov\_ftn\])) to Bayesian kernel regression (BKR), i.e. GPR with Laplacian-based covariance matrix. GPR has an advantage in terms of accuracy, and BKR has an advantage in terms of computational complexity. When we compare the formulations of means of target values in two models, in the one in BKR (Eq (\[eqn:result\])) only a limited correlation among data points is expressed, but in the one in GPR (Eq (\[eq:transd\_lap\_mean\])) much more complicated correlation among data points is expressed. So regression performance in GPR should be better than the one in BKR. However, unllike in GPR, BKR has advantages that it does not require the computation of matrix inversion, which is the computational bottleneck of GPR, for the mean and variance of target value (See Eq (\[eqn:result\]) vs. Eq (\[eq:transd\_lap\_mean\]), (\[eq:transd\_lap\_std\])), and that it requires the computation of matrix determinant but not matrix inversion for the log evidence (See Eq (\[eqn:loglikeli\]) vs. Eq (\[eqn:bkr\_log\_ev\])). (It requires the computation of matrix inversion for the derivative of the log evidence.)

There are some possible future works on the proposed methods. For Bayesian kernel regression, we have had a new bandwidth selection scheme. This scheme can be applied to many application areas including financial engineering and image processing. [@ait1998nonparametric; @ait2000nonparametric] used kernel regression with a heuristic bandwidth selection method in financial engineering. It should be interesting to use Bayesian kernel regression to their tasks, and to compare two results. Data-adapted extensions of kernel regression such as bilateral filtering in image processing can be extended to Bayesian models like in the proposed Bayesian kernel regression so that it may have bandwidth or hyperparameter selection scheme. As future works for Bayesian mutual $k$-NN regression, it may be possible to build Bayesian mutual $k$-NN classifiers, or related models.

[^1]: H.-C. Kim is with R$^2$ Research, Seoul, South Korea. E-mail: hckim.sr[@]{}gmail.com

[^2]: To our knowledge, no Bayesian approach for either $k$-NN regression or mutual $k$-NN regression has been proposed.

[^3]: Some functions in Econometrics package in Octave was used for kernel regression and bandwidth selection by leave-one-out cross-validation

[^4]: Some functions in Econometrics package in Octave was modified and used for the multiple bandwidth scheme in the traditional kernel regression.
