Recursive Multikernel Filters Exploiting Nonlinear Temporal Structure

TL;DR

This paper introduces a recursive multikernel filtering method that adaptively combines multiple kernels to better capture nonlinear temporal structures in data, outperforming classical approaches in various tasks.

Contribution

It proposes a novel recursive multikernel filtering approach that adaptively learns to combine kernels, capturing complex temporal dependencies more effectively than previous methods.

Findings

Outperforms classical methods in artificial and real tasks

Automatically learns to process nonlinear temporal information

Effective in both batch and online settings

Abstract

In kernel methods, temporal information on the data is commonly included by using time-delayed embeddings as inputs. Recently, an alternative formulation was proposed by defining a gamma-filter explicitly in a reproducing kernel Hilbert space, giving rise to a complex model where multiple kernels operate on different temporal combinations of the input signal. In the original formulation, the kernels are then simply combined to obtain a single kernel matrix (for instance by averaging), which provides computational benefits but discards important information on the temporal structure of the signal. Inspired by works on multiple kernel learning, we overcome this drawback by considering the different kernels separately. We propose an efficient strategy to adaptively combine and select these kernels during the training phase. The resulting batch and online algorithms automatically learn to process highly nonlinear temporal information extracted from the input signal, which is implicitly encoded in the kernel values. We evaluate our proposal on several artificial and real tasks, showing that it can outperform classical approaches both in batch and online settings.

Tables

Table 1. TABLE I: Experimental comparison of several kernels in the batch case, including a standard RBF kernel with time-delayed embeddings, the two composite recursive kernels from [ 13 ] , and the proposed recursive multikernel (RMK) strategies.

	Standard kernel	Composite recursive kernel		Proposed recursive MK
Dataset	RBF	Average	Symmetric	SimpleMKL	Stacking	Sparse Stacking
MG30	$-18.99$ dB	$-23.26$ dB	$-22.21$ dB	$-23.42$ dB	$-20.05$ dB	$-19.92$ dB
Narendra	$-14.81$ dB	$-15.72$ dB	$-15.09$ dB	$-15.29$ dB	$-17.01$ dB	$-16.78$ dB
Wiener	$-18.20$ dB	$-17.58$ dB	$-17.23$ dB	$-17.56$ dB	$-18.19$ dB	$-18.19$ dB
EEG	$-5.92$ dB	$-3.69$ dB	$-3.60$ dB	$-6.23$ dB	$-6.94$ dB	$-7.69$ dB
Respiratory	$-18.53$ dB	$-14.92$ dB	$-12.58$ dB	$-16.82$ dB	$-18.84$ dB	$-18.64$ dB
EUR-USD	$-25.10$ dB	$-23.24$ dB	$-21.47$ dB	—	$-25.83$ dB	$-25.68$ dB

Table 2. TABLE II: Results for adaptive filtering with different kernels.

Dataset	KLMS	KLMS (RMK)
MG30	$-12.50$ dB	$-17.99$ dB
Narendra	$-5.65$ dB	$-14.55$ dB
Wiener	$-13.42$ dB	$-13.68$ dB
EEG	$-5.32$ dB	$-7.33$ dB
Respiratory	$-12.56$ dB	$-11.53$ dB
EUR-USD	$-21.35$ dB	$-25.13$ dB

Equations

$$y_n$$

$$x_n^i$$

$$y_n$$

$${\varphi }_n^i$$

$$y_n$$

$${\varphi }_n^i$$

$$w^i=\sum _{m=1}^N{\beta }_m^i{\varphi }_m^i,$$

$${\hat{y}}_n=\sum _{i=1}^P\sum _{m=1}^N{\beta }_m^i⟨{\varphi }_m^i,{\varphi }_n^i⟩=\sum _{i=1}^P\sum _{m=1}^N{\beta }_m^i{\kappa }^i(m,n).$$

$${\kappa }^i(m,n)=\{\begin{array}{ll}\kappa (x_m,x_n),& i=1\\ {\bar{\mu }}^2{\kappa }^i(m-1,n-1)& \\ +{\mu }^2{\kappa }^{i-1}(m-1,n-1)& \\ +{\mu }^2\sum _{j=2}^{m-1}{\bar{\mu }}^{j-1}{\kappa }^{i-1}(m-j,n-1)& \\ +{\mu }^2\sum _{j=2}^{n-1}{\bar{\mu }}^{j-1}{\kappa }^{i-1}(m-1,n-j),& i>1\end{array}$$

$$f^i(\mathbf{x})={\mathbf{y}}^T{\left({\mathbf{K}}^i+c\mathbf{I}\right)}^{-1}{\mathbf{\kappa }}^i(x),$$

$$\underset{{\alpha }^1,\dots ,{\alpha }^P}{\min }\left\{\frac{1}{2}\sum _{n=1}^N{\left(y_n-\sum _{i=1}^P{\alpha }^i{f}^i(x_n)\right)}^2\right\},$$

$$\mathbf{\alpha }={\mathbf{F}}^{-1}\mathbf{y},$$

$$f_n^i=f_{n-1}^i+\eta \left[y_n-f_{n-1}^i(x_n)\right]{\kappa }^i(\cdot ,x_n),$$

$${\mathbf{\alpha }}_n={\mathbf{\alpha }}_{n-1}+\nu \left(y_n-\sum _{i=1}^P{\alpha }^i{f}_n^i(x_n)\right){\mathbf{f}}_n(x_n),$$

$${\mathbf{r}}_n^i(m)={\mu }^2\sum _{j=2}^{n-1}(1-\mu {)}^{j-1}{\kappa }^i(m-1,n-j).$$

$${\mathbf{r}}_n^i(m)=(1-\mu )\left({\mathbf{r}}_{n-1}^i(m)+{\mu }^2{\kappa }^i(m-1,n-j-1)\right).$$

$$y_n=0.3y_{n-1}+0.6y_{n-2}+f(e_n),$$

$$f(e)=0.6\sin \left(\pi e\right)+0.3\sin \left(3\pi e\right)+0.1\sin \left(5\pi e\right)$$

$$x_n=b{x}_{n-1}+\sqrt{1-b^2}{e}_x,$$

$$\text{nMSE}=10{\log }_{10}\left({\text{E}}^2/{\hat{\sigma }}_y^2\right),$$

Full text

I Introduction

In recent years, kernel adaptive filters (KAF) have become a popular approach for online machine learning and time-series prediction, thanks to numerous theoretical and practical advances [1, 2, 3]. Differently from deep recurrent neural networks [4], whose training is generally formulated in batch fashion, KAFs can be trained efficiently with a single pass over the training data. Popular examples of KAFs include kernel least-mean-square [5, 2] and kernel recursive least-squares [6, 1]. When operating on temporal data, most KAFs apply common kernel functions, e.g. Gaussian or polynomial, to time-delay embeddings of the input data. Choosing a specific embedding is not trivial in general, and it might be suboptimal in the case of non-stationary signals. These are known problems in other kernel methods as well, such as support vector machines (SVMs) [7] and kernel ridge regression (KRR).

In order to overcome these limitations, several authors have proposed kernelized extensions of classical recursive models, including the recurrent least-squares SVM [8], the autoregressive and moving average (ARMA) SVM [9], the kernel machine and space projection (KMSP) method [10], the kernel $\gamma$-filter [11], and, more recently, the kernel adaptive ARMA algorithm [12]. All of these works share a common methodology, which is composed of three major steps: (i) define a proper state-space model (SSM) in the input space; (ii) map the input and/or the state of the model to a high-dimensional feature map corresponding to a properly defined reproducing kernel Hilbert space (RKHS); and (iii) solve the resulting model by substituting all dot products with evaluations of the associated kernel function according to the so-called “kernel trick”. Although this is a powerful methodology, it is hard to obtain insights from the model and, more importantly, selecting a proper embedding remains a crucial problem.

In this paper, we focus on the alternative methodology recently proposed in [13]. The basic idea consists in defining the SSM explicitly in a proper RKHS, where samples might not correspond one-to-one with the original inputs. In particular, it is possible to define a proper reproducing theorem such that the model can be expressed as a summation of kernel values, which can be computed recursively at each time instant. For the specific case of the $\gamma$-filter, the resulting model is particularly appealing because each “tap” in the RKHS corresponds to a filtering operation on a different time-scale of the original input [13]. Unlike previously proposed recursive kernels, e.g. [14], this class of kernels was found to work robustly in many problems, including time-series prediction and array processing. Some theoretical aspects of a related class of kernels were investigated independently in [15].

Here, we are interested in extending this methodology by focussing on a shortcoming of the original model. In particular, standard kernel methods consider a single kernel value for each input datum, while in this case, we obtain several kernel values per input. In [13], this problem was side-stepped by either averaging the values, or by computing inner products. This approach lacks in flexibility, however, because it implicitly assumes that all time-scales are equally important. In this paper, we put forth the idea of considering each kernel value as coming from a different kernel function, and to apply proper adaptive strategies to learn the dependency with respect to each of them. In the literature, the idea of adaptively combining kernel functions goes under the name of multiple kernel (MK) learning. MK algorithms originated in the SVM literature [16], and their benefits have also been proven by a number of authors for KAFs, including the MK normalized LMS [17], the mixture KMLS [18], the doubly regularized MKLMS [19], and Cartesian HYPASS [3].

In order to test the feasibility of our idea, we focus on a simple strategy in which multiple kernel estimators are linearly combined following a stacking-like [20] algorithm. This allows us to test the same formulation in both batch and online settings. In the experimental section, we show that by combining our procedure with the recursive kernel of [13], we obtain a performance that is comparable to or better than competing approaches. To this end, we evaluate several benchmarks using both artificial and real-world datasets.

The rest of the paper is organized as follows. Section II describes the RKHS $\gamma$-filter from [13], which is extended through the proposed MK strategy in Section III. We briefly consider computational considerations of the recursive kernel evaluation in Section IV. Experiments are presented in Section V, before giving some conclusive remarks in Section VI.

II Recursive $\gamma$-filtering in RKHS

Let us denote by $(x_{n},y_{n})$ a generic input-output pair observed at time instant $n$. We assume these data to be generated by the following nonlinear model

[TABLE]

where $x_{n}^{i}$ is the input signal at the $i$-th filter tap, $\langle\cdot\rangle$ denotes the inner product, $w^{i}$ are the filter weights, $f(\cdot),g(\cdot)$ are smooth nonlinear functions, and $e_{x},e_{y}$ represent state and output noise respectively. The main idea introduced in [13] is to model a similar process, defined instead in a proper Hilbert space $\mathcal{H}$

[TABLE]

where $w^{i},\phi^{i}_{n}$ are now samples in the (possibly infinite-dimensional) $\mathcal{H}$, and $e_{\phi}$, $e_{y}$ represent the state and output noise. Differently from previous works, it is not required for $\phi_{n}^{i}$ to have $x_{n}^{i}$ as its preimage, in order to provide more flexibility to the model. Proving a representer theorem in general is not trivial, except for specific instantiations of Eq. (2). A particularly interesting case arises by assuming that the filtering operation in Eq. (2) is a $\gamma$-filter [21, 11] given by (omitting noise for simplicity)

[TABLE]

where $\psi(x_{n})$ is some nonlinear transformation in $\mathcal{H}$, $P$ is the filter length controlling the memory depth, and $0<\mu\leq 1$ is a free parameter controlling stability. It is interesting to observe that each “tap” in the Hilbert space is defined by a recursive equation, so as to work over different temporal combinations of the (nonlinearly transformed) original input sequence. In [13], it is proved that, for a given sequence of length $N$, the filter weights $w^{i}$ can be expressed as

[TABLE]

for some coefficients $\beta_{1}^{i},\ldots,\beta^{i}_{N}\in\mathbb{R}$. After substituting (4) in (3), and making use of the kernel definition $\kappa^{i}(m,n)=\langle\phi_{m}^{i},\phi_{n}^{i}\rangle_{\mathcal{H}}$, we obtain

[TABLE]

Again, it is worth underlining that, in general, $\kappa^{i}(m,n)\neq\kappa^{i}(x_{m}^{i},x_{n}^{i})$. In particular, due to the recursive model definition, the kernel itself can be defined recursively, and the closed-form expression is given by111Note that the original formula in [13, Eq. (13)] was missing a summation from $j=2$ to $n-1$. The correct formula is given by Eq. (6).

[TABLE]

where $\bar{\mu}=1-\mu$, and $\kappa(x_{m},x_{n})$ is the classical (scalar) kernel function corresponding to $\langle\psi(x_{m}),\psi(x_{n})\rangle_{\mathcal{H}}$. The calculation of the associated kernel matrix is illustrated in Fig. 1. More generally, we can replace $x_{m}$ and $x_{n}$ with the corresponding time-delayed embeddings $\bm{x}_{m}$ and $\bm{x}_{n}$ to obtain a more expressive model with similar computational complexity.

The recursive kernel $\gamma$-filter generalizes several known models, such as the classical $\gamma$-filter, recursive auto-regressive filters, and several others (see [13, Section III-C]). Differently from standard KRR and KAF algorithms, which require one kernel value for each time instant, the model structure in (5) requires $P$ kernel values, each of which can be seen as operating at a different time-scale. The filtering approach proposed in [13] did not adapt all coefficients corresponding to all kernel values, but employed a simple combination of the kernels instead. One such combination is the composite kernel $\kappa(m,n)=\frac{1}{P}\sum_{i=1}^{P}\kappa^{i}(m,n)$. However, this approach has a drawback, namely, it assumes that all kernels are equally significant, and it may lose important information when combining them. In the next section, we will describe a more principled approach to combine the different kernels.

III Proposed adaptive multi-kernel approach

We now describe an adaptive formulation that automatically weights the different kernels. The proposed algorithm is relatively inexpensive, and it can be implemented easily in both batch and online settings. Note, however, that nothing prevents the use of more advanced multi-kernel or ensemble strategies to further exploit the multi-kernel structure.

Let us consider the batch case first. Given $N$ training samples $(x_{n},y_{n})$, denote by $\mathbf{y}$ the vector of all $N$ outputs, and by $\mathbf{K}^{i}$ the kernel matrix corresponding to the $i$-th tap in (6). A set of $P$ KRR models is trained as

[TABLE]

where $\bm{\kappa}^{i}(x)=\left[\kappa^{i}(x,x_{1}),\ldots,\kappa^{i}(x,x_{N})\right]$ and $c$ is a regularization constant. We combine the basic models as $h(\mathbf{x})=\sum_{i=1}^{P}\alpha^{i}f^{i}(x)$, where the coefficients $\bm{\alpha}=[\alpha^{1},\ldots,\alpha^{P}]^{T}$ are found by minimizing

[TABLE]

which has an immediate closed form solution

[TABLE]

where $\lfloor\mathbf{F}\rfloor_{ni}=f^{i}(x_{n})$. This formulation is a basic form of what is known as “stacking” in the machine learning literature [20, 22]. When $P\ll N$, Eq. (8) does not require regularization. Otherwise, Eq. (8) can be replaced by a more general leave-one-out strategy as in the original stacking problem [20]. More generally, we can include several constraints and/or regularization terms to Eq. (8) in order to force a specific structure on $\bm{\alpha}$, such as the requirement to lie in a $P$-simplex (similar to an adaptive combination of filters [23]), or impose an $\ell_{1}$-norm regularization to remove unnecessary lags.

A similar stacking strategy can be applied in the online case. In particular, given the new datum $(x_{n},y_{n})$, we first update $P$ KAFs in parallel, for instance employing $P$ KLMS filters [5]

[TABLE]

where $\eta$ is the step-size. Then, we update the current estimate $\bm{\alpha}_{n-1}$ of the weighting coefficients following an instantaneous descent on (8)

[TABLE]

where $\nu$ is a step-size parameter and $\lfloor\bm{f}_{n}(\cdot)\rfloor_{i}=f_{n}^{i}(\cdot)$.

Note that KLMS algorithms require calculating arrays of kernel evaluations in each iteration, which is feasible by applying the recursive formula discussed in the next section. KRLS algorithms, on the other hand, require calculating entire kernel matrices, which is less obvious in recursive settings. This, and other concepts, such as sparsity, require more investigation.

IV Efficient computation of the recursive kernel

We now briefly outline an implementation to compute the recursive kernel using fast operations on column vectors.

The slowest operations in Eq. (6) are the summations of the third and fourth term. In order to speed up the calculation of the last term, we define the column vector $\bm{r}^{i}_{n}$ with elements

[TABLE]

This vector can be obtained recursively by observing that

[TABLE]

A similar recursive calculation is not possible for the third term in Eq. (6), though this term can be obtained efficiently from the elements of the convolution $\bm{\mu}_{mem}\ast\bm{k}^{i}_{n}$, where $\bm{k}^{i}_{n}(m)=\kappa^{i}(m,n-1)$ and $\bm{\mu}_{mem}=\mu^{2}[1-\mu,(1-\mu)^{2},...,(1-\mu)^{N}]^{\top}$.

Fig. 2 compares the computation times of the recursive kernel matrix with $P=5$, for an explicit implementation of Eq. (6) and the discussed efficient implementation, both in Matlab R2015a on an Intel Core i7 PC with $3.4$ GHz processor.

V Experimental results

The proposed batch and online algorithms are evaluated on six benchmark problems, three of which are defined on artificially generated datasets and the rest on real-world data.

The first artificial dataset is the Mackey-Glass time-series (with delay $30$) taken from [13], denoted as MG30, on which we perform one-step-ahead prediction. The second task is a noisy version of a nonlinear prediction problem introduced in [24], denoted as “Narendra”, which is defined by

[TABLE]

where the unknown function $f(\cdot)$ has the form

[TABLE]

and $e_{n}=\sin\left((1+a)\omega_{0}n\right)$, $\omega_{0}=2\pi/250$, $a$ is uniformly distributed in the interval $\left[0.1,2.9\right]$, and we set $y_{-1}=y_{0}=1$. We additionally add Gaussian noise with variance $0.1$ to the desired output during training.

The third task on artificial data is the identification of the nonlinear Wiener model described in [25], where the input is generated according to

[TABLE]

with $x_{0}$ randomly generated according to a uniform distribution, $e_{x}$ is Gaussian noise with variance $0.1$, and we set $b=0.8$. The output is given by first applying a linear filter to an embedding of the last $8$ inputs, and then applying a soft nonlinearity on the resulting scalar value.

The first problem on real-world data considers the $4$-step ahead prediction of the EEG dataset from [13], extracted from the MIT-BIH Polysomnographic Database222https://www.physionet.org/physiobank/database/slpdb/. We then consider the prediction of a respiratory motion trace recorded at the Georgetown University Hospital333http://signals.rob.uni-luebeck.de/, originally described in [26]. Finally, the EUR-USD dataset contains the EUR vs. USD exchange rates in minute intervals, taken on the days January 2nd and 5th of 2009. The task is $2$-step ahead prediction.

V-A Batch experiments

In all batch experiments, we use $200$ elements for training and $1000$ separate elements for testing, except for EUR-USD, where we use $1440$ samples for training and $1370$ for testing, respectively. Parameters are fine-tuned following the same grid-search procedure described in [13], optimizing over a third, independent validation set.

We evaluate several KRR models, trained using (i) a standard RBF kernel with time-delayed embeddings on input, (ii) the recursive kernel with the two composition strategies described in [13] (averaging and symmetrization), and (iii) the stacking procedure described in Section III, both with $\ell_{2}$ and $\ell_{1}$ regularization. Additionally, we consider an extension of our idea by training a support vector regression model with the SimpleMKL algorithm [27], which finds a new kernel via an adaptive combination of the base kernel matrices. Denoting by $\text{E}^{2}$ the mean-squared error (MSE) over the test set, we evaluate the models using a normalized MSE on a logarithmic scale,

[TABLE]

where $\hat{\sigma}_{y}^{2}$ is the empirical variance of the output computed over the test set.

The results for the different benchmarks are given in Table I, where the best result for each dataset is highlighted with a bold font. We observe that the proposed stacking procedures are outperforming the standard KRR models in $4$ out of $6$ benchmarks, while SimpleMKL achieves slightly better performance in the Mackey-Glass task. To confirm the results, we employ the corrected Friedman test described in [28], according to which the performance of the algorithms are statistically different with a confidence value of $\alpha=0.05$. A set of Nemenyi post-hoc tests shows that the performance of the base stacking procedure is statistically better than the standard kernel and the two composite recursive kernels, while the sparse stacking procedure is statistically better than the standard kernel and the symmetric recursive kernel, only.

V-B Online experiments

In a second set of experiments we evaluate the online performance of the described KLMS algorithm with recursive multikernel, and classical KLMS [5], on the six benchmark problems. Both algorithms use the same kernel parameter and learning rate. Table II lists the nMSE results obtained after convergence, indicating clear benefits of the RMK strategy.

Finally, in Fig. 3 we reproduce the online learning experiment from [5, Sec. 2.11.2]. In this experiment, a binary signal is fed into a nonlinear communications channel, and the goal is to estimate the correct input signal given its output, i.e. to construct a channel equalizer. KLMS with RMK employs $P=5$ lags and $\mu=0.9$, and the same kernel parameter and learning rate as standard KLMS. While both nonlinear algorithms converge to a similar MSE, the RMK algorithm enjoys a much faster convergence rate.

VI Conclusions

In this paper, we proposed a novel approach for exploiting multiple time-scale information in kernel regression and filtering problems by combining a previously introduced $\gamma$-filter (defined in a proper Hilbert space), with an adaptive strategy for combining kernel functions. In this way, the algorithm automatically adapts to the most significant time-scales of the original input signal. Additionally, the kernel functions are particularly suitable for online processing because they can be recursively computed from previous values, and we briefly discussed on their efficient implementation.

Experimental simulations show that the algorithm has similar or better performance on a wide range of tasks, when compared to several alternative strategies such as time-delayed embeddings of the input. In future works, we plan to further extend our idea by leveraging upon the recent literature on MK filters, and by exploring nonlinear combinations of the different kernels.

Acknowledgment

S. Van Vaerenbergh is supported by the Spanish Ministry of Economy and Competitiveness (under project TEC2014-57402-JIN). S. Scardapane is supported in part by Italian MIUR, “Progetti di Ricerca di Rilevante Interesse Nazionale”, GAUChO project, under Grant 2015YPXH4W_004.