Stable spline identification of linear systems under missing data

TL;DR

This paper introduces a nonparametric Bayesian approach using stable spline kernels for identifying linear systems with missing data, improving model prediction accuracy over traditional methods.

Contribution

It presents a novel class of radial basis function kernels derived from stable spline kernels, enabling effective system identification with incomplete data.

Findings

The new kernel approach outperforms standard parametric methods in prediction accuracy.

The method is applicable to both discrete and continuous-time systems.

Numerical experiments validate the effectiveness of the proposed approach.

Abstract

A different route to identification of time-invariant linear systems has been recently proposed which does not require committing to a specific parametric model structure. Impulse responses are described in a nonparametric Bayesian framework as zero-mean Gaussian processes. Their covariances are given by the so-called stable spline kernels encoding information on regularity and BIBO stability. In this paper, we demonstrate that these kernels also lead to a new family of radial basis functions kernels suitable to model system components subject to disturbances given by filtered white noise. This novel class, in cooperation with the stable spline kernels, paves the way to a new approach to solve missing data problems in both discrete and continuous-time settings. Numerical experiments show that the new technique may return models more predictive than those obtained by standard parametric Prediction Error Methods, also when these latter exploit the full data set.

Equations

$${\mathbf{y}}_o=[y(t_1),\dots ,y(t_n){]}^{\top }.$$

$$y(t)=\left\{\sum _{\ell =1}^{p-1}(u_{\ell }\otimes f_{\ell })(t)\right\}+(e\otimes f_p)(t)$$

$$D_t=\left\{y(x),u_{\ell }(x):x\le t,\ell =1,\dots ,p-1\right\},x\in \mathbb{Z}\vee \mathbb{R}$$

$$\tilde{y}(t)=\mathcal{F}(D_t)$$

$$\mathcal{F}(D_t)=\left\{\sum _{\ell =1}^{p-1}(u_{\ell }\otimes g_{\ell })(t)\right\}+(y\otimes g_p)(t)$$

$$\left[\begin{array}[]{c}{\bf y}_{o}\\ {\bf y}_{m}\end{array}\right]=\left[\begin{array}[]{c}\Phi_{o}\\ \Phi_{m}\end{array}\right]\,{\bf f}+\left[\begin{array}[]{c}{\boldsymbol{\xi}}\\ {\boldsymbol{\xi}}_{m}\end{array}\right]$$

$${\mathbf{R}}_m:=\mathbb{E}[{\mathbf{\xi }}_m{\mathbf{\xi }}^{\top }],\mathbf{R}:=\mathbb{E}[\mathbf{\xi }{\mathbf{\xi }}^{\top }]$$

$$\begin{array}[]{rcl}\hat{\bf y}_{m}&:=&\hat{\mathbb{E}}\left[{\bf y}_{m}|{\bf y}_{o}\right]=\Phi_{m}\hat{\mathbb{E}}\left[{\bf f}|{\bf y}_{o}\right]+\hat{\mathbb{E}}\left[{\boldsymbol{\xi}}_{m}|{\bf y}_{o}\right]\\ &=&\left(\Phi_{m}\Sigma_{f}\Phi_{o}^{\top}+\mathbf{R}_{m}\right)\left(\Phi_{o}\Sigma_{f}\Phi_{o}^{\top}+\mathbf{R}\right)^{-1}{\bf y}_{o}\end{array}$$

$$\mathbb{V}[f_{\ell }(t)f_{\ell }(s)]=\lambda K_{\ell }(t,s),\ell =1,2,\dots ,p$$

$$W_q(s,t)={\int }_0^1{G}_q(s,u)G_q(t,u)du$$

$$\\ G_{q}(r,u)=\frac{(r-u)_{+}^{m-1}}{(m-1)!},\qquad(u)_{+}=\left\{\begin{array}[]{cl}u&\mbox{if}~{}u\geq 0\\ 0&\mbox{otherwise}\end{array}\right.$$

$$K(s,t):=W_q(e^{-\beta s},e^{-\beta t}),m=1,2,\dots .$$

$$K(s,t)=\frac{e^{-\beta (s+t)}{e}^{-\beta \max (s,t)}}{2}-\frac{e^{-3\beta \max (s,t)}}{6}$$

$$\xi (t):=(e\otimes f_p)(t)$$

$$R(s,t)$$

$$h(x)$$

$$h(x)=\frac{e^{-\beta |x|}}{\beta }$$

$$h(x)=\frac{3e^{-2\beta |x|}-e^{-3\beta |x|}}{18\beta }$$

$$K(s,t):=W_q(e^{-\beta s},e^{-\beta t}),(s,t)\in \mathbb{N}\times \mathbb{N}$$

$$R(s,t)$$

$$h(x)$$

$$h(x)=\frac{e^{-\beta (|x|+1)}}{1-e^{-\beta }}$$

$$h(x)=\frac{3e^{-2\beta |x|}-e^{-3\beta |x|}}{6}\frac{e^{-3\beta }}{1-e^{-3\beta }}$$

$$f_{\ell }=g_{\ell }\otimes h_{\ell }$$

$$H(z)=\frac{z^2}{z^2+\varphi z+\phi }.$$

$$[\mathbf{R}{]}_{ij}=R(t_i,t_j)$$

$$P(x,y)=\sum _{\ell =1}^{p-1}\left(u_{\ell }\otimes h_x\right)(y)$$

$$h_x(t)=(u_{\ell }(\cdot )\otimes K(t,\cdot ))(x)$$

$$[\mathbf{P}{]}_{ij}=P(t_i,t_j)$$

$$[{\mathbf{K}}_{\ell }^{\infty }{]}_{ij}=K_{\ell }(i,j),(i,j)\in \mathbb{N}\times \mathbb{N}$$

$$P(i,j)=\sum _{\ell =1}^{p-1}{\mathbf{U}}_{\ell }^i{\mathbf{K}}_{\ell }^{\infty }{\left({\mathbf{U}}_{\ell }^j\right)}^{\top }$$

Full text

1 Introduction

The main approach to identification of time-invariant linear discrete-time models is given by parametric Prediction Error Methods (PEM), see [1]. As a rule, model complexity is unknown and model-order selection is a key ingredient of the identification process. Models of different order are identified from data and compared resorting to either complexity measures such as AIC or cross validation, see e.g. [2].

Recent work has proposed an alternative nonparametric approach that focuses on the direct identification of the impulse response [3, 4, 37]. Following the Gaussian regression framework [5], the unknown impulse response is modeled as a Gaussian process whose autocovariance encodes the available prior knowledge. Of particular interest is a class of autocovariances, named stable spline kernels [3, 6, 34], encoding system exponential stability. More precisely, the associated Gaussian process is the $q$-fold integration of white noise subject to an exponential time transformation. The derivation of the first-order stable spline kernel via “deterministic” arguments can be found in in [7]. See also [8] for a survey and [26, 27, 28, 35, 36, 31] for even more sophisticated models. The connection with the concept of stable Reproducing kernel Hilbert spaces is also described in [33, 30, 24, 25].

A significant advantage of the stable spline kernel is that it depends on few hyperparameters that can be estimated from data e.g. via marginal likelihood maximization [32, 39, 38]. Once these hyperparameters have been fixed, system identification boils down to a convex problem and the impulse response estimate can be obtained in closed-form. This approach has been proved to be competitive with respect to established identification methods such as PEM and subspace methods.

Within this framework for system identification, the new contributions obtained in this paper are the following ones. First, we use stable spline kernels to derive a new class of radial basis functions (RBF) kernels tailored to describe system components fed with disturbances. Next, we show that the synergic use of this new class and of the stable spline kernels paves the way to a new nonparametric solution of the missing data problem.

It is useful now to recall that identification problems under missing observations have attracted enormous interest in the control as well as in the statistics communities since, at least, the second half of the last century. The number of contributions is hard to be exhaustively surveyed. However, it is interesting to notice that solutions have been developed mainly in parametric settings. For instance, one possible finite-dimensional approach to deal with missing data is to compute the marginal likelihood, i.e. marginalizing the likelihood function over missing data111Note that this depends on the model parameters and should not to be confused with the marginal density in the Bayesian setup which instead depends on the model only through its prior.. Perhaps the first notable attempts in this direction are due to [9, 10], who compute the exact likelihood for an ARMA model under missing observation using the Kalman filter; more efficient algorithms are described in [11] and approximated versions for the AR case are discussed in [12]. However, the marginal likelihood is not straightforward to compute and identification requires the use of costly approaches such as stochastic simulation techniques, e.g. Markov chain Monte Carlo [13], or non linear optimization in high-dimensional domains, likely to converge only to a local minimum. In [14], estimators relying upon the method of moments as well as frequency domain approaches based on the periodogram are also discussed. They can be used as initialization for non-linear optimization procedures, see also [15] and references therein.

A different approach to deal with missing data is to “fill-in” somehow the missing data, a procedure usually referred to as imputation in the statistics literature, and then solve a “complete data” problem. The EM algorithm introduced in the celebrated paper [16] can be utilized to tackle missing data problems. Even though the authors of [16] did not specifically address ARMA estimation, their general framework can be applied to this problem without major difficulties, e.g. see [17] where identification of ARX systems with also missing input data is discussed. Interestingly, the common feature of all the methods so far discussed is that a parametric structure has to be postulated. Thus, model order selection has to be performed on top of the identification procedure using e.g. AIC or BIC [1, 18]. Hence, imputation has to be repeated many times, one for each postulated model order, leading to a high computational complexity.

The novelty in this work is that we also use imputation but adopting stable spline kernels in cooperation with the new class of RBF kernels. The result is a new nonparametric framework where the linear minimum variance estimator of the missing observations can be explicitly worked out. Such estimator contains very few unknown parameters that can be determined from data. In particular, instead of solving many high-dimensional nonlinear optimization problems, as in the parametric setting, imputation is reduced to a single optimization in a two or four-dimensional space at the most. Once the hyperparameter vector is achieved, all the missing observations can be computed in closed form. An additional advantage of the new technique is that it can tackle both discrete-time and continuous-time identification problems. Also, in comparison with [19], the approach here developed allows to use deterministic optimization techniques for performing imputation in place of more costly approaches like stochastic simulation via Markov Chain Monte Carlo.

We also report numerical experiments involving discrete-time ARMAX models. Results reveal that, in many cases of interest, the proposed algorithm may return models that have a better predictive capability on new data than those obtained by standard parametric Prediction Error Methods, also when these latter exploit the full data set.

The paper is organized as follows. Section 2 reports the statement of the problem while our Bayesian model is described in Section 3. In Section 4 we obtain a new class of RBF kernels suitable to describe system components subject to disturbances. In Section 5, we derive the minimum variance linear estimator of the missing observations. This forms the basis of a new nonparametric procedure which we call stable spline imputation. Section 6.1 reports numerical experiments regarding the identification of discrete-time ARMAX models in presence of missing output samples. Conclusions then end the paper while the Appendix gathers the proofs of some technical results.

2 Problem statement

2.1 Notation

We use $\mathcal{X}$ to denote either the set of natural numbers ${\mathbb{N}}$ or the positive real axis ${\mathbb{R}}^{+}$. The $\ell$-th system impulse response is denoted by $f_{\ell}:\mathcal{X}\rightarrow{\mathbb{R}}$ and is assumed to be stable, i.e. absolute summable on $\mathcal{X}$. The $\ell$-th observable and deterministic input is instead denoted by $u_{\ell}$. It is a scalar function defined on ${\mathbb{Z}}$ or ${\mathbb{R}}$, depending if a discrete or a continuous-time problem is considered. All the results present in the paper hold also assuming that the inputs are stochastic processes, independent of the noise entering the system. Without loss of generality, we will consider only MISO systems, fed with $p-1$ observable inputs and one disturbance $e$.

The system output is denoted by $y$ and is a function whose domain can be either ${\mathbb{Z}}$ or ${\mathbb{R}}$. The set of sampling instants where output data are collected is $\{t_{i}\}_{i=1}^{n}$. All the vectors will be column vectors and, in particular, the one containing the observed output measurements is

[TABLE]

The symbol $\otimes$ indicates convolution in discrete or continuous-time. In particular, the expression $\left(q\otimes w\right)(t)$ is the convolution between the functions $q$ and $w$ evaluated at $t$. When involving functions defined only on $\mathcal{X}$, e.g. impulse responses, the operation $\otimes$ is well defined assuming that such functions are null outside $\mathcal{X}$.

2.2 Identification under missing output observations

The measurements model is

[TABLE]

where all the $f_{\ell}$ are unknown and $e$ is the innovation process, i.e. zero-mean white noise of unit variance representing the unobservable system input. For $t\in{\mathbb{Z}}$ or $t\in{\mathbb{R}}$, let

[TABLE]

Then, the predictor model $\tilde{y}(t)$ associated with (1) is

[TABLE]

for a suitable functional $\mathcal{F}$. For instance, $\tilde{y}(t)$ can be the estimator of $y(t+\Delta)$ based on past input-output data up to instant $t$. In particular, in discrete-time, one can recover the familiar one-step ahead linear predictor setting $\Delta=1$ and

[TABLE]

where the $g_{\ell}$ are the predictor impulse responses.

When using PEM, first, the functional $\mathcal{F}$ has to be determined from data, e.g. the predictor can be searched within the class (4), exploiting a finite-dimensional parametrization for $g_{\ell}$ [1] or the nonparametric approach described in [4]. Once the $g_{\ell}$ are known, the system impulse responses $f_{\ell}$ can be computed. However, this paradigm can be utilized only if, for a certain $t$, the set (2) is completely available. If this is not the case, a missing data problem arises.

Our problem is thus to develop a new approach that, starting from $\mathbf{y}_{o}$, provides an estimate of all the output samples in $D_{t}$ that are missing.

2.3 Bayesian Missing data problems in linear models

In the sequel, given a random vector $q$, ${\mathbb{E}}[q]$ denotes its mean while, given $q$ and $w$, $\hat{\mathbb{E}}\left[q|w\right]$ denotes the best linear mean squared estimator of $q$ given $w$. Below, we also use notation of ordinary algebra to handle infinite-dimensional objects.

Let us now consider model (1) and assume, for simplicity of exposition, that $t\in{\mathbb{Z}}$. One has

[TABLE]

where ${\bf y}_{m}$ is the vector with the stacked missing output data, $\Phi_{o},\Phi_{m}$ are suitable matrices containing input (past) data, ${\bf f}$ is an infinite-dimensional vector containing the impulse response coefficients $f_{\ell}(t)$ (as $t\in{\mathbb{Z}}^{+}$ and $\ell\in\{1,2,\ldots,p-1\}$ ). Finally, ${\boldsymbol{\xi}}$ and ${\boldsymbol{\xi}_{m}}$ are column vectors whose entries are suitably stacked “noise” components from $\xi(t)=(e\otimes f_{p})(t)$ stemming from the second term on the right hand side of (1).

Let us now assume that ${\bf f}$ is a zero mean random variable with covariance matrix $\Sigma_{f}$. Assume also that ${\boldsymbol{\xi}}$ and ${\boldsymbol{\xi}}_{m}$ are zero mean random vectors, uncorrelated from ${\bf f}$ and with covariance matrices

[TABLE]

It is a standard result from optimal estimation [20] that the best linear predictor of ${\bf y}_{m}$ given ${\bf y}_{o}$ is given by

[TABLE]

The above equation thus allows to compute in closed form the best linear mean squared estimator of the missing data as a function of the joint (second order) statistics given by the covariances $\Sigma_{f}$, $\mathbf{R}_{m}$ and $\mathbf{R}$.

Eq. (7) will form the basis of our “data imputation” mechanism; in order to do so, we have to specify the statistical description (prior) of both the “deterministic” impulse response ${\bf f}$ as well as of the “noise” part $\xi(t)$, which will be done in the forthcoming Sections. Note that the impulse response vector ${\bf f}$ lives, in principle, in an infinite dimensional space and, as such, the role of the prior is essential for the solution to (5) to be well posed.

3 Bayesian description of the system identification problem

Our main objective is the introduction of a new kernel to model stationary disturbances and the use of a linear minimum variance estimator for the imputation of missing data. For this purpose, a Bayesian description of the system identification problem is now introduced.

Below, given the random vectors $q$ and $w$, we define ${\mathbb{V}}(q,w)={\mathbb{E}}\left[(q-{\mathbb{E}}[q])(w-{\mathbb{E}}[w])^{\top}\right]$ and ${\mathbb{V}}(q)={\mathbb{E}}\left[(q-{\mathbb{E}}[q])(q-{\mathbb{E}}[q])^{\top}\right]$. The Bayesian network in Fig. 1 describes our system identification problem, using as starting point the measurements model (1). Solid and dotted nodes represent, respectively, random and deterministic variables, with arrows to denote stochastic relationships.

Differently from the classical parametric approaches for system identification, the network models the impulse responses $f_{\ell}$ as stochastic processes. In particular, under the framework of Gaussian regression [5], each impulse response is interpreted as the realization of a nonstationary and zero-mean Gaussian process with covariance, also called kernel, defined by

[TABLE]

In (8), $K_{\ell}$ is a symmetric, continuous and positive-definite function $\mathcal{X}\times\mathcal{X}\rightarrow{\mathbb{R}}$ that depends on an unknown hyperparameter vector $\theta$ (also discussed in the next section), while $\lambda\in{\mathbb{R}}$ plays the role of a scale factor common to all the $f_{\ell}$. Looking at the top of the network, one can see that both $\theta$ and $\lambda$ are deterministic quantities. They are contained in a node connected with all the $f_{\ell}$ since it determines the impulse responses statistics.

Notice also that all the nodes $f_{\ell}$ and $e$ are not connected to each other since are all assumed mutually independent random processes. The presence in Fig. 1 of a super-node $\xi$ gathering $f_{p}$ and $e$, and describing the noise part in (5), is instrumental to solving the missing output data problem via a suitable convexification; its role will be elucidated in the next section.

Finally, the node $y$ is the output sequence. The network connections illustrate that it is determined by the impulse responses $f_{\ell}$, by the innovation sequence $e$ and by the deterministic node $u$, which gathers all the observable inputs $u_{\ell}$.

4 Kernels for linear system identification

4.1 Stable spline kernels

In the literature on Gaussian regression, the adopted priors usually reflect only knowledge on the smoothness of the unknown function. One popular approach is to model it as the $q$-fold integral of white Gaussian noise. The resulting covariance becomes proportional to

[TABLE]

where

[TABLE]

This class contains the so-called spline kernels and underlies the Bayesian interpretation of $q$-th order smoothing splines, see [21] for details. One can see that the kernel $W$ does not account for impulse response stability: the variance of $f$ increases as time progresses.

In [3, 6] new kernels for linear system identification have been introduced to include the knowledge on smoothness and exponential BIBO stability. This is obtained by an exponential time transformation regulated by the hyperparameter $\beta>0$ which establishes how fast the variance of $f$ goes to zero. This leads to the following class of kernels parametrized by $q$:

[TABLE]

The choice $q=1$ leads to $K(s,t)=e^{-\beta\max{(s,t)}}$. When $q=2$ one instead obtains the stable spline kernel originally introduced in [3], i.e.

[TABLE]

4.2 RBF kernels for linear system identification

According to our Bayesian paradigm, the impulse responses $f_{\ell}$ are modeled as Gaussian processes with covariance proportional to the stable spline kernel (11). Now, let us focus on the last component of the model (1). It involves $f_{p}$ and the noise $e$, and is given by

[TABLE]

It is easy to see that $\xi$ is a (non Gaussian) zero-mean stationary stochastic process. As already mentioned at the beginning of the previous section, our aim is to define linear minimum variance estimators of the missing data. Hence, now we just need the second-order statistics of $\xi$. They are derived in the following proposition whose proof is not reported since it relies upon simple computations.

Proposition 1**.**

Consider the Bayesian network in Fig. 1. Let $f_{p}$ be a continuous-time zero-mean Gaussian process on ${\mathbb{R}}^{+}$ with covariance $\lambda K$, where $K$ is the stable spline kernel (10). It comes that $\xi$ in (12) is a zero-mean stationary stochastic process on ${\mathbb{R}}$, with covariance $\lambda R$ where

[TABLE]

In particular, for $q=1$ one has

[TABLE]

while $q=2$ leads to

[TABLE]

Eq. (13) thus provides a new class of RBF kernels for identification, useful to describe system components subject to disturbances. Remarkably, the simplest element, obtained with $q=1$ and reported in (14), corresponds to the familiar Laplace kernel.

4.3 Stable spline and RBF kernels in discrete-time

It is straightforward to extend the results obtained in the previous two subsections to the discrete-time context. For what concerns the stable spline kernels (10), one can just consider their sampled versions

[TABLE]

Then, starting from (16), the discrete-time versions of the RBF kernels (13) become

[TABLE]

In particular, for $q=1$ one has

[TABLE]

while $q=2$, that corresponds to using (11), leads to

[TABLE]

4.4 Enriching the kernels and the ARMAX case

In some circumstances, it can be useful to add to the stable spline kernels (10) and (16) some components able to capture dynamics which are hardly represented by smooth processes, e.g. high-frequency poles. As described in [4], this goal can be obtained modeling each $f_{\ell}$ as

[TABLE]

where $g_{\ell}$ is a zero-mean Gaussian process, with covariance proportional to the stable spline kernel, while $h_{\ell}$ is a low-order parametric impulse response. Notice also that the definition of the new RBF kernels remains the same, except that $K$ in (13) and (17) has to be replaced by the stable spline kernel enriched with the parametric component.

For discrete-time ARMAX models, it is useful to set $h=h_{1}=\ldots=h_{p}$. In what follows, the zeta transform of $h$ is then given by

[TABLE]

The overall model so depends on four hyper-parameters: $\phi$ and $\varphi$ which carry the information on the poles common to the impulse responses, the variance decay rate $\beta$ and the scale factor $\lambda$.

5 Stable spline imputation

5.1 Notation

In order to introduce the new imputation procedure, first we need to set up some additional notation.

Given the RBF kernel $R$ and the sampling instants $\{t_{i}\}_{i=1}^{n}$, $\mathbf{R}\in\mathbb{R}^{n\times n}$ is a positive semidefinite matrix, already introduced in (6), that we call RBF kernel matrix, whose $(i,j)$ entry is

[TABLE]

Given $K$ and the observable system inputs, the output kernel $P:\mathcal{X}\times\mathcal{X}\rightarrow{\mathbb{R}}$ is defined, for every $x,y\in\mathcal{X}$, by

[TABLE]

where $h_{x}(t)$ is a function, parametrized by $x$, defined $\forall t$ by

[TABLE]

Notice that, when performing the outer convolution, $h_{x}$ is thought of as a function parametrized by $x$. Furthermore, $\mathbf{P}\in\mathbb{R}^{n\times n}$ denotes the output kernel matrix whose $(i,j)$ entry is

[TABLE]

Remark 2**.**

When working in discrete-time, i.e. $\mathcal{X}={\mathbb{N}}$, $P(i,j)$ admits a simple expression using a matrix vector notation. In fact, let $\mathbf{K}^{\infty}_{\ell}$ denote the infinite-dimensional kernel matrix whose $(i,j)$ entry is

[TABLE]

Then, one obtains

[TABLE]

where

[TABLE]

Above, as in Section 2.3, notation of ordinary algebra has been adopted to handle infinite-dimensional objects.

5.2 Minimum variance linear estimator of the missing data

The next proposition provides the minimum variance linear estimator and the posterior covariance of the system output for known $\theta$ (see Appendix for the proof).

Proposition 3**.**

Consider the Bayesian model displayed in Fig. 1, where each $f_{\ell}$ is a zero-mean stochastic process of covariance $\lambda K$, with $K$ the stable spline kernel (10). Then, for known $\theta$ and arbitrary time instant $t$, the minimum variance linear estimator of $y(t)$ is

[TABLE]

where the $\hat{c}_{i}$ are the components of the vector $\hat{\mathbf{c}}=\left(\mathbf{P}+\mathbf{R}\right)^{-1}\mathbf{y}_{o}.$ Finally, the posterior covariance of $y(t)$ given $\mathbf{y}_{o}$ is

[TABLE]

*where $a_{t}=\left[P(t,t_{1})+R(t,t_{1})\quad\ldots\quad P(t,t_{n})+R(t,t_{n})\right].$ *

Eq. (24) thus makes available all the estimate of the system output in closed form. One can also see that the system output estimator has the structure of a particular regularization network [22]. It is a sum of $n$ basis functions with expansions coefficients obtained by solving a linear system of equations. Each basis function is the sum of the output kernel section $P(\cdot,t_{i})$, coming from the stable spline kernels convoluted with the system inputs, and of the RBF kernel section $R(\cdot,t_{i})$. Notice also that the estimate does not depend on the scale factor $\lambda$.

5.3 Stable spline imputation

In real applications, the estimator (24) can not be directly applied since it depends on the unknown vector $\theta$ entering the kernels $P$ and $R$. This problem can be faced exploiting the Bayesian framework underlying the stable spline estimator. In particular, $\theta$ can be estimated by optimizing the marginal likelihood, i.e. the joint density of $\mathbf{y}_{o}$ and the impulse responses $f_{\ell}$, where the impulse responses are integrated out. Adopting a Gaussian approximation for the disturbance, then using the same arguments adopted in Appendix of [4], the estimate of $\theta$ becomes

[TABLE]

Our numerical procedure, namely stable spline imputation, able to return any missing output sample is then given below.

Algorithm 4** (Stable Spline Imputation).**

The input to this algorithm includes the system inputs $u_{\ell}$, the sampling instants $\{t_{i}\}_{i=1}^{n}$, and the measurement vector $\mathbf{y}_{o}$. The output of this algorithm is the estimate $\hat{y}(t)$, where $t$ denotes the time instant where the system output needs to be estimated. The steps are as follows:

• •

*Compute the hyperparameter vector via marginal likelihood optimization as described in (26). *

• •

Set the hyperparameter vector $\theta$ to $\hat{\theta}$ and return the following estimate of $y(t)$:

[TABLE]

where the $\hat{c}_{i}$ are the components of the vector

[TABLE]

Usual parametric approaches to imputation must solve several non convex optimization problems, possibly in high-dimensional spaces. This e.g. happens adopting rational transfer functions of unknown order so that it is necessary consider various model parameterizations, possibly involving vectors with tens/hundreds of components. In the stable spline imputation procedure, model order is instead encoded in the hyperparameter vector $\theta$. Hence, only $J$ must be optimized. Any evaluation of such objective requires $O(nq^{2})$ operations using $q$ to denote the number of estimated coefficients of the one-step ahead predictor, e.g. see [8] for details on marginal likelihood computation. Once $\theta$ is known, the solution is then available in closed form. When ARMAX models are e.g. considered, $\theta$ contains only $(\beta,\lambda)$ or, at the most, also the other two parameters $(\phi,\varphi)$ describing two poles. Thus, a two or a four-dimensional space needs to be explored so that, even if the objective is non-convex, grid methods could be also adopted to mitigate the risk of local minima.

6 Numerical experiments

6.1 Identification of discrete-time ARMAX models

Let us now consider a Monte Carlo study of 500 runs. During each run a discrete-time ARMAX linear system with $3$ observable inputs is randomly generated as follows:

• •

each rational transfer function is given the same order, randomly chosen in $\{1,2,...,30\}$;

• •

the polynomials defining the model are generated using the MATLAB function drmodel.m. The system and the predictor poles are restricted to have modulus less than 0.95 with the ratio between the sum of the $\ell_{2}$ norms of $f_{\ell}$, $\ell=1,2,3$, and of $f_{4}$ falling in $[1,5]$ (drmodel.m is repeatedly called at any run until such requirements are fulfilled). The delay between all the inputs and the output is always equal to 1.

Independent realizations of unit variance white noises are used as inputs. Output data are collected after getting rid of initial conditions to define a training set of 300 input-output pairs and a test set of size 1000. To simulate a missing data problem, at every run the training set is reduced by randomly discarding an output value with probability 0.25.

6.2 Performance measures

Two performance measures will be used to compare the performance of different estimators. The first index is the quality in the reconstruction of the missing data contained in the vector ${\bf y}_{m}$. It is given by the Coefficient of Determination $COD_{miss}$, computed at any run as

[TABLE]

where $\|\cdot\|$ is the Euclidean norm and $\bar{{\bf y}}_{m}$ indicates the average value of the components of ${\bf y}_{m}$.

The second index measures the ability of the estimated model in predicting the test set, as a function of the prediction horizon $k$. In particular, we use the $k$-step ahead Coefficient of Determination, denoted by $COD_{k}$ for $k=1,2,\ldots,20$, and computed at any run as follows:

[TABLE]

where ${\bf y}_{test}$ is the vector containing the outputs in the test set, whose sample mean is denoted by $\bar{{\bf y}}_{test}$, while the components of the vector $\hat{\bf y}_{test}^{k}$ are the $k$-step ahead predictions computed using the estimated model. The average of the values of $COD_{k}$ obtained after the 500 Monte Carlo runs is then denoted by $\overline{COD}_{k}$.

6.3 The adopted estimators

The following 6 estimators are used:

• •

PEM+Oracle (missing): This algorithm computes the PEM (ML in the Gaussian case) estimator of the system parameters; the Prediction Error cost is computed using the Kalman filter (see, e.g. [9]). Note that, alternatively, the same result would have been obtained formulating it as a missing data problem by adding the missing observations as unknowns and then using the EM algorithm, which iterates between computing conditional expectation of the conditional log likelihood over the missing data for fixed model parameters and maximizing of the expected log likelihood [16, 17]. The procedure is repeated for all ARMAX models with the three polynomials with the same degree ranging from 1 to 20 (increasing this number does not improve the results); the model order selection is then performed by an oracle which, at each run, maximizes $COD_{k}$ and $COD_{miss}$. Notice that the chosen orders depend on the target, i.e. prediction quality on a certain horizon $k$ or quality in the reconstruction of the missing data. The same order selection procedure is used also when the other PEM-based approaches described below are used. The missing data are then estimated using a Rauch-Tung-Striebel smoother based on the estimated model parameters, which are the maximum a posteriori estimates of the missing observations conditionally on the estimated model and the observed data.

• •

PEM+Oracle (full): the same as above, except that the estimator uses all the 300 identification data and the pem.m function of the MATLAB System Identification Toolbox is used to estimate the model.

• •

PEM+BIC (full): the same as above, except that BIC is used for model order selection.

• •

PEM+AICC (full): the same as above except that the corrected version of Akaike criterion (AICC) [23] is used for model order selection.

• •

SS (full): this is the Stable Spline estimator described in [4] to which the reader is referred for all the details. Here, we just recall that the stable spline kernel of order $q=2$ enriched with a parametric part describing two poles is used. Hence, the dimension of the hyperparameter vector is 4 and its components are estimated via marginal likelihood optimization. For computational reasons, the number of estimated coefficients of the one-step predictor impulse responses is set to 100. To form the regression matrices, entries depending on samples values at time instants $t<t_{1}$ are set to zero, even if data are not generated starting from null initial conditions.

• •

SS imputation+SS: this estimator has to identify the system using the reduced training set. First, the stable spline imputation procedure defined by Algorithm 4 is used to recover the missing output values. Then, the ARMAX model is estimated using the Stable Spline estimator fed with the union of the available measurements and the estimated outputs as returned by Algorithm 4. Also in this case, regression matrices’ entries depending on samples values at time instants $t<t_{1}$ are set to zero.

The information that the delay between all the inputs and the output is equal to 1 is provided to all the system identification algorithms listed above. Notice that all the estimators exploit the full data set of 300 samples, except PEM+Oracle (missing), and SS imputation+SS which, on average, use only 225 output measurements.

6.4 SS imputation+SS vs PEM+Oracle (missing)

We start comparing the performance of PEM+Oracle (missing), and SS imputation+SS which are the two estimators which have to handle missing data situations.

The top panel of Fig. 2 displays results regarding the prediction capability of the estimated models. In particular, the abscissae and the ordinates of the 2000 points contained in the figure correspond to the values of $\{COD_{k}\}_{k=1}^{20}$ returned, respectively, by PEM+Oracle (missing) and SS imputation+SS after the first 100 runs. It turns out the the predictive performance of SS imputation+SS is superior to that of the oracle-based procedure in almost $87\%$ of the cases. This result is remarkable since the oracle uses additional information not available to SS imputation+SS to select that model order (function of $k$) which maximizes $COD_{k}$ or $COD_{miss}$.

The bottom panel of Fig. 2 compares the quality in the reconstruction of ${\bf{y}}_{m}$ by reporting the MATLAB boxplots of the values of $COD_{miss}$ returned by the two estimators after the 500 runs. One can see that the performance of the stable spline imputation (which is implementable in real applications) is very similar to that of PEM+Oracle (missing) (which is not implementable in practice).

6.5 SS imputation+SS vs estimators using the full data set

We now compare the performance of SS imputation+SS with that of all the other estimators exploiting the full data set. The top panel of Fig. 3 displays $\overline{COD}_{k}$ achieved by the 5 estimators after the 500 runs, as a function of the prediction horizon. The bottom panel display the boxplots of the 500 values of $COD_{5}$. The mean performance of the new estimator SS imputation+SS, is comparable to that of SS (full) which provides result similar of those of PEM+Oracle (full).

6.6 Variants of the experiment

We have also considered two variants of the experiment. In the first one, in place of white noises, the observable system inputs are given by low pass signals. In particular, the inputs are independent realizations of white noises filtered by a strictly proper second-order rational transfer function randomly generated at every run (the same mechanism used to generate the system impulse responses is used). In the second one, the system and the predictor poles are restricted to have modulus less than 0.999 (in place of 0.95) and, in addition, the unknown system impulse responses have been enriched by adding to them a couple of lowly damped poles. In particular, each transfer function is that obtained by the MATLAB generator multiplied by

[TABLE]

where $b=0.999$ while $a$ is, at every run, a different realization from a uniform distribution on $[-1,1]$. The left and right panels of Fig. 4 displays $\overline{COD}_{k}$ following the same rationale adopted in Fig. 2. As a matter of fact, in both the variants, the stable spline estimators still outperform the classical system identification approaches also when the latter exploit the full data set.

7 Conclusions

In this paper we have shown that stable spline kernels can be used also to derive a new class of RBF covariances useful to model that part of the system output due to disturbances. The stable spline and the new RBF kernels lead to a new solution of the missing data problem based on a new imputation procedure, namely stable spline imputation. It returns all the missing output samples just solving one low-dimensional optimization problem.

The new technique has been used to identify discrete-time ARMAX models under missing data. In many cases of interest the new stable spline imputation followed by the stable spline estimator developed in [4] may return models more predictive than those obtained by standard parametric PEM, also when the latter have access to the full data set.

Appendix: Proof of Proposition 3

The proof relies upon well known results regarding the estimation of stochastic processes, e.g. see [20]. The minimum variance linear estimator of $y(t)$ given $\mathbf{y}_{o}$ is

[TABLE]

Now, recall that the innovation $e$ and the $\{f_{\ell}\}_{\ell=1}^{p}$ are all assumed mutually independent. This implies also the independence of $\xi$ and $\{f_{\ell}\}_{\ell=1}^{p-1}$. Thus, we obtain ${\mathbb{V}}(\mathbf{y}_{o})=\lambda\left(\mathbf{P}+\mathbf{R}\right)$ and one also has

[TABLE]

so that (24) is obtained. Similar arguments lead to (25) and this completes the proof.

of Proposition 3.