# A recursive parameter identification algorithm for nonlinear errors-in-variables models

Hugo Koide, Jérémy Vayssettes, Guillaume Mercère

► To cite this version:

Hugo Koide, Jérémy Vayssettes, Guillaume Mercère. A recursive parameter identification algorithm for nonlinear errors-in-variables models. IFAC Journal of Systems and Control, 2026, 35, pp.100381. 10.1016/j.ifacsc.2026.100381. hal-05511384

HAL Id: hal-05511384

https://hal.science/hal-05511384v1

Submitted on 15 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Full length article

# A recursive parameter identification algorithm for nonlinear errors-in-variables models☆,☆☆

Hugo Koide a,b, Jérémy Vayssettes b, Guillaume Mercère a,*

## ARTICLE INFO

**Article history:**

Received 27 October 2025

Received in revised form 30 December 2025

Accepted 27 January 2026

Available online 2 February 2026

**Keywords:**

Time varying system identification

Errors-in-variables

Recursive parameter identification

Nonlinear programming

## ABSTRACT

Recursive identification of model parameters from data is an essential tool in adaptive signal processing and adaptive control applications. For general nonlinear regression problems where errors are present in both dependent and independent variables, the errors-in-variables (EIV) model is frequently employed to avoid a biased estimation of the system parameters. This work presents a recursive identification algorithm for nonlinear EIV models which is derived from an offline sequential quadratic programming solution. A proof of convergence is provided and the performance of the proposed algorithm is illustrated on benchmark simulation examples. The results demonstrate the numerical efficiency of the recursive algorithm, as well as its ability to track time-varying system parameters.

© 2026 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license

(http://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Nonlinear model structures are ubiquitous in science and engineering applications, where the problem of estimating model parameters from data is often referred to as nonlinear regression or data fitting (Hansen et al., 2013; Seber & Wild, 2003). Assuming the number of data samples, N, is sufficiently large, the nonlinear regression problem boils down to solving an overdetermined system of nonlinear equations,

$$ \mathbf{y}_i \approx f(\mathbf{x}_i, \boldsymbol{\theta}^*), \quad i \in \{1, \dots, N\}, \qquad (1) $$

where $\mathbf{y}_i \in \mathbb{R}^{n_y}$ and $\mathbf{x}_i \in \mathbb{R}^{n_x}$ are the output (dependent) and input (independent) variables for a multiple-input multiple-output (MIMO) system, respectively, $\boldsymbol{\theta}^* \in \mathbb{R}^{n_\theta}$ is the vector of unknown true system parameters contained in the parameter space $\Theta \subset \mathbb{R}^{n_\theta}$, and where $\mathbf{f}: \mathbb{R}^{n_x} \times \mathbb{R}^{n_\theta} \to \mathbb{R}^{n_y}$ is a sufficiently smooth nonlinear function in $\boldsymbol{\theta} \in \Theta$. The inconsistencies in (1) (represented by an approximative relationship) typically arise from inaccuracies in the chosen model or due to errors in the

observed data, such that there is no exact solution for $\boldsymbol{\theta}^*$ that satisfies the equations in (1). In the well-known nonlinear least squares (NLS) approach (Björk, 1996), the parameters in (1) are estimated by considering the presence of additive errors on the output variables $\mathbf{y}_i$ exclusively, in what is generally referred to as an output error model (Ljung, 1999).

In many practical applications, however, the input variables $\mathbf{x}_i$ are measured, and may also be subject to error. It is well known that in the case of errors in both dependent and independent variables, the NLS approach is ill-advised, resulting in a biased estimation of the system parameters (Ljung, 1999). To address this issue, a useful extension of the output error model for general regression problems is the errors-in-variables (EIV) model,

$$ \begin{gathered} \mathbf{y}_i = f(\mathbf{x}_i, \mathbf{e}_{\mathbf{x}_i}, \boldsymbol{\theta}^*) + \mathbf{e}_{\mathbf{y}_i} \quad \text{with} \quad \mathbf{y}_i = \mathbf{y}_i^* + \mathbf{e}_{\mathbf{y}_i}, \\ \mathbf{x}_i = \mathbf{x}_i^* + \mathbf{e}_{\mathbf{x}_i}, \end{gathered} \qquad (2) $$

where $\mathbf{e}_{\mathbf{x}_i} \in \mathbb{R}^{n_x}$ and $\mathbf{e}_{\mathbf{y}_i} \in \mathbb{R}^{n_y}$ are the errors acting on the true data $\mathbf{x}_i^* \in \mathbb{R}^{n_x}$ and $\mathbf{y}_i^* \in \mathbb{R}^{n_y}$, respectively, and where $\mathbf{f}: \mathbb{R}^{n_x} \times \mathbb{R}^{n_x} \times \mathbb{R}^{n_\theta} \to \mathbb{R}^{n_y}$. EIV models have been extensively studied in the fields of system identification, statistics, econometrics, and signal processing, owing to the pervasiveness of errors in real-world data (Söderström, 2018; Wolter & Fuller, 1982).

Assuming access to the data samples $\{(\mathbf{x}_i, \mathbf{y}_i)\}_{i=1}^N$, the input/output data can be stacked in column-vector form

$$ \mathbf{x} = \begin{bmatrix} \mathbf{x}_1 \\ \vdots \\ \mathbf{x}_N \end{bmatrix} \in \mathbb{R}^{Nn_x}, \quad \mathbf{y} = \begin{bmatrix} \mathbf{y}_1 \\ \vdots \\ \mathbf{y}_N \end{bmatrix} \in \mathbb{R}^{Nn_y}. \qquad (3) $$

<footnote>a Université de Poitiers, ISAE-ENSMA Poitiers, LIAS, Poitiers, France</footnote>

<footnote>b Michelin Centre de Recherche de Ladoux, Rue Bleue, Cébazat, 63118, France</footnote>

<footnote>☆ This article is part of a Special issue entitled: 'IFAC WC 2026 - TC1.1' published in IFAC Journal of Systems and Control.</footnote>

<footnote>☆☆ This work was partially supported by I-TireLab (2025-2029), a joint public-private laboratory focused on modeling, estimation, and control for connected tires and rubber manufacturing processes. The laboratory is supported by Michelin, Université de Poitiers, Univ. Grenoble Alpes, France, Grenoble INP, and the CNRS.</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: hugo.koide@univ-poitiers.fr (H. Koide), jeremy.vayssettes@ michelin.com (J. Vayssettes), guillaume.mercere@univ-poitiers.fr (G. Mercère).</footnote>
The parameters in (2) can subsequently be estimated by solving
the following nonlinear total least squares (NTLS) problem,

$$
\min_{\theta, e_y, e_x} \frac{1}{2} \|e_y\|^2_{\mathbf{w}_y} + \frac{1}{2} \|e_x\|^2_{\mathbf{w}_x} \quad (4)
$$

subject to $\mathbf{y} = \mathbf{f}(\mathbf{x}, \mathbf{e}_x, \theta) + \mathbf{e}_y$,

where $\mathbf{W}_y = \text{diag}(\mathbf{W}_{y1}, \dots, \mathbf{W}_{yN}) \in \mathbb{R}^{N_{ny} \times N_{ny}}$ and $\mathbf{W}_x = \text{diag}(\mathbf{W}_{x1}, \dots, \mathbf{W}_{xN}) \in \mathbb{R}^{N_{nx} \times N_{nx}}$ are block diagonal matrices containing the symmetric positive definite weighting matrices $\{\mathbf{W}_{y1}\}_{i=1}^N$ and $\{\mathbf{W}_{x1}\}_{i=1}^N$, respectively, and where $\mathbf{f}: \mathbb{R}^{N_{nx}} \times \mathbb{R}^{N_{ny}} \times \mathbb{R}^{N_{\theta}} \to \mathbb{R}^{N_{ny}}$. In the nonlinear regression literature, the minimization problem (4) is also known under the rubric of orthogonal distance regression (ODR) (Boggs et al., 1987).

Offline solutions to the ODR problem date as far back as Ad-
cock (1878) for the case of **f** being linear in **θ**, and Deming (1943)
for nonlinear systems. In Schwetlick and Tiller (1985) and Boggs
et al. (1987), Gauss-Newton algorithms with trust region strate-
gies are developed for solving an unconstrained reformulation of
(4). In Tjoa and Biegler (1992) a decoupled sequential quadratic
programming (SQP) algorithm is developed for solving the con-
strained optimization problem in (4). A commonality between the
aforementioned offline solutions is that they are limited by the
additional degrees of freedom in e<sub>x</sub> ∈ ℝ<sup>N<sub>nx</sub></sup> and e<sub>y</sub> ∈ ℝ<sup>N<sub>ny</sub></sup>, mean-
ing that the estimation problem scales linearly in the number of
data samples N. Offline estimation of the parameters therefore
becomes prohibitively expensive for large datasets, and methods
to reduce the computational complexity of EIV estimation strate-
gies are consequently proposed in Schwetlick and Tiller (1989)
and Tjoa and Biegler (1992).

For adaptive signal processing and adaptive control applica-
tions, a recursive (online) identification algorithm for nonlinear
EIV models is therefore highly desirable, not only to process data
sequentially and allow tracking of time-varying parameters, but
also to avoid the dimensionality issues associated with offline
variants. Regarding linear EIV models, Davila (2002) and Rhode
et al. (2014) proposed recursive total least squares (TLS) algo-
rithms based on eigenvector tracking, and in Koide et al. (2024)
a recursive TLS algorithm is developed based on a constrained
optimization formulation of the linear EIV problem. To the best
of the authors' knowledge, however, no recursive identification
algorithm for general nonlinear EIV models has been proposed in
the literature.

To this end, this work presents a recursive nonlinear total least
squares algorithm which is derived from a SQP solution adapted
from Tjoa and Biegler (1992). A proof of convergence is provided
using the theoretical tools from Ljung (1977), Ljung and Söder-
ström (1983) and the performance of the proposed algorithm is
evaluated on benchmark simulation examples. The remainder of
this paper is organized as follows: Section 2 presents the offline
SQP algorithm for solving the NTLS problem, Section 3 derives
the recursive identification algorithm, including a variant for
efficient implementation, algorithm modifications, and a proof of
convergence, Section 4 illustrates the convergence characteristics,
efficiency, and tracking ability of the proposed algorithm through
simulation experiments, and Section 5 concludes this work.

2. Offline NTLS

Let us begin by noticing that the NTLS problem in (4) is
a quadratic optimization problem with nonlinear equality con-
straints and can therefore be solved by leveraging sequential
quadratic programming techniques which benefit from strong
local convergence properties (Nocedal & Wright, 2006; Tjoa &
Biegler, 1992).

With this in mind, let us define the following quantities

$$
\psi = [\theta^{\top} \quad \mathbf{e}_{x}^{\top} \quad \mathbf{e}_{y}^{\top}]^{\top} \in \mathbb{R}^{(n_{\theta} + N(n_{x} + n_{y}))}, \quad (5a)
$$

$$
h(\psi) = \frac{1}{2} \|\mathbf{e}_y\|^2_{\mathbf{w}_y} + \frac{1}{2} \|\mathbf{e}_x\|^2_{\mathbf{w}_x}, \quad (5b)
$$

$$
\mathbf{g}(\psi) = \mathbf{y} - \mathbf{e}_y - \mathbf{f}(\mathbf{x}, \mathbf{e}_x, \theta). \quad (5c)
$$

The NTLS problem can then be reformulated as the following
nonlinear equality constrained minimization problem

$$
\min_{\psi} h(\psi) \quad \text{such that} \quad \mathbf{g}(\psi) = \mathbf{0}, \tag{6}
$$

and the Lagrangian function associated with this minimization
problem takes the form

$$
\mathcal{L}(\psi, \lambda) = \frac{1}{2} \mathbf{e}_y^\top \mathbf{W}_y \mathbf{e}_y + \frac{1}{2} \mathbf{e}_x^\top \mathbf{W}_x \mathbf{e}_x + \lambda^\top (\mathbf{y} - \mathbf{e}_y - \mathbf{f}(\mathbf{x}, \mathbf{e}_x, \theta)), \quad (7)
$$

with $\lambda \in \mathbb{R}^{N_{ny}}$ a vector of Lagrange multipliers. Assuming access to $\{\hat{\psi}^{(j)}, \hat{\lambda}^{(j)}\}$ at iteration $j$, the SQP algorithm generates the next iterates $\{\hat{\psi}^{(j+1)}, \hat{\lambda}^{(j+1)}\}$ via the scheme

$$
\begin{bmatrix} \hat{\psi}^{(j+1)} \\ \hat{\lambda}^{(j+1)} \end{bmatrix} = \begin{bmatrix} \hat{\psi}^{(j)} \\ \hat{\lambda}^{(j)} \end{bmatrix} + \alpha^{(j)} \begin{bmatrix} \hat{\mathbf{p}}_{\psi}^{(j)} \\ \hat{\mathbf{\mu}}^{(j)} - \hat{\lambda}^{(j)} \end{bmatrix} \quad (8)
$$

where $\hat{\mathbf{p}}_{\psi}^{(j)} \in \mathbb{R}^{(n_{\theta}+N(n_{x}+n_{y}))}$ and $\hat{\mathbf{\mu}}^{(j)} \in \mathbb{R}^{N_{ny}}$ generate the search directions for the iterates $\hat{\psi}^{(j)}$ and $\hat{\lambda}^{(j)}$, respectively, and where $\alpha^{(j)} \in (0, 1]$ is a scalar step length. The SQP method then boils down to solving the following quadratic programming subproblem at each iteration

$$
\begin{bmatrix} \nabla_{\psi}^2 \mathcal{L}(\hat{\psi}^{(j)}, \hat{\lambda}^{(j)}) & \mathbf{J}(\hat{\psi}^{(j)}) \\ \mathbf{J}^{\top}(\hat{\psi}^{(j)}) & \mathbf{0} \end{bmatrix} \begin{bmatrix} \hat{\mathbf{p}}_{\psi}^{(j)} \\ \hat{\mathbf{\mu}}^{(j)} \end{bmatrix} = - \begin{bmatrix} \nabla_{\psi} h(\hat{\psi}^{(j)}) \\ \mathbf{g}(\hat{\psi}^{(j)}) \end{bmatrix}, \quad (9)
$$

where the gradient and Hessian of the Lagrangian function w.r.t. $\psi$
are given by

$$
\nabla_{\psi} \mathcal{L}(\psi, \lambda) = \begin{bmatrix} \mathbf{0} \\ \mathbf{W}_x \mathbf{e}_x \\ \mathbf{W}_y \mathbf{e}_y \end{bmatrix} + \begin{bmatrix} -\mathbf{A}^{\top} \lambda \\ -\mathbf{B}^{\top} \lambda \\ -\lambda \end{bmatrix}, \quad (10a)
$$

$$
\nabla_{\psi}^2 \mathcal{L}(\psi, \lambda) = \begin{bmatrix} \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{W}_x & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{W}_y \end{bmatrix} + \sum_{i=1}^{N_{ny}} \lambda_i \begin{bmatrix} -\nabla_{\theta}^2 f_i & -\nabla_{\theta e_x}^2 f_i & \mathbf{0} \\ -\nabla_{e_x \theta}^2 f_i & -\nabla_{e_x e_x}^2 f_i & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} \end{bmatrix}, \quad (10b)
$$

and where the Jacobians $\mathbf{A} \in \mathbb{R}^{N_{ny} \times n_{\theta}}$, $\mathbf{B} \in \mathbb{R}^{N_{ny} \times N_{nx}}$ of $\mathbf{f}$ are
defined as

$$
\mathbf{A} = \begin{bmatrix} \frac{\partial \mathbf{f}_1(\mathbf{x}, \mathbf{e}_x, \theta)}{\partial \theta^\top} \\ \vdots \\ \frac{\partial \mathbf{f}_N(\mathbf{x}, \mathbf{e}_x, \theta)}{\partial \theta^\top} \end{bmatrix}, \quad \mathbf{B} = \begin{bmatrix} \frac{\partial \mathbf{f}_1(\mathbf{x}, \mathbf{e}_x, \theta)}{\partial \mathbf{e}_x^\top} \\ \vdots \\ \frac{\partial \mathbf{f}_N(\mathbf{x}, \mathbf{e}_x, \theta)}{\partial \mathbf{e}_x^\top} \end{bmatrix}. \quad (11)
$$

The Jacobian of the constraints $\mathbf{g}(\psi)$ and the gradient of the
cost function $h(\psi)$ in (9) are then given by

$$
\mathbf{J}(\psi) = \begin{bmatrix} -\mathbf{A}^{\top} \\ -\mathbf{B}^{\top} \\ -\mathbf{I}_{N_{ny}} \end{bmatrix}, \quad \nabla_{\psi} h(\psi) = \begin{bmatrix} \mathbf{0} \\ \mathbf{W}_x \mathbf{e}_x \\ \mathbf{W}_y \mathbf{e}_y \end{bmatrix}, \quad (12)
$$

where $\mathbf{I}_{N_{ny}}$ is the $N_{ny} \times N_{ny}$ identity matrix. Let us then notice that, according to the first order necessary conditions for optimality (Nocedal & Wright, 2006),

$$
\nabla_{\psi} \mathcal{L}(\hat{\psi}^{(j)}, \hat{\lambda}^{(j)}) = \mathbf{0}, \quad (13)
$$

and thus from (10a) the Lagrange multipliers should satisfy $\hat{\lambda}^{(j)} = \mathbf{0}$. The quadratic subproblem at each iteration can then be reduced to solving the following system of equations (dropping the
iteration index for clarity)

$$ \begin{bmatrix} \mathbf{0} & \mathbf{0} & \mathbf{0} & -\mathbf{A}^{\top} \\ \mathbf{0} & \mathbf{W}_{\mathbf{x}} & \mathbf{0} & -\mathbf{B}^{\top} \\ \mathbf{0} & \mathbf{0} & \mathbf{W}_{\mathbf{y}} & -\mathbf{I}_{N_{\mathbf{ny}}} \\ -\mathbf{A} & -\mathbf{B} & -\mathbf{I}_{N_{\mathbf{ny}}} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \hat{\mathbf{p}}_{\theta} \\ \hat{\mathbf{p}}_{e_x} \\ \hat{\mathbf{p}}_{e_y} \\ \hat{\mathbf{\mu}} \end{bmatrix} = - \begin{bmatrix} \mathbf{0} \\ \mathbf{W}_{\mathbf{x}}\hat{\mathbf{e}}_{x} \\ \mathbf{W}_{\mathbf{y}}\hat{\mathbf{e}}_{y} \\ \mathbf{y} - \hat{\mathbf{y}} - \hat{\mathbf{e}}_{y} \end{bmatrix}, \quad (14) $$

where $\hat{\mathbf{y}} := \mathbf{f}(\mathbf{x}, \hat{\mathbf{e}}_x, \hat{\theta})$. Solving the second and third rows of (14) yields

$$ \hat{\mathbf{p}}_{e_x} = \mathbf{W}_{\mathbf{x}}^{-1}\mathbf{B}^{\top}\hat{\mathbf{\mu}} - \hat{\mathbf{e}}_{x}, \quad (15a) $$

$$ \hat{\mathbf{p}}_{e_y} = \mathbf{W}_{\mathbf{y}}^{-1}\hat{\mathbf{\mu}} - \hat{\mathbf{e}}_{y}, \quad (15b) $$

and from the fourth row of (14),

$$ \hat{\mathbf{\mu}} = \mathbf{M}^{-1}(\mathbf{y} - \hat{\mathbf{y}} + \mathbf{B}\hat{\mathbf{e}}_{\mathbf{x}} - \mathbf{A}\hat{\mathbf{p}}_{\theta}), \quad (16a) $$

$$ \mathbf{M} = \mathbf{W}_{\mathbf{y}}^{-1} + \mathbf{B}\mathbf{W}_{\mathbf{x}}^{-1}\mathbf{B}^{\top}. \quad (16b) $$

Finally, the first row of (14) gives the parameter search direction as

$$ \hat{\mathbf{p}}_{\theta} = (\mathbf{A}^{\top}\mathbf{M}^{-1}\mathbf{A})^{-1}\mathbf{A}^{\top}\mathbf{M}^{-1}(\mathbf{y} - \hat{\mathbf{y}} + \mathbf{B}\hat{\mathbf{e}}_{\mathbf{x}}). \quad (17) $$

The iterates at $j+1$ can then be updated according to (8).

## 2.1. Numerical complexity

A notable disadvantage of the offline NTLS algorithm described by (15)-(17) is that it suffers from the curse of dimensionality. Since the dimension of the optimization variable $\psi$ scales linearly with the number of data points $N$, the algorithm requires $N$-dimensional matrix operations, which is computationally expensive for large datasets.

Computing the matrix inverses in (15)-(17) can be done in an efficient manner by leveraging the block-diagonal structure of the weighting matrices $\mathbf{W}_y$ and $\mathbf{W}_x$ and of the Jacobian matrix $\mathbf{B}$. Indeed, we have that

$$ \mathbf{B} = \begin{bmatrix} \mathbf{B}_1 & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{B}_2 & \cdots & \mathbf{0} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{B}_N \end{bmatrix}, \quad (18) $$

with $\mathbf{B}_i = \partial \mathbf{f}_i(\mathbf{x}_i, \mathbf{e}_{x_i}, \theta) / \partial \mathbf{e}_{x_i} \in \mathbb{R}^{n_y \times n_x}$ for $i \in \{1, \dots, N\}$. This further leads to a block-diagonal structure for the matrix $\mathbf{M}$, i.e.,

$$ \mathbf{M} = \begin{bmatrix} \mathbf{M}_1 & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{M}_2 & \cdots & \mathbf{0} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{M}_N \end{bmatrix}, \quad (19) $$

with $\mathbf{M}_i = \mathbf{W}_{\mathbf{y}_i}^{-1} + \mathbf{B}_i \mathbf{W}_{\mathbf{x}_i}^{-1} \mathbf{B}_i^{\top} \in \mathbb{R}^{n_y \times n_y}$ for $i \in \{1, \dots, N\}$. Computing the inverse of $\mathbf{W}_y$, $\mathbf{W}_x$, and $\mathbf{M}$ then only requires the inversion of $N n_y \times n_y$ and $n_x \times n_x$ matrices, reducing the computational complexity from $\mathcal{O}(N^3 \max(n_y^3, n_x^3))$ to $\mathcal{O}(N \max(n_y^3, n_x^3))$ operations. Note that the final inversion of the matrix $(\mathbf{A}^{\top}\mathbf{M}^{-1}\mathbf{A})$ has complexity $\mathcal{O}(n_y^3)$, which is cheap to compute.

The algorithm remains limited, however, by the $N$-dimensional matrix multiplications in (15)-(17) which remain of computational complexity $\mathcal{O}(N^3)$. A recursive solution description then becomes important for memory-scarce and real-time applications with limited computational resources.

## 3. Online NTLS

Recursive identification algorithms provide a computationally efficient way of inferring the model parameters online, as new data becomes available, which is essential in adaptive filtering

and adaptive control applications which may require tracking of time-varying parameters. In the following section, a recursive algorithm is derived from the offline SQP solution developed in Section 2. An efficient implementation strategy is proposed along with important algorithm modifications for practical implementation purposes, and a proof of convergence is provided based on the analysis of Ljung (1977) for general recursive identification algorithms.

### 3.1. Recursive identification algorithm

Let us begin by noticing that the parameter search direction (17) can be expressed as

$$ \hat{\mathbf{p}}_{\theta_k} = \left[ \sum_{i=1}^{k} \mathbf{A}_i^{\top} \mathbf{M}_i^{-1} \mathbf{A}_i \right]^{-1} \sum_{i=1}^{k} \mathbf{A}_i^{\top} \mathbf{M}_i^{-1} (\mathbf{y}_i - \hat{\mathbf{y}}_i + \mathbf{B}_i \hat{\mathbf{e}}_{x_i}), \quad (20) $$

where data are accumulated up until the discrete index $k$. The weighted information matrix can then be defined as

$$ \mathbf{R}_k = \sum_{i=1}^{k} \mathbf{A}_i^{\top} \mathbf{M}_i^{-1} \mathbf{A}_i \quad (21) $$

such that

$$ \sum_{i=1}^{k-1} \mathbf{A}_i^{\top} \mathbf{M}_i^{-1} (\mathbf{y}_i - \hat{\mathbf{y}}_i + \mathbf{B}_i \hat{\mathbf{e}}_{x_i}) = \mathbf{R}_{k-1} \hat{\mathbf{p}}_{\theta_{k-1}}. \quad (22) $$

and where

$$ \mathbf{R}_{k-1} = \mathbf{R}_k - \mathbf{A}_k^{\top} \mathbf{M}_k^{-1} \mathbf{A}_k. \quad (23) $$

The parameter search direction can then be rearranged as

$$ \hat{\mathbf{p}}_{\theta_k} = \mathbf{R}_k^{-1} \left[ \sum_{i=1}^{k-1} \mathbf{A}_i^{\top} \mathbf{M}_i^{-1} (\mathbf{y}_i - \hat{\mathbf{y}}_i + \mathbf{B}_i \hat{\mathbf{e}}_{x_i}) + \mathbf{A}_k^{\top} \mathbf{M}_k^{-1} (\mathbf{y}_k - \hat{\mathbf{y}}_k + \mathbf{B}_k \hat{\mathbf{e}}_{x_k}) \right], \quad (24) $$

which can be combined with (22) and (23) such that

$$ \hat{\mathbf{p}}_{\theta_k} = \mathbf{R}_k^{-1} \left[ (\mathbf{R}_k - \mathbf{A}_k^{\top} \mathbf{M}_k^{-1} \mathbf{A}_k) \hat{\mathbf{p}}_{\theta_{k-1}} + \mathbf{A}_k^{\top} \mathbf{M}_k^{-1} (\mathbf{y}_k - \hat{\mathbf{y}}_k + \mathbf{B}_k \hat{\mathbf{e}}_{x_k}) \right], \quad (25) $$

and thus

$$ \hat{\mathbf{p}}_{\theta_k} = \hat{\mathbf{p}}_{\theta_{k-1}} + \mathbf{R}_k^{-1} \mathbf{A}_k^{\top} \mathbf{M}_k^{-1} (\mathbf{y}_k - \hat{\mathbf{y}}_k + \mathbf{B}_k \hat{\mathbf{e}}_{x_k} - \mathbf{A}_k \hat{\mathbf{p}}_{\theta_{k-1}}). \quad (26) $$

In order to render the update (26) fully recursive, let us assume a single optimization step is performed each time a new data sample $\{\mathbf{x}_k, \mathbf{y}_k\}$ becomes available. That is, for each $k \in \mathbb{Z}$, the optimization algorithm is initialized at $\hat{\theta}_{k-1}$, and a single iteration (26) is performed such that $\hat{\mathbf{p}}_{\theta_{k-1}} = \mathbf{0}$. The full recursive algorithm is then given by

$$ \zeta_k = \mathbf{y}_k - \hat{\mathbf{y}}_k + \mathbf{B}_k \hat{\mathbf{e}}_{x_{k-1}}, \quad (27a) $$

$$ \hat{\mathbf{p}}_{\theta_k} = \gamma_k \mathbf{R}_k^{-1} \mathbf{A}_k^{\top} \mathbf{M}_k^{-1} \zeta_k, \quad (27b) $$

$$ \hat{\mathbf{\mu}}_k = \mathbf{M}_k^{-1}(\mathbf{y}_k - \hat{\mathbf{y}}_k + \mathbf{B}_k \hat{\mathbf{e}}_{x_{k-1}} - \mathbf{A}_k \hat{\mathbf{p}}_{\theta_k}), \quad (27c) $$

$$ \hat{\mathbf{p}}_{e_{x_k}} = \mathbf{W}_{x_k}^{-1} \mathbf{B}_k^{\top} \hat{\mathbf{\mu}}_k - \hat{\mathbf{e}}_{x_{k-1}}, \quad (27d) $$

$$ \hat{\mathbf{p}}_{e_{y_k}} = \mathbf{W}_{y_k}^{-1} \hat{\mathbf{\mu}}_k - \hat{\mathbf{e}}_{y_{k-1}}, \quad (27e) $$

$$ \begin{bmatrix} \hat{\theta}_k \\ \hat{\mathbf{e}}_{x_k} \\ \hat{\mathbf{e}}_{y_k} \\ \hat{\mathbf{\lambda}}_k \end{bmatrix} = \begin{bmatrix} \hat{\theta}_{k-1} \\ \hat{\mathbf{e}}_{x_{k-1}} \\ \hat{\mathbf{e}}_{y_{k-1}} \\ \hat{\mathbf{\lambda}}_{k-1} \end{bmatrix} + \alpha_k \begin{bmatrix} \hat{\mathbf{p}}_{\theta_k} \\ \hat{\mathbf{p}}_{e_{x_k}} \\ \hat{\mathbf{p}}_{e_{y_k}} \\ \hat{\mathbf{\lambda}}_{k-1} \end{bmatrix}, \quad (27f) $$
<table><caption>**Table 1**
Algorithm floating-point operations.</caption><thead><tr><th>Algorithm</th><th>FLOPS</th></tr></thead><tbody><tr><td>NTLS (Offline)</td><td>N<sup>3</sup>[4n<sub>θ</sub><sup>2</sup>n<sub>y</sub> + 4n<sub>y</sub><sup>2</sup>n<sub>θ</sub> + 2n<sub>x</sub><sup>2</sup>n<sub>x</sub> + 2n<sub>x</sub><sup>2</sup>n<sub>y</sub>]<br>+N<sup>2</sup>[4n<sub>y</sub><sup>2</sup>+4n<sub>x</sub>n<sub>y</sub>+2n<sub>x</sub><sup>2</sup>]+N[n<sub>y</sub><sup>3</sup>+2n<sub>y</sub>n<sub>θ</sub>]</td></tr><tr><td>NTLS (Online)</td><td>2n<sub>y</sub><sup>3</sup>+n<sub>y</sub><sup>2</sup>+2n<sub>x</sub><sup>2</sup>n<sub>y</sub>+n<sub>y</sub><sup>2</sup>(4n<sub>x</sub>+4n<sub>θ</sub>)+8n<sub>y</sub><sup>2</sup>n<sub>y</sub></td></tr></tbody></table>

where $\gamma_k$ is introduced as a sequence of positive scalars tending to zero (a gain typically chosen to be 1/k) and where $\hat{y}_k = f(x_k, \hat{e}_{x_{k-1}}, \hat{\theta}_{k-1})$. The information matrix can then be updated according to (Ljung & Söderström, 1983)

$$ \mathbf{R}_k = \mathbf{R}_{k-1} + \gamma_k (\mathbf{A}_k^\top \mathbf{M}_k^{-1} \mathbf{A}_k - \mathbf{R}_{k-1}), \quad (28) $$

A standard way to improve the numerical efficiency of (27) is to use the Woodbury matrix inversion lemma (Zhang, 2017) to update the inverse of the information matrix directly, such that (28) becomes

$$ \mathbf{P}_k = \frac{1}{\kappa_k} (\mathbf{P}_{k-1} - \mathbf{L}_k \mathbf{A}_k \mathbf{P}_{k-1}), \quad (29a) $$

$$ \mathbf{L}_k = \mathbf{P}_{k-1} \mathbf{A}_k^\top (\kappa_k \mathbf{M}_k + \mathbf{A}_k \mathbf{P}_{k-1} \mathbf{A}_k^\top)^{-1}, \quad (29b) $$

where we introduce the covariance matrix $\mathbf{P}_k := \gamma_k \mathbf{R}_k^{-1}$ and the estimator gain $\mathbf{L}_k \in \mathbb{R}^{n_\theta \times n_y}$, and where $\kappa_k \in (0, 1]$ is a forgetting factor that serves to discount past data, and is related to the gain sequence $\gamma_k$ via

$$ \kappa_k = \frac{\gamma_{k-1}}{\gamma_k} (1 - \gamma_k). \quad (30) $$

Note that the choosing the gain $\gamma_k = 1/k$ is equivalent to setting $\kappa_k = 1$ and that setting $\kappa_k < 1$ allows for tracking of time-varying parameters via exponential forgetting of past data. The parameters can then be updated as

$$ \hat{\theta}_k = \hat{\theta}_{k-1} + \alpha_k \mathbf{L}_k \zeta_k. \quad (31) $$

The computational complexities of the offline NTLS and online NTLS algorithms for a single iteration/recursion are summarized in Table 1, where the highest order floating point operations (FLOPS) are reported. Note that the FLOP count for the offline algorithm accounts for the efficient matrix inversions discussed in Section 2. The recursive NTLS algorithm has overall complexity $O(\max(n_x^3, n_y^3))$, whereas the batch NTLS algorithm has overall complexity $O(N^3 \max(n_x^2, n_y^2, n_\theta^2))$. It is straightforward to see that $N \gg \max(n_x, n_y)$ for large sample sizes, such that the recursive formulation becomes significantly more efficient.

An important modification to improve the convergence properties of recursive identification algorithms is to consider the case where the Jacobian $\mathbf{A}_k$ is not full rank (e.g. due to lack of excitation in the data). For such ill-conditioned problems, the information matrix may become quasi-singular, and its inversion will be numerically unstable. To avoid this issue, a Levenberg-Marquardt (LM) regularization term can be added to the information matrix to ensure it is always invertible (Ljung & Söderström, 1983). The LM modification

$$ \bar{\mathbf{R}}_k = \mathbf{R}_{k-1} + \gamma_k (\mathbf{A}_k^\top \mathbf{M}_k^{-1} \mathbf{A}_k - \mathbf{R}_{k-1}), \quad (32a) $$

$$ \begin{aligned} \mathbf{R}_k &= \bar{\mathbf{R}}_k + \delta_k \mathbf{I}_{n_\theta}, \\ &= (\mathbf{I}_{n_\theta} + \delta_k \bar{\mathbf{R}}_k^{-1}) \bar{\mathbf{R}}_k, \end{aligned} \quad (32b) $$

where $\delta_k$ is a small positive scalar, can then be used, where the inversion of (32b) leads to

$$ \mathbf{P}_k = \bar{\mathbf{P}}_k (\mathbf{I}_{n_\theta} + \delta_k \bar{\mathbf{P}}_k)^{-1}, \quad (33) $$

and where (33) can be used as an additional update step after (29a).

Finally, to ensure convergence of the algorithm, it is essential to project the parameter estimates into a stability region $D_M \subset \mathbb{R}^{n_\theta}$ consistent with the model set $M$ such that

$$ D_M = \{ \theta \in \mathbb{R}^{n_\theta} | \mathcal{M}(\theta) \text{ is stable} \}. \quad (34) $$

A heuristic projection of the form

$$ \hat{\theta}_k = [\hat{\theta}_{k-1} + \alpha_k \mathbf{L}_k \zeta_k]_{D_M}, \quad (35) $$

$$ [b]_{D_M} = \begin{cases} b, & \text{if } b \in D_M, \\ \text{a value interior to } D_M, & \text{if } b \notin D_M, \end{cases} \quad (36) $$

can thus be included as a robust alternative to (31). Since the algorithm (27) is nothing more than a recursive implementation of a nonlinear optimization algorithm, the problem of converging to local minima is naturally inherited. If sufficient prior knowledge of the parameter space is available, global convergence can be ensured by initializing the algorithm within the domain of attraction of the global minimum.

## 3.2. Convergence analysis

The convergence analysis is based on the method for the analysis of recursive stochastic approximation algorithms developed in Ljung (1977), which was later adapted to more general recursive identification algorithms in Ljung and Söderström (1983). The analysis relies on relating the asymptotic trajectories of the recursive algorithm to the solutions of an associated ordinary differential equation (ODE) that would be obtained if the sequence $\{\hat{\theta}_k\}$ were replaced by a constant parameter $\theta$.

In Ljung and Söderström (1983), the authors adapted the analysis from Ljung (1977) to prove the convergence for a general class of recursive Gauss-Newton (GN) algorithms. The general GN algorithm is one which considers the quadratic criterion

$$ \ell(\theta_k) = \frac{1}{2} \varepsilon_k(\theta_k)^\top \mathbf{\Lambda}_k^{-1} \varepsilon_k(\theta_k), \quad (37) $$

where $\varepsilon_k(\theta_k) := \hat{y}_k - \hat{y}_k(\theta_k)$, where $\hat{y}_k(\theta_k)$ is the model prediction, and where $\mathbf{\Lambda}_k$ is a positive definite weighting matrix. The GN search direction $\hat{\mathbf{p}}_{\theta_k}^{\text{GN}}$ at time $k$ is then found as the solution to the least squares problem

$$ \min_{\mathbf{p}_{\theta_k}} \frac{1}{2} \left\| \mathbf{J}_k(\hat{\theta}_k) \mathbf{p}_{\theta_k} + \varepsilon_k(\hat{\theta}_k) \right\|_{\mathbf{\Lambda}_k^{-1}}^2, \quad (38) $$

where $\mathbf{J}_k(\hat{\theta}_k)$ is the Jacobian of the prediction error with respect to the parameters. The resulting recursion for the general GN algorithm studied in Ljung and Söderström (1983, Chapter 4) is then given by

$$ \hat{\mathbf{p}}_{\theta_k}^{\text{GN}} = \gamma_k \mathbf{R}_k^{-1} \mathbf{J}_k \mathbf{\Lambda}_k^{-1} \varepsilon_k, \quad (39a) $$

$$ \mathbf{R}_k = \mathbf{R}_{k-1} + \gamma_k (\mathbf{J}_k^\top \mathbf{\Lambda}_k^{-1} \mathbf{J}_k - \mathbf{R}_{k-1}), \quad (39b) $$

$$ \hat{\theta}_k = \hat{\theta}_{k-1} + \hat{\mathbf{p}}_{\theta_k}^{\text{GN}}, \quad (39c) $$

where $\gamma_k$ is, as before, a sequence of positive scalars tending to zero. Notice the resemblance between (39) and the recursions in the NTLS algorithm (27). This is no coincidence, since the recursive NTLS algorithm (27) developed from the SQP method is found to coincide exactly with the GN solution for the criterion function

$$ \ell(\boldsymbol{\phi}_k) = \frac{1}{2} \boldsymbol{g}_k(\boldsymbol{\phi}_k)^{\top} \mathbf{W}_k \boldsymbol{g}_k(\boldsymbol{\phi}_k), \quad (40) $$

where $\boldsymbol{\phi}_k = [\boldsymbol{\theta}_k^\top \mathbf{e}_{x_k}^\top]^\top$, and where

$$ \boldsymbol{g}_k(\boldsymbol{\phi}_k) = \left[ \mathbf{y}_k - \mathbf{f}_k(\mathbf{x}_k, \mathbf{e}_{x_k}, \boldsymbol{\theta}_k) \quad \mathbf{e}_{x_k} \right]. \quad (41) $$

This result is formalized in the following Lemma.
**Lemma 1.** Consider the nonlinear least squares problem

$$ \min_{\theta_k, e_{x_k}} \frac{1}{2} \left\| \left[ \mathbf{y}_k - \mathbf{f}_k(\mathbf{x}_k, \mathbf{e}_{x_k}, \theta_k) \right] \right\|_{\mathbf{w}_k}^2 = \min_{\phi_k} \frac{1}{2} \left\| \mathbf{g}_k(\phi_k) \right\|_{\mathbf{w}_k}^2, \quad (42) $$

where $\phi_k = [\theta_k^\top \ e_{x_k}^\top]^\top \in \mathbb{R}^{n_\theta+n_x}$, $g_k: \mathbb{R}^{n_\theta+n_x} \to \mathbb{R}^{n_y+n_x}$, and where

$$ \mathbf{W}_k = \begin{bmatrix} \mathbf{W}_{y_k} & \mathbf{0} \\ \mathbf{0} & \mathbf{W}_{x_k} \end{bmatrix} \quad (43) $$

is composed of the symmetric positive definite weighting matrices $\mathbf{W}_{y_k} \in \mathbb{R}^{n_y \times n_y}$ and $\mathbf{W}_{x_k} \in \mathbb{R}^{n_x \times n_x}$, respectively. Let us define the Jacobian matrix

$$ \mathbf{J}_k(\boldsymbol{\phi}_k) = \frac{\partial \mathbf{g}_k(\boldsymbol{\phi}_k)}{\partial \boldsymbol{\phi}_k^\top} = \begin{bmatrix} -\mathbf{A}_k & -\mathbf{B}_k \\ \mathbf{0} & \mathbf{I}_{n_x} \end{bmatrix} \quad (44) $$

where

$$ \mathbf{A}_k = \frac{\partial \mathbf{f}_k(\mathbf{x}_k, \mathbf{e}_{x_k}, \theta_k)}{\partial \theta_k^\top}, \quad \mathbf{B}_k = \frac{\partial \mathbf{f}_k(\mathbf{x}_k, \mathbf{e}_{x_k}, \theta_k)}{\partial \mathbf{e}_{x_k}^\top}. \quad (45) $$

Then, the Gauss-Newton search direction $\hat{\mathbf{p}}_{\phi_k}$ for the problem (42) at $\hat{\mathbf{p}}_k$ is the solution to the least squares problem

$$ \min_{\mathbf{p}_{\phi_k}} \frac{1}{2} \left\| \mathbf{J}_k(\hat{\mathbf{p}}_{\phi_k}) \mathbf{p}_{\phi_k} + \mathbf{g}_k(\hat{\mathbf{p}}_{\phi_k}) \right\|_{\mathbf{w}_k}^2, \quad (46) $$

where $\mathbf{p}_{\phi_k} = [\mathbf{p}_{\theta_k}^\top \ \mathbf{p}_{e_{x_k}}^\top]^\top$ and is given by

$$ \hat{\mathbf{p}}_{\theta_k} = (\mathbf{A}_k^\top \mathbf{M}_k^{-1} \mathbf{A}_k)^{-1} \mathbf{A}_k^\top \mathbf{M}_k^{-1} (\mathbf{y}_k - \hat{\mathbf{y}}_k + \mathbf{B}_k \hat{\mathbf{e}}_{x_k}), \quad (47) $$

$$ \hat{\mathbf{p}}_{e_{x_k}} = \mathbf{W}_{x_k}^{-1} \mathbf{B}_k^\top \mathbf{M}_k^{-1} (\mathbf{y}_k - \hat{\mathbf{y}}_k + \mathbf{B}_k \hat{\mathbf{e}}_{x_k} - \mathbf{A}_k \hat{\mathbf{p}}_{\theta_k}) - \hat{\mathbf{e}}_{x_k}, \quad (48) $$

**Proof.** See Appendix.

By observing that the resulting GN search directions (47)-(48) coincide with those generated by the SQP algorithm (15)-(17) with the chosen step length $\alpha_k = 1$, it follows that the convergence analysis of the recursive NTLS algorithm can be carried out using the results from Ljung and Söderström (1983) for the convergence of general recursive GN algorithms.

Let us then consider the following recursive NTLS algorithm,

$$ \zeta_k = \mathbf{y}_k - \hat{\mathbf{y}}_k + \mathbf{B}_k \hat{\mathbf{e}}_{x_k-1}, \quad (49a) $$

$$ \mathbf{M}_k = \mathbf{W}_{y_k}^{-1} + \mathbf{B}_k \mathbf{W}_{x_k}^{-1} \mathbf{B}_k^\top, \quad (49b) $$

$$ \bar{\mathbf{R}}_k = \mathbf{R}_{k-1} + \gamma_k (\mathbf{A}_k^\top \mathbf{M}_k^{-1} \mathbf{A}_k - \mathbf{R}_{k-1}), \quad (49c) $$

$$ \mathbf{R}_k = \begin{cases} \bar{\mathbf{R}}_k, & \text{if } \bar{\mathbf{R}}_k \ge \delta \mathbf{I}_{n_\theta}, \\ \bar{\mathbf{R}}_k + \delta_k \mathbf{I}_{n_\theta}, & \text{otherwise,} \end{cases} \quad (49d) $$

$$ \hat{\mathbf{p}}_{\theta_k} = \gamma_k \mathbf{R}_k^{-1} \mathbf{A}_k^\top \mathbf{M}_k^{-1} \zeta_k, \quad (49e) $$

$$ \hat{\mathbf{u}}_k = \mathbf{M}_k^{-1}(\mathbf{y}_k - \hat{\mathbf{y}}_k + \mathbf{B}_k \hat{\mathbf{e}}_{x_k-1} - \mathbf{A}_k \hat{\mathbf{p}}_{\theta_k}), \quad (49f) $$

$$ \hat{\mathbf{p}}_{e_{x_k}} = \mathbf{W}_{x_k}^{-1} \mathbf{B}_k^\top \hat{\mathbf{u}}_k - \hat{\mathbf{e}}_{x_k-1}, \quad (49g) $$

$$ \hat{\mathbf{p}}_{e_{y_k}} = \mathbf{W}_{y_k}^{-1} \hat{\mathbf{u}}_k - \hat{\mathbf{e}}_{y_k-1}, \quad (49h) $$

$$ \begin{bmatrix} \hat{\mathbf{\theta}}_k \\ \hat{\mathbf{e}}_{x_k} \\ \hat{\mathbf{e}}_{y_k} \\ \hat{\mathbf{\lambda}}_k \end{bmatrix} = \begin{bmatrix} \hat{\mathbf{\theta}}_{k-1} \\ \hat{\mathbf{e}}_{x_k-1} \\ \hat{\mathbf{e}}_{y_k-1} \\ \hat{\mathbf{\lambda}}_{k-1} \end{bmatrix} + \alpha_k \begin{bmatrix} \hat{\mathbf{p}}_{\theta_k} \\ \hat{\mathbf{p}}_{e_{x_k}} \\ \hat{\mathbf{p}}_{e_{y_k}} \\ \hat{\mathbf{\mu}}_k - \hat{\mathbf{\lambda}}_{k-1} \end{bmatrix} \quad (49i) $$

$$ \xi_{k+1} = A(\hat{\mathbf{\psi}}_k)\xi_k + B(\hat{\mathbf{\psi}}_k)z_k, \quad (49j) $$

$$ \begin{bmatrix} \hat{\mathbf{y}}_{k+1} \\ \mathrm{vec}(\mathbf{J}_{k+1}^\top) \end{bmatrix} = C(\hat{\mathbf{\psi}}_k)\xi_{k+1}, \quad (49k) $$

where $\psi_k = [\theta_k^\top \ e_{x_k}^\top \ e_{y_k}^\top]^\top$, where a LM regularization step is included in (49d) and where $z_k = [\mathbf{x}_k^\top \ \mathbf{y}_k^\top]^\top$ is the vector of input-output data. Eqs. (49j)-(49k) imply a (non-restrictive) condition that the predictions $\{\hat{\mathbf{y}}_k\}$ and the Jacobians $\{\mathbf{J}_k\}$ can be obtained

by filtering the data $\{z_k\}$ through a finite-dimensional linear filter governed by the parameter-dependent matrices $A(\psi)$, $B(\psi)$, and $C(\psi)$ via the internal state vector $\xi_k$. The convergence properties of (49) can then be analyzed by considering the following regularity conditions, where the compact notation $| \cdot |$ denotes a fixed matrix/vector norm.

**M1.** $D_M$ is a compact set of $\mathbb{R}^{n_\theta+n_x+n_y}$ such that $\psi \in D_M \Rightarrow A(\psi)$ has all eigenvalues strictly inside the unit circle.

**M2.** The matrices $A(\psi)$, $B(\psi)$, and $C(\psi)$ are continuously differentiable with respect to $\psi$ for $\psi \in D_M$.

**Cr1.** The function $\mathbf{H}_k(\mathbf{R}, \psi, \mathbf{g}, \mathbf{J}) = \mathbf{A}_k^\top \mathbf{M}_k^{-1} \mathbf{A}_k - \mathbf{R}_{k-1}$ is differentiable w.r.t. $\mathbf{R}, \psi, \mathbf{g}$, and $\mathbf{J}$, such that, for some $C < \infty$,

$$ \begin{aligned} |\mathbf{H}_k(\mathbf{R}, \psi, \mathbf{g}, \mathbf{J})| &\le C(1 + |\mathbf{g}|^2 + |\mathbf{J}|^2 + |\mathbf{R}|), \\ |\mathbf{H}_{k,R}(\mathbf{R}, \psi, \mathbf{g}, \mathbf{J})| &+ |\mathbf{H}_{k,\psi}(\psi, \mathbf{g}, \mathbf{J})| \le C(1 + |\mathbf{g}|^2 + |\mathbf{J}|^2), \\ |\mathbf{H}_{k,g}(\mathbf{R}, \psi, \mathbf{g}, \mathbf{J})| &+ |\mathbf{H}_{k,J}(\mathbf{R}, \psi, \mathbf{g}, \mathbf{J})| \le C(1 + |\mathbf{g}| + |\mathbf{J}|) \end{aligned} $$

for $\psi \in D_M$.

**Cr2.** The function $\ell_k(\psi, \mathbf{g}) = \frac{1}{2}\mathbf{g}_k(\psi)^\top \mathbf{W}_k \mathbf{g}_k(\psi)$ is twice continuously differentiable w.r.t. $\psi$ and $\mathbf{g}$ such that, for some $C < \infty$,

$$ \begin{aligned} |\ell_{k,\psi}(\psi, \mathbf{g})| + |\ell_{k,\psi\psi}(\psi, \mathbf{g})| &\le C(1 + |\mathbf{g}|)^2, \\ |\ell_{k,g}(\psi, \mathbf{g})| + |\ell_{k,g\psi}(\psi, \mathbf{g})| &\le C(1 + |\mathbf{g}|), \\ |\ell_{k,gg}(\psi, \mathbf{g})| &\le C, \end{aligned} $$

for $\psi \in D_M$.

**R1.** The generation of the matrix $\mathbf{R}_k$ by (49d) is such that $\mathbf{R}_k$ is symmetric and $\mathbf{R}_k \ge \delta \mathbf{I}$ for some $\delta > 0$.

**G1.** $\lim_{t \to \infty} \gamma_k = \mu > 0$.

**A1.** The limits

$$ \frac{1}{2} \bar{E}[\mathbf{g}_k^\top(\psi) \mathbf{W}_k \mathbf{g}_k(\psi)] = \bar{V}(\psi), \quad (50) $$

$$ \bar{E}[\mathbf{J}_k^\top(\psi) \mathbf{W}_k \mathbf{g}_k(\psi)] = \mathbf{h}(\psi), \quad (51) $$

$$ \bar{E}[\mathbf{J}_k^\top(\psi) \mathbf{W}_k \mathbf{J}_k(\psi)] = \mathbf{G}(\psi). \quad (52) $$

exist for all $\psi \in D_M$, where the limit in expectation is defined as $\bar{E}[f_k(a)] := \lim_{N \to \infty} \frac{1}{N} \sum_{k=1}^N \bar{E}[f_k(a)]$.

**S1.** For each $k, k_0, k \ge k_0$, there exists a random vector $\tilde{\mathbf{z}}_k$ that belongs to the $\sigma$-algebra generated by $\{z_i\}_{i=0}^{k_0}$ but is independent of $\{z_i\}_{i=0}^{k_0}$, such that,

$$ \bar{E}[|\mathbf{z}_k - \tilde{\mathbf{z}}_k|^4] < C \omega^{k-k_0}, \quad C < \infty, \quad \omega < 1. $$

Note that the condition **A1** implies that the data sequence $\{z_k\}$ is asymptotically mean stationary and **S1** that the data generation is exponentially stable. The following theorem can then be stated.

**Theorem 1.** Consider the algorithm (49). Assume that it includes a projection to keep $\hat{\mathbf{y}}_k \in D_M$, and to assure a bounded subsequence of $\{\xi_k\}$. Assume that **M1**, **M2**, **Cr1**, **Cr2**, **A1**, **S1**, **R1**, and **G1** hold. Then $\{\hat{\mathbf{y}}_k\}$ w.p.1 converges either to the set

$$ D_c = \{ \psi \mid \frac{d}{d\psi} \bar{V}(\psi) = 0 \}, $$

or to the boundary of $D_M$ as $k \to \infty$.

**Proof.** See Theorem 4.3 in Ljung and Söderström (1983).

## 4. Simulation experiments

The convergence properties of the proposed algorithm (49) are first illustrated on the example nonlinear EIV regression problem

$$ y = f(\mathbf{x}, \theta) = x_1 \theta_1^{1/2} + x_2 \theta_2^{1/2}. \quad (53) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Trajectories for the recursive NTLS algorithm with multiple initializations around the global minimum for the MISO model $y = x_1\theta_1^{1/2} + x_2\theta_2^{1/2}$. The performance of the algorithm is tested against 25 dB, 20 dB, and 15 dB SNR measurement noise conditions from left to right.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Computation time for a single iteration of batch size N for the offline algorithm vs N individual recursions for the online algorithm, executed on an Intel Core i7-10610U CPU.</figcaption></figure>

In Fig. 1 the performance of the recursive algorithm is tested for six different parametric initializations within the domain of attraction of the function's global minimum, with the gain sequence $\gamma_k = 1/k$, step length $\alpha_k = 0.5$, the LM parameter $\delta_k = 1e^{-4}$, and the information matrix initialized with $R_0 = I_{n_0}$. The input/output data are generated with heterogeneous additive Gaussian white noise (AGWN) sequences. The robustness of the algorithm is tested under increasingly difficult noise conditions with signal to noise ratio (SNR) levels at 25 dB, 20 dB, and 15 dB from left to right, respectively. In each case, the inverse of the algorithm weighting matrices are tuned with the noise covariance matrices

$$W_{y_k}^{-1} = \Sigma_{y_k} = \sigma_y^2, \quad W_{x_k}^{-1} = \Sigma_{x_k} = \begin{bmatrix} \sigma_{x_1}^2 & 0 \\ 0 & \sigma_{x_2}^2 \end{bmatrix}, \qquad (54)$$

where $\sigma_y^2$, $\sigma_{x_1}^2$, and $\sigma_{x_2}^2$ are the variances of the AGWN sequences given to the output and input data, respectively. The results in Fig. 1 demonstrate that the algorithm converges to the global minimum from all initializations, regardless of the noise level.

In order to highlight the improved computational efficiency of the recursive algorithm, the execution time for N individual recursions is compared to that of a single iteration with batch size N for the equivalent offline algorithm detailed in Section 2. The execution time against increasing sample size N is illustrated in Fig. 2. At $N = 1e^4$, a single iteration of the offline algorithm takes approximately 60 s whereas the online algorithm takes <3 s to process the same amount of data sequentially on an Intel Core i7-10610U CPU. Naturally, it would take more than a single iteration for the offline algorithm to converge, further exacerbating this difference.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Tracking ability of the recursive algorithm for time-varying system parameters for the SISO model $y = x\theta^2$. The algorithms are tested with different forgetting factors $\kappa$ with 20 dB SNR noise on the data.</figcaption></figure>

Finally, the tracking ability of the recursive algorithm is illustrated on a simple nonlinear single-input single-output (SISO) model $y = x\theta^2$ (directly estimating $\theta$ rather than $\theta^2$ introduces a nonlinearity in the model). Three time-varying parameter scenarios are considered: (i) a linear ramp between $k = 2e^3$ and $k = 8e^3$, (ii) a triangular pulse between $k = 2e^3$ and $k = 8e^3$, and (iii) a step function at $k = 5e^3$. The input/output data are generated with 20 dB SNR AGWN. The results are illustrated in Fig. 3, where the tracking ability of the algorithm is compared for
three different forgetting factors $κ = \{0.9, 0.99, 0.999\}$. Notably, for sharp changes in the parameter values, superior performance is observed with smaller forgetting factors at around $κ = 0.9$. The results show that the recursive algorithm is well suited for tracking time-varying system parameters, which presents a significant advantage over the equivalent batch algorithm.

## 5. Conclusion

In this work, a recursive identification algorithm for nonlinear EIV models is presented, and is derived from an offline sequential quadratic programming solution. Several algorithmic modifications are proposed to enhance the numerical efficiency and robustness of the algorithm for practical implementation purposes. A proof of convergence is provided using the theoretical tools from Ljung and Söderström (1983) for the analysis of recursive stochastic Gauss–Newton algorithms. The performance of the proposed algorithm is then illustrated in simulation on some toy nonlinear EIV regression problems, including for the tracking of time-varying system parameters. The results demonstrate the improved numerical efficiency and tracking ability of the recursive algorithm compared to its offline counterpart.

### CRediT authorship contribution statement

**Hugo Koide:** Writing – original draft, Software, Project administration, Methodology, Formal analysis, Data curation, Conceptualization. **Jérémy Vayssettes:** Writing – review & editing, Validation, Supervision, Project administration, Funding acquisition, Formal analysis. **Guillaume Mercère:** Writing – review & editing, Validation, Supervision, Project administration, Methodology, Funding acquisition, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Appendix. Proof of Lemma 1

**Proof.** The normal equations associated with (46) are given by

$$ 
\begin{bmatrix} A_k^\top W_{yk} A_k & A_k^\top W_{yk} B_k \\ B_k^\top W_{yk} A_k & B_k^\top W_{yk} B_k + W_{xk} \end{bmatrix} \begin{bmatrix} \hat{\rho}_{\theta_k} \\ \hat{\rho}_{e_{xk}} \end{bmatrix} = - \begin{bmatrix} -A_k^\top W_{yk}(\mathbf{y}_k - \hat{\mathbf{y}}_k) \\ -B_k^\top W_{yk}(\mathbf{y}_k - \hat{\mathbf{y}}_k) + W_{xk} \hat{\mathbf{e}}_{xk} \end{bmatrix}. \quad (A.1) $$

Starting from the second row of (A.1), we have

$$ \hat{\mathbf{p}}_{e_{xk}} = (\mathbf{W}_{xk} + \mathbf{B}_k^\top \mathbf{W}_{yk} \mathbf{B}_k)^{-1} (\mathbf{B}_k^\top \mathbf{W}_{yk} (\mathbf{y}_k - \hat{\mathbf{y}}_k) - \mathbf{W}_{xk} \hat{\mathbf{e}}_{xk} - \mathbf{B}_k^\top \mathbf{W}_{yk} \mathbf{A}_k \hat{\mathbf{p}}_{\theta_k}). \quad (A.2) $$

The first row of (A.1) then yields

$$ (A_k^\top W_{yk} A_k) \hat{\mathbf{p}}_{\theta_k} = A_k^\top W_{yk} (\mathbf{y}_k - \hat{\mathbf{y}}_k) - A_k^\top W_{yk} \mathbf{B}_k \hat{\mathbf{p}}_{e_{xk}}. \quad (A.3) $$

Making use of the matrix inversion identity (Zhang, 2017)

$$ (\mathbf{D} + \mathbf{U}\mathbf{C}^{-1}\mathbf{V})^{-1} = \mathbf{D}^{-1} - \mathbf{D}^{-1}\mathbf{U}(\mathbf{C} + \mathbf{V}\mathbf{D}^{-1}\mathbf{U})^{-1}\mathbf{V}\mathbf{D}^{-1}, \quad (A.4) $$

and substituting (A.2) in (A.3) results in

$$ 
\begin{aligned} (A_k^\top (W_{yk}^{-1} + \mathbf{B}_k W_{xk}^{-1} \mathbf{B}_k^\top)^{-1} A_k) \hat{\mathbf{p}}_{\theta_k} & = A_k^\top (W_{yk}(\mathbf{y}_k - \hat{\mathbf{y}}_k) - W_{yk} \mathbf{B}_k (W_{xk} + \mathbf{B}_k^\top W_{yk} \mathbf{B}_k)^{-1} \\ & \times (\mathbf{B}_k^\top W_{yk}(\mathbf{y}_k - \hat{\mathbf{y}}_k) - W_{xk} \hat{\mathbf{e}}_{xk})). \end{aligned} \quad (A.5) $$

Then, via the matrix push-through identity (Zhang, 2017)

$$ \mathbf{D}^{-1}\mathbf{U}(\mathbf{C}^{-1} + \mathbf{V}\mathbf{D}^{-1}\mathbf{U})^{-1} = (\mathbf{D} + \mathbf{U}\mathbf{C}\mathbf{V})^{-1}\mathbf{U}\mathbf{C}, \quad (A.6) $$

we have that

$$ 
\begin{aligned} (A_k^\top (W_{yk}^{-1} + \mathbf{B}_k W_{xk}^{-1} \mathbf{B}_k^\top)^{-1} A_k) \hat{\mathbf{p}}_{\theta_k} & = A_k^\top (W_{yk}^{-1} + \mathbf{B}_k W_{xk}^{-1} \mathbf{B}_k^\top)^{-1} ((\mathbf{y}_k - \hat{\mathbf{y}}_k) + \mathbf{B}_k \hat{\mathbf{e}}_{xk}). \end{aligned} \quad (A.7) $$

The parameter search direction can then be found via

$$ \hat{\mathbf{p}}_{\theta_k} = (A_k^\top M_k^{-1} A_k)^{-1} A_k^\top M_k^{-1} (\mathbf{y}_k - \hat{\mathbf{y}}_k + \mathbf{B}_k \hat{\mathbf{e}}_{xk}), \quad (A.8) $$

where

$$ \mathbf{M}_k = \mathbf{W}_{yk}^{-1} + \mathbf{B}_k \mathbf{W}_{xk}^{-1} \mathbf{B}_k^\top, \quad (A.9) $$

which is the desired result for (47). Then, the expression (A.2) can be reformulated as

$$ 
\begin{aligned} \hat{\mathbf{p}}_{e_{xk}} = & (\mathbf{W}_{xk} + \mathbf{B}_k^\top \mathbf{W}_{yk} \mathbf{B}_k)^{-1} \mathbf{B}_k^\top \mathbf{W}_{yk} \\ & (\mathbf{y}_k - \hat{\mathbf{y}}_k + \mathbf{B}_k \hat{\mathbf{e}}_{xk} - A_k \hat{\mathbf{p}}_{\theta_k}) - \hat{\mathbf{e}}_{xk}. \end{aligned} \quad (A.10) $$

Let us notice that

$$ 
\begin{aligned} & (\mathbf{W}_{xk} + \mathbf{B}_k^\top \mathbf{W}_{yk} \mathbf{B}_k)^{-1} \\ & = (\mathbf{W}_{xk} + \mathbf{B}_k^\top \mathbf{W}_{yk} \mathbf{B}_k)^{-1} (\mathbf{W}_{xk} + \mathbf{B}_k^\top \mathbf{W}_{yk} \mathbf{B}_k - \mathbf{B}_k^\top \mathbf{W}_{yk} \mathbf{B}_k) \mathbf{W}_{xk}^{-1}, \\ & = \mathbf{W}_{xk}^{-1} - (\mathbf{W}_{xk} + \mathbf{B}_k^\top \mathbf{W}_{yk} \mathbf{B}_k)^{-1} \mathbf{B}_k^\top \mathbf{W}_{yk} \mathbf{B}_k \mathbf{W}_{xk}^{-1}. \end{aligned} \quad (A.11) $$

Applying the push-through identity (A.6) twice to (A.11) and substituting the result in (A.10) then gives

$$ 
\begin{aligned} \hat{\mathbf{p}}_{e_{xk}} = & \mathbf{W}_{xk}^{-1} \mathbf{B}_k^\top (\mathbf{W}_{yk} + \mathbf{W}_{yk} \mathbf{B}_k) (\mathbf{W}_{xk} + \mathbf{B}_k^\top \mathbf{W}_{yk} \mathbf{B}_k)^{-1} \\ & \times \mathbf{B}_k^\top \mathbf{W}_{yk} 1tt(\mathbf{y}_k - \hat{\mathbf{y}}_k + \mathbf{B}_k \hat{\mathbf{e}}_{xk} - A_k \hat{\mathbf{p}}_{\theta_k}) - \hat{\mathbf{e}}_{xk}. \end{aligned} \quad (A.12) $$

Finally, applying the matrix inversion identity (A.4) to (A.12) results in

$$ \hat{\mathbf{p}}_{e_{xk}} = \mathbf{W}_{xk}^{-1} \mathbf{B}_k^\top \mathbf{M}_k^{-1} (\mathbf{y}_k - \hat{\mathbf{y}}_k + \mathbf{B}_k \hat{\mathbf{e}}_{xk} - A_k \hat{\mathbf{p}}_{\theta_k}) - \hat{\mathbf{e}}_{xk}, \quad (A.13) $$

which is the desired result for (48). This concludes the proof.

## Data availability

The data that has been used is confidential.

## References

Adcock, R. J. (1878). A problem in least squares. *The Analyst*, 5(2), 53–54.

Björk, A. (1996). *Numerical methods for least squares problems*. SIAM.

Boggs, P. T., Byrd, R. H., & Schnabel, R. B. (1987). A stable and efficient algorithm for nonlinear orthogonal distance regression. *SIAM Journal on Scientific and Statistical Computing*, 8(6), 1052–1078.

Davila, C. E. (2002). An efficient recursive total least squares algorithm for FIR adaptive filtering. *IEEE Transactions on Signal Processing*, 42(2), 268–280.

Deming, W. E. (1943). Statistical adjustment of data. *Wiley*.

Hansen, P. C., Pereyra, V., & Scherer, G. (2013). Least squares data fitting with applications. JHU Press.

Koide, H., Vayssettes, J., & Mercère, G. (2024). A robust and regularized algorithm for recursive total least squares estimation. *IEEE Control Systems Letters*, 8, 1006–1011.

Ljung, L. (1977). Analysis of recursive stochastic algorithms. *IEEE Transactions on Automatic Control*, 22(4), 551–575.

Ljung, L. (1999). System Identification: Theory for the User, In *Prentice hall information and system sciences series*, (2nd ed.) New Jersey: Prentice-Hall, Inc.

Ljung, L., & Söderström, T. (1983). *Theory and practice of recursive identification*. MIT Press.

Nocedal, J., & Wright, S. J. (2006). *Numerical optimization* (2nd ed.) Springer.

Rhode, S., Usevich, K., Markovsky, I., & Gauterin, F. (2014). A recursive restricted total least-squares algorithm. *IEEE Transactions on Signal Processing*, 62(21), 5652–5662.

Schwetlick, H., & Tiller, V. (1985). Numerical methods for estimating parameters in nonlinear models with errors in the variables. *Technometrics*, 27(1), 17–24.
Schwetlick, H., & Tiller, V. (1989). Nonstandard scaling matrices for trust region Gauss–Newton methods. *SIAM Journal on Scientific and Statistical Computing*, **10**(4), 654–670.

Seber, G., & Wild, C. J. (2003). *Nonlinear regression*. John Wiley & Sons.

Söderström, T. (2018). *Errors-in-variables methods in system identification*. Springer.

Tjoa, I. B., & Biegler, I. (1992). Reduced successive quadratic programming strategy for errors-in-variables estimation. *Computers & Chemical Engineering*, **16**(6), 523–533.

Wolter, K. M., & Fuller, W. A. (1982). Estimation of nonlinear errors-in-variables models. *The Annals of Statistics*, **539–548**.

Zhang, X. D. (2017). *Matrix analysis and applications*. Cambridge University Press.