# On Kernel Design for Regularized Non-Causal System Identification

Xiaozhu Fang<EMAIL_ADDRESS>Tianshi Chen<EMAIL_ADDRESS>School of Data Science and Shenzhen Research Institute of Big Data, The
Chinese University of Hong Kong, Shenzhen 518172, China

###### Abstract

Through one decade’s development, the kernel-based regularization method (KRM)
has become a complement to the classical maximum likelihood/prediction error
method and an emerging new system identification paradigm. One recent example
is its application in the non-causal system identification, and the key issue
lies in the design and analysis of kernels for non-causal systems. In this
paper, we develop systematic ways to deal with this issue. In particular, we
first introduce the guidelines for kernel design and then extend the system
theoretic framework to design the so-called non-causal simulation-induced
(NCSI) kernel, and we also study its structural properties, including
stability and semiseparability. Finally, we consider some special cases of the
NCSI kernel and show their advantage over the existing kernels through
numerical simulations.

###### keywords:

Non-causal system identification, Kernel-based regularization method, Kernel
design.

††thanks: A prelimnary version of this paper was presented at the 41st Chinese
Control Conference, Heifei, China, July 25-27, 2022. Corresponding author
Tianshi Chen.
This work was funded by NSFC under contract No. 62273287, the Shenzhen Science
and Technology Innovation Council under contract No. JCYJ20220530143418040,
and the Thousand Youth Talents Plan funded by the central government of China.

,

##  1 Introduction

In the control community, system identification is the term for the area of
building up mathematical models of dynamical systems based on the measured
input and output data, and has a history of nearly 60 years [1]. In the first
50 years, the major advance of system identification is on the maximum
likelihood/prediction error method (ML/PEM) (and its asymptotic analysis),
which is often called the classical system identification, e.g., [2, 3, 4, 5].
In the last decade, the major advance is on the kernel-based regularization
method (KRM), e.g., the survey papers [6, 7, 8] and the book [9]. In contrast
with ML/PEM, KRM has the feature that it finds a systematic way to engage the
prior knowledge of the underlying system to be identified in the system
identification loop, in particular in the selection of both the model
structure and the model complexity. The carrier of the prior knowledge is the
so-called kernel, which determines the model structure, and its parameter
(called hyper-parameters) determines the model complexity, which can often be
tuned in a continuous way, e.g., [6, 7, 8, 10]. Hence, the kernel plays a
fundamental role and its design is a key issue. A couple of kernel design
methods and many kernels have been proposed, e.g., [11, 10, 12, 13]. Now KRM
has become not only a complement to the classical ML/PEM [2, 14], but also an
emerging new system identification paradigm, which is often called the
regularized system identification [8, 9].

One recent example is the application of KRM in the non-causal system
identification [15], where the causal tuned-correlated (TC) kernel proposed in
[14] was first extended to the non-causal case, and then KRM was used to
identify the non-causal impulse response and achieved satisfying results. The
success of KRM in [15] is mainly because KRM can engage, through the designed
kernel, the prior knowledge on the non-causal system to be identified, e.g.,
the stability and smoothness of the non-causal impulse response. This success
confirms the efficacy and also motivates us to further explore the potential
of KRM in the non-causal system identification. Then, it is worth noting that
all kernels, except the non-causal TC kernel proposed in [15], are causal
kernels designed for causal systems, and the design and analysis of non-causal
kernels for non-causal systems has not been studied in a systematic way
before.

In this paper, we focus on this problem. First, we introduce the guidelines
for the non-causal kernel design and then extend the system theoretic
framework proposed in [10] to design the so-called non-causal simulation-
induced (NCSI) kernel. In particular, the NCSI kernel employs the
multiplicative uncertainty configuration in robust control, e.g., [16], where
the nominal model is used to embed the prior knowledge, the uncertainty is
modelled by a Gaussian process, and the overall model is simulated with an
impulsive input to get the NCSI kernel. Then, we study the structural
properties of the NCSI kernel, including stability and semiseparability.
Finally, we consider some special parameterizations of the NCSI kernel and
show that they give better model estimates/tracking performance than the ones
proposed in [15] through numerical simulations. In contrast with [15], this
paper has the following contributions:

  1. (i)

we develop a systematic framework to design NCSI kernels for non-causal
systems, and moreover, the prior knowledge embedded in NCSI kernels has clear
physical interpretation;

  2. (ii)

we study the structural properties of NCSI kernels, including stability and
semiseparability;

  3. (iii)

we test the designed NCSI kernels by Monte Carlo simulations with randomly
generated systems and moreover, the simulation results show that the designed
NCSI kernels perform (in terms of average accuracy and robustness)
significantly better.

The remaining part of this paper is organized as follows. In Section 2, we
first introduce some background materials and then the problem statement. In
Section 3.1, we introduce the guidelines for non-causal kernel design. In
Section 3, we design the NCSI kernel followed by some analysis and examples.
In Section 4, we run numerical simulations to illustrate the efficacy of the
designed kernels. Finally, we conclude in Section 5. The proofs of theoretical
results are deferred to the appendix.

##  2 Background and Problem Statement

In this section, we first introduce some background materials and then the
problem statement.

###  2.1 Non-Causal System Identification

We consider linear time-invariant (LTI), discrete-time, non-causal and
bounded-input-bounded-output (BIBO) stable systems described by

$\displaystyle y(t)=G(q)u(t)+v(t),t=1,2,\cdots,N,$ (1)

where $t$ is the time instant, $u(t),y(t)\in\mathbb{R}$ are the input and
output of the system, respectively, $v(t)$ is the measurement noise and
assumed to be a white noise with mean 0 and variance $\sigma^{2}$, $q$ is the
forward-shift operator, i.e., $qu(t)=u(t+1)$, and $G(q)$ the transfer function
of the system. Moreover, $G(q)$ can be represented as

$\displaystyle G(q)=\sum\limits_{k=-\infty}^{\infty}g^{0}(k)q^{-k},$ (2)

where $\\{g^{0}(k)\\}_{k\in\mathbb{Z}}$ is called the non-causal impulse
response of $G(q)$, and clearly, $G(q)$ is BIBO stable if its non-causal
impulse response is absolutely summable. The non-causal system identification
problem is to estimate $\\{g^{0}(k)\\}_{k\in\mathbb{Z}}$ or simply $g^{0}$
below as well as possible based on the input-output data
$\\{u(t),y(t)\\}_{t=1}^{N}$.

It is worth to stress that the non-causal system (1) naturally arises from the
inverse model control, where the true system is unknown, stable, non-minimum
phase and has no zeros on the unit circle, e.g., [15] for more details.
Besides, the non-causal systems also arise from other contexts, e.g., the
continuous-time system identification [17], and unstable system identification
[18].

###  2.2 Kernel-Based Regularization Method

The non-causal system identification problem of (1) can be handled by the
kernel-based regularization method (KRM), which relies on a positive
semidefinite kernel
$k(t,s;\eta):\mathbb{Z}\times\mathbb{Z}\rightarrow\mathbb{R}$ with
$\eta\in\Omega\subset{\mathbb{R}}^{p}$ being the hyper-parameter and
$p\in\mathbb{N}$ being its dimension. In particular, KRM searches for a
regularized least squares (RLS) estimate $\hat{g}^{\text{R}}$ of $g^{0}$ in
the reproducing kernel Hilbert space (RKHS) $\mathcal{H}_{k}$ associated with
$k(t,s;\eta)$:

$\hat{g}^{\text{R}}=\underset{g\in\mathcal{H}_{k}}{\text{argmin}}\sum\limits_{t=1}^{N}(y(t)-\sum\limits_{k=-\infty}^{\infty}g(k)u(t-k))^{2}+\gamma||g||^{2}_{\mathcal{H}_{k}},$
(3)

where $||\cdot||_{\mathcal{H}_{k}}$ is the norm of $\mathcal{H}_{k}$, and
$\gamma>0$ is the regularization parameter.

The design of kernels is a core issue for KRM and is referred to as the
parameterization of the kernel $k(t,s;\eta)$ through the hyper-parameter
$\eta$ by encoding the prior knowledge of the underlying system to be
identified. The design of causal kernels, i.e.,
$k(t,s;\eta):\mathbb{N}\times\mathbb{N}\rightarrow\mathbb{R}$, for causal
impulse responses has been well studied, e.g., [10], and many causal kernels
have been introduced, e.g., the diagonal-correlated (DC) kernel and the
tuned/correlated (TC) kernel [14]:

$\displaystyle k^{\text{DC}}(t,s;\eta)=c\lambda^{(t+s)/2}\rho^{|t-s|},$ (4a)
$\displaystyle\eta=[c,\lambda,\rho],\quad c\geq 0,\quad
0\leq\lambda<1,\quad|\rho|\leq 1,$ $\displaystyle
k^{\text{TC}}(t,s;\eta)=c\min\\{\lambda^{t},\lambda^{s}\\},$ (4b)
$\displaystyle\eta=[c,\lambda],\quad c\geq 0,\quad 0\leq\lambda<1.$

In contrast, the design of non-causal kernels, i.e.,
$k(t,s;\eta):\mathbb{Z}\times\mathbb{Z}\rightarrow\mathbb{R}$, for non-causal
impulse responses has not been studied before until recently in [15]. Therein,
a so-called non-causal TC (NCTC) kernel was introduced and takes the following
form

$\displaystyle k^{\text{NC-TC}}(t,s;\eta)=c\min\\{b(t),b(s))\\},$ (5)
$\displaystyle\qquad b(t)=\left\\{\begin{array}[]{ll}\lambda_{c}^{t},&t\geq
0,\\\ \lambda_{a}^{-t},&t<0,\end{array}\right.$ (8)
$\displaystyle\eta=[c,\lambda_{c},\lambda_{a}],\ c\geq 0,\quad
0\leq\lambda_{c}<1,0\leq\lambda_{a}<1,$

where $c$ is a scale factor, $\lambda_{c},\lambda_{a}$ describe the decay
rates for the causal part, i.e., $t\geq 0$ and the anti-causal part, i.e.,
$t<0$, respectively.

###  2.3 Problem Statement

As shown in [15], KRM with the proposed kernels, e.g., (5), can achieve
satisfying results for the non-causal system identification, which motivates
us to further explore the potential of KRM. In particular, we aim to address
the following questions in relation to the kernel design:

  * 1)

what guidelines should be followed when designing non-causal kernels for non-
causal impulse responses;

  * 2)

how to design non-causal kernels in a systematic way following those
guidelines;

  * 3)

is it possible to design non-causal kernels that give even better performance
than the ones proposed in [15], e.g., (5).

##  3 Non-Causal Kernel Design

The idea is to extend the results for causal kernels to non-causal kernels.

###  3.1 Guidelines for Non-Causal Kernel Design

Following this idea, we first derive the optimal kernel in the non-causal
context and then introduce the guidelines for non-causal kernel design.

######  Proposition 3.1 (Optimal Kernel, [14, 6]).

Let $\bar{g}^{0}$ and $\bar{\hat{g}}$ denote any finite dimensional vector
obtained by sampling $g^{0}(t)$ and its estimate $\hat{g}(t)$ at the same but
arbitrary sampling time instants in $\mathbb{Z}$. The optimal kernel is
defined as

$\displaystyle k^{\text{opt}}(t,s)=g^{0}(t)g^{0}(s),\quad t,s\in\mathbb{Z},$
(9)

which minimizes the MSE matrix

$\displaystyle
MSE(k(t,s))=\mathbb{E}[(\bar{\hat{g}}-\bar{g}^{0})(\bar{\hat{g}}-\bar{g}^{0})^{T}],$
(10)

in the sense that $MSE(k(t,s))-MSE(k^{\text{opt}}(t,s))$ is positive
semidefinite for any positive semidefinite kernel $k(t,s)$, where $\mathbb{E}$
denotes the mathematical expectation.

Similar to the causal case, we have the following guidelines for non-causal
kernel design: first, let the kernel mimic the behavior of (9), and moreover,
the prior knowledge of ${g}^{0}$ should be used in the kernel design [14, 10];
second, let the kernel have some special structure that can ease the
computation of KRM [19, 20]. For non-causal systems, the most common prior
knowledge is the BIBO stability, and besides, some other common prior
knowledge will be introduced later in Section 3.5.

###  3.2 Non-Causal Simulation-Induced Kernel

Following the guidelines mentioned above, we propose to design kernels for
non-causal impulse responses from a system theory perspective. As will be seen
shortly, although the same idea has been used in [10] to design kernels for
causal impulse responses, some challenges arise due to the non-causal nature
of the problem.

Before proceeding to the details, it is worth sketching the idea of our kernel
design method first. Following the first guideline for kernel design, i.e., to
mimic the optimal kernel (9), it is natural to embed the prior knowledge of
$G(q)$ from a system theoretic perspective, e.g., properness, stability,
dominant dynamics, into a stochastic process $g(t),t\in\mathbb{Z}$, and then
design the kernel as follows

$\displaystyle k(t,s)=\mathbb{C}ov(g(t),g(s)),\ t,s\in\mathbb{Z},$ (11)

where $\mathbb{C}ov(\cdot,\cdot)$ denotes the covariance between two random
variables. Following this idea, the question now becomes “how to embed the
prior knowledge of $G(q)$ into a stochastic process $g(t)$”?

One way is to employ the multiplicative uncertainty framework in robust
control, e.g., [16], as follows:

$\displaystyle G(q)=G_{0}(q)(1+G_{\Delta}(q))$ (12)

where $G_{0}(q)$ and $G_{\Delta}(q)$ are called the nominal model and
uncertainty, respectively. For the nominal model $G_{0}(q)$, we often
postulate a classical parameterization $G_{0}(q,\theta)$ in terms of $\theta$
to embed the prior knowledge of $G(q)$. For example, if the prior knowledge is
that $G(q)$ is proper, BIBO stable and its dominant dynamics is over-damped,
then we can simply choose $G_{0}(q,\theta)$ to be a first order system and
constrain the values of $\theta$ such that the pole of $G_{0}(q,\theta)$ is
within the unit circle. For the uncertainty $G_{\Delta}(q)$, we have to
suppose as little prior knowledge as possible, and for stable system
identification, the only prior knowledge imposed on $G_{\Delta}(q)$ should be
the stability, and in this case, we should design a zero mean Gaussian process
to embed such prior knowledge and model the impulse response of
$G_{\Delta}(q)$, then the task of kernel design would be done. To be more
specific, we draw the block diagram of (12) in Fig. 1, where $\delta(t)$ is a
unit impulsive input signal, $\bar{u}(t),g(t)$ are the input and output of
$G_{0}(q)$, respectively. Then clearly, if it is possible to construct a zero
mean Gaussian process to model the impulse response $h(t)$ of $G_{\Delta}(q)$
in Fig. 1, then $\bar{u}(t)$ and in turn $g(t)$ would be both Gaussian
processes. In this way, we finish the design of the kernel (11), and call (11)
the non-causal simulation-induced (NCSI) kernel.

$G_{\Delta}(q)$$+$$G_{0}(q)$$\delta(t)$$h(t)$$\bar{u}(t)$$g(t)$

Figure 1: The block diagram for the multiplicative uncertainty.

####  3.2.1 The Nominal Model

We will show how to postulate a classical parameterization $G_{0}(q,\theta)$
of $G_{0}(q)$ in terms of $\theta$ to embed the prior knowledge of $G(q)$,
e.g., properness, stability, dominant dynamics. In this section, we will only
consider the most fundamental prior knowledge that $G(q)$ is proper and BIBO
stable, and later in Section 3.5, we will consider some more specific prior
knowledge.

As well known from e.g., [21, 22], proper and BIBO stable $G_{0}(q)$ has the
following state-space model realization:

$\displaystyle x_{c}(t+1)=A_{c}x_{c}(t)+B_{c}\bar{u}(t),$ (13a) $\displaystyle
G_{0}(q,\theta):\quad$ $\displaystyle
x_{a}(t+1)=A_{a}x_{a}(t)+B_{a}\bar{u}(t),$ (13b) $\displaystyle
g(t)=C_{c}x_{c}(t)+C_{a}x_{a}(t)+D\bar{u}(t),$ (13c)

where $A_{c},B_{c},C_{c},A_{a},B_{a},C_{a},D$ have compatible dimensions,
$t\in\mathbb{Z}$ is the time index, $x_{c}(t)$ and $x_{a}(t)$ are the states
for the causal part and anti-causal part, respectively, such that all stable
poles are contained in $A_{c}$ and all unstable poles are contained in
$A_{a}$. Let $\theta=\\{A_{c},B_{c},C_{c},A_{a},B_{a},C_{a},D\\}$. Then the
state-space model (13) is a classical parameterization $G_{0}(q,\theta)$ of
$G_{0}(q)$.

####  3.2.2 The Uncertainty

Noting that for stable system identification, the only prior knowledge imposed
on $G_{\Delta}(q)$ should be its stability, we will show how to design a zero
mean Gaussian process (or equivalently a kernel) to embed such prior knowledge
and model the impulse response $h(t)$ of $G_{\Delta}(q)$. To this goal, since
there is no prior knowledge about the smoothness of $h(t)$, it is natural to
model $h(t)$ by the following Gaussian process:

$\displaystyle G_{\Delta}(q):\quad h(t)=b(t)w(t),\quad t\in\mathbb{Z},$ (14)

where $w(t)$ is a white Gaussian noise with zero mean and unit variance, and
$\\{b(t)\\}_{t\in\mathbb{Z}}\in\ell_{1}(\mathbb{Z})$. Clearly, $h(t)$ is a
zero mean Gaussian process with a non-causal kernel (covariance function)
$b(t)b(s)\delta_{t,s}$ with $t,s\in\mathbb{Z}$, where $\delta_{t,s}$ is the
Kronecker delta function, and which is an extension of the causal kernel [10,
eq. (21)].

####  3.2.3 The NCSI Kernel

Now we put (13) and (14) together and embed them in (12) to get the NCSI
kernel (11).

######  Lemma 3.1.

Consider the following state-space model

$\displaystyle x_{c}(t+1)=A_{c}x_{c}(t)+B_{c}b(t)w(t),$ (15a) $\displaystyle
x_{a}(t+1)=A_{a}x_{a}(t)+B_{a}b(t)w(t),$ (15b) $\displaystyle
g(t)=C_{c}x_{c}(t)+C_{a}x_{a}(t)+Db(t)w(t),$ (15c) $\displaystyle
t=-M_{c},-M_{c}+1,\cdots,\ M_{a},\text{for }M_{c},M_{a}\in\mathbb{N},$ (15d)

where $-M_{c}$ and $M_{a}$ are the starting time for the causal part and anti-
causal part, respectively. Assume that the initial states $x_{c}(-M_{c})$ for
the causal part and $x_{a}(M_{a})$ for the anti-causal part are such that
$[x_{c}(-M_{c})^{T},x_{a}(M_{a})^{T}]^{T}$ has zero mean and bounded
covariance matrix. Then for any finite number of sampling time instants
$t_{1},t_{2},\cdots,t_{m}\in\mathbb{Z}$,
$\lim_{M_{c}\to\infty}\lim_{M_{a}\to\infty}[g(t_{1}),\cdots,g(t_{m})]^{T}$
converges in distribution to a Gaussian random vector with mean zero and
covariance matrix, whose $(t,s)$th element is defined through

$\displaystyle k^{\text{NCSI}}(t,s)$
$\displaystyle=\sum\limits_{k=-\infty}^{\infty}b^{2}(k)g_{0}(t-k)g_{0}(s-k),$
(16a) $\displaystyle t,s\in\\{t_{1},t_{2},\cdots,t_{m}\\},$ where
$\\{g_{0}(t)\\}_{t\in\mathbb{Z}}$ denotes the non-causal impulse response of
the nominal model (13),i.e., $\displaystyle
g_{0}(t)=\left\\{\begin{array}[]{ll}C_{c}A_{c}^{t-1}B_{c},&\quad t\geq 1\\\
D-C_{a}A^{-1}_{a}B_{a},&\quad t=0\\\ -C_{a}A_{a}^{t-1}B_{a},&\quad
t\leq-1.\end{array}\right.,$ (16e)

######  Remark 3.1.

It is worth stressing that although the system theoretic perspective for the
design of causal kernels proposed in [10] is extended for that of non-causal
ones here, due to the difference in the nature of causal and non-causal
systems, the design of non-causal kernel has the combination of the following
issues need to be addressed:

  * •

the assumption on $M_{c}$ and $M_{a}$: shall they be treated as hyper-
parameters or chosen to be such that $M_{c}=+\infty$ and $M_{a}=+\infty$?

  * •

the assumption on $x_{c}(-M_{c})$ and $x_{a}(M_{a})$: shall they be assumed to
be a Gaussian random variable with zero mean and a parameterized covariance
matrix? and shall the mutual or partial independence be assumed on
$x_{c}(-M_{c})$, $x_{a}(M_{a})$ and $w(t)$?

  * •

the assumption on the uncertainty: the general guideline is that the
assumption on the uncertainty should be made as little as possible (because it
is uncertainty), but shall the uncertainty be assumed to be causal and stable,
or anti-causal and stable, or non-causal and stable?

Based on some analysis and moreover, Monte Carlo simulations, we have
solutions to the above issues, which are summarized in Lemma 4.1.

######  Remark 3.2.

Lemma 3.1 shows that as the starting time $M_{c}$ (actually $-M_{c}$) for the
causal part and $M_{a}$ for the anti-causal part go to infinity, then the
initial states $x_{c}(-M_{c})$ for the causal part and $x_{a}(M_{a})$ for the
anti-causal part will asymptotically have no influence on the NCSI kernel
(16). This is different from the causal case, where the initial state at
starting time $t=0$ has influence on the SI kernel, e.g., [10, Eq. (22)-(23)].
On the other hand, (16) shows that the NCSI kernel can be interpreted as an
infinite sum of rank-1 kernels, $g_{0}(t-k)g_{0}(s-k)$, weighted by
$b^{2}(k)$, see Fig. 2.

Figure 2: Interpretation of the NCSI kernel (16) as the infinite sum of rank-1
kernels, $g_{0}(t-k)g_{0}(s-k)$, weighted by $b^{2}(k)$, where we choose
$g_{0}(t)=10\exp(-|t|)$, and $b(k)=\exp(-|k|)$.

###  3.3 Stability of NCSI Kernel

As can be seen in the following theorem, the NCSI kernel (16) is stable in the
sense that its associated RKHS $\mathcal{H}_{k}\subset\ell_{1}(\mathbb{Z})$,
see [23, 6], where $\ell_{1}(\mathbb{Z})$ is set of absolutely summable
sequences of real numbers, indexed by $\mathbb{Z}$.

######  Theorem 3.1.

If $\\{g_{0}(t)\\}_{t\in\mathbb{Z}}\in\ell_{1}(\mathbb{Z})$ and
$\\{b(t)\\}_{t\in\mathbb{Z}}\in\ell_{1}(\mathbb{Z})$, then the NCSI kernel
(16) is stable.

######  Remark 3.3.

The technique used in the proof of Theorem 3.1 can be used to prove the
stability of the causal SI kernel proposed in [10], resulting that the
sufficient condition [10, eq. (29a)] can be replaced by a weaker one
$\\{b(t)\\}_{t\in\mathbb{N}}\in\ell_{1}(\mathbb{N})$.

###  3.4 Semiseparability of NCSI Kernel

It was shown in [20] that, if the kernel for regularized impulse response
estimation has some special rank structures, e.g., semiseparability, [24], the
computational complexity of KRM can be reduced. For example, the widely used
DC and TC kernels (4) are both semiseparable, and the simulation-induced
kernel for causal impulse response estimation is also semiseparable [20].
Interestingly, we will show below that, the NCSI kernel (16) also has some
special rank structure and in fact, it is semiseparable plus diagonal.

######  Definition 3.1.

A positive semidefinite kernel
$k(t,s):\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ is said to be
extended $p$-semiseparable plus diagonal, if there exist
$\mu_{i},\nu_{i}:\mathbb{Z}\rightarrow\mathbb{R}$, $i=1,\cdots,p$, such that

$\displaystyle
k(t,s)=\left\\{\begin{array}[]{lll}\sum_{i=1}^{p}\mu_{i}(t)\nu_{i}(s),&\quad
t>s\\\ d(t),&\quad t=s\\\ \sum_{i=1}^{p}\nu_{i}(t)\mu_{i}(s),&\quad
t<s\end{array}\right.$ (20)

where $\mathcal{X}$ represents either $\mathbb{N}$ or $\mathbb{Z}$,
$p\in\mathbb{N}$, and $d(t)\in\mathbb{R}$ with $t\in\mathcal{X}$. Moreover, if
$d(t)$ further satisfies

$\displaystyle d(t)=\sum_{i=1}^{p}\mu_{i}(t)\nu_{i}(t),$ (21)

then $k(t,s)$ is said to be extended $p$-semiseparable.

Then we show that the NCSI kernel (16) is extended $p$-semiseparable plus
diagonal.

######  Theorem 3.2.

Consider the NCSI kernel (16) and assume that $A_{c}$ is non-singular. Then
the NCSI kernel (16) is extended $\bar{p}$-semiseparable plus diagonal with
$\bar{p}\in\mathbb{N}$ and $\bar{p}\leq p$, where $p$ is the dimension of
$A_{c}$ in (15a) plus the dimension of $A_{a}$ in (15b).

######  Corollary 3.1.

Consider the NCSI kernel (16) and assume that $A_{c}$ is non-singular. Then
the NCSI kernel (16) is extended $\bar{p}$-semiseparable, if

$\displaystyle D=0,\quad\text{ or }\quad
D-C_{a}A_{a}^{-1}B_{a}=C_{c}A_{c}^{-1}B_{c}$

######  Remark 3.4.

One may wonder what would happen for the case when $A_{c}$ is singular. In
this case, the NCSI kernel (16) is not semiseparable any more, but its
corresponding kernel matrix may still have some rank-structure properties.
Since this case is less interesting, the corresponding details will not be
given here.

As shown in [25], there are many efficient algorithms on matrix computations
for extended $p$-semiseparable plus diagonal kernels. We will illustrate this
with a concrete example in the next section.

###  3.5 Embedding More Prior Knowledge

In this section, we try to embed some more specific prior knowledge of $G(q)$
into the NCSI kernel (16).

####  3.5.1 Over-Damped Dominant Dynamics

First, we assume that the prior knowledge is that the dominant dynamics of
$G(q)$ is over-damped and the function $b(t)$ in the model of the uncertainty
$G_{\Delta}(q)$ in (14) is an exponential decay one. Therefore, the state-
space model (13) of $G_{0}(q)$ is parameterized as follows:

$\displaystyle A_{c}=B_{c}=a_{c},\quad A_{a}^{-1}=B_{a}=a_{a},$ $\displaystyle
C_{c}=c_{c},\quad-C_{a}A_{a}^{-2}=c_{a},\quad D-C_{a}A^{-1}_{a}B_{a}=c_{0},$
$\displaystyle-1<a_{c},a_{a}<1,\ -\infty<c_{c},c_{0},c_{a}<\infty,$ such that
$\displaystyle g_{0}(t)=\left\\{\begin{array}[]{ll}c_{c}a_{c}^{t},&\quad
t>0\\\ c_{0},&\quad t=0\\\ c_{a}a_{a}^{-t},&\quad t<0\end{array}\right.,$
(22d) and moreover, $b(t)$ is parameterized as follows: $\displaystyle
b(t)=\left\\{\begin{array}[]{ll}\sigma_{c}\lambda_{c}^{t/2},&\quad t>0\\\
\sigma_{0},&\quad t=0\\\ \sigma_{a}\lambda_{a}^{-t/2},&\quad
t<0\end{array}\right.,$ (22h) $\displaystyle 0<\lambda_{c},\lambda_{a}<1,\
\sigma_{c},\sigma_{a}\geq 0.$ As a result, we obtain the following non-causal
kernel $\displaystyle k^{\text{NCSI-
FO}}(t,s;\eta)=\text{\eqref{eq:si_expre}}\text{ with
}\eqref{eq:ncsi_fo1}\text{ and }\eqref{eq:ncsi_fo2},$ (22i)
$\displaystyle\eta=[a_{c},a_{a},c_{c},c_{0},c_{a},\lambda_{c},\lambda_{a},\sigma_{c},\sigma_{a},\sigma_{0}].$

Since the nominal model is a first order system for both the causal and anti-
causal part, the kernel (22) is called the first order NCSI (NCSI-FO) kernel.
The NCSI-FO kernel (22) has a closed-form expression that can be found in
Appendix A.4. As shown in Remark A.1, one of the scaling hyper-parameters
$c_{c},c_{0},c_{a},\sigma_{c},\sigma_{0},\sigma_{a}$ is redundant and thus,
$\sigma_{0}$ is fixed to be 1 and not treated as a hyper-parameter hereafter.

######  Example 3.1.

It follows from Theorem 3.2 that, the NCSI-FO kernel (22) is an extended
$2$-semiseparable plus diagonal kernel. Define the kernel matrix
$K^{\text{NCSI-FO}}_{ij}=k^{\text{NCSI-FO}}(i,j;\eta_{0})$,
$i,j=-n/2,\cdots,n/2$, with $n=1000,\cdots,8000$. Then Fig. 3 shows the time
for computing the Cholesky factor of $K^{\text{NCSI-FO}}$ by using chol in
MATLAB (in blue) and by exploiting its high-order semiseparable structure with
Algorithm 3 in [25] (in red).

Figure 3: The time cost for computing the Cholesky factor of $K^{\text{NCSI-
FO}}$ by using chol in MATLAB (in blue) and by exploiting its high-order
semiseparable structure with Algorithm 3 in [25] (in red).

####  3.5.2 Mirrored Poles

Second, we further assume that there is an extra prior knowledge that the
dominant dynamics of $G(q)$ has a pair of poles mirrored with respect to the
unit circle. As mentioned in [15, Sec. 4.3], mechanical structures may have
such prior knowledge. In this case, it is natural to further let $a_{c}=a_{a}$
in the NCSI-FO kernel (22), leading to the following kernel

$\displaystyle k^{\text{NCSI-
FO}^{\text{mp}}}(t,s;\eta)=\eqref{eq:kernel_si_9}\text{ with }a_{c}=a_{a},$
(23)
$\displaystyle\eta=[a_{c},c_{c},c_{0},c_{a},\lambda_{c},\lambda_{a},\sigma_{c},\sigma_{a},\sigma_{0}].$

###  3.6 Kernel Simplification

Since the NCSI-FO kernel (22) has 9 hyper-parameters, its hyper-parameter
estimation may have some difficulties, e.g., its initialization is harder, its
local minima issues may be more severe, and it may take more time to find the
optimal solution. Hence, following the second guideline for kernel design,
i.e., to ease the computation of KRM, we simplify the NCSI-FO kernel (22)
below by linking some hyper-parameters together.

There are many ways to link the hyper-parameters. First, we can link the decay
rate of the nominal model with that of the uncertainty as follows:

$\displaystyle k^{\text{NCSI-DC}}(t,s;\eta)=k^{\text{NCSI-FO}}(t,s;\eta)\text{
with: }$ $\displaystyle\qquad a_{c}=\rho_{c}\lambda_{c}^{1/2},\
a_{a}=\rho_{a}\lambda_{a}^{1/2},\ -1\leq\rho_{c},\rho_{a}\leq 1,$
$\displaystyle\qquad\sigma_{c}=\sqrt{1-\rho_{c}^{2}},\
\sigma_{a}=\sqrt{1-\rho_{a}^{2}},$ (24)
$\displaystyle\eta=[\lambda_{c},\lambda_{a},c_{c},c_{0},c_{a},\rho_{c},\rho_{a}].$

The kernel (24) is called the NCSI-DC kernel, because it has a nice property
similar to the DC kernel (4a), that is, the kernel (24) reduces to a rank-1
kernel, i.e., $k(t,s)=g_{0}(t)g_{0}(s)$, for $\rho_{c}=\rho_{a}=1$, and a
diagonal kernel, i.e., $k(t,s)=b(t)b(s)\delta_{t,s}$, for
$\rho_{c}=\rho_{a}=0$, which shows different degrees of balance between the
nominal model $G_{0}(q)$ (described by $g_{0}(t)$) and the uncertainty
$G_{\Delta}(q)$ (described by $b(t)$) by varying $\rho_{c}$ and $\rho_{a}$.
Moreover, similar to the relation between the TC kernel (4b) and DC kernel
(4a), we can obtain another special case of the NSCI-DC kernel (24) as follows

$\displaystyle k^{\text{NCSI-TC}}(t,s;\eta)=k^{\text{NCSI-DC}}(t,s;\eta)\text{
with: }$ $\displaystyle\qquad\rho_{c}=\lambda_{c}^{1/2},\
\rho_{a}=\lambda_{a}^{1/2},$ (25)
$\displaystyle\eta=[\lambda_{c},\lambda_{a},c_{c},c_{0},c_{a}].$

which is called the NCSI-TC kernel.

######  Remark 3.5.

The $\text{NCSI-FO}^{\text{mp}}$ kernel (23) can also be seen as being
obtained from the NCSI-FO kernel (22) by hyper-parameter simplification, i.e.,
by letting $a_{c}=a_{a}$. The difference between the $\text{NCSI-
FO}^{\text{mp}}$ kernel (23) and the NSCI-DC kernel (24) (also the NSCI-TC
kernel (25)) is that the extra prior knowledge embedded in the former due to
the hyper-parameter simplification $a_{c}=a_{a}$ has a clear physical
interpretation, i.e., a pair of poles mirrored with respect to the unit
circle, but not in the latter.

######  Remark 3.6.

Despite having the same suffix ’TC’, as illustrated in Fig. 4, the NCSI-TC
kernel (25) has very different behavior from the NC-TC kernel (5) and its
block diagonal variant in [15],

$\displaystyle k^{\text{NCBD-
TC}}(t,s;\eta)=\left\\{\begin{array}[]{ll}c\min\\{\lambda_{c}^{t},\lambda_{c}^{s}\\},&\text{if
}t\geq 0,\ s\geq 0\\\ c\min\\{\lambda_{a}^{-t},\lambda_{a}^{-s}\\},&\text{if
}t<0,\ s<0\\\ 0,&\text{otherwise}.\end{array}\right.,$ (29)
$\displaystyle\eta=[c,\lambda_{c},\lambda_{a}],$

which is called the NCBD-TC kernel below. Similar to the construction of the
NCBD-TC kernel, we can also construct the NCBD-DC kernel below

$\displaystyle k^{\text{NCBD-
DC}}(t,s;\eta)=\left\\{\begin{array}[]{ll}c_{c}^{2}\lambda_{c}^{(t+s)/2}\rho_{c}^{|t-s|},&\text{if
}t\geq 1,s\geq 1\\\ c_{0}^{2},&\text{if }t=0,s=0\\\
c_{a}^{2}\lambda_{a}^{-(t+s)/2}\rho_{a}^{|t-s|},&\text{if }t\leq-1,s\leq-1\\\
0,&\text{otherwise}\end{array}\right.,$ (34)
$\displaystyle\eta=[\lambda_{c},\lambda_{a},c_{c},c_{0},c_{a},\rho_{c},\rho_{a}].$

Figure 4: The illustrations of three TC-like non-causal kernels with
$\lambda_{c}=0.9$, $\lambda_{a}=0.8$, $c=1$ for the NC-TC kernel (5) and its
block diagonal variant, the NCBD-TC kernel (29), and $\lambda_{c}=0.9$,
$\lambda_{a}=0.8$, $c_{c}=-1$, $c_{0}=c_{a}=1$ for the NCSI-TC kernel (25).
One significant difference of the NCSI-TC kernel (25) is its flexibility of
choosing negative off-diagonal blocks, which partially accounts for its
advantages over the NC-TC kernel (5) and the NCBD-TC kernel (29).

##  4 Numerical Simulation

In this section, we run simulations to test the kernels proposed in Sections
3.5 and 3.6, and the ones in [15].

###  4.1 Test Systems and Data-Bank

As mentioned in the last of Section 2.1, system (1) often arises from the
inverse model control, e.g., [15]. More specifically, if the control plant in
the feedforward control is denoted by $P(q)$ and assumed to be unknown,
stable, non-minimum phase and has no zeros on the unit circle, then its non-
causal inverse $G(q)$ can be put in the form of (1), e.g., [15, Theorem 1].
Moreover, it should be noted that the input $u(t)$ and output $y(t)$ of $G(q)$
correspond to the output and input of $P(q)$, respectively.

In what follows, we will generate four data banks D1-D4, associated with
different types of test systems. We first introduce four types of test systems
$P(q)$, then generate the input-output data for each test system $P(q)$, and
finally the input-output data for the non-causal inverse of $P(q)$, i.e.,
$G(q)$ in (1).

First, we generate the test systems in D1-D4.

  * •

Data-bank D1 contains a single non-minimum phase system studied in [15]

$\begin{array}[]{cc}P(q)=\frac{1.550(q^{2}-2.035q+1.052)(q^{2}-1.844q+0.9391)}{q^{2}(q-0.9514)(q-0.9511)},\end{array}$
(35)

with two stable zeros at $0.922\pm 0.298i$ and two unstable zeros at $1.018\pm
0.126i$.

  * •

Data-bank D2 contains 1000 randomly generated systems $P(q)$. We generate each
system as follows. First, a 30th order continuous-time system is generated
using function m=rss(30) in MATLAB. Then system m is sampled at 100 times of
its bandwidth by function md=c2d(m,2pi/(100$\cdot$bandwidth(m))) in MATLAB.
Finally, if the system md satisfies the following two conditions:

    1. (a)

md is a stable, bi-proper and non-minimum phase system,

    2. (b)

md has all poles and zeros within $\\{z\in\mathbb{C}:|z|<0.96\text{ or
}|z|>1.04\\}$,

then md is saved as one of 1000 systems in D2.

  * •

Data-bank D3 contains 1000 randomly generated systems $P(q)$. We generate each
system as follows. For each test system md in D2, we replace the two zeros
(can be two real zeros or a pair of complex conjugate zeros) having the
smallest magnitude with two new zeros sampled from the uniform distributions,
$\mathcal{U}[0.9,0.8]$ and $\mathcal{U}[1.1,1.2]$, respectively. This new
system is then saved as one of 1000 test systems in D3 and has the feature
that its non-causal inverse $G(q)$ has dominant over-damped dynamics.

Data-bank D4 contains 1000 randomly generated systems. Each test system takes
the form of

$\displaystyle
P(q)=\frac{(q-z_{1})(q-1/z_{1})(q-z_{2})(q-1/z_{2})}{(q-p_{1})(q-p_{2})(q-p_{3})(q-p_{4})},$

where $z_{1}=0.9$, and $z_{2},p_{1},p_{2},p_{3},p_{4}$ are all sampled from
the uniform distribution $\mathcal{U}[0,0.9]$, and thus its non-causal inverse
$G(q)$ has mirrored poles.

For illustration, some true non-causal impulse responses $g^{0}(t)$,
$t\in\mathbb{Z}$, of the non-causal inverse $G(q)$ of $P(q)$, are shown in
Fig. 5. The basic prior knowledge that the test systems contain in each data
banks of D2-D4 are summarized below follows:

  * D2:

BIBO stability and properness;

  * D3:

BIBO stability, properness, and over-damped dominant dynamics;

  * D4:

BIBO stability, properness, and over-damped dominant dynamics with mirrored
poles.

Second, we generate the input-output data in D1-D4. Each of D1-D4 contains
1000 data sets $\\{u(t),y(t)\\}^{N_{max}}_{t=1}$, where $N_{max}$ is 2000 for
D1 and 700 for D2-D4, respectively. Moreover, we generate 1000 data sets for
the single test system (35) in D1, but one data set for each test system in
D2-D4. More specifically, the input signal $u(t)$ in (1) is simulated by
MATLAB function: u(t)=lsim(md,$y_{0}(t)$) from $t=-n_{c}$ to
$t=N_{max}+n_{a}$, where $y_{0}(t)$ is a white Gaussian noise with unit
variance. Then the output signal $y(t)$ in (1) is obtained by perturbing
$y_{0}(t)$ with an additive white Gaussian noise $v(t)$ with mean zero and the
variance one tenth of the variance of $y_{0}(t)$. Finally,
$\\{y(t),u(t)\\}_{t=1}^{N_{max}}$ is saved as one test input-output data for
the identification of (1).

Figure 5: The illustrations of the true non-causal impulse response $g^{0}$:
the single system (35) in D1, and five system examples in D2-D4.

###  4.2 Simulation Setup

We first truncate the non-causal impulse response (2) at sufficiently large
orders $n_{c}$ and $n_{a}$, and obtain the so-called non-causal finite impulse
response (FIR) model

$\displaystyle y(t)=\sum\limits_{\tau=-n_{a}}^{n_{c}}g(\tau)u(t-\tau)+v(t),$
(36)

which can be rewritten into the matrix-vector form as

$\displaystyle Y=\Psi\theta+V,$ (37)

where

$\displaystyle Y=[y(n_{c}+1)\ y(n_{c}+2)\ \cdots\ y(N-n_{a})]^{\text{T}},$
$\displaystyle V=[v(n_{c}+1)\ v(n_{c}+2)\ \cdots\ v(N-n_{a})]^{\text{T}},$
$\displaystyle\theta=[g(-n_{a})\ g(-n_{a}+1)\ \cdots\ g(n_{c})]^{T},$
$\displaystyle\Psi=\left[\begin{array}[]{cccc}u(1+n_{c}+n_{a})&u(n_{c}+n_{a})&\cdots&u(1)\\\
u(2+n_{c}+n_{a})&u(1+n_{c}+n_{a})&\cdots&u(2)\\\
\vdots&\vdots&\ddots&\vdots\\\ u(N)&u(N-1)&\cdots&u(N-n_{c}-n_{a})\\\
\end{array}\right].$

Accordingly, the regularized least square (RLS) estimate (3) can be rewritten
as

$\displaystyle\hat{\theta}^{\text{R}}=\underset{\theta\in\mathbb{R}^{n}}{\text{argmin}}||Y-\Psi\theta||_{2}^{2}+\sigma^{2}\theta^{T}K(\eta)^{-1}\theta,$
(38)

where $\hat{\theta}^{\text{R}}=[\hat{g}^{\text{R}}(-n_{a})\
\hat{g}^{\text{R}}(-n_{a}+1)\ \cdots\ \hat{g}^{\text{R}}(n_{c})]^{T}$,
$\sigma^{2}$ is the noise variance of $v(t)$, and
$K(\eta)\in\mathbb{R}^{n\times n}$ is a positive semidefinite matrix with its
$(i,j)$th entry $K_{i,j}(\eta)$ corresponding to $k(i,j;\eta)$ in (3).

The non-causal FIR model orders, $n_{c}$ and $n_{a}$, in (36) are both chosen
to be 150 for D1 and 50 for D2-D4, respectively. For each data bank, we will
use data sets with the different sample size $N$ with $N\leq N_{max}$, and
test the following kernels:

  * •

the NC-TC kernel (5),

  * •

the NCBD-TC kernel (29),

  * •

the NCBD-DC kernel (34),

  * •

the NCBD-TC${}^{\text{mp}}$ kernel, (29) with $\lambda_{c}=\lambda_{a}$,

  * •

the NCBD-DC${}^{\text{mp}}$ kernel, (34) with $\lambda_{c}=\lambda_{a}$,

  * •

the NCSI-TC kernel (25),

  * •

the NCSI-DC kernel (24),

  * •

the NCSI-FO kernel (22),

  * •

the NCSI-FO${}^{\text{mp}}$ kernel (23).

######  Remark 4.1.

The NCBD-DC kernel (34) is proposed by mimicking the structure of the NCBD-TC
kernel (29). It is worth noting that the NCBD-DC kernel (34) and the NCSI-DC
kernel (24) have the same number of hyper-parameters.

The hyper-parameters are estimated using the empirical Bayes (EB) method,
i.e.,

$\displaystyle\hat{\eta}=\underset{\eta\in\Omega_{k}}{\text{argmin}}\\{Y(\Psi
K(\eta)\Psi^{T}+\sigma^{2}I_{N-n})^{-1}Y+$ $\displaystyle\qquad\log\det(\Psi
K(\eta)\Psi^{T}+\sigma^{2}I_{N-n})\\},$ (39)

where $\sigma^{2}$ is estimated from the sample variance of the estimated non-
causal FIR model (36) using the LS method. The optimization (39) is
implemented by the MATLAB function MultiStart with fmincon and the number of
runs being ten times the number of hyper-parameters. Substituting (39) and the
estimate of $\sigma^{2}$ into (38), we obtain the estimate
$\hat{\theta}^{\text{R}}$.

The performance of the estimate $\hat{\theta}^{\text{R}}$ is evaluated by the
model fit [2] and the tracking error [15]. More specifically, the model fit
(FIT) is defined by

$\displaystyle\text{FIT}=100\left(1-\left[\frac{\sum^{n_{c}}_{k=-n_{a}}|g^{0}(k)-\hat{g}^{\text{R}}(k)|^{2}}{\sum^{n_{c}}_{k=-n_{a}}|g^{0}(k)-\bar{g}^{0}|^{2}}\right]^{\frac{1}{2}}\right),$
$\displaystyle\bar{g}^{0}=\frac{1}{n}\sum^{n_{c}}_{k=-n_{a}}g^{0}(k),$

where $\\{g^{0}(k)\\}_{k\in\mathbb{Z}}$ is the true non-causal impulse
response of $P^{-1}(q)$, and the tracking error (ERR) is defined by

$\displaystyle\text{ERR}=\sqrt{\sum_{k=1}^{L}e(k)^{2}/L},$ $\displaystyle
e(t)=r(t)-P(q)\sum\limits_{k=-n_{a}}^{n_{c}}\hat{g}^{\text{R}}(k)q^{-k}r(t),$
$\displaystyle t=1-2n_{c},\cdots,L+n_{a},$

where the test length $L=1000$, and the test reference signal $r(t)$ is a
white Gaussian noise with unit variance.

Figure 6: Box plots of FITs and ERRs for the data sets in D1-D3.  Figure 7:
Estimates $\hat{\theta}^{\text{R}}$ for one data set in D1 with $N=600$, and
the true non-causal impulse response $\tilde{g}_{0}$ (dot line). The estimates
$\hat{\theta}^{\text{R}}$ with the NCSI kernels (solid lines) show a better
model fit than the existing non-causal kernels (dash lines).

Figure 8: Box plots of FITs and ERRs for the data sets in D4.  Table 1:
Average FITs and ERRs, and standard deviation of FITs and ERRs for the data
sets in D1-D3.

avg. FIT (standard deviation)
---
Data-bank | N | NC-TC | NCBD-TC | NCBD-DC | NCSI-TC | NCSI-DC | NCSI-FO
D1 | 600 | 12.98( 7.62) | 27.06(10.47) | 30.51(11.98) | 47.64( 8.98) |  52.97( 7.77) | 52.69( 8.03)
D1 | 2000 | 49.48(10.11) | 60.41( 4.06) | 62.78( 3.80) | 76.02( 3.48) |  76.96( 3.25) | 76.85( 3.33)
D2 | 200 | 76.22(16.61) | 77.49(16.50) | 78.62(16.09) | 78.78(15.15) | 81.09(13.07) |  81.18(13.14)
D2 | 700 | 88.54(10.87) | 88.98(11.09) | 89.46(11.06) | 89.87(10.38) | 91.16( 9.65) |  91.28( 9.61)
D3 | 200 | 47.58(23.40) | 55.39(23.65) | 59.58(25.26) | 63.37(20.25) | 65.16(21.29) |  67.20(21.45)
D3 | 700 | 74.76(16.61) | 79.26(16.40) | 81.73(16.54) | 83.13(14.93) | 84.45(15.17) |  85.41(15.25)
avg. ERR (standard deviation)
Data-bank | N | NC-TC | NCBD-TC | NCBD-DC | NCSI-TC | NCSI-DC | NCSI-FO
D1 | 600 | 0.329(0.030) | 0.291(0.036) | 0.287(0.042) | 0.205(0.025) | 0.189(0.031) |  0.188(0.029)
D1 | 2000 | 0.165(0.020) | 0.137(0.012) | 0.134(0.011) | 0.093(0.008) | 0.087(0.007) |  0.087(0.007)
D2 | 200 | 0.194(0.069) | 0.181(0.067) | 0.171(0.067) | 0.175(0.062) |  0.163(0.058) | 0.165(0.057)
D2 | 700 | 0.083(0.027) | 0.078(0.027) | 0.073(0.027) | 0.074(0.024) | 0.068(0.022) |  0.067(0.022)
D3 | 200 | 0.280(0.073) | 0.240(0.062) | 0.219(0.072) | 0.211(0.053) | 0.198(0.060) |  0.189(0.063)
D3 | 700 | 0.124(0.025) | 0.102(0.024) | 0.091(0.025) | 0.089(0.019) | 0.081(0.020) |  0.077(0.021)

Table 2: Average FITs and ERRs and standard deviation of FITs and ERRs for the
data sets in D4.

avg. FIT (standard deviation)
---
Data-bank | N | NCBD-TC | NCBD-TC${}^{\text{mp}}$ | NCBD-DC | NCBD-DC${}^{\text{mp}}$ | NCSI-FO | NCSI-FO${}^{\text{mp}}$
D4 | 200 | 34.26(37.70) | 32.00(39.83) | 49.17(37.61) | 51.86(33.74) | 64.33(30.08) |  69.31(27.35)
D4 | 700 | 61.80(27.65) | 64.99(27.23) | 78.15(20.33) | 77.50(20.65) | 83.84(16.88) |  85.32(17.13)
avg. ERR (standard deviation)
Data-bank | N | NCBD-TC | NCBD-TC${}^{\text{mp}}$ | NCBD-DC | NCBD-DC${}^{\text{mp}}$ | NCSI-FO | NCSI-FO${}^{\text{mp}}$
D4 | 200 | 0.208(0.065) | 0.236(0.067) | 0.178(0.063) | 0.170(0.053) | 0.140(0.053) |  0.132(0.049)
D4 | 700 | 0.102(0.034) | 0.106(0.032) | 0.074(0.026) | 0.0761(0.026) | 0.060(0.023) |  0.058(0.023)

###  4.3 Results, Observations and Take-home Messages

First, we consider the data-banks D1-D3. The box plots of FITs and ERRs are
shown in Fig. 6, and the average FITs and ERRs and their standard deviations
are shown in Table 1. For illustration, the estimates $\hat{\theta}^{R}$ for
one data set in D1 with $N=600$ are shown in Fig. 7. Based on these results,
we have the following observations:

  * •

The kernels more to the left in general give better (and thus the NCSI-FO
kernel (22) in general gives the best) performance in terms of the average FIT
and ERR, and robustness of the FITs and ERRs through the distribution and also
the standard deviation of the FITs and ERRs. In particular, for data-banks
D1-D2, the NCSI-DC kernel (24) and NCSI-FO kernel (22) have very close
performance

Second, we consider the data-bank D4. The box plots of FITs and ERRs are shown
in Fig. 8, and the average FITs and ERRs and their standard deviations are
shown in Table 2. Then we have the following observations:

  * •

in contrast with the NCBD-TC kernel (29) (resp. the NCBD-DC kernel (34)), the
NCBD-TC${}^{\text{mp}}$ kernel (resp. the NCBD-DC${}^{\text{mp}}$ kernel) does
not give consistent improvement in the performance in terms of the average FIT
and ERR, and robustness of the FITs and ERRs through the distribution and also
the standard deviation of the FITs and ERRs; on the other hand, in contrast
with the NCSI-FO kernel (22), the NCSI-FO${}^{\text{mp}}$ kernel (23) gives
consistent improvement in the performance in terms of the average FIT and
robustness of the FITs and ERRs through the distribution of the FITs and ERRs,
but a bit worse standard deviation.

Finally, we have the following take-home messages:

  * •

in general, the NCSI-FO kernel (22) is recommended to be used, and if
computation issues are concerned, then the NCSI-DC kernel (24) can be tried as
an alternative;

  * •

if there is extra prior knowledge (e.g., mirrored poles), then it is suggested
to further explore the structure of the NCSI-FO kernel (22) and design a more
specific kernel (e.g., the NCSI-FO${}^{\text{mp}}$ kernel) to embed such prior
knowledge

##  5 Conclusion

In this paper, we studied the non-causal system identification problem by
using the kernel-based regularization method. To tackle the key difficulty of
non-causal kernel design, we first introduced the guidelines, then extended
the system theoretic framework to design the so-called non-causal simulation-
induced kernel and studied its properties, including stability and
semiseparability. The numerical simulation results showed that the proposed
kernels can give better model estimates/tracking performance than the state-
of-art kernels. The importance of our obtained results lies in that we develop
a systematic framework to design kernels for non-causal systems, the prior
knowledge embedded in the designed kernels have clear physical interpretation,
and the designed kernels are semiseparable and thus can be used to develop
efficient implementations.

## Appendix

###  A.1 Proof of Lemma 3.1

We will prove Lemma 3.1 in two steps. First, it is straightforward to show
that, given $x_{c}(-M_{c})$, $x_{a}(M_{a})$, $g(t)$ takes the form of

$\displaystyle g(t)=\left\\{\begin{array}[]{ll}&g_{ini}(t)+g_{imp}(t),\text{
if }t,s=-M_{c},\cdots,M_{a}\\\ &0,\quad\text{otherwise}\end{array}\right.,$
(A.3) $\displaystyle
g_{ini}(t)=C_{c}A_{c}^{t+M_{c}}x_{c}(-M_{c})+C_{a}A_{a}^{t-M_{a}}x_{a}(M_{a}),$
$\displaystyle g_{imp}(t)=\sum\limits_{k=-M_{c}}^{M_{a}-1}g_{0}(t-k)b(k)w(k),$

where $g_{0}(t)$ is defined in (16e), and moreover,
$\\{g_{0}(t)\\}_{t\in\mathbb{Z}}\in\ell_{1}(\mathbb{Z})$.

Second, we let $Z=[g(t_{1}),\cdots,g(t_{m})]^{T}$ and then decompose $Z$ as
$Z=Z_{imp}+Z_{ini}$, where

$\displaystyle Z_{ini}=[g_{ini}(t_{1}),\cdots,g_{ini}(t_{m})]^{\text{T}},$
$\displaystyle Z_{imp}=[g_{imp}(t_{1}),\cdots,g_{imp}(t_{m})]^{\text{T}},$
$\displaystyle-M_{c}\leq t_{1},t_{2},\cdots,t_{m}\leq M_{a}.$

Then we show that as ${M_{c}\to\infty}$ and ${M_{a}\to\infty}$, $Z_{ini}$
converges in probability to a zero vector and $Z_{imp}$ converges in
distribution to a Gaussian random vector.

On the one hand, for $Z_{ini}$ and for any $\epsilon>0$, we have

$\displaystyle P$ $\displaystyle\bigg{(}\|Z_{ini}\|_{2}>m\epsilon\bigg{)}$
$\displaystyle\leq P\bigg{(}|g_{ini}(t_{1})|>\epsilon\text{ or }\cdots\text{
or }|g_{ini}(t_{m})|>\epsilon\bigg{)}$
$\displaystyle\leq\sum\limits_{i=1}^{m}P\bigg{(}|g_{ini}(t_{i})|>\epsilon\bigg{)}\leq\sum\limits_{i=1}^{m}\frac{\mathbb{C}ov(g_{ini}(t_{i}),g_{ini}(t_{i}))}{\epsilon^{2}}$

where $\|\cdot\|_{2}$ is the Euclidean norm and the last inequality is due to
the Chebyshev’s inequality. Under the assumption that
$x_{ini}(M_{c},M_{a})=[x_{c}(-M_{c})^{T},x_{a}(M_{a})^{T}]^{T}$ has zero mean
and bounded covariance matrix $\Sigma_{ini}(M_{c},M_{a})$, we have

$\displaystyle\mathbb{C}ov(g_{ini}(t_{i}),g_{ini}(t_{i}))=[C_{c}A_{c}^{t_{i}+M_{c}}\
C_{a}A_{a}^{t_{i}-M_{a}}]\cdots$ $\displaystyle\quad\Sigma_{ini}(M_{c},M_{a})\
[C_{c}A_{c}^{t_{i}+M_{c}}\ C_{a}A_{a}^{t_{i}-M_{a}}]^{T}\rightarrow 0$

as $M_{c},M_{a}\rightarrow\infty$, and moreover,

$\displaystyle\lim\limits_{M_{c}\rightarrow\infty}\lim\limits_{M_{a}\rightarrow\infty}P\bigg{(}\|Z_{ini}\|_{2}>m\epsilon\bigg{)}\leq
0,$ (A.4)

implying that $Z_{ini}$ converges in probability to a zero vector.

On the other hand, for $Z_{imp}$, we consider its characteristic function
given by

$\displaystyle\varphi(v)=\mathbb{E}[\exp(jv^{T}Z_{imp})]=\exp(-\frac{1}{2}v^{T}\Sigma\
v)\text{ with }$
$\displaystyle(\Sigma)_{i,j}=\mathbb{C}ov(g_{imp}(t_{i}),g_{imp}(t_{j})),\
1\leq i,j\leq m,$
$\displaystyle=\sum\limits_{k=-M_{c}}^{M_{a}-1}b^{2}(k)g_{0}(t_{i}-k)g_{0}(t_{j}-k).$

Then we have

$\displaystyle\lim\limits_{M_{c}\rightarrow\infty}\lim\limits_{M_{a}\rightarrow\infty}\varphi(v)=\tilde{\varphi}(v)=\exp(-\frac{1}{2}v^{T}\
\tilde{\Sigma}\ v),$ (A.5) $\displaystyle\text{ with
}(\tilde{\Sigma})_{i,j}=\sum\limits_{k=-\infty}^{\infty}b^{2}(k)g_{0}(t_{i}-k)g_{0}(t_{j}-k),$

where $(\tilde{\Sigma})_{i,j}$ exists for all $i,j=1,\cdots,m$ due to the
boundedness of $b^{2}(t)$ and
$\\{g_{0}(t)\\}_{t\in\mathbb{Z}}\in\ell_{1}(\mathbb{Z})$. It then follows from
the Levy’s continuity theorem, e.g., [26, p. 204], that (A.5) leads to

$\displaystyle
Z_{imp}\overset{d}{\longrightarrow}\mathcal{N}(0,\tilde{\Sigma}),\quad\text{
as }M_{c},M_{a}\rightarrow\infty,$ (A.6)

where $\overset{d}{\longrightarrow}$ denotes the convergence in distribution.

Finally, from the multivariate version of the Slutsky’s theorem, e.g., [27,
Thm. 3.4.3], (A.4) and (A.6), it holds that
$\lim_{M_{c}\to\infty}\lim_{M_{a}\to\infty}Z$ converges in distribution to a
Gaussian vector with mean zero and covariance matrix $\tilde{\Sigma}$. This
completes the proof.

###  A.2 Proof of Theorem 3.1

Note that

$\displaystyle\sum\limits_{s=-\infty}^{\infty}\sum\limits_{t=-\infty}^{\infty}|k^{\text{NCSI}}(t,s)|$
$\displaystyle=\sum\limits_{s=-\infty}^{\infty}\sum\limits_{t=-\infty}^{\infty}|\sum\limits_{k=-\infty}^{\infty}b^{2}(k)g_{0}(t-k)g_{0}(s-k)|$
$\displaystyle\leq\sum\limits_{k=-\infty}^{\infty}|b^{2}(k)|\sum\limits_{s=-\infty}^{\infty}\sum\limits_{t=-\infty}^{\infty}|g_{0}(t-k)g_{0}(s-k)|$
$\displaystyle=\sum\limits_{k=-\infty}^{\infty}b^{2}(k)\sum\limits_{s=-\infty}^{\infty}\sum\limits_{t=-\infty}^{\infty}|g_{0}(t)g_{0}(s)|$
$\displaystyle\leq\left(\sum\limits_{k=-\infty}^{\infty}b^{2}(k)\right)\left(\sum\limits_{k=-\infty}^{\infty}|g_{0}(k)|\right)^{2}<\infty,$
where the last inequality is true for
$\\{g_{0}(t)\\}_{t\in\mathbb{Z}}\in\ell_{1}(\mathbb{Z})$ and
$\\{b(t)\\}_{t\in\mathbb{Z}}$
$\in\ell_{1}(\mathbb{Z})\subset\ell_{2}(\mathbb{Z})$. Then by [28, Lemma 3] we
complete the proof.

###  A.3 Proof of Theorem 3.2 and Corollary 3.1

Before proceeding to the details, we introduce a lemma.

######  Lemma A.1.

Given three matrices $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times
1}$ and $C\in\mathbb{R}^{1\times n}$ with $A$ being non-singular, there exist
two real-valued functions $h_{j}(t),\hbar_{j}(t),t\in\mathbb{Z}$ such that

$\displaystyle CA^{t-\tau}B=\sum\limits_{j=1}^{n}h_{j}(t)\hbar_{j}(\tau),\quad
t,\tau\in\mathbb{Z},$

Proof: The proof can be found in that of [20, Theorem 1].

First, the non-causal impulse response of the nominal model (13) takes the
form of (16e). Then we have

$\displaystyle
g_{0}(t-\tau)=\left\\{\begin{array}[]{lll}&C_{c}A^{t-\tau-1}_{c}B_{c},&\quad
t-\tau\geq 1\\\ &D-C_{a}A^{-1}_{a}B_{a},&\quad t-\tau=0\\\
&-C_{a}A^{t-\tau-1}_{a}B_{a},&\quad t-\tau\leq-1\end{array}\right.$

where both $A_{c}$ and $A_{a}$ are non-singular. By Lemma A.1, there exist
four real-valued functions $h_{c,j}(t)$, $\hbar_{c,j}(t)$, $h_{a,j}(t)$,
$\hbar_{a,j}(t)$, $t\in\mathbb{Z}$ such that

$\displaystyle
C_{c}A^{t-\tau-1}_{c}B_{c}=\sum\limits_{j=1}^{n_{d,c}}h_{c,j}(t)\hbar_{c,j}(\tau),\quad
t,\tau\in\mathbb{Z},$ $\displaystyle-
C_{a}A^{t-\tau-1}_{a}B_{a}=\sum\limits_{j=1}^{n_{d,a}}h_{a,j}(t)\hbar_{a,j}(\tau),\quad
t,\tau\in\mathbb{Z},$

where $n_{d,c},n_{d,a}\in\mathbb{N}$ are the dimension of $A_{c}$ and $A_{a}$,
respectively. Referring to Definition 3.1, we first write down the generator
functions of (16), denoted by, $\mu_{j}(t),\nu_{j}(s)$, $t,s\in\mathbb{Z}$,
$j=1,\cdots,n_{d,c}+n_{d,a}$. For $j=1,\cdots,n_{d,c}$, we have

$\displaystyle\begin{array}[]{lll}&\mu_{j}(t)=h_{c,j}(t),\\\
&\nu_{j}(s)=M_{1,j}(s)+M_{2,j}(s)+M_{3,j}(s),\end{array}$ $\displaystyle
M_{1,j}(s)=\sum\limits_{j^{\prime}=1}^{n_{d,c}}h_{c,j^{\prime}}(s)\sum\limits_{i=-\infty}^{s-1}b^{2}(i)\hbar_{c,j}(i)\hbar_{c,j^{\prime}}(i),$
$\displaystyle M_{2,j}(s)=\hbar_{c,j}(s)b^{2}(s)(D-C_{a}A_{a}^{-1}B_{a}),$
$\displaystyle
M_{3,j}(s)\\!=\\!\left\\{\begin{array}[]{lll}&\sum\limits_{j^{\prime}=1}^{n_{d,a}}h_{a,j^{\prime}}(s)\sum\limits_{i=s+1}^{0}b^{2}(i)\hbar_{c,j}(i)\hbar_{a,j^{\prime}}(i),&\text{
if }s\leq-1\\\ &0,&\text{ if }s=0\\\
&-\sum\limits_{j^{\prime}=1}^{n_{d,a}}h_{a,j^{\prime}}(s)\sum\limits_{i=1}^{s}b^{2}(i)\hbar_{c,j}(i)\hbar_{a,j^{\prime}}(i),&\text{
if }s\geq 1\end{array}\right.$

and for $j=1,\cdots,n_{d,a}$, we have

$\displaystyle\begin{array}[]{lll}&\mu_{n_{d,c}+j}(t)=N_{1,j}(t)+N_{2,j}(t)+N_{3,j}(t),\\\
&\nu_{n_{d,c}+j}(s)=h_{a,j}(s),\end{array}$ $\displaystyle
N_{1,j}(t)=\sum\limits_{j^{\prime}=1}^{n_{d,a}}h_{a,j^{\prime}}(t)\sum\limits^{\infty}_{i=t+1}b^{2}(i)\hbar_{a,j}(i)\hbar_{a,j^{\prime}}(i),$
$\displaystyle N_{2,j}(t)=\hbar_{a,j}(t)b^{2}(t)(D-C_{a}A_{a}^{-1}B_{a}),$
$\displaystyle
N_{3,j}(t)\\!=\\!\left\\{\begin{array}[]{lll}&\sum\limits_{j^{\prime}=1}^{n_{d,c}}h_{c,j^{\prime}}(t)\sum\limits_{i=1}^{t-1}b^{2}(i)\hbar_{a,j}(i)\hbar_{c,j^{\prime}}(i),&\text{
if }t\geq 2\\\ &0,&\text{ if }t=1\\\
&-\sum\limits_{j^{\prime}=1}^{n_{d,c}}h_{c,j^{\prime}}(t)\sum\limits_{i=t}^{0}b^{2}(i)\hbar_{a,j}(i)\hbar_{c,j^{\prime}}(i),&\text{
if }t\leq 0\end{array}\right.$

To prove Theorem 3.2, we need to show
$k^{\text{NCSI}}(t,s)=\sum_{j=1}^{n_{d,c}+n_{d,a}}\mu_{j}(t)\nu_{j}(s)$, for
$t>s$. From (16), we have

$\displaystyle
k^{\text{NCSI}}(t,s)=\sum\limits_{i=-\infty}^{s-1}b^{2}(i)g_{0}(t-i)g_{0}(s-i)$
$\displaystyle+\sum\limits_{i=s}^{t}b^{2}(i)g_{0}(t-i)g_{0}(s-i)+\sum\limits_{i={t+1}}^{\infty}b^{2}(i)g_{0}(t-i)g_{0}(s-i).$

where

$\displaystyle\sum\limits_{i=-\infty}^{s-1}b^{2}(i)g_{0}(t-i)g_{0}(s-i)$
$\displaystyle=\sum\limits_{j=1}^{n_{d,c}}h_{c,j}(t)\left(\sum\limits_{j^{\prime}=1}^{n_{d,c}}h_{c,j^{\prime}}(s)\sum\limits_{i=-\infty}^{s-1}b^{2}(i)\hbar_{c,j}(i)\hbar_{c,j^{\prime}}(i)\right)$
$\displaystyle=\sum\limits_{j=1}^{n_{d,c}}\mu_{j}(t)M_{1,j}(s),$
$\displaystyle\sum\limits_{i=t+1}^{\infty}b^{2}(i)g_{0}(t-i)g_{0}(s-i)$
$\displaystyle=\sum\limits_{j=1}^{n_{d,a}}h_{a,j}(s)\left(\sum\limits_{j^{\prime}=1}^{n_{d,a}}h_{a,j^{\prime}}(t)\sum\limits^{\infty}_{i=t+1}b^{2}(i)\hbar_{a,j}(i)\hbar_{a,j^{\prime}}(i)\right)$
$\displaystyle=\sum\limits_{j=1}^{n_{d,a}}N_{1,j}(t)\nu_{n_{d,c}+j}(s),$
$\displaystyle\sum\limits_{i=s}^{t}b^{2}(i)g_{0}(t-i)g_{0}(s-i)=b^{2}(s)g_{0}(t-s)g_{0}(0)$
$\displaystyle\quad+\sum\limits_{i=s+1}^{t-1}b^{2}(i)g_{0}(t-i)g_{0}(s-i)+b^{2}(t)g_{0}(0)g_{0}(s-t)$
$\displaystyle=\sum\limits_{j=1}^{n_{d,c}}\mu_{j}(t)M_{2,j}(s)+\sum\limits_{j=1}^{n_{d,c}}\mu_{j}(t)M_{3,j}(s)$
$\displaystyle\quad+\sum\limits_{j=1}^{n_{d,a}}N_{3,j}(t)\nu_{n_{d,c}+j}(s)+\sum\limits_{j=1}^{n_{d,a}}N_{2,j}(t)\nu_{n_{d,c}+j}(s).$

Then note that, if any generator functions $\mu_{j}(t),\nu_{j}(s)$,
$j=1,\cdots,n_{d,c}+n_{d,a}$ vanish, the semiseparable rank will be smaller
than $n_{d,c}+n_{d,a}$. This completes the proof of Theorem 3.2.

To prove Corollary 3.1, we need to show
$k^{\text{NCSI}}(t,t)=\sum_{j=1}^{n_{d,c}+n_{d,a}}\mu_{j}(t)\nu_{j}(t)$. On
one hand, we have from (16)

$\displaystyle
k^{\text{NCSI}}(t,t)=\sum\limits_{i=-\infty}^{t-1}b^{2}(i)g_{0}(t-i)g_{0}(t-i)$
$\displaystyle+b^{2}(t)g_{0}(0)g_{0}(0)+\sum\limits_{i={t+1}}^{\infty}b^{2}(i)g_{0}(t-i)g_{0}(t-i).$
$\displaystyle=\sum\limits_{j=1}^{n_{d,c}}\mu_{j}(t)M_{1,j}(t)+b^{2}(t)(D-C_{a}A^{-1}_{a}B_{a})^{2}$
$\displaystyle+\sum\limits_{j=1}^{n_{d,a}}N_{1,j}(t)\nu_{j+n_{d,c}}(t)$ (A.9)

On the other hand, we have

$\displaystyle\sum\limits_{j=1}^{n_{d,c}+n_{d,a}}\mu_{j}(t)\nu_{j}(t)=\sum\limits_{j=1}^{n_{d,c}}\mu_{j}(t)M_{1,j}(t)+\sum\limits_{j=1}^{n_{d,c}}\mu_{j}(t)M_{2,j}(t)$
$\displaystyle+\sum\limits_{j=1}^{n_{d,c}}\mu_{j}(t)M_{3,j}(t)+\sum\limits_{j=1}^{n_{d,a}}N_{3,j}(t)\nu_{j+n_{d,c}}(t)$
$\displaystyle+\sum\limits_{j=1}^{n_{d,a}}N_{2,j}(t)\nu_{j+n_{d,c}}(t)+\sum\limits_{j=1}^{n_{d,a}}N_{1,j}(t)\nu_{j+n_{d,c}}(t)$
$\displaystyle=\sum\limits_{j=1}^{n_{d,c}}\mu_{j}(t)M_{1,j}(t)+b^{2}(t)(C_{c}A_{c}^{-1}B_{c})(D-C_{a}A_{a}^{-1}B_{a})$
$\displaystyle-b^{2}(t)(C_{c}A_{c}^{-1}B_{c})(-C_{a}A_{a}^{-1}B_{a})$
$\displaystyle+b^{2}(t)(-C_{a}A_{a}^{-1}B_{a})(D-C_{a}A_{a}^{-1}B_{a})$
$\displaystyle+\sum\limits_{j=1}^{n_{d,a}}N_{1,j}(t)\nu_{j+n_{d,c}}(t)$ (A.10)

Clearly, (A.9) and (A.10) are not equal in general unless

$\displaystyle D(D-C_{c}A_{c}^{-1}B_{c}-C_{a}A_{a}^{-1}B_{a})=0.$

This completes the proof of Corollary 3.1.

###  A.4 Closed Form of the NCSI-FO Kernel (22)

To avoid the abuse of notation, we let $k(t,s;\eta)$ represent $k^{\text{NCSI-
FO}}(t,s;\eta)$ in this section. Note that the kernel is symmetric, so we only
consider $k(t,s;\eta)$ with $t\geq s$. It follows from (16) that

$\displaystyle
k(t,s;\eta)=k_{c}(t,s;\eta)+k_{\delta}(t,s;\eta)+k_{a}(t,s;\eta),$
$\displaystyle
k_{c}(t,s;\eta)=\sum\limits_{k=1}^{\infty}\sigma^{2}_{c}\lambda_{c}^{k}g_{0}(t-k;\eta)g_{0}(s-k;\eta),$
$\displaystyle k_{\delta}(t,s;\eta)=\sigma^{2}_{0}g_{0}(t;\eta)g_{0}(s;\eta),$
$\displaystyle
k_{a}(t,s;\eta)=\sum\limits_{k=-\infty}^{-1}\sigma^{2}_{a}\lambda_{a}^{-k}g_{0}(t-k;\eta)g_{0}(s-k;\eta).$

Substituting (22d) into $g_{0}(\cdot)$, we obtain the closed form as follows.

$\displaystyle
k_{\delta}(t,s;\eta)=\left\\{\begin{array}[]{lll}\sigma_{0}^{2}c_{c}^{2}a_{c}^{t+s},&t\geq
s>0\\\ \sigma_{0}^{2}c_{c}c_{0}a_{c}^{t},&t>s=0\\\
\sigma_{0}^{2}c_{c}c_{a}a_{c}^{t}a_{a}^{-s},&t>0>s\\\
\sigma_{0}^{2}c_{0}c_{a}a_{a}^{-s},&t=0>s\\\
\sigma_{0}^{2}c_{a}^{2}a_{a}^{-t-s},&0>t\geq s\\\
\sigma_{0}^{2}c_{0}^{2},&t=s=0\end{array}\right.,$ $\displaystyle
k_{c}(t,s;\eta)=\left\\{\begin{array}[]{lll}&\sum\limits_{i=1}^{5}K_{i}(t,s),&t>s>1\\\
&\sum\limits_{i=1}^{4}K_{i}(t,s),&t>s=1\\\
&\sum\limits_{i=1}^{3}K_{i}(t,s),&t>1>s\\\
&\sum\limits_{i=1}^{2}K_{i}(t,s),&t=1>s\\\
&\sum\limits_{i=1}^{1}K_{i}(t,s),&1>t>s\\\
&K_{1}(t,t)+K_{5}(t,t)+K_{6}(t,t),&t=s>1\\\ &K_{1}(t,t)+K_{6}(t,t),&t=s=1\\\
&K_{1}(t,t),&1>t=s\\\ \end{array}\right.,$

with

$\displaystyle
K_{1}(t,s)=\sigma_{c}^{2}c_{a}^{2}a_{a}^{-t-s}\frac{(\lambda_{c}a_{a}^{2})^{\max\\{1,t+1\\}}}{1-\lambda_{c}a_{a}^{2}},$
$\displaystyle K_{2}(t,s)=\sigma_{c}^{2}c_{0}c_{a}a_{a}^{t-s}\lambda_{c}^{t},$
$\displaystyle
K_{3}(t,s)=\sigma_{c}^{2}c_{c}c_{a}a_{c}^{t}a_{a}^{-s}\frac{(\lambda_{c}a_{a}a_{c}^{-1})^{\max\\{s+1,1\\}}-(\lambda_{c}a_{a}a_{c}^{-1})^{t}}{1-\lambda_{c}a_{a}a_{c}^{-1}},$
$\displaystyle K_{4}(t,s)=\sigma_{c}^{2}c_{0}c_{c}a_{c}^{t-s}\lambda_{c}^{s},$
$\displaystyle
K_{5}(t,s)=\sigma_{c}^{2}c_{c}^{2}a_{c}^{t+s}\frac{\lambda_{c}a_{c}^{-2}-(\lambda_{c}a_{c}^{-2})^{s}}{1-\lambda_{c}a_{c}^{-2}},$
$\displaystyle K_{6}(t,t)=\sigma_{c}^{2}c_{0}^{2}\lambda_{c}^{t}.$

Note that $k_{a}(t,s;\eta)$ can be obtained in the same way as
$k_{c}(t,s;\eta)$ because

$\displaystyle k_{a}(t,s;\eta)=k_{c}(-t,-s;\eta^{*}),$
$\displaystyle\eta=[a_{c},a_{a},c_{c},c_{0},c_{a},\lambda_{c},\lambda_{a},\sigma_{c},\sigma_{a},\sigma_{0}],$
$\displaystyle\eta^{*}=[a_{a},a_{c},c_{a},c_{0},c_{c},\lambda_{a},\lambda_{c},\sigma_{a},\sigma_{c},\sigma_{0}].$

######  Remark A.1.

It is shown that $k_{c}(t,s)$, $k_{\delta}(t,s)$ and $k_{a}(t,s)$ all have the
product $\sigma^{2}_{(\cdot)}c_{(\cdot)}c^{\prime}_{(\cdot)}$ with
$\sigma_{(\cdot)}\in\\{\sigma_{c},\sigma_{a},\sigma_{0}$} and
$c_{(\cdot)},c^{\prime}_{(\cdot)}\in\\{c_{c},c_{a},c_{0}$}. It follows that

$\displaystyle k^{\text{NCSI-FO}}(t,s;\eta_{1})=k^{\text{NCSI-
FO}}(t,s;\eta_{2}),$
$\displaystyle\eta_{1}=[a_{c},a_{a},c_{c},c_{0},c_{a},\lambda_{c},\lambda_{a},\sigma_{c},\sigma_{a},\sigma_{0}],$
$\displaystyle\eta_{2}=[a_{c},a_{a},c_{c}\\!\cdot\\!\sigma_{0},c_{0}\\!\cdot\\!\sigma_{0},c_{a}\\!\cdot\\!\sigma_{0},\lambda_{c},\lambda_{a},\frac{\sigma_{c}}{\sigma_{0}},\frac{\sigma_{a}}{\sigma_{0}},1],$

This leads to the identifiability issue in the hyper-parameter estimation.
Thus, $\sigma_{0}$ is always fixed to be 1 and not treated as the hyper-
parameter.

## References

  * [1] L. A. Zadeh.  On the identification problem.  IRE Transactions on Circuit Theory, 3:277–281, 1956.
  * [2] L. Ljung.  System Identification - Theory for the User.  Prentice-Hall, Upper Saddle River, N.J., 2nd edition, 1999.
  * [3] T. Söderström and P. Stoica.  System Identification.  Prentice-Hall Int., London, 1989.
  * [4] R. Pintelon and J. Schoukens.  System Identification – A Frequency Domain Approach.  IEEE Press, New York, 2001.
  * [5] L. Ljung.  The System Identification Toolbox: The Manual.  The MathWorks Inc. 1st edition 1986, 4th edition 1995, Natick, MA, 1995\.
  * [6] G. Pillonetto, F. Dinuzzo, T. Chen, G. De Nicolao, and L. Ljung.  Kernel methods in system identification, machine learning and function estimation: A survey.  Automatica, 50(3):657–682, 2014.
  * [7] A. Chiuso.  Regularization and Bayesian learning in dynamical systems: Past, present and future.  Annual Reviews in Control, 41:24 – 38, 2016.
  * [8] L. Ljung, T. Chen, and B. Mu.  A shift in paradigm for system identification.  International Journal of Control, 93(2):173–180, 2020.
  * [9] G. Pillonetto, T. Chen, A. Chiusoand, G. De Nicolao, and L. Ljung.  Regularized System Identification: Learning Dynamic Models from Data.  Springer Nature, 2022.
  * [10] T. Chen.  On kernel design for regularized lti system identification.  Automatica, 90:109–122, 2018.
  * [11] G. Pillonetto and G. De Nicolao.  A new kernel-based approach for linear system identification.  Automatica, 46(1):81–93, 2010.
  * [12] Mattia Zorzi and Alessandro Chiuso.  The harmonic analysis of kernel functions.  Automatica, 94:125–137, 2018.
  * [13] A. Marconato, M. Schoukens, and J. Schoukens.  Filter-based regularisation for impulse response modelling.  IET Control Theory $\&$ Applications, 11:194–204, 2016.
  * [14] T. Chen, H. Ohlsson, and L. Ljung.  On the estimation of transfer functions, regularizations and gaussian processes—revisited.  Automatica, 48(8):1525–1535, 2012.
  * [15] L. Blanken and T. Oomen.  Kernel-based identification of non-causal systems with application to inverse model control.  Automatica, 114, 2020.
  * [16] K. Zhou, J. C. Doyle, and K. Glover.  Robust and optimal control.  Prentice-Hall, Inc., 1996.
  * [17] R. A. González, C. R. Rojas, and H. Hjalmarsson.  Non-causal regularized least-squares for continuous-time system identification with band-limited input excitations.  In 2021 60th IEEE Conference on Decision and Control, pages 114–119, 2021.
  * [18] Y. Fujimoto.  Kernel-based regularization for unstable systems.  In 2022 IEEE 61st Conference on Decision and Control (CDC), pages 209–214. IEEE, 2022.
  * [19] T. Chen, M. S. Andersen, L. Ljung, A. Chiuso, and G. Pillonetto.  System identification via sparse multiple kernel-based regularization using sequential convex optimization techniques.  IEEE Transactions on Automatic Control, 59(11):2933–2945, 2014\.
  * [20] T. Chen and M. S. Andersen.  On semiseparable kernels and efficient implementation for regularized system identification and function estimation.  Automatica, 132:109682, 2021.
  * [21] M. Verhaegen.  A subspace model identification solution to the identification of mixed causal, anti-causal lti systems.  SIAM Journal on Matrix Analysis and Applications, 17(2):332–347, 1996.
  * [22] J. van Zundert and T. Oomen.  On inversion-based approaches for feedforward and ilc.  Mechatronics, 50:282–291, 2018.
  * [23] C. Carmeli, E. De Vito, and A. Toigo.  Vector valued reproducing kernel hilbert spaces of integrable functions and mercer theorem.  Analysis and Applications, 4(04):377–408, 2006.
  * [24] R. Vandebril, M. van Van Barel, and N. Mastronardi.  Matrix Computations and Semiseparable Matrices: Linear Systems.  Johns Hopkins University Press, 2008.
  * [25] M. S. Andersen and T. Chen.  Smoothing splines and rank structured matrix: Revisiting the spline kernel.  SIAM Journal on Matrix Analysis and Applications, 41(2):389–412, 2020.
  * [26] M. Loève.  Probability theory, Vol. I.  Springer, 1977.
  * [27] P.K. Sen and J.M. Singer.  Large sample methods in statistics: an introduction with applications, volume 25.  CRC press, 1994.
  * [28] F. Dinuzzo.  Kernels for linear time invariant system identification.  SIAM Journal on Control and Optimization, 53(5):3299–3317, 2015\.