# An Online Algorithm for Contrastive Principal Component Analysis

###### Abstract

Finding informative low-dimensional representations that can be computed
efficiently in large datasets is an important problem in data analysis.
Recently, contrastive Principal Component Analysis (cPCA) was proposed as a
more informative generalization of PCA that takes advantage of contrastive
learning. However, the performance of cPCA is sensitive to hyper-parameter
choice and there is currently no online algorithm for implementing cPCA. Here,
we introduce a modified cPCA method, which we denote cPCA$\\!\ast$, that is
more interpretable and less sensitive to the choice of hyper-parameter. We
derive an online algorithm for cPCA$\\!\ast$ and show that it maps onto a
neural network with local learning rules, so it can potentially be implemented
in energy efficient neuromorphic hardware. We evaluate the performance of our
online algorithm on real datasets and highlight the differences and
similarities with the original formulation.

Index Terms—  Contrastive principal component analysis, online algorithm,
neural network, local learning rules

##  1 Introduction

Principal Component Analysis (PCA) is a long-standing pillar of dimensionality
reduction that is often a first step in data analysis. As a spectral method
without hyper-parameters, PCA is robust and interpretable. Furthermore, there
exist efficient algorithms for computing the principal subspace projections
[1], including online algorithms that map onto networks with local learning
rules and can potentially be implemented in energy efficient neuromorphic
hardware [2].

However, PCA projects a dataset onto the directions of highest variance,
without regard to the source of this variability. In particular, if there are
directions with high variability due to noise, PCA projects onto those
directions, potentially at the cost of more informative subspaces. Recently,
contrastive PCA (cPCA) was proposed as a generalization of PCA that can
account for such scenarios [3, 4]. By using a background dataset comprised of
domain-relevant ‘negative samples’ (e.g., a music studio recording sans
music), cPCA can distinguish directions of interest from those that have high
variance due to background noise. For example, in a drug trial, cPCA picks an
informative subspace by discounting the directions of high variance in the
control (or placebo) group. This method has proved useful for both analysis
and visualization in domains where a background dataset is available [5, 6,
7]. It has also been applied to the image and video compression domain [8].

While cPCA has been shown to find more informative projections than PCA [3],
this improvement has come at the cost of interpretability and efficiency. In
particular, cPCA is highly sensitive to a hyper-parameter that sets the
relative contribution of the positive and negative samples. This renders the
algorithm less robust and also less interpretable than PCA. Furthermore,
online PCA algorithms are not readily adapted to perform cPCA, so cPCA is less
amenable to large datasets or in scenarios with limited amount of compute.

In this work, we introduce cPCA$\\!\ast$, a more robust variation of cPCA. We
show cPCA$\\!\ast$ is less sensitive to hyper-parameter choice and, for
specific generative examples, can be interpreted as finding the subspace that
maximizes the signal-to-noise ratio. We derive an online algorithm for
cPCA$\\!\ast$ and map it onto a neural network with local learning rules which
can be implemented on efficient neuromorphic chips. In summary,

  * •

We introduce cPCA$\\!\ast$, an interpretable contrastive PCA method, and
derive an online algorithm that maps onto a neural network with local learning
rules.

  * •

We empirically show that cPCA$\\!\ast$ is more robust to the choice of hyper-
parameter than cPCA.

##  2 cPCA: original and our method

Given a dataset consisting of positive and negative samples, we want to
identify directions that have large variance in the positive samples but small
variance in the negative samples.

To be precise, let $d\geq 2$ and $({\bf x}_{1},\delta_{1}),\dots,({\bf
x}_{T},\delta_{T})\in\mathbb{R}^{d}\times\\{0,1\\}$ be a sequence of inputs.
The input ${\bf x}_{t}$ is a feature vector, which can either be a positive
sample or a negative sample. The scalar $\delta_{t}$ is an indicator variable
such that $\delta_{t}=1$ for positive samples and $\delta_{t}=0$ for negative
samples. We define covariance matrices of the positive and negative samples,
as follows:

$\displaystyle{\bf C}_{(+)}:=\langle{\bf x}_{t}{\bf
x}_{t}^{\top}|\delta_{t}=1\rangle_{t},$ $\displaystyle{\bf
C}_{(-)}:=\langle{\bf x}_{t}{\bf x}_{t}^{\top}|\delta_{t}=0\rangle_{t}.$

cPCA finds directions (i.e., unit vectors ${\bf v}\in\mathbb{R}^{d}$) that
simultaneously maximize ${\bf v}^{\top}{\bf C}_{(+)}{\bf v}$ and minimize
${\bf v}^{\top}{\bf C}_{(-)}{\bf v}$.

###  2.1 Contrastive PCA.

Abid et al. [3, 4] proposed projecting the positive samples onto the top
$k$-dimensional eigen-subspace of the matrix difference

$\displaystyle{\bf A}_{\alpha}:=(1-\alpha){\bf C}_{(+)}-\alpha{\bf
C}_{(-)},\qquad\alpha\in[0,1],$ (1)

where $1\leq k<d$ and $\alpha$ is the contrast parameter.111Abid et al. [4]
find the top eigen-subspace of ${\bf C}_{(+)}-\alpha{\bf C}_{(-)}$ for
$\alpha\geq 0$. We reparametrize the contrast parameter for a more direct
comparison with our method.

To gain intuition, consider the special case where positive samples correspond
to noisy signal and negative samples are pure noise; that is,

$\displaystyle{\bf x}_{t}=\delta_{t}\times\text{signal}+\text{noise}.$ (2)

If $\delta_{t}=1$, then ${\bf x}_{t}$ is signal plus noise, whereas if
$\delta_{t}=0$, then ${\bf x}_{t}$ is just noise. Under the assumption that
the signal and noise are decorrelated, the covariance of the signal is equal
to $2{\bf A}_{1/2}={\bf C}_{(+)}-{\bf C}_{(-)}$, so the directions that
maximize the variance of the signal are the top eigenvectors of ${\bf
A}_{1/2}$.

The contrast parameter $\alpha$ represents the trade-off between maximizing
the target variance and minimizing the background variance. When $\alpha=0$,
cPCA reduces to PCA applied to the positive samples. As $\alpha\to 1$,
directions that reduce the variance of the negative samples become more
optimal and the contrastive principal subspace is driven towards the minor
subspace of the negative samples. Therefore, each value of $\alpha$ yields a
direction with a different trade-off between positive and negative samples
variance.

In practice, the performance of cPCA is sensitive to the contrast parameter
$\alpha$ (see Fig. 2). While the generative model suggests that it should
perform optimally when $\alpha=1/2$, in practice, cPCA often performs better
for $\alpha>1/2$ (Fig. 2) and there is no theory for why this is true.
Furthermore, since ${\bf A}_{\alpha}$ is not a covariance matrix for
$\alpha\in(0,1)$, online PCA algorithms [1, 2] cannot be adapted to cPCA.

###  2.2 Our method: Contrastive PCA$\\!\ast$.

Our method, cPCA$\\!\ast$, is based on projecting the positive samples onto
the top $k$-dimensional eigen-subspace of the generalized eigenvalue problem

$\displaystyle{\bf C}_{(+)}{\bf v}=\lambda{\bf B}_{\beta}{\bf v},$ (3)

where

$\displaystyle{\bf B}_{\beta}:=(1-\beta){\bf I}_{d}+\beta{\bf
C}_{(-)},\qquad\beta\in[0,1],$ (4)

and the hyper-parameter $\beta$ is our contrast parameter.

To motivate cPCA$\\!\ast$, consider the generative model in Eq. (2), under the
assumption that the signal and noise are uncorrelated. Let ${\bf
C}_{\text{signal}}={\bf C}_{(+)}-{\bf C}_{(-)}$ and ${\bf
C}_{\text{noise}}={\bf B}_{1}={\bf C}_{(-)}$. Then the direction (i.e., unit
vector ${\bf v}$) that maximizes the signal-to-noise ratio is the unit vector
${\bf v}$ that maximizes the ratio

$\displaystyle\frac{{\bf v}^{\top}{\bf C}_{\text{signal}}{\bf v}}{{\bf
v}^{\top}{\bf C}_{\text{noise}}{\bf v}}=\frac{{\bf v}^{\top}{\bf C}_{(+)}{\bf
v}}{{\bf v}^{\top}{\bf B}_{1}{\bf v}}-1.$

The optimal ${\bf v}$ corresponds to the top eigenvector of the generalized
eigenvalue problem (3) with $\beta=1$.

As with the cPCA contrast parameter $\alpha$, when $\beta=0$, the problem
reduces to PCA applied to positive samples. As $\alpha\to 1$ in cPCA, the
projection subspace converges to the negative sample minor subspace, whereas
as $\beta\to 1$ in cPCA$\\!\ast$, the projection converges to maximize the
ratio of the variance of the positive sample projections and the variance of
the negative sample projections. For $\beta\in(0,1)$, our method is
interpretable as maximizing the ratio of the variance of the positive sample
projections and the variance of the negative sample projections when the
negative samples have been conditioned with independent isotropic noise with
variance $\tfrac{1-\beta}{\beta}$.

##  3 Online cPCA$\\!\ast$ algorithm

To derive an online algorithm, we concatenate the sequence of products
$\delta_{1}{\bf x}_{1},\dots,\delta_{T}{\bf x}_{T}\in\mathbb{R}^{d}$ into the
data matrix:

$\displaystyle{\bf X}_{(+)}:=[\delta_{1}{\bf x}_{1},\dots,\delta_{T}{\bf
x}_{T}]\in\mathbb{R}^{d\times T}.$

Since $\delta_{t}{\bf x}_{t}$ is zero for negative samples and equal to ${\bf
x}_{t}$ for positive samples, the covariance of ${\bf X}_{(+)}$ is a scalar
multiple of the covariance of the positive samples ${\bf C}_{(+)}$.

Similarity matching objective. Our starting point is the similarity matching
objective:

$\displaystyle\min_{{\bf Z}\in\mathbb{R}^{k\times T}}\frac{1}{T^{2}}\|{\bf
Z}^{\top}{\bf Z}-{\bf X}_{(+)}^{\top}{\bf B}_{\beta}^{-1}{\bf
X}_{(+)}\|_{F}^{2}.$ (5)

The similarity matching objective, which is a special case of a cost function
from multidimensional scaling [9], minimizes the difference in similarity
(measured using inner products) between the outputs ${\bf Z}$ and the positive
samples ${\bf X}_{(+)}$ normalized by ${\bf B}_{\beta}$. The optimal solution
${\bf Z}$ of equation (5) is the projection of ${\bf B}_{\beta}^{-1/2}{\bf
X}_{(+)}$ onto its $k$-dimensional principal subspace [9], which is equal to
the projection of ${\bf X}_{(+)}$ onto the top $k$-dimensional eigen-subspace
of the generalized eigenvalue problem (3).

Legendre transforms. Directly optimizing Eq. (5) does not result in an online
algorithm. Rather, we follow the approach in [10, 11] for deriving an online
algorithm. Expanding the square in Eq. (5) and dropping terms that do not
depend on ${\bf Z}$ yields the minimization problem

$\displaystyle\min_{{\bf Z}\in\mathbb{R}^{k\times
T}}\frac{1}{T^{2}}\operatorname{Tr}\left({\bf Z}^{\top}{\bf Z}{\bf
Z}^{\top}{\bf Z}-2{\bf Z}^{\top}{\bf Z}{\bf X}_{(+)}^{\top}{\bf
B}_{\beta}^{-1}{\bf X}_{(+)}\right).$

We substitute in with the Legendre transforms:

$\displaystyle\frac{1}{T^{2}}\operatorname{Tr}\left({\bf Z}^{\top}{\bf Z}{\bf
Z}^{\top}{\bf Z}\right)=\max_{{\bf
M}\in\mathbb{S}_{++}^{k}}\frac{2}{T}\operatorname{Tr}\left({\bf Z}^{\top}{\bf
M}{\bf Z}\right)-\operatorname{Tr}\left({\bf M}^{2}\right),$

where $\mathbb{S}_{++}^{k}$ is the set of $k\times k$ positive definite
matrices, and

$\displaystyle\frac{1}{T^{2}}\operatorname{Tr}\left({\bf Z}^{\top}{\bf Z}{\bf
X}_{(+)}^{\top}{\bf B}_{\beta}^{-1}{\bf X}_{(+)}\right)$
$\displaystyle\qquad=\max_{{\bf W}\in\mathbb{R}^{k\times
d}}\frac{2}{T}\operatorname{Tr}\left({\bf Z}^{\top}{\bf W}{\bf
X}_{(+)}\right)-\operatorname{Tr}\left({\bf W}{\bf B}_{\beta}{\bf
W}^{\top}\right),$

which are respectively optimized at ${\bf M}^{\ast}=\frac{1}{T}{\bf Z}{\bf
Z}^{\top}$ and ${\bf W}^{\ast}=\frac{1}{T}{\bf Z}{\bf X}_{(+)}^{\top}{\bf
B}_{\beta}^{-1}$. After substitution, we have

$\displaystyle\min_{{\bf Z}\in\mathbb{R}^{k\times T}}\min_{{\bf
W}\in\mathbb{R}^{k\times d}}\max_{{\bf M}\in\mathbb{S}_{++}^{k}}L({\bf W},{\bf
M},{\bf Z}),$

where

$\displaystyle L({\bf W},{\bf M},{\bf Z})$
$\displaystyle:=\frac{1}{T}\operatorname{Tr}(2{\bf Z}^{\top}{\bf M}{\bf
Z}-4{\bf Z}{\bf W}{\bf X}_{(+)})$ $\displaystyle\qquad-\operatorname{Tr}({\bf
M}^{2}+2{\bf W}{\bf B}_{\beta}{\bf W}^{\top}).$

Since the objective satisfies the saddle-point property with respect to ${\bf
M}$ and ${\bf Z}$, we can interchange the order of optimization:

$\displaystyle\min_{{\bf W}\in\mathbb{R}^{k\times d}}\max_{{\bf
M}\in\mathbb{S}_{++}^{k}}\min_{{\bf Z}\in\mathbb{R}^{k\times T}}L({\bf W},{\bf
M},{\bf Z}),$ (6)

Offline algorithm. Before deriving an online algorithm, it is instructive to
first solve Eq. (6) in the offline setting where we have access to the dataset
${\bf X}_{(+)}$ and the matrix ${\bf B}_{\beta}$. In this setting, we first
minimize $L({\bf W},{\bf M},{\bf Z})$ with respect to ${\bf Z}$:

$\displaystyle{\bf Z}={\bf M}^{-1}{\bf W}{\bf X}_{(+)}.$

We then take gradient descent-ascent steps:

$\displaystyle{\bf W}$ $\displaystyle\leftarrow{\bf
W}+2\eta\left(\frac{1}{T}{\bf Z}{\bf X}_{(+)}^{\top}-{\bf W}{\bf
B}_{\beta}\right)$ (7) $\displaystyle{\bf M}$ $\displaystyle\leftarrow{\bf
M}+\frac{\eta}{\tau}\left(\frac{1}{T}{\bf Z}{\bf Z}^{\top}-{\bf M}\right).$
(8)

Here $\tau>0$ denotes the ratio between the updates to ${\bf W}$ and ${\bf M}$
(which we set to 1 in our experiments) and $\eta\in(0,\tau)$ is the step size.

Online algorithm. To solve the minimax objective (6) in the online setting, we
take stochastic gradient descent-ascent steps. At each time step $t$, we first
minimize over the output ${\bf z}_{t}$ by running the continuous dynamics to
equilibrium:

$\displaystyle\dot{\bf z}_{t}(\gamma)=\delta_{t}{\bf c}_{t}-{\bf M}{\bf
z}_{t}(\gamma)\qquad\Rightarrow\qquad{\bf z}_{t}=\delta_{t}{\bf M}^{-1}{\bf
c}_{t},$ (9)

where we have defined ${\bf c}_{t}:={\bf W}{\bf x}_{t}$. Since $\delta_{t}=0$
when ${\bf x}_{t}$ is a negative sample, the algorithm only produces a non-
trivial output when presented with a positive sample. Next, we take a
stochastic gradient descent-ascent step with respect to $({\bf W},{\bf M})$.
We can replace the averages $\tfrac{1}{T}{\bf Z}{\bf X}_{(+)}^{\top}$ and
$\tfrac{1}{T}{\bf Z}{\bf Z}^{\top}$ in equations (7) and (8) with the rank-1
approximations ${\bf z}_{t}{\bf x}_{t}^{\top}$ and ${\bf z}_{t}{\bf
z}_{t}^{\top}$, respectively. To estimate ${\bf B}_{\beta}$ in the online
setting, we note that

$\displaystyle{\bf C}_{(-)}=\langle{\bf x}_{t}{\bf
x}_{t}^{\top}|\delta_{t}=0\rangle_{t}=\frac{\langle(1-\delta_{t}){\bf
x}_{t}{\bf x}_{t}^{\top}\rangle_{t}}{\langle 1-\delta_{t}\rangle_{t}}.$

Therefore, we can replace ${\bf B}_{\beta}$ with the online estimate

$\displaystyle\beta\frac{1-\delta_{t}}{p_{t}}{\bf x}_{t}{\bf
x}_{t}^{\top}+(1-\beta){\bf I}_{d},$

where $p_{t}$ denotes a running estimate of $\langle 1-\delta_{t}\rangle_{t}$
that is initialized at $p_{0}=0.5$ and updated at each time step, as follows:

$\displaystyle p_{t}=p_{t-1}+\frac{1}{t}(1-\delta_{t}-p_{t-1}).$

These substitutions result in the online updates

$\displaystyle{\bf W}$ $\displaystyle\leftarrow{\bf W}+2\eta\left({\bf
z}_{t}-\beta\frac{1-\delta_{t}}{p_{t}}{\bf c}_{t}\right){\bf
x}_{t}^{\top}-2\eta(1-\beta){\bf W}$ (10) $\displaystyle{\bf M}$
$\displaystyle\leftarrow{\bf M}+\frac{\eta}{\tau}\left({\bf z}_{t}{\bf
z}_{t}^{\top}-{\bf M}\right).$ (11)

We can map these onto a neural network, Fig. 1. Assuming that each neuron has
access to the running average $p_{t}$, the learning rules are local in the
sense that the update to a synaptic weight depends only on variables that are
available in the pre- and postsynaptic neurons.

Fig. 1: Neural network implementation of our online cPCA$\\!\ast$ algorithm.
At each $t$, the inputs $(x_{t}^{1},\dots,x_{t}^{d})$ are projected onto the
feedforward weights ${\bf W}$ to obtain $(c_{t}^{1},\dots,c_{t}^{k})$. Lateral
weights $-{\bf M}$ connect the output neurons. The neural dynamics in Eq. (9)
are run to equilibrium and the output of the network is
$(z_{t}^{1},\dots,z_{t}^{k})$. The synaptic weights ${\bf W}$ and $-{\bf M}$
are then updated according to Eqs. (10) and (11).

##  4 Numerical experiments

###  4.1 Comparison of cPCA and cPCA$\\!\ast$

We compare cPCA and cPCA$\\!\ast$ on two naturalistic datasets used in [3] as
well as a synthetic dataset. The results are summarized in Fig. 3.

Artificial dataset. The artificial dataset consists of 200 30-dimensional
samples. The negative samples are pure Gaussian noise while the positive
samples exhibit a bimodal distribution along some of the dimensions. Due to
high variance in the noisy directions, PCA is ineffective for identifying the
clusters.

Noisy Digits dataset. In the first naturalistic dataset, the positive samples
consist of 5,000 synthetic images generated by randomly superimposing images
of handwritten digits 0 and 1 from MNIST dataset [12] on top of natural images
taken from the UPenn Natural Image Database [13]. The negative samples consist
of natural images. We apply cPCA and cPCA$\\!\ast$ to find 2-dimensional
projections, for varying hyper-parameter choices, Figs. 2 and 3. To measure
how well-separated the projections of the noisy digits 0 and 1 are, we plot
the symmetrized KL divergence of the distributions of the projections of the
noisy 0’s and the distributions of the projections of the noisy 1’s. We see
that our method separates the digits for a much larger range of the hyper-
parameter (Fig. 2).

Mouse dataset. The second naturalistic dataset consists of protein expression
measurements $(d=77)$ of mice that have received shock therapy [14], some of
which develop Down Syndrome. We apply cPCA and cPCA$\\!\ast$ (with $k=2$)
using negative samples that consist of protein expression measurements from a
set of mice that have not been exposed to shock therapy and measure the degree
of linear separability between mice that do not have Down Syndrome and mice
that have Down Syndrome, Fig. 3.

Fig. 2: cPCA$\\!\ast$ is more robust than cPCA to the choice of hyper-
parameters. Separation of noisy MNIST digits 0 and 1 with (top) cPCA (top) vs.
(middle) cPCA$\\!\ast$. Bottom: Quantitative measure of separation of the
MNIST digits vs. value of hyper-parameter ($\alpha$ for cPCA and $\beta$ for
cPCA$\\!\ast$).  Fig. 3: Range of hyper-parameters ($\alpha$ for cPCA and
$\beta$ for cPCA$\\!\ast$) that lead to good performance (greater than 90%
classification accuracy using linear discriminant analysis) is broader in
cPCA$\\!\ast$.

###  4.2 Evaluation of online cPCA$\\!\ast$

Fig. 4 shows the convergence of the online cPCA$\\!\ast$ algorithm on the
mouse dataset. We see that cPCA$\\!\ast$ converges to the optimal projector.
Here, the learning rate is $\eta=0.003$ and the ratio between gradient descent
and ascent steps is $\tau=1$.

Fig. 4: Convergence of the online cPCA$\\!\ast$ algorithm on the mouse
dataset. The $y$-axis shows the alignment of the projector found via the
online algorithm to the optimal projector, computed by taking the trace of the
product of the two projectors. Shaded region shows standard deviation over 5
runs of the experiment.

##  5 Summary

In this work, we introduced cPCA$\\!\ast$, a modified contrastive Principal
Component Analysis method. We showed that this method is interpretable as
maximizing the signal to noise ratio and leads to an online algorithm which
can be mapped onto a neural network with local learning rules. In terms of the
scope of applicability, cPCA$\\!\ast$ has the same requirements as cPCA: as a
contrastive algorithm it needs a relevant background dataset. However, we hope
that the derivation of this online algorithm with more robust hyper-parameter
sensitivity will broaden the possible use-cases of this algorithm to larger
datasets.

## References

  * [1] Zeyuan Allen-Zhu and Yuanzhi Li,  “First efficient convergence for streaming k-PCA: a global, gap-free, and near-optimal rate,”  in 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS). IEEE, 2017, pp. 487–492.
  * [2] Cengiz Pehlevan, Tao Hu, and Dmitri B Chklovskii,  “A Hebbian/anti-Hebbian neural network for linear subspace learning: A derivation from multidimensional scaling of streaming data,”  Neural Computation, vol. 27, no. 7, pp. 1461–1495, 2015.
  * [3] Abubakar Abid, Vivek Kumar Bagaria, Martin Jinye Zhang, and James Y. Zou,  “Contrastive principal component analysis,”  arXiv preprint arXiv:1709.06716, 2017.
  * [4] Abubakar Abid, Martin J Zhang, Vivek K Bagaria, and James Zou,  “Exploring patterns enriched in a dataset with contrastive principal component analysis,”  Nature Communications, vol. 9, no. 1, pp. 1–7, 2018.
  * [5] Imrul Kaish, Jakir Hossain, Evangelos Papalexakis, and Jia Chen,  “COVID-19 or flu? Discriminative knowledge discovery of COVID-19 symptoms from Google Trends data,”  4th International Workshop on Epidemiology meets Data Mining and Knowledge discovery, 2021.
  * [6] Micol Marchetti-Bowick,  Structured Sparse Regression Methods for Learning from High-Dimensional Genomic Data,  Ph.D. thesis, Carnegie Mellon University, 2020.
  * [7] Takanori Fujiwara, Oh-Hyun Kwon, and Kwan-Liu Ma,  “Supporting analysis of dimensionality reduction results with contrastive learning,”  IEEE Transactions on Visualization and Computer Graphics, vol. 26, no. 1, pp. 45–55, 2020.
  * [8] Junfeng Xiao and Di Zhang,  “Full-reference image/video quality assessment algorithms based on contrastive principal component analysis,”  in 2022 7th International Conference on Image, Vision and Computing (ICIVC), 2022, pp. 648–653.
  * [9] Michael AA Cox and Trevor F Cox,  “Multidimensional scaling,”  in Handbook of Data Visualization, pp. 315–347. Springer, 2008\.
  * [10] David Lipshutz, Charles Windolf, Siavash Golkar, and Dmitri B Chklovskii,  “A biologically plausible neural network for slow feature analysis,”  Advances in Neural Information Processing Systems, vol. 33, pp. 14986–14996, 2020.
  * [11] David Lipshutz, Yanis Bahroun, Siavash Golkar, Anirvan M Sengupta, and Dmitri B Chklovskii,  “A biologically plausible neural network for multichannel canonical correlation analysis,”  Neural Computation, vol. 33, no. 9, pp. 2309–2352, 2021.
  * [12] Yann LeCun,  “The MNIST database of handwritten digits,”  http://yann.lecun.com/exdb/mnist/, 1998.
  * [13] Gašper Tkačik, Patrick Garrigan, Charles Ratliff, Grega Milčinski, Jennifer M Klein, Lucia H Seyfarth, Peter Sterling, David H Brainard, and Vijay Balasubramanian,  “Natural images from the birthplace of the human eye,”  PLOS One, vol. 6, no. 6, pp. e20409, 2011.
  * [14] Clara Higuera, Katheleen J Gardiner, and Krzysztof J Cios,  “Self-organizing feature maps identify proteins critical to learning in a mouse model of down syndrome,”  PLOS One, vol. 10, no. 6, pp. e0129126, 2015.


# Training Recurrent Neural Networks by Sequential Least Squares and the
Alternating Direction Method of Multipliers

Alberto Bemporad The author is with the IMT School for Advanced Studies,
Piazza San Francesco 19, Lucca, Italy. Email<EMAIL_ADDRESS>This paper was partially supported by the Italian Ministry of University and
Research under the PRIN’17 project “Data-driven learning of constrained
control systems”, contract no. 2017J89ARP.

###### Abstract

This paper proposes a novel algorithm for training recurrent neural network
models of nonlinear dynamical systems from an input/output training dataset.
Arbitrary convex and twice-differentiable loss functions and regularization
terms are handled by sequential least squares and either a line-search (LS) or
a trust-region method of Levenberg-Marquardt (LM) type for ensuring
convergence. In addition, to handle non-smooth regularization terms such as
$\ell_{1}$, $\ell_{0}$, and group-Lasso regularizers, as well as to impose
possibly non-convex constraints such as integer and mixed-integer constraints,
we combine sequential least squares with the alternating direction method of
multipliers (ADMM). We call the resulting algorithm NAILS (nonconvex ADMM
iterations and least squares) in the case line search (LS) is used, or NAILM
if a trust-region method (LM) is employed instead. The training method, which
is also applicable to feedforward neural networks as a special case, is tested
in three nonlinear system identification problems.

Keywords: Recurrent neural networks; nonlinear system identification;
nonlinear least-squares; generalized Gauss-Newton methods; Levenberg-Marquardt
algorithm; alternating direction method of multipliers; non-smooth loss
functions.

##  1 Introduction

In the last few years we have witnessed a tremendous attention towards machine
learning methods in several fields, including control. The use of neural
networks (NNs) for modeling dynamical systems, already popular in the nineties
[10, 28], is flourishing again, mainly due to the wide availability of
excellent open-source packages for automatic differentiation (AD). AD greatly
simplifies setting up and solving the training problem via numerical
optimization, in which the unknowns are the weights and bias terms of the NN.
While the attention has been mainly on the use of _feedforward_ NNs to model
the output function of a dynamical system in autoregressive form with
exogenous inputs (NNARX), _recurrent_ neural networks (RNNs) often better
capture the dynamics of the system, due to the fact that they are
intrinsically state-space models.

While minimizing the one-step-ahead output prediction error of the NNARX model
leads to a (nonconvex) unconstrained optimization problem with a separable
objective function, training RNNs is more difficult. This is because the
hidden states and parameters of the network are coupled by nonlinear equality
constraints. The training problem is usually solved by gradient descent
methods in which backpropagation through time (or its approximate truncated
version [40]) is used to evaluate derivatives. To enable the use of minibatch
stochastic gradient descent (SGD) methods, in [3] the authors propose to split
the experiment in multiple sections where the initial state of each section is
written as a NN function of a finite set of past inputs and outputs.

The main issue of gradient-descent methods is their well-known slow-rate of
convergence. To improve speed and quality of training, as well as to be able
to train RNN models in real-time as new input/output data become available,
training methods based on extended Kalman filtering (EKF) [25, 36, 15] or
recursive least squares (RLS) [41] were proposed based on minimizing the mean-
squared error (MSE) loss. These results were recently extended in [4] to
handle rather arbitrary convex and twice-differentiable loss terms when
training RNNs whose state-update and output functions are multi-layer NNs,
showing a drastic improvement in terms of computation efficiency and quality
of the trained model with respect to gradient-descent methods. For the special
case of shallow RNNs in which the states are known functions of the output
data and MSE loss, a combination of gradient descent and recursive least
squares was proposed in [24] to optimize the sequence of outputs and weights,
respectively.

In this paper we propose an offline training method based on two main
contributions. First, to handle arbitrary convex and twice-differentiable loss
terms, we propose the use of sequential least squares, constructed by
linearizing the RNN dynamics successively and taking a quadratic approximation
of the loss to minimize. An advantage of this approach, compared to more
classical backpropagation [38] is that the required Jacobian matrices, whose
evaluation is usually the most expensive part of the training algorithm, can
be computed in parallel, rather than sequentially, to propagate the gradients.
To force the objective function to decrease monotonically, we consider two
alternative ways: a line-search (LS) method and a trust-region method, the
latter in the classical Levenberg-Marquardt (LM) setting [13, 20, 33].

A second contribution of this paper is to combine sequential least squares
with nonconvex alternating direction method of multipliers (ADMM) iterations
to handle non-smooth and possibly non-convex regularization terms. ADMM has
been proposed in the past for training NNs, mainly feedforward NNs [31] but
also recurrent ones [30]. However, existing methods entirely rely on ADMM to
solve the learning problem, such as to gain efficiency by parallelizing
numerical operations. In this paper, ADMM is combined with the smooth
nonlinear optimizer based on sequential LS mentioned above to handle non-
smooth terms, therefore allowing handling sparsification (via $\ell_{1}$,
$\ell_{0}$, or group-Lasso penalties) and to impose (mixed-)integer
constraints, such as for quantizing the network coefficients. In particular,
we show that we can use ADMM and group Lasso to select the number of states to
include in the nonlinear RNN dynamics, therefore allowing a tradeoff between
model complexity and quality of fit, a fundamental aspect for the
identification of control-oriented black-box models. See also [14] for the use
of nonconvex ADMM to handle combinatorial constraints in combination with RNN
training.

We call the overall algorithm NAILS (nonconvex ADMM iterations and least
squares) when line search (LS) is used, or NAILM if a trust-region approach of
Levenberg-Marquardt (LM) type is employed instead. We show that NAILS/NAILM is
computationally more efficient, provides better-quality solutions than
classical gradient-descent methods, and is competitive with respect to EKF-
based approaches. The algorithm is also applicable to the simpler case of
training feedforward neural networks.

The paper is organized as follows. After formulating the training problem in
Section 2, we present the sequential least-squares approach to solve the
smooth version of the problem in Section 3 and extend the method to handle
non-smooth/non-convex regularization terms by ADMM iterations in Section 4.
The performance and versatility of the proposed approach is shown in three
nonlinear identification problems in Section 5.

##  2 Problem formulation

Given a set of input/output training data $(u_{0},y_{0})$, $\ldots$,
$(u_{N-1},y_{N-1})$, $u_{k}\in{\mathbb{R}}^{n_{u}}$,
$y_{k}\in{\mathbb{R}}^{n_{y}}$, we want to identify a nonlinear state-space
model in recurrent neural network (RNN) form

$x_{k+1}=f_{x}(x_{k},u_{k},\theta_{x}),\quad\hat{y}_{k}=f_{y}(x_{k},u_{k},\theta_{y})$
(1)

where $k$ denotes the sample instant, $x_{k}\in{\mathbb{R}}^{n_{x}}$ is the
vector of hidden states,
$f_{x}:{\mathbb{R}}^{n_{x}}\times{\mathbb{R}}^{n_{u}}\to{\mathbb{R}}^{n_{x}}$
and
$f_{y}:{\mathbb{R}}^{n_{x}}\times{\mathbb{R}}^{n_{u}}\to{\mathbb{R}}^{n_{y}}$
are feedforward neural networks parameterized by weight/bias terms
$\theta_{x}\in{\mathbb{R}}^{n_{\theta x}}$ and
$\theta_{y}\in{\mathbb{R}}^{n_{\theta y}}$, respectively, as follows [4]:

$\left\\{\begin{array}[]{rcl}v_{1k}^{x}&=&A_{1}^{x}\left[\begin{smallmatrix}x_{k}\\\
u_{k}\end{smallmatrix}\right]+b_{1}^{x}\\\
v_{2k}^{x}&=&A_{2}^{x}f_{1}^{x}(v_{1k}^{x})+b_{2}^{x}\\\ \vdots&&\vdots\\\
v_{L_{x}k}^{x}&=&A_{L_{x}}f_{L_{x}-1}^{x}(v_{(L_{x}-1)k}^{x})+b_{L_{x}}^{x}\\\
x_{k+1}&=&v_{L_{x}k}^{x}\end{array}\right.$ (2a)
$\left\\{\begin{array}[]{rcl}v_{1k}^{y}&=&A_{1}^{y}\left[\begin{smallmatrix}x_{k}\\\
u_{k}\end{smallmatrix}\right]+b_{1}^{y}\\\
v_{2k}^{y}&=&A_{2}^{y}f_{1}^{y}(v_{1k}^{y})+b_{2}^{y}\\\ \vdots&&\vdots\\\
v_{L_{y}k}^{y}&=&A_{L_{y}}^{y}f_{L_{y}-1}^{y}(v_{(L_{y}-1)k}^{y})+b_{L_{y}}^{y}\\\
\hat{y}_{k}&=&f_{L_{y}}^{y}(v_{L_{y}k}^{y})\end{array}\right.$ (2b)

where $L_{x}-1\geq 0$ and $L_{y}-1\geq 0$ are the number of hidden layers of
the state-update and output functions, respectively,
$v_{i}^{x}\in{\mathbb{R}}^{n_{i}^{x}}$, $i=1,\ldots L_{x}-1$,
$v_{L_{x}}^{x}\in{\mathbb{R}}^{n_{x}}$, and
$v_{i}^{y}\in{\mathbb{R}}^{n_{i}^{y}}$, $i=1,\ldots L_{y}$, the values of the
corresponding inner layers,
$f_{i}^{x}:{\mathbb{R}}^{n_{i}^{x}}\to{\mathbb{R}}^{n_{i+1}^{x}}$, $i=1,\ldots
L_{x}-1$, $f_{i}^{y}:{\mathbb{R}}^{n_{i}^{y}}\to{\mathbb{R}}^{n_{i+1}^{y}}$,
$i=1,\ldots L_{y}-1$ the corresponding activation functions, and
$f_{L_{y}}^{y}:{\mathbb{R}}^{n_{L_{y}}^{y}}\to{\mathbb{R}}^{n_{y}}$ the output
function, e.g., $f_{L_{y}}^{y}(v_{L_{y}}^{y})=v_{L_{y}}^{y}$ to model
numerical outputs, or
$[f_{L_{y}}^{y}(v_{L_{y}}^{y})]_{i}=1/\left(1+e^{[v_{L_{y}}^{y}]_{i}}\right)$,
$i=1,\ldots,n_{y}$, for binary outputs. Note that the RNN (1) can be made
strictly causal by simply setting $v_{1}^{y}(k)=A_{1}^{y}x_{k}+b_{1}^{y}$ in
(2b).

The training problem we want to solve is:

$\displaystyle\min_{z}$ $\displaystyle
r(x_{0},\theta_{x},\theta_{y})+g(\theta_{x},\theta_{y})+\displaystyle{\sum_{k=0}^{N-1}\ell(y_{k},f_{y}(x_{k},u_{k},\theta_{y}))}$
s.t. $\displaystyle x_{k+1}=f_{x}(x_{k},u_{k},\theta_{x}),\ k=0,\ldots,N-2$
(3b)

where $z=[x_{0}^{\prime}\ \ldots\ x_{N-1}^{\prime}\ \theta_{x}^{\prime}\
\theta_{y}^{\prime}]^{\prime}\in{\mathbb{R}}^{n}$ is the overall optimization
vector, $n=Nn_{x}+n_{\theta_{x}}+n_{\theta_{y}}$,
$\ell:{\mathbb{R}}^{n_{y}}\times{\mathbb{R}}^{n_{y}}\to{\mathbb{R}}$ is a loss
function that we assume strongly convex and twice differentiable with respect
to its second argument,
$r:{\mathbb{R}}^{n_{x}}\times{\mathbb{R}}^{n_{\theta_{x}}}\times{\mathbb{R}}^{n_{\theta_{y}}}\to{\mathbb{R}}$
is a regularization term (such as an $\ell_{2}$-regularizer) that we also
assume strongly convex and twice differentiable, and
$g:{\mathbb{R}}^{n_{\theta_{x}}}\times{\mathbb{R}}^{n_{\theta_{x}}}\to{\mathbb{R}}\cup\\{+\infty\\}$
is a possibly non-smooth and non-convex regularization term. For simplicity,
in the following we assume that $r$ is separable with respect to $x_{0}$,
$\theta_{x}$, and $\theta_{y}$, i.e.,
$r(x_{0},\theta_{x},\theta_{y})=r_{x}(x_{0})+r_{\theta}^{x}(\theta_{x})+r_{\theta}^{y}(\theta_{y})$,
and that in case of multiple outputs $n_{y}>1$, the loss function $\ell$ is
also separable, i.e.,
$\ell(y,\hat{y})=\sum_{i=1}^{n_{y}}\ell_{i}(y_{ki},f_{yi}(x_{k},u_{k},\theta_{y}))$,
where the subscript $i$ denotes the $i$th component of the output signal
$y_{k}$ and the corresponding loss function.

Problem (3) can be rewritten as the following unconstrained nonlinear
programming problem in _condensed form_

$\min_{x_{0},\theta_{x},\theta_{y}}V(x_{0},\theta_{x},\theta_{y})+g(\theta_{x},\theta_{y})$
(4)

where we get

$V(x_{0},\theta_{x},\theta_{y})=r(x_{0},\theta_{x},\theta_{y})+\sum_{k=0}^{N-1}\ell(y_{k},f_{y}(x_{k},u_{k},\theta_{y}))$
(5)

by replacing $x_{k+1}=f_{x}(x_{k},u_{k},\theta_{x})$ iteratively. Note that
Problem (4) can be highly nonconvex and present several local minima, see,
e.g., recent studies reported in [27].

###  2.1 Multiple training traces and initial-state encoder

For simplicity of notation, in this paper we focus on training a RNN model
(2b) based on a single input/output trace. The extension to $M$ multiple
experiments is trivial and can be achieved by introducing an initial state
$x_{0,j}$ per trace and minimizing the sum of the corresponding loss functions
and the regularization terms with respect to
$(x_{0,1},\ldots,x_{0,M},\theta_{x},\theta_{y})$. In alternative, we can
parameterize

$x_{0,i}=f_{x0}(v_{0,i},\theta_{x0})$ (6)

where $v_{0}\in{\mathbb{R}}^{n_{v}}$ is a measured vector available at time 0,
such as a collection of $n_{a}$ past outputs and $n_{b}$ past inputs as
suggested in [21, 3], and/or other measurable values that are known to affect
the initial state of the system we want to model. The initial-state encoder
function
$f_{x0}:{\mathbb{R}}^{n_{v}}\times{\mathbb{R}}^{n_{\theta_{x0}}}\to{\mathbb{R}}^{n_{x}}$
is a feedforward neural network defined as in (2b), parameterized by the
vector $\theta_{x0}\in{\mathbb{R}}^{n_{\theta_{x0}}}$ to be learned jointly
with $\theta_{x}$, $\theta_{y}$. The regularization term in this case is
defined as $r(\theta_{x0},\theta_{x},\theta_{y})$.

##  3 Sequential linear least squares

We first handle the smooth case $g(\theta_{x},\theta_{y})=0$. Let
$(x_{0}^{0},\theta_{x}^{0},\theta_{y}^{0})$ be an initial guess. By simulating
the RNN model (1) we get the corresponding sequence of hidden states
$x_{k}^{0}$, $k=1,\ldots,N-1$, and predicted outputs $\hat{y}_{k}^{0}$. At a
generic iteration $h$ of the algorithm proposed next, we assume that the
nominal trajectory is feasible, i.e.,
$x_{k+1}^{h}=f_{x}(x_{k}^{h},u_{k},\theta_{x}^{h})$, $k=0,\ldots,N-2$. We want
to find updates $x_{k}^{h+1}$, $\theta_{x}^{h+1}$, $\theta_{y}^{h+1}$ by
solving a least-squares approximation of problem (3) as described in the next
paragraphs.

###  3.1 Linearization of the dynamic constraints

To get a quadratic approximation of problem (4) we first linearize the RNN
dynamic constraints (3b) around $x_{k}^{h}$ by approximating
$x_{k+1}^{h}+p_{x_{k+1}}$ $=$
$f_{x}(x_{k}^{h}+p_{x_{k}},u_{k},\theta_{x}^{h}+p_{\theta_{x}})$ $\approx$
$f_{x}(x_{k}^{h},u_{k},\theta_{x}^{h})$ $+$
$\nabla_{x}f_{x}(x_{k}^{h},u_{k},\theta_{x}^{h})p_{x_{k}}$ $+$
$\nabla_{\theta_{x}}f_{x}(x_{k}^{h},u_{k},\theta_{x}^{h})p_{\theta_{x}}$.
Since $x_{k+1}^{h}=f_{x}(x_{k}^{h},u_{k},\theta_{x}^{h})$, we can write the
linearized state-update equations as

$p_{x_{k+1}}=\nabla_{x}f_{x}(x_{k}^{h},u_{k},\theta_{x}^{h})p_{x_{k}}+\nabla_{\theta_{x}}f_{x}(x_{k}^{h},u_{k},\theta_{x}^{h})p_{\theta_{x}}$
(7)

The dynamics (7) can be condensed to

$p_{x_{k}}=M_{kx}p_{x_{0}}+M_{k\theta_{x}}p_{\theta_{x}}$ (8)

where matrices $M_{kx}\in{\mathbb{R}}^{n_{x}\times n_{x}}$ and
$M_{k\theta_{x}}\in{\mathbb{R}}^{n_{x}\times n_{\theta_{x}}}$,
$k=0,\ldots,N-2$, are recursively defined as follows:

$\displaystyle M_{0x}$ $\displaystyle=$ $\displaystyle I,\quad
M_{0\theta_{x}}=0$ (9) $\displaystyle M_{(k+1)x}$ $\displaystyle=$
$\displaystyle\nabla_{x}f_{x}(x_{k}^{h},u_{k},\theta_{x}^{h})M_{kx}$
$\displaystyle M_{(k+1)\theta_{x}}$ $\displaystyle=$
$\displaystyle\nabla_{x}f_{x}(x_{k}^{h},u_{k},\theta_{x}^{h})M_{k\theta_{x}}+\nabla_{\theta_{x}}f_{x}(x_{k}^{h},u_{k},\theta_{x}^{h})$

In case $x_{0}$ is parameterized as in (6), by approximating

$x_{0}^{h}\approx
f_{x0}(v_{0},\theta_{x0}^{h})+\nabla_{\theta_{x0}}f_{x0}(v_{0},\theta_{x0}^{h})p_{\theta_{x0}}$

where $p_{\theta_{x0}}$ is the increment of $\theta_{x0}$ to be computed, it
is enough to replace $p_{x_{0}}$ with $p_{\theta_{x0}}$ in (8) and
$M_{0x}=\nabla_{\theta_{x0}}f_{x0}(v_{0},\theta_{x0}^{h})$ in (9).

###  3.2 Quadratic approximation of the cost function

Let us now find a quadratic approximation of the cost function (5). First,
regarding the regularization term $r$, we take the 2nd-order Taylor expansions

$\displaystyle r_{x}(x_{0}+p_{x_{0}})$ $\displaystyle\approx$
$\displaystyle\frac{1}{2}p_{x_{0}}^{\prime}\nabla^{2}r_{x}(x_{0}^{h})p_{x_{0}}+\nabla
r_{x}(x_{0}^{h})p_{x_{0}}+r(x_{0}^{h})$ $\displaystyle
r_{\theta}^{x}(\theta_{x}^{h}+p_{\theta_{x}})$ $\displaystyle\approx$
$\displaystyle\frac{1}{2}p_{\theta_{x}}^{\prime}\nabla^{2}r_{\theta}^{x}(\theta_{x}^{h})p_{\theta_{x}}+\nabla
r_{\theta}^{x}(\theta_{x}^{h})p_{\theta_{x}}+r_{\theta}^{x}(\theta_{x}^{h})$
$\displaystyle r_{\theta}^{y}(\theta_{y}^{h}+p_{\theta_{y}})$
$\displaystyle\approx$
$\displaystyle\frac{1}{2}p_{\theta_{y}}^{\prime}\nabla^{2}r_{\theta}^{y}(\theta_{y}^{h})p_{\theta_{y}}+\nabla
r_{\theta}^{y}(\theta_{y}^{h})p_{\theta_{y}}+r_{\theta}^{y}(\theta_{y}^{h})$

By neglecting the constant terms, minimizing the 2nd-order Taylor expansion of
$r(x_{0}+p_{x_{0}},\theta_{x}+p_{\theta_{x}},\theta_{y}+p_{\theta_{y}})$ is
equivalent to minimizing the least-squares term

$\frac{1}{2}\left\|\left[\begin{smallmatrix}L_{r0}p_{x0}+(L_{r0}^{\prime})^{-1}\nabla_{x}r_{x}(x_{0}^{h})\\\
L_{rx}p_{\theta_{x}}+(L_{rx}^{\prime})^{-1}\nabla
r_{\theta}^{x}(\theta_{x}^{h})\\\
L_{ry}p_{\theta_{y}}+(L_{ry}^{\prime})^{-1}\nabla
r_{\theta}^{y}(\theta_{y}^{h})\end{smallmatrix}\right]\right\|_{2}^{2}$ (10)

where $L_{r0}^{\prime}L_{r0}=\nabla^{2}r_{x}(x_{0})$,
$L_{rx}^{\prime}L_{rx}=\nabla^{2}r_{\theta}^{x}(\theta_{x})$, and
$L_{ry}^{\prime}L_{ry}=\nabla^{2}r_{\theta}^{y}(\theta_{y})$ are Cholesky
factorizations of the corresponding Hessian matrices. Note that in the case of
standard $\ell_{2}$-regularization
$r_{x}(x_{0})=\frac{\rho_{x}}{2}\|x_{0}\|_{2}^{2}$,
$r_{\theta}^{x}(\theta_{x})=\frac{\rho_{\theta}}{2}\|\theta_{x}\|_{2}^{2}$,
$r_{\theta}^{y}(\theta_{y})=\frac{\rho_{\theta}}{2}\|\theta_{y}\|_{2}^{2}$, we
simply have $L_{r0}=\sqrt{\rho_{x}}I$, $L_{rx}=\sqrt{\rho_{\theta}}I$,
$L_{ry}=\sqrt{\rho_{\theta}}I$, $\nabla_{x}r_{x}(x_{0})=\sqrt{\rho_{x}}x_{0}$,
$r_{\theta}^{x}(\theta_{x})=\sqrt{\rho_{\theta}}\theta_{x}$,
$r_{\theta}^{y}(\theta_{y})=\sqrt{\rho_{\theta}}\theta_{y}$, and hence (10)
becomes

$\frac{1}{2}\left\|\left[\begin{smallmatrix}\sqrt{\rho_{x}}p_{x0}+\sqrt{\rho_{x}}x_{0}^{h}\\\
\sqrt{\rho_{\theta}}p_{\theta_{x}}+\sqrt{\rho_{\theta}}\theta_{x}^{h}\\\
\sqrt{\rho_{\theta}}p_{\theta_{y}}+\sqrt{\rho_{\theta}}\theta_{y}^{h}\end{smallmatrix}\right]\right\|_{2}^{2}$
(11)

Regarding the loss terms penalizing the $i$th output-prediction error at step
$k$, we have that

$\displaystyle\nabla_{x_{k}}\ell_{i}(y_{ki},f_{yi}(x_{k}^{h},u_{k},\theta_{y}^{h}))$
$\displaystyle=$
$\displaystyle\nabla_{x_{k}}f_{yi}(x_{k}^{h},u_{k},\theta_{y}^{h})\ell^{\prime}_{ik}$
$\displaystyle\nabla_{\theta_{y}}\ell_{i}(y_{ki},f_{yi}(x_{k}^{h},u_{k},\theta_{y}^{h}))$
$\displaystyle=$
$\displaystyle\nabla_{\theta_{y}}f_{yi}(x_{k}^{h},u_{k},\theta_{y}^{h})\ell^{\prime}_{ik}$

where
$\ell^{\prime}_{ik}\triangleq\frac{d\ell_{i}(y_{k},f_{yi}(x_{k}^{h},u_{k},\theta_{y}^{h}))}{d\hat{y}_{i}}$,
$\ell^{\prime}_{ik}\in{\mathbb{R}}$, and hence, by differentiating again with
respect to $\left[\begin{smallmatrix}x_{k}\\\
\theta_{y}\end{smallmatrix}\right]$ and setting
$\ell^{\prime\prime}_{ik}\triangleq\frac{d^{2}\ell_{i}(y_{k},f_{yi}(x_{k}^{h},u_{k},\theta_{y}^{h})}{d\hat{y}_{i}^{2}}$,
$\ell^{\prime\prime}_{ik}\in{\mathbb{R}}$, $\ell^{\prime\prime}_{ik}>0$, we
get

$\nabla^{2}_{\left[\begin{smallmatrix}x_{k}\\\
\theta_{y}\end{smallmatrix}\right]\left[\begin{smallmatrix}x_{k}\\\
\theta_{y}\end{smallmatrix}\right]}\ell_{i}(y_{ki},f_{yi}(x_{k}^{h},u_{k},\theta_{y}^{h}))=(\nabla_{\left[\begin{smallmatrix}x_{k}\\\
\theta_{y}\end{smallmatrix}\right]}f_{yi})(\nabla_{\left[\begin{smallmatrix}x_{k}\\\
\theta_{y}\end{smallmatrix}\right]}f_{yi})^{\prime}\ell^{\prime\prime}_{ik}+\nabla^{2}_{\left[\begin{smallmatrix}x_{k}\\\
\theta_{y}\end{smallmatrix}\right]\left[\begin{smallmatrix}x_{k}\\\
\theta_{y}\end{smallmatrix}\right]}f_{yi}(x_{k}^{h},u_{k},\theta_{y}^{h})\ell^{\prime}_{ik}$

(we have omitted the dependence of $\ell^{\prime}$, $\ell^{\prime\prime}$ on
the iteration $h$ for simplicity). By neglecting the Hessian of the output
function $f_{yi}$ (as it is common in Gauss-Newton methods), the minimization
of
$\ell_{i}(y_{ki},f_{yi}(x_{k}^{h}+p_{x_{k}},u_{k},\theta_{x}^{h}+p_{\theta_{x}}))$
can be approximated by the minimization of

$\frac{1}{2}\left[\begin{smallmatrix}p_{x_{k}}\\\
p_{\theta_{x}}\end{smallmatrix}\right]^{\prime}\ell^{\prime\prime}_{ik}(\nabla_{\left[\begin{smallmatrix}x_{k}\\\
\theta_{y}\end{smallmatrix}\right]}f_{yi})(\nabla_{\left[\begin{smallmatrix}x_{k}\\\
\theta_{y}\end{smallmatrix}\right]}f_{yi})^{\prime}\left[\begin{smallmatrix}p_{x_{k}}\\\
p_{\theta_{y}}\end{smallmatrix}\right]+\ell^{\prime}_{ik}\left[\begin{smallmatrix}\nabla_{x_{k}}f_{yi}\\\
\nabla_{\theta_{y}}f_{yi}\end{smallmatrix}\right]^{\prime}\left[\begin{smallmatrix}p_{x_{k}}\\\
p_{\theta_{y}}\end{smallmatrix}\right]$

By strong convexity of the loss function $\ell$, this can be further rewritten
as the minimization of the least-squares term

$\frac{1}{2}\left\|(\ell^{\prime\prime}_{ik})^{\frac{1}{2}}(\nabla_{\left[\begin{smallmatrix}x_{k}\\\
\theta_{y}\end{smallmatrix}\right]}f_{yi})^{\prime}\left[\begin{smallmatrix}p_{x_{k}}\\\
p_{\theta_{y}}\end{smallmatrix}\right]+(\ell^{\prime\prime}_{ik})^{-\frac{1}{2}}\ell^{\prime}_{ik}\right\|_{2}^{2}$
(12)

that, by exploiting (8), is equivalent to

$\displaystyle\frac{1}{2}\left\|(\ell^{\prime\prime}_{ik})^{\frac{1}{2}}\left[\begin{array}[]{ccc}(\nabla_{x_{k}}f_{yi})^{\prime}M_{kx}&(\nabla_{x_{k}}f_{yi})^{\prime}M_{k\theta_{x}}&(\nabla_{\theta_{y}}f_{yi})^{\prime}\end{array}\right]\left[\begin{smallmatrix}p_{x_{k}}\\\
p_{\theta_{x}}\\\
p_{\theta_{y}}\end{smallmatrix}\right]+(\ell^{\prime\prime}_{ik})^{-\frac{1}{2}}\ell^{\prime}_{ik}\right\|_{2}^{2}$
(14)

Finally, the minimization of the quadratic approximation of (5) can be recast
as the linear LS problem

$p^{h}=\arg\min_{p}\frac{1}{2}\|A^{h}p-b^{h}\|_{2}^{2}$ (15)

where $p=\left[\begin{smallmatrix}p_{x_{0}}\\\ p_{\theta_{x}}\\\
p_{\theta_{y}}\end{smallmatrix}\right]$, $A^{h}$ and $b^{h}$ are constructed
by collecting the least-squares terms from (11) and (14),
$A^{h}\in{\mathbb{R}}^{n_{A}\times n_{p}}$, $b^{h}\in{\mathbb{R}}^{n_{A}}$,
$n_{A}=Nn_{y}+n_{x}+n_{\theta_{x}}+n_{\theta_{y}}$, and
$n_{p}=n_{x}+n_{\theta_{x}}+n_{\theta_{y}}$.

We remark that most of the computation effort is usually spent in computing
the Jacobian matrices $\nabla_{\left[\begin{smallmatrix}x_{k}\\\
\theta_{x}\end{smallmatrix}\right]}f_{x}(x_{k}^{h},u_{k},\theta_{x}^{h})$ and
$\nabla_{\left[\begin{smallmatrix}x_{k}\\\
\theta_{y}\end{smallmatrix}\right]}f_{yi}(x_{k}^{h},u_{k},\theta_{y}^{h})$.
Given the current nominal trajectory $\\{x_{k}^{h}\\}$, such computations can
be completely _parallelized_. Hence, we are not forced to use _sequential_
computations as in backpropagation [38] to solve for $p^{h}$ in (15), as
typically done in classical training methods based on gradient-descent.

Note also that for training _feedforward_ neural networks we set $n_{x}=0$,
$v_{1}^{y}(k)=A_{1}^{y}u_{k}+b_{1}^{y}$ in (2b), and only optimize with
respect to $\theta_{y}$. In this case, the operations performed at each
iteration $h$ to predict $\hat{y}_{k}$ and compute $V^{h}$ can be also
parallelized.

For the sake of comparison, a gradient-descent method would set

$\displaystyle p^{h}$ $\displaystyle=$ $\displaystyle-\nabla
V(x_{0}^{h},\theta_{x}^{h},\theta_{y}^{h})$ $\displaystyle=$
$\displaystyle\left[\begin{smallmatrix}\nabla_{x_{0}}r_{x}(x_{0}^{h})\\\
\nabla_{\theta_{x}}r_{\theta}^{x}(\theta_{x}^{h})\\\
\nabla_{\theta_{y}}r_{\theta}^{y}(\theta_{y}^{h})\end{smallmatrix}\right]+\sum_{k=0}^{N-1}\sum_{i=1}^{n_{y}}\ell^{\prime}_{ik}\left[\begin{smallmatrix}\nabla_{x_{k}}f_{yi}(x_{k}^{h},\theta_{y}^{h})^{\prime}M_{kx}\\\
\nabla_{x_{k}}f_{yi}(x_{k}^{h},\theta_{y}^{h})^{\prime}M_{k\theta_{x}}\\\
\nabla_{\theta_{y}}f_{yi}(x_{k}^{h},\theta_{y}^{h})\end{smallmatrix}\right]$

since, from (9), $\nabla_{x_{0}}x_{k}=M_{kx}$ and
$\nabla_{\theta_{x}}x_{k}=M_{k\theta_{x}}$.

####  3.2.1 Sequential unconstrained quadratic programming

The formulation (15) requires storing matrix $A^{h}$ and vector $b^{h}$. In
order to decrease memory requirements, especially in the case of a large
number $N$ of training data, an alternative method is to incrementally
construct the matrices
$H^{h}=(A^{h})^{\prime}A^{h}\in{\mathbb{R}}^{n_{p}\times n_{p}}$ and
$c^{h}=-(A^{h})^{\prime}b^{h}\in{\mathbb{R}}^{n_{p}}$ of the normal equation
$H^{h}p^{h}+c^{h}=0$ associated with (15). Although this approach might be
computationally more expensive, it only requires storing
$\frac{n_{p}^{2}+3n_{p}}{2}$ numbers. A further advantage of solving the
normal equation is that the convexity of $\ell$ with respect to $\hat{y}$ and
of $r(x,\theta_{x},\theta_{y})$ do not have to be strong, as instead required
to express the quadratic approximation of $f$ as a squared norm from (14) and
(10), respectively, as long as the computed $H^{h}$ is nonsingular.

####  3.2.2 Recursive least-squares

A further alternative to precomputing and storing the matrices $A^{h},b^{h}$
is to exploit (14) and use the following recursive least-square (RLS)
formulation of (15)

$\displaystyle P_{k}$ $\displaystyle=$ $\displaystyle
P_{k-1}-\frac{P_{k-1}\phi_{k}\phi_{k}^{\prime}P_{k-1}}{1+\phi_{k}^{\prime}P_{k-1}\phi_{k}}$
$\displaystyle p_{k}$ $\displaystyle=$ $\displaystyle
p_{k-1}+P_{k}\phi_{k}(\eta_{k}-p_{k-1}^{\prime}\phi_{k})$ (17)

for $k=0,\ldots,N-1$ where
$\phi_{k}=\sqrt{\ell^{\prime\prime}_{ik}}\left[\begin{smallmatrix}M^{\prime}_{kx}\nabla_{x_{k}}f_{yi}\\\
M^{\prime}_{k\theta_{x}}\nabla_{x_{k}}f_{yi}\\\
\nabla_{\theta_{y}}f_{yi}\end{smallmatrix}\right]$,
$\eta_{k}=\ell^{\prime}_{ik}/\sqrt{\ell^{\prime\prime}_{ik}}$. The RLS
iterations in (17) provide the optimal solution $p^{h}\triangleq p_{N-1}$ (as
well as $P_{N-1}=((A^{h})^{\prime}A^{h})^{-1})$ when the matrix
$P_{k}\in{\mathbb{R}}^{n_{p}\times n_{p}}$ and vector $p_{k}$ are initialized,
respectively, as

$P_{-1}=\left[\begin{smallmatrix}\frac{1}{\rho_{x}}I_{n_{x}}&0&0\\\
0&\frac{1}{\rho_{\theta}}I_{n_{\theta_{x}}}&0\\\
0&0&\frac{1}{\rho_{\theta}}I_{n_{\theta_{x}}}\end{smallmatrix}\right],\quad
p_{-1}=\left[\begin{smallmatrix}x_{0}^{h}\\\ \theta_{y}^{h}\\\
\theta_{x}^{h}\end{smallmatrix}\right]$ (18)

Different algorithms were proposed in the literature to implement RLS
efficiently, such as the inverse QR method proposed in [1], in which a
triangular matrix $R_{k}$ is updated recursively starting from
$R_{-1}=(P_{-1})^{\frac{1}{2}}$.

###  3.3 Line search

To enforce the decrease of the cost function (LABEL:eq:training-cost-reg)
across iterations $h=0,1,\ldots$, a possibility is to update the solution as

$\left[\begin{smallmatrix}x_{0}^{h+1}\\\ \theta_{x}^{h+1}\\\
\theta_{y}^{h+1}\end{smallmatrix}\right]=\left[\begin{smallmatrix}x_{0}^{h}\\\
\theta_{x}^{h}\\\
\theta_{y}^{h}\end{smallmatrix}\right]+\alpha_{h}\left[\begin{smallmatrix}p^{h}_{x_{0}}\\\
p^{h}_{\theta_{x}}\\\ p^{h}_{\theta_{y}}\end{smallmatrix}\right]$ (19)

and set $x_{k+1}^{h+1}=f_{x}(x_{k}^{h+1},u_{k},\theta_{x}^{h+1})$,
$k=1,\ldots,N-1$. We choose the step-size $\alpha_{h}$ by finding the largest
$\alpha$ satisfying the Armijo condition

$V^{h+1}\leq V^{h}+c_{1}\alpha(\nabla V^{h})^{\prime}p^{h}$ (20)

where $V^{h+1}$, $V^{h}$ are the costs defined as in (3) associated with the
corresponding combinations of $(x_{0},\theta_{x},\theta_{y})$, the gradient is
taken with respect to $\left[\begin{smallmatrix}x_{0}\\\ \theta_{x}\\\
\theta_{y}\end{smallmatrix}\right]$, and $c_{1}$ is a constant (e.g.,
$c_{1}=10^{-4}$) [23, Chapter 3.1]. In this paper we use the common approach
of starting with $\alpha=1$ and then change $\alpha\leftarrow\sigma\alpha$,
$\sigma\in(0,1)$ until (20) is satisfied. We set a bound $N_{\sigma}$ on the
number of line-search steps that can be performed, setting $\alpha=0$ in the
case of line-search failure. Since $V(x_{0}+\alpha
p_{x_{0}}^{h},\theta_{x}+\alpha p_{\theta_{x}}^{h},\theta_{y}+\alpha
p_{\theta_{y}}^{h})\approx
V(x_{0},\theta_{x}^{h},\theta_{y}^{h})+\frac{1}{2}\|A^{h}\alpha
p^{h}-b^{h}\|_{2}^{2}-\frac{1}{2}\|b^{h}\|_{2}^{2}$ for small values of
$\alpha$, the directional derivative in (20) can be computed from
$A^{h},b^{h}$:

$\begin{array}[]{rcl}(\nabla
V^{h})^{\prime}p^{h}&=&\underset{\alpha\rightarrow
0}{\lim}\frac{V(x_{0}+\alpha p_{x_{0}}^{h},\theta_{x}+\alpha
p_{\theta_{x}}^{h},\theta_{y}+\alpha
p_{\theta_{y}}^{h})-V(x_{0},\theta_{x}^{h},\theta_{y}^{h})}{\alpha}\\\
&=&\underset{\alpha\rightarrow 0}{\lim}\frac{1}{2}\frac{\|A^{h}\alpha
p^{h}-b^{h}\|_{2}^{2}-\|b^{h}\|_{2}^{2}}{\alpha}=-(b^{h})^{\prime}A^{h}p^{h}\end{array}$

Note that when the loss function
$\ell(y,\hat{y})=\frac{1}{2}\|y-\hat{y}\|_{2}^{2}$ and the regularization
terms $r_{x},r_{\theta_{y}},r_{\theta_{x}}$ are quadratic, the loss
$V(x_{0},\theta_{x}^{h},\theta_{y}^{h})=\frac{1}{2}(b^{h})^{\prime}b^{h}$ by
construction.

The overall nonlinear least-squares solution method for training recurrent
neural networks described in the previous sections is summarized in Algorithm
1.

Algorithm 1 belongs to the class of Generalized Gauss-Newton (GGN) methods
[22] applied to solve problem (5), with the addition of line search. In fact,
following the notation in [22], we can write
$V(x_{0},\theta_{x},\theta_{y})=\Phi_{0}(F(w))$ where
$w\triangleq[x_{0}^{\prime}\ {\theta_{x}}^{\prime}\ {\theta_{y}}]^{\prime}$ is
the optimization vector, function
$\Phi_{0}:{\mathbb{R}}^{N_{\phi}}\to{\mathbb{R}}$, $N_{\phi}\triangleq
n_{x}+n_{\theta_{x}}+n_{\theta_{y}}+Nn_{y}$, is given by the composition
$\Phi_{0}=\Phi_{+}\circ\Phi_{\ell}$ of the summation function
$\Phi_{+}:{\mathbb{R}}^{N_{\phi}}\to{\mathbb{R}}$,
$\Phi_{+}(x)=\sum_{i=1}^{N_{\phi}}x_{i}$, with the vectorized loss
$\Phi_{\ell}:{\mathbb{R}}^{N_{\phi}}\to{\mathbb{R}}^{N_{\phi}}$,
$\Phi_{\ell}([x_{0}^{\prime}\ \theta_{x}^{\prime}\ \theta_{y}^{\prime}\
\hat{y}_{0}^{\prime}\ \ldots\ \hat{y}_{N-1}^{\prime}]^{\prime})=[r_{x}(x_{0})\
r_{\theta}^{x}(\theta_{x})\ r_{\theta}^{y}(\theta_{y})\
\ell(y_{0},\hat{y}_{0})^{\prime}\
\ldots\ell(y_{N-1},\hat{y}_{N-1})^{\prime}]^{\prime}$. By assumption, all the
components of $\Phi_{\ell}$ are strongly convex functions, and since
$\Phi_{+}$ is also strongly convex, $\Phi_{0}$ is strongly convex with respect
to its argument $F(w)$. As in GGN methods, Algorithm 1 takes a second-order
approximation of $\Phi_{0}$ and a linearization of $F(w)$ during the
iterations. In [22], local convergence under the full step $\alpha_{h}\equiv
1$, $\forall h\geq 0$, is shown under additional assumptions on $\Phi_{0}$ and
the neglected second-order derivatives of $F$, while in Algorithm 1 we force
monotonic decrease of the loss function $V$ at each iteration $h$ by line
search, hence obtaining a “damped” GGN algorithm. Note that the latter may
fail ($\alpha_{h}=0$) in the case $p^{h}$ is not a descent direction or
$N_{\sigma}$ is not large enough. In such a case, Algorithm 1 would stop
immediately as $V^{h}-V^{h-1}=0$.

Algorithm 1 Sequential least squares with line search.

Input: Training dataset $\\{u_{k},y_{k}\\}_{k=0}^{N-1}$, initial guess
$(x_{0}^{0}$, $\theta_{x}^{0}$, $\theta_{y}^{0})$, maximum number $E$ of
epochs, tolerance $\epsilon_{V}>0$, line-search parameters $c_{1}>0$,
$\sigma\in(0,1)$, $N_{\sigma}\geq 1$.

  1. 1.

Simulate (1) to get $x_{k+1}^{0}$, $k=0,\ldots,N-1$, and the initial cost
$V^{0}$ from (5);

  2. 2.

$h\leftarrow 0$;

  3. 3.

repeat:

    1. 0..1.

Construct and solve the linear least-squares problem (15) to get a descent
direction $p^{h}$;

    2. 0..2.

Find step-length $\alpha_{h}\in[0,1]$ by line-search using Armijo rule (20)
and get the new cost $V^{h+1}$;

    3. 0..3.

Update $x_{0}^{h+1},\theta_{x}^{h+1},\theta_{y}^{h+1}$ as in (19);

    4. 0..4.

$h\leftarrow h+1$;

  4. 4.

until $V^{h}-V^{h-1}\leq\epsilon_{V}$ or $h\geq E$;

  5. 5.

end.

Output: RNN parameters $\theta_{x}^{h},\theta_{y}^{h}$ and initial hidden
state $x_{0}^{h}$.

###  3.4 Initialization

In this paper we initialize $x_{0}^{0}=0$, set zero initial bias terms, and
draw the remaining components of $\theta_{x}^{0}$ and $\theta_{y}^{0}$ from
the normal distribution with zero mean and standard deviation defined as in
[7], further multiplied by the quantity $\sigma_{0}\leq 1$. The rationale to
select small initial values of $\theta_{x}^{0}$, $\theta_{y}^{0}$ is to avoid
that the initial dynamics $f_{x}$ used to compute $x_{k}^{0}$ for
$k=1,\ldots,N-1$ are unstable. Then, since the overall loss $V^{h}$ is
decreasing, the divergence of state-trajectories is unlikely to occur at
subsequent iterations. Note that this makes the elimination of the state
variables $x_{1},\ldots,x_{N-1}$ in the condensing approach (4) suitable for
the current RNN training setting, contrarily to solving nonlinear MPC problems
in which the initial state $x_{0}$ and the model coefficients $\theta_{x}$,
$\theta_{y}$ are fixed and may therefore lead to excite unstable linearized
model responses, with consequent numerical issues [5]. Note also that ideas to
enforce open-loop stability could be introduced in our setting by adopting
ideas as in [12], where the authors propose to also learn a Lyapunov function
while learning the model.

###  3.5 A Levenberg-Marquardt variant

Trust-region methods are an alternative family of approaches to line-search
for forcing monotonicity of the objective function. In solving nonlinear
least-squares problems, the well-known Levenberg-Marquardt (LM) algorithm [13,
20] is often used due to its good performance. In LM methods, (15) is changed
to

$p^{h}=\arg\min_{p}\frac{1}{2}\|A^{h}p-b^{h}\|_{2}^{2}+\frac{1}{2}\lambda_{h}\|p\|_{2}^{2}$
(21)

where the regularization term $\lambda_{h}$ is tuned at each iteration $h$ to
meet the condition $V^{h}\leq V^{h-1}$. It is easy to show that tuning
$\lambda_{h}$ is equivalent to tuning the trust-region radius, see e.g. [23,
Chapter 10]. Clearly, setting $\lambda_{h}=0$ corresponds to the full
Generalized Gauss-Newton step, as in line-search when $\alpha_{h}=1$, and
$\lambda_{h}\rightarrow\infty$ corresponds to $p^{h}\rightarrow 0$
($\alpha_{h}\rightarrow 0$). Several methods exist to update $\lambda_{h}$
(and possibly generalizing $\|p\|_{2}^{2}$ to $\|D^{h}p\|_{2}^{2}$ for some
weight matrix $D^{h}$) [33]. In this paper, we take the simplest approach
suggested in [33, Section 2.1] and, starting from $\lambda_{h}=\lambda_{h-1}$,
keep multiplying $\lambda_{h}\leftarrow c_{2}\lambda_{h}$ until the condition
$V^{h}\leq V^{h-1}$ is met, then reduce
$\lambda_{h}\leftarrow\lambda_{h}/c_{3}$ before proceeding to step $h+1$. The
approach is summarized in Algorithm 2.

Algorithm 2 Sequential least squares, Levenberg-Marquardt variant.

Input: Training dataset $\\{u_{k},y_{k}\\}_{k=0}^{N-1}$, initial guess
$(x_{0}^{0}$, $\theta_{x}^{0}$, $\theta_{y}^{0})$, maximum number $E$ of
epochs, tolerance $\epsilon_{V}>0$, initial value $\lambda_{0}>0$, LM
parameters $c_{2},c_{3}>1$, number $N_{\lambda}\geq 1$ of search iterations.

  1. 1.

Simulate (1), get $x_{k+1}^{0}$, $k=0,\ldots,N-1$, and initial cost $V^{0}$
from (5);

  2. 2.

$h\leftarrow 0$;

  3. 3.

repeat:

    1. 5..1.

Construct the least-squares problem (15);

    2. 5..2.

$j\leftarrow 0$;

    3. 5..3.

repeat:

      1. 5..3.1.

Solve the regularized linear least-squares problem (21) and get a candidate
$p^{h}$;

      2. 5..3.2.

Update $x_{0}^{h+1},\theta_{x}^{h+1},\theta_{y}^{h+1}$ as in (19) when
$\alpha_{h}=1$;

      3. 5..3.3.

Get new cost $V^{h+1}$;

      4. 5..3.4.

if $V^{h+1}\leq V^{h}$ decrease $\lambda_{h}\leftarrow\lambda_{h}/c_{3}$;

      5. 5..3.5.

else increase $\lambda_{h}\leftarrow c_{2}\lambda_{h}$;

      6. 5..3.6.

set $j\leftarrow j+1$;

    4. 5..4.

until $V^{h+1}\leq V^{h}$ or $j=N_{\lambda}$;

    5. 5..5.

$h\leftarrow h+1$;

    6. 5..6.

if $V^{h}>V^{h-1}$ go to Step 5.;

  4. 4.

until $V^{h}-V^{h-1}\leq\epsilon_{V}$, or $h\geq E$;

  5. 5.

end.

Output: RNN parameters $\theta_{x}^{h},\theta_{y}^{h}$ and initial hidden
state $x_{0}^{h}$.

Similarly to Algorithm 1, in Algorithm 2 we have set a bound $N_{\lambda}$ on
the number of least-squares problems that can be solved at each epoch $h$. In
the case an improvement $V^{h}\leq V^{h-1}$ has not been found, Algorithm 2
stops due to Step 5.5..6.

The main difference between performing line-search and computing instead LM
iterations is that at each epoch $h$ the former only solves the linear least-
squares problems (15), the latter requires solving up to $N_{\lambda}$
problems (21). Note that, however, both only construct $A^{h}$ and $b^{h}$
once at each iteration $h$, which, as remarked earlier, is usually the
dominant computation effort due to the evaluation of the required Jacobian
matrices $\nabla_{x}f_{x}$, $\nabla_{\theta_{x}}f_{x}$, $\nabla_{x}f_{y}$,
$\nabla_{\theta_{y}}f_{y}$.

##  4 Non-smooth regularization

We now want to solve the full problem (3) by considering also the additional
non-smooth regularization term $g(\theta_{x},\theta_{y})$. Consider again the
condensed loss function (5) and rewrite problem (3) according to the following
split

$\displaystyle\min_{x_{0},\theta_{x},\theta_{y},\nu_{x},\nu_{y}}$
$\displaystyle V(x_{0},\theta_{x},\theta_{y})+g(\nu_{x},\nu_{y})$
$\displaystyle\mathop{\rm s.t.}\nolimits$
$\displaystyle\theta_{x}=\nu_{x},\quad\theta_{y}=\nu_{y}$ (22)

We solve problem (22) by executing the following scaled ADMM iterations

$\displaystyle\left[\begin{smallmatrix}x_{0}(t+1)\\\ \theta_{x}(t+1)\\\
\theta_{y}(t+1)\end{smallmatrix}\right]$ $\displaystyle=$
$\displaystyle\arg\min_{x_{0},\theta_{x},\theta_{y}}V(x_{0},\theta_{x},\theta_{y})+\frac{\rho}{2}\left\|\left[\begin{smallmatrix}\theta_{x}-\nu_{x}(t)+w_{x}(t)\\\
\theta_{y}-\nu_{y}(t)+w_{y}(t)\end{smallmatrix}\right]\right\|_{2}^{2}$
$\displaystyle\left[\begin{smallmatrix}\nu_{x}(t+1)\\\
\nu_{y}(t+1)\end{smallmatrix}\right]$ $\displaystyle=$
$\displaystyle\arg\min_{\nu_{x},\nu_{y}}\frac{1}{\rho}g(\nu_{x},\nu_{y})+\frac{1}{2}\left\|\left[\begin{smallmatrix}\nu_{x}-\theta_{x}(t+1)-w_{x}(t)\\\
\nu_{y}-\theta_{y}(t+1)-w_{y}(t)\end{smallmatrix}\right]\right\|_{2}^{2}$
(23b) $\displaystyle=$ $\displaystyle\mathop{\rm
prox}\nolimits_{\frac{1}{\rho}g}(\theta_{x}(t+1)+w_{x}(t),\theta_{y}(t+1)+w_{y}(t))$
$\displaystyle\left[\begin{smallmatrix}w_{x}(t+1)\\\
w_{y}(t+1)\end{smallmatrix}\right]$ $\displaystyle=$
$\displaystyle\left[\begin{smallmatrix}w_{x}^{h}+\theta_{x}(t+1)-\nu_{x}(t+1)\\\
w_{y}^{h}+\theta_{y}(t+1)-\nu_{y}(t+1)\end{smallmatrix}\right]$ (23c)

for $t=0,\ldots,N_{\rm ADMM}-1$, where “$\mathop{\rm prox}\nolimits$” denotes
the proximal operator. Problem (LABEL:eq:ADMM-x) is solved by running either
Algorithm 1 or 2 with initial condition $x_{0}^{0}=x_{0}(t)$,
$\theta_{x}^{0}=\theta_{x}(t)$, $\theta_{y}^{0}=\theta_{y}(t)$. Note that when
solving the least-squares problem (15) or (21), $A^{h},b^{h}$ are augmented to
include the additional $\ell_{2}$-regularization term
$\frac{1}{2}\left\|\sqrt{\rho}\left[\begin{smallmatrix}p_{\theta_{x}}\\\
p_{\theta_{y}}\end{smallmatrix}\right]-\sqrt{\rho}\left[\begin{smallmatrix}\nu_{x}(t)-\theta_{x}^{h}-w_{x}(t)\\\
\nu_{y}(t)-\theta_{y}^{h}-w_{y}(t)\end{smallmatrix}\right]\right\|_{2}^{2}$
introduced in the ADMM iteration (LABEL:eq:ADMM-x). Note also that for $t>0$
Step 5. of Algorithms 1 and 2 can be omitted, since the initial simulated
trajectory $x_{k+1}^{0}$ is available from the last execution of the ADMM step
(LABEL:eq:ADMM-x), $\forall k=0,\ldots,N-2$, and $V_{0}$ can be obtained by
adjusting the cost computed at the last iteration of the previous execution of
Algorithm 1 or 2 to take into account that $\nu_{x}$, $\nu_{y}$, $w_{x}$,
$w_{y}$ have possibly changed their value.

Our numerical experiments show that setting $E=1$ is a good choice, which
corresponds to just computing one sequential LS iteration to solve
(LABEL:eq:ADMM-x) approximately. The overall algorithm, that we call NAILS
(nonconvex ADMM iterations and least squares with line search) in the case
Algorithm 1 is used, or NAILM (nonconvex ADMM iterations and least squares
with Levenberg-Marquardt variant) when instead Algorithm 2 is used, is
summarized in Algorithm 3. Note that NAILS/NAILM processes the training
dataset for a maximum of $E\cdot N_{\rm ADMM}$ epochs.

Algorithm 3 Nonconvex ADMM iterations and least squares, either with a line-
search (NAILS) or Levenberg-Marquardt (NAILM) approach

Input: Training dataset $\\{u_{k},y_{k}\\}_{k=0}^{N-1}$, initial guess
$x_{0}^{0}$, $\theta_{x}^{0}$, $\theta_{y}^{0}$, number $N_{\rm ADMM}$ of ADMM
iterations, ADMM parameter $\rho>0$, maximum number $E$ of epochs, tolerance
$\epsilon_{V}>0$; line-search parameters $c_{1}>0$, $\sigma\in(0,1)$, and
$N_{\sigma}>1$ (NAILS), or LM parameters $\lambda_{0}>0$, $c_{2},c_{3}>1$, and
$N_{\lambda}>0$ (NAILM).

  1. 1.

$\left[\begin{smallmatrix}\nu_{x}^{0}\\\
\nu_{y}^{0}\end{smallmatrix}\right]\leftarrow\left[\begin{smallmatrix}\theta_{x}^{0}\\\
\theta_{y}^{0}\end{smallmatrix}\right]$,
$\left[\begin{smallmatrix}w_{x}^{0}\\\
w_{y}^{0}\end{smallmatrix}\right]\leftarrow 0$;

  2. 2.

for $t=0,1,\ldots,N_{\rm ADMM}-1$ do:

    1. 5..1.

run Algorithm 1 or 2 to update $x_{0}(t+1)$, $\theta_{x}(t+1)$, and
$\theta_{y}(t+1)$ as in (LABEL:eq:ADMM-x);

    2. 5..2.

update $\left[\begin{smallmatrix}\nu_{x}(t+1)\\\
\nu_{y}(t+1)\end{smallmatrix}\right]$ as in (23b) and
$\left[\begin{smallmatrix}w_{x}(t+1)\\\ w_{y}(t+1)\end{smallmatrix}\right]$ as
in (23c);

    3. 5..3.

$t\leftarrow t+1$;

  3. 3.

end.

Output: RNN parameters $\theta_{x}=\nu_{x}(t),\theta_{y}=\nu_{y}(t)$ and
initial hidden state $x_{0}(t)$.

In NAILS, if RLS are used in alternative to forming and solving (15), since

$\arg\min_{\alpha\in{\mathbb{R}}}\frac{1}{2}(\rho_{1}(\alpha-\alpha_{1})^{2}+\rho_{2}(\alpha-\alpha_{2})^{2})=\arg\min_{\alpha\in{\mathbb{R}}}\frac{\rho_{1}+\rho_{2}}{2}\left(\alpha-\frac{\rho_{1}\alpha_{1}+\rho_{2}\alpha_{2}}{\rho_{1}+\rho_{2}}\right)^{2}$

for any $\rho_{1},\rho_{2}>0$ and $\alpha_{1},\alpha_{2}\in{\mathbb{R}}$, we
can initialize

$P_{-1}=\left[\begin{smallmatrix}\frac{1}{\rho_{x}}I_{n_{x}}&0&0\\\
0&\frac{1}{\rho_{\theta}+\rho}I_{n_{\theta_{x}}}&0\\\
0&0&\frac{1}{\rho_{\theta}+\rho}I_{n_{\theta_{y}}}\end{smallmatrix}\right],\
p_{-1}=\left[\begin{smallmatrix}x_{0}^{h}\\\
\frac{\rho_{\theta}\theta_{x}^{h}+\rho(\nu_{x}(t)-\theta_{x}^{h}-w_{x}(t))}{\rho_{\theta}+\rho}\\\
\frac{\rho_{\theta}\theta_{y}^{h}+\rho(\nu_{y}(t)-\theta_{y}^{h}-w_{y}(t))}{\rho_{\theta}+\rho}\end{smallmatrix}\right]$

The final model parameters are given by $\theta_{x}=\nu_{x}(t)$,
$\theta_{y}=\nu_{y}(t)$, $x_{0}(t)$, where $t=N_{\rm ADMM}$ is the last ADMM
iteration executed. In case $x_{0}$ also needs to be regularized by a non-
smooth penalty, we can extend the approach to include a further split
$\nu_{x_{0}}=x_{0}$.

Algorithm 3 can be modified to return the initial state and model obtained
after each ADMM step that provide the best loss $\sum\ell(y_{k},\hat{y}_{k})$
observed during the execution of the algorithm on training data or, if
available, on a separate validation dataset.

Note that the iterations (23) are not guaranteed to converge to a global
optimum of the problem. Moreover, it is well known that in some cases
nonconvex ADMM iterations may even diverge [9]. The reader is referred to [9,
8, 37, 32] for convergence conditions and/or alternative ADMM formulations.

We finally remark that, in the absence of the non-smooth regularization term
$g$, Algorithm 3 simply reduces to a smooth GNN method running with tolerance
$\epsilon_{V}$ for a maximum number $E$ of epochs.

###  4.1 $\ell_{1}$-regularization

A typical instance of regularization is the $\ell_{1}$-penalty

$g(\theta_{x},\theta_{y})=\tau_{x}\|\theta_{x}\|_{1}+\tau_{y}\|\theta_{y}\|_{1}$
(24)

to attempt pruning an overly-parameterized network structure, with
$\tau_{x},\tau_{y}\geq 0$. In this case, the update in (23b) simply becomes

$\begin{array}[]{rcl}\nu_{x}(t+1)&=&S_{\frac{\tau_{x}}{\rho}}(\theta_{x}(t+1)+w_{x}^{h})\\\
\nu_{y}(t+1)&=&S_{\frac{\tau_{y}}{\rho}}(\theta_{y}(t+1)+w_{y}^{h})\end{array}$
(25)

where $S_{\alpha}$ is the soft-threshold operator

$[S_{\alpha}(w)]_{i}=\left\\{\begin{array}[]{rll}w_{i}+\alpha&\mbox{if}&w_{i}\leq-\alpha\\\
0&\mbox{if}&|w_{i}|\leq\alpha\\\
w_{i}-\alpha&\mbox{if}&w_{i}\geq\alpha\end{array}\right.$

###  4.2 Group-Lasso regularization and model reduction

The complexity of the structure of the networks $f_{x}$, $f_{y}$ can be
reduced by using group-Lasso penalties [42] that attempt zeroing entire
subsets of parameters. In particular, the order $n_{x}$ of the RNN can be
penalized by creating $n_{x}$ groups $\theta_{i}^{g}$ of variables,
$\theta_{i}^{g}\in{\mathbb{R}}^{n_{g}}$,
$n_{g}=n_{1}^{x}+n_{1}^{y}+n_{L_{x}-1}^{x}+1$, $i=1,\ldots,n_{x}$, each one
stacking the $i$th columns of $A_{1}^{x}$ and $A_{1}^{y}$, the $i$th row of
$A_{L_{x}}$, and the $i$th entry of $b_{L_{x}}$. In this case, in (22) we only
introduce the splitting $\theta_{i}^{g}=\nu_{i}^{g}$ on those variables and
penalize

$g(\nu_{i}^{g})=\tau_{g}\sum_{i=1}^{n_{x}}\|\nu_{i}^{g}\|_{2}$

where $\tau_{g}\geq 0$ is the group-Lasso penalty. Then, vectors
$\nu_{i}^{g}(t+1)$ are updated as in (25) by replacing
$S_{\frac{\tau_{g}}{\rho}}$ with the block soft thresholding operator
$S_{\frac{\tau_{g}}{\rho}}^{g}:{\mathbb{R}}^{n_{g}}\to{\mathbb{R}}^{n_{g}}$
defined as

$S_{\alpha}^{g}(w)=\left\\{\begin{array}[]{rll}(1-\frac{\alpha}{\|w\|_{2}})w&\mbox{if}&\|w\|_{2}>\alpha\\\
0&\mbox{if}&\|w\|_{2}\leq\alpha\end{array}\right.$

The larger the penalty $\tau_{g}$, the larger is expected the number of state
indices $i$ such that the $i$th column of $A_{1}^{x}$ and $A_{1}^{y}$, the
$i$th row of $A_{L_{x}}$, and the $i$th entry of $b_{L_{x}}$ are all zero. As
a consequence, the corresponding $i$th hidden state has no influence in the
model.

###  4.3 $\ell_{0}$-regularization

The $\ell_{0}$-regularization term

$g(\theta_{x},\theta_{y})=\tau_{x}\|\theta_{x}\|_{0}+\tau_{y}\|\theta_{y}\|_{0}$
(26)

with $\tau_{x}^{0},\tau_{y}^{0}\geq 0$ also admits the explicit proximal
“hard-thresholding” operator

$[\mathop{\rm
prox}\nolimits_{\frac{\tau}{\rho}\|x\|_{0}}]_{i}=\left\\{\begin{array}[]{rll}0&\mbox{if}&x_{i}^{2}<\frac{2\rho}{\tau}\\\
x_{i}&\mbox{if}&x_{i}^{2}\geq\frac{2\rho}{\tau}\end{array}\right.$

and can be used in alternative to (24).

###  4.4 Integer constraints

Constraints on $\theta_{x}$, $\theta_{y}$ that can expressed as mixed-integer
linear inequalities

$A_{\rm MI}\left[\begin{smallmatrix}\theta_{x}\\\
\theta_{y}\end{smallmatrix}\right]\leq b_{\rm
MI},\quad{\theta_{x}}_{i},{\theta_{y}}_{j}\in\\{0,1\\},\ \forall i\in I_{x},\
j\in I_{y}$ (27)

where $I_{x}\subseteq\\{1,\ldots,n_{\theta_{x}}\\}$,
$I_{y}\subseteq\\{1,\ldots,n_{\theta_{y}}\\}$, $A_{\rm
MI}\in{\mathbb{R}}^{n_{\rm MI}\times(n_{\theta_{x}}+n_{\theta_{y}})}$, $b_{\rm
MI}\in{\mathbb{R}}^{n_{\rm MI}}$, $n_{\rm MI}\geq 0$, can be handled by
setting $g$ as the indicator function

$g(\theta_{x},\theta_{y})=\left\\{\begin{array}[]{rl}0&\mbox{if~{}\eqref{eq:MI-
constraints} is satisfied}\\\ +\infty&\mbox{otherwise}\end{array}\right.$

In this case, $\nu_{x}^{h+1},\nu_{y}^{h+1}$ can be computed from (23b) by
solving the mixed-integer quadratic programming (MIQP) problem

$\begin{array}[]{rl}\min_{\nu_{x},\nu_{y}}&\frac{1}{2}\left\|\left[\begin{smallmatrix}\nu_{x}-\theta_{x}(t+1)-w_{x}(t)\\\
\nu_{y}-\theta_{y}(t+1)-w_{y}(t)\end{smallmatrix}\right]\right\|_{2}^{2}\\\
\mathop{\rm s.t.}\nolimits&A_{\rm MI}\left[\begin{smallmatrix}\theta_{x}\\\
\theta_{y}\end{smallmatrix}\right]\leq b_{\rm MI}\\\
&{\theta_{x}}_{i},{\theta_{y}}_{j}\in\\{0,1\\}\ \forall i\in I_{x},\ j\in
I_{y}\end{array}$ (28)

Note that in the special case $n_{\rm MI}=0$, (28) has the following explicit
solution [29]

$\begin{array}[]{rcl}\nu_{xi}(t+1)&=&\mathop{\rm
round}\nolimits({\theta_{x}}_{i}(t+1)+w_{xi}(t)),\ i\in I_{x}\\\
\nu_{yj}(t+1)&=&\mathop{\rm round}\nolimits({\theta_{x}}_{j}(t+1)+w_{xj}(t)),\
j\in I_{y}\end{array}$ (29)

where $\mathop{\rm round}\nolimits(\alpha)=0$ if $\alpha<0.5$, or $1$
otherwise.

Binary constraints can be extended to more general quantization constraints
${\theta_{x}}_{i},{\theta_{y}}_{j}\in{\mathcal{Q}}\triangleq\\{q_{1},\ldots,q_{n_{Q}}\\}\subset{\mathbb{R}}$,
$\forall i\in I_{x}$, $\forall j\in I_{y}$. In fact, these can be handled as
in (29) by setting $\nu_{xi}(t+1)$ equal to the value in ${\mathcal{Q}}$ that
is closest to ${\theta_{x}}_{i}(t+1)+w_{xi}^{h}$, and similarly for
$\nu_{yj}(t+1)$.

Figure 1: MSE loss (5) evaluated on training data as a function of training time on different runs $L_{x}-1$ | number of hidden layers in $f_{x}$
---|---
$(n_{1}^{x},\ldots,n_{L_{x}-1}^{x})$ | number of neurons per layer in $f_{x}$
$n_{\theta_{x}}$ | dimension of parameter vector $\theta_{x}$
$L_{y}-1$ | number of hidden layers in $f_{y}$
$(n_{1}^{y},\ldots,n_{L_{y}-1}^{y})$ | number of neurons per layer in $f_{y}$
$f_{L_{y}}^{y}$ | nonlinear output function in $f_{y}$
$n_{\theta_{y}}$ | dimension of parameter vector $\theta_{y}$
$\sigma_{0}$ | scaling factor used in random initialization [7]
$\rho_{x}$ | penalty on $\|x_{0}\|_{2}^{2}$ (or on $\|\theta_{x0}\|_{2}^{2}$)
$\rho_{\theta}$ | penalty on $\|\theta_{x}\|_{2}^{2}$, $\|\theta_{y}\|_{2}^{2}$
$\tau_{x}$ ($\tau_{y}$) | penalty on $\|\theta_{x}\|_{1}$ ($\|\theta_{y}\|_{1}$)
$\tau_{g}$ | group-Lasso penalty
$\tau_{x}^{0}$ ($\tau_{y}^{0}$) | penalty on $\|\theta_{x}\|_{0}$ ($\|\theta_{y}\|_{0}$)
$E$ | max number of epochs when solving
| the nonlinear LS problem
$\epsilon_{V}$ | optimality threshold
$c_{1},\sigma,N_{\sigma}$ | line-search parameters
$c_{2},c_{3},N_{\lambda},\lambda_{0}$ | LM parameters
$N_{\rm ADMM}$ | number of ADMM iterations
$\rho$ | ADMM parameter
Table 1: Hyper-parameters of Algorithm 3 (NAILS/NAILM)

##  5 Numerical experiments

We test NAILS and NAILM on three nonlinear system identification problems. The
hyper-parameters of the Algorithm 3 are summarized in Table 1. All
computations have been carried out in MATLABR2022b on an Apple M1 Max CPU,
with the library CasADi [2] used to generate Jacobian matrices via automatic
differentiation. The scaling factor $\sigma_{0}=0.15$ is used to initialize
the weights of the neural networks in all the experiments. Unless the state
encoder (6) is used, after training a model the best initial condition $x_{0}$
is computed by solving the small-scale non-convex simulation-error
minimization problem

$\min_{x_{0}}\rho_{x}\|x_{0}\|_{2}^{2}+\sum_{k=0}^{N_{0}}\ell(y_{k},\hat{y}_{k})$
(30)

by using the particle swarm optimizer (PSO) PSwarm [34] with initial
population of $2n_{x}$ samples, each component of $x_{0}$ constrained in
$[-3,3]$, and with $N_{0}=100$.

###  5.1 RNN training with smooth quadratic loss

We first test NAILS/NAILM against the EKF approach [4] and gradient-descent
based on the AMSGrad algorithm [26], for which we compute the gradient as in
(3.2), on real-world data from the nonlinear magneto-rheological fluid damper
benchmark proposed in [35]. We use $N=2000$ data samples for training and
$1499$ samples for testing, obtained from the System Identification (SYS-ID)
Toolbox for MATLAB R2021a [16], where a nonlinear autoregressive (NLARX) model
identification algorithm is employed for black-box nonlinear modeling.

We consider a RNN model (2b) with $n_{x}=4$ hidden states and shallow state-
update and output network functions ($L_{x}=L_{y}=2$) with
$n_{1}^{x}=n_{1}^{y}=4$ neurons, hyperbolic-tangent activation functions
$f_{1}^{x},f_{1}^{y}$, and linear output function $f_{2}^{y}$, parameterized
by $\theta_{x}\in{\mathbb{R}}^{44}$, $\theta_{y}\in{\mathbb{R}}^{29}$. The CPU
time to solve (30) is approximately 30 ms in all tests.

In (3) we use the quadratic loss
$\ell(y_{k},\hat{y}_{k})=\frac{1}{N}\|y_{k}-\hat{y}_{k}\|_{2}^{2}$ and
quadratic regularization
$r(x_{0},\theta_{x},\theta_{y})=\frac{1}{2}(\rho_{x}\|x_{0}\|_{2}^{2}+\rho_{\theta}\|\theta_{x}\|_{2}^{2}+\rho_{\theta}\|\theta_{x}\|_{2}^{2})$
with $\rho_{x}=1$, $\rho_{\theta}=0.1$. As $g(\theta_{x},\theta_{y})=0$, NAILS
and NAILM only execute Algorithm 1 once with parameters $\epsilon_{V}=10^{-6}$
and, respectively, $c_{1}=10^{-4}$, $\sigma=0.5$, $N_{\sigma}=20$, and
$\lambda_{0}=100$, $c_{2}=1.5$, $c_{3}=5$, $N_{\lambda}=20$. The initial
condition is $x_{0}^{0}=0$ and the components of $\theta_{x}^{0}$,
$\theta_{y}^{0}$ are randomly generated as described in Section 3.4.

Table 2 shows the mean and standard deviation (computed over 20 experiments,
each one starting from a different random set of weights and zero bias terms)
of the final best fit rate BFR = $100(1-\|Y-\hat{Y}\|_{2}/\|Y-\bar{y}\|_{2})$,
where $Y$ is the vector of measured output samples, $\hat{Y}$ the vector of
output samples simulated by the identified model fed in open-loop with the
input data, and $\bar{y}$ is the mean of $Y$, achieved on training and test
data.

BFR | training | test
---|---|---
NAILS | 94.41 (0.27) | 89.35 (2.63)
NAILM | 94.07 (0.38) | 89.64 (2.30)
EKF | 91.41 (0.70) | 87.17 (3.06)
SGD | 84.69 (0.15) | 80.56 (0.18)
Table 2: Fluid damper benchmark: mean (standard deviation) of BFR on training
and test data

For illustration, the evolution of the MSE loss
$\frac{1}{2N}\sum_{k=0}^{N-1}(y_{k}-\hat{y}_{k})^{2}$ as a function of
training time is shown in Figure 1 for each run. The average CPU time per
epoch spent by the training algorithm is 54.3 ms (NAILS), 72.2 ms (NAILM),
253.3 ms (EKF), and 34.6 ms (AMSGrad).

It is apparent that NAILS and NAILM obtain similar quality of fit and are
computationally comparable. They get better fit results than EKF on average
and, not surprisingly, outperform gradient descent due to their second-order
approximation nature. On the other hand, EKF converges more quickly to a good
quality of fit, as the model parameters $\theta_{x}$, $\theta_{y}$ are updated
_within_ each epoch, rather than _after_ each epoch is processed.

####  5.1.1 Silverbox benchmark problem

We test the proposed training algorithm on the Silverbox dataset, a popular
benchmark for nonlinear system identification, for which we refer the reader
to [39] for a detailed description. The first 40000 samples are used as the
test set, the remaining samples, which correspond to a sequence of 10
different random odd multi-sine excitations, are manually split in $M=10$
different training traces of about 8600 samples each (the intervals between
each different excitation are removed from the training set, as they bring no
information).

We consider a RNN model (2b) with $n_{x}=8$ hidden states, no I/O feedthrough
($A_{1}^{y}\in{\mathbb{R}}^{n_{1}^{y}\times n_{x}}$ in (2b)), and state-update
and output network functions with three hidden layers ($L_{x}=L_{y}=4$) with
$8$ neurons each, a neural network model $f_{x0}$ (6) with two hidden layers
of $4$ neurons each mapping the vector $v_{0}$ of past $n_{a}=8$ outputs and
$n_{b}=8$ inputs to the initial state $x_{0}$, hyperbolic-tangent activation
functions, and linear output functions. The total number of parameters is
$n_{\theta_{x}}+n_{\theta_{y}}+n_{\theta_{x0}}=296+225+128=649$, that we train
using NAILM111The resulting model can be retrieved at
http://cse.lab.imtlucca.it/~bemporad/shared/silverbox/rnn888.zip on $E=150$
epochs with $\rho_{\theta}=0.01$ and regularization
$\frac{0.1}{2}\|\theta_{x0}\|_{2}^{2}$ on the parameters defining $f_{x0}$.
For comparison, we consider the autoregressive models considered in [18], in
particular the ARX model ml and NLARX models ms, mlc, ms, and ms8c50, that we
trained on the same dataset by using the SYS-ID Toolbox [16] using the same
commands reported in [18]. Note that the latter two models contain explicitly
the nonlinear term $y_{k-1}^{3}$, in accordance with the physics-based model
of the Silverbox system, which is an electronic implementation of the Duffing
oscillator. The results of the comparison are shown in Table 3, which also
includes the results reported in the recent paper [3] and those obtained by
the LSTM model proposed in [17]. The results reproduced here differ from those
reported originally in the papers, that we show included in brackets (the
results originally reported in [17] have been omitted, as they were obtained
on a reduced subset of test and training data).

identification method | RMSE (mV) | BFR
---|---|---
ARX (ml) [18] | 16.29 [4.40] | 69.22 [73.79]
NLARX (ms) [18] |  08.42 [4.20] | 83.67 [92.06]
NLARX (mlc) [18] |  01.75 [1.70] | 96.67 [96.79]
NLARX (ms8c50) [18] |  01.05 [0.30] | 98.01 [99.43]
Recurrent LSTM model [17] |  02.20 [5.90] | 95.83 [88.90]
State-space encoder [3] ($n_{x}=4$) |  00.00 [1.40] |  00.00 [97.35]
NAILM (this paper) | 0.35 [0.00] | 99.33 [00.00]
Table 3: Silverbox benchmark: RMSE and BFR obtained by different methods. The
numbers reported in brackets refer to the results reported in the
corresponding original papers.

###  5.2 RNN training with $\ell_{1}$-regularization

We consider again the fluid damper benchmark dataset and add an
$\ell_{1}$-regularization term $g(\theta_{x},\theta_{y})$ as in (24) and use
NAILS to train the same RNN model structure (2b), for different values of
$\tau_{x}=\tau_{y}\triangleq\tau$. Figure 2 shows the resulting BRF. Figure 2
also shows the percentage of zero coefficients in the resulting parameter
vector $\left[\begin{smallmatrix}\theta_{x}\\\
\theta_{y}\end{smallmatrix}\right]\in{\mathbb{R}}^{73}$, i.e., its sparsity.

Figure 2: Fluid damper benchmark: BRF and sparsity of $\theta_{x},\theta_{y}$
for different $\ell_{1}$-regularization coefficients $\tau$ (mean results over
20 runs from different random initial weights)

The benefit of using $\ell_{1}$-regularization to reduce the complexity of the
model without sacrificing fit quality is clear. In fact, the BRF value on test
data remains roughly constant until $\tau\approx 0.1$, which corresponds to
$\approx 30\%$ zero coefficients in the model.

In order to compare NAILS/NAILM with other state-of-the-art training
algorithms that can deal with $\ell_{1}$-regularization terms (Adam [11],
AMSGrad [26], DiffGrad [6]), we report in Table 4 the results obtained by
solving the training problem for $\tau=0.2$ using the EKF approach [4] and
some of the SGD algorithms mostly used in machine-learning packages.
NAILS/NAILM is run for $E=250$, SGD for $E=2000$, and EKF for $E=50$ epochs.
While all methods get a similar BFR, NAILS and NAILM are more effective in
both sparsifying the model-parameter vector and execution time.

training | BFR | BFR | sparsity | CPU
---|---|---|---|---
algorithm | training | test | % | time (s)
NAILS | 91.00 (1.66) | 87.71 (2.67) | 65.1 (6.5) | 11.4
NAILM | 91.32 (1.19) | 87.80 (1.86) | 64.1 (7.4) | 11.7
EKF | 89.27 (1.48) | 86.67 (2.71) | 47.9 (9.1) | 13.2
AMSGrad | 91.04 (0.47) | 88.32 (0.80) | 16.8 (7.1) | 64.0
Adam | 90.47 (0.34) | 87.79 (0.44) |  08.3 (3.5) | 63.9
DiffGrad | 90.05 (0.64) | 87.34 (1.14) |  07.4 (4.5) | 63.9
Table 4: Fluid damper benchmark: mean BRF, sparsity of $[\theta_{x}^{\prime}\
\theta_{y}^{\prime}]^{\prime}$ (with standard deviation in parentheses), and
mean CPU time given by different solution methods when solving
$\ell_{1}$-regularized training for $\tau=0.2$ (results obtained over 20 runs
from different random initial weights)

For the Silverbox benchmark problem, following the analysis performed in [19],
Figure 3 shows the obtained RMSE on test data as a function of the number of
nonzero parameters in the model that are obtained by running NAILM to train
the same RNN model structure defined in Section 5.1.1 under an additional
$\ell_{1}$-regularization term with weight $\tau$ between $10^{-4}$ and 20.
The figure shows the effectiveness of NAILM in trading off model simplicity vs
quality of fit.

Figure 3: Silverbox benchmark: RMSE vs number of parameters. Models ml, ms,
mlc, and ms8c50 [18], LSTM model [17], and RNN models generated by NAILM for
values of the $\ell_{1}$-regularization parameter $\tau$ between 0 and 20.

###  5.3 RNN training with quantization

Consider again the fluid damper benchmark dataset and let $\theta_{x}$,
$\theta_{y}$ be only allowed to take values in the finite set ${\mathcal{Q}}$
of the multiples of $0.1$ between $-0.5$ and $0.5$. We also change the
activation function of the neurons to the leaky-ReLU function
$f(x)=\max\\{x,0\\}+0.1\min\\{x,0\\}$, so that the evaluation of the resulting
model amounts to extremely simple arithmetic operations, and increase the
number of hidden neurons to $n_{1}^{x}=n_{1}^{y}=6$. We train the model using
NAILS to handle the non-smooth and nonconvex regularization term
$g(\theta_{x},\theta_{y})=0$ if $\theta_{x}\in{\mathcal{Q}}^{n_{\theta_{x}}}$,
$\theta_{y}\in{\mathcal{Q}}^{n_{\theta_{y}}}$,
$g(\theta_{x},\theta_{y})=+\infty$ otherwise, with $E=1$, $N_{\rm ADMM}=200$,
$\rho=10$, without changing the remaining parameters. For comparison, we solve
the same problem without non-smooth regularization term $g$ by running
Algorithm 1 for $E=N_{\rm ADMM}$ epochs and then quantize the components of
the resulting vectors $\theta_{x}$, $\theta_{y}$ to their closest values in
${\mathcal{Q}}$.

The mean and standard deviation of the BRF and corresponding CPU time obtained
over 20 runs from different random initial weights is reported in Table 5. It
is apparent that quantizing a posteriori the parameters of a trained RNN leads
to much poorer results.

BFR | NAILS | quantization
---|---|---
training | 84.36 (3.00) | 17.64 (6.41)
test | 78.43 (4.38) | 12.79 (7.38)
CPU time | 12.04 (0.54) |  07.19 (2.82)
Table 5: Fluid damper benchmark with quantized model coefficients: mean BRF
and CPU time (standard deviation) over 20 runs

###  5.4 Group-Lasso penalty and model-order reduction

We use the group-Lasso penalty $\tau_{g}$ described in Section 4.2 to control
the model order $n_{x}$ when training a model on the fluid damper benchmark
dataset. We consider the same settings as in Section 5.1, except for the
model-order $n_{x}=8$, $n_{1}^{x}=n_{1}^{y}=6$ neurons in the hidden layers of
the state-update and output functions, and set the ADMM parameter
$\rho=10\tau_{g}$. The mean BRF results over 20 different runs of the NAILS
algorithm obtained for different values of $\tau_{g}$ are reported in Figure
4, along with the resulting mean model order obtained. From the figure, one
can see that $n_{x}=3$ is a good candidate model order to provide the best BRF
on test data.

Figure 4: BRF and resulting model order $n_{x}$ for different group-Lasso
regularization coefficients $\tau_{g}$ (mean results over 20 runs from
different random initial weights)

###  5.5 RNN training with smooth non-quadratic loss

To test the effectiveness of the proposed approach in handling non-quadratic
strongly convex and smooth loss terms, we consider $2000$ input/output pairs
generated by the following nonlinear system with binary outputs

$\begin{array}[]{rcl}x(k+1)&=&\left[\begin{smallmatrix}.8&.2&-.1\\\ 0&.9&.1\\\
.1&-.1&.7\end{smallmatrix}\right]\mathop{\rm
diag}\nolimits(x(k))(0.9+0.1\sin(x(k)))\\\\[5.0pt]
&&+\left[\begin{smallmatrix}-1\\\ .5\\\
1\end{smallmatrix}\right]u(k)(1-u^{3}(k))+\xi(k)\\\\[5.0pt]
y(k)&=&\left\\{\begin{array}[]{ll}1&\mbox{if}\
{\scriptsize\left[\begin{array}[]{c}-2\\\\[-2.10002pt] 1.5\\\\[-2.10002pt]
0.5\end{array}\right]}^{\prime}(x(k)+\frac{1}{3}x^{3}(k))-{\scriptsize{4}}+\zeta(k)\geq
0\\\ 0&\mbox{otherwise}\end{array}\right.\end{array}$

from $x(0)=0$, with the values of the input $u(k)$ changed with 90%
probability from step $k$ to $k+1$ with a new value drawn from the uniform
distribution on $[0,1]$. We consider independent noise signals
$\xi_{i}(k),\zeta(k)\sim{\mathcal{N}}(0,\sigma_{n}^{2})$, $i=1,2,3$, for
different values of the standard deviation $\sigma_{n}$. The first $N=1000$
samples are used for training, the rest for testing. NAILS and NAILM are used
to train a RNN model (2b) with no feedthrough, one hidden layer
($L_{x}=L_{y}=2$) with $n_{1}^{x}=n_{1}^{y}=5$ neurons and hyperbolic-tangent
activation function, and sigmoid output function
$f_{2}^{y}(y)=1/(1+e^{-A_{2}^{y}[x^{\prime}(k)\ u(k)]^{\prime}-b_{2}^{y}})$.
The resulting model-parameter vectors $\theta_{x}\in{\mathbb{R}}^{43}$ and
$\theta_{y}\in{\mathbb{R}}^{26}$ are trained using the modified cross-entropy
loss
$\ell(y(k),\hat{y})=\sum_{i=1}^{n_{y}}-y_{i}(k)\log(\epsilon+\hat{y}_{i})-(1-y_{i}(k))\log(1+\epsilon-\hat{y}_{i})$
[4], where $\hat{y}=f_{y}(x_{k},\theta_{y})$ and $\epsilon=10^{-4}$, quadratic
regularization with $\rho_{x}=0.1$, $\rho_{\theta}=0.01$, optimality tolerance
$\epsilon_{V}=10^{-6}$, parameters $c_{1}=10^{-4}$, $\sigma=0.5$,
$N_{\sigma}=10$ (NAILS), and $\lambda_{0}=100$, $c_{2}=1.5$, $c_{3}=5$,
$N_{\lambda}=30$ (NAILM), for maximum $N_{\rm ADMM}=150$ epochs.

Table 6 shows the final accuracy (%) achieved on training and test data for
different levels $\sigma_{n}$ of noise, averaged over 20 runs from different
initial values of the model parameters. The average CPU time to converge is
5.2 s (NAILS) and 6.2 s (NAILM).

$\sigma_{n}$ |  | accuracy (%) | accuracy (%)
---|---|---|---
|  | training data | test data
0.00 | NAILS | 99.16 (0.6) | 96.87 (1.0)
| NAILM | 99.30 (1.0) | 96.66 (1.3)
0.05 | NAILS | 98.45 (0.7) | 94.53 (1.2)
| NAILM | 98.13 (0.3) | 94.73 (1.9)
0.20 | NAILS | 86.03 (4.1) | 83.32 (5.5)
| NAILM | 86.33 (1.5) | 85.52 (1.7)
Table 6: Nonlinear system with binary outputs: mean final accuracy (%)
(standard deviation) over 20 runs on training and test data

##  6 Conclusions

We have proposed a training algorithm for RNNs that is computationally
efficient, provides very good quality solutions, and can handle rather general
loss and regularization terms. Although we focused on recurrent neural
networks of the form (2b), the algorithm is applicable to learn parametric
models with rather arbitrary structures, including other black-box state-space
models (such as LSTM models), gray-box and white-box models, and static
models.

Future research will be devoted to further analyze the convergence properties
of Algorithms 1 and 2 within the Generalized Gauss-Newton framework, and to
establish conditions on the ADMM parameter $\rho$, number of epochs $E$ used
when solving (LABEL:eq:ADMM-x), the model structure to learn, and the
regularization function $g$ that guarantee convergence of the nonconvex ADMM
iterations.

##  7 Acknowledgement

The author thanks M. Schoukens for exchanging ideas on the setup of the
Silverbox benchmark problem.

## References

  * [1] S.T. Alexander and A.L. Ghirnikar.  A method for recursive least squares filtering based upon an inverse QR decomposition.  IEEE Trans. Signal Processing, 41(1):20–30, 1993.
  * [2] J.A.E. Andersson, J. Gillis, G. Horn, J.B. Rawlings, and M. Diehl.  CasADi – A software framework for nonlinear optimization and optimal control.  Mathematical Programming Computation, 11(1):1–36, 2019.
  * [3] G. Beintema, R. Toth, and M. Schoukens.  Nonlinear state-space identification using deep encoder networks.  In Proc. Machine Learning Research, volume 144, pages 241–250, 2021\.
  * [4] A. Bemporad.  Recurrent neural network training with convex loss and regularization functions by extended Kalman filtering.  2021\.  Submitted for publication. Available on https://arxiv.org/abs/2111.02673.
  * [5] A. Bemporad and G. Cimini.  Variable elimination in model predictive control based on K-SVD and QR factorization.  IEEE Trans. Automatic Control, 2022.  In press. Also available on arXiv at http://arxiv.org/abs/2012.10423.
  * [6] S.H. Dubey, S. Chakraborty, S.K. Roy, S. Mukherjee, S.K. Singh, and B.B. Chaudhuri.  diffGrad: an optimization method for convolutional neural networks.  IEEE Transactions on Neural Networks and Learning Systems, 31(11):4500–4511, 2019.
  * [7] X. Glorot and Y. Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In Proc. 13${}^{t}h$ Int. Conf. Artificial Intelligence and Statistics, pages 249–256, 2010.
  * [8] M. Hong, Z.-Q. Luo, and M. Razaviyayn.  Convergence analysis of alternating direction method of multipliers for a family of nonconvex problems.  SIAM Journal on Optimization, 26(1):337–364, 2016.
  * [9] B. Houska, J. Frasch, and M. Diehl.  An augmented Lagrangian based algorithm for distributed nonconvex optimization.  SIAM Journal on Optimization, 26(2):1101–1127, 2016.
  * [10] K.J. Hunt, D. Sbarbaro, R. Żbikowski, and P.J. Gawthrop.  Neural networks for control systems – A survey.  Automatica, 28(6):1083–1112, 1992.
  * [11] D.P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [12] J.Z. Kolter and G. Manek.  Learning stable deep dynamics models.  Advances in neural information processing systems, 32, 2019.
  * [13] K. Levenberg.  A method for the solution of certain non-linear problems in least squares.  Quarterly of applied mathematics, 2(2):164–168, 1944.
  * [14] Z. Li, C. Ding, S. Wang, W. Wen, Y. Zhuo, C. Kiu, Q. Qiu, W. Xu, X. Lin, X. Qian, and Y. Wang.  E-RNN: Design optimization for efficient recurrent neural networks in FPGAs.  In IEEE Int. Symp. High Performance Computer Architecture, pages 69–80, 2019.
  * [15] M.M. Livstone, J.A. Farrell, and W.L. Baker.  A computationally efficient algorithm for training recurrent connectionist networks.  In Proc. American Control Conference, pages 555–561, 1992.
  * [16] L. Ljung.  System Identification Toolbox for MATLAB.  The Mathworks, Inc., 2001.  https://www.mathworks.com/help/ident.
  * [17] L. Ljung, C. Andersson, K. Tiels, and T.B. Schön.  Deep learning and system identification.  IFAC-PapersOnLine, 53(2):1175–1181, 2020.
  * [18] L. Ljung, Q. Zhang, P. Lindskog, and A. Juditski.  Estimation of grey box and black box models for non-linear circuit data.  IFAC Proceedings Volumes, 37(13):399–404, 2004.
  * [19] A. Marconato, M. Schoukens, Y. Rolain, and J. Schoukens.  Study of the effective number of parameters in nonlinear identification benchmarks.  In Proc. 52nd IEEE Conf. Dec. Contr., pages 4308–4313, 2013\.
  * [20] D.W. Marquardt.  An algorithm for least-squares estimation of nonlinear parameters.  Journal of the society for Industrial and Applied Mathematics, 11(2):431–441, 1963.
  * [21] D. Masti and A. Bemporad.  Learning nonlinear state-space models using autoencoders.  Automatica, 129:109666, 2021.
  * [22] F. Messerer, K. Baumgärtner, and M. Diehl.  Survey of sequential convex programming and generalized Gauss-Newton methods.  ESAIM. Proceedings and Surveys, 71:64, 2021.
  * [23] J. Nocedal and S.J. Wright.  Numerical Optimization.  Springer, 2nd edition, 2006.
  * [24] R. Parisi, E.D. Di Claudio, A. Rapagnetta, and G. Orlandi.  Recursive least squares approach to learning in recurrent neural networks.  In Proc. Int. Conf. on Neural Networks, volume 2, pages 1350–1354, 1996.
  * [25] G.V. Puskorius and L.A. Feldkamp.  Neurocontrol of nonlinear dynamical systems with Kalman filter trained recurrent networks.  IEEE Transactions on Neural Networks, 5(2):279–297, 1994.
  * [26] S.J. Reddi, S. Kale, and S. Kumar.  On the convergence of adam and beyond.  arXiv preprint arXiv:1904.09237, 2019.
  * [27] A.H. Ribeiro, K. Tiels, J. Umenberger, T.B. Schön, and L.A. Aguirre.  On the smoothness of nonlinear system identification.  Automatica, 121:109158, 2020.
  * [28] J.A.K. Suykens, J.P.L. Vandewalle, and B.L.R. De Moor.  Artificial neural networks for modelling and control of non-linear systems.  Springer Science & Business Media, 1995.
  * [29] R. Takapoui, N. Moehle, S. Boyd, and A. Bemporad.  A simple effective heuristic for embedded mixed-integer quadratic programming.  Int. J. Control, 93(1):2–12, 2020.
  * [30] Y. Tang, Z. Kan, D. Sun, L. Qiao, K. Xiao, Z. Lai, and D. Li.  ADMMiRNN: Training RNN with stable convergence via an efficient ADMM approach.  In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pages 3–18, 2020.
  * [31] G. Taylor, R. Burmeister, Z. Xu, B. Singh, A. Patel, and T. Goldstein.  Training neural networks without gradients: A scalable admm approach.  In International Conference on Machine Learning, pages 2722–2731, 2016.
  * [32] A. Themelis and P. Patrinos.  Douglas–Rachford splitting and ADMM for nonconvex optimization: Tight convergence results.  SIAM Journal on Optimization, 30(1):149–181, 2020.
  * [33] M.K. Transtrum and J.P. Sethna.  Improvements to the Levenberg-Marquardt algorithm for nonlinear least-squares minimization.  arXiv preprint arXiv:1201.5885, 2012.
  * [34] A.I.F. Vaz and L.N. Vicente.  PSwarm: A hybrid solver for linearly constrained global derivative-free optimization.  Optimization Methods and Software, 24:669–685, 2009.  http://www.norg.uminho.pt/aivaz/pswarm/.
  * [35] J. Wang, A. Sano, T. Chen, and B. Huang.  Identification of Hammerstein systems without explicit parameterisation of non-linearity.  International Journal of Control, 82(5):937–952, 2009.
  * [36] X. Wang and Y. Huang.  Convergence study in extended Kalman filter-based training of recurrent neural networks.  IEEE Transactions on Neural Networks, 22(4):588–600, 2011.
  * [37] Y. Wang, W. Yin, and J. Zeng.  Global convergence of ADMM in nonconvex nonsmooth optimization.  Journal of Scientific Computing, 78(1):29–63, 2019.
  * [38] P.J. Werbos.  Backpropagation through time: what it does and how to do it.  Proceedings of the IEEE, 78(10):1550–1560, 1990.
  * [39] T. Wigren and J. Schoukens.  Three free data sets for development and benchmarking in nonlinear system identification.  In European Control Conference, pages 2933–2938, Zürich, Switzerland, 2013.
  * [40] R.J. Williams and J. Peng.  An efficient gradient-based algorithm for on-line training of recurrent network trajectories.  Neural computation, 2(4):490–501, 1990.
  * [41] Q. Xu, K. Krishnamurthy, B. McMillin, and W. Lu.  A recursive least squares training algorithm for multilayer recurrent neural networks.  In Proc. American Control Conference, volume 2, pages 1712–1716, 1994.
  * [42] M. Yuan and Y. Lin.  Model selection and estimation in regression with grouped variables.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(1):49–67, 2006.