Recursive nonlinear-system identification using latent variables

I Introduction

In this paper we consider the problem of learning a nonlinear dynamical system model with multiple outputs $y (t)$ and multiple inputs $u (t)$ (when they exist). Generally this identification problem can be tackled using different model structures, with the class of linear models being arguably the most well studied in engineering, statistics and econometrics [27, 15, 6, 2, 7].

Linear models are often used even when the system is known to be nonlinear [10, 25]. However certain nonlinearities, such as saturations, cannot always be neglected. In such cases using block-oriented models is a popular approach to capture static nonlinearities [11]

. Recently, such models have been given semiparametric formulations and identified using machine learning methods, cf.

[22, 23]. To model nonlinear dynamics a common approach is to use Narmax models [26, 5].

In this paper we are interested in recursive identification methods [18]. In cases where the model structure is linear in the parameters, recursive least-squares can be applied. For certain models with nonlinear parameters, the extended recursive least-squares has been used [8]. Another popular approach is the recursive prediction error method which has been developed, e.g., for Wiener models, Hammerstein models, and polynomial state-space models [33, 19, 31].

Nonparametric models are often based on weighted sums of the observed data [24]. The weights vary for each predicted output and the number of weights increases with each observed datapoint. The weights are typically obtained in a batch manner; in [1, 4] they are computed recursively but must be recomputed for each new prediction of the output.

For many nonlinear systems, however, linear models work well as an initial approximation. The strategy in [21] exploits this fact by first finding the best linear approximation using a frequency domain approach. Then, starting from this approximation, a nonlinear polynomial state-space model is fitted by solving a nonconvex problem. This two-step method cannot be readily implemented recursively and it requires input signals with appropriate frequency domain properties.

In this paper, we start from a nominal model structure. This class can be based on insights about the system, e.g. that linear model structures can approximate a system around an operating point. Given a record of past outputs, $y (t)$ and inputs $u (t)$ , that is,

D_{t} ≜ {(y (1), u (1)), \dots, (y (t), u (t))},

a nominal model yields a predicted output $y_{0} (t + 1)$ which differs from the output $y (t + 1)$ . The resulting prediction error is denoted $ε (t + 1)$ [16]. By characterizing the nominal prediction errors in a data-driven manner, we aim to develop a refined predictor model of the system. Thus we integrate classic and data-driven system modeling approaches in a natural way.

The general model class and problem formulation are introduced in Section II. Then in Section III we apply the principle of maximum likelihood to derive a statistically motivated learning criterion. In Section IV this nonconvex criterion is minimized using a majorization-minimization approach that gives rise to a convex user-parameter free method. We derive a computionally efficient recursive algorithm for solving the convex problem, which can be applied to large datasets as well as online learning scenarios. In Section V, we evaluate the proposed method using both synthetic and real data examples.

In a nutshell, the contribution of the paper is a modelling approach and identification method for nonlinear multiple input-multiple output systems that:

explicitly separates modeling based on application-specific insights from general data-driven modelling,
circumvents the choice of regularization parameters and initialization points,
learns parsimonious predictor models,
admits a computationally efficient implementation.

Notation: $E_{i, j}$ denotes the $i j$ th standard basis matrix. $\otimes$ and $⊙$ denote the Kronecker and Hadamard products, respectively. $vec (\cdot)$

is the vectorization operation.

∥ x ∥_{2}

∥ x ∥_{1}

and

∥ X ∥_{W} = \sqrt{tr {X^{⊤} W X}}

, where

W ≻ 0

, denote

ℓ_{2}

ℓ_{1}

- and weighted norms, respectively. The Moore-Penrose pseudoinverse of

X

is denoted

X^{†}

Remark 1.

An implementation of the proposed method is available at https://github.com/magni84/lava.

Ii Problem formulation

Given $D_{t - 1}$ , the $n_{y}$ -dimensional output of a system can always be written as

y (t) = y_{0} (t) + ε (t),

(1)

where $y_{0} (t)$ is any one-step-ahead predictor based on a nominal model. Here we consider nominal models on the form

y_{0} (t) = Θ φ (t),

(2)

where the $p \times 1$ vector $φ (t)$ is a given function of $D_{t - 1}$ and $Θ$ denotes the unknown parameters.

Remark 2.

A typical example of $φ (t)$ is

φ (t) = [y^{⊤} (t - 1) \dots y^{⊤} (t - n_{a}) u^{⊤} (t - 1) \dots u^{⊤} (t - n_{b}) 1]^{⊤},

(3)

in which case the nominal predictor is linear in the data and therefore captures the linear system dynamics. Nonlinearities can be incorporated if such are known about the system, in which case $φ (t)$ will be nonlinear in the data.

The popular Arx model structure, for instance, can be cast into the framework (1) and (2) by assuming that the nominal prediction error $ε (t)$

is a white noise process

[27, 15]. For certain systems, (2) may accurately describe the dynamics of the system around its operation point and consequently the white noise assumption on

ε (t)

may be a reasonable approximation. However, this ceases to be the case even for systems with weak nonlinearities, cf. [10].

Next, we develop a data-driven model of the prediction errors $ε (t)$ in (1), conditioned on the past data $D_{t - 1}$ . Specifically, we assume the conditional model

ε (t) | D_{t - 1} \sim N (Z γ (t), Σ),

(4)

where $Z$ is an $n_{y} \times q$ matrix of unknown latent variables, $Σ$ is an unknown covariance matrix, and the $q \times 1$ vector $γ (t)$ is any given function of $D_{t - 1}$ . This is a fairly general model structure that can capture correlated data-dependent nominal prediction errors.

Note that when $Z \equiv 0$ , the prediction errors are temporally white and the nominal model (2) captures all relevant system dynamics. The latent variable is modeled as random here. Before data collection, we assume $Z$ to have mean $0$

as we have no reason to depart from the nominal model assumptions until after observing data. Using a Gaussian distribution, we thus get

vec (Z) \sim N (0, D),

(5)

where $D$ is an unknown covariance matrix.

Our goal here is to identify a refined predictor of the form

^y (t) = ˆ Θ φ (t)      {^y}_{0} (t) + ˆ Z γ (t)     ^ε (t),

(6)

from a data set $D_{t - 1}$

, by maximizing the likelihood function. The first term is an estimate of the nominal predictor model while the second term tries to capture structure in the data that is not taken into account by the nominal model. Note that when

ˆ Z

is sparse we obtain a parsimonious predictor model.

Remark 3.

The model (1)-(4) implies that we can write the output in the equivalent form

y (t) = Θ φ (t) + Z γ (t) + v (t),

where $v (t)$ is a white process with covariance $Σ$ . In order to formulate a flexible data-driven error model (4), we overparametrize it using a high-dimensional $γ (t)$ . In this case, regularization of $Z$ is desirable, which is achieved by (5). Note that $D$ and $Σ$ are both unknown. Estimating these covariance matrices corresponds to using a data-adaptive regularization, as we show in subsequent sections.

Remark 4.

The nonlinear function $γ (t)$ of $D_{t - 1}$ can be seen as a basis expansion which is chosen to yield a flexible model structure of the errors. In the examples below we will use the Laplace operator basis functions [28]. Other possible choices include the polynomial basis functions, Fourier basis functions, wavelets, etc. [26, 15, 32].

Remark 5.

In (6), $^y (t)$ is a one-step-ahead predictor. However, the framework can be readily applied to $k$ -step-ahead prediction where $φ (t)$ and $γ (t)$ depend on $y (1), \dots, y (t - k)$ .

Iii Latent variable framework

Given a record of $N$ data samples, $D_{N}$ , our goal is to estimate $Θ$ and $Z$ to form the refined predictor (6). In Section III-A, we employ the maximum likelihood approach based on the likelihood function $p (Y | Θ, D, Σ)$ , which requires the estimation of nuisance parameters $D$ and $Σ$ . For notational simplicity, we write the parameters as $Ω = {Θ, D, Σ}$ and in Section III-B we show how an estimator of $Z$ is obtained as a function of $Ω$ and $D_{N}$ .

Iii-a Parameter estimation

We write the output samples in matrix form as

Y = [\begin{matrix} y (1) & \dots & y (N) \end{matrix}] \in R^{n_{y} \times N} .

In order to obtain maximum likelihood estimates of $Ω$ , we first derive the likelihood function by marginalizing out the latent variables from the data distribution:

p (Y | Ω) = \int p (Y | Ω, Z) p (Z) d Z,

(7)

where the data distribution $p (Y | Ω, Z)$ and $p (Z)$ are given by (4) and (5), respectively.

To obtain a closed-form expression of (7), we begin by constructing the regression matrices

	$Φ$	$= [\begin{matrix} φ (1) & \dots & φ (N) \end{matrix}] \in R^{p \times N},$
	$Γ$	$= [\begin{matrix} γ (1) & \dots & γ (N) \end{matrix}] \in R^{q \times N} .$

It is shown in Appendix A that (7) can be written as

p (Y | Ω) = \frac{1}{(2 π)^{N n_{y}} | R |} exp (- \frac{1}{2} ∥ y - F θ ∥_{R^{- 1}}^{2}),

(8)

where

y = vec (Y), θ = vec (Θ), z = vec (Z),

(9)

are vectorized variables, and

	$F$	$= Φ^{⊤} \otimes I_{n_{y}}, G = Γ^{⊤} \otimes I_{n_{y}},$		(10)
	$R$	$≜ G D G^{⊤} + I_{N} \otimes Σ .$		(11)

Note that (8) is not a Gaussian distribution in general, since $R$ may be a function of $Y$ . It follows that maximizing (8) is equivalent to solving

min Ω V (Ω),

(12)

where

V (Ω) ≜ ∥ y - F θ ∥_{R^{- 1}}^{2} + ln | R |

(13)

and $y - F θ = vec (Y - Θ Φ) = vec ([ε (1) \dots ε (N)])$ is nothing but the vector of nominal prediction errors.

Iii-B Latent variable estimation

Next, we turn to the latent variable $Z$ which is used to model the nominal prediction error $ε (t)$ in (4). As we show in Appendix A, the conditional distribution of $Z$ is Gaussian and can be written as

p (Z | Ω, Y) = \frac{1}{\sqrt{(2 π)^{n_{y} q} | Σ_{z} |}} exp (- \frac{1}{2} ∥ z - ζ ∥_{Σ_{z}^{- 1}}^{2}),

(14)

with conditional mean

ζ = D G^{⊤} R^{- 1} (y - F θ),

(15)

and covariance matrix

Σ_{z} = {(D^{- 1} + G^{⊤} (I_{N} \otimes Σ^{- 1}) G)}_{N}^{- 1} .

(16)

An estimate $ˆ Z$ is then given by evaluating the conditional (vectorized) mean (15) at the optimal estimate $ˆ Ω$ obtained via (12).

Iv Majorization-minimization approach

The quantities in the refined predictor (6) are readily obtained from the solution of (12). In general, $V (Ω)$ may have local minima and (12) must be tackled using computationally efficient iterative methods to find the optimum. The obtained estimates will then depend on the choice of initial point $Ω_{0}$ . Such methods includes the majorization-minimization approach [12, 35]

, which in turn include Expectation Maximization

[9] as a special case.

The majorization-minimization approach is based on finding a majorizing function $V^{'} (Ω | ˜ Ω)$ with the following properties:

V (Ω) \leq V^{'} (Ω | ˜ Ω), \forall Ω

(17)

with equality when $Ω = ˜ Ω$ . The key is to find a majorizing function that is simpler to minimize than $V (Ω)$ . Let $Ω_{k + 1}$ denote the minimizer of $V^{'} (Ω | Ω_{k})$ . Then

\begin{matrix} V (Ω_{k + 1}) & \leq V^{'} (Ω_{k + 1} | Ω_{k}) \leq V^{'} (Ω_{k} | Ω_{k}) = V (Ω_{k}) . \end{matrix}

(18)

This property leads directly to an iterative scheme that decreases $V (Ω)$ monotonically, starting from an initial estimate $Ω_{0}$ .

Given the overparameterized error model (4), it is natural to initialize at points in the parameter space which correspond to the nominal predictor model structure (2). That is,

Ω_{0} = {Θ_{0}, 0, Σ_{0}}, where Σ_{0} ≻ 0,

(19)

at which $ˆ Z \equiv 0$ .

Iv-a Convex majorization

For a parsimonious parameterization and computationally advantageous formulation, we consider a diagonal structure of the covariance matrices in (4), i.e., we let

Σ = diag (σ_{1}, σ_{2}, \dots, σ_{n_{y}}),

(20)

and we let $D_{i} = diag (d_{i, 1}, \dots, d_{i, q})$ denote the covariance matrix corresponding to the $i$ th row of $Z$ , so that

D = n_{y} \sum i = 1 D_{i} \otimes E_{i, i} .

(21)

We begin by majorizing (13) via linearization of the concave term $ln | R |$ :

\begin{matrix} ln | R | & \leq ln | ˜ R | - tr {{˜ R}^{- 1} (I_{N} \otimes ˜ Σ)} - tr {G^{⊤} {˜ R}^{- 1} G ˜ D} + tr {{˜ R}^{- 1} (I_{N} \otimes Σ)} + tr {G^{⊤} {˜ R}^{- 1} G D}, \end{matrix}

(22)

where $˜ D$ and $˜ Σ$ are arbitrary diagonal covariance matrices and $˜ R$ is obtained by inserting $˜ D$ and $˜ Σ$ into (11). The right-hand side of the inequality above is a majorizing tangent plane to $ln | R |$ , see Appendix B. The use of (22) yields a convex majorizing function of $V (Ω)$ in (12):

\begin{matrix} V^{'} (Ω | ˜ Ω) & = ∥ y - F θ ∥_{R^{- 1}}^{2} + tr {{˜ R}^{- 1} (I_{N} \otimes Σ)} + tr {G^{⊤} {˜ R}^{- 1} G D} + ˜ K, \end{matrix}

(23)

where $˜ K = ln | ˜ R | - tr {{˜ R}^{- 1} (I_{N} \otimes ˜ Σ)} - tr {G^{⊤} {˜ R}^{- 1} G ˜ D}$ is a constant. To derive a computationally efficient algorithm for minimizing (23), the following theorem is useful:

Theorem 1.

The majorizing function (23) can also be written as

V^{'} (Ω | ˜ Ω) = min Z V^{'} (Ω | Z, ˜ Ω)

(24)

where

\begin{matrix} V^{'} (Ω | Z, ˜ Ω) & = ∥ Y - Θ Φ - Z Γ ∥_{Σ^{- 1}}^{2} + ∥ vec (Z) ∥_{D^{- 1}}^{2} + tr {{˜ R}^{- 1} (I_{N} \otimes Σ)} + tr {G^{⊤} {˜ R}^{- 1} G D} + ˜ K . \end{matrix}

(25)

The minimizing $Z$ is given by the conditional mean (15).

Proof.

The problem in (24) has a minimizing $Z$ which, after vectorization, equals $ζ$ in (15). Inserting the minimizing $Z$ into (25) and using (41) yields (23). ∎

Remark 6.

The augmented form in (23), enables us to solve for the nuisance parameters $D$ and $Σ$ in closed-form and also yields the conditional mean $ˆ Z$ as a by-product.

To prepare for the minimization of the function (25) we write the matrix quantities using variables that denote the $i$ th rows of the following matrices:

Y = ⎡ ⎢ ⎢ ⎢ ⎣ \begin{matrix} ⋮ y_{i}^{⊤} ⋮ \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎦, Θ = ⎡ ⎢ ⎢ ⎢ ⎣ \begin{matrix} ⋮ θ_{i}^{⊤} ⋮ \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎦, Z = ⎡ ⎢ ⎢ ⎢ ⎣ \begin{matrix} ⋮ z_{i}^{⊤} ⋮ \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎦, Γ = ⎡ ⎢ ⎢ ⎢ ⎣ \begin{matrix} ⋮ γ_{i}^{⊤} ⋮ \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎦ .

Theorem 2.

After concentrating out the nuisance parameters, the minimizing arguments $Θ$ and $Z$ of the function (25) are obtained by solving the following convex problem:

min Θ, Z n_{y} \sum i = 1 (∥ y_{i} - Φ^{⊤} θ_{i} - Γ^{⊤} z_{i} ∥_{2} + ∥ w_{i} ⊙ z_{i} ∥_{1})

(26)

where

	$w_{i}$	$= {[\begin{matrix} w_{i, 1} & \dots & w_{i, q} \end{matrix}]}^{⊤}, w_{i, j} = \sqrt{\frac{γ_{j}^{⊤} {˜ R}_{i}^{- 1} γ_{j}}{tr {{˜ R}_{i}^{- 1}}}}$		(27)
	${˜ R}_{i}$	$= Γ {˜ D}_{i} Γ^{⊤} + {˜ σ}_{i} I_{N} .$		(28)

The closed-form expression for the minimizing nuisance parameters (20) and (21) are given by

{^σ}_{i}

=∥yi−Φ⊤θi−Γ⊤zi∥2√tr{˜R−1i},^di,j=|zi,j|√γ⊤j˜R−1iγj.

(29)

Proof.

See Appendix C. ∎

Remark 7.

Problem (26) contains a data-adaptive regularizing term which typically leads to parsimonious estimates of $Z$ , cf. [30].

Remark 8.

Majorizing at a nominal predictor model, i.e. $˜ Ω = {˜ Θ, 0, ˜ Σ}$ as discussed above, yields ${˜ R}_{i} = {˜ σ}_{i} I_{N}$ and the weights are given by

w_{i, j} = \frac{∥ γ_{j} ∥_{2}}{\sqrt{N}} .

(30)

Then problem (26) and consequently the minimization of (24) becomes invariant with respect to $˜ Σ$ .

The iterative scheme (18) is executed by initializing at $˜ Ω := Ω_{0}$ and solving (26). The procedure is then repeated by updating the majorization point using the new estimate $˜ Ω := ˆ Ω$ . It follows that the estimates will converge to a local minima of (13). The following theorem establishes that the local minima found, and thus the resulting predictor (6), is invariant to $Ω_{0}$ in the form (19).

Theorem 3.

All initial points $Ω_{0}$ in the form (19) result in the same sequence of minimizers $({ˆ Θ}_{k}, {ˆ Z}_{k})$ of (26), for all $k > 0$ . Moreover, the sequence $({ˆ D}_{k}, {ˆ Σ}_{k})$ converges to a unique point when $k \to \infty$ .

Proof.

See Appendix D. ∎

Remark 9.

As a result we may initialize the scheme (18) at $Ω_{0} = {0, 0, I_{n_{y}}}$ . This obviates the need for carefully selecting an initialization point, which would be needed in e.g. the Expectation Maximization algorithm.

Iv-B Recursive computation

We now show that the convex problem (26) can be solved recursively, for each new sample $y (t)$ and $u (t)$ .

Iv-B1 Computing $ˆ Θ$

If we fix $Z$ and only solve for $Θ$ , the solution is given by

(31)

where

¯ ¯¯¯ ¯ Θ = Y Φ^{†} % a n d H^{⊤} = Γ Φ^{†} .

Note that both $¯ ¯¯¯ ¯ Θ$ and $H$ are independent of $Z$ , and that they can be computed for each sample $t$ using a standard recursive least-squares (LS) algorithm:

$¯ ¯¯¯ ¯ Θ (t)$	$= ¯ ¯¯¯ ¯ Θ (t - 1) + (y (t) - ¯ ¯¯¯ ¯ Θ (t - 1) φ (t)) φ^{⊤} (t) P (t)$	(32)
$H (t)$	$= H (t - 1) + P (t) φ (t) (γ^{⊤} (t) - φ^{⊤} (t) H (t - 1))$	(33)
$P (t)$	$= P (t - 1) - \frac{P (t - 1) φ (t) φ^{⊤} (t) P (t - 1)}{1 + φ^{⊤} (t) P (t - 1) φ (t)} .$	(34)

Remark 10.

Natural initial values are $¯ ¯¯¯ ¯ Θ (0) = 0$ and $H (0) = 0$ . The matrix $Φ^{†}$ equals $Φ^{⊤} (Φ Φ^{⊤})^{- 1}$ when $t \geq p$ samples yield a full-rank matrix $Φ$ . The matrix $P (t)$ converges to $(Φ Φ^{⊤})^{- 1}$ . A common choice for the initial value of $P (t)$ is $P (0) = c I$ , where a larger constant $c > 0$ leads to a faster convergence of (34), cf. [27, 29].

Iv-B2 Computing $ˆ Z$

Using (31), we concentrate out $Θ$ from (26) to obtain

V^{'} (Z) = n_{y} \sum i = 1 (∥ ξ_{i} - (Γ^{⊤} - Φ^{⊤} H) z_{i} ∥_{2} + ∥ w_{i} ⊙ z_{i} ∥_{1})

where

ξ_{i} = y_{i} - Φ^{⊤} {¯ ¯ ¯ θ}_{i} .

In Appendix E it is shown how the minimum of $V^{'} (Z)$ can be found via cyclic minimization with respect to the elements of $Z$ , similar to what has been done in [36] in a simpler case. This iterative procedure is implemented using recursively computed quantities and produces an estimate $ˆ Z (t)$ at sample $t$ .

Iv-B3 Summary of the algorithm

The algorithm computes $ˆ Θ (t)$ and $ˆ Z (t)$ recursively by means of the following steps at each sample $t$ :

Compute $¯ ¯¯¯ ¯ Θ (t)$ and $H (t)$ , using (32)-(34).
Compute $ˆ Z (t)$ via the cyclic minimization of $V^{'} (Z)$ . Cycle through all elements $L$ times.
Compute $ˆ Θ (t) = ¯ ¯¯¯ ¯ Θ (t) - ˆ Z (t) H^{⊤} (t)$

The estimates are initialized as $ˆ Θ (0) = 0$ and $ˆ Z (0) = 0$ . In practice, small $L$ works well since we cycle $L$ times through all elements of $Z$ for each new data sample. The computational details are given in Algorithm 1 in Appendix E, which can be readily implemented e.g. in Matlab.

V Numerical experiments

In this section we evaluate the proposed method and compare it with two alternative identification methods.

V-a Identification methods and experimental setup

The numerical experiments were conducted as follows. Three methods have been used: LS identification of affine Arx, Narx using wavelet networks (Wave for short), and the latent variable method (Lava) presented in this paper. From our numerical experiments we found that performing even only one iteration of the majorization-minimization algorithm produces good results, and doing so leads to a computationally efficient recursive implementation (which we denote Lava-R for recursive).

For each method the function $φ (t)$ is taken as the linear regressor in (3). Then the dimension of $φ (t)$ equals $p = n_{y} n_{a} + n_{u} n_{b} + 1$ . For affine ARX, the model is given by

{^y}_{A R X} (t) = ¯ ¯¯¯ ¯ Θ φ (t),

where $¯ ¯¯¯ ¯ Θ$ is estimated using recursive least squares [27]. Note that in Lava-R, $¯ ¯¯¯ ¯ Θ$ is computed as a byproduct (32).

For the wavelet network, nlarx in the System Identification Toolbox for Matlab was used, with the number of nonlinear units automatically detected [17].

For Lava-R, the model is given by (6) and $ˆ Θ, ˆ Z$ are found by the minimization of (26) using ${˜ Ω}_{0} = {0, 0, I_{n_{y}}}$ . The minimization is performed using the recursive algorithm in Section IV-B3 with $L = 5$ . The nonlinear function $γ (t)$ of the data $D_{t - 1}$ can be chosen to be a set of basis functions evaluated at $φ (t)$ . Then $Z γ (t)$ can be seen as a truncated basis expansion of some nonlinear function. In the numerical examples, $γ (t)$ uses the Laplace basis expansion due to its good approximation properties [28]. Each element in the expansion is given by

γ_{k_{1}, \dots, k_{p - 1}} (t) = p - 1 \prod i = 1 \frac{1}{\sqrt{ℓ_{i}}} sin (\frac{π k_{i} (φ_{i} (t) + ℓ_{i})}{2 ℓ_{i}}),

(35)

where $ℓ_{i}$ are the boundaries of the inputs and outputs for each channel and $k_{i} = 1, \dots, M$ are the indices for each element of $γ (t)$ . Then the dimension of

γ (t) = [γ_{1, \dots, 1} (t) \dots γ_{p - 1, \dots, p - 1} (t)]^{⊤}

equals $q = M^{p - 1}$ , where $M$ is a user parameter which determines the resolution of the basis.

Finally, an important part of the identification setup is the choice of input signal. For a nonlinear system it is important to excite the system both in frequency and in amplitude. For linear models a commonly used input signal is a pseudorandom binary sequence (PRBS), which is a signal that shifts between two levels in a certain fashion. One reason for using PRBS is that it has good correlation properties [27]

. Hence, PRBS excites the system well in frequency, but poorly in amplitude. A remedy to the poor amplitude excitation is to multiply each interval of constant signal level with a random factor that is uniformly distributed on some interval

[0, A]

, cf. [34]. Hence, if the PRBS takes the values

- 1

and

1

, then the resulting sequence will contain constant intervals with random amplitudes between

- A

and

A

. We denote such a random amplitude sequence

RS (A)

where

A

is the maximum amplitude.

V-B Performance metric

For the examples considered here the system does not necessarily belong to the model class, and thus there is no true parameter vector to compare with the estimated parameters. Hence, the different methods will instead be evaluated with respect to the simulated model output ${^y}_{s} (t)$ . For Lava-R

	$ˆ φ (t) = [{^y}_{s}^{⊤} (t - 1) \dots {^y}_{s}^{⊤} (t - n_{a}) u^{⊤} (t - 1) \dots u^{⊤} (t - n_{b}) 1]^{⊤} .$
	${ˆ y}_{s} (t) = ˆ Θ ˆ φ (t) + ˆ Z ˆ γ (t)$

and $ˆ γ (t)$ is computed as (35) with $φ (t)$ replaced by $^φ (t)$ .

The performance can then be evaluated using the root mean squared error (RMSE) for each output channel $i$ ,

{RMSE}_{i} = \sqrt{\frac{1}{T} T \sum t = 1 E [∥ y_{i} (t) - {^y}_{s, i} (t) ∥_{2}^{2}]} .

The expectations are computed using 100 Monte Carlo simulations on validation data.

For the dataset collected from a real system, it is not possible to evaluate the expectation in the RMSE formula. For such sets we use the fit of the data, i.e.,

{FIT}_{i} = 100 (1 - \frac{∥ y_{i} - {^y}_{s, i} ∥_{2}}{} ∥ y_{i} - {¯ y}_{i} 1 ∥_{2}),

where ${^y}_{s, i}$ contains the simulated outputs for channel $i$ , ${¯ y}_{i}$ is the empirical mean of $y_{i}$ and $1$ is a vector of ones. Hence, FIT compares the simulated output errors with those obtained using the empirical mean as the model output.

V-C System with saturation

Consider the following state-space model,

$x_{1} (t + 1)$	$= {sat}_{2} [0.9 x_{1} (t) + 0.1 u_{1} (t)]$	(36)
$x_{2} (t + 1)$	$= 0.08 x_{1} (t) + 0.9 x_{2} (t) + 0.6 u_{2} (t)$	(37)
$y (t)$	$= x (t) + e (t) .$	(38)

where $x (t) = {[\begin{matrix} x_{1} (t) & x_{2} (t) \end{matrix}]}^{⊤}$ and

{sat}_{a} (x) = {\begin{matrix} x & if | x | < a sign (x) a & if | x | \geq a \end{matrix} .

A block-diagram for the above system is shown in Fig. 1. The measurement noise $e (t)$ was chosen as a white Gaussian process with covariance matrix $σ I$ where $σ = 2.5 \cdot 10^{- 3}$ .

Data was collected from the system using an $RS (A)$ input signal for several different amplitudes $A$ . The identification was performed using $n_{a} = 1$ , $n_{b} = 1$ , $M = 4$ , and $N = 1000$ data samples. This means that $p = 5$ and $q = 256$ , and therefore there are 10 parameters in $Θ$ and $512$ in $Z$ .

Fig. 1: A block diagram of the system used in Example V-C.

Note that, for sufficiently low amplitudes $A$ , $x_{1} (t)$ will be smaller than the saturation level $a = 2$ for all $t$ , and thus the system will behave as a linear system. However, when $A$ increases, the saturation will affect the system output more and more. The RMSE was computed for eight different amplitudes $A$ , and the result is shown in Fig. 2. It can be seen that for small amplitudes, when the system is effectively linear, the Arx model gives a marginally better result than Lava-R and Wave. However, as the amplitude is increased, the nonlinear effects become more important, and Lava-R outperforms both Wave and Arx models.

Fig. 2: The RMSE for Example (V-C) computed for different input amplitudes, using Lava-R (solid), affine Arx (dashed) and Wave (dash-dotted).

V-D Water tank

In this example a real cascade tank process is studied. It consists of two tanks mounted on top of each other, with free outlets. The top tank is fed with water by a pump. The input signal is the voltage applied to the pump, and the output consists of the water level in the two tanks. The setup is described in more detail in [34]. The data set consists of 2500 samples collected every five seconds. The first 1250 samples where used for identification, and the last 1250 samples for validation.

Fig. 3: The output in Example V-D (blue), plotted together with the output of the model identified by Lava-R (red). The system was identified using the first 1250 data samples. The validation set consisted of the remaining 1250 samples.

The identification was performed using $n_{a} = 2$ , $n_{b} = 2$ and $M = 3$ . With two outputs, this results in 14 parameters in $Θ$ and 1458 parameters in $Z$ . Lava-R found a model with only 37 nonzero parameters in $Z$ , and the simulated output together with the measured output are shown in Fig. 3. The FIT values, computed on the validation data are shown in Table I. It can be seen that an affine ARX model gives a good fit, but also that using Lava-R the FIT measure can be improved significantly. In this example, Wave did not perform very well.

	Lava-R	Wave	Arx
Upper tank	$91.6 %$	$79.2 %$	$84.9 %$
Lower tank	$90.8 %$	$76.9 %$	$78.6 %$

TABLE I: FIT for Example V-D.

V-E Pick-and-place machine

In the final example, a real pick-and-place machine is studied. This machine is used to place electronic components on a circuit board, and is described in detail in [13]. This system exhibits saturation, different modes, and other nonlinearities. The data used here are from a real physical process, and were also used in e.g. [3, 14, 20]. The data set consists of a 15s recording of the single input $u (t)$ and the vertical position of the mounting head $y (t)$ . The data was sampled at 50 Hz, and the first 8s ( $N = 400$ ) were used for identification and the last 7s for validation.

The identification was performed using $n_{a} = 2$ , $n_{b} = 2$ and $M = 6$ . For the SISO system considered here, this results in 5 parameters in $Θ$ and 1296 parameters in $Z$ . Lava-R found a model with 33 of the parameters in $Z$ being nonzero, the output of which is shown in Fig. 4.

Fig. 4: The output in Example V-E (blue), plotted together with the output of the model identified by Lava-R (red). The system was identified using the first part of the data, while the validation set consisted of the remaining samples indicated after the dashed line.

	Lava-R	Wave	Arx
FIT	$83.2 %$	$78.2 %$	$73.1 %$

TABLE II: FIT for Example V-E.

The FIT values, computed on the validation data, for Lava-R, Wave and affine Arx are shown in Table II. Lava-R outperforms Narx using wavelet networks, and both are better than Arx.

Vi Conclusion

We have developed a method for learning nonlinear systems with multiple outputs and inputs. We began by modelling the errors of a nominal predictor using a latent variable formulation. The nominal predictor could for instance be a linear approximation of the system but could also include known nonlinearities. A learning criterion was derived based on the principle of maximum likelihood, which obviates the tuning of regularization parameters. The criterion is then minimized using a majorization-minimization approach. Specifically, we derived a convex user-parameter free formulation, which led to a computationally efficient recursive algorithm that can be applied to large datasets as well as online learning problems.

The method introduced in this paper learns parsimonious predictor models and captures nonlinear system dynamics. This was illustrated via synthetic as well as real data examples. As shown in these examples a recursive implementation of the proposed method was capable of outperforming a batch method using a Narx model with a wavelet network.

Appendix A Derivation of the distributions (8) and (14)

We start by computing $p (Y | Ω)$ given in (7). The function $p (Y | Ω, Z)$ can be found from (1)–(4

) and the chain rule:

\begin{matrix} p (Y | Ω, Z) & = N \prod t = 1 p_{ε} (y (t) - Θ φ (t) | D_{t - 1}, Z), \end{matrix}

(39)

where we have neglected initial conditions [27]. Since

\begin{matrix} p_{ε} (y (t) - Θ φ (t) | D_{t - 1}, Z) \propto exp (- \frac{1}{2} ∥ y (t) - Θ φ (t) - Z γ (t) ∥_{Σ^{- 1}}^{2}), \end{matrix}

it follows that

\begin{matrix} p (Y | Ω, Z) = \frac{1}{\sqrt{(2 π)^{n_{y} N} | Σ |^{N}}} exp (- \frac{1}{2} ∥ Y - Θ Φ - Z Γ ∥_{Σ^{- 1}}^{2}) . \end{matrix}

(40)

Using the vectorized variable in (9)-(10) we can see that

vec (Θ Φ) = F θ and vec (Z Γ) = G z .

and thus,

∥ Y - Θ Φ - Z Γ ∥_{Σ^{- 1}}^{2} = ∥ y - F θ - G z ∥_{I_{N} \otimes Σ^{- 1}}^{2} .

Next, we note that the following useful equality holds:

∥ y - F θ - G z ∥_{I_{N} \otimes Σ^{- 1}}^{2} + ∥ z ∥_{D^{- 1}}^{2} = ∥ y - F θ ∥_{R^{- 1}}^{2} + ∥ z - ζ ∥_{Σ_{z}^{- 1}}^{2}

(41)

where $R$ is given by (11), $ζ$ by (15), and $Σ_{z}$ by (16). To see that the equality holds, expand the norms on both sides of (41) and apply the matrix inversion lemma.

The sought-after distribution $p (Y | Ω)$ is given by (7). By using (41) it follows that

p (Y | Ω, Z) p (Z) \propto exp (- \frac{1}{2} ∥ y - F θ ∥_{R^{- 1}}^{2}) exp (- \frac{1}{2} ∥ z - ζ ∥_{Σ_{z}^{- 1}}^{2}) .

(42)

with the normalization constant $((2 π)^{n_{y} (N + q)} | Σ |^{N} | D |)^{- 1 / 2}$ . Noting that

\int exp (- \frac{1}{2} ∥ z - ζ ∥_{Σ_{z}^{- 1}}^{2}) d Z = \sqrt{(2 π)^{n_{y} q} | Σ_{z} |}

it can be seen that

p (Y | Ω) = \frac{1}{\sqrt{(2 π)^{N n_{y}} | R |}} exp (- \frac{1}{2} ∥ y - F θ ∥_{R^{- 1}}^{2}),

(43)

which proves (8). To obtain an expression for $p (Z | Ω, Y)$ simply insert (42) and (43) into Bayes’ rule to get (14).

Appendix B Derivation of the majorizing tangent plane (22)

The first-order Taylor expansion of the log-determinant can be written as

\begin{matrix} ln | R | & ≃ ln | ˜ R | + (\partial_{σ} ln | R |) |_{R = ˜ R} (σ - ˜ σ) + (\partial_{d} ln | R |_{R = ˜ R}) (d - ~ d) \end{matrix}

where $σ$ is the vector of diagonal elements in $Σ$ and $d$ contains the diagonal elements in $D$ .

For the derivatives with respect to $d$ we have

\frac{\partial}{\partial d_{i, j}} ln | R | = tr (R^{- 1} \frac{\partial R}{\partial d_{i, j}}) = tr (G^{⊤} R^{- 1} G \frac{\partial D}{\partial d_{i, j}}) .

Note that

n_{y} \sum i = 1 q \sum j = 1 d_{i, j} \frac{\partial D}{\partial d_{i, j}} = D

so Hence

\partial_{d} ln | R |_{R = ˜ R} (d - ~ d) = tr (G^{⊤} {˜ R}^{- 1} G (D - ˜ D)) .

In the same way

\partial_{σ} ln | R |_{R = ˜ R} (σ - ~ σ) = tr ({˜ R}^{- 1} (I_{N} \otimes (Σ - ˜ Σ)) .

Since $ln | R |$ is concave in $σ$ and $d$ , it follows that

ln | R | \leq ˜ K + tr (G^{⊤} {˜ R}^{- 1} G D) + tr ({˜ R}^{- 1} (I_{N} \otimes Σ))

(44)

where

˜ K = ln | ˜ R | - tr (G^{⊤} {˜ R}^{- 1} G ˜ D) - tr ({˜ R}^{- 1} (I_{N} \otimes ˜ Σ)) .

Appendix C Proof of Theorem 2

It follows from (21) that

R = n_{y} \sum i = 1 (R_{i} \otimes E_{i, i})

where $R_{i} = Γ^{⊤} D_{i} Γ + σ_{i} I_{N}$ . Hence,

R^{- 1} = n_{y} \sum i = 1 (R_{i}^{- 1} \otimes E_{i, i}) .

Thus, we can rewrite (25) as (to within an additive constant):

V^{'} (Ω | Z, ˜ Ω) = n_{y} \sum i = 1 (\frac{1}{σ_{i}} ∥ {¯ y}_{i} ∥_{2}^{2} + ∥ z_{i} ∥_{D_{i}^{- 1}}^{2} + σ_{i} tr ({˜ R}_{i}^{- 1}) + tr (Γ {˜ R}_{i}^{- 1} Γ^{⊤} D_{i})) .

(45)

where ${¯ y}_{i} = y_{i} - Φ^{⊤} θ_{i} - Γ^{⊤} z_{i}$ .

We next derive analytical expressions for the $Σ$ and $D$ that minimize $V^{'} (Ω | Z, ˜ Ω)$ . Note that

\frac{\partial}{\partial σ_{i}} V^{'} (Ω | Z, ˜ Ω) = - \frac{1}{σ_{i}^{2}} ∥ {¯ y}_{i} ∥_{2}^{2} + tr ({˜ R}_{i}^{- 1}),

and setting the derivative to zero yields the estimate (29). In the same way it can be seen that the minimum of $d_{i, j}$ is attained at (29).

Inserting ${^σ}_{i}$ and ${^d}_{i, j}$ into (45), we see that we can find the minimizing $Θ$ and $Z$ by minimizing

	$2$	$n_{y} \sum i = 1 (\sqrt{tr ({˜ R}_{i}^{- 1})} ∥ y_{i} - Φ^{⊤} θ_{i} - Γ^{⊤} z_{i} ∥_{2} +$
		$q \sum j = 1 \| z_{i, j} \| \sqrt{{¯ ¯ ¯ γ}_{j}^{⊤} {˜ R}_{i}^{- 1} {¯ ¯ ¯ γ}_{j}}) .$

Since term $i$ in the above sum is invariant with respect to $θ_{k}$ and $z_{k}$ for $k \neq i$ , we can divide term $i$ by $2 \sqrt{tr ({˜ R}_{i}^{- 1})}$ , and see that minimizing the criterion above is equivalent to (26).

Appendix D Proof of Theorem 3

Initializing (18) at ${¯ ¯¯¯ ¯ Ω}_{0} = {0, 0, I_{n_{y}}}$ and $Ω_{0} = {Θ_{0}, 0, Σ_{0}}$ where $Σ_{0} = diag (σ_{1}^{(0)}, \dots, σ_{n_{y}}^{(0)})$ , produces two sequences denoted ${¯ ¯¯¯ ¯ Ω}_{k} = {{¯ ¯¯¯ ¯ Θ}_{k}, {¯ ¯¯¯ ¯ D}_{k}, {¯ ¯¯¯ ¯ Σ}_{k}}$ and $Ω_{k} = {Θ_{k}, D_{k}, Σ_{k}}$ for $k > 0$ , respectively. This results also in sequences ${¯ ¯¯ ¯ Z}_{k}$ and $Z_{k}$ . The theorem states that:

Θ_{k} = {¯ ¯¯¯ ¯ Θ}_{k} and Z_{k} = {¯ ¯¯ ¯ Z}_{k}

(46)

D_{k} - {¯ ¯¯¯ ¯ D}_{k} \to 0 and Σ_{k} - {¯ ¯¯¯ ¯ Σ}_{k} \to 0 .

(47)

We now show the stronger result that the covariance matrices converge as

D_{i}^{(k)} = c_{i}^{(k)} {¯ ¯¯¯ ¯ D}_{i}^{(k)}, σ_{i}^{(k)} = c_{i}^{(k)} {¯ ¯ ¯ σ}_{i}^{(k)}, \forall k > 0,

(48)

where $c_{i}^{(k)} = (σ_{i}^{(0)})^{\frac{1}{2^{k}}}$ . Note that $c_{i}^{(k)} \to 1$ as $k \to \infty$ . Hence (48) implies (47). We prove (48) and (46) by induction. That (48) and (46) holds for $k = 1$ follows directly from Theorem 2. Now assume that (48) holds for some $k \geq 1$ . Let

	${¯ ¯¯¯ ¯ R}_{i}$	$= Γ {¯ ¯¯¯ ¯ D}_{i}^{(k)} Γ^{⊤} + {¯ ¯ ¯ σ}_{i}^{(k)} I_{N},$

Recursive nonlinear-system identification using latent variables

Authors

Variational Autoencoders and Nonlinear ICA: A Unifying Framework

Invariant Gaussian Process Latent Variable Models and Application in Causal Discovery

Reduced-dimensional Monte Carlo Maximum Likelihood for Latent Gaussian Random Field Models

Identification of complex systems in the basis of wavelets

Restricted maximum-likelihood method for learning latent variance components in gene expression data with known and unknown confounders

Learning latent variable structured prediction models with Gaussian perturbations

Robust Nonlinear Component Estimation with Tikhonov Regularization

This week in AI

I Introduction

Remark 1.

Ii Problem formulation

Remark 2.

Remark 3.

Remark 4.

Remark 5.

Iii Latent variable framework

Iii-a Parameter estimation

Iii-B Latent variable estimation

Iv Majorization-minimization approach

Iv-a Convex majorization

Theorem 1.

Proof.

Remark 6.

Theorem 2.

Proof.

Remark 7.

Remark 8.

Theorem 3.

Proof.

Remark 9.

Iv-B Recursive computation

Iv-B1 Computing ˆΘ

Remark 10.

Iv-B2 Computing ˆZ

Iv-B3 Summary of the algorithm

V Numerical experiments

V-a Identification methods and experimental setup

V-B Performance metric

V-C System with saturation

V-D Water tank

V-E Pick-and-place machine

Vi Conclusion

Appendix A Derivation of the distributions (8) and (14)

Appendix B Derivation of the majorizing tangent plane (22)

Appendix C Proof of Theorem 2

Appendix D Proof of Theorem 3

Iv-B1 Computing $ˆ Θ$

Iv-B2 Computing $ˆ Z$