An Incremental Approach to Online Dynamic Mode Decomposition for Time-Varying Systems with Applications to EEG Data Modeling

TL;DR

This paper introduces incremental algorithms for online dynamic mode decomposition tailored to time-varying systems, with applications to EEG data modeling, enabling real-time analysis of evolving dynamics.

Contribution

It develops novel incremental DMD algorithms for systems with time-varying dynamics, extending existing methods to handle high-dimensional, non-stationary data efficiently.

Findings

Effective reconstruction and prediction of EEG data.

Algorithms successfully handle singular data matrices.

Demonstrated applicability to both autonomous and non-autonomous systems.

Abstract

Dynamic Mode Decomposition (DMD) is a data-driven technique to identify a low dimensional linear time invariant dynamics underlying high-dimensional data. For systems in which such underlying low-dimensional dynamics is time-varying, a time-invariant approximation of such dynamics computed through standard DMD techniques may not be appropriate. We focus on DMD techniques for such time-varying systems and develop incremental algorithms for systems without and with exogenous control inputs. We build upon the work in [35] to scenarios in which high dimensional data are governed by low dimensional time-varying dynamics. We consider two classes of algorithms that rely on (i) a discount factor on previous observations, and (ii) a sliding window of observations. Our algorithms leverage existing techniques for incremental singular value decomposition and allow us to determine an appropriately reduced model at each time and are applicable even if data matrix is singular. We apply the developed algorithms for autonomous systems to Electroencephalographic (EEG) data and demonstrate their effectiveness in terms of reconstruction and prediction. Our algorithms for non-autonomous systems are illustrated using randomly generated linear time-varying systems.

Equations

$${\mathbf{x}}_{k+1}$$

$$\mathbf{X}=\left[\begin{array}{cccc}{\mathbf{x}}_1& {\mathbf{x}}_2& \cdots & {\mathbf{x}}_m\end{array}\right],\mathbf{Y}=\left[\begin{array}{cccc}{\mathbf{x}}_2& {\mathbf{x}}_3& \cdots & {\mathbf{x}}_{m+1}\end{array}\right].$$

$${\mathbf{x}}_{k+1}=\bar{\mathbf{A}}{\mathbf{x}}_k,$$

$$\bar{\mathbf{A}}=\mathbf{Y}{\bar{\mathbf{V}}}_x{\bar{\mathbf{\Sigma }}}_x^{-1}{\bar{\mathbf{U}}}_x^{\ast }.$$

$$\tilde{\mathbf{A}}={\bar{\mathbf{U}}}_x^{\ast }\bar{\mathbf{A}}{\bar{\mathbf{U}}}_x={\bar{\mathbf{U}}}_x^{\ast }\mathbf{Y}{\bar{\mathbf{V}}}_x{\bar{\mathbf{\Sigma }}}_x^{-1},$$

$${\tilde{\mathbf{x}}}_{k+1}=\tilde{\mathbf{A}}{\tilde{\mathbf{x}}}_k.$$

$$\tilde{\mathbf{A}}\mathbf{W}=\mathbf{W\Lambda }.$$

$${\mathbf{x}}_{k+1}$$

$${\mathbf{x}}_{k+1}=\mathbf{A}{\mathbf{x}}_k+\mathbf{B}{\mathbf{\gamma }}_k=\mathbf{G}\left[\begin{array}{c}{\mathbf{x}}_k\\ {\mathbf{\gamma }}_k\end{array}\right],$$

$$\mathbf{G}=\mathbf{Y}{\mathbf{V}}^c({\mathbf{\Sigma }}^c{)}^{-1}{\mathbf{U}}^{c\ast }.$$

$${\mathbf{A}}^c=\mathbf{Y}{\mathbf{V}}^c({\mathbf{\Sigma }}^c{)}^{-1}{\mathbf{U}}^{c_a\ast }\text{and}{\mathbf{B}}^c=\mathbf{Y}{\mathbf{V}}^c({\mathbf{\Sigma }}^c{)}^{-1}{\mathbf{U}}^{c_b\ast }.$$

$${\tilde{\mathbf{x}}}_{k+1}={\tilde{\mathbf{A}}}^c{\tilde{\mathbf{x}}}_k+{\tilde{\mathbf{B}}}^c{\mathbf{\gamma }}_k,$$

$${\tilde{\mathbf{A}}}^c={\bar{\mathbf{U}}}^{\ast }\mathbf{Y}{\bar{\mathbf{V}}}^c({\bar{\mathbf{\Sigma }}}^c{)}^{-1}{\bar{\mathbf{U}}}^{c_a\ast }\bar{\mathbf{U}}\text{and}{\tilde{\mathbf{B}}}^c={\bar{\mathbf{U}}}^{\ast }\mathbf{Y}{\mathbf{V}}^c({\mathbf{\Sigma }}^c{)}^{-1}{\mathbf{U}}^{c_b\ast }.$$

$${\mathbf{x}}_{k+1}$$

$${\mathbf{x}}_{k+1}={\mathbf{A}}_k{\mathbf{x}}_k,$$

$$\mathbf{J}({\mathbf{A}}_k)=\sum _{j=k-w+1}^k{\Vert {\rho }^{(k-j)}({\mathbf{y}}_j-{\mathbf{A}}_k{\mathbf{x}}_j)\Vert }^2,$$

$${\mathbf{x}}_{k+1}$$

$${\mathbf{x}}_{k+1}={\mathbf{A}}_k^c{\mathbf{x}}_k+{\mathbf{B}}_k^c{\mathbf{\gamma }}_k,$$

$${\mathbf{J}}^c({\mathbf{A}}_k^c,{\mathbf{B}}_k^c)=\sum _{j=k-w+1}^k{\Vert {\rho }^{(k-j)}\left({\mathbf{y}}_j-{\mathbf{A}}_k^c{\mathbf{x}}_j-{\mathbf{B}}_k^c{\mathbf{\gamma }}_j\right)\Vert }^2,$$

$${\mathbf{x}}_{k+1}^c={\mathbf{G}}_k{\mathbf{x}}_k^c,$$

$${\mathbf{J}}^c({\mathbf{G}}_k)=\sum _{j=k-w+1}^k{\Vert {\rho }^{(k-j)}\left({\mathbf{y}}_j-{\mathbf{G}}_k{\mathbf{x}}_j^c\right)\Vert }^2.$$

$${\mathbf{X}}_k$$

$${\mathbf{Y}}_k$$

$${\mathbf{A}}_k={\mathbf{Y}}_k{\mathbf{X}}_k^{\ast }{\left({\mathbf{X}}_k{\mathbf{X}}_k^{\ast }\right)}^{-1}={\mathbf{Y}}_k{\mathbf{V}}_{x_k}{\mathbf{\Sigma }}_{x_k}^{-1}{\mathbf{U}}_{x_k}^{\ast },$$

$${\bar{\mathbf{A}}}_k={\mathbf{Y}}_k{\bar{\mathbf{V}}}_{x_k}{\bar{\mathbf{\Sigma }}}_{x_k}^{-1}{\bar{\mathbf{U}}}_{x_k}^{\ast },$$

$${\mathbf{X}}_k^c$$

$${\mathbf{Y}}_k$$

$${\mathbf{G}}_k={\mathbf{Y}}_k{\mathbf{X}}_k^{c\ast }{\left({\mathbf{X}}_k^c{\mathbf{X}}_k^{c\ast }\right)}^{-1}={\mathbf{Y}}_k{\mathbf{V}}_{x_k}^c({\mathbf{\Sigma }}_{x_k}^c{)}^{-1}{\mathbf{U}}_{x_k}^{c\ast },$$

$${\bar{\mathbf{G}}}_k={\mathbf{Y}}_k^c{\bar{\mathbf{V}}}_{x_k}^c({\bar{\mathbf{\Sigma }}}_{x_k}^c{)}^{-1}{\bar{\mathbf{U}}}_{x_k}^{c\ast }$$

$${\bar{\mathbf{A}}}_k^c={\mathbf{Y}}_k{\bar{\mathbf{V}}}_{x_k}^c({\bar{\mathbf{\Sigma }}}_{x_k}^c{)}^{-1}{\bar{\mathbf{U}}}_{x_k}^{c_a\ast },{\bar{\mathbf{B}}}_k^c={\mathbf{Y}}_k{\bar{\mathbf{V}}}_{x_k}^c({\bar{\mathbf{\Sigma }}}_{x_k}^c{)}^{-1}{\bar{\mathbf{U}}}_{x_k}^{c_b\ast }.$$

$${\mathbf{X}}_k=\left[\begin{array}{cccc}{\rho }^{k-1}{\mathbf{x}}_1& {\rho }^{k-2}{\mathbf{x}}_2& \cdots & {\mathbf{x}}_k\end{array}\right],$$

Full text

Keywords: data-driven dynamics, dynamic mode decomposition, time-varying systems, incremental algorithms

1 Introduction

Emergence of low cost sensors and their widespread deployment has led to an unprecedented amount of data. Extraction of actionable information from such plethora of data remains a challenge. It is often the case that a low dimensional dynamical system governs the high-dimensional spatiotemporal sensory data. Dynamic Mode Decomposition (DMD) has emerged as a popular data-driven technique to efficiently compute such low dimensional dynamics [30, 19, 33]. One of the attractive features of the DMD approach is that it requires minimal assumptions on the data. Furthermore, its strong connections with the Koopman operator [29, 23, 24, 25, 17] makes it further appealing and theoretically grounded. The efficacy and simplicity of the DMD has inspired its application in a wide range of areas from fluid dynamics [30, 19] to video processing [10, 32] to epidemiology [28] to neuroscience [3].

For high-dimensional data that is generated by an underlying low dimensional time-varying dynamics, the standard DMD approach may be applicable only locally in time and appropriate time-varying DMD operators can be computed by discounting old data or considering a sliding window of observations [35]. Even in scenarios when underlying dynamics is time-invariant and nonlinear, such time-varying linear approximations have been shown to be very effective [7]. In this paper, we develop incremental methods for the computation of the DMD for time-varying systems and demonstrate the utility of the developed algorithm using Electroencephalographic (EEG) data.

Incremental approaches to DMD refer to methods that allow for efficiently computing the DMD when data is provided sequentially instead of being provided as a batch. The sequential arrival might be due to the inherent nature of the application or it may be used for computational efficiency even if all the data is accessible. One such technique is streaming DMD [13, 12] that first computes the projection of new data on the current orthogonal basis for the range of the data matrix and compares the norm of the difference between the projection and the new data with a threshold. If the difference is larger than the threshold, then it appends the basis with an additional element. Incremental total DMD introduced in [21] applies incremental Singular Value Decomposition (SVD) [1, 2, 26] for incremental computation of dynamic modes and subsequent identification of dominant modes. Unlike the streaming DMD and the total DMD, our focus in this paper is on systems with time-varying dynamics.

Within the context of time-varying systems, an incremental approach is the Online DMD  [35], which adapts the sequential least squares method [16, 15] to incrementally update the DMD computation. The online DMD incrementally computes the DMD operator using the Sherman-Morrison identity [31]. Authors of [35] consider two strategies in their DMD computation to account for time-varying systems: (i) discounting the old data, and (ii) using a sliding window of observations. They also consider DMD for systems with control input [27, 18]. In particular, they extend the DMD with control approach proposed in [27] to develop an online DMD with control algorithm for time-varying system.

In this paper, we build upon the work in [35], and develop incremental SVD based techniques for the computation of the online DMD for time-varying systems. The proposed approach provides access to singular values of the data matrix at each time and hence allows it to determine a reduced order model. For the scenarios in which the high-dimensional data is governed by a low-dimensional dynamics, such reduced order models provide better future prediction accuracy. Furthermore, for scenarios where the Sherman-Morrison identity cannot be applied, e.g., if the data matrix remains singular even after sufficiently long time, the access to singular values allows for efficient computation of appropriate pseudo-inverse. The major contributions of this work are fourfold

1. (i).

We leverage incremental SVD techniques to develop two algorithms for the computation of the online DMD for time-varying systems that rely on discounting old data and a sliding window of data, respectively; 2. (ii).

We extend both these algorithms to develop incremental SVD based algorithms for the computation of the online DMD with control input for time-varying systems; 3. (iii).

We apply our algorithms for the case without control input on an EEG dataset and show their efficacy in reconstructing and predicting error related potential, a slow cortical potential seen in the EEG signal that is elicited by an unexpected outcome; 4. (iv).

We apply our algorithms for the case with control input on randomly generated linear time-varying dynamical systems, and demonstrate their efficacy.

The remainder of the paper is organized as follows. In Section 2, we present some background on the DMD for time-invariant and time-varying systems. In Section 3, we present some background for the incremental SVD algorithms. In Section 4, we leverage the incremental SVD algorithms to develop novel incremental DMD algorithms. We apply these algorithms to EEG data in Section 5. We present incremental SVD-based algorithms for online DMD with control and illustrate its utility using numerical examples in Section 6. Finally, we conclude in Section 7.

Notation: We denote the set of real and complex matrices of size $n\times m$ by $\mathbb{R}^{n\times m}$ and $\mathbb{C}^{n\times m}$, respectively. The set of real and complex $n$-dimensional vectors are denoted by $\mathbb{R}^{n}$ and $\mathbb{C}^{n}$, respectively. We denote matrices with bold upper case letters, vectors with bold lower case letters, and scalar with normal lower case letters. We denote the matrix transpose and the matrix conjugate transpose by $(\cdot)^{\text{T}}$ and $(\cdot)^{*}$, respectively. The kernel and range of a matrix $\mathbf{X}$ are denoted by $\text{Ker}\mathbf{X}$, and $\text{Ran}\mathbf{X}$, respectively.

2 Dynamic Mode Decomposition for Time-Varying Systems

In this section, we recall the DMD setup and describe the time-varying DMD problem that we study in this paper.

2.1 Dynamic Mode Decomposition

Consider the following discrete-time system:

[TABLE]

where $\mathbf{x}_{k}\in\mathbb{R}^{n}$ is a high-dimensional state vector ($n\gg 1$) sampled at $t_{k}=k\Delta t$, $k\in\{1,\dots,(m+1)\}$, and $\mathbf{f}$ is an unknown map which describes the evolution of the state vector between two subsequent sampling times. Suppose that the evolution of the high-dimensional state $\mathbf{x}$ is governed by some underlying low-dimensional dynamics. Then, the DMD computes a data-driven linear approximation to the system (1) as follows [19, 33].

Consider a collection of $(m+1)$ sequential measurements arranged in the following two datasets:

[TABLE]

Let the projection of $\mathbf{X}$ onto its leading $r$ singular vectors be $\bar{\mathbf{X}}=\bar{\mathbf{U}}_{x}\bar{\mathbf{\Sigma}}_{x}\bar{\mathbf{V}}_{x}^{*}$ , where $r$ is the targeted dimension of the underlying low-dimensional dynamical system, $\bar{\mathbf{U}}_{x}\in\mathbb{C}^{n\times r}$, $\bar{\mathbf{V}}_{x}\in\mathbb{C}^{r\times m}$, and $\bar{\mathbf{\Sigma}}_{x}=\text{diag}\{\sigma_{1},\dots,\sigma_{r}\}\in\mathbb{C}^{r\times r}$. Then, the DMD framework approximates the dynamics (1) by

[TABLE]

where $\bar{\mathbf{A}}\in\mathbb{R}^{n\times n}$ is called the DMD operator and is given by

[TABLE]

Let the projection of the state $\mathbf{x}$ and the DMD operator $\bar{\mathbf{A}}$ onto the space spanned by leading $r$ singular vectors of $\mathbf{X}$ be $\tilde{\mathbf{x}}_{k}=\bar{\mathbf{U}}^{\ast}\mathbf{x}_{k}$ and

[TABLE]

respectively. Then, the approximation to the underlying $r$-dimensional dynamics is

[TABLE]

Let $\mathbf{\Lambda}$ be the diagonal matrix of the eigenvalues of $\tilde{\mathbf{A}}$ and $\mathbf{W}$ be the matrix of the associated eigenvectors such that

[TABLE]

Then, two types of DMD modes can be identified using the eigen-decomposition in (4): the projected DMD modes $\hat{\mathbf{\Phi}}=\bar{\mathbf{U}}_{x}\mathbf{W}$ and the exact DMD modes $\mathbf{\Phi}=\mathbf{Y}\bar{\mathbf{V}}_{x}\bar{{\mathbf{\Sigma}}}_{x}^{-1}\mathbf{W}$ [33]. Either of these DMD modes and the associated eigenvalues describe the evolution of the high-dimensional system using low-dimensional dynamics (3).

2.2 Dynamic Mode Decomposition with Control

Consider the following non-autonomous discrete-time system:

[TABLE]

where $\mathbf{x}_{k}\in\mathbb{R}^{n}$ is a high-dimensional state vector ($n\gg 1$) sampled at $t_{k}=k\Delta t$, $k\in\{1,\dots,(m+1)\}$, and $f^{c}$ is an unknown time-varying map which describes the evolution of the state vector between two subsequent sampling times, and $\boldsymbol{\gamma}_{k}\in\mathbb{R}^{l}$ is the exogenous input. Then, the Dynamic Mode Decomposition with Control (DMDc) computes a data-driven linear approximation to the system (5) as follows [27].

Consider a collection of $m$ sequential measurements $\mathbf{x}_{i}\in\mathbb{R}^{n}$ and the associated exogeneous inputs $\boldsymbol{\gamma}_{i}\in\mathbb{R}^{l}$, $i\in\{1,\dots,m\}$. DMDc approximates the non-autonomous dynamics underlying these measurements by

[TABLE]

where $\mathbf{A}\in\mathbb{R}^{n\times n}$ is called the DMD operator, $\mathbf{B}\in\mathbb{R}^{n\times l}$ is called the input matrix, and $\mathbf{G}=\begin{bmatrix}\mathbf{A}&\mathbf{B}\end{bmatrix}\in\mathbb{R}^{n\times(n+l)}$. In particular, the DMDc algorithm constructs matrices $\mathbf{X}=\begin{bmatrix}\mathbf{x}_{1}&\cdots&\mathbf{x}_{m-1}\end{bmatrix}$, $\mathbf{Y}=\begin{bmatrix}\mathbf{x}_{2}&\cdots&\mathbf{x}_{m}\end{bmatrix}$, and $\mathbf{\Gamma}=\begin{bmatrix}\boldsymbol{\gamma}_{1}&\cdots&\boldsymbol{\gamma}_{m-1}\end{bmatrix}$. Let $\mathbf{X}^{c}=\begin{bmatrix}\mathbf{X}\\ \mathbf{\Gamma}\end{bmatrix}\in\mathbb{R}^{(n+l)\times(m-1)}$, and $\mathbf{X}^{c}=\mathbf{U}^{c}\mathbf{\Sigma}^{c}\mathbf{V}^{c*}$ be its singular value decomposition, where $\mathbf{U}^{c}\in\mathbb{C}^{(n+l)\times(n+l)}$, $\mathbf{\Sigma}^{c}\in\mathbb{C}^{(n+l)\times(n+l)}$, and $\mathbf{V}^{c}\in\mathbb{C}^{(m-1)\times(n+l)}$.

The DMDc algorithm estimates $\mathbf{G}$ by

[TABLE]

Finally, $\mathbf{U}^{c*}$ can be written in a partitioned form of $\mathbf{U}^{c*}=\begin{bmatrix}\mathbf{U}^{c_{a}*}&\mathbf{U}^{c_{b}*}\end{bmatrix}$, where $\mathbf{U}^{c_{a}}\in\mathbb{C}^{n\times(n+l)}$ and $\mathbf{U}^{c_{b}}\in\mathbb{C}^{l\times(n+l)}$. Then, the DMDc algorithm computes the matrices $\mathbf{A}^{c}$ and $\mathbf{B}^{c}$ given by

[TABLE]

Similar to the case of the DMD, the DMDc algorithm also enables identification of the low-dimensional non-autonomous system underlying the high-dimensional measurements. Assume that the matrix $\mathbf{X}^{c}$ can be approximated by its projection onto the leading $p$ singular vectors by $\bar{\mathbf{X}}^{c}=\bar{\mathbf{U}}^{c}\bar{\mathbf{\Sigma}}^{c}\bar{\mathbf{V}}^{c*}$, where $\bar{\mathbf{U}}^{c}\in\mathbb{C}^{(n+l)\times p}$, $\bar{\mathbf{\Sigma}}^{c}\in\mathbb{C}^{p\times p}$, and $\bar{\mathbf{V}}^{c}\in\mathbb{C}^{(m-1)\times p}$, and the matrix $\mathbf{X}$ can be approximated by its projection onto the leading $r\leq p$ singular vectors by $\bar{\mathbf{X}}=\bar{\mathbf{U}}\bar{\mathbf{\Sigma}}\bar{\mathbf{V}}^{*}$, where $\bar{\mathbf{U}}\in\mathbb{C}^{n\times r}$, $\bar{\mathbf{\Sigma}}\in\mathbb{C}^{r\times r}$, and $\bar{\mathbf{V}}\in\mathbb{C}^{(m-1)\times r}$. Then, a reduced order model can be represented as follows:

[TABLE]

where the lower-dimensional matrices $\tilde{\mathbf{A}}^{c}\in\mathbf{R}^{r\times r}$ and $\tilde{\mathbf{B}}^{c}\in\mathbf{R}^{r\times l}$ can be calculated by:

[TABLE]

2.3 Problem Formulation

In this paper, we study incremental algorithms for the time-varying DMD and DMDc problems defined below.

2.3.1 Time-varying DMD

Consider the time-varying discrete time system of the form

[TABLE]

where $\mathbf{x}_{k}\in\mathbb{R}^{n}$ is a high-dimensional state vector ($n\gg 1$) sampled at $t_{k}=k\Delta t$, $k\in\{1,\dots,(m+1)\}$, and $\mathbf{f}$ is an unknown time-varying map which describes the evolution of the state vector between two subsequent sampling times. Similar to the standard DMD setup, we assume that the high-dimensional dynamics (6) is generated by a low-dimensional time-varying dynamics.

Suppose that, at sampling time $t_{k+1}$, we have access to a collection of $(w+1)$ sequential measurements $\begin{bmatrix}\mathbf{x}_{k-w+1}&\mathbf{x}_{k-w+2}&\mathbf{x}_{k-w+3}&\cdots&\mathbf{x}_{k+1}\end{bmatrix}$. Our objective is to design efficient computational techniques to identify the time-varying DMD operator $\mathbf{A}_{k}$ that approximates the system (6) by the following linear time-varying system

[TABLE]

such that the cost function

[TABLE]

is minimized, where $\rho\in(0,1]$ is a discounting factor.

2.3.2 Time-varying DMDc

Consider the following variant of the system in (6)

[TABLE]

where $f^{c}$ is an unknown time-varying map which describes the evolution of the state vector between two subsequent sampling times, and $\boldsymbol{\gamma}_{k}\in\mathbb{R}^{l}$ is the exogenous input.

Assume that, at sampling time $\mathbf{t}_{k+1}$, two collections of sequential measurements are available, one for the system states $\begin{bmatrix}\mathbf{x}_{k-w+1}&\mathbf{x}_{k-w+2}&\mathbf{x}_{k-w+3}&\cdots&\mathbf{x}_{k+1}\end{bmatrix}$, and the other one for the exogenous inputs $\begin{bmatrix}\boldsymbol{\gamma}_{k-w+1}&\boldsymbol{\gamma}_{k-w+2}&\boldsymbol{\gamma}_{k-w+3}&\cdots&\boldsymbol{\gamma}_{k+1}\end{bmatrix}$. Our objective is to design efficient computational techniques to identify the time-varying DMDc operator $\mathbf{A}^{c}_{k}$ and input matrix $\mathbf{B}^{c}_{k}$ that approximates the system (9) by the following linear time-varying system

[TABLE]

such that the cost function

[TABLE]

is minimized, where $\rho\in(0,1]$ is a discounting factor.

Let $\mathbf{G}^{*}_{k}=\begin{bmatrix}\mathbf{A}_{k}^{c*}&\mathbf{B}_{k}^{c*}\end{bmatrix}$, and $\mathbf{x}^{c*}_{j}=\begin{bmatrix}\mathbf{x}_{j}^{*}&\boldsymbol{\gamma}_{j}^{*}\end{bmatrix}$. Then, the linear model (10) can be written as

[TABLE]

and the cost function (11) can be equivalently written as

[TABLE]

Online DMD approach [35] also minimizes the cost functions (8) and (11), and computes the time-varying DMD and DMDc operators. In addition to the computation of these operators, the computational techniques developed in this paper also enable efficient computation of (i) the associated DMD eigenvalues and modes, and (ii) the linear time-varying systems that approximate the low-dimensional dynamics underlying (6) and (9), respectively.

2.4 Block Computation of Time-varying DMD and DMDc

We refer to a computation that requires all the data until the sampling time $t_{k}$ to compute the desired solution at time $t_{k}$ as a block computation. In contrast, an incremental computation uses only the data at time $t_{k}$ and the solution at time $t_{k-1}$ to compute the desired solution.

In this section, we describe block computation for time-varying DMD and DMDc in Lemma 1 and 2, respectively. These lemmas are immediate consequence of least square solutions of a linear system of equations [22, Section 5.13] and we omit their proofs. We will use these block computations to derive the incremental computations later in the paper.

Lemma 1**.**

Consider a sequence of $(w+1)$ measurements $\{\mathbf{x}_{k-w+1},\ldots,\mathbf{x}_{k+1}\}$ at sampling time $t_{k+1}$ that is arranged in the following two matrices

[TABLE]

where $\mathbf{x}_{j}\in\mathbb{R}^{n}$, for each $j\in\{k-w+1,\cdots,k+1\}$, and $\mathbf{y}_{j}=\mathbf{x}_{j+1}$. Then, the following statements hold for the cost function (8)

1. (i).

if $\text{Ker}\mathbf{X}_{k}^{*}$ is trivial, then the unique minimizer of (8) is

[TABLE]

where $\mathbf{U}_{x_{k}}\in\mathbb{C}^{n\times n}$ is a unitary matrix, $\mathbf{\Sigma}_{x_{k}}=\text{diag}\{\sigma_{1},\dots,\sigma_{n}\}\in\mathbb{C}^{n\times n}$, $\mathbf{V}_{x_{k}}\in\mathbb{C}^{k\times n}$ has orthonormal columns, and are defined by the reduced SVD of $\mathbf{X}_{k}=\mathbf{U}_{x_{k}}\mathbf{\Sigma}_{x_{k}}\mathbf{V}_{x_{k}}^{*}$. 2. (ii).

if $\text{Ker}\mathbf{X}_{k}^{*}$ is non-trivial, then there exists infinitely many minimizers of (8), and the unique minimizer with the smallest induced two-norm is

[TABLE]

where $\bar{\mathbf{U}}_{x_{k}}\in\mathbb{C}^{n\times r}$, $\bar{\mathbf{\Sigma}}_{x_{k}}=\text{diag}\{\sigma_{1},\dots,\sigma_{r}\}\in\mathbb{C}^{r\times r}$, $\bar{\mathbf{V}}_{x_{k}}\in\mathbb{C}^{k\times r}$, and $r=n-\text{dim}\left(\text{Ker}\mathbf{X}_{k}^{*}\right)$, are defined by the reduced SVD of $\mathbf{X}_{k}\approx\bar{\mathbf{U}}_{x_{k}}\bar{\mathbf{\Sigma}}_{x_{k}}\bar{\mathbf{V}}_{x_{k}}^{*}$.

Lemma 2**.**

Consider two sequences of $(w+1)$ measurements $\{\mathbf{x}_{k-w+1},\ldots,\mathbf{x}_{k+1}\}$, and $\{\boldsymbol{\gamma}_{k-w+1},\cdots,\boldsymbol{\gamma}_{k+1}\}$ at sampling time $t_{k+1}$ that are arranged in the following two matrices

[TABLE]

where $\mathbf{y}_{j}=\mathbf{x}_{j+1}$, $\mathbf{x}^{c*}_{j}=\begin{bmatrix}\mathbf{x}_{j}^{*}&\boldsymbol{\gamma}_{j}^{*}\end{bmatrix}\in\mathbf{R}^{(n+l)}$. Then the following statements hold for the cost function in (13)

1. (i).

if $\text{Ker}\mathbf{X}_{k}^{c*}$ is trivial, then the unique minimizer of (11) is

[TABLE]

where $\mathbf{U}^{c}_{x_{k}}\in\mathbb{C}^{(n+l)\times(n+l)}$ is a unitary matrix, $\mathbf{\Sigma}^{c}_{x_{k}}=\text{diag}\{\sigma_{1},\dots,\sigma_{n+l}\}\in\mathbb{C}^{(n+l)\times(n+l)}$, $\mathbf{V}^{c}_{x_{k}}\in\mathbb{C}^{k\times(n+l)}$ has orthonormal columns, and are defined by the reduced SVD of $\mathbf{X}^{c}_{k}=\mathbf{U}^{c}_{x_{k}}\mathbf{\Sigma}^{c}_{x_{k}}\mathbf{V}_{x_{k}}^{c*}$. 2. (ii).

if $\text{Ker}\mathbf{X}_{k}^{*}$ is non-trivial, then there exists infinitely many minimizers of (11), and the unique minimizer with the smallest induced two-norm is

[TABLE]

where $\bar{\mathbf{U}}^{c}_{x_{k}}\in\mathbb{C}^{(n+l)\times r}$, $\bar{\mathbf{\Sigma}}^{c}_{x_{k}}=\text{diag}\{\sigma_{1},\dots,\sigma_{r}\}\in\mathbb{C}^{r\times r}$, $\bar{\mathbf{V}}^{c}_{x_{k}}\in\mathbb{C}^{k\times r}$, and $r=(n+l)-\text{dim}\left(\text{Ker}\mathbf{X}_{k}^{c*}\right)$, are defined by the reduced SVD of $\mathbf{X}^{c}_{k}\approx\bar{\mathbf{U}}^{c}_{x_{k}}\bar{\mathbf{\Sigma}}^{c}_{x_{k}}\bar{\mathbf{V}}_{x_{k}}^{c*}$.

Note that using the partitioned form of $\mathbf{U}_{k}^{c*}=\begin{bmatrix}\mathbf{U}_{k}^{c_{a}*}&\mathbf{U}_{k}^{c_{b}*}\end{bmatrix}$, where $\mathbf{U}_{k}^{c_{a}}\in\mathbb{C}^{n\times r}$ and $\mathbf{U}_{k}^{c_{b}}\in\mathbb{C}^{l\times r}$, the matrices $\bar{\mathbf{A}}_{k}^{c}$ and $\bar{\mathbf{B}}_{k}^{c}$ associated with (19) can be obtained as follows

[TABLE]

Similar to the online DMD [35], we focus on two special cases of the cost functions (8) and (13), namely, the weighted cost functions and the windowed cost functions. The weighted cost functions consider a gradual elimination for the old measurements by assigning the discounting factor to be $\rho\in(0,1)$ and consider all the available data, i.e., they select $w=k$. The windowed cost functions consider a sharp cut-off window and use a non-discounted window of recent measurements, i.e., they select $\rho=1$ and $w<k$. In the remainder of the paper, we will refer to the DMD and the DMDc operators/matrices obtained using the weighted (resp., windowed) cost functions by the weighted (resp., windowed) DMD and DMDc operator/matrices and denote them by $\mathbf{A}_{k}^{\rho}$ (resp., $\mathbf{A}_{k}^{w}$) and $(\mathbf{A}_{k}^{\rho_{c}},\mathbf{B}_{k}^{\rho_{c}})$ (resp., $(\mathbf{A}_{k}^{w_{c}},\mathbf{B}_{k}^{w_{c}})$), respectively.

3 Background on Incremental SVD Algorithms

At the heart of the DMD framework is the utilization of the SVD to identify an invariant subspace for the low-dimensional dynamics underlying the high-dimensional data. The techniques proposed in this paper for the computation of the time-varying DMD operator require access to the SVD of certain time-varying data matrix at each time. For efficient computation of these SVDs, we resort to incremental SVD techniques proposed in [1, 2, 26]. In this section, we present these incremental SVD techniques. The presentation below is adapted from [1] in order to facilitate better exposition in Section 4.

3.1 Weighted Incremental SVD Algorithm

Consider a weighted dataset

[TABLE]

where $0\leq\rho\leq 1$ is a scalar discounting factor and $k\geq n$. Given an additional datum $\mathbf{x}_{k+1}$, let

[TABLE]

be the incremented dataset. Provided that the SVD of $\mathbf{X}_{k}$ is known, the SVD of $\mathbf{X}_{k+1}$ can be incrementally computed using the following proposition, which is proved in Appendix A.

Proposition 1**.**

Let at sampling time $t_{k}$ the SVD of the dataset in (21) be defined by $~{}\boldsymbol{X}_{k}=\mathbf{U}_{x_{k}}\mathbf{\Sigma}_{x_{k}}\mathbf{V}_{x_{k}}^{*}$, where ${\mathbf{U}_{x_{k}}}\in\mathbb{C}^{n\times n}$ is a unitary matrix, ${\mathbf{V}_{x_{k}}}\in\mathbb{C}^{k\times n}$ has orthonormal columns, and $\mathbf{\Sigma}_{x_{k}}=\text{diag}\{\sigma_{1},\dots,\sigma_{n}\}\in\mathbb{C}^{n\times n}$ is a diagonal matrix. Assume that at sampling time $t_{k+1}$, a new datum $\mathbf{x}_{k+1}$ is accessed. Then, the SVD for the dataset in (22), defined by $~{}\mathbf{X}_{k+1}=\mathbf{U}_{x_{k+1}}\mathbf{\Sigma}_{x_{k+1}}\mathbf{V}^{*}_{x_{k+1}}$, is given by :

[TABLE]

where $\mathbf{U}_{s_{k}}$, $\mathbf{\Sigma}_{s_{k}}$, and $\mathbf{V}_{s_{k}}=\begin{bmatrix}\mathbf{V}_{s_{k,1}}\\ \mathbf{v}_{s_{k,2}}\end{bmatrix}$, with $\mathbf{v}_{s_{k,2}}\in\mathbb{C}^{1\times n}$, are defined by the following SVD

[TABLE]

The above incremental SVD computation requires the computation of the SVD of matrix $\mathbf{S}_{k}\in\mathbb{C}^{n\times(n+1)}$. The size of this matrix is smaller than the size of dataset $\mathbf{X}_{k+1}$ and it has the so-called broken arrow structure which enables efficient computation of its SVD [11, 14, 6].

3.2 Windowed Incremental SVD Algorithm

Consider a windowed dataset

[TABLE]

where $n$ is the dimension of measurements and $w$ is the length of the time-window of the desired measurements. Define $q$ as follows

[TABLE]

Consider another dataset

[TABLE]

Then, the incremental SVD computes the SVD of $\boldsymbol{\chi}_{k+1}=\mathbf{U}_{{\chi}_{k+1}}\mathbf{\Sigma}_{{\chi}_{k+1}}\mathbf{V}_{{\chi}_{k+1}}^{*}$ as follows. First, it eliminates $\mathbf{x}_{k-w+1}$ from the dataset to obtain decremented dataset and the associated SVD

[TABLE]

Then, it increments $\boldsymbol{\acute{\chi}}_{k}$ with $\mathbf{x}_{k+1}$ to obtain $\boldsymbol{\acute{\chi}}_{k+1}$ and the associated SVD. The SVD of the dataset $\boldsymbol{\acute{\chi}}_{k}$ can be incrementally computed using the following proposition, which is proved in Appendix B.

Proposition 2**.**

Let at sampling time $t_{k}$, the SVD of the dataset in (24) be defined by $~{}\boldsymbol{\chi}_{k}=\mathbf{U}_{\chi_{k}}\mathbf{\Sigma}_{\chi_{k}}{\mathbf{V}_{\chi_{k}}^{*}}$, where $\mathbf{U}_{\chi_{k}}\in\mathbb{C}^{{n}\times{q}}$, $\mathbf{\Sigma}_{\chi_{k}}=\begin{bmatrix}\text{diag}\{\sigma_{1},\dots,\sigma_{q}\}&\mathbf{0}_{q\times(w-q)}\end{bmatrix}\in\mathbb{C}^{{q}\times{w}}$, and $\mathbf{V}_{\chi_{k}}\in\mathbb{C}^{{w}\times{w}}$ is a unitary matrix. Then, the SVD of the dataset in (27) defined as $~{}\boldsymbol{\acute{\chi}}_{k}=\mathbf{U}_{\boldsymbol{\acute{\chi}}_{k}}\mathbf{\Sigma}_{\boldsymbol{\acute{\chi}}_{k}}\mathbf{V}_{\boldsymbol{\acute{\chi}}_{k}}^{*}$, is given by :

[TABLE]

where $\mathbf{U}_{\acute{s}_{k}}$, $\mathbf{V}_{\acute{s}_{k}}$, and $\mathbf{\Sigma}_{\acute{s}_{k}}$ are defined by the following SVD

[TABLE]

and $\mathbf{V}_{\chi_{k,2}}$ is the submatrix of $\mathbf{V}_{\boldsymbol{{\chi}}_{k}}$ obtained after removing its first row.

In the incremental step, the new datum $\mathbf{x}_{k+1}$ is appended to $\boldsymbol{\acute{\chi}}_{k}$ and the SVD of the resulting matrix $\boldsymbol{\chi}_{k+1}$ can be obtained using the procedure presented in Section 3.1. In particular, let $\mathbf{\hat{S}}_{k}\triangleq\begin{bmatrix}{\mathbf{\Sigma}}_{\acute{\chi}_{k}}&\mathbf{U}_{\acute{\chi}_{k}}^{*}\mathbf{x}_{k+1}\end{bmatrix}=\mathbf{U}_{\hat{s}_{k}}\mathbf{\Sigma}_{\hat{s}_{k}}\mathbf{V}^{*}_{\hat{s}_{k}}$, and using equation (23), the SVD of $\boldsymbol{\chi}_{k+1}=\mathbf{U}_{\chi_{k+1}}\mathbf{\Sigma}_{\chi_{k+1}}\mathbf{V}_{\chi_{k+1}}^{*}$ is

[TABLE]

The windowed incremental SVD computation requires the computation of the SVD of $\mathbf{\acute{S}}_{k}\in\mathbb{R}^{q\times w}$. In general, $\mathbf{\acute{S}}_{k}$ may not possess any special structure. However, if $n\gg w$, in which case $q=w$, the above incremental update maybe a lot cheaper than the computation of the SVD of an $n\times w$ matrix. Also, in this context, the above procedure helps in term of storage as it does not require the large matrix $\boldsymbol{\chi}_{k+1}$ to be stored and uses only the first and last column of this matrix along with the SVD at the previous iteration.

4 Incremental DMD Algorithms

In this section, we present two algorithms, namely weighted incremental DMD and windowed incremental DMD, respectively, for incremental computation of a time-varying lower-dimensional DMD operator. Similar to the existing approaches in the literature, the proposed algorithms perform incremental computations either by assigning decaying weights to the received data [16, 35], or by using a sliding window of the received data [10, 20, 35]. The recursive updates derived in this section are similar to those in sequential least square method [16, Appendix 8C], recursive least square method [15, Section 9.4], as well as online DMD [35].

4.1 Weighted Incremental DMD Algorithm

In this algorithm, the time-varying DMD operator is estimated by assigning a decaying weight to past measurements in order to gradually discount their representation as newer measurements become available. Assume that at sampling time $t_{k+1}$, we have access to the following weighted datasets:

[TABLE]

where $\mathbf{y}_{k}=\mathbf{x}_{k+1}$, for each $k\in\mathbb{N}$. Assume that at sampling time $t_{k+2}$, the datasets in (30) are updated with a pair of measurements $(x_{k+1},y_{k+1})$ such that

[TABLE]

Suppose that at sampling time $t_{k+1}$, the DMD operator $\mathbf{A}^{\rho}_{k}$ and the SVD of $\mathbf{X}_{k}$ are known, then at sampling time $t_{k+2}$, the DMD operator can be updated with the new pair of measurements $(\mathbf{x}_{k+1},\mathbf{y}_{k+1})$ using the following theorem, which is proved in Appendix C.

Theorem 1**.**

Let at sampling time $t_{k+1}$ the SVD of $\mathbf{X}_{k}=\mathbf{U}_{x_{k}}\mathbf{\Sigma}_{x_{k}}\mathbf{V}_{x_{k}}^{*}$ be known, and let the DMD operator $\mathbf{A}^{\rho}_{k}\in\mathbb{R}^{n\times n}$ minimize the cost function in (8) with $\mathbf{Y}_{k}$ and $\mathbf{X}_{k}$ given in (30). Assume that at sampling time $t_{k+2}$, a new pair of measurements $(\mathbf{x}_{k+1},\mathbf{y}_{k+1})$ is used to incrementally compute the SVD of $\mathbf{X}_{k+1}$ using Proposition 1. Then, the DMD operator $\mathbf{A}^{\rho}_{k+1}\in\mathbb{R}^{n\times n}$ which minimizes the cost function in (8) is

[TABLE]

If $\mathbf{X}_{k+1}$ is well-approximated by its projection onto its leading $r$ singular vectors given by $\bar{\mathbf{X}}_{k+1}=\bar{\mathbf{U}}_{x_{k+1}}\bar{\mathbf{\Sigma}}_{x_{k+1}}\bar{\mathbf{V}}^{*}_{x_{k+1}}$, then the DMD operator update in (32) takes the form

[TABLE]

Moreover, the DMD operator $\bar{\mathbf{A}}^{\rho}_{k+1}\in\mathbb{R}^{n\times n}$ can be projected onto a subspace spanned by the leading $r$ left singular vectors of $\mathbf{X}_{k+1}$ to obtain a lower-dimensional DMD operator $\tilde{\mathbf{A}}^{\rho}_{k+1}\in\mathbb{R}^{r\times r}$ given by

[TABLE]

The update in equation (32) requires only the current DMD operator, measurement pair $(\mathbf{x}_{k+1},\mathbf{y}_{k+1})$, and the incremental SVD update. Specifically, these updates do not require the data matrix to be stored.

4.2 Windowed Incremental DMD Algorithm

We now focus on windowed incremental DMD algorithm in which a sliding-window of $w$ most recent measurements is used to estimate the time-varying DMD operator. Assume that we have access to the following windowed datasets at sampling time $t_{k+1}$:

[TABLE]

Accordingly, at sampling time $t_{k+2}$, we have access to

[TABLE]

Suppose that at sampling time $t_{k+1}$, the DMD operator $\mathbf{A}^{w}_{k}$ and the SVD of $\boldsymbol{\chi}_{k}$ are known, then at sampling time $t_{k+2}$, the DMD operator can be updated with the new pair of measurements $(\mathbf{x}_{k+1},\mathbf{y}_{k+1})$ using the following theorem, which is proved in Appendix D.

Theorem 2**.**

Let at sampling time $t_{k+1}$ the SVD of $\boldsymbol{\chi}_{k}=\mathbf{U}_{\chi_{k}}\mathbf{\Sigma}_{\chi_{k}}\mathbf{V}^{*}_{\chi_{k}}$ be known, and let the DMD operator $\mathbf{A}^{w}_{k}\in\mathbb{R}^{n\times n}$ minimize the cost function in (8) with $\boldsymbol{\Upsilon}_{k}$ and $\boldsymbol{\chi}_{k}$ given in (35). Assume that at sampling time $t_{k+2}$, a new pair of measurements $(\mathbf{x}_{k+1},\mathbf{y}_{k+1})$ is used to incrementally compute the SVD of $\boldsymbol{\chi}_{k+1}$ using Proposition 2 and equation (3.2). Then, the DMD operator $\mathbf{A}^{w}_{k+1}\in\mathbb{R}^{n\times n}$ which minimizes the cost function in (8) is

[TABLE]

If $\boldsymbol{\chi}_{k+1}$ is well-approximated by its projection onto its leading $r$ singular vectors given by $\bar{\boldsymbol{\chi}}_{k+1}=\bar{\mathbf{U}}_{\chi_{k+1}}\bar{\mathbf{\Sigma}}_{\chi_{k+1}}\bar{\mathbf{V}}^{*}_{\chi_{k+1}}$, then the DMD operator takes the form

[TABLE]

Moreover, the DMD operator $\bar{\mathbf{A}}^{w}_{k+1}\in\mathbb{R}^{n\times n}$ can be projected onto a subspace spanned by the leading $r$ left singular vectors of $\boldsymbol{\chi}_{k+1}$ to obtain a lower-dimensional DMD operator ${\tilde{\mathbf{A}}}^{w}_{k+1}\in\mathbb{R}^{r\times r}$ given by

[TABLE]

Similar to the weighted incremental DMD, the update in equation (37) requires only the current DMD operator, measurement pair $(\mathbf{x}_{k+1},\mathbf{y}_{k+1})$, and the windowed incremental SVD update.

For a fixed dimension $r$ of the reduced order system, the computational complexity of rank one incremental SVD update is $\mathcal{O}(nr+r^{3})$ [2, Section 4.3]. In comparison, computational complexities of each iteration in the streaming DMD and the online DMD are $\mathcal{O}(nr^{2})$ and $\mathcal{O}(n^{2})$, respectively. When data matrices are non-singular, both online and incremental DMD approaches can be used to compute time-varying DMD operators. In such cases, if a reduced order DMD operator of dimension $r$ is of interest, then the proposed incremental DMD approach has a better computational complexity, when $r<n^{2/3}$.

5 Dynamic Mode Decomposition for Error-Related Potentials in EEG data

Error‐ Related Potentials (ErrPs) are slow cortical potentials seen in the EEG signal of human subjects that are elicited by an unexpected (erroneous) outcome. For example, within the context of human-machine interaction, such signals are observed when the machine takes an unexpected action [9, 4, 5]. In this section, we illustrate the efficacy of the proposed DMD algorithms on EEG recordings taken from Monitoring Error-Related Potentials database [4, 5]. These EEG recordings were collected from six subjects while they were observing the movement of a cursor between two labeled targets on a screen. The subjects had no control on the motion of the cursor. However, they had a priori knowledge about the intended directions of movement. The cursor movements were set to elicit ErrPs by generating two types of events, namely, correct events and erroneous events. The correct events correspond to the cursor movement in the intended direction, while the erroneous events correspond to the cursor movement in any unintended direction. EEG signals were recorded at sampling rate of $512$ Hz, using Biosemi ActiveTwo system with $64$ electrodes distributed according to the standard $10/20$ international system as it is shown in Figure 1.

A two step standard EEG preprocessing for ErrPs which consists of Common Average Reference (CAR) filtering, and then $(1-10)$ Hz Band Pass (BP) filtering is performed on the data. Segments of EEG signals, from all the participants, corresponding to the correct and erroneous events were extracted into two separate groups. Then, the segments within each group were averaged to obtain a signal each for the correct and erroneous events. Finally, the ErrPs are calculated by subtracting the signal for correct event from the signal for the erroneous event. Figure 2 shows the ERP signals and topographical views corresponding to the correct and erroneous events, and the calculated ErrPs. EEGLAB was used for the processing of EEG signals and generating the topographical views for the related brain activities [8]. Note that for consistency with the literature, we investigate the averaged ERP signal. However, the confidence regions around the averaged signal shown in Figure 2 reveal that despite the variability in these signals, the key features are also seen in the signal associated with a single event. Thus, the following analysis can be performed on a single subject real time data as well.

As reported in [4, 5], the ErrPs are characterized by a sequence of three peaks/troughs after the onset of the event in the EEG signal measured at channel FCz (see Figure 1). A sample ErrP is shown in Figure 2c where the first small positive peak is observed at 200 ms, followed by a negative trough at 260 ms, and finally a larger positive peak at 330 ms.

Two observations can be made by comparing the temporal behavior and the topographical views during the correct and erroneous events in Figures 2a and 2b, respectively. First, the ERP signal during the correct event tends to have broader peaks in comparison with the ERP signal during erroneous events, which indicates that it has slower dynamics, i.e., has smaller growth or decay rates. Second, by comparing the topographical views during the three characterizing peaks at FCz channel, it can be noticed that the ERP peaks during erroneous events have higher amplitudes than the peaks during correct events in both negative and positive directions. Note that the colormap bar values are different for the topographical views for the correct ERP, the erroneous ERP, and the ErrP in Figure 2.

As shown in Figure 2, the EEG patterns show coherent activities within the brain suggesting the possibility of a low-dimensional dynamics underlying this high dimensional data. Furthermore, these patterns change significantly with time and can possibly be better explained using a time-varying underlying model. To investigate these hypotheses, we apply the incremental DMD algorithms to the EEG recordings corresponding to correct and erroneous events. The weighted incremental DMD algorithm is applied with discounting factors $\rho=\left\{0.1,0.2,0.4,0.8\right\}$, and the windowed incremental DMD algorithm with a time-window of width $w=512$ samples. All singular values greater than $\sigma_{\text{thr}}=\{0.01,0.001\}$ are used for computing the reduced order DMD operator. An initial window of $512$ samples before the event occurrence is used to initialized the DMD models for both algorithms.

We compare the proposed algorithms with the online DMD approach proposed in [35]. To this end, we use the same parameters (discounting factor and width of time-window) in the online DMD algorithms as the incremental DMD algorithms. Recall that the online DMD algorithms rely on the Sherman-Morrison identity to compute the DMD operator and do not have access to the (time-varying) singular values of the data matrix and hence, the reduction of the DMD operator to leading singular vectors does not apply to online DMD algorithms.

In order to compute the time-varying DMD operator, we first compute a standard DMD operator using EEG data from a one-second window ($512$ samples) just before the correct or erroneous event. Subsequently, the EEG measurements are streamed with a rate of $1$ $\text{sample}/\text{iteration}$ to update the initial DMD operator. At each iteration, both weighted and windowed versions of the incremental DMD and online DMD algorithms are used to predict $64$ future samples of the EEG signal at channel FCz. We compare the performance of these algorithms using normalized Root Mean Square (RMS) prediction error, denoted as ${e}_{\text{nrms}}(k)$ at each $t_{k}\in\{0,\cdots,0.6\}$ using the following formula:

[TABLE]

where $\hat{\mathbf{y}}_{i}$ and $\mathbf{y}_{i}$ are the estimated and recorded values, respectively, $\overline{\mathbf{y}}_{k}=\max\limits_{k+1\leq i\leq k+64}\mathbf{y}_{i}$, and $\underline{\mathbf{y}}_{k}=\min\limits_{k+1\leq i\leq k+64}\mathbf{y}_{i}$.

Figure 3 shows the mean of the normalized RMS prediction error computed over all iterations as well as the associated 95% confidence sets for weighted and windowed incremental DMD algorithms and different choices of parameters. The performance of the incremental DMD does not appear to vary much with the weighting factor and the two choices of the threshold on singular values. The performance of the windowed incremental DMD algorithm is similar to the weighted incremental DMD. However, incremental DMD algorithms seem to outperform online DMD algorithms. This suggests that having a small threshold $\sigma_{\textup{thr}}$ is more beneficial than having no threshold as in the case of online DMD algorithms. Figure 13 in Appendix E presents the evolution of the normalized RMS prediction error for each event.

In order to gain further insight into the influence of different thresholds $\sigma_{\text{thr}}$ and weighting factors $\rho$, we investigate how accurately the incremental DMD predicts the signature ERPs at channel FCz. To this end, we updated the initial DMD operator using incremental DMD algorithms until the first peak at 200 msec. We then used the updated lower-dimensional DMD operator to predict the future ERP until the third peak at 330 msec. The predicted signal and the associated normalized RMS are shown in Figures 4 and 5, respectively. For the erroneous event, the predicted signal from incremental DMD algorithms capture the trend of the original signal, while for the correct event, most incremental DMD algorithms fail to capture the final dip in the original signal. The only exception is the windowed incremental DMD with $\sigma_{\textup{thr}}=0.001$. For $\sigma_{\textup{thr}}=0.01$, the number of DMD modes used by the incremental DMD algorithms range from 26 to 31 for correct events and 28 to 32 for erroneous events. Similarly, for $\sigma_{\textup{thr}}=0.001$, the number of DMD modes used by the incremental DMD algorithms range from 26 to 31 for correct events and 28 to 33 for erroneous events. In comparison, the predicted signal from the online DMD algorithms performs poorly. The key reason for this poor performance, is that the data matrix obtained after prepossessing is not full rank even if all the measurements are included. Thus, the Sherman-Morrison update cannot be applied in the standard form. Recall that the online DMD algorithm starts by initializing the DMD operator $\mathbf{A}_{k}$ and the inverse of covariance matrix $\mathbf{P}_{k}=\left(\mathbf{X}_{k}\mathbf{X}_{k}^{\text{T}}\right)^{-1}$. For initialization of online DMD algorithms, we adopted the suggestion in [35]: in the event that $\mathbf{X}_{k}\mathbf{X}_{k}^{\text{T}}$ is not full rank, $\mathbf{A}_{k}$ can be initialized by an $n\times n$ zero matrix, and $\mathbf{P}_{k}$ can be initialized by $\mathbf{P}_{\text{init}}=\alpha\text{I}$ where $\alpha$ is a large positive scalar.

During the correct events, there appears to be an interplay between the threshold $\sigma_{\text{thr}}$ and weighting factor $\rho$. As shown in Figure 5 for correct events, the normalized RMS decreases with decreasing the weighting factor until $\rho=0.2$ and increases with further decrease in $\rho$. For $\rho=0.2$, the higher dimension of the reduced system ($\sigma_{\textup{thr}}=0.001$) leads to higher prediction error compared with the lower dimension reduced system ($\sigma_{\textup{thr}}=0.01$). This suggests that for a weighting factor that leads to the smallest prediction error, including too many modes is not beneficial for prediction performance, since it leads to over-fitting. A similar interplay is not seen for the erroneous event in Figure 5. However, if we choose a smaller initial window to initialize the DMD model, we observe similar effects as shown in Figure 14 in Appendix F. Similar effects are observed for the windowed DMD as shown in Figure 15 in Appendix F.

To further compare the online DMD with the incremental DMD, we took the raw EEG data, i.e., data without any preprocessing (CAR and BP filtering). In this case, the data matrix is well-conditioned and the Sherman-Morrison identity is well defined. The performance of online DMD and the incremental DMD in terms of predicting future signal is shown in Figure 6. In this case, both incremental and online DMD algorithms have similar poor performance. This suggests that the difference in the performance of these algorithms is primarily due to ill-conditioned dataset. This also highlights the utility of applying appropriate band-pass filtering to the raw EEG data. Without such filtering, the key activities seem to be lost in the background noise resulting in a poor prediction performance.

Topographical views for the real part of the four dominant DMD modes at 200 msec obtained from the incremental DMD algorithms for the correct and erroneous events are shown in Figure 7 and Figure 8, respectively. A comparison between the DMD modes during the erroneous and correct events shows that there is a stronger activity in the frontal lobe during the erroneous events especially around the FCz channel. This is consistent with the topographical views from EEG recordings shown in Figures 2a and 2b that show stronger activity in the frontal lobe during an erroneous event. This illustrates that a few principal DMD modes are able to capture the dominant activity in the brain during these experiments.

Figure 9 shows the logarithm of eigenvalues of the DMD operator at 200 msec associated with different thresholds $\sigma_{\text{thr}}$, weight factor $\rho$ and time-window width. For smaller value of $\sigma_{\textup{thr}}$, the range of eigenvalues in the left half plane is bigger. This is consistent with the fact that smaller $\sigma_{\text{thr}}$ implies a larger reduced order model in which some states converge to zero much faster than other states. For the erroneous event, some eigenvalues are spread in the right half plane, which suggest a locally unstable dynamics underlie the evolution of ERP signal during erroneous trials. This is also consistent with the faster dynamics (sharper peaks) in the EEG signal during the erroneous event (see Figure 2).

We now investigate the effectiveness of the incremental DMD algorithms in terms of reconstructing the ERP signals during the correct and erroneous events at FCz channel. To this end, we reconstruct the signal starting from the onset of the events until the end of the third peak at $330$ msec. The reconstructed signals and the normalized RMS error for the reconstruction are shown in Figures 10 and 11, respectively. These results indicate that the incremental DMD models incur smaller reconstruction error than the online DMD models. This result is counter-intuitive, since online DMD uses more modes and should lead to smaller reconstruction error. This discrepancy happens due to the fact that the data matrix never achieves full rank and the online DMD is not able to overcome the heuristic initialization discussed above. There is no significant difference in the reconstruction error for incremental DMD models with different $\rho$ and $\sigma_{\text{thr}}$ values as shown in Figure 10a.

Finally, we summarize the above investigation into the utility of incremental and online DMD in modeling EEG data. It appears that the pre-processing of the EEG data that requires domain-specific knowledge such as the band of frequency in which the event of interest is observed is vital to obtain sensible DMD-based data-driven models. However, after such pre-processing, the data matrix may become ill-conditioned and under such scenarios incremental DMD techniques proposed in this paper appear more promising that the online DMD techniques. There exists a trade-off between the threshold on the singular values used for model reduction and the window-width or the discount factor used in the incremental DMD algorithms. We also observed that the principal DMD modes are consistent with the activity in the brain during the studied experiments. This suggests that the incremental DMD algorithms can compute an efficient basis for describing the evolution of the EEG activity. Finally, we observed that under erroneous events, some DMD eigenvalues moved towards the right-half complex plane suggesting that certain events trigger (small time) unstable dynamics.

6 Incremental Dynamic Mode Decomposition for Systems with Control Input

In this section, we extended the incremental DMD algorithms to the case of non-autonomous dynamical systems. We first recall the Dynamic Mode Decomposition with Control (DMDc) algorithm proposed in [27] that estimates time-invariant non-autonomous system dynamics underlying high-dimensional data.

For the scenarios in which the non-autonomous system underlying the measurements is time-varying, the DMDc algorithm has been extended to the online DMDc algorithm [35]. The online DMDc algorithm operates similarly to the online DMD algorithm and yields time-varying matrices $\mathbf{A}^{c}_{k}\in\mathbb{R}^{n\times n}$ and $\mathbf{B}^{c}_{k}\in\mathbb{R}^{n\times l}$, $k\in\{1,2,\cdots\}$. Similar to the online DMD algorithm, the DMDc algorithm also does not allow for incremental computation of the lower dimensional non-autonomous system underlying high dimensional measurements. In the following, we extend the DMDc algorithm to the weighted and windowed incremental DMDc algorithms that enable us to obtain a lower-dimensional time-varying approximation for the underlying dynamics.

6.1 Weighted Incremental DMDc Algorithm

Assume that at sampling time $t_{k+1}$, the measurements and the associated exogeneous inputs are arranged in the following datasets

[TABLE]

respectively, where $\mathbf{X}_{k},\mathbf{Y}_{k}\in\mathbb{R}^{n\times k}$, $\mathbf{\Gamma}_{k}\in\mathbb{R}^{l\times k}$, and $\mathbf{X}_{k}^{c}\in\mathbb{R}^{(n+l)\times k}$. At sampling time $t_{k+1}$, the datasets are updated with a new set of measurements $\left(\mathbf{x}_{k+1},\mathbf{y}_{k+1},\boldsymbol{\gamma}_{k+1}\right)$ which results in the following

[TABLE]

Suppose that the weighted DMDc operator and input matrix at $t_{k+1}$, and the SVD of both $\mathbf{X}_{k}$ and $\mathbf{X}_{k}^{c}$ are known, then according to the new incoming measurements at $t_{k+2}$, the updates for the weighted DMDc operator and input matrix can be obtained using the following theorem, which is proved in Appendix G.

Theorem 3**.**

Let at sampling time $t_{k+1}$ the SVD of $\mathbf{X}^{c}_{k}=\mathbf{U}^{c}_{x_{k}}\mathbf{\Sigma}^{c}_{x_{k}}\mathbf{V}^{c*}_{x_{k}}$ and $\mathbf{X}_{k}=\mathbf{U}_{x_{k}}\mathbf{\Sigma}_{x_{k}}\mathbf{V}^{*}_{x_{k}}$ be known. Assume that the pair $(\mathbf{A}^{\rho_{c}}_{k},\mathbf{B}^{\rho_{c}}_{k})$ minimizes the cost function in (11) with $\mathbf{X}_{k}$ and $\mathbf{\Gamma}_{k}$ defined in (41). Consider that at sampling time $t_{k+2}$, a new triplet of measurements $(\mathbf{x}_{k+1},\mathbf{y}_{k+1},\boldsymbol{\gamma}_{k+1})$ is used to incrementally compute the SVD of both $\mathbf{X}_{k+1}^{c}$ and $\mathbf{X}_{k+1}$ using Proposition 1. Then, the pair $(\mathbf{A}^{\rho_{c}}_{k+1},\mathbf{B}^{\rho_{c}}_{k+1})$ which minimizes the cost function in (11) is

[TABLE]

Assume that $\mathbf{X}^{c}_{k+1}$ can be well-approximated by its projection onto its leading $p$ singular vectors such that $\mathbf{X}^{c}_{k+1}\approx\bar{\mathbf{X}}^{c}_{k+1}=\bar{\mathbf{U}}^{c}_{x_{k+1}}\bar{\mathbf{\Sigma}}^{c}_{x_{k+1}}\bar{\mathbf{V}}^{c*}_{x_{k+1}}$, where $\bar{\mathbf{U}}_{x_{k+1}}^{c}\in\mathbb{C}^{(n+l)\times p}$, $\bar{\mathbf{\Sigma}}_{x_{k+1}}^{c}\in\mathbb{C}^{p\times p}$ ,and $\bar{\mathbf{V}}_{x_{k+1}}^{c}\in\mathbb{C}^{(k+1)\times p}$. Let $\bar{\mathbf{U}}_{x_{k+1}}^{c*}$ be partitioned such that $\bar{\mathbf{U}}_{x_{k+1}}^{c*}=[\bar{\mathbf{U}}^{c_{a}*}_{x_{k+1}}\bar{\mathbf{U}}^{c_{b}*}_{x_{k+1}}]$, where $\bar{\mathbf{U}}^{c_{a}}_{x_{k+1}}\in\mathbb{C}^{n\times p}$ and $\bar{\mathbf{U}}^{c_{b}}_{x_{k+1}}\in\mathbb{C}^{l\times p}$. Then, the update in (43) yields the following updates for $\bar{\mathbf{A}}^{\rho_{c}}_{k+1}\in\mathbb{R}^{n\times n}$ and $\bar{\mathbf{B}}^{\rho_{c}}_{k+1}\in\mathbb{R}^{n\times l}$:

[TABLE]

Assume that $\mathbf{X}_{k+1}$ is well-approximated by its projection onto its leading $r$ singular vectors such that $\mathbf{X}_{k+1}\approx\bar{\mathbf{X}}_{k+1}=\bar{\mathbf{U}}_{x_{k+1}}\bar{\mathbf{\Sigma}}_{x_{k+1}}\bar{\mathbf{V}}^{*}_{x_{k+1}}$ with $r\leq p$. Then a reduced order model can be represented as:

[TABLE]

where the lower-dimensional approximation $\tilde{\mathbf{A}}^{\rho_{c}}_{k+1}\in\mathbb{R}^{r\times r}$ and $\tilde{\mathbf{B}}^{\rho_{c}}_{k+1}\in\mathbb{R}^{r\times l}$ can obtained by projecting $\bar{\mathbf{A}}^{\rho_{c}}_{k+1}$ and $\bar{\mathbf{B}}^{\rho_{c}}_{k+1}$ onto a subspace spanned by the columns of $\bar{\mathbf{U}}_{x_{k+1}}$ as follows:

[TABLE]

The update in equation (43) requires only the current DMD operator $(\mathbf{A}^{\rho_{c}}_{k},\mathbf{B}^{\rho_{c}}_{k})$, measurements set of $(\mathbf{y}_{k+1},\mathbf{x}_{k+1},{\gamma}_{k+1})$, and the incremental SVD update. Specifically, these updates do not require the large data matrix to be stored.

6.2 Windowed Incremental DMDc Algorithm

Assume that at sampling time $t_{k+1}$, the past $w$ states and control input measurements are arranged in the following datasets

[TABLE]

where $\boldsymbol{\chi}_{k},\boldsymbol{\Upsilon}_{k}\in\mathbb{R}^{n\times w}$, $\mathbf{\Gamma}_{l}\in\mathbb{R}^{l\times w}$, and $\boldsymbol{\chi}^{c}_{k}\in\mathbb{R}^{(n+l)\times w}$, respectively. At sampling time $t_{k+1}$, the datasets are updated by adding a new set of measurements $\left(\mathbf{x}_{k+1},\mathbf{y}_{k+1},\boldsymbol{\gamma}_{k+1}\right)$, and removing the oldest set of measurements $\left(\mathbf{x}_{k-w+1},\mathbf{y}_{k-w+1},\boldsymbol{\gamma}_{k-w+1}\right)$, such that:

[TABLE]

Suppose that the windowed DMDc operator and input matrix at $t_{k+1}$, and the SVD of both $\boldsymbol{\chi}_{k}$ and $\boldsymbol{\chi}^{c}_{k}$ are known. Then, according to the new incoming measurements at $t_{k+2}$, the updates for the windowed DMDc operator and input matrix can be obtained using the following theorem, which is proved in Appendix H.

Theorem 4**.**

Let at sampling time $t_{k+1}$ the SVD of $\boldsymbol{\chi}_{k}=\mathbf{U}_{\chi_{k}}\mathbf{\Sigma}_{\chi_{k}}\mathbf{V}^{*}_{\chi_{k}}$ and $\boldsymbol{\chi}_{k}^{c}=\mathbf{U}^{c}_{\chi_{k}}\mathbf{\Sigma}^{c}_{\chi_{k}}\mathbf{V}^{c*}_{\chi_{k}}$ be known, and let the pair $(\mathbf{A}^{w_{c}}_{k},\mathbf{B}^{w_{c}}_{k})$ minimize the cost function in (11) with $\boldsymbol{\Upsilon}_{k}$ and $\boldsymbol{\chi}_{k}$ given in (44). Assume that at sampling time $t_{k+2}$, a new triplet of measurements $(\mathbf{x}_{k+1},\mathbf{y}_{k+1},\boldsymbol{\gamma}_{k+1})$ is used to incrementally compute the SVD of $\boldsymbol{\chi}^{c}_{k+1}$ and $\boldsymbol{\chi}_{k+1}$ using Proposition 2 and equation (3.2). Then, the pair $(\mathbf{A}^{w_{c}}_{k+1},\mathbf{B}^{w_{c}}_{k+1})$ which minimizes the cost function in (11) is

[TABLE]

If $\boldsymbol{\chi}^{c}_{k+1}$ can be well-approximated by its projection onto its leading $p$ singular vectors such that $\bar{\boldsymbol{\chi}}^{c}_{k+1}=\bar{\mathbf{U}}^{c}_{\chi_{k+1}}\bar{\mathbf{\Sigma}}^{c}_{\chi_{k+1}}\bar{\mathbf{V}}^{c*}_{\chi_{k+1}}$, where $\bar{\mathbf{U}}_{\chi_{k+1}}^{c}\in\mathbb{C}^{(n+l)\times p}$, $\bar{\mathbf{\Sigma}}_{\chi_{k+1}}^{c}\in\mathbb{C}^{p\times p}$, and $\bar{\mathbf{V}}_{\chi_{k+1}}^{c}\in\mathbb{C}^{w\times p}$, then equation (46) yields the following updates for $\bar{\mathbf{A}}^{w_{c}}_{k+1}\in\mathbb{R}^{n\times n}$ and $\bar{\mathbf{B}}^{w_{c}}_{k+1}\in\mathbb{R}^{n\times n}$:

[TABLE]

where $\bar{\mathbf{U}}^{c_{a}}_{\chi_{k+1}}\in\mathbb{C}^{n\times p}$ and $\bar{\mathbf{U}}^{c_{b}}_{\chi_{k+1}}\in\mathbb{C}^{l\times p}$ are defined such that $\bar{\mathbf{U}}^{c*}_{\chi_{k+1}}=[\bar{\mathbf{U}}^{c_{a}*}_{\chi_{k+1}}\bar{\mathbf{U}}^{c_{b}*}_{\chi_{k+1}}]$.

Assume that $\boldsymbol{\chi}_{k+1}$ can be projected onto its leading $r$ singular vectors such that $\bar{\boldsymbol{\chi}}_{k+1}=\bar{\mathbf{U}}_{\chi_{k+1}}\bar{\mathbf{\Sigma}}_{\chi_{k+1}}\bar{\mathbf{V}}^{*}_{\chi_{k+1}}$ with $r\leq p$. Then a reduced order model can be represented as:

[TABLE]

where the lower-dimensional approximation of $\tilde{\mathbf{A}}^{w_{c}}_{k+1}\in\mathbb{R}^{r\times r}$ and $\tilde{\mathbf{B}}^{w_{c}}_{k+1}\in\mathbb{R}^{r\times l}$ can obtained by projecting $\bar{\mathbf{A}}^{w_{c}}_{k+1}$ and $\bar{\mathbf{B}}^{w_{c}}_{k+1}$ onto a subspace spanned by the columns of $\bar{\mathbf{U}}_{\chi_{k+1}}$ as follows:

[TABLE]

For simplicity of exposition, in Section 6.1 and 6.2, we assumed that the SVD of the data matrix $\mathbf{X}_{k}~{}(\boldsymbol{\chi}_{k})$ and the augmented data matrix $\mathbf{X}^{c}_{k}~{}(\boldsymbol{\chi}^{c}_{k})$ is known. However, if the dimension of the input $\boldsymbol{\gamma}_{k}$ is small, or if the large dimension of the data matrix leads to storage concerns for two SVDs, it may be beneficial to store only the SVD of $\mathbf{X}_{k}~{}(\boldsymbol{\chi}_{k})$ and subsequently, use incremental SVD updates to first compute SVD of $\mathbf{X}_{k+1}~{}(\boldsymbol{\chi}_{k+1})$ and then compute the SVD of $\mathbf{X}^{c}_{k+1}~{}(\boldsymbol{\chi}^{c}_{k+1})$. The former can be accomplished using the incremental SVD updates described in Section 3, while the latter update is described in Appendix I.

6.3 Numerical Illustration for the Incremental DMDc Algorithm

In this section, we illustrate the efficiency of incremental DMDc algorithms. To this end, we generate a time-varying linear dynamical system as follows. We generate random linear discrete time invariant system matrices $\mathbf{A}\in\mathbb{R}^{n\times n}$ and $\mathbf{B}\in\mathbb{R}^{n\times l}$ using MATLAB function drss. These system matrices were used to generate the following time-varying linear system:

[TABLE]

where $\mathbf{A}_{k}=(1+\epsilon\sin{\omega k})\mathbf{A}$, and $\mathbf{B}_{k}=(1+\epsilon\sin{\omega k})\mathbf{B}$ and $\boldsymbol{\gamma}_{k}\in\mathbb{R}^{l}$ is selected as an i.i.d. sequence of standard Gaussian random vector. The system in (47) is used to run a simulation for $m\gg n$ time steps, to obtain two datasets $\mathbf{X},\mathbf{Y}\in\mathbb{R}^{n\times(m-1)}$, with $\omega=1$, $\epsilon=0.001$, $m=200$, $n=20$, and $l=2$. Initial models are estimated using DMD and DMDc algorithms using initial windows of $\mathbf{X}=\begin{bmatrix}\mathbf{x}_{1},\cdots,\mathbf{x}_{40}\end{bmatrix},\mathbf{Y}=\begin{bmatrix}\mathbf{y}_{1},\cdots,\mathbf{y}_{40}\end{bmatrix}$, and $\mathbf{\Gamma}=\begin{bmatrix}\boldsymbol{\gamma}_{1},\cdots,\boldsymbol{\gamma}_{40}\end{bmatrix}$. The weighted ($\rho=0.9$) and windowed ($w=40$ samples) incremental DMD and DMDc algorithms are applied at each iteration on the generated data $\mathbf{x}_{k},\mathbf{y}_{k},\boldsymbol{\gamma}_{k}$, for $k\in\{41,\cdots,200\}$, to update the initial DMD model and DMDc model. The two models are also used to predict $10$ future state vectors at each iteration. The Frobenius norm of the prediction error for incremental DMDc (blue line) and the incremental DMD (red line) are shown in Figure 12. The incremental DMDc models have higher prediction accuracy because of its ability to characterize the relationship between the states and the control input which is vital for any predictive model. The windowed incremental DMDc algorithm appears to have smaller prediction error than the weighted incremental DMDc algorithm.

7 Conclusion and Future Direction

In this paper, we developed algorithms for incremental computation of time-varying dynamic mode decomposition for autonomous and non-autonomous systems. In contrast to existing algorithms, these algorithms rely on incremental singular value decomposition to update singular values of the data matrix and allow for computation of a reduced order model at each time step. These algorithms are particularly useful for the cases in which the data matrix is singular and incremental matrix inversion based algorithms cannot be applied. We applied the proposed algorithms to an EEG dataset associated with error-related potentials and showed the efficacy of the algorithms in terms of predicting the future EEG signal. We also illustrated that the principal DMD modes obtained were consistent with the EEG activity seen in the brain.

There are several interesting future research directions for this work. First, the proposed approach can be adapted to apply the incremental updates on a nonlinear mapping of the measurements. This could be done in the spirit of extended dynamic mode decomposition proposed in [34]. Such extension may further improve the predictive power of the computed models. Second, it would be interesting to conduct human-in-the-loop experiments in which the EEG data is used in real-time to computed a DMD-based model of human performance and the control is designed to improve the performance.

Code

The MATLAB code used to generate the numerical results presented in this paper is available at

https://github.com/MSU-dcypherlab/Incremental-DMD-for-EEG-Data.git.

Appendix A Proof of Proposition 1

The incremented dataset $\mathbf{X}_{k+1}$ can be represented in term of $\mathbf{X}_{k}$ using the following additive updating formula

[TABLE]

where $\mathbf{z}_{k+1}=\begin{bmatrix}0&0&\cdots&1\end{bmatrix}^{\text{T}}\in\mathbb{R}^{k+1}$.

Let the SVD of $\mathbf{S}_{k}\triangleq\begin{bmatrix}\mathbf{\Sigma}_{x_{k}}&\rho^{-1}\mathbf{U}_{x_{k}}^{*}\mathbf{x}_{k+1}\end{bmatrix}=\mathbf{U}_{s_{k}}\mathbf{\Sigma}_{s_{k}}\mathbf{V}^{*}_{s_{k}}$ with $\mathbf{U}_{s_{k}}\in\mathbb{C}^{n\times n}$, $\mathbf{\Sigma}_{s_{k}}\in\mathbb{C}^{n\times n}$, and $\mathbf{V}_{s_{k}}\in\mathbb{C}^{(n+1)\times n}$. Then, $\mathbf{V}_{s_{k}}$ can be partitioned as $\mathbf{V}_{s_{k}}=\begin{bmatrix}\mathbf{V}_{s_{k,1}}\\ \mathbf{v}_{s_{k,2}}\end{bmatrix}$ such that $\mathbf{v}_{s_{k,2}}\in\mathbb{C}^{1\times n}$ is orthogonal to $\mathbf{V}_{s_{k,1}}\in\mathbb{C}^{n\times n}$, and $\mathbf{X}_{k+1}$ in (22) can be written as:

[TABLE]

Thus, the SVD of $\mathbf{X}_{k+1}=\mathbf{U}_{x_{k+1}}\mathbf{\Sigma}_{x_{k+1}}\mathbf{V}^{*}_{x_{k+1}}$ is given by

[TABLE]

Appendix B Proof of Proposition 2

Consider the updated dataset

[TABLE]

where $\boldsymbol{\acute{\tilde{\chi}}}_{k}\triangleq\begin{bmatrix}0&\mathbf{x}_{k-w+2}&\cdots&\mathbf{x}_{k}\end{bmatrix}$ is the $\boldsymbol{\acute{\chi}}_{k}$ dataset padded with $n$-dimensional zero column, and $\mathbf{z}_{1}=\begin{bmatrix}1&0&\cdots&0\end{bmatrix}^{\text{T}}\in\mathbb{R}^{w}$.

Let the SVD of $\mathbf{\acute{S}}_{k}\triangleq\mathbf{\Sigma}_{\chi_{k}}-\mathbf{U}_{\chi_{k}}^{*}\mathbf{x}_{k-w+1}\mathbf{z}_{1}^{\text{T}}\mathbf{V}_{\chi_{k}}$ be

[TABLE]

where $\mathbf{U}_{\acute{s}_{k}}\in\mathbb{C}^{q\times q}$, ${\mathbf{\Sigma}}_{\acute{s}_{k}}\in\mathbb{C}^{q\times q}$, $\mathbf{V}_{\acute{s}_{k}}\in\mathbb{C}^{w\times q}$ and $\tilde{\mathbf{V}}_{\acute{s}_{k}}\in\mathbb{C}^{w\times(w-q)}$. Let $\mathbf{V}_{\chi_{k}}$ be partitioned into $\mathbf{V}_{\chi_{k}}=\begin{bmatrix}\mathbf{v}_{\chi_{k,1}}\\ \mathbf{V}_{\chi_{k,2}}\end{bmatrix}$, where $\mathbf{v}_{\chi_{k,1}}\in\mathbb{C}^{1\times w}$ and $\mathbf{V}_{\chi_{k,2}}\in\mathbb{C}^{(w-1)\times w}$.

Substituting the SVD from equation (50) into equation(49), the dataset $\mathbf{\acute{{\tilde{X}}}}_{k}$ can be expressed as:

[TABLE]

Since, $\boldsymbol{\acute{\tilde{\chi}}}_{k}=\begin{bmatrix}\mathbf{0}&\boldsymbol{\acute{{\chi}}}_{k}\end{bmatrix}$, it follows that $\mathbf{U}_{\chi_{k}}\mathbf{U}_{\acute{s}_{k}}{\mathbf{\Sigma}}_{\acute{s}_{k}}\mathbf{V}_{\acute{s}_{k}}^{*}\mathbf{v}_{\chi_{k,1}}^{*}=\mathbf{0}$ and $\boldsymbol{\acute{{\chi}}}_{k}=\mathbf{U}_{\chi_{k}}\mathbf{U}_{\acute{s}_{k}}{\mathbf{\Sigma}}_{\acute{s}_{k}}\mathbf{V}_{\acute{s}_{k}}^{*}\mathbf{V}_{\chi_{k,2}}^{*}$. Thus, the SVD of $\boldsymbol{\acute{\chi}}_{k}=\mathbf{U}_{\acute{\chi}_{k}}{\mathbf{\Sigma}}_{\acute{\chi}_{k}}{\mathbf{V}}_{\acute{\chi}_{k}}^{*}$ is

[TABLE]

A noteworthy property of the above decomposition that we will use in the later developments is that $\mathbf{v}_{\chi_{k,1}}\mathbf{V}_{\acute{s}_{k}}=\mathbf{0}$. To establish this property, we note that since $\mathbf{U}_{\chi_{k}}\mathbf{\acute{S}}_{k}\mathbf{v}_{\chi_{k,1}}^{*}=\mathbf{0}$, $\mathbf{v}^{*}_{\chi_{k,1}}\in\text{Ker}(\mathbf{U}_{\chi_{k}}\mathbf{\acute{S}}_{k})$. Furthermore, since $\mathbf{U}_{\chi_{k}}$ is full rank, then $\mathbf{v}^{*}_{\chi_{k,1}}\in\text{Ker}(\mathbf{\acute{S}}_{k})=\text{span}(\tilde{\mathbf{V}}_{\acute{s}_{k}})$, and consequently

[TABLE]

Appendix C Proof of Theorem 1

Let the SVD of $\mathbf{X}_{k+1}=\mathbf{U}_{x_{k+1}}{\mathbf{\Sigma}_{x_{k+1}}}\mathbf{V}^{*}_{x_{k+1}}$ be calculated by the incremental SVD update (23). Then, the time-varying DMD operator (15) at sampling time $t_{k+1}$ is:

[TABLE]

Substituting the incrementally computed $\mathbf{V}_{k+1}$ from (23) into the DMD operator (53) yields:

[TABLE]

where $I_{n}$ is the identity matrix of order $n$, and $\mathbf{A}^{\rho}_{k}$ has been added and subtracted in (54). The term $\left(\mathbf{U}_{x_{k+1}}\mathbf{\Sigma}_{x_{k+1}}-\rho\mathbf{U}_{x_{k}}\mathbf{\Sigma}_{x_{k}}\mathbf{V}_{s_{k,1}}\right)$ can be rewritten as

[TABLE]

where (56) is obtained by substituting incremental update for $\mathbf{V}_{x_{k+1}}$ from (23), and (57) is obtained by using the fact that the rows in $\mathbf{V}_{s_{k}}=\begin{bmatrix}\mathbf{V}_{s_{k,1}}\\ \mathbf{v}_{s_{k,2}}\end{bmatrix}$ are orthogonal to each other. Substituting equation (58) into the equation (55):

[TABLE]

Appendix D Proof of Theorem 2

Let the SVD of $\boldsymbol{\chi}_{k+1}=\mathbf{U}_{\chi_{k+1}}\mathbf{\Sigma}_{\chi_{k+1}}\mathbf{V}^{*}_{\chi_{k+1}}$ be calculated by the windowed SVD update (3.2). Then, the time-varying DMD operator (15) at sampling time $t_{k+1}$ is

[TABLE]

Substituting the incrementally computed $\mathbf{V}_{k+1}$ from (3.2) into the DMD operator (59) yields:

[TABLE]

where equation (D) holds because $\mathbf{v}_{\chi_{k,1}}\mathbf{V}_{\acute{s}_{k}}=\mathbf{0}$ from equation (52).

After following the same steps presented in the proof of Theorem 1, equation (61) becomes:

[TABLE]

The term $\left(\mathbf{U}_{\chi_{k+1}}\mathbf{\Sigma}_{\chi_{k+1}}-\mathbf{U}_{\chi_{k}}\mathbf{\Sigma}_{\chi_{k}}\mathbf{V}_{\acute{s}_{k}}\mathbf{V}_{\hat{s}_{k,1}}\right)$ can be rewritten as

[TABLE]

where equation (63) is obtained after substituting the windowed incremental update formula from (3.2), and (64) is obtained after substituting $\mathbf{V}_{\acute{s}_{k}}\mathbf{V}_{\acute{s}_{k}}^{*}=\mathbf{V}_{\hat{s}_{k,1}}\mathbf{V}_{\hat{s}_{k,1}}^{*}=\mathbf{I}$ and $\mathbf{V}_{\hat{s}_{k,1}}\mathbf{v}_{\hat{s}_{k,2}}^{*}=\mathbf{0}$. The substitution of (65) in (62) yields:

[TABLE]

Appendix E EEG Prediction Error for Incremental DMD and Online DMD Algorithms

The figures below show the normalized RMS for the predicted EEG signal at channel FCz during correct and erroneous events using incremental DMD with $\sigma_{thr}=\{0.01,0.001\}$, $\rho=\{0.1,0.2,0.4,0.8\}$ and $w=512$, and online DMD with $\rho=\{0.1,0.2,0.4,0.8\}$ and $w=512$.

Appendix F ERP Prediction Error for Incremental DMD Algorithm with Different Initial Window Sizes

Figure 14 show the normalized RMS for the prediction error of the ERP signal during the correct and erroneous events using weighted incremental DMD with initial window size of $w_{0}=128$ samples, and Figure 15 shows the normalized RMS for the same error using windowed incremental DMD with $w=\{64,128,256\}$.

Appendix G Proof of Theorem 3

Consider the SVD of $\mathbf{X}^{c}_{k}=\mathbf{U}^{c}_{x_{k}}\mathbf{\Sigma}^{c}_{x_{k}}\mathbf{V}^{c*}_{x_{k}}$, where $\mathbf{U}_{x_{k}}^{c}\in\mathbb{C}^{(n+l)\times(n+l)}$, $\mathbf{\Sigma}_{x_{k}}^{c}\in\mathbb{C}^{(n+l)\times(n+l)}$, and $\mathbf{V}_{x_{k}}^{c}\in\mathbb{C}^{k\times(n+l)}$, and the SVD of $\mathbf{S}^{c}_{k}\triangleq\begin{bmatrix}\mathbf{\Sigma}^{c}_{x_{k}}&\rho^{-1}\mathbf{U}_{x_{k}}^{c*}\mathbf{x}^{c}_{k+1}\end{bmatrix}=\mathbf{U}^{c}_{s_{k}}\mathbf{\Sigma}_{s_{k}}^{c}\begin{bmatrix}\mathbf{V}^{c*}_{s_{k,1}}&\mathbf{v}^{c*}_{s_{k,2}}\end{bmatrix}$. Then, the weighted incremental SVD presented in Proposition 1 can be applied to obtain the SVD of $\mathbf{X}^{c}_{k+1}=\mathbf{U}^{c}_{x_{k+1}}\mathbf{\Sigma}^{c}_{x_{k+1}}\mathbf{V}^{c*}_{x_{k+1}}$ as follows

[TABLE]

Given the SVD of $\mathbf{X}_{k}=\mathbf{U}_{x_{k}}\mathbf{\Sigma}_{x_{k}}\mathbf{V}^{*}_{x_{k}}$, the updated SVD for $\mathbf{X}_{x_{k+1}}$ is given by (23). Following the same steps presented in the proof of Theorem 1 with the SVD updates in (66), we can obtain the following weighted incremental DMDc update

[TABLE]

where $\mathbf{G}^{\rho_{c}}_{k}=\begin{bmatrix}\mathbf{A}^{\rho_{c}}_{k}&\mathbf{B}^{\rho_{c}}_{k}\end{bmatrix}$. Let ${\mathbf{U}}_{x_{k+1}}^{c*}$ be partitioned such that ${\mathbf{U}}_{x_{k+1}}^{c*}=[{\mathbf{U}}^{c_{a}*}_{x_{k+1}}{\mathbf{U}}^{c_{b}*}_{x_{k+1}}]$, where ${\mathbf{U}}^{c_{a}}_{x_{k+1}}\in\mathbb{C}^{n\times n}$ and ${\mathbf{U}}^{c_{b}}_{x_{k+1}}\in\mathbb{C}^{l\times n}$. Then, the update in (67) yields the following updates for ${\mathbf{A}}^{\rho_{c}}_{k+1}\in\mathbb{R}^{n\times n}$ and ${\mathbf{B}}^{\rho_{c}}_{k+1}\in\mathbb{R}^{n\times l}$

[TABLE]

Appendix H Proof of Theorem 4

Let $q_{c}$ be defined by

[TABLE]

Consider the SVD of $\boldsymbol{\chi}_{k}^{c}=\mathbf{U}^{c}_{\chi_{k}}\mathbf{\Sigma}^{c}_{\chi_{k}}\mathbf{V}^{c*}_{\chi_{k}}$, where $\mathbf{U}_{\chi_{k}}^{c}\in\mathbb{C}^{(n+l)\times q_{c}}$, $\mathbf{\Sigma}_{\chi_{k}}^{c}\in\mathbb{C}^{q_{c}\times w}$, and $\mathbf{V}_{\chi_{k}}^{c}\in\mathbb{C}^{w\times w}$. Define the matrices $\mathbf{\acute{S}}^{c}_{k}$ and $\mathbf{\hat{S}}^{c}_{k}$ with their associated SVDs by

[TABLE]

Then, the windowed incremental SVD presented in Proposition 2 and equation (3.2) can be applied to obtain the SVD of $\boldsymbol{\chi}^{c}_{k+1}=\mathbf{U}_{\chi_{k+1}}^{c}\mathbf{\Sigma}_{\chi_{k+1}}^{c}\mathbf{V}^{c*}_{\chi_{k+1}}$ as follows

[TABLE]

Similarly, given the SVD of $\boldsymbol{\chi}_{k}=\mathbf{U}_{\chi_{k}}\mathbf{\Sigma}_{\chi_{k}}\mathbf{V}^{*}_{\chi_{k}}$, the SVD of $\boldsymbol{\chi}_{k+1}=\mathbf{U}_{\chi_{k+1}}\mathbf{\Sigma}_{\chi_{k+1}}\mathbf{V}^{*}_{\chi_{k+1}}$ can be obtained using the updates in (3.2).

The same steps presented in the proof of Theorem 2 can be followed with the SVD update in (70) to obtain the following update

[TABLE]

Let ${\mathbf{U}}_{\chi_{k+1}}^{c*}$ be partitioned such that ${\mathbf{U}}_{\chi_{k+1}}^{c*}=[{\mathbf{U}}^{c_{a}*}_{\chi_{k+1}}{\mathbf{U}}^{c_{b}*}_{\chi_{k+1}}]$, where ${\mathbf{U}}^{c_{a}}_{\chi_{k+1}}\in\mathbb{C}^{n\times n}$ and ${\mathbf{U}}^{c_{b}}_{\chi_{k+1}}\in\mathbb{C}^{l\times n}$. Then, the updates in (71) yields the following updates for ${\mathbf{A}}^{w_{c}}_{k+1}\in\mathbb{R}^{n\times n}$ and ${\mathbf{B}}^{w_{c}}_{k+1}\in\mathbb{R}^{n\times l}$

[TABLE]

Appendix I Incremental Singular Value Decomposition for Data Matrix Augmented with Exogenous Input

In this section, we present an incremental update to compute the SVD of the data matrix augmented with the exogeneous input. This update can be beneficial in case of low-dimensional control input or limited storage capacity. Assume that at $t_{k}$, the state measurements are arranged in a matrix denoted by $\mathbf{X}_{k}\in\mathbb{R}^{n\times q_{m}}$, and the exogeneous input measurements are arranged in a matrix denoted by $\mathbf{\Gamma}_{k}\in\mathbb{R}^{l\times q_{m}}$, where $(q_{m}=k)$ for weighted incremental SVD and $(q_{m}=w)$ for windowed incremental SVD. Both $\mathbf{X}_{k}$ and $\mathbf{\Gamma}_{k}$ are arranged in an augmented matrix denoted by $\mathbf{X}^{c}_{k}\in\mathbb{R}^{(n+l)\times q_{m}}$ as follows:

[TABLE]

where for weighted incremental SVD we have $\mathbf{\Gamma}_{k}=\begin{bmatrix}\rho^{k-1}\boldsymbol{\gamma}_{1}&\cdots&\boldsymbol{\gamma}_{k}\end{bmatrix}$ and $\mathbf{X}_{k}=\begin{bmatrix}\rho^{k-1}\mathbf{x}_{1}&\cdots&\mathbf{x}_{k}\end{bmatrix}$, and for windowed incremental SVD we have $\mathbf{\Gamma}_{k}=\begin{bmatrix}\boldsymbol{\gamma}_{k-w+1}&\cdots&\boldsymbol{\gamma}_{k}\end{bmatrix}$ and $\mathbf{X}_{k}=\begin{bmatrix}\mathbf{x}_{k-w+1}&\cdots&\mathbf{x}_{k}\end{bmatrix}$. Let $q_{n}\triangleq n$ for weighted incremental SVD and $q_{n}\triangleq q$ for windowed incremental SVD, where $q$ is given in (25). Then, the SVD of $\mathbf{X}^{c}_{k}$ can be calculated by updating the SVD of $\mathbf{X}_{k}$ using the following proposition.

Proposition 3**.**

Let at sampling time $t_{k}$, the SVD of $\mathbf{X}_{k}=\mathbf{U}_{k}\mathbf{\Sigma}_{k}\mathbf{V}_{k}^{*}$ be known, where $\mathbf{U}_{k}\in\mathbb{C}^{n\times q_{n}}$, $\mathbf{\Sigma}_{k}\in\mathbb{C}^{q_{n}\times q_{m}}$, and ${\mathbf{V}_{k}}\in\mathbb{C}^{q_{m}\times q_{m}}$. Then, the SVD of the dataset in (72), defined by $\mathbf{X}^{c}_{k}=\mathbf{U}^{c}_{k}\mathbf{\Sigma}^{c}_{k}\mathbf{V}^{c*}_{k}$, is

[TABLE]

where $\mathbf{U}^{c}_{R}=\begin{bmatrix}\mathbf{U}^{c}_{R_{1}}\\ \mathbf{U}^{c}_{R_{2}}\end{bmatrix}$, $\mathbf{\Sigma}^{c}_{R}$, and ${\mathbf{V}}^{c}_{R}$ are given by the following SVD:

[TABLE]

Proof.

The new dataset matrix $\mathbf{X}^{c}_{k}$ can be represented in term of $\mathbf{X}_{k}$ using the following additive update formula:

[TABLE]

where $\mathbf{Z}_{n+1}=\begin{bmatrix}\mathbf{0}&\mathbf{0}&\dots&\mathbf{I}_{l}\end{bmatrix}^{\text{T}}\in\mathbb{R}^{(n+l)\times l}$, where $\mathbf{I}_{l}$ is the identity matrix of order $l$. Then, the new dataset in (74) can be written as:

[TABLE]

Define $\mathbf{R}^{c}\triangleq\begin{bmatrix}\mathbf{\Sigma}_{k}\\ \left({\mathbf{V}^{*}_{k}}\mathbf{\Gamma}_{k}\right)^{*}\end{bmatrix}$, which has the SVD of $\mathbf{R}^{c}=\mathbf{U}^{c}_{R}\mathbf{\Sigma}^{c}_{R}{\mathbf{V}}^{{c_{l}}*}_{{R}}$, where $\mathbf{U}^{c}_{R}=\begin{bmatrix}\mathbf{U}^{c}_{R_{1}}\\ \mathbf{U}^{c}_{R_{2}}\end{bmatrix}\in\mathbb{C}^{(n+l)\times q_{n}}$, $\mathbf{U}^{c}_{R_{1}}\in\mathbb{C}^{n\times q_{n}}$, $\mathbf{U}^{c}_{R_{2}}\in\mathbb{C}^{l\times q_{n}}$, $\mathbf{\Sigma}^{c}_{R}\in\mathbb{C}^{q_{n}\times q_{m}}$, and $\mathbf{V}^{c}_{R}\in\mathbb{C}^{q_{m}\times q_{m}}$. Then, after substituting the SVD of matrix $\mathbf{R}^{c}$ in (75), the updated dataset becomes:

[TABLE]

Finally the SVD factors for $\mathbf{X}^{c}_{k}$ is defined by:

[TABLE]

∎

For $l=1$, the matrix $\mathbf{R}$ has a broken-arrow structure and its SVD can be computed efficiently [2]. For $l>1$, the matrix is still sparse but does not carry the broken-arrow structure. In this case, either the above procedure can be applied recursively with one input added at a time or problem specific sparse SVD solvers can be employed.