# Decoupling P-NARX models using filtered CPD

Jan Decuyper  David Westwick  Kiana Karami  Johan Schoukens  Department of
Engineering Technology, Vrije Universiteit Brussel, Pleinlaan 2, 1050
Brussels, Belgium (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Electrical and Computer Engineering, University of Calgary,
Calgary, Canada (e-mail<EMAIL_ADDRESS> Department of Electrical
Engineering and Electrical Engineering Technology, Penn state Harrisburg,
Middletown, PA, USA, 17057 (e-mail<EMAIL_ADDRESS> Department of Electrical
Engineering, Eindhoven University of Technology, Eindhoven, The Netherlands

###### Abstract

Nonlinear Auto-Regressive eXogenous input (NARX) models are a popular class of
nonlinear dynamical models. Often a polynomial basis expansion is used to
describe the internal multivariate nonlinear mapping (P-NARX). Resorting to
fixed basis functions is convenient since it results in a closed form solution
of the estimation problem. The drawback, however, is that the predefined basis
does not necessarily lead to a sparse representation of the relationship,
typically resulting in very large numbers of parameters. So-called decoupling
techniques were specifically designed to reduce large multivariate functions.
It was found that, often, a more efficient parameterisation can be retrieved
by rotating towards a new basis. Characteristic to the decoupled structure is
that, expressed in the new basis, the relationship is structured such that
only single-input single-output nonlinear functions are required. Classical
decoupling techniques are unfit to deal with the case of single-output NARX
models. In this work, this limitation is overcome by adopting the filtered CPD
decoupling method of Decuyper et al. (2021b). The approach is illustrated on
data from the Sliverbox benchmark: measurement data from an electronic circuit
implementation of a forced Duffing oscillator.

###### keywords:

Polynomial-NARX, decoupling, filtered CPD, model reduction, nonlinear system
identification

††thanks: This work was supported by the Flemish fund for scientific research
FWO under license number G0068.18N.

##  1 Introduction

Nonlinear auto-regressive exogenous input (NARX) models have been extensively
used to describe nonlinear systems. It is a black-box system identification
technique which has proven to be useful in a wide range of applications (Chan
et al., 2015; Zhao et al., 2013). NARX models describe dynamical nonlinear
behaviour by relating the current output sample to both past output samples,
and current and past input samples. Defining

$\bm{x}=\\{u(t),u(t-1),\dots,u(t-n_{u}),\\\
y(t-1),\dots,y(t-n_{y})\\},\quad\quad$ (1)

a general single-output NARX model is described by (Billings, 2013)

$y(t)=F(\bm{x})+e(t),$ (2)

where $F:\mathbb{R}^{n_{u}+n_{y}+1}\rightarrow\mathbb{R}$ is a static
multiple-input single-output (MISO) nonlinear function and $e(t)$ is an
equation error which is assumed to be a sequence of independent identically
distributed (IID) random variables. The function $F$ may be described by any
class of functions, e.g. it can be a neural network or a wavelet network.
Often, however, a basis expansion is preferred. In that case, a direct
estimate of the model parameters follows from linear regression by using the
measured outputs in the regressor (minimising the equation error). A popular
choice is the polynomial basis, leading to so-called Polynomial-NARX or P-NARX
models.

The ease of identification is however countered by a number of fundamental
disadvantages. P-NARX models are notorious for the number of parameters that
have to be identified, growing both with the number of past inputs and past
outputs, $n=n_{u}+n_{y}+1$, and combinatorially with the nonlinear degree $d$.
Moreover, being a black-box method, the obtained model is often very hard to
interpret. General practice is to use regularisation to steer the optimisation
towards meaningful terms and hence reduce the number of parameters. Often,
however, no sparse representation can be obtained using the classical monomial
basis. The present work promotes the use of decoupling techniques, which were
designed to translate the relationship into a more favourable basis.

Similar issues have already been addressed in the context of nonlinear state-
space models, where typically a multiple-input multiple-output (MIMO)
polynomial is used to describe the nonlinearity. It was found that by
_decoupling_ the multivariate polynomial into a number of univariate
polynomials, both insight and model reduction could be obtained (Decuyper et
al., 2021a).

In general, decoupling techniques aim at transforming generic multivariate
nonlinear functions into decoupled functions. The decoupled structure is
characterised by the fact that the nonlinear relationship is described by a
number of univariate functions of intermediate variables. Given a generic
nonlinear function

$\bm{q}=\bm{f}(\bm{p})$ (3)

with $\bm{q}\in\mathbb{R}^{m}$ and $\bm{p}\in\mathbb{R}^{n}$, the idea is to
introduce an appropriate linear transformation of $\bm{p}$, denoted $\bm{V}$,
such that in this alternative basis, univariate functions may be used to
describe the nonlinear mapping. The decoupled function is then of the
following form

$\bm{f}(\bm{p})=\bm{W}\bm{g}(\bm{V}^{\text{T}}\bm{p})$ (4)

where the $i$th function is $g_{i}(z_{i})$ with
$z_{i}=\bm{v}_{i}^{\text{T}}\bm{p}$, emphasising that all functions are
strictly univariate. The number of allowed univariate functions, denoted $r$,
is a user choice which can be used to control the model complexity. In some
cases prior knowledge of the system may dictate that a certain number of
nonlinear components drive the nonlinearity. Whenever no such knowledge is
available, a scan over $r$ is performed. The number of univariate functions
plays a crucial role since it will determine whether the implied equivalence
of Eq. (4) exists. A second linear transformation $\bm{W}$, maps the function
back onto the outputs. The linear transformations then have the following
dimensions: $\bm{V}\in\mathbb{R}^{n\times r}$ and
$\bm{W}\in\mathbb{R}^{m\times r}$.

The original decoupling procedure requires the function to be of the MIMO-type
(Dreesen et al., 2015), i.e $m>1$, excluding the class of single-output P-NARX
models. The method is based on the canonical polyadic decomposition (CPD) of a
3-way tensor, constructed out of evaluations of the Jacobian matrix along a
number of operating points (see Section 3.2). It exploits the uniqueness
properties of the CPD in order to retrieve estimates of the univariate
functions $g_{i}$. In Westwick et al. (2018), the issue with single-output
functions was circumvented by resorting to evaluations of the Hessian, again
retrieving a 3-way tensor. In practice, however, the uniqueness of the CPD is
often not guaranteed (Decuyper et al., 2019), leading to very noisy estimates
of $g_{i}$ or its derivatives. As a solution, Karami et al. (2021) proposed to
add polynomial constraints when factoring the Hessian. Accurate decoupled
models could be obtained, although at the price of computing the Hessian.

In this work, the filtered CPD approach of Decuyper et al. (2021b) will be
used. It acts on the basis of first order derivative information, avoiding the
computation of the Hessian, and resorts to non-parametric soft constraints to
ensure meaningful estimates of $g_{i}$ are obtained.

##  2  Filtered CPD decoupling

The filtered CPD approach (F-CPD) links the original function to a decoupled
function on the basis of its first order derivative information. The method
relies on the underlying diagonal structure of the Jacobian (Eq. (5)), which
follows from the use of univariate functions $g_{i}$. Denoting the left and
right hand side Jacobian of Eq. (4) by $\bm{J}$ and $\bm{J}^{\prime}$,
respectively, we have that

$\bm{J}^{\prime}=\bm{W}\operatorname{diag}\left(\begin{bmatrix}h_{1}(z_{1})\quad\cdots\quad
h_{r}(z_{r})\end{bmatrix}\right)\bm{V}^{\text{T}}$ (5)

in which case $h_{i}(z_{i})\coloneqq\frac{dg_{i}(z_{i})}{dz_{i}}$ represents
the derivative of the univariate function $g_{i}(z_{i})$ with respect to its
argument. Evaluating this Jacobian in $N$ operating points, i.e. for
$\\{\bm{p}[1],\cdots,\bm{p}[N]\\}$, and collecting the evaluations in a third
dimension, expands the object into a three-way array
$\mathcal{J}^{\prime}\in\mathbb{R}^{n\times m\times N}$. This is illustrated
graphically in Fig. 1. Notice that the diagonal plane,
$\bm{H}\in\mathbb{R}^{N\times r}$, stores evaluations of $h_{i}$, i.e the
derivative of the functions $g_{i}$. Given the diagonal form, the collection
of Jacobians may be written as a sum of $r$ outer products (or rank-one
terms). Element-wise we have that,

$\mathcal{J}^{\prime}_{\left[s,k,l\right]}=\sum_{i=1}^{r}w_{si}~{}v_{ki}~{}h_{li},$
(6)

$\text{for}~{}s=1,\dots,m~{}\text{and}~{}k=1,\dots,n~{}\text{and}~{}l=1,\dots,N$.
A sum of rank-one terms defines a diagonal tensor decomposition (Kolda and
Bader, 2009). The latter is illustrated on the right in Fig. 1.

Dreesen et al. (2015) found that the underlying diagonal form of
$\mathcal{J}^{\prime}$ can be exploited in the decoupling process. It was
suggested to construct a third order tensor, $\mathcal{J}$, out of evaluations
of the Jacobian of the known function, $\bm{f}(\bm{p})$, and compute a
diagonal decomposition such that $\mathcal{J}\approx\mathcal{J}^{\prime}$. In
the filtered CPD approach, finite difference filters are introduced into the
decomposition. This allows for the Jocobian tensor to be decomposed into the
more convenient factors $\\{\bm{W},\bm{V},\bm{G}\\}$, where $\bm{G}$ stores
evaluations of the univariate functions $g_{i}$.

Figure 1: Centre: a collection of evaluations of the Jacobian of the decoupled
function (Eq. (4)), stacked in the third dimension. Left: the corresponding
third order tensor, Right: extracting the central diagonal plane reveals a
diagonal tensor decomposition.

The method of Decuyper et al. (2021b) can be summarised in three steps:

  1. 1.

Collect the Jacobian matrices of the known function, $\bm{J}$, and stack them
into a three-way array, i.e. the Jacobian tensor
$\mathcal{J}\in\mathbb{R}^{n\times m\times N}$.

  2. 2.

Factor $\mathcal{J}$ into $\\{\bm{W},\bm{V},\bm{G}\\}$ by computing a filtered
diagonal tensor decomposition (F-CPD).

  3. 3.

Retrieve the functions, $g_{i}$, by parametrising the nonparametric estimates
stored in $\bm{G}$.

At the core of the algorithm lies an alternating least-squares routine which
iteratively updates $\bm{W},\bm{V}$, and $\bm{G}$. Starting from a random
initialisation111The retrieved local minimum depends on the initialisation
point. In this work, however, no significant impact of the initialisation was
observed when tackling the benchmark problem of Section 3. the following norm
is progressively minimised

$\underset{\bm{W},\bm{V},\bm{G}}{\operatorname{arg~{}min}}~{}\lVert\mathcal{J}-\llbracket\bm{W},\bm{V},\mathcal{F}_{C}(\bm{V})\circ\bm{G}\rrbracket\rVert_{F}^{2},$
(7)

where the shorthand notation
$\mathcal{J}^{\prime}=\llbracket\bm{W},\bm{V},\mathcal{F}_{C}(\bm{V})\circ\bm{G}\rrbracket$
is used and ‘$\circ$’ defines a matrix-column product, illustrated by Eq. (9).
In this formulation, $\mathcal{F}_{C}\in\mathbb{R}^{N\times N\times r}$ stores
a collection of _finite difference filters_ , with the $i$th filter

$\bm{F}_{C_{i}}\coloneqq\mathcal{F}_{C_{[:,:,i]}}.$ (8)

A finite difference filter is a matrix which upon multiplication with a vector
of function evaluations, returns a finite difference approximation. Recall
that the $i$th column of $\bm{H}$ stores evaluations of the derivative of
$g_{i}$ with respect to $z_{i}=\bm{v}_{i}^{\text{T}}\bm{p}$, resulting in the
dependency $\mathcal{F}_{C}(\bm{V})$. We may then express $\bm{h}_{i}$, as the
finite difference of the $i$th column of $\bm{G}$,

$\bm{h}_{i}\coloneqq\bm{F}_{C_{i}}\bm{g}_{i}.$ (9)

The subscript $C$ denotes that a central differencing scheme is used.

Eq. (7) can be broken down into 3 separate objectives, i.e. Eq. (10), (11),
and (12). Denoting the matricisations of $\mathcal{J}$ along its rows, columns
and tubes by $\bm{J}_{(1)},\bm{J}_{(2)},\bm{J}_{(3)}$, respectively, and using
‘$\odot$’ to denote the Khatri-Rao product, an update of $\bm{W}$ is found
from

$\underset{\bm{W}}{\operatorname{arg~{}min}}~{}\|\bm{J}_{(1)}-\bm{W}((\mathcal{F}_{C}(\bm{V})\circ\bm{G})\odot\bm{V})^{\text{T}}\|_{F}^{2}$
(10)

Given that $\bm{W}$ appears linearly in the objective, a closed-form update
formula may be obtained.

What is imperative is to retrieve meaningful estimates of the functions
$g_{i}$ along the columns of $\bm{G}$. Originally one relied on the uniqueness
properties of the CPD to meet this requirement (Dreesen et al., 2015). As was
mentioned earlier, it has been shown that this uniqueness is often not
guaranteed, leading to very noisy estimate.

The F-CPD method, on the other hand, relies on finite difference filters to
steer the decomposition towards meaningful factors. What is required in
practice is to promote smoothness on the estimates stored in $\bm{G}$. Given
that both $\bm{V}$ (via $z_{i}=\bm{v}_{i}^{\text{T}}\bm{p}$) and $\bm{G}$
contribute to the smoothness of the estimates of $g_{i}(z_{i}$), a smoothness
objective should be reflected in both their cost functions. Regularisation is
used to penalise noisy estimates. An update of $\bm{V}$ is found from

$\begin{split}\underset{\bm{V}}{\operatorname{arg~{}min}}&~{}\|\bm{J}_{(2)}-\bm{V}((\mathcal{F}_{C}(\bm{V})\circ\bm{G})\odot\bm{W})^{\text{T}}\|_{F}^{2}\\\
&+\lambda\|\left(\mathcal{F}_{L}(\bm{V})\circ\bm{G}\right)-\left(\mathcal{F}_{R}(\bm{V})\circ\bm{G}\right)\|_{F}^{2}\end{split}$
(11)

where $\lambda$ is a hyperparameter which balances both objectives. The
additional term penalises divergent results from a left ($\mathcal{F}_{L}$)
and a right ($\mathcal{F}_{R}$) finite difference filtering operation,
ultimately steering the optimisation towards smooth solutions. Given that
$\bm{V}$ appears nonlinearly in Eq. (11), nonlinear optimisation is required
when computing an update.

In analogy with Eq. (11), the update formula of $\bm{G}$ is also found from a
joint objective function.

$\begin{split}\underset{\bm{G}}{\operatorname{arg~{}min}}&~{}\|\bm{J}_{(3)}-(\mathcal{F}_{C}(\bm{V})\circ\bm{G})(\bm{V}\odot\bm{W})^{\text{T}}\|_{F}^{2}\\\
&+\lambda\|\left(\mathcal{F}_{L}(\bm{V})\circ\bm{G}\right)-\left(\mathcal{F}_{R}(\bm{V})\circ\bm{G}\right)\|_{F}^{2}\end{split}$
(12)

It can be shown that $\bm{G}$ appears linearly in Eq. (12).

The smoothness objective ensures that the functions may be parameterised using
an appropriate basis expansion. What basis expansion to use depends on the
application and can be freely chosen by the user.

The computational cost of the algorithm is quadratic in both $N$ and $r$. More
efficient implementations are the subject of future study. For now, the
benchmark problem of Section 3 could be solved in a computing time in the
order of minutes. The present algorithm can be summarised by the following
pseudocode.

F-CPD Algorithm
---
Construct Jacobian tensor $\mathcal{J}$
Randomly select $N$ points from the input space
$\left\\{\bm{p}[k]\right\\}_{k=1}^{N}$
Compute $\bm{J}_{k}$ on the operating points
$\left\\{\bm{p}[k]\right\\}_{k=1}^{N}$
Stack the Jacobian matrices $\mathcal{J}[:,:,k]\coloneqq\bm{J}_{k}$
Factor $\mathcal{J}$ into $\left\\{\bm{W},\bm{V},\bm{G}\right\\}$
Initialise $\bm{W},\bm{V},\bm{G}$ for chosen value of $r$
Repeat ALS routine
update $\bm{W},\bm{V},\bm{G}$ via (10) to (12)
Until maximum number of iterations
Parameterise $\bm{G}$
obtain $g_{i}(z_{i})$ from appropriate basis expansion of $\bm{G}$

The filtered CPD approach is a powerful tool since it no longer relies on the
uniqueness properties of the CPD. As a result, the number of univariate
functions, $r$, to be used in the decoupled function, has become a design
choice.

Moreover, also multiple-input single-output functions, may be decoupled using
filtered CPD. MISO functions result in Jacobian matrices rather than tensors,
preventing the possibility of exploiting the uniqueness properties of the CPD.
F-CPD therefore enables the decoupling of single-output NARX models, and by
extent of general MISO functions.

##  3 Benchmark problem

The decoupling of single-output P-NARX models will be demonstrated on a
nonlinear benchmark data set of the forced duffing oscillator. The data is
obtained from an electrical implementation of a mechanically resonating system
involving a moving mass $m$, a viscous damping $c$ and a nonlinear spring
$k(y(t))$. The analogue electrical circuitry generates data close to but not
exactly equal to the idealised representation given by the nonlinear ordinary
differential equation (ODE)

$m\ddot{y}(t)+c\dot{y}(t)+k(y(t))y(t)=u(t),$ (13)

where the presumed displacement, $y(t)$, is considered the output and the
presumed force, $u(t)$, is considered the input. Overdots denote the
derivative with respect to time. The static position-dependent stiffness is
given by

$k(y(t))=\alpha+\beta y^{2}(t),$ (14)

which can be interpreted as a cubic hardening spring.

The training data consists of 9 realisations of a random-phase odd multisine.
The period of the multisine is ${}^{1}\\!/_{f_{0}}$ with
$f_{0}={}^{f_{s}}\\!/_{8192}$ Hz and $f_{s}\approx 610$ Hz. The number of
excited harmonics is $L=1342$ resulting in an $f_{\text{max}}\approx 200$ Hz.
Each multisine realisation is given a unique set of phases that are
independent and uniformly distributed in $[0,2\pi[$. The signal to noise ratio
(SNR) at the output is estimated at approximately 40 dB. This is measurement
noise (or sensor noise). It is high levels of measurement noise can
deteriorate the performance of NARX models (Schoukens et al., 2021). It is
important to stress that this inherent sensitivity is not removed by replacing
the internal function by a decoupled one.

As validation data a filtered Gaussian noise sequence of the same band width
and with a linearly increasing amplitude is used.

The data are part of three benchmark data sets for nonlinear system
identification described in Wigren and Schoukens (2013). The approach will
consist of 3 steps:

  1. 1.

Identify a P-NARX model on the basis of the training data.

  2. 2.

Decouple the NARX model using the F-CPD technique.

  3. 3.

Use a final optimisation to minimise the simulation error of the decoupled
model.

###  3.1 Reference P-NARX model

Using the System Identification toolbox in MATLAB a P-NARX model is estimated
with the following properties: $n_{u}=1$, $n_{y}=3$, and $d=3$. All cross-term
monomials where included leading to a model with 55 parameters. We will denote
the P-NARX model by
$f(\bm{x}):\mathbb{R}^{n_{u}+n_{y}+1}\rightarrow\mathbb{R}$, with $\bm{x}$ as
defined in Eq. (1), emphasising that the model is obtained by minimising the
equation error (focus on prediction). As performance metric a relative root-
mean-squared simulation error is used

$e_{\text{rms}}=\frac{\sqrt{\frac{1}{N}\sum_{k=1}^{N}\left(y[k]-y_{s}[k]\right)^{2}}}{\sqrt{\frac{1}{N}\sum_{k=1}^{N}y[k]^{2}}}\times
100$ (15)

where $y_{s}$ denotes the simulated output and $N$ is the record length. The
estimation process returns an accurate P-NARX model yielding a
$e_{\text{rms}}$ of 1.01% when simulating the validation data.

###  3.2 Decoupled P-NARX model

Given that using the F-CPD method, $r$, has become a design choice, we are
able to scan over $r$ and study the performance of the obtained decoupled
models. Additionally, a scan over the hyperparameter $\lambda$ is required.
Besides the simulation error (Eq. (15)) we will also introduce a function
approximation error,

$e_{f}=\frac{\sqrt{\frac{1}{N}\sum_{k=1}^{N}\left(f(\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{x}_{o}[k]\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0})-f_{d}(\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{x}_{o}[k]\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0})\right)^{2}}}{\sqrt{\frac{1}{N}\sum_{k=1}^{N}f(\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{x}_{o}[k]\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0})^{2}}}\times
100$ (16)

in which $f_{d}$ represents the decoupled polynomial function,

$f_{d}=\bm{w}^{\text{T}}\bm{g}\left(\bm{V}^{\text{T}}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{x}_{o}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\right),$
(17)

and $\bm{x}_{o}$ are the operating points. In this case $N$ refers to the
number of operating points for which the decomposition is computed.

In order to ensure that the operating points cover the region of interest,
they are selected on the basis of the training data. Simulating the training
data using the P-NARX estimate $f$ returns the collection
$\bm{X}_{T}=\\{\bm{x}_{T}[k]\\}_{k=0}^{N_{T}}$, with $\bm{x}_{T}$ containing
the simulated output samples

$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{x}_{T}=\\{u(t),u(t-1),\dots,u(t-n_{u}),\\\
\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}y_{s}(t-1),\dots,y_{s}(t-n_{y})\\},\quad\quad$
(18)

and $N_{T}$ denoting the training record length. The operating points are then
drawn from the joint normal distribution, inferred from $\bm{X}_{T}$. The
number of operating points to use can be considered a hyperparameter. In this
work $N=200$ points is used. In correspondence to the reference model, $f$,
the functions $g_{i}$ of the decoupled models will be parametrised using third
order polynomials, $d=3$.

We will study the decoupled models following the grid defined by $r=1,\dots,6$
and $\lambda=10^{l}$, with $l=-1,\dots,5$. For every value of $r$ the
decoupled model yielding the lowest function approximation error, $e_{f}$, is
selected. This will correspond to the value of $\lambda$ which results in the
appropriate balance between the tensor approximation objective and the
smoothness objective (Eq. (11) and (12)). The overview of the function
approximation error is depicted in Fig. 2. Accurate function approximations,
$f_{d}\approx f$, could be obtained yielding errors below 1% for values of
$r>1$.

Figure 2: Relative function approximation error computed on the operating
points.

In a final step, nonlinear optimisation is used to fine-tune the decoupled
models on the basis of their training data simulation error. A Levenberg-
Marquardt algorithm (Levenberg, 1944) is used to minimise the output error
objective given by

$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\underset{\bm{w},\bm{V},\bm{\theta}}{\operatorname{arg~{}min}}\sum_{k=1}^{N_{T}}\left(y[k]-\bm{w}^{\text{T}}\bm{g}(\bm{V}^{\text{T}}\bm{x}_{s}[k],\bm{\theta})\right)^{2},$
(19)

with $y$ the true output, $\bm{x}_{s}$ containing simulated output samples of
the decoupled model (similar to Eq. (18)), and $\bm{\theta}$ storing the
coefficients of the third order polynomials in $\bm{g}$.

Fig. 3 illustrates the performance of the optimised model set when simulating
the validation data. It is clear that decoupled models for which $r\geq 4$
perform equally well as the reference model, i.e $\approx 1$% (indicated by
the red line). Models containing less univariate mapping functions miss the
required complexity to reproduce the data up to such precision. For some
applications, however, the performance of the $r=2$ and $r=3$ model may still
be acceptable. The simulation performance of the $r=4$ model is illustrated in
Fig. 4.

Figure 3: Black markers: Relative simulation error of the decoupled P-NARX
models, computed on the validation data. Red: error obtained from the
reference P-NARX model. Figure 4: Blue: validation output data corresponding
to a filtered Gaussian noise sequence. The amplitude extends beyond the
training data amplitude resulting in an extrapolation of the model. Red:
simulation error of the $r=4$ decoupled model.

Noticing that $\bm{w}$ merely serves as a scaling on the functions $\bm{g}$,
the coefficients can be incorporated in $\bm{\theta}$. Doing so results in a
total parameter count of 36 for the $r=4$ model, compared to the 55 parameters
required in the standard basis expansion of the reference model. The
univariate functions of the decoupled model are depicted in Fig. 5. From
visual inspection of the functions one may conclude that the system behaviour
is dominantly cubic. Additionally the SISO functions can easily be monitored
to flag extrapolation (Karami et al., 2019). It should be noted that in this
case, even after decoupling, physical interpretability remains hard. An
important aspect is believed to be the choice of regressors, i.e.
$n_{u},n_{y}$. The use of shifted input samples contributes to retrieving a
non-physical model, in this case of the forced Duffing oscillator. One such
mechanism is the introduction of sampling zeros (Goodwin et al., 2013).

Figure 5: Univariate functions of the $r=4$ decoupled P-NARX model. The
functions are evaluated over the validation range.

##  4 Comparing polynomial constraints to F-CPD

In this section the decoupling technique of Karami et al. (2021), which will
be referred to as the _the structured Hessian method_ , is compared to the
F-CPD method. Following the author’s notation, the decoupled model may be
written as

$\displaystyle f_{d}=c_{0}+\sum_{i=1}^{r}g_{i}(\bm{v}_{i}^{\text{T}}\bm{x})$
(20a) $\displaystyle g_{i}(z_{i})=\sum_{j=1}^{d}c_{j,i}z_{i}^{j},$ (20b)

which corresponds to Eq. (17) where $\bm{w}$ is incorporated in the vector of
coefficients,

$\bm{c}=[c_{0}~{}c_{1,1}~{}\dots~{}c_{d,1}~{}c_{1,2}~{}\dots~{}c_{d,r}]^{\text{T}}.$
(21)

Computing the Hessian of the output of the decoupled function with respect to
the input variables generates a matrix for each operating point $\bm{x}_{o}$

$\bm{H}^{\prime}_{y}(\bm{x}_{o})=\sum_{i=1}^{r}g_{i}^{\prime\prime}(\bm{v}_{i}^{\text{T}}\bm{x}_{o})\bm{v}_{i}\bm{v}_{i}^{\text{T}}$
(22)

where $g_{i}^{\prime\prime}$ denotes the second derivative of $g_{i}$ with
respect to $z_{i}$. The collection of Hessian matrices may again be stacked
into a three-way array whose entries are given by

$\mathcal{H}^{\prime}_{[j,k,l]}=\sum_{i=1}^{r}v_{ij}~{}v_{ik}~{}g_{il}^{\prime\prime}$
(23)

$\text{for}~{}j=1,\dots,n~{}\text{and}~{}k=1,\dots,n~{}\text{and}~{}l=1,\dots,N$.

The idea is to construct the Hessian tensor $\mathcal{H}$ out of evaluations
of the known function and use the CPD to factor it into
$\mathcal{H}^{\prime}=\llbracket\bm{V},\bm{V},\bm{G}^{\prime\prime}\rrbracket$,
such that $\mathcal{H}\approx\mathcal{H}^{\prime}$. In this case
$\bm{G}^{\prime\prime}=[\bm{g}_{1}^{\prime\prime}\dots\bm{g}_{r}^{\prime\prime}]$
is a matrix storing evaluations of the second derivate of $\bm{g}$.

From the previous it is clear that relying on the uniqueness properties of the
CPD alone does not guarantee to obtain accurate estimates of
$g_{i}^{\prime\prime}$. It was therefore proposed not to solve for
$\bm{G}^{\prime\prime}$ directly, but to formulate the matrix factor into a
polynomial form and solve for the polynomial coefficients instead. The
optimisation problem to be solved can be formulated as

$\underset{\bm{V},\bm{c}}{\operatorname{arg~{}min}}\lVert\mathcal{H}-\llbracket\bm{V},\bm{V},\bm{G}^{\prime\prime}_{P}(\bm{c},\bm{V})\rrbracket\rVert_{F}^{2},$
(24)

where the columns of $\bm{G}^{\prime\prime}_{P}$ are expressed as polynomials
using the coefficients $\bm{c}$ and the Vandermonde matrices based on all the
$z_{i}$ (hence the dependence on $\bm{V}$ given
$z_{i}=\bm{v}_{i}^{\text{T}}\bm{x}$). The optimisation problem is solved using
a quasi-Newton algorithm.

The obtained estimate of $\bm{V}$ is then used to initialise a final
optimisation based on the simulated training output of the decoupled model

$\underset{\bm{V},\bm{c}}{\operatorname{arg~{}min}}\sum_{k=1}^{N_{T}}\left(y[k]-y_{s}[k](\bm{V},\bm{c})\right)^{2}.$
(25)

Starting from the reference model of Section 3.1, the structured Hessian
method is used to obtain a decoupled P-NARX model. The exact same data and
operating points are used for fair comparison.

A decoupled P-NARX model with $r=4$ functions is obtained. Also here
polynomials of the third degree are used. The obtained model is of an
identical architecture as the $r=4$ model obtained in Section 3.2, i.e. when
considering the $\bm{w}$ vector to be incorporated in $\bm{g}$. In Fig. 6, the
performance of the decoupled model is compared to the models obtained in
Section 3.2. The decoupled model has a slightly higher error of
$e_{\text{rms}}=1.78$%.

Figure 6: Black markers: Relative simulation error of the decoupled P-NARX
models, computed on the validation data. Red: error obtained from the
reference P-NARX model. Blue: error of the $r=4$ model obtained from the
structured Hessian method.

A number of advantages and disadvantages of both methods are listed:

  * •

The structured Hessian method requires the expensive computation of the
Hessian.

  * •

The F-CPD method requires a search over the hyperparameter $\lambda$.

  * •

The structured Hessian method enforces polynomial constraints on one of the
factors. This may be inappropriate when facing non-polynomial nonlinearities.

  * •

The structured Hessian method retrieves an estimate of the factor $\bm{V}$
from the constrained CPD of the Hessian. The Hessian, however, no longer
contains information on the linear part of the NARX model.

##  5 Conclusion

In this work a filtered CPD is used to decouple single-output polynomial NARX
models. Polynomial NARX models typically require a large number of parameters
in their description. For decoupled structures, the number of parameters grows
linearly with the degree, resulting in a substantial model reduction. It was
shown that the filtered CPD method no longer relies on the uniqueness
properties of the CPD. As a result, the number of univariate functions in the
decoupled structure has become a design choice. The method is illustrated on
the benchmark problem of the forced Duffing oscillator and compared to the
results obtained from the structured Hessian method.

## References

  * Billings (2013) Billings, S. (2013).  _Nonlinear System Identification: NARMAX Methods in the Time, Frequency and Spatio-Temporal Domains_.  Wiley.
  * Chan et al. (2015) Chan, R., Yuen, J., Lee, E., and Arashpour, M. (2015).  Application of nonlinear-autoregressive-exogenous model to predict the hysteretic behaviour of passive control systems.  _Eng. Struct._ , 85, 1–10.
  * Decuyper et al. (2019) Decuyper, J., Dreesen, P., Schoukens, J., Runacres, M., and Tiels, K. (2019).  Decoupling multivariate polynomials for nonlinear state-space models.  _IEEE Control Systems Letters_ , 3(3), 745–750.
  * Decuyper et al. (2021a) Decuyper, J., Tiels, K., Runacres, M.C., and Schoukens, J. (2021a).  Retrieving highly structured models starting from black-box nonlinear state-space models using polynomial decoupling.  _Mechanical Systems And Signal Processing_ , 146.
  * Decuyper et al. (2021b) Decuyper, J., Tiels, K., Weiland, S., and Schoukens, J. (2021b).  Decoupling multivariate functions using a non-parametric Filtered-CPD approach.  In _IFAC Symposium on System Identification: learning models for decision and control_. Padova, Italy.
  * Dreesen et al. (2015) Dreesen, P., Ishteva, M., and Schoukens, J. (2015).  Decoupling multivariate polynomials using first-order information.  _SIAM Journal on Matrix Analysis and Applications_ , 36(2), 864–879.
  * Goodwin et al. (2013) Goodwin, G.C., Aguero, J.C., Cea Garridos, M.E., Salgado, M.E., and Yuz, J.I. (2013).  Sampling and sampled-data models: The interface between the continuous world and digital algorithms.  _IEEE Control Systems Magazine_ , 33(5), 34–53.  10.1109/MCS.2013.2270403.
  * Karami et al. (2021) Karami, K., Westwick, D., and Schoukens, J. (2021).  Applying polynomial decoupling methods to the polynomial NARX model.  _Mechanical Systems And Signal Processing_ , 148.
  * Karami et al. (2019) Karami, K., Westwick, D., and Schoukens, J. (2019).  Identification of decoupled polynomial narx model using simulation error minimization* supported by the natural sciences and engineering research council of canada (nserc) through grant rgpin/06464-2015.  In _2019 American Control Conference (ACC)_ , 4362–4367.  10.23919/ACC.2019.8815172.
  * Kolda and Bader (2009) Kolda, T.G. and Bader, B.W. (2009).  Tensor decomposition and applications.  _SIAM Rev._ , 51(3), 455–500.
  * Levenberg (1944) Levenberg, K. (1944).  A method for the solution of certain non-linear problems in least squares.  _Quarterly of Applied Mathematics_ , 2(2), 164–168.
  * Schoukens et al. (2021) Schoukens, J., Westwick, D., Ljung, L., and Dobrowiecki, T. (2021).  Nonlinear system identification with dominating output noise - a case study on the Silverbox.  In _IFAC Symposium on System Identification: learning models for decision and control_.
  * Westwick et al. (2018) Westwick, D., Hollander, G., Karami, K., and Schoukens, J. (2018).  Using decoupling methods to reduce polynomial NARX models.  In _18th IFAC Symposium on System Identification_. Stockholm, Sweden.
  * Wigren and Schoukens (2013) Wigren, T. and Schoukens, J. (2013).  Three free data sets for development and benchmarking in nonlinear system identification.  In _European Control Conference (ECC)_ , 2933–2938. Zurich, Switzerland.
  * Zhao et al. (2013) Zhao, Y., Billings, S., Wei, H., He, F., and Sarrigiannis, P. (2013).  A new NARX-based Granger linear and nonlinear casual influence detection method with application to EEG data.  _J. Neurosci. Methods_ , 212(1), 79–86.