The Harmonic Analysis of Kernel Functions

TL;DR

This paper introduces a harmonic analysis approach to kernel functions in Bayesian system identification, aiding kernel design and reducing computational complexity.

Contribution

It presents a novel harmonic analysis method for non-stationary Gaussian process kernels, improving kernel design and computational efficiency.

Findings

Harmonic analysis provides insights into kernel structure.

The approach enables efficient kernel approximation.

Results suggest improved computational performance.

Abstract

Kernel-based methods have been recently introduced for linear system identification as an alternative to parametric prediction error methods. Adopting the Bayesian perspective, the impulse response is modeled as a non-stationary Gaussian process with zero mean and with a certain kernel (i.e. covariance) function. Choosing the kernel is one of the most challenging and important issues. In the present paper we introduce the harmonic analysis of this non-stationary process, and argue that this is an important tool which helps in designing such kernel. Furthermore, this analysis suggests also an effective way to approximate the kernel, which allows to reduce the computational burden of the identification procedure.

Figures

Equations

$$y(t)=G(z)u(t)+e(t),t\in \mathbb{N}$$

$$G(z)$$

$${\mathcal{Z}}^N=\{u(1),y(1)\dots u(N),y(N)\}.$$

$$g_a(t)=\sum _{l=1}^N|c_l|e^{{\alpha }_{l}t}\cos ({\omega }_{l}t+\angle c_l),t\in {\mathbb{R}}_{+}$$

$$\mathbb{E}[c_l\overline{c_{l^{\prime }}}]={\varphi }_l{\delta }_{l-l^{\prime }}$$

$$\mathbb{E}[c_l{c}_{l^{\prime }}]=0$$

$$g_a(t)=\sum _{i=1}^{N_{\alpha }}\sum _{k=1}^{N_{\omega }}|c_{ik}|e^{{\alpha }_{i}t}\cos ({\omega }_{k}t+\angle c_{ik})$$

$$\mathbb{E}[c_{ik}\overline{c_{i^{\prime }{k}^{\prime }}}]={\varphi }_{ik}{\delta }_{i-i^{\prime }}{\delta }_{k-k^{\prime }}$$

$$\mathbb{E}[c_{ik}{c}_{i^{\prime }{k}^{\prime }}]=0$$

$$g(t)={\int }_{-\infty }^0{\int }_{-\infty }^{\infty }|c(\alpha ,\omega )|e^{\alpha t}\cos (\omega t+\angle c(\alpha ,\omega ))\mathrm{d}\omega \mathrm{d}\alpha$$

$$\mathbb{E}[c(\alpha ,\omega )\overline{c({\alpha }^{\prime },{\omega }^{\prime })}]=\varphi (\alpha ,\omega )\delta (\alpha -{\alpha }^{\prime })\delta (\omega -{\omega }^{\prime })$$

$$\mathbb{E}[c(\alpha ,\omega )c({\alpha }^{\prime },{\omega }^{\prime })]=0$$

$$K(t,s)=\frac{1}{2}{\int }_{-\infty }^0{\int }_{-\infty }^{\infty }\varphi (\alpha ,\omega )e^{\alpha (t+s)}\cos (\omega (t-s))\mathrm{d}\omega \mathrm{d}\alpha .$$

$$\mathbb{E}\left[g(t{)}^2\right]=K(t,t)=\frac{1}{2}{\int }_{-\infty }^0{\int }_{-\infty }^{\infty }\varphi (\alpha ,\omega )e^{2\alpha t}\mathrm{d}\alpha \mathrm{d}\omega .$$

$$\mathbb{E}\left[g(0{)}^2\right]=K(t,t)=\frac{1}{2}{\int }_{-\infty }^0{\int }_{-\infty }^{\infty }\varphi (\alpha ,\omega )\mathrm{d}\alpha \mathrm{d}\omega .$$

$$K(t,s)=K_1(t-s)=\frac{1}{2}{\int }_{-\infty }^{\infty }{\varphi }_1(\omega )\cos (\omega (t-s))\mathrm{d}\omega$$

$$g(t)={\int }_{-\infty }^{\infty }|c(\omega )|\cos (\omega t+\angle c(\omega ))\mathrm{d}\omega$$

$$g(t)={\int }_{-\infty }^{\infty }{e}^{j\omega t}\mathrm{d}C(\omega ).$$

$${\varphi }_1(\omega )={\tilde{\varphi }}_1(\omega )+{\tilde{\varphi }}_1(-\omega )$$

$${\tilde{\varphi }}_L(\omega )$$

$${\tilde{\varphi }}_C(\omega )$$

$${\tilde{\varphi }}_G(\omega )$$

$$K_L(t-s)$$

$$K_C(t-s)$$

$$K_G(t-s)$$

$$K(t,s)=K_1(t-s)K_2\left(t+s\right)$$

$$K_2\left(t+s\right)={\int }_{-\infty }^0{\varphi }_2(\alpha )e^{\alpha (t+s)}\mathrm{d}\alpha .$$

$$K(t,s)=e^{{\alpha }_0(t+s)}{K}_1(t-s).$$

$$\varphi (\alpha ,\omega )=2\frac{-\alpha }{\pi [{\omega }^2+{\alpha }^2]}{1}_{[{\alpha }_m/2,{\alpha }_M/2]}(\alpha )$$

$$\displaystyle\mathbf{1}_{[\alpha_{m}/2,\alpha_{M}/2]}(\alpha)=\left\{\begin{array}[]{cc}1,&\alpha\in[\alpha_{m}/2,\alpha_{M}/2]\\ 0,&\hbox{otherwise}\end{array}\right.$$

$$K(t,s)=$$

$$\times {\int }_{-\infty }^{\infty }\frac{-\alpha }{\pi [{\omega }^2+{\alpha }^2]}\cos (\omega (t-s))\mathrm{d}\omega \mathrm{d}\alpha$$

$$={\int }_{{\alpha }_m/2}^{{\alpha }_M/2}{e}^{\alpha (t+s)}{e}^{\alpha |t-s|}\mathrm{d}\alpha$$

$$={\int }_{{\alpha }_m/2}^{{\alpha }_M/2}{e}^{2\alpha \max \{t,s\}}\mathrm{d}\alpha$$

$$=\frac{e^{{\alpha }_M\max \{t,s\}}-e^{{\alpha }_m\max \{t,s\}}}{2\max \{t,s\}}$$

$$\varphi (\alpha ,\omega )={\varphi }_1(\omega ;\alpha )1_{[{\alpha }_m,{\alpha }_M]}(\alpha ),$$

Full text

††thanks: This work has been partially supported by the FIRB project “Learning meets time” (RBFR12M3AC) funded by MIUR.

,

1 Introduction

Building upon the theory of reproducing kernel Hilbert spaces and statistical learning, kernel-based methods for linear system identification have been recently introduced in the system identification literature, see Pillonetto & De Nicolao (2010); Pillonetto et al. (2011); Chiuso & Pillonetto (2012); Chen et al. (2012); Pillonetto et al. (2014); Zorzi & Chiuso (2017, 2015); Chiuso (2016); Fraccaroli et al. (2015).

These methods, framed in the context of Prediction Error Minimization, differ from classical parametric methods Ljung (1999); Söderström & Stoica (1989), in that models are searched for in possibly infinite dimensional model classes, described by Reproducing Kernel Hilbert Spaces (RKHS). Equivalently, in a Bayesian framework, models are described assigning as prior a Gaussian distribution; estimation is then performed following the prescription of Bayesian Statistics, combining the “prior” information with the data in the posteriors distribution. Choosing the covariance function of the prior distribution, or equivalently the Kernel defining the RKHS, is one of the most challenging and important issues. For instance the prior distribution could reflect the fact that the system is Bounded Input Bounded Output (BIBO) stable, its impulse response possibly smooth and so on (Pillonetto et al., 2014).

Within this framework, the purpose of this paper is to discuss the properties of certain kernel choices from the point of view of Harmonic Analysis of stationary processes. The latter is well defined for stationary processes (Lindquist & Picci, 2015, Chapter 3). In particular, it defines as Power Spectral Density (PSD) the function describing how the statistical power is distributed over the frequency domain. In this paper, we extend this analysis for a particular class of non-stationary processes modeling impulse responses of marginally stable systems. Accordingly, we define as Generalized Power Spectral Density (GPSD) the function describing how the statistical power is distributed over the decay rate-frequency domain.

Under the assumption that the prior density is Gaussian, the probability density function (PDF) of the prior is linked to the GPSD. The main difference is that while the former is defined over an infinite dimensional space (the underlying RKHS ${\mathcal{H}}_{K}$), the latter is defined over the bidimensional decay rate-frequency space. As a consequence, the latter is simple to depict but also to interpret from an engineering point of view. We show experimentally that, over the class of second-order linear systems, the two provide similar information. This class is important because: 1) it contains the simplest systems that exhibit oscillations and overshoot; 2) second order systems are building block of higher order systems and, as such, understanding second order systems helps understanding higher ones. Furthermore, for a special class of exponentially convex locally stationary processes (ECLS) typically used in system identification (Chen & Ljung, 2015b, a), it is possible to provide (i) a clear link between the GPSD and the Fourier transform of the exponentially modulated sample trajectory and (ii) characterize the posterior mean in terms of the GPSD. As a consequence, it is possible to outline a simple procedure for the design of the kernel, through the GPSD.

Another important aspect in kernel-based methods is to reduce the computational burden (Chen & Ljung, 2013). This task can be accomplished by approximating the kernel functions through eigen-decomposition (Carli et al., 2012) or random features, (Rahimi & Recht, 2007) techniques. However, these methods can be applied only to special kernel functions. We show that, the GPSD provides a general procedure to approximate a wide class kernel functions.

The outline of the paper follows. In Section 2 we review the Gaussian process regression framework used for kernel-based methods. In Section 3 we present the harmonic representation of the kernel function of continuos time non-stationary processes modeling impulse responses of marginally stable systems and in Section 4 the corresponding discrete time version. In Section 5 we show the relation between the GPSD and the probability density function of the prior over the class of second-order linear systems. In Section 6 we characterize the posterior mean in terms of the GPSD for a special class of ECLS processes. Section 7 regards kernel approximation using the GPSD. Finally, conclusions are drawn in Section 8. In order to streamline the presentation all proofs are deferred to the Appendix.

Notation: $\mathbb{N}$ denotes the set of the natural numbers, $\mathbb{Z}$ is the set of the integer numbers, $\mathbb{R}_{+}$ the set of the nonnegative real numbers, $\mathbb{R}_{-}$ the set of the negative real numbers, and $\mathbb{C}$ is the set of the complex numbers. Given $x\in\mathbb{C}$, $|x|$ denotes its absolute value, $\angle{x}$ denotes its phase and $\overline{x}$ denotes its conjugate. Given $A\in\mathbb{C}^{m\times n}$, $A^{*}$ denotes its transposed conjugate. $\mathbb{E}[\cdot]$ denotes the expectation operator.

2 System Identification and Gaussian Process Regression

For convenience in what follows we consider a discrete time, stable and linear time-invariant (LTI) single input-single output (SISO) in OE form:

[TABLE]

where $z^{-1}$ denotes the backward shift operator; $u(t)$ is the input; $y(t)$ is the noisy output; $e(t)$ is zero mean white noise, that is $\mathbb{E}[e(t)e(s)]=\sigma^{2}\delta_{t-s}$ where $\delta_{k}$ denotes the Kronecker delta function. The transfer function $G(z)$ is stable and strictly causal, i.e. $G(\infty)=0$. Expanding $G(z)$ in $z^{-1}$ we obtain the impulse response of the linear system

[TABLE]

The system identification problem can be frased as that of estimating the impulse response $\{g(t)\}_{t\in\mathbb{N}}$, from the given data record

[TABLE]

In the Gaussian process regression framework, $g(t)$ is modeled as a zero-mean discrete time Gaussian process with kernel (covariance) function $K(t,s):=\mathbb{E}[g(t)g(s)]$. The minimum variance estimator of $g(t)$ is given by its posterior mean given $\mathcal{Z}^{N}$, (Pillonetto & De Nicolao, 2010). It is clear that the posterior highly depends on the kernel functions. Accordingly, the most challenging part of this system identification procedure is to design $K$ so that the posterior has some desired properties.

Similarly, in the continuous time case, $\{g(t)\}_{t\in\mathbb{R}_{+}}$ is a zero-mean continuous time Gaussian process with kernel function $K(t,s):=\mathbb{E}[g(t)g(s)]$, with $t,s\in\mathbb{R}_{+}$. In what follows Gaussian processes (both discrete time and continuous time) are always understood with zero-mean.

3 Harmonic Analysis: continuous time case

It is well known that the impulse response of a finite dimensional LTI stable (or marginally stable) system can be written as a linear combination of decaying sinusoids111For simplicity we exclude here the case of multiple eigenvalues. (i.e. modes)

[TABLE]

where $\alpha_{l}\in\mathbb{R}_{-}\cup\{0\}$ and $\omega_{l}\in\mathbb{R}$ are, respectively, the decay rate and the angular frequency of the $l$-th damped oscillation, and $c_{l}\in\mathbb{C}$. Adopting the Bayesian perspective, $g(t)$ is modeled as a Gaussian process where the coefficients $c_{l}$ are zero mean complex Gaussian random variables such that

[TABLE]

with $l,l^{\prime}=1\ldots N$ and $\phi_{l}\geq 0$. For convenience, we rewrite (2) as

[TABLE]

that is $(\alpha_{i},\omega_{k})$ belongs to a $N_{\alpha}\times N_{\omega}$ grid contained in $\mathbb{R}_{-}\cup\{0\}\times\mathbb{R}$ and $c_{ik}$ is a complex Gaussian random variable such that

[TABLE]

with $\phi_{ik}\geq 0$. It is then natural to generalize (3) as an infinite “dense” sum of decaying sinusoids 222Strictly speaking the integral over $\alpha$ should be understood in an open interval of the form $(-\infty,\epsilon)$ with $\epsilon>0$ and $\epsilon\rightarrow 0$.:

[TABLE]

where $c(\alpha,\omega)$ is a bidimensional complex Gaussian process333To be precise, we should work with the “generalized spectral measure” $C(\omega,\alpha)$, a Gaussian process with orthogonal increments; formally $dC(\omega,\alpha)=c(\omega,\alpha)d\omega d\alpha$., hereafter called generalized Fourier transform of $g(t)$, such that

[TABLE]

where $\phi(\alpha,\omega)$ is a nonnegative function on $\mathbb{R}_{-}\cup\{0\}\times\mathbb{R}$ such that $\phi(\alpha,\omega)=\phi(\alpha,-\omega)$ and $\delta(\cdot)$ denotes the Dirac delta function.

Proposition 1**.**

Let $K(t,s)$ be the kernel function of $g(t)$ in (4) then,

[TABLE]

Formula (5) is the harmonic representation of the covariance function of the non-stationary process (4). We refer to $\phi(\alpha,\omega)$ as generalized power spectral density (GPSD) of $g(t)$. The latter describes how the “statistical power” of $g(t)$ (which depends on $t$) is distributed over the decay rate-angular frequency space $\mathbb{R}_{-}\cup\{0\}\times\mathbb{R}$ according to:

[TABLE]

In particular, when $t=0$, the exponential term $e^{2\alpha t}$ disappears, so that the statistical power of $g(0)$ is obtained as in the stationary case:

[TABLE]

It is worth noting that the GPSD can be understood as a function in $\mathbb{C}$, that is $\phi(s)$ with $s=\alpha+j\omega$, and its support is the left half-plane. Now we show that the harmonic representation (5) describes many kernel functions used in system identification.

3.1 Stationary kernels

The special case of stationary processes is recaptured when $\phi(\alpha,\omega)=\delta(\alpha)\phi_{1}(\omega)$, with $\phi_{1}(\omega)=\phi_{1}(-\omega)$ nonnegative function; in fact, under this assumption, we have

[TABLE]

which is a stationary kernel and $\phi_{1}(\omega)$ is the corresponding power spectral density (PSD). Note that, (7) is the usual harmonic representation of a stationary covariance function which is also exploited in spectral estimation problems (Zorzi, 2015b, a). In this case the stationary process $g(t)$ is an infinite “dense” sum of sinusoids

[TABLE]

which, more rigorously, should be written in terms of the spectral measure $C(\omega)$

[TABLE]

Since $\phi_{1}(\omega)$ is an even function, we can rewrite it as

[TABLE]

with $\tilde{\phi}_{1}(\omega)$ nonnegative function. For instance, choosing $\tilde{\phi}_{1}(\omega)$ as one of the following:

[TABLE]

we obtain

[TABLE]

where $\omega_{0}$ denotes the angular frequency for which $\tilde{\phi}$ takes the maximum and $\beta$ is proportional to the bandwidth. Setting $\omega_{0}=0$ we obtain, respectively, the Laplacian, Cauchy and Gaussian kernel (Rasmussen & Williams, 2006). In particular, the latter is widely used in robotics for the identification of the inverse dynamic (Romeres et al., 2016).

3.2 Exponential Convex Locally Stationary (ECLS) kernels

A generalization of stationary kernels, introduced by Silverman (1957), is the so-called class of Exponentially Convex Locally Stationary (ECLS) kernels; this is obtained postulating a separable structure for the GPSD $\phi(\alpha,\omega)=\phi_{1}(\omega)\phi_{2}(\alpha)$, with $\phi_{1}(\omega)$ and $\phi_{2}(\alpha)$ nonnegative functions, so that the kernel $K(t,s)$ inherits the decomposition

[TABLE]

where $K_{1}$ has been defined in (7) and

[TABLE]

In the case that $\phi_{2}\left(\alpha\right)=\delta(\alpha-\alpha_{0})$ with $\alpha_{0}\in\mathbb{R}_{-}$ we obtain the ECLS kernel (Chen & Ljung, 2015b, 2016) :

[TABLE]

Specializing $K_{1}(t-s)=K_{L}(t-s)$ with $\omega_{0}=0$, in (9) we obtain the so-called diagonal-correlated (DC) kernel. Furthermore, for suitable choices of $\phi_{1}(\omega)$ in (9), one can obtain the stable-spline (SS), the diagonal (DI) and the tuned-correlated kernel (TC), see Chen et al. (2012).

3.3 Integrated kernels

Consider the GPSD

[TABLE]

where

[TABLE]

with $\alpha_{m}<\alpha_{M}<0$. Then, the corresponding kernel function is

[TABLE]

which is similar to the integrated TC kernel introduced in Pillonetto et al. (2016) 444See Remark 3 for more details.. In general, taking

[TABLE]

where we made explicit the dependence of $\phi_{1}$ upon $\alpha$, we have

[TABLE]

where $K_{1}(t-s;\alpha)$ is the stationary kernel corresponding to $\phi_{1}(\omega;\alpha)$. Notice that, kernel (14) is obtained by integrating the ECLS kernel $e^{\alpha(t+s)}K_{1}(t-s;\alpha)$ over the interval $[\alpha_{m},\alpha_{M}]$, which justifies the name “integrated”. Another possible way to construct an integrated kernel is choosing

[TABLE]

where $\phi_{1}$ does not depend on $\alpha$. Then, the corresponding kernel is

[TABLE]

which is an ECLS kernel with

[TABLE]

4 Harmonic analysis: discrete time case

Following the same argumentations of Section 3, a Gaussian process describing a discrete time causal impulse response can be understood as

[TABLE]

where $c(\lambda,\vartheta)$ is the generalized Fourier transform of $g(t)$, $\lambda$ is the decay rate and $\vartheta$ is the normalized angular frequency. Moreover, the kernel function of $g(t)$ admits the following harmonic representation

[TABLE]

and $\phi(\lambda,\vartheta)=\phi(\lambda,-\vartheta)$ is the GPSD of $g(t)$. The latter is a nonnegative function over the decay rate-normalized angular frequency space $[0,1]\times[-\pi,\pi]$. Also in this case the GPSD can be understood as a function in $\mathbb{C}$, that is $\phi(z)$ with $z=\lambda e^{j\vartheta}$, and its support is the unit circle.

In system identification, it is usual to design a kernel function for a continuous time Gaussian process $g_{c}(t)$, $t\in\mathbb{R}_{+}$. Then, the “discrete time” kernel is obtained by sampling the “continuous time” kernel with a certain sampling time $T$. The latter corresponds to the discrete time Gaussian process $g_{d}(k)$, $k\in\mathbb{N}$, obtained by sampling $g_{c}(t)$.

Proposition 2**.**

Let $\phi_{c}(\alpha,\omega)$ and $\phi_{d}(\lambda,\vartheta)$ be GPSD of $g_{c}(t)$, $t\in\mathbb{R}_{+}$, and $g_{d}(k)$, $k\in\mathbb{N}$, respectively. Then,

[TABLE]

Accordingly, if the continuous time GPSD is such that

[TABLE]

then its discretized version is such that

[TABLE]

Remark 3**.**

Discretizing (3.3) with $T=1$, we obtain

[TABLE]

with $0<\lambda_{m}<\lambda_{M}<1$ and $t,s\in\mathbb{N}$. However, the integrated TC kernel derived in Pillonetto et al. (2016) is slightly different:

[TABLE]

Indeed, the latter has been derived by discretizing the TC kernel, and then the integration along the decay rate has been performed in the discrete domain. Note that, the integration along the decay rate in $\bar{K}_{iTC}$is uniform, while, in $K_{iTC}$ such integration is warped according to (16).

4.1 Filtered kernels

In order to account for high frequency components of predictor impulse responses, Pillonetto et al. (2011) have introduced a class of priors obtained as filtered versions of stable spline kernels, using second order filters of the form:

[TABLE]

with $|\rho_{0}|<1$ and $\vartheta_{0}\in[-\pi,\pi]$. The latter filter is fed by a Gaussian process $\tilde{g}(t)$ with kernel function $\tilde{K}(t,s)$, for instance in the class of “stable-spline” kernels (see Pillonetto et al. (2010)), which have most of the statistical power (6) concentrated around $\vartheta=0$; in this paper $\tilde{K}$ is chosen as TC kernel. $F(z)$ plays the role of a shaping filter, which concentrates the statistical power of the stationary part around $\vartheta=\vartheta_{0}$. It is not difficult to see that the kernel function of $g(t)$ is

[TABLE]

where $\{f_{l}\}_{l\in\mathbb{N}}$ is the impulse response of the filter $F(z)$, i.e. $F(z)=\sum_{s=0}^{\infty}f_{s}z^{-s}$. Assuming that $\tilde{K}(t,s)$ admits the harmonic representation

[TABLE]

then we have

[TABLE]

where

[TABLE]

Accordingly, (19) is the representation of the kernel function of $g(t)$ according to the basis functions in (20). Accordingly, the effect of filtering is to change the term $\lambda^{t+s}\cos(\vartheta(t-s))$ with (20).

5 Probability Density Function of Gaussian Processes and their GPSD

Given a discrete time Gaussian process $g(t)$ with kernel $K(t,s)$, in this Section we shall study the relation between the associated Gaussian probability distribution and the Generalized Spectral Density introduced in (5). For simplicity we consider the discrete time case with sampling time $T=1$. In what follows, this process is understood as an infinite dimensional vector $g=[\,g(0)\,g(1)\,\ldots\,]^{T}$ and the corresponding kernel as an infinite dimensional matrix $K$. Then, the PDF of the process is

[TABLE]

where “$\propto$” stands for “proportional to”. Practically one can evaluate the r.h.s. of (21) with respect to the impulse response $g$, which, with some abuse of notation, represents a realization of the process, truncating it according to its practical length and taking the corresponding finite dimensional sub-matrix555It can be shown that this makes sense provided $g$ belongs to the RKHS ${\mathcal{H}}$ with kernel $K$. in $K$. In this Section we consider a particular class of candidate impulse responses,

[TABLE]

where

[TABLE]

Notice that each model in this class is uniquely characterized by $p$, accordingly this class is isomorphic to the upper part of the open unit circle.

Thus, we consider the PDF (21) conditionally to the event $g\in{\mathcal{G}}_{p}$,

[TABLE]

As we shall see, there is a close connection between the GPSD of the process $g(t)$ and $p(g|g\in{\mathcal{G}}_{p})$, suggesting that indeed GPSDs whose stationary part concentrate energy around specific frequency bands are well suited to describe second order systems with modes in the same band (i.e. the phase of $p$ in (22)). Similarly, the same applies to the decay rate, which relates to $|p|$ in (22).

5.1 ECLS kernels

We consider the ECLS kernels

[TABLE]

where $\alpha_{0}=-0.1$, $\beta=0.1$ and $\omega_{0}=3\pi/5$. In Figure 1 (top) we show the corresponding PSDs of the stationary part.

Then, we discretize these kernels with $T=1$. The corresponding GPSDs have as support a circle with radius $\lambda_{0}=e^{\tilde{\alpha}_{0}T}=0.9$ centered in zero, see left picture of Figure 2. Since condition (17) holds, the shape of the discretized versions essentially reflect the continuous time version in Figure 1 (top).

This means that the statistical power at time $t=0$ is concentrated in a neighborhood of $z_{0}=\lambda_{0}e^{j\vartheta_{0}}$, with $\vartheta_{0}=\omega_{0}$, along the normalized angular frequency domain. In particular, the one of $K^{\hbox{\tiny ECLS}}_{C}$ is more concentrated in $z_{0}$ than the one of $K^{\hbox{\tiny ECLS}}_{L}$, and the latter is more concentrated than the one of $K^{\hbox{\tiny ECLS}}_{G}$. Finally we also consider the filtered kernel with $\rho_{0}=0.93$ and $\vartheta_{0}=3\pi/5$. The corresponding PDFs are depicted in Figure 3(a).

For all these kernels the impulse responses taking high probability have a pole in a neighborhood of $p_{0}=\lambda_{0}e^{j\vartheta_{0}}$. Consistently with the GPSDs, the neighborhood of high-probability of $K_{G}^{ECLS}$ is more spread along $\vartheta$ than the one of $K_{L}^{ECLS}$, and the latter is more spread along $\vartheta$ than the one of $K_{C}^{ECLS}$. It is worth noting that for the filtered kernel there is a neighborhood of high-probability around a pole with phase close to zero. Indeed, process $\tilde{g}(t)$ at the input of $F(z)$ gives high-probability to impulse responses with pole whose phase is close to zero, and $F(z)$ does not sufficiently attenuate these impulse responses at the output.

Decreasing $\beta$ to $0.05$, the statistical power at time $t=0$ is more concentrated around $z_{0}$ along the normalized angular frequency domain, see Figure 1 (bottom). The corresponding PDFs are depicted in Figure 3(b). Consistently with the GPSDs, the neighborhood of high-probability is more squeezed along $\vartheta$. Regarding the filtered kernel, $\rho$ has been increased to $0.96$; then the bandwidth of $F(z)$ is decreased. As a consequence, the neighborhood of high-probability around the pole with phase close to zero disappeared because $F(z)$ now drastically attenuate these impulse responses at the output.

5.2 Integrated kernels

We consider the integrated versions of (5.1)

[TABLE]

where $\mathcal{T}\in\{L,C,G\}$, $\alpha_{m}=-0.92$, $\alpha_{M}=-0.1$, $\beta=0.05$ and $\omega_{0}=3\pi/5$. Similarly as before, we discretize these kernels with $T=1$. The corresponding GPSDs have as support an annulus inside the unit circle, where the lower radius is $\lambda_{m}=e^{\alpha_{m}}=0.4$ and the upper radius $\lambda_{M}=e^{\alpha_{M}}=0.9$, see Figure 2 (right). If we develop any circle in the support, then we find the function of Figure 1 (bottom) up to a scaling factor. Therefore, the statistical power at time $t=0$ is concentrated in a neighborhood of $z_{0}$ which now is spread also along the decay rate interval $[\lambda_{m},\lambda_{M}]$. The corresponding PDFs are depicted in Figure 3 (c). Consistently with the GPSDs, the PDFs are such that the neighborhood of high-probability probability around $p_{0}$ is now also spread over the decay rate domain. Note that, the integrated version of the filtered kernel is obtained by feeding $F(z)$, here $\rho_{0}=0.96$, with process $\tilde{g}(t)$ with kernel the integrated TC kernel (18).

6 GPSD design for a special class of ECLS kernels

In order to study in more detail the relation among the (sample) properties of Bayes estimators, the frequency domain description of the unknown impulse responses (Fourier transform) and the properties of the kernel, we shall now specialize to a particular class of priors on $g(t)$, admitting an ECLS kernel as in (9). Under this restriction, the generalized Fourier transform takes the form

[TABLE]

where $c_{1}(\omega)$ is a complex Gaussian process such that

[TABLE]

and $\phi_{1}(\omega)\geq 0$ denotes its PSD. Indeed,

[TABLE]

and $\phi(\alpha,\omega)=\phi_{1}(\omega)\delta(\alpha-\alpha_{0})$. By (4) and (24), it is not difficult to see that

[TABLE]

where $G_{\alpha_{0}}(\omega):=(c_{1}(\omega)+\overline{c_{1}(-\omega)})/2$. Accordingly, $G_{\alpha_{0}}(\omega)$ is the Fourier transform of $\pi^{-1}e^{-\alpha_{0}t}g(t)$:

[TABLE]

The next Proposition characterizes the posterior mean of $g(t)$ in terms of the PSD $\phi_{1}(\omega)$.

Proposition 4**.**

Consider the continuos time model

[TABLE]

$u(t)$* is the (known) measured input.

Let $y^{N}=[\,y(1)\,\ldots\,y(N)\,]^{T}$ be sampled noisy measurements 666For simplicity, here we assume that the sampling time is $T=1$.*

[TABLE]

where $e(k)$, $k=1,...,N$ are i.i.d. zero mean Gaussian with variance $\sigma^{2}$, independent of $g(t)$. Assuming the prior distribution on $g(t)$ is Gaussian with kernel (9), the posterior mean $\mathbb{E}[g(t)|y^{N}]$ of $g(t)$ given $y^{N}=[\,y(1)\,\ldots\,y(N)\,]^{T}$ satisfies

[TABLE]

where

[TABLE]

is the posterior mean of $G_{\alpha_{0}}(\omega)$ with

[TABLE]

Proposition 4 can be also adapted to the discrete time case as follows:

Proposition 5**.**

Consider a discrete time process $g(t)$ with ECLS kernel $K(t,s)=\lambda_{0}^{t+s}K_{1}(t-s)$, with $0<\lambda_{0}<1$ and $t,s\in\mathbb{N}$. Then,

[TABLE]

Consider the discrete time OE model

[TABLE]

where $u(t)$ is the measured input and $e(t)$ is zero mean white Gaussian noise with variance $\sigma^{2}$, independent of $g(t)$ and assume $\phi_{1}(\vartheta)$ is the PSD of $K_{1}(t-s)$. Then the posterior mean $\mathbb{E}[g(t)|y^{N}]$ of $g(t)$ given $y^{N}=[\,y(1)\,\ldots\,y(N)\,]^{T}$ is given by

[TABLE]

where

[TABLE]

is the posterior mean of $G_{\lambda_{0}}(\vartheta)$ with

[TABLE]

For simplicity in what follows we consider the discrete time case, but the same observations hold for the continuous time case. Proposition 5 shows that the absolute value of the posterior mean of $G_{\lambda_{0}}(\vartheta)$ is proportional to $\phi_{1}(\vartheta)$ (frequency wise). To understand better this fact, we consider a data record $\mathcal{Z}^{N}$ of length $N=500$ generated by the discrete time model (1) with transfer function having poles in $0.936$, $-0.45\pm 0.8$, $-0.25\pm 0.85$ and zeros in $0.16$, $-0.8\pm 0.4$. Since the dominant pole is $0.936$, we take $\lambda_{0}=0.94$.

In this way $\pi^{-1}\lambda_{0}^{-t}g(t)$ admits Fourier transform $G_{\lambda_{0}}(\vartheta)$, see Figure 4.

We consider the posterior mean estimators of $g(t)$ with three different sampled version of kernel (9) with $T=1$. In particular we choose $\alpha_{0}=\log(\lambda_{0})$, and $\phi_{1}(\omega)$ takes one of the following three shapes:

• •

$\phi_{1}(\omega)=\phi_{L0}(\omega)=2\beta/(\pi(\beta^{2}+\omega^{2}))$

• •

$\phi_{1}(\omega)=\phi_{L}(\omega)=\tilde{\phi}_{L}(\omega)+\tilde{\phi}_{L}(-\omega)$ with $\tilde{\phi}_{L}$ as in (3.1)

• •

$\phi_{1}(\omega)=\phi_{M}(\omega)=\phi_{L0}(\omega)+\phi_{L}(\omega)$ (a mixture kernel inspired by (Chen et al., 2014)).

We shall denote, correspondingly, with $g_{L0}(t)$, $g_{L}(t)$ and $g_{M}(t)$ the three estimators.

The hyperparameters of the kernel are estimated from the data by minimizing the negative log-likelihood (e.g. the ones for $g_{M}(t)$ are $\beta$, $\omega_{0}$ and the scaling factor). In Figure 5,

the Fourier transform of the estimates, as well as the absolute errors are depicted. It is clear that $g_{M}(t)$ provides the best approximation of $g(t)$. In particular, compared to $g_{L0}(t)$ and $g_{L}(t)$, it improves the approximation at low ($\vartheta\simeq 0$), medium ($\vartheta\simeq 1$) and high ($\vartheta\simeq\pi$) normalized angular frequencies. In Figure 6,

the PSDs $\phi_{L0}(\vartheta)$, $\phi_{L}(\vartheta)$ and $\phi_{M}(\vartheta)$ of the stationary part of the three discretized kernels are depicted. Note that $\phi_{L0}(\vartheta)$, $\phi_{L}(\vartheta)$ and $\phi_{M}(\vartheta)$ are the periodic repetition (up to a scaling factor) of $\phi_{L0}(\omega)$, $\phi_{L}(\omega)$ and $\phi_{M}(\omega)$, respectively, with period $2\pi$. It can be noticed that, only $\phi_{M}$ follows the shape of $|G_{\lambda_{0}}(\vartheta)|$. This confirms the intuition that, if the PSD of the stationary part of the kernel has a similar shape of $|G_{\lambda_{0}}(\vartheta)|$, then we expect the corresponding estimate of $g(t)$ is good. Hence this provides guidelines as to how the (stationary part of the) kernel should be designed; in particular, it should mimic the frequency response function $\mathrm{G}_{\lambda_{0}}^{\dagger}(\vartheta)$ of $\pi^{-1}\lambda_{0}^{-t}\mathrm{g}^{\dagger}(t)$ where $\mathrm{g}^{\dagger}(t)$ is the “true” impulse response function.

Note that that optimality of the stationary part of the kernel is coupled with the choice of the decay $\lambda_{0}$; in practice, estimating the hyperparameters using the marginal likelihood estimator allows to optimize jointly $\lambda_{0}$ and the stationary part $\phi_{1}(\omega)$.

Note also that, for the discussion above to make sense, $\lambda_{0}^{-t}\mathrm{g}^{\dagger}(t)$ should admit a Fourier transform, which imposes constraints on $\lambda_{0}$. To the purpose of illustration, let us postulate that $\mathrm{g}^{\dagger}(t)=\sum_{k=1}^{\infty}a_{i}e^{\alpha_{i}t}\cos(\delta_{i}t)$ with $\alpha_{i}\in\mathbb{R}_{-}$, i.e. $\mathrm{g}^{\dagger}(t)$ is a sum of damped sinusoids, then (26)-(27) hold if and only if

[TABLE]

For instance, the TC kernel, i.e. stable spline kernel of order one, is defined as $K_{TC}(t,s)=e^{-\gamma\max\{t,s\}}$ with $\gamma>0$ which can be written in the form (9) choosing $\alpha_{0}=-\gamma/2$ and $\phi_{L}(\omega)$ in (3.1) with $\beta=\gamma/2$, therefore condition (32) becomes $\alpha_{max}<-\gamma/2$. Furthermore, the SS kernel, i.e. stable spline kernel of order two, is defined

[TABLE]

we can rewrite it as (9) by choosing $\alpha_{0}=-3\gamma/2$ and

[TABLE]

therefore condition (32) becomes $\alpha_{max}<3\gamma/2$. It is interesting to note that the two conditions above coincide with the conditions derived in Pillonetto et al. (2010) for the posterior mean estimator to be statistically consistent, namely:

[TABLE]

where $m$ is the order of the stable spline kernel. Indeed, if (26)-(27) hold then the hypothesis space is endowed by a probability density which is strictly positive at the “true” system $g^{\dagger}(t)$; under this condition it can be proved that the posterior mean estimator almost surely converges to $\mathrm{g}^{\dagger}(t)$. A similar reasoning can be applied to the discrete time case.

7 Kernel approximation

Kernel approximation is widely used in machine learning and system identification to reduce the computational burden. Next, we show that the GPSD represents a powerful tool for this problem. To this purpose, note that the process $g_{a}(t)$ in (3) can be understood as an approximation of process $g(t)$ in (4). In particular, its kernel function is

[TABLE]

We define $z_{\alpha,\omega}(t)=e^{\alpha t}[\,\cos(\omega t)\;\,\sin(\omega t)]^{\top}$, then

[TABLE]

The computational burden of the identification procedure based on $K_{a}(t,s)$ is related to the number of points $N_{\alpha}$ and $N_{\omega}$, which can thus be chosen to trade off kernel approximation and computational cost. For fixed $N_{\alpha}$ and $N_{\omega}$, the quality of the approximation $K_{a}(t,s)\approx K(t,s)$ depends on the choice of the $\phi_{ik}$’s, which will be now discussed. First of all, let us observe that (6) can be approximated as follows:

[TABLE]

where $\Delta\alpha_{i}=\alpha_{i+1}-\alpha_{i}$, $\Delta\omega_{k}=\omega_{k+1}-\omega_{k}$ with $\alpha_{i+1}>\alpha_{i}$, $\omega_{k+1}>\omega_{k}$. On the other hand, we have

[TABLE]

Matching (34) and (35), we obtain

[TABLE]

An alternative way to approximate $K(t,s)$ takes inspiration from the random features approach for stationary kernels (Rahimi & Recht, 2007). Observe that (5) can be rewritten as

[TABLE]

where $\mathbb{E}[\cdot]$ is the expectation operator taken with respect to the PDF $\tilde{p}(\alpha,\omega)=K(0,0)^{-1}\phi(\alpha,\omega)$. Then, we can approximate (36) with

[TABLE]

where $\phi_{ik}$’s are drawn from $\tilde{p}(\alpha,\omega)$ and the constant $\gamma$ satisfies $\gamma=\frac{1}{2}\sum_{i=1}^{N_{\alpha}}\sum_{k=1}^{N_{\omega}}\phi_{ik}$.

7.1 Numerical experiments

Data set. We generate 1000 discrete time SISO systems of order $30$. The poles are randomly generated as follows: 75% have phase randomly generated over an interval of size $\pi/6$ centered in $\vartheta_{0}\sim{\mathcal{U}}[\pi/4,3/4\pi]$ and absolute value $\sim{\mathcal{U}}[0.8,0.95]$; the remaining poles are generated uniformly inside the closed unit disc of radius $[0,0.95]$. For each system a data set of 230 points is obtained feeding the linear system with zero mean, unit variance, white Gaussian noise and corrupting the output with additive zero mean white Gaussian noise so as to guarantee a signal to noise ratio equal to $10$.

Simulation setup and results. We consider model (1) with $G(z)=\sum_{t=1}^{n}g(t)z^{-t}$ where $n=100$ is the practical length. We consider several estimators of $g(t)$ which differ on the choice of the kernel describing the prior distribution on $g(t)$. In particular we shall use the following subscripts:

• •

$L0$ for estimator which uses the kernel $K^{INT}_{L}$ with $\omega_{0}=0$;

• •

$L0A$ as above but with $K^{INT}_{L}$ approximated using (33) with $N_{\omega}=5$ and $N_{\alpha}=3$;

• •

$L$ for estimator which uses the kernel $K^{INT}_{L}$;

• •

$LA$ as above but with $K^{INT}_{L}$ approximated using (33) with $N_{\omega}=5$ and $N_{\alpha}=3$;

The subscripts $G0$, $G0A$, $G$, $GA$ and $C0$, $C0A$, $C$, $CA$ will have the same meaning as above but with respect to kernel $K^{INT}_{G}$ and $K^{INT}_{C}$ respectively. Note that all hyperparameters (e.g. $\alpha_{m}$, $\alpha_{M}$, $\beta$, $\omega_{0}$ and the scaling factor for $K^{INT}_{L}$) are estimated from the data by minimizing the negative log-likelihood. To evaluate the various kernels, the impulse response estimates $\hat{g}(t)$, $t=1\ldots n$, are compared to the true one, i.e. $g(t)$, by the average fit

[TABLE]

The distribution of the fits are shown by box-plots in Figure 7: $L0$ and $C0$ estimators are outperformed by their approximated versions. In the remaining cases the approximated versions provides similar performance to the exact version.

8 Conclusions

In this paper we have introduced the harmonic representation of kernel functions used in system identification. In doing that, we have introduced the GPSD which represents the generalization of the power spectral density used for the harmonic representation of stationary kernels. We have showed by simulation that the GPSD carries similar information that we can find in the PDF of the process over the class of second-order systems. Moreover, we have characterized the posterior mean in terms of the GPSD for a special class of ECLS kernels. Finally, we have showed that the GPSD provides a powerful tool to approximate kernels function, and thus to reduce the computational burden in the system identification procedure.

Appendix

Proof of Proposition 1

Rewriting

[TABLE]

then we have

[TABLE]

∎

Proof of Proposition 2

Let $K(t,s)$ be the kernel of $g_{c}(t)$. Define

[TABLE]

Then, by (5) we have

[TABLE]

accordingly $H_{s}(\omega)$ is the Fourier transform of $h_{s}(t)$. Now, we consider $h_{s}^{d}(k):=h_{s}(kT)$, $k\in\mathbb{N}$, which is the sampled version of $h_{s}(t)$ and $T$ is the sampling time. Therefore,

[TABLE]

where $\vartheta=\omega T$ is the normalized angular frequency and

[TABLE]

Setting $n\in\mathbb{N}$, we have

[TABLE]

where we used the substitution $\lambda=e^{\alpha T}$. Substituting (Proof of Proposition 2) in (38) and in view of (37), we obtain

[TABLE]

In view of the discrete time harmonic representation (15), we conclude that $\phi_{d}$ is the GPSD of the sampled kernel $K(kT,nT)$ corresponding to process $g_{d}$.∎

Proof of Proposition 4

By (26), we trivially have (29). Then,

[TABLE]

where

[TABLE]

By (28), we have

[TABLE]

where we exploited relation (26) and

[TABLE]

Let $k,n=1\ldots N$. It follows that

[TABLE]

where we exploited (6) and the fact that $U_{\alpha_{0},t}(\omega)$ is an Hermitian function. Moreover,

[TABLE]

Substituting (26) in (Proof of Proposition 4) and using (6), we obtain

[TABLE]

where we exploited the fact that $\phi_{1}(\omega)=\phi_{1}(-\omega)$. Since $U_{\alpha_{0}}(\omega)=[\,U_{\alpha_{0},1}(\omega)\,\ldots\,\,U_{\alpha_{0},N}(\omega)]^{T}$, we have that $\mathbb{E}[c_{1}(\omega)(y^{N})^{T}]=\phi_{1}(\omega)U_{\alpha_{0}}(\omega)^{*}/2$ and

[TABLE]

where $V_{\alpha_{0}}$ has been defined in (4). In similar way, it can be proven that

[TABLE]

Finally, substituting (42) and (43) in (Proof of Proposition 4) we obtain (30). ∎

Proof of Proposition 5

The proof is similar to the one of Proposition 4. ∎