Exploring Dynamical Systems With DMD: Part 1

Posted by Kristian M. Eschenburg on May 22, 2018 · 10 mins read

In the next two posts, I want to talk briefly about an algorithm called Dynamic Mode Decomposition (DMD). DMD is a spatiotemporal modal decomposition technique that can be used to identify spatial patterns in a signal (modes), along with the time course of these spatial patterns (dynamics). As such, the algorithm assumes that the input data has a both a spatial and a temporal component. We are interested in modeling how the system evolves over time.

If you’d like to find more information about DMD, (Schmid, 2010; Tu, Rowley, Luchtenburg, Brunton, & Kutz, 2014) are good references. Likewise, if you’d like to follow along with the code for the following analysis, see my repo. For a more in-depth analysis that applies DMD to brain activity in the resting brain, see this pre-print of a paper my colleagues and I wrote (Kunert-Graf et al., 2018), along with the code used for our analysis.

The DMD Algorithm

Let’s assume that you’ve taken measurements from specific points in space for time points, where for now we assume that . For now, we’ll assume that the sampling frequency, , is stable across the entire experiment. We define our entire data matrix as

We are interested in solving for the matrix, , such that

Given our full data matrix , we can define two matrices and such that

so that we can write

If $n$ is small, this is relatively easy to compute – however, if is large, as is the case when modeling temporal dynamics in resting-state MRI, it would be computationally inefficient to compute A directly. To alleviate this, we can make use of the Singular Value Decomposition (SVD) of our predictor matrix . We define the SVD of as

as well as the Moore-Penrose psuedo-inverse of as

such that we can write

Additionally, if we assume that , then we can use the truncated SVD such that

and

As it stands now, we still compute an matrix. However, because we have a potentially low-rank system, we can apply a Similarity Transformation to in order to reduce its dimensionality, without changing its spectrum. Using our spatial singular vectors $U$, we define

where . If we consider the above SVD, we see that $U$ is the matrix of left singular vectors, an orthogonal basis that spans , which is an r-dimensional subspace of . Thus, the similarity transform represents a mapping . We now have a reduced-dimensional representation of our linear operator, from which we can compute the spatial modes and dynamic behavior of each mode. First, however, because of the notion of variance captured by the singular values of our original predictor matrix, we weight by the singular values as

such that our computed spatial modes have been weighted by the amount they contribute to our measured signal. We can now compute the eigendecomposition of as

where the eigenvectors are the reduced-dimension representations of our spatial modes, and the eigenvalues capture the dynamic behavior of our spatial modes. Because our original data matrix had spatial dimension and our eigenvectors have dimension , we need to up-project our eigenvectors to compute the final spatial modes, via

From the SVD of our prediction matrix , the matrix is the matrix of right singular vectors, an orthogonal basis spanning the space of (i.e. $r$ basis vectors spanning the space of the measured time courses). Thus, we see that represents a linear combination of the temporal basis vectors (a mapping from ) for each eigenvector of , weighted by the corresponding singular value (that acts to normalize the spatial mode amplitudes). Finally, we see that computes how much of each temporal basis vector is present in the measured time course at each point in space.

Because we are modeling a dynamical system, we can compute the continuous time dynamics of our system using our spatial modes and eigenvalues as

where is a growth-decay constant and is the frequency of oscillation of the spatial mode . We can compute these two constants as

So, we can see that DMD linearizes our measured time series, by fitting what can be analogized to a “global” regression. That is, instead of computing how a single time point predicts the next time point, which could readily be solved using the simple Normal equations, DMD computes how a matrix of time points predicts another matrix of time points that is shifted one unit of time into the future. To this extent, DMD minimizes the Frobenius norm of

However, rather than explicitly computing the matrix , DMD computes the eigenvectors and eigenvalues of , by utilizing the Singular Value Decomposition, along with a Similarity Transformation, in order to generate a reduced-dimensional representation of .

This spectral decomposition of our linear operator is of particular importance, because it sheds light on the fact the DMD models the temporal dynamics of our system using a Fourier basis. Each spatial mode is represented by a particular Fourier frequency along and growth-decay constant that determines the future behavior of our spatial mode. Additionally, the Fourier basis also determines what sorts of time series can be modeled using DMD – time series that are expected to have sinusoidal behavior will be more reliably modeled using DMD, whereas signals that show abrupt spike patterns might be more difficult to model.

  1. Schmid, P. (2010). Dynamic Mode Decomposition of Numerical and Experimental Data. Journal of Fluid Mechanics, 656, 5–28.
  2. Tu, J., Rowley, C., Luchtenburg, D., Brunton, S., & Kutz, N. (2014). On Dynamic Mode Decomposition: Theory and Applications. Journal of Computational Dynamics.
  3. Kunert-Graf, J., Eschenburg, K., Galas, D. J., Kutz, N., Rane, S. D., & Brunton, B. W. (2018). Extracting Time-Resolved Resting State Networks Using Dynamic Mode Decomposition. BioRxiv.