# Estimation of Switched Markov
Polynomial NARX models

Alessandro Brusaferria,b, Matteo Matteuccib, Stefano Spinellia,b aCNR-
Institute of Intelligent Industrial Technologies and Systems for Advanced
Manufacturing, Milan, Italy
bPolitecnico di Milano - Department of Electronics, Informatics and
Bioengineering, Milan, Italy
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This work targets the identification of a class of models for hybrid dynamical
systems characterized by nonlinear autoregressive exogenous (NARX) components,
with finite-dimensional polynomial expansions, and by a Markovian switching
mechanism. The estimation of the model parameters is performed under a
probabilistic framework via Expectation Maximization, including submodel
coefficients, hidden state values and transition probabilities. Discrete mode
classification and NARX regression tasks are disentangled within the
iterations. Soft-labels are assigned to latent states on the trajectories by
averaging over the state posteriors and updated using the parametrization
obtained from the previous maximization phase. Then, NARXs parameters are
repeatedly fitted by solving weighted regression subproblems through a
cyclical coordinate descent approach with coordinate-wise minimization.
Moreover, we investigate a two stage selection scheme, based on a $\ell
1$-norm bridge estimation followed by hard-thresholding, to achieve
parsimonious models through selection of the polynomial expansion. The
proposed approach is demonstrated on a SMNARX problem composed by three
nonlinear sub-models with specific regressors.

###### Index Terms:

Hybrid systems identification, NARX, Markov model, Expectation Maximization,
Two-stage selection, $\ell 1$-norm

##  I Introduction

Hybrid systems provide an integrated framework to model complex interactions
between discrete and continuous behaviors, as it is the case for Cyber
Physical Systems [1]. Piecewise and Switched model forms are widely considered
for identification purpose. While the former characterizes the discrete modes
by partitioning the input/state space, the latter considers arbitrarily
switching patterns, both in terms of the transition instants and discrete
state values [2].

A substantial amount of research works has been dedicated over the past
decades to the special subclass of affine dynamics. A non-exhaustive list of
the most popular approaches includes sparse regression methods, sum-of-norms
regularization, recursive methods, clustering based techniques, the bounded
error approach and the mathematical programming methods (see e.g,[2] and
references therein).

Despite being crucial to address complex applications, far fewer studies have
been dedicated to the general case of Switched nonlinear dynamical systems
[3]. A first attempt has been performed in [4] where kernel regression and
Support Vector Machines are employed, but resulting feasible just for rather
limited dataset sizes. Hence, several extensions have been investigated,
including differentiable approximations using the products of error estimator,
Least-Squares Support Vector Machines, Feature Vector Selection, Kernel
Principal Component Regression and Reduced KPCR by Incomplete Cholesky
Decomposition [5]. Authors in [6] proposed an extension of the error
sparsification approach to the nonlinear case. Then, a convex relaxation of
the $l0$-norm is employed, providing optimality guarantees for noiseless data.
Robust sparsity is introduced in [7] to address common noisy conditions,
followed by the iterative solution of Support Vector Regression subproblems.

While aforementioned studies exploit nonparametric techniques (i.e, kernels),
a parametric perspective is proposed in [3] by adopting finite-dimensional
parameterized polynomial expansions, enabling the characterization of a broad
range of nonlinear systems using a small number of parameters while improving
model interpretability. Such representations suffer curse of dimensionality,
therefore the integration of model structure selection (MSS) facilities is
required for their application in practice. An iterative method is developed
to alleviate the combinatorial complexity of the resulting mixed integer
optimization problem, where mode assignment and MSS are tuned through a
randomized strategy, followed by refinement of the switching instants.

Despite the specific class of model adopted to represent the continuous
dynamics, a common feature of the reported approaches resides in the treatment
of the switching mechanism. Indeed, the observations are typically segmented -
by inferring the hidden active states - to identify the submodels, but the
stochastic information behind the transitions is neither reconstruct nor
exploited [8]. In various application fields, as e.g, in the operation of
industrial processes, mode transitions follows certain patterns that can be
exploited to improve predictions [9]. To this end, the discrete state
evolution can be modeled by Markov chains, leading to Jump Markov Systems
[10]. As for the piecewise and switching models above, most studies in this
context have focused the linear subclasses, i.e, Jump Markov Linear Systems,
jump BoxJenkins and Switched Markov AutoRegressive eXogenous systems (see e.g,
[8],[9],[10]). Jump Markov Nonlinear Systems are investigated in [11],
proposing a recursive maximum likelihood approach based on a Rao-Blackwellized
particle filter.

In this work, we investigate a nonlinear extension of the Switched Markov
Autoregressive eXogenous systems. Specifically, we consider a parametric
approach to represent the continuous dynamics by including polynomial
expansions, thus obtaining a Switched Markov Nonlinear ARX system (labeled
SMNARX hereafter for notation simplicity). The estimation of the overall model
parameters, including NARXs, hidden state values and transition probabilities,
is performed under a probabilistic framework via Expectation Maximization.
Then, mode-wise regressors are sparsified within the disentangled weighted
regression subproblems, where the soft-assignments to the latent states are
iteratively updated through the posterior estimations until convergence.
Several MSS techniques have been proposed in the literature for NARXs,
including heuristic search, mathematical programming and sparsity promoting
terms in overparametrized models (see e.g, [12] and references therein). Here,
we focus on a two-step approach including $\ell 1$-norm relaxation of the
$\ell 0$-norm followed by thresholding.

The paper is structured as follows: Section 2 details the developed SMNARX;
Section 3 reports the developed algorithm for parameters estimation; Section 4
describes the adopted case study and summarizes the results achieved.

##  II SMNARX identification problem

Consider the class of hybrid dynamical systems defined in discrete time as
follows:

$\displaystyle\begin{dcases}y_{k}&=g_{1}(x_{k};\theta_{1})+e_{k},\ \
\text{when:}\ z_{k}=1\\\ &...\\\ y_{k}&=g_{S}(x_{k};\theta_{S})+e_{k},\
\text{when:}\ z_{k}=S\\\ \end{dcases}$ (1)

where $x_{k}=[y_{k-1},...,y_{k-n_{a}},u^{T}_{k-1},...,u^{T}_{k-n_{b}}]^{T}$ is
a vector of past system outputs $y_{k}\in\mathbb{R}$ and exogenous inputs
$u_{k}\in\mathbb{R}^{q}$ at discrete time $k\in\mathbb{Z}$, up to lags
$n_{a},n_{b}\in\mathbb{Z}^{+}$, and $e_{k}\in\mathbb{R}$ is an additive noise
term. Each discrete mode $s\in{1,...,S}$ assumes a specific deterministic
nonlinear mapping defined by a linear combination of basis functions
$\left\\{\varphi_{i}(x_{k})\right\\}_{i=1}^{n}$:

$g_{s}(x_{k};\theta_{s})=\sum_{i=1}^{n}\theta_{s,i}\varphi_{i}(x_{k})$ (2)

As introduced above, in this work we adopt a polynomial NARX form in the
states, hence the regressor vector $\varphi_{i}(.)$ is constituted by
monomials of $x_{k}$ up to a fixed degree $n_{d}\in\mathbb{Z}^{+}$.
Considering a sequence of $N$ observed outputs
$\mathcal{Y}_{N}=\left\\{y_{k}\right\\}_{k=1}^{N}$ and related lagged vectors
$\mathcal{X}_{N}=\left\\{x_{k}\right\\}_{k=1}^{N}$, and following the first-
order Markovian assumption over the latent modes:

$p(z_{k}|z_{k-1},z_{k-2},...,z_{1},\mathcal{Y}_{N-1},\mathcal{X}_{N-1})=p(z_{k}|z_{k-1})$
(3)

the joint probability over the path of the hidden states
$\mathcal{Z}_{N}=\left\\{z_{k}\right\\}_{k=1}^{N}$ factorizes as follows:

$p(\mathcal{Z}_{N})=p(z_{N},...,z_{k},...,z_{1})=\prod_{k=2}^{N}p(z_{k}|z_{k-1})p(z_{1})$
(4)

while the joint mode/observation probability at $k$ results:

$p(y_{k},z_{k}|\mathcal{Z}_{k-1},\mathcal{X}_{k})=p(y_{k}|z_{k},x_{k})p(z_{k}|z_{k-1})$
(5)

The conditional likelihood of the observation sequence $\mathcal{Y}_{N}$ is
obtained through marginalization over the hidden-state path:

$\small
p_{\Theta}(\mathcal{Y}_{N}|\mathcal{X}_{N})=\sum_{\mathcal{Z}_{N}}\prod_{k=2}^{N}p_{\theta}(y_{k}|x_{k},z_{k})p_{\mathcal{A}}(z_{k}|z_{k-1})p_{\theta}(y_{1}|x_{1},z_{1})\Pi$
(6)

where $p_{\mathcal{A}}(z_{k}|z_{k-1})$ represents components of the transition
probability matrix
$\mathcal{A}=\left\\{a^{i,j}:=p(z_{k}=j|z_{k-1}=i)\right\\}_{i,j=1}^{S}$,
$\Pi=\left\\{\pi^{i}\right\\}_{i=1}^{S}$, $\sum_{i=1}^{S}\pi^{i}=1$ the
initial state probabilities and $p_{\theta}(y_{k}|x_{k},z_{k})$ the
$\theta$-parametrized emission probabilities defined hereafter.
$\Theta=\left\\{\theta,\mathcal{A},\Pi\right\\}$ summarizes the set of SMNARX
model parameters to be estimated. The likelihoods of the output observations
$y_{k}$ are assumed independent across the time instants given the state and
the conditioning inputs at $k$. The shape of the emission distributions in the
discrete states depends on the considered noise form. In this work, we follow
the common Gaussian assumption $y_{k}|z_{k}=s,x_{k}\ \sim
N(.;g_{s}(x_{k};\theta_{s}),\sigma^{2})$, obtaining a set of $S$ subsystems
densities:

$p(y_{k}|x_{k},z_{k}=s,\theta_{s},\sigma)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{1}{2\sigma^{2}}\left(y_{k}-g_{s}(x_{k};\theta_{s})\right)^{2}}$
(7)

where the NARX predictors provide the conditional means and
$\theta=\\{\left\\{\theta_{s}\\}_{s=1}^{S},\sigma\right\\}$ summarizes the
submodels parameters.

It is worth noting that, beside being a class of models for hybrid systems,
SMNARX are a particular form of conditional probability Mixtures where
component responsibilities depend on the assignments in the preceding time
instant, referred to as Hidden Markov Experts in the time series analysis
field [13]. Besides, they constitute an extension of the Hidden Markov Models,
widely exploited in the speech recognition context, including conditioned
output distribution [14]. If the Markovian conditioning over the mode
switching mechanism is removed - hence loosing the strength of SMNARX reported
above - the active state is estimated using the normalized emission
probabilities. In fact, such simplification leads to conventional Mixture of
Experts previously employed for PWARX and timed switching ARX systems
identification [9].

The overall SMNARX identification given a collection of input/output data
$\mathcal{D}=\left\\{(x_{k},y_{k})\right\\}_{k=1}^{N}$ is a highly difficult
task involving discrete mode inference over the observed trajectories, the
estimation of the transition probabilities and NARX parameters
$\left\\{\theta_{s}\right\\}_{s=1}^{S}$, the selection of the components in
the regressor vectors $\left\\{\varphi_{i}(.)\right\\}_{s=1}^{S}$, as well as
model orders $n_{a},n_{b}$ and the number of modes $S$. In this work, we
follow the common assumption of predetermined number of modes and maximum
order/lag, often inferred in practice by preliminary data analysis, cross
validation and order selection techniques [2]. Despite such simplification,
the identification problem is still challenging, mainly due to the unknown
active states and proper regressors shapes. To tackle this issues, we
implemented the estimation method reported in the next section.

##  III Parameters estimation method

###  III-A Expectation Maximization approach

While a brute force maximization of the conditional likelihood
$p_{\Theta}(\mathcal{Y}_{N}|\mathcal{X}_{N})$ would be impracticable for
realistic problem sizes, a much more efficient solution can be achieved under
the Expectation Maximization (EM) framework. Constituting a general purpose
technique for parameters estimation in case of missing data, EM can be
employed to address the latent modes of switching Markovian models.
Specifically, the likelihood optimization problem is approached by iterating
over the following Expectation and Maximization steps:

$\displaystyle\begin{aligned} E&-step:\\\
&\mathcal{Q}\left(\Theta,\Theta^{(t)}\right)=\mathbb{E}_{\mathcal{Z_{N}}}\left[p(\mathcal{Y}_{N},\mathcal{Z}_{N}|\mathcal{X}_{N},\Theta)|\mathcal{Y}_{N},\mathcal{X}_{N},\Theta^{(t)}\right]\\\
M&-step:\end{aligned}$ $\Theta^{(t+1)}=\underset{\Theta}{\small{argmax}}\
\mathcal{Q}\left(\Theta,\Theta^{(t)}\right)$ (8)

where $\mathcal{Q}(\Theta,\hat{\Theta})$ is a deterministic auxiliary function
representing the expectation over the latent states of the complete data
likelihood (i.e, assuming the latter to be given) determined using the model
parameters from the previous iteration $t$. Hence, the $\mathcal{Q}$-function
evaluation involves the computation of the posterior distribution of the
latent states, as shown here:

$\mathcal{Q}(\Theta,\Theta^{(t)})=\sum_{\mathcal{Z_{N}}}p(\mathcal{Z}_{N}|\mathcal{Y}_{N},\mathcal{X}_{N},\Theta^{(t)})p(\mathcal{Y}_{N},\mathcal{Z}_{N}|\mathcal{X}_{N},\Theta)$

Practically, it works as a kind of ”approximated oracle”: soft-labels are
assigned to latent states over the trajectories by averaging over the
posterior, updated using the parametrization obtained from the previous
maximization phase, starting from initial guesses[14]. The algorithm iterates
until convergence, tested by stopping conditions on the likelihood or maximum
number of iterations. As common for EM, several initialization are required,
by sampling different values of the parameters to tackle eventual convergence
to poor local minima. Notably, EM has interesting convergence properties,
holding also when a closed form solution is not available for the M-step [15].
In this case, an increase of the $\mathcal{Q}$-function is required in each
iteration, leading to the Generalized EM algorithm (GEM).

###  III-B Computation of the latent state posterior in the E-step

As the probability graph of the SMNARX has a tree structure similar to Hidden
Markov Models (and Hidden Markov Experts), the posterior can be computed
efficiently through a slight modification of a widely adopted message passing
algorithm, often referred to as forward-backward or Baum-Welch algorithm
[15],[13]. Specifically, two auxiliary variables are introduced to recursively
perform the estimation step, namely the forward and backward variables. By
introducing a specific term for the probability of being in state $s$ at $k$
after seeing the first t observations (given the model parameters):

$p(z_{k}=i,y_{k},...,y_{1}):=\alpha_{k}^{i}$ (9)

the probability of the overall sequences of observations is obtained by
marginalizing the states occurring at the final step:

$p(\mathcal{Y}_{N})=\sum_{i=1}^{S}p(z_{N}=i,y_{N},..,y_{k},..,y_{1})=\sum_{i=1}^{S}\alpha_{N}^{i}$
(10)

We hided explicit references to the conditioning parameters to lighten
notation. $\alpha_{k}^{i}$, typically referred to as the forward variable in
the HMM context, can be computed efficiently using a recursive procedure,
starting from initial assignments $\alpha_{1}^{i}=\pi^{i}b_{1}^{i}$, through
the following factorization:

$\displaystyle\begin{aligned}
\alpha_{k}^{j}=&\sum_{i=1}^{S}[p(z_{k-1}=i,y_{k-1},...,y_{1})p(z_{k}=j|z_{k-1}=i)\\\
&*p(y_{k}|x_{k},z_{k}=j)]=\sum_{i=1}^{S}\left(\alpha_{k-1}^{i}a^{i,j}\right)b_{k}^{j}\end{aligned}$
(11)

where $b_{k}^{j}=p(y_{k}|x_{k},z_{k}=j)$ represents the emission probability
given the active state $z_{k}=j$. It is worth noting that such operations
result linear in time, as opposed to the exponential increase with the
sequence length given by full paths coverage. The sequences are processed in
the opposite direction by the backward variables, representing the joint
probability of the observations from $k+1$, given the active state at $k$:

$p(y_{N},...,y_{k+1}|z_{k}=i):=\beta_{k}^{i}$ (12)

computed through the following recursion, starting from $N$:

$\displaystyle\begin{aligned}
\small\beta_{k}^{i}=&\sum_{j=1}^{S}[p(z_{k+1}=j|z_{k}=i)p(y_{k+1}|x_{k+1},z_{k+1}=j)\\\
&*p(y_{N},...,y_{k+2}|z_{k+1}=j)]=\sum_{j=1}^{S}a^{i,j}b_{k+1}^{j}\beta_{k+1}^{j}\end{aligned}$
(13)

Hence, the overall observation path (summarized by
$\alpha_{k}^{s},\beta_{k}^{s}$) is employed to compute the state posterior at
time $k$:

$p(z_{k}=j|\mathcal{Y}_{N})=\frac{\alpha_{k}^{j}\beta_{k}^{j}}{\sum_{s=1}^{S}\alpha_{k}^{s}\beta_{k}^{s}}:=\gamma_{k}^{j}$
(14)

while the joint posterior distribution of two successive latent states over
the path is obtained by:

$p(z_{k}=i,z_{k+1}=j|\mathcal{Y}_{N})=\frac{\alpha_{k}^{i}a^{i,j}b_{k+1}^{j}\beta_{k+1}^{j}}{\sum_{s=1}^{S}\alpha_{k}^{s}\beta_{k}^{s}}:=\xi_{k}^{i,j}$
(15)

As the computations of the auxiliary variables for recursive estimation
include also future observations over the training sequence (by
$\beta_{k+1}^{j}$), the state posterior $\gamma_{k}^{j}$ cannot be considered
for prediction in testing conditions. To this end, the following predictive
mode probabilities can be adopted:

$p(z_{k}=j|\mathcal{Y}_{k-1})=\frac{\sum_{i=1}^{S}\alpha_{k-1}^{i}a^{i,j}}{\sum_{i=1}^{S}\sum_{j=1}^{S}\alpha_{k-1}^{i}a^{i,j}}:=f_{k}^{j}$
(16)

inferring the latent state at $k$ by iterating the forward variable till $k-1$
(giving the modes probabilities in the previous instant), times the consequent
transitions. The predicted means are then obtained by combining the modes
through $f_{k}^{s}$:

$\hat{y}_{k}=\sum_{s=1}^{S}f_{k}^{s}\sum_{i=1}^{n}\theta_{s,i}\varphi_{i}(x_{k})$
(17)

On the other hand, the posteriors $\gamma_{k}^{s},\xi_{k}^{i,j}$ provided by
the forward-backward procedure are employed to update the parameters in the
M-step, as detailed in the next subsection.

###  III-C Model parameters updates in the M-step

The probability matrix is computed by the expected number of transitions from
state i to state j on the overall the expected number of transitions from
state i.

$\hat{a}^{i,j}=\frac{\sum_{n=1}^{N-1}\xi_{k}^{i,j}}{\sum_{k=1}^{N-1}\sum_{s=1}^{S}\xi_{k}^{i,s}}$
(18)

while the initial probability parameter is assigned to
$\pi^{i}=\gamma_{1}^{i},\ \forall i$. The variance parameter is computed by:

$\hat{\sigma^{2}}=\sum_{s=1}^{S}\sum_{k=1}^{N}\gamma_{k}^{s}\left(y_{k}-\sum_{i=1}^{n}\theta_{s,i}\varphi_{i}(x_{k})\right)^{2}/\sum_{s=1}^{S}\sum_{k=1}^{N}\gamma_{k}^{s}$
(19)

The update of the NARXs parameters result in the following set of subproblems,
each weighted by the related sample-wise soft-assighments to the modes
$\gamma_{k}^{s}$ during the iterations:

$\mathcal{L}_{s}(\theta_{s})=\sum_{k=1}^{N}\gamma_{k}^{s}(y_{k}-\sum_{i=1}^{n}\theta_{s,i}\varphi_{i}(x_{k}))^{2},\
\forall s\in S$ (20)

Notably, discrete mode classification and NARX regression tasks are
disentangled within the EM iterations.

As reported in Section I, automated mechanisms has to be introduced to select
the most suitable NARX components (in terms of contribution to the prediction
performances), addressing the curse of dimensionality of polynomial expansions
as well as the increased overfitting potential of overparametrized models.
Hence, the $\theta_{s},\forall s\in S$ are assumed to attain specific sparse
patterns (i.e, most coefficients equal zero) as the corresponding basis terms
do not provide sensible contributions to the predictions, resulting in
parsimonious models.

To this end, the best subset estimation approach considers a $l0$-norm on the
parameters vector, thus constraining the regression fitting with at most
$\kappa$ non-vanishing components. Considering the specific characteristics of
the SMNARX model, this leads to a set of $S$ sub-problems of the form:

$\underset{\theta_{s}}{\small{min}}\ \mathcal{L}_{s}(\theta_{s})\ s.t.\
||\theta_{s}||_{0}\leq\kappa;\
||\theta_{s}||_{0}=\sum_{i=1}^{n}\mathbf{1}\left\\{\theta_{s,i}\neq
0\right\\}$ (21)

A sensible amount of studies has been dedicated to the achievement of scalable
and computational efficient solution of such NP-hard problems [16]. Popular
techniques include greedy forward/backward stepwise selection [17], randomized
algorithms [3], Mixed Integer Programming based methods [18] and convex
relaxations through bridge regression [19].

In this work we focus on the latter class, leaving the investigation of
alternative approaches to future developments. Specifically, we consider the
closest convex surrogate of the $l0$-norm, i.e, the $l1$-norm (often referred
to as least absolute shrinkage and selection operator - LASSO), widely
employed in the machine learning field to introduce sparsity priors and
regularization in high dimensional settings [20]. Adopting the conventional
Lagrangian form, the discrete mode weighted regressions are augmented by
sparsity promoting terms:

$\underset{\theta_{s}}{\small{min}}\
\sum_{k=1}^{N}\gamma_{k}^{s}\left(y_{k}-\sum_{i=1}^{n}\theta_{s,i}\varphi_{i}(x_{k})\right)^{2}+\lambda\sum_{i=1}^{n}|\theta_{s,i}|$
(22)

where $\lambda\geq 0$ is a tuning hyperparameter weighting the regularization.
Practically, such penalized regression form encourages sparse solutions while
performing parameters estimation, resulting in a soft-thresholding behavior,
as opposed to the hard-thresholds of the $\ell 0$-norm based subset
estimation. While LASSO can theoretically uncover sparse structures under
specific strong conditions, it suffers several limitations in finite sample
settings due to its intrinsic downward bias caused by the shrinkage action
[16]. Hence, several variants have been investigated to obtain weaker
requirements than the irrepresentable condition of the conventional LASSO for
variable selection consistency [19]. In this work, we focus on the two stage
selection scheme developed in [21], post-processing the parameters vectors
computed by the $\ell 1$-regularized regressions component-wise by the hard-
thresholding function (i.e., $\ell 1$2S):

$\small\bar{\theta}_{s}=\eta_{0}(\hat{\theta}_{s};\upsilon^{2}/2),\text{with:
}\eta_{0}(x;\varUpsilon)=x\mathbf{1}\left\\{|x|\geq\sqrt{2\varUpsilon}\right\\}$
(23)

where $\hat{\theta}_{s}$ represent the parameters estimate obtained during
each EM iteration by the solution of (21) and $\upsilon\geq 0$ is a threshold
hyperparameter. Hence, sign consistency can be achieved under weaker
requirements, while the first stage is exploited to provide good
approximations in the $\ell 2$-sense. A detailed comparison of two-stage
variable selection techniques to LASSO under various Signal to Noise ratios is
reported in [19], showing that the optimal configuration of the regularization
weight in the bridge regression step is the same for estimation and selection.
The sparse patterns are then recovered under the assumption that proper
nonzero coefficients are sufficiently large (in absolute value), controlling
the discovery trade-off by tuning $\upsilon$ [21]. Within the different
approaches that can be exploited for hyperparameters calibration (e.g,
combining the two-step procedure with the knockoff framework), in this work we
adopt the decreasing search procedure deployed in [19] with configurable
window/grid size through cross-validation. Considering the specific
characteristics of the SMNARX identification problem at hand, requiring both
parameter estimation, selection and discrete mode inference, we found useful
experimentally (as shown in the next section) to perform a burn-in period
using only the $\ell 1$-norm to get rid of the soft-assignments performed
during the first parts of the EM iterations. Hence, the hard-thresholding step
is introduced once the likelihood reach stable conditions, controlled by a
parametrized tolerance. The $\ell 1$ action is then modulated by the selected
features. To solve the weighted regressions, we adopted cyclical coordinate-
wise minimization [22].

##  IV Numerical results

To the best of our knowledge, this work constitute the first investigation of
Switched Markov Polynomial NARX systems in the literature, and a benchmark
problem for this systems has not been proposed yet. Hence, we constructed the
following system by merging the characteristics of the SNARX problem deployed
in [3] (aimed to investigate switched polynomial NARX) with the SMARX system
considered in [9]. The resulting SMNARX system to be identified is composed by
three different nonlinear sub-models with specific regressors:

$\displaystyle\begin{dcases}y_{k}&=0.5y_{k-1}+0.8u_{k-2}+u_{k-1}^{2}-0.3y_{k-2}^{2}+e_{k}\\\
y_{k}&=0.2y_{k-1}^{3}-0.5y_{k-2}-0.7y_{k-2}u_{k-2}^{2}+0.6u_{k-2}^{2}+e_{k}\\\
y_{k}&=0.5y_{k-2}-0.4y_{k-1}+0.2u_{k-1}-0.4u_{k-3}y_{k-1}+e_{k}\end{dcases}$

with $e_{k}\sim N(0,0.1^{2})$ and input sampled from a uniform distribution
$u_{k}\sim U[-1,1]$. The transition probability matrix governing the Markovian
discrete mode switching has been defined as
$A=[[0.98,0.02,0],[0,0.98,0.02],[0.02,0,0.98]]$.

We generated a sequence of 12000 samples, thus including a sensible number
($>$ 300) of state switches in the identification set. The first 10000 has
been devoted to training, the successive 1000 to validation and the last 1000
constitute the test set. Following [3], we set the orders to oversized values
$n_{a},n_{b}=4,n_{d}=3$ resulting in 165 regressors. Then, we structured the
training sequence into mini-batches of length 200.

Both the overall model architecture and the EM-based estimation algorithm have
been developed in Numpy-1.18. The weighted least squares subproblems,
including the $\ell 1$-norm and thresholding steps are implemented using the
LinearModel and SelectFromModel classes of Scikit-learn-0.23.

To avoid convergence to poor local minimum (where the states collapse to
common, i.e., average, modes), we found useful experimentally to initialize
$\gamma_{k}^{s}\sim U[0.31,0.35],\forall s,k$, while constraining a valid
distribution parametrization. Moreover, the EM is executed (for each
hyperparameters configuration) by starting from 10 random initialization and
then selecting the run with the higher likelihood. The tolerances on the
likelihood improvements have been set to $1e^{-2}$ and $1e^{-6}$ to control
burn-in and training stops respectively, including a maximum number of
iterations of 100. Prediction accuracy is measured by the conventional Root
Mean Squared Error (RMSE). Using the decreasing grid search procedure in
cross-validation over a window $[1e^{-6},1e^{1}]$, $\lambda=5e^{-4}$ resulted
a reasonable choice for test purpose. Table II reports the obtained validation
RMSE for the most significant subset of $\lambda$ in the search window. The
threshold have been set to $\upsilon=5e^{-2}$.

To evaluate the accuracy of the SMNARX reconstruction we employed the
performance measures adopted in [23][9]:

$\mathcal{F}_{\theta}=\frac{1}{S}\sum_{s=1}^{S}\left(1-\frac{\lVert\theta_{s}-\hat{\theta}_{s}\lVert}{\lVert\theta_{s}\lVert}\right),\
\ \mathcal{F}_{A}=1-\frac{\lVert\hat{A}-A\lVert}{\lVert A\lVert}$ (24)

namely the parameters and transition probability indexes. To investigate the
mode estimation performance, we adopted the index

$\mathcal{F}_{s}=1/N_{t}\sum_{k=1}^{N_{t}}\mathbf{1}_{\hat{s}_{k}=s_{k}}$

where $\hat{s}_{k}$ and $s_{k}$ represents the inferred and true state at $k$
while $N_{t}$ the train/test size, where the subscript $s_{tr}(s_{te})$ in
Table II refers to training (test) set.
As the order of the discrete states in the model do not necessarily match that
of the target system, post-processing is required before evaluation [24]. As
suggested in [23], without loss of generality, we reordered the sub-models
estimate by the Euclidean norm to the target parameters. Following [25],
statistics are computed over 100 runs of the estimation algorithm on randomly
sampled sequences using the same hyperparameters. Under the adopted
configuration, convergence have been observed in approximately 20 iterations
on average. The mean computation time of one iteration on a CPU-i7-2.5-GHz-
RAM-8Gb for this case study has been 0.8 seconds.

TABLE I: Validation RMSE over different $\lambda$

$\lambda$ | 1e-1 | 5e-2 | 1e-2 | 5e-3 | 1e-3 | 5e-4 | 1e-4
---|---|---|---|---|---|---|---
| 0.5013 | 0.3382 | 0.2007 | 0.1757 | 0.1692 | 0.1691 | 0.1691

TABLE II: Performance indexes and parameters statistics (100 runs)

| $\mathcal{F}_{s_{tr}}$ | $\mathcal{F}_{s_{te}}$ | $\mathcal{F}_{A}$ | $\mathcal{F}_{\theta}$ | $n_{feat}$
---|---|---|---|---|---
SMNARX-EM | 0.995 | 0.966 | 0.995 | 0.941 | 165
SMNARX-EM-$\ell 1$ | 0.996 | 0.967 | 0.995 | 0.969 | 165
SMNARX-EM-$\ell 1$2S | 0.996 | 0.967 | 0.995 | 0.990 | 4

$\ell 1$2S | $Mean$ | $Std$
---|---|---
$\theta_{1}$ | [0.500, 0.800, 0.999, -0.299] | [2.3e-3, 3.1e-3, 4.2e-3, 1.9e-3]
$\theta_{2}$ | [0.200, -0.500, -0.698, 0.599] | [5.3e-3, 6.0e-3, 1.4e-2, 4.0e-3]
$\theta_{3}$ | [0.498, -0.400, 0.200, -0.399] | [6.7e-3, 7.7e-3, 2.9e-3, 9.3e-3]
$\sigma^{2}$ | 0.010 | 1.3e-4

Figure 1: Predicted (blue) vs True (orange) modes over the test set

Table II reports the measured indexes and the submodels parameters statistics
computed over 100 runs, comparing the integration of the shrinkage action and
the thresholding step. Notably, the EM-based estimation technique achieved a
proper classification of the active states and switching patterns, as depicted
also by Figure 1, reporting the predicted modes over the test set.
Furthermore, the integration of the two stage selection utility enabled the
extraction of the correct set of features ($n_{feat}$) - for each submodel -
from the sparsified patterns provided by the $\ell 1$-norm bridge estimator,
and the consequent improvement of NARXs submodel parameters fitting
($\mathcal{F}_{\theta}$). To provide further insights into such behavior, we
report in Figure 2 the plots of the coefficients paths of the NARX over the
iteration before/after the threshold insertion on the conventional LASSO,
starting from random initialization. Visibly, the base $\ell 1$-norm
sparsifier, employed during the burn-in phase, provides foreseeable
approximations [21], where false positive terms have been assigned with a very
small parameter value. Once the threshold phase is activated, such minor
components are incrementally removed, while the coefficients of the remaining
parameters are slightly adjusted accordingly, until convergence.

Clearly, a proper hyperparameter configuration is crucial to achieve such
behavior. While the adopted cross-validation based procedure provided suitable
results on this case study, this topic would require further investigation.
Chiefly, enhanced heuristics/search techniques might further speed-up tuning,
fostering scalability to increasing problem complexity.

Figure 2: NARXs parameters during iterations for each mode

##  V Conclusion and Next Steps

In this work we investigated the identification of Switched Markov Nonlinear
ARX system with parameterized polynomial expansions. To tackle the challenging
brute force likelihood maximization, estimation is performed using an
Expectation Maximization algorithm, disentangling discrete mode classification
and NARXs regression tasks. Sparse features patters are enforced through a two
stage approach including a $\ell 1$-norm based bridge estimation. Numerical
analysis has been performed on a benchmark SMNARX problem, showing the
capability to achieved a proper identification of the active states, switching
patterns, and submodel parameters.

Next developments will include the integration of enhanced search techniques
for hyperparameters tuning, the investigation of further model structure
selection approaches, modes/order selection methods, and the application to a
real case study.

## References

  * [1] Y. Yuan, X. Tang, W. Zhou, W. Pan, X. Li, H.-T. Zhang, H. Ding, and J. Goncalves, “Data driven discovery of cyber physical systems,” Nature Communications, vol. 10, no. 1, p. 4894, 2019.
  * [2] F. Lauer and G. Bloch, Hybrid system identification: Theory and algorithms for learning switching models, vol. 478.  Springer, 2018.
  * [3] F. Bianchi, M. Prandini, and L. Piroddi, “A randomized two-stage iterative method for switched nonlinear systems identification,” Nonlinear Analysis: Hybrid Systems, vol. 35, p. 100818, 2020.
  * [4] F. Lauer and G. Bloch, “Switched and piecewise nonlinear hybrid system identification,” in Hybrid Systems: Computation and Control (M. Egerstedt and B. Mishra, eds.), 2008.
  * [5] V. L. Le, G. Bloch, and F. Lauer, “Reduced-size kernel models for nonlinear hybrid system identification,” IEEE Transactions on Neural Networks, vol. 22, no. 12, pp. 2398–2405, 2011.
  * [6] L. Bako, K. Boukharouba, and S. Lecoeuche, “An ℓ0−ℓ1 norm based optimization procedure for the identification of switched nonlinear systems,” in 49th IEEE Conf. on Decision and Control (CDC), 2010.
  * [7] V. L. Le, F. Lauer, L. Bako, and G. Bloch, “Learning nonlinear hybrid systems: From sparse optimization to support vector regression,” in Proc. of the 16th Int. Conf. on Hybrid Systems: Computation and Control, HSCC ’13, (New York, NY, USA), p. 33–42, 2013.
  * [8] D. Piga, V. Breschi, and A. Bemporad, “Estimation of jump box–jenkins models,” Automatica, vol. 120, p. 109126, 2020.
  * [9] X. Jin and B. Huang, “Identification of switched markov autoregressive exogenous systems with hidden switching state,” Automatica, vol. 48, no. 2, pp. 436 – 441, 2012.
  * [10] A. Garulli, S. Paoletti, and A. Vicino, “A survey on switched and piecewise affine system identification,” 16th IFAC Symposium on System Identification, vol. 45, no. 16, 2012.
  * [11] E. Özkan, F. Lindsten, C. Fritsche, and F. Gustafsson, “Recursive maximum likelihood identification of jump markov nonlinear systems,” IEEE Transactions on Signal Processing, vol. 63, no. 3, pp. 754–765, 2015.
  * [12] G. C. Calafiore, L. E. Ghaoui, and C. Novara, “Sparse identification of polynomial and posynomial models,” IFAC Proceedings Volumes, vol. 47, no. 3, pp. 3238 – 3243, 2014.  19th IFAC World Congress.
  * [13] A. S. Weigend and S. Shi, “Hidden markov experts,” in Quantitative Analysis In Financial Markets, ch. 2, pp. 35–70, 2001.
  * [14] Y. Bengio and P. Frasconi, “An input output hmm architecture,” Advances in Neural Information Processing Systems, vol. 7, 12 1995.
  * [15] C. M. Bishop, Pattern Recognition and Machine Learning (Information Science and Statistics).  Berlin, Heidelberg: Springer-Verlag, 2006.
  * [16] M. Yin, N. Ho, B. Yan, and M. Zhou, “Probabilistic best subset selection by gradient-based optimization,” arXiv:2006.06448, June2020.
  * [17] T. Hastie, R. Tibshirani, and R. J. Tibshirani, “Extended comparisons of best subset selection, forward stepwise selection, and the lasso,” arXiv:1707.08692v2, 2017.
  * [18] D. Bertsimas, A. King, and R. Mazumder, “Best subset selection via a modern optimization lens,” The Annals of Statistics, vol. 44, 07 2015.
  * [19] S. Wang, H. Weng, and A. Maleki, “Which bridge estimator is optimal for variable selection?,” CoRR, vol. abs/1705.08617, 2017.
  * [20] I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning.  MIT Press, 2016.  http://www.deeplearningbook.org.
  * [21] N. Meinshausen and B. Yu, “Lasso-type recovery of sparse representations for high-dimensional data,” The Annals of Statistics, vol. 37, 2008.
  * [22] R. Tibshirani, T. Hastie, and J. Friedman, “Regularized paths for generalized linear models via coordinate descent,” Journal of Statistical Software, vol. 33, 02 2010.
  * [23] G. Pillonetto, “A new kernel-based approach to hybrid system identification,” Automatica, vol. 70, pp. 21 – 31, 2016.
  * [24] S. Fruhwirth-Schnatter, G. Celeux, and C. P. Robert, Handbook of Mixture Analysis.  Chapman and Hall/CRC, 2019.
  * [25] A. Hartmann, J. M. Lemos, R. S. Costa, J. Xavier, and S. Vinga, “Identification of switched arx models via convex optimization and expectation maximization,” Journ. of Process Control, vol. 28, 2015.