A novel Multiplicative Polynomial Kernel for Volterra series identification

TL;DR

This paper introduces a new kernel-based regularization method for Volterra series identification, improving the selection of influential monomials and enhancing prediction accuracy in nonlinear system modeling.

Contribution

It proposes a novel multiplicative polynomial kernel for Volterra models, with parameter estimation via marginal likelihood, outperforming existing methods in monomial selection and prediction.

Findings

Better monomial selection improves model accuracy

Numerical experiments demonstrate superior prediction capabilities

Kernel approach effectively handles curse of dimensionality

Abstract

Volterra series are especially useful for nonlinear system identification, also thanks to their capability to approximate a broad range of input-output maps. However, their identification from a finite set of data is hard, due to the curse of dimensionality. Recent approaches have shown how regularized kernel-based methods can be useful for this task. In this paper, we propose a new regularization network for Volterra models identification. It relies on a new kernel given by the product of basic building blocks. Each block contains some unknown parameters that can be estimated from data using marginal likelihood optimization. In comparison with other algorithms proposed in the literature, numerical experiments show that our approach allows to better select the monomials that really influence the system output, much increasing the prediction capability of the model.

Tables

Table 1. Table 1: One-step-ahead (Pred.) and simulation (Sim.) performance of the PK and MPK based estimators obtained in the Silverbox test dataset.

	Pred. (Fit%)	Sim. (Fit%)	Sim. (RMSE [mV])
PK	$97.79$	$81.13$	$17.2213$
MPK	$99.70$	$98.67$	$0.8862$

Equations

$$y_k$$

$$H_i({\mathbf{u}}_k)=$$

$$H_i({\mathbf{u}}_k)={\mathbf{\varphi }}_i^T\left({\mathbf{u}}_k\right){\mathbf{w}}_i\text{,}$$

$$y_k=\mathbf{\varphi }{\left({\mathbf{u}}_k\right)}^T\mathbf{w}+e_k\text{,}$$

$$\mathbf{\varphi }\left({\mathbf{u}}_k\right)={\left[h_0{\mathbf{\varphi }}_1^T({\mathbf{u}}_k)\dots {\mathbf{\varphi }}_r^T({\mathbf{u}}_k)\right]}^T\text{,}$$

$$\mathbf{w}\left({\mathbf{u}}_k\right)={\left[w_0{\mathbf{w}}_1^T\dots {\mathbf{w}}_r^T\right]}^T\text{.}$$

$$\underset{f\in H}{\mathrm{argmin}}\sum _{t=1}^T{\left(y_t-z({\mathbf{u}}_t)\right)}^2+{\gamma }^2||f|{|}_H^2\text{,}$$

$$\hat{f}({\mathbf{u}}_k)={\hat{z}}_k=\sum _{t=1}^T{\alpha }_{t}k({\mathbf{u}}_k,{\mathbf{u}}_t)\text{,}$$

$$k({\mathbf{u}}_k,{\mathbf{u}}_j)=\sum _q{\lambda }_q{\varphi }_q({\mathbf{u}}_k){\varphi }_q({\mathbf{u}}_j)\text{,}$$

$$f({\mathbf{u}}_k)=\sum _q{c}_q{\varphi }_q({\mathbf{u}}_k)\text{,}$$

$$||f|{|}_H^2=\sum _q\frac{c_q^2}{{\lambda }_q}\text{.}$$

$$k^{(r)}({\mathbf{u}}_k,{\mathbf{u}}_j)={\left(1+{\mathbf{u}}_k^T{\mathbf{u}}_j\right)}^r\text{,}$$

$$f({\mathbf{u}}_k)=u_k^3+u_k^2{u}_{k-1}+0.5\text{.}$$

$$k^{(3)}({\mathbf{u}}_i,{\mathbf{u}}_j)$$

$$+3u_i^2{u}_{i-1}{u}_j^2{u}_{j-1}+3u_i{u}_{i-1}^2{u}_j{u}_{j-1}^2$$

$$+3u_i^2{u}_j^2+3u_{i-1}^2{u}_{j-1}^2+6u_i{u}_{i-1}{u}_j{u}_{j-1}$$

$$+3u_i{u}_j+3u_{i-1}{u}_{j-1}+1\text{.}$$

$${\lambda }_{3,0}={\lambda }_{0,3}=1\text{, }{\lambda }_{2,1}={\lambda }_{1,2}=3\text{,}$$

$${\lambda }_{2,0}={\lambda }_{0,2}=3\text{, }{\lambda }_{1,1}=6\text{,}$$

$${\lambda }_{1,0}={\lambda }_{0,1}=3\text{, }{\lambda }_{0,0}=1\text{.}$$

$$k^{(r)}({\mathbf{u}}_k,{\mathbf{u}}_j)=\prod _{i=1}^r\left({\sigma }_{0_i}+({\mathbf{u}}_k{)}^T{\Sigma }_i{\mathbf{u}}_j\right)\text{,}$$

$${\lambda }_{3,0}=\prod _{j=1}^3{\sigma }_0^{(j)}\text{,}{\lambda }_{0,3}=\prod _{j=1}^3{\sigma }_1^{(j)}\text{,}$$

$${\lambda }_{2,1}=\sum _{j=1}^3{\sigma }_1^{(j)}\prod _{l\ne j}{\sigma }_0^{(l)}\text{,}{\lambda }_{1,2}=\sum _{j=1}^3{\sigma }_0^{(j)}\prod _{l\ne j}{\sigma }_1^{(l)}\text{,}$$

$${\lambda }_{2,0}=\sum _{j=1}^3{\sigma }_{0_j}\prod _{l\ne j}{\sigma }_0^{(l)}\text{,}{\lambda }_{0,2}=\prod _{j=1}^3{\sigma }_{0_j}\prod _{l\ne m}{\sigma }_1^{(l)}\text{,}$$

$${\lambda }_{1,1}=\sum _{j=1}^3{\sigma }_{0_j}\sum _{l_1\ne l_2\ne j}{\sigma }_0^{(l_1)}{\sigma }_1^{(l_2)}\text{,}$$

$${\lambda }_{1,0}=\sum _{j=1}^3{\sigma }_0^{(j)}\prod _{l\ne j}{\sigma }_{0_l}\text{,}$$

$${\lambda }_{0,0}=\prod _{j=1}^3{\sigma }_{0_j}\text{,}$$

$${\sigma }_{0_1}=1\text{ , }\left[{\sigma }_0^{(1)}{\sigma }_1^{(1)}\right]=\left[11\right]\text{,}$$

$${\sigma }_{0_2}=1\text{ , }\left[{\sigma }_0^{(2)}{\sigma }_1^{(2)}\right]=\left[10\right]\text{,}$$

$${\sigma }_{0_3}=1\text{ , }\left[{\sigma }_0^{(3)}{\sigma }_1^{(3)}\right]=\left[10\right]\text{,}$$

$${\lambda }_{3,0}=1\text{, }{\lambda }_{2,1}=1\text{, }{\lambda }_{0,3}={\lambda }_{1,2}=0\text{,}$$

$${\lambda }_{2,0}=2\text{, }{\lambda }_{0,2}=0\text{, }{\lambda }_{1,1}=2\text{,}$$

$${\lambda }_{1,0}=2\text{, }{\lambda }_{0,1}=1\text{, }{\lambda }_{0,0}=1\text{.}$$

$${\sigma }_{0_1}=1\text{ , }\left[{\sigma }_0^{(1)}{\sigma }_1^{(1)}\right]=\left[10\right]\text{,}$$

$${\sigma }_{0_2}=1\text{ , }\left[{\sigma }_0^{(2)}{\sigma }_1^{(2)}\right]=\left[11\right]\text{,}$$

$${\sigma }_{0_3}=1\text{ , }\left[{\sigma }_0^{(3)}{\sigma }_1^{(3)}\right]=\left[10\right]\text{.}$$

$${\Sigma }_r$$

$${\Sigma }_i$$

Full text

1 Introduction

In many real world applications, linear models are not able to adequately describe dynamic systems. This can be due to the presence of saturations, quantizers or static nonlinearities at the input and/or the output Ljung (1999)[Section 5]. Even if some insight on the nonlinearities can be available, the formulation of parametric models from finite data records is a difficult task Haber and Unbehauen (1990); Lind and Ljung (2008); Sjöberg et al. (1995). In particular, nonlinear system identification is often seen as an extended parametric regression where the choice of regressors and basis functions plays a crucial role. In this context, Volterra series are especially useful since they can represent a broad range of nonlinear systems Rugh (1980); Boyd and Chua (1985); Cheng et al. (2017). When working in discrete-time, such models correspond to Taylor expansions of the input-output map. Indeed, a truncated Volterra series describes the system as the sum of all the possible monomials in the past inputs and outputs, up to a certain order. The problem is however the curse of dimensionality: the number of monomials grows quickly w.r.t. the polynomial degree and the system memory (given e.g. by the number of past input values that determine the output). Thus, a careful selection of the relevant components to be included in the model is crucial to control the complexity of the estimator, a problem known as regressors selection.

Suboptimal solutions are often searched through greedy approaches like forward/backward subset selection, see for instance Chen et al. (1989); Billings et al. (1989). These regressor selection methods have however difficulties in handling high-dimensional regression spaces. An interesting option is joint estimation and variable selection. This can be performed using e.g. the $\ell_{1}$-norm regularizer which leads to the famous LASSO Tibshirani (1996).

More recent approaches proposed to deal with the dimensionality of Volterra models can be found in Birpoutsoukis et al. (2017); Stoddard et al. (2017). In particular, in Birpoutsoukis et al. (2017), inspired by ideas developed for linear system identification in Pillonetto and De Nicolao (2010), the authors proposed a regularization strategy suitable for Volterra series with smooth exponential decay.

An alternative route to the approaches mentioned above is the use of kernel-methods, which lead e.g. to the so called regularization networks Poggio and Girosi (1990). Here, an unknown function is determined as the minimizer of an objective that is sum of two terms: a quadratic loss and regularizer defined by a positive definite kernel. The choice of the kernel has a major effect on the quality of the estimate since it encodes the expected properties of the function to reconstruct. Just looking at the function to reconstruct as the unknown system (input-output map), in recent years kernel-based approaches have been widely exploited also for nonlinear system identification and prediction, see for instance Espinoza et al. (2005); Hall et al. (2012).

Another popular model is the polynomial kernel, which has a deep connection with Volterra series. In fact, it implicitly encodes all the monomials up to the desired degree $r$, a kernel parameter tunable by the user. Regularization networks for efficient Volterra identification that exploit this kernel can be found in Franz and Schölkopf (2006).

The approach described in this paper is based on a new polynomial kernel, named Multiplicative Polynomial kernel (MPK). Similarly to the polynomial kernel, the MPK encodes all the monomials up to degree $r$, but it is defined by the product of $r$ linear kernels, each one equipped with a distinct set of hyperparameters. The MPK has some important features w.r.t. the polynomial kernel used in Franz and Schölkopf (2006). As already said, the polynomial kernel depends only on the polynomial degree $r$ and it encodes a number of monomials rapidly increasing with $r$ and the system memory. When plugged in a regularization network, it induces a penalty that cannot promote any sparsity in the solution. On the contrary, the MPK is equipped with a augmented set of hyperparameters, which allows promoting the monomials w.r.t. their maximum relative degree, improving regularization performance. Tests performed both in simulation and with data coming from a real system show that the MPK hyperparameters can be tuned via marginal likelihood optimization or cross-validation.

The paper is organized as follows. In Section 2 we provide a brief overview on Volterra series and the main identification approaches adopted. In Section 3 we highlight some critic aspects of the standard polynomial kernel, and we introduce our kernel function, the Multiplicative Polynomial Kernel, highlighting its regularization capabilities. Finally, in Section 4 we report numerical results, in which we compare performance of the proposed kernel and the standard polynomial kernel.

2 BACKGROUND

2.1 Volterra series

Let $u_{k}$ and $z_{k}$ be the one dimensional input and output signals at time $k$. When modeling the system response with a discrete time Volterra series of order $r$ the noisy output $y_{k}$ is assumed to be the sum of measurement noise and $r+1$ contributions acting on the lagged inputs $u_{k},\,u_{k-1},\,u_{k-2}\,\dots$. Assume that the system has finite memory $m$, and define the input vector $\boldsymbol{u}_{k}=[u_{k}\,\dots\,u_{k-m}]$. Then we have

[TABLE]

where $e_{k}\sim N\left(0,\sigma_{n}^{2}\right)$ is the measurement noise, $h_{0}$ is a constant accounting for the zero-order Volterra contribution, while the $H_{i}$ are the i-th Volterra contributions. In particular, each $H_{i}$ is defined as the convolution between $\boldsymbol{u}_{k}$ and a Volterra map $h_{i}$, namely,

[TABLE]

In this paper, we consider symmetric Volterra kernels, i.e., given a set of lags $\tau_{1},\dots,\tau_{i}$, the value of $h_{i}$ is independent on the lags order. For instance, with $i=2$ we have $h_{2}\left(\tau_{1},\tau_{2}\right)=h_{2}\left(\tau_{2},\tau_{1}\right)$.

Alternatively, each $H_{i}$ term in (1) can be rewritten more compactly as an inner product. Let $\boldsymbol{\phi}_{i}(\boldsymbol{u}_{k})$ be a vector collecting all the distinct monomials with degree $i$ in the components of $\boldsymbol{u}_{k}$. We have

[TABLE]

where $\boldsymbol{w}_{i}$ is the vector collecting the distinct $h_{i}$ coefficients ordered in accordance with $\boldsymbol{\phi}_{i}$, and opportunely scaled to account for repetitions due to symmetry. More precisely, the $\boldsymbol{w}_{i}$ entry associated to the monomial $\prod_{j=0}^{m}u_{k-j}^{d_{j}}$ is the product between the correspondent $h_{i}$ coefficient scaled by the multinomial coefficient $\binom{i}{d_{0},\dots,d_{m}}$. Based on the above definitions, (1) can be rewritten as

[TABLE]

where

[TABLE]

The authors in Birpoutsoukis et al. (2017) have proposed to learn the input-output relations by directly estimating the elements of $\boldsymbol{w}$. This estimation is performed by solving a least square problem defined by (3), given a data set of input-output measurements $D=\left\{\left(\boldsymbol{u}_{k},y_{k}\right),\,k=1,\dots T\right\}$.

It is worth stressing that the applicability of algorithms based on a least-square approach is strongly limited by the high computational and memory requirements related to the dimension of $\boldsymbol{w}$. Indeed, the number of Volterra coefficients grows rapidly with the system memory $m$ and the Volterra order $r$. More precisely, assuming that the Volterra maps are symmetric, we have that $\boldsymbol{w}_{i}$ is composed of $N_{i}=\binom{m+i}{i}$ elements, leading to a total number of $N=1+\sum_{i=1}^{r}N_{i}$ parameters to be estimated.

2.2 Polynomial kernel and Volterra series

An alternative solution to accomplish the Volterra series identification has been proposed in Franz and Schölkopf (2006). Instead of formulating the identification problem directly w.r.t. $\boldsymbol{w}$, the authors rely on kernel based techniques. The input-output map $f$ is assumed to belong to a reproducing kernel Hilbert space (RKHS) Scholkopf and Smola (2001), defined by a kernel function $k(\boldsymbol{u}_{k},\boldsymbol{u}_{j})$. Given an input-output dataset $D$ like that previously introduced, $\hat{f}$, the estimate of $f$, is obtained solving the following problem,

[TABLE]

where the first term of the loss function accounts for the adherence to experimental data, while the second is the regularization term, given by the squared RKHS norm of $f$; the balance between these two contributions can be controlled tuning the hyperparameter $\gamma$. According to the representer theorem, $\hat{f}$ is expressed in closed form as

[TABLE]

where $\boldsymbol{\alpha}=\left[\alpha_{1},\dots,\alpha_{T}\right]^{T}$ is equal to $(K+\gamma^{2}I_{T})^{-1}\boldsymbol{y}$, $\boldsymbol{y}=\left[y_{1},\dots,y_{T}\right]^{T}$ denotes the vector containing all the output measurements, and $K$ is the Kernel matrix, i.e. its $(k,j)$ entry is $K_{k,j}=k(\boldsymbol{u}_{k},\boldsymbol{u}_{j})$.

For our future use, it is also useful recalling the following fundamental facts regarding RKHS theory. Under mild assumptions, a kernel function admits an expansion (possibly infinite) in terms of basis functions $\phi_{q}$, namely,

[TABLE]

where $\lambda_{q}$ are positive scalars. It can then be proved that any function in the RKHS induced by the above kernel has the representation

[TABLE]

for suitable coefficients $c_{q}$. In addition, if all the basis functions $\phi_{q}$ are linearly independent, one also has

[TABLE]

This last relation shows how the $\lambda_{q}$ coefficients are related to each $\phi_{q}$ in determining the regularization term present in (6). In particular, small values of $\lambda_{q}$ entail high penalization of $\phi_{q}$.

As far as the kernel function is concerned, in Franz and Schölkopf (2006) the authors considered the polynomial kernel. In particular, we discuss the inhomogeneous polynomial kernel, defined as

[TABLE]

where $k$ is a tunable hyperparameter corresponding to the degree of the polynomial kernel. As showed in Scholkopf and Smola (2001), the polynomial kernel in (11) admits an expansion in the monomials in $\boldsymbol{u}_{k}$, with degree up to $k$. Namely, referring to (8) and (10), we have that the $\phi_{q}$ corresponds to the elements of the $\boldsymbol{\phi}$ vector defined in (4); accordingly, $1/\lambda_{q}$ then defines the penalty assigned to the relative monomial.

We conclude this subsection with a computational note. The computation of $\hat{f}$ involves the inversion of a $T\times T$ matrix, see (7). The number of operations so scales with the cube of the data set size, and there is no direct dependence on $N$, the dimensions of $\boldsymbol{\phi}$, allowing the use of high-order Volterra models.

3 PROPOSED KERNEL

The Volterra series learning strategy we propose is based on a novel polynomial kernel, called Multiplicative Polynomial Kernel (MPK). Compared to the standard polynomial kernel reported in (11), our kernel is equipped with a set of parameters which allows to assign suitable priors to the different basis functions of the RKHS, thus leading to better performance in terms of estimation and generalization. Before describing the proposed kernel function we highlight some critical issues of standard inhomogeneous polynomial kernel.

3.1 Penalties induced by (11)

As stated in Rasmussen and Williams (2006) (Chapter 4.2.2), polynomial kernels are not widely used in regression problems, since they are prone to overfitting, in particular, in presence of high dimensional inputs and when the degree is greater than two. Indeed, in the kernel formulation given in (11) there are not parameters that allow to weigh differently the monomials composing the RKHS.

To clarify this concept we consider a simple example, a third order Volterra series with $m=1$ defined as follow

[TABLE]

We compute the $\lambda_{q}$ obtained with the standard polynomial kernel expanding (11), and comparing the result with (8). For the sake of clarity, we will denote with $\lambda_{d_{0},\dots,d_{m}}$ the penalty coefficient associated to the monomial $\prod_{\tau=0}^{m}u^{d_{\tau}}_{k-\tau}$. The kernel expansion is

[TABLE]

Then, by inspection, we obtain

[TABLE]

These values show that (11) assigns penalties based on the monomial degree, and promoting mixed terms. This trend might not be representative of the Volterra kernel, leading to the need of more training data to obtain accurate estimates. For instance, consider the function reported in (12). It is evident that the $\lambda$ values obtained with (11) do not describe properly the contributions of the different monomials, since, for example, the higher values of $\lambda$ are assigned to monomials that are not present.

3.2 Multiplicative Polynomial Kernel

The kernel function we propose to model the $r$ order Volterra series is given by the product of $r$ linear kernels, and it is formally defined as

[TABLE]

where the matrices $\Sigma_{i}\in\mathbb{R}^{(m+1)\times(m+1)}$ are diagonal. In particular, for each $i$ we have $\Sigma_{i}=diag\left(\left[\sigma^{(i)}_{0},\dots,\sigma^{(i)}_{m}\right]\right)$, with the diagonal elements greater of equal than zero.

Exploiting the kernel properties it can be easily shown that the function defined in (13) is a well-defined kernel function, since it is the product of several valid kernel functions, see Rasmussen and Williams (2006).

3.3 Penalties induced by the MPK

In this subsection we analyze the advantages of the proposed kernel function, focusing in the role played by the kernel parameters. To this aim, we consider the example analyzed in the previous subsection, that is, the identification of the input-output behavior of a Volterra series with $r=3$ and $m=1$. Starting from the kernel definition given in (13), through standard algebraic computations, we can derive the penalties coefficients as functions of the MPK parameters. In particular, the penalties assigned to monomials of degree three are

[TABLE]

the ones assigned to monomials of degree two are

[TABLE]

and, finally, the ones assigned to monomials of degree one and zero are

[TABLE]

Some interesting insights can be obtained from the previous penalties expressions. Notice that the MPK parameters allow penalizing the monomials w.r.t. their relative degree. For instance, consider the penalties assigned to monomials that contains $u_{k}$, namely, $\lambda_{3,0}$, $\lambda_{2,1}$, $\lambda_{1,2}$, $\lambda_{2,0}$, $\lambda_{1,1}$, and $\lambda_{1,0}$, as function of $\sigma^{(j)}_{0}$, with $j=1,2,3$. Analyzing these expression we can appreciate that to promote monomials in which $u_{k}$ appears with degree i at least $i$ of the $\lambda^{(j)}_{0}$ parameters need to be significantly greater than zero. On the contrary, to penalize monomials in which $u_{k}$ appears with relative degree greater than $i$ we need to set to zero $r-i$ of the $\sigma^{(j)}_{0}$ elements. Thus, referring to the test function in (12), where the maximum relative degree of $u_{k}$ and $u_{k-1}$ are, respectively, $3$ and $1$, we can exclude part of the monomials that do not influence the input-output defining

[TABLE]

which determine the following penalties,

[TABLE]

The hyperparameters tuning can be accomplished relying on empirical methods, like cross validation, or optimizing a given loss function, like Marginal Likelihood (ML).

3.4 Parametrization of the $\Sigma_{i}$ matrices

Notice that the MPK is given by the product of $r$ equal blocks. Consequently, permuting the $\Sigma_{i}$ matrices we obtain different configurations of the MPK hyperparameters associated to the same set of penalties. For instance, referring to the test case analyzed in the previous subsection, the configuration reported in (14) and the following configuration are associated to the same set of penalties, i.e. (15),

[TABLE]

When optimizing the hyperparameters by ML, this fact could lead to undesired behaviors, due to the presence of several local maxima. To avoid such behaviors, we propose a iterative parametrization of the $\Sigma_{i}$ diagonal elements. More specifically, the $\Sigma_{i}$, are defined by a backward iteration as follows,

[TABLE]

where the $a^{(i)}_{j}$ elements are greater of equal than zero. We conclude by emphasizing the relation between the proposed re-parameterization and the relative degree with which each term appears. Notice that increasing $a^{(i)}_{j}$ we promote simultaneously all the monomials in which $u_{k-j}$ appears with relative degree up to $i$. Moreover, to penalize monomials in which $u_{k-j}$ appears with relative degree greater than $i$ we need to set to zero the $a^{(l)}_{j}$ with $l>i$.

4 EXPERIMENTAL RESULTS

4.1 Simulated environment

In this set of experiments we test the performance of the MPK in a simulated environment, the benchmark system introduced in Spinelli et al. (2005), i.e. a third order Volterra series described by the following equation

[TABLE]

The MPK-based estimator is compared with the one based on the polynomial kernel reported in (11), hereafter denoted with PK; both for Pk and MPK we considered $r=3$. Input signals are $1000$ samples obtained from a realization of Gaussian noise. Concerning $m^{tr}_{u}$, $m^{ts}_{u}$, $\sigma^{tr}_{u}$ and $\sigma^{ts}_{u}$, respectively, the input mean and standard deviation of the training and test samples, we consider four different scenarios:

• •

Experiment 1: $m^{tr}_{u}=m^{ts}_{u}=0$, $\sigma^{tr}_{u}=\sigma^{ts}_{u}=4$;

• •

Experiment 2: $m^{tr}_{u}=m^{ts}_{u}=0$, $\sigma^{tr}_{u}=\sigma^{ts}_{u}=2$;

• •

Experiment 3: $m^{tr}_{u}=-12$, $m^{ts}_{u}=12$, $\sigma^{tr}_{u}=\sigma^{ts}_{u}=4$;

• •

Experiment 4: $m^{tr}_{u}=-12$, $m^{ts}_{u}=12$, $\sigma^{tr}_{u}=\sigma^{ts}_{u}=2$.

In all the experiments the noise standard deviation is $\sigma_{n}=4$. The PK and MPK hyperparameters have been trained optimizing the ML of the training samples. As concerns the optimization, we used standard gradient descent algorithm, with adaptive learning rate. The two estimators are implemented in PyTorch Paszke et al. (2017), to exploit automatic differentiation functionalities.

The four experiments can be grouped in two sets. Generalization properties are stressed more in Experiment 3 and Experiment 4 since $m^{tr}_{u}\neq m^{ts}_{u}$, and hence the training and test input signals are significantly different with each other; in particular, the mean values are such that with high probability the test inputs are outside the $3\sigma$ of training inputs distribution. In Experiment 1 (resp. Experiment 3) and Experiment 2 (resp. Experiment 4) we considered different values of the input standard deviation in order to analyze the estimators behaviors with different system excitations.

For each experiment we performed a Monte Carlo of 100 simulations. In each simulation the same training and test data sets have been used to implement and test the MPK and PK based estimators. Results are reported in Figure 1. Performance is measured by the percentage fit (Fit%), defined as

[TABLE]

where $\boldsymbol{z}$ and $\hat{\boldsymbol{z}}$ are the vectors collecting the real and estimated system output, while $\bar{z}$ is the mean of $\boldsymbol{z}$.

In Figure 1, we visualized the results obtained through some boxplots. The estimator based on the MPK outperforms the standard polynomial kernel, since in all the tests the Fit% obtained with MPK are higher than the ones obtained with PK. Besides improving the estimation accuracy, the MPK parametrization improves also the generalization performance. Indeed, comparing results obtained in Experiment 1 and Experiment 3, we can appreciate how the penalties learned by the MPK estimator by optimizing ML provides robustness to variations of the inputs distribution; more specifically, the MPK performance decreases less than the PK performance when the test input locations that are far from the training inputs.

As far as variations of the system excitation, comparing results obtained in Experiment 3 and Experiment 4, we can observe how not sufficiently exciting training samples can lead to a bad identification of the MPK parameters. Notice how from Experiment 3 to Experiment 4, the variance of MPK based estimator grows up more that the one of the PK based, highlighting the importance of using sufficiently exciting input trajectories.

4.2 Identification of the Silverbox system

In this subsection we test the estimators based on MPK and PK with data collected on a real system, the Silverbox Wigren and Schoukens (2013). The Silverbox is an electrical system that simulates a mass-spring-damper system. The spring exhibits non-linear behaviors, and the system is described by the following differential equation

[TABLE]

where $z$ and $u$ are, respectively, the mass displacement and the input force applied to the mass, while $d$, $k_{1}$, and $k_{3}$ are the parameters of the damper and the non-linear spring.

We used the MPK and PK based estimators to learn the evolution of the mass displacement, modeling the unknown target function with a third order Volterra series. We considered as input of the model the past $m=5$ $u$ and $y$, namely, at time $k$, the model input and output are, respectively, $\left[u_{k}\,\dots u_{k-m}\,z_{k-1}\,\dots\,z_{k-m}\right]$ and $z_{k}$. The original training and test dataset are publicly available111http://www.nonlinearbenchmark.org/. The training dataset accounts approximately for $80000$, obtained inputting to the system an odd random phase multisine signal, while the test set accounts for approximately $40000$ samples, collected exciting the system with filtered Gaussian noise.

To further stress generalization properties, we derived and trained the estimators using just the first $200$ samples of the training dataset. Besides testing the one-step-ahead prediction performance, we measured also the simulation performance: given an initial state of the system we simulate its evolution evaluating iteratively the function learned for the one-step-ahead problem, inputting to the estimator the past predicted output instead of the measured output. In this context, we noticed that optimizing the hyperparameters through cross-validation increases the simulation performance. To deal with the considerable number of parameters, we adopted a gradient-based strategy. Specifically, we randomly select $5$ partitions of the training data. Each partition accounts for two sets, composed of $100$ samples. As loss loss function we considered the sum of the mean squared errors (MSEs) in validation. Namely, for each partition, we derive the estimator based on the first set, and we compute the MSE of the second set. The loss function is the sum of the MSEs. The kernel hyperparameters are updated minimizing the loss though a gradient-based procedure, iterated until convergence of the loss.

Performance is reported in Table 1. As before, we compare the one-step-ahead performance using the Fit%. Simulation performance are measured both with Fit% and root mean squared error (RMSE). The MPK outperforms the standard polynomial kernel, both in one-step-ahead prediction and simulation. The gap is particularly evident in simulation, where MPK significantly outperforms PK; see Figure 2 to compare the evolution of the simulation errors. Despite we used just the first $200$ training samples to derive the model, and hyperparameters optimization was focused on optimizing the one-step-ahead performance, the MPK performance is close ot the best results obtained in this benchmark.

5 CONCLUSIONS

In this paper we have introduced a novel polynomial kernel, name Multiplicative Polynomial Kernel. Compared to the standard polynomial kernel, the MPK is equipped with a set of parameters that allows to better select the monomials that really influence the system output. As proven by numerical results, this fact entails a performance improvement, both in terms of accuracy and generalization properties.