Neural-Kernelized Conditional Density Estimation

Suppose that we are given $T$ pairs of data samples drawn from the joint distribution with density $p(\bm{y},\bm{x})$:

where $\bm{x}(t)$ and $\bm{y}(t)$ denote the $t$-th observations of $\bm{x}$ and $\bm{y}$, respectively, and $d_{\mathrm{x}}$ and $d_{\mathrm{y}}$ denote the dimensions of $\bm{x}$ and $\bm{y}$, respectively. For example, in supervised learning, $\bm{y}$ would be the output data. We will also consider the case of unsupervised learning, where following Hyvärinen and Morioka (2016), $\bm{y}$ could a time index (more on this below). Our primary interest is to estimate the logarithmic conditional density of $\bm{y}$ given $\bm{x}$ up to the partition function from $\mathcal{D}$:

where $q$ is an unnnormalized conditional density, and $Z(\bm{x})$ is the partition function.

In this paper, we propose a function approximation framework where we estimate $\log q(\bm{y}|\bm{x})$ under the following form:

where $\bm{w}(\bm{y})=(w_{1}(\bm{y}),\dots,w_{d}(\bm{y}))^{\top}$ and $\bm{h}(\bm{x})=(h_{1}(\bm{x}),\dots,h_{d}(\bm{x}))^{\top}$ are vector-valued functions.

The model (3) in our framework may seem to be of a restricted form, and unable to express general continuous functions of $\bm{y}$ and $\bm{x}$. However, we next prove that the model can actually approximate a wide range of continuous functions. To this end, we employ the well-known concept of universal approximation capability:

Let us express two sets of functions on domains $\bm{y}\in\mathcal{Y}$ and $\bm{x}\in\mathcal{X}$ by $w_{i}\in\mathcal{F}_{\mathrm{y}}(\mathcal{Y})$ and $h_{i}\in\mathcal{F}_{\mathrm{x}}(\mathcal{X})$ for all $i$, respectively. Eq.(3) indicates that $\log q(\bm{y}|\bm{x})$ is estimated as an element in the following set of functions:

$\mathcal{F}_{\mathrm{y}}(\mathcal{Y})\odot\mathcal{F}_{\mathrm{x}}(\mathcal{X})$ does not necessarily have the universal approximation capability to $\mathcal{C}(\mathcal{Y}\times\mathcal{X})$ in general. Thus, it is important to understand under what conditions $\mathcal{F}_{\mathrm{y}}(\mathcal{Y})\odot\mathcal{F}_{\mathrm{x}}(\mathcal{X})$ has the universal approximation capability. The following proposition states sufficient conditions.

We omit the proof because it is essentially the same as Waegeman et al. (2012, Thereom II.I).

Proposition 2 implies that the model (3) can approximate a wide range of continuous functions on $\mathcal{Y}\times\mathcal{X}$ when both $\mathcal{F}_{\mathrm{y}}(\mathcal{Y})$ and $\mathcal{F}_{\mathrm{x}}(\mathcal{X})$ have the universal approximation capability. Thus, we should use function approximators for $w_{i}$ and $h_{i}$ such that both of them have the universal approximation capability.

It is well-known that the function set of feedforward neural networks have the universal approximation capability (Hornik, 1991). Another example of function sets with the universal approximation capability is reproducing kernel Hilbert spaces (RKHSs) with universal kernels whose example is the Gaussian kernel (Micchelli et al., 2006). Motivated by these facts and Proposition 2, we employ both neural networks and RKHSs in the learning algorithm below.

To estimate $\log q(\bm{y}|\bm{x})$, we employ score matching (Hyvärinen, 2005). It is based on fitting a model $g(\bm{y},\bm{x})$ to $\log q(\bm{y}|\bm{x})$ under the Fisher divergence (Cox, 1985; Sasaki et al., 2014; Sriperumbudur et al., 2017), defined as

where $\nabla_{\bm{y}}:=\left(\frac{\partial}{\partial y_{1}},\dots,\frac{\partial}{\partial y_{d_{\mathrm{y}}}}\right)$. A useful property is that the minimizer of $J(g)$ coincides with $\log p(\bm{y}|\bm{x})$ up to the partition function, which is stated by the following proposition:

The proof can be seen in the supplementary material.

Proposition 3 implies that minimising $J$ asymptotically yields a consistent estimator. However, in practice, we need to estimate $J$ from finite samples. Following score matching (Hyvärinen, 2005) and further developments by Arbel and Gretton (2018), we can derive an easily computable version of $J$ under some mild assumptions as follows:

where $C=\frac{1}{2}\iint\|\nabla_{\bm{y}}\log q(\bm{y}|\bm{x})\|^{2}p(\bm{y},\bm{x})\mathrm{d}\bm{y}\mathrm{d}\bm{x}$, and we used the relation $\nabla_{\bm{y}}\log q(\bm{y}|\bm{x})=\nabla_{\bm{y}}\log p(\bm{y}|\bm{x})$ and applied the integration by parts on the second line under a mild assumption that $\lim_{|y_{j}|\rightarrow\infty}\frac{\partial}{\partial y_{j}}g(\bm{y},\bm{x})p(\bm{y},\bm{x})=0$ for all $\bm{x}$ and $j$. After substituting the model $\bm{w}^{\top}\bm{h}$ into $g$, the empirical version of $\tilde{J}$ can be obtained as

Thus, an estimator can be obtained by minimising $\widehat{J}$ with respect to $w_{i}$ and $h_{i}$ for all $i$.

Motivated by Proposition 2, we model $h_{i}(\bm{x})$ by a feedforward neural network (fNN), $h_{i}(\bm{x};\bm{\theta}_{i})$ where $\bm{\theta}_{i}$ denotes a vector of parameters in the fNN. On the other hand, modelling $w_{i}(\bm{y})$ by fNNs is problematic because $\widehat{J}$ includes the partial derivatives of $w_{i}(\bm{y})$, and computation of the (particularly second-order) derivatives are highly complicated in learning fNNs (Martens et al., 2012). Alternatively, following Proposition 2, we estimate $w_{i}$ as an element in an RKHS with a universal kernel, which has the universal approximation capability (Micchelli et al., 2006). This approach greatly simplifies the derivative computation in $\widehat{J}$ because of the representer theorem: Inspired by the representer theorem for derivatives (Zhou, 2008)¹¹ 1 If we exactly follow the representer theorem in Zhou (2008), $w_{i}$ should have the following form: $w_{i}(\bm{y})=\sum_{t=1}^{T}\sum_{j=1}^{d_{\mathrm{y}}}\alpha_{i,(j-1)T+t}\frac{\partial}{\partial y_{j}^{\prime}}k(\bm{y},\bm{y}^{\prime})\Bigr|_{\bm{y}^{\prime}=\bm{y}(t)}+\beta_{i,(j-1)T+t}\frac{\partial^{2}}{\partial y_{j}^{\prime 2}}k(\bm{y},\bm{y}^{\prime})\Bigr|_{\bm{y}^{\prime}=\bm{y}(t)}$, where $\beta_{i,j}$ denotes coefficients. To decrease the computational costs, we simplify the model by setting all $\beta_{i,j}=0$ and using only a randomly chosen subset of $\{\bm{y}(t)\}_{t=1}^{T}$ as the centers in the kernel function. The same simplification under the Fisher divergence is previously performed and theoretically shown not to incur the performance degeneration under some conditions (Sutherland et al., 2018)., we employ the following linear-in-parameter model for $w_{i}$:

where $k$ and $\alpha_{i,j}$ are the kernel function and coefficients respectively. To decrease the computation cost, we follow the Nystöm approximation (Rudi et al., 2015) and use $\{\tilde{\bm{y}}(b)\}_{b=1}^{B}$ as the center points in $k$ by randomly choosing $B(\ll T)$ points from $\{\bm{y}(t)\}_{t=1}^{T}$.

By substituting an fNN $h_{i}(\bm{x};\bm{\theta}_{i})$ and the linear-in-parameter model $w_{i}(\bm{y};\bm{\alpha}_{i})$ into $\widehat{J}(\{w_{i},h_{i}\}_{i=1}^{d})$, the optimal parameters are obtained as

Setting $\widehat{\bm{w}}(\bm{y}):=(w_{1}(\bm{y};\widehat{\bm{\alpha}}_{1}),\dots,w_{d}(\bm{y};\widehat{\bm{\alpha}}_{d}))^{\top}$ and $\widehat{\bm{h}}(\bm{x}):=(h_{1}(\bm{x};\widehat{\bm{\theta}}_{1}),\dots,h_{d}(\bm{x};\widehat{\bm{\theta}}_{d}))^{\top}$, the estimator is finally given by

We call this method the neural-kernelized conditional density estimator (NKC).

We report the results on the following three methods:

• •

  NKC (Proposed method): $\bm{h}$ was modelled by a feedforward neural network with three layers: The numbers of hidden units in the three layers were $100,50$ and $d$, where $d$ was either $3$ or $5$ (we performed preliminary experiments for $d=1$ as well, but the results were often worse than $d=3,5$). The activation functions were all ReLU. For $\bm{w}$, we employed the Gaussian kernel, and fixed the number of the kernel centers at $B=100$. To perform stochastic gradient descent, RMSprop (Hinton et al., 2012) was used with $128$ minibatches. When learning $\bm{h}$ and $\bm{w}$, model selection was also performed for the learning rate and the width parameter in the Gaussian kernel with a early stopping technique: The datasets were first divided into the training (80%) and validation (20%) datasets, and then $\bm{h}$ and $\bm{w}$ were learned using the training datasets for $100$ epochs. Finally, the best epoch result was chosen as the final result using the validation dataset in terms of $\widehat{J}$. Other details are given in the supplementary material.

• •

  Conditional variational autoencoder (CVAE) (Sohn et al., 2015): This is a conditional density estimator based on variational methods, shown for comparison. In CVAE, three conditional densities related to latent variables, which are called recognition, prior and generation networks, were modelled by parametric Gaussian densities: the means and the logarithms of standard deviations in the Gaussian densities of the recognition and generation networks were expressed by neural networks, while the prior network was modelled by the normal density independent of $\bm{x}$ as in (Kingma et al., 2014; Sohn et al., 2015). The neural networks had two layers with $50$ hidden units. The activation function was ReLU in the middle layer, while no activation function was applied in the final layer. The dimensionality of latent variables was $50$. The same optimisation procedures as in NKC was applied. It should be noted that the number of parameters in CVAE was a bit, but not much, larger than in NKC (the exact numbers of parameters are given in the supplementary material).

• •

  LSCDE (Sugiyama et al., 2010)²² 2 http://www.ms.k.u-tokyo.ac.jp/software.html#LSCDE: This is a least-squares estimation method for conditional densities based on the Gaussian kernel, shown for comparison. In the default setting, LSCDE only uses $\min(100,T)$ centers in the Gaussian kernel (Sugiyama et al., 2010), but to apply LSCDE to a large amount of and high-dimensional data, we set $\min(1000,T)$ centers, which were randomly chosen from the whole data. The bandwidth in the Gaussian kernel, as well as the regularisation parameters, were determined by five-fold cross-validation.

The performance of each method was evaluated by the log-likelihood using test data. To estimate the log-likelihood in NKC, we estimated the partition function by importance sampling, while the log-likelihood in CVAE is estimated as in Sohn et al. (2015, Eq.(6)).

Conditional density estimation is particularly useful compared to simple conditional mean estimation in the case where conditional densities are multimodal. Therefore, we first numerically demonstrate the performance of NKC on multimodal density estimation. To this end, we first generated data samples of $\bm{x}$ from the normal density, and then data samples of $y$ was drawn from a mixture of three Gaussians: $p(y|\bm{x})=\sum_{k=1}^{3}\frac{c_{k}}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{(y-\mu_{k}(\bm{x}))^{2}}{2\sigma^{2}}\right)$ where $d_{\mathrm{y}}=1$, $c_{k}=1/3$ for all $k$ and $\sigma=0.25$. The dimension $d_{\mathrm{x}}$ took the value $50$, $100$, $200$ or $300$. Here, $\mu_{k}(\bm{x})$ were modelled by a three layer neural network with leaky ReLU where the number of hidden units is $d_{\mathrm{x}}$, $d_{\mathrm{x}}$ and $1$, and the parameter values were randomly determined in each run. The total number of samples was fixed at $T=30,000$.

The results are presented in the left panel of Table 1. As $d_{\mathrm{x}}$ increases, the performance of LSCDE gets worse, as would be expected from a method based on Gaussian kernels. On the other hand, NKC and CVAE work well even in $d_{\mathrm{x}}=100$. When $d_{\mathrm{x}}=200,300$, NKC performs better than CVAE. Thus, NKC should be a useful method for multimodal density estimation with high-dimensional data. Regarding the computational time, LSCDE was the fastest because the optimal solution in LSCDE is analytically computed. On the other hand, NKC and CVAE had comparable time if any model selection was not used.

Second, we evaluate the performance of NKC on benchmark datasets. The datasets were the publicly available³³ 3 http://www.dcc.fc.up.pt/~ltorgo/Regression/DataSets.html datasets Ailerons $(d_{\mathrm{x}}=26,T=13,750)$, Friedman $(d_{\mathrm{x}}=10,T=40,768)$, Kinematics$(d_{\mathrm{x}}=8,T=8,192)$, and Pumadyn$(d_{\mathrm{x}}=32,T=8,192)$. Before performing conditional density estimation, we first standardised the datasets by subtracting the means and dividing each variable by its standard deviation. Then, in each dataset, 10% of the data samples was used for test and the remaining data samples were for training.

The results are summarised in the right panel of Table 1, showing that NKC compares favourably with CVAE, and clearly beats LSCDE.