Dynamic Mode Decomposition
The Dynamic Mode Decomposition (DMD) is a relatively recent mathematical innovation that, among other things, allows us to solve or approximate dynamical systems in terms of coherent structures that grow, decay, and/ or oscillate in time. We refer to the coherent structures as DMD modes. Each DMD mode has corresponding time dynamics defined in terms of a single eigenvalue.
In other words, the DMD converts a dynamical system into a superposition of modes whose dynamics are governed by eigenvalues.
Amazingly, although the mathematical procedure for identifying the DMD modes and eigenvalues is purely linear, the system itself can be nonlinear! I won’t go into it here but there are sound theoretical underpinnings for the claim that a nonlinear system can be described by a set of mode and eigenvalue pairs. Read up on the Koopman operator and the DMD’s connection to it for more information ^{1}^{2}^{3}.
Not only is the DMD a useful diagnostic tool for analyzing the inner workings of a system, but it can also be used to predict the future state of the system. All that is needed are the modes and the eigenvalues. With very little effort, the modes and eigenvalues can be combined to produce a function that approximates the system state at any moment in time.
In this post, I start by briefly introducing the mathematical theory of the DMD. Then I go over a number of examples with annotated Python code. The examples are just toy problems, but they succeed in illustrating how to use the DMD in several different situations, including:

Given: a 1dimensional scalar function evolving in time.

Given: a set of trajectories in 3 dimensions produced by an unknown vector field.

Given: a set of 3dimensional vectors sampled from an unknown vector field.

Given: a 2dimensional scalar function in spherical coordinates evolving in time.
It is important to fully understand the limitations of any modeling strategy. Thus, I also talk a bit about when and how the DMD can fail. Finally, I conclude the post with a mention of notable extensions to the DMD and a summary.
Formal Definition
Consider a set of
in which each
If we define an operator
where
The definition above, from Tu et al. ^{2}, is known as the exact DMD. It is currently the most general definition and can be applied to any dataset that satisfies the given requirements. In this post, we are mostly interested in the cases where
Clearly,
There are a number of theorems that go along with this definition of the DMD ^{2}. One of the more useful theorems states that
DMD Algorithm
At first glance, the task of finding the eigendecomposition of pinv
and eig
methods from Numpy or MATLAB will do the trick. The problem comes when
Fortunately, the problem can be broken down into smaller pieces with the help of the exact DMD algorithm.

Compute the SVD of
$X$ , optionally performing a lowrank truncation at the same time:$$X=U\Sigma V^*$$ 
Compute
$\tilde A$ , the projection of the full matrix$A$ onto$U$ :$$\tilde A=U^* AU=U^*YV\Sigma^{1}$$ 
Compute the eigenvalues
$\lambda_i$ and eigenvectors$w_i$ of$\tilde A$ :$$\tilde AW=W\Lambda$$ 
Reconstruct the eigendecomposition of
$A$ from$W$ and$\Lambda$ . The eigenvalues for$A$ are equivalent to those of$\tilde A$ . The eigenvectors of$A$ are given by the columns of$\Phi$ :$$A\Phi=\Phi\Lambda,\quad \Phi=YV\Sigma^{1}W$$
A more indepth explanations of the algorithm’s derivation can be found in the references ^{1}^{2}. It also might be of theoretical interest to note that
Now let’s go through the algorithm stepbystep in Python. Start by installing and importing all necessary packages.
import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from numpy import dot, multiply, diag, power from numpy import pi, exp, sin, cos, cosh, tanh, real, imag from numpy.linalg import inv, eig, pinv from scipy.linalg import svd, svdvals from scipy.integrate import odeint, ode, complex_ode from warnings import warn
Let’s generate some play data. Keep in mind that in practice, one doesn’t necessarily know the governing equations for a data. Here we’re just inventing some equations to create a dataset. Forget they ever existed once the data has been generated.
# define time and space domains x = np.linspace(10, 10, 100) t = np.linspace(0, 6*pi, 80) dt = t[2]  t[1] Xm,Tm = np.meshgrid(x, t) # create three spatiotemporal patterns f1 = multiply(200.2*power(Xm, 2), exp((2.3j)*Tm)) f2 = multiply(Xm, exp(0.6j*Tm)) f3 = multiply(5*multiply(1/cosh(Xm/2), tanh(Xm/2)), 2*exp((0.1+2.8j)*Tm)) # combine signals and make data matrix D = (f1 + f2 + f3).T # create DMD inputoutput matrices X = D[:,:1] Y = D[:,1:]
Now compute the SVD of
# SVD of input matrix U2,Sig2,Vh2 = svd(X, False)
Given the singular values shown above, we conclude that the data has three modes of any significant interest. Therefore, we truncate the SVD to only include those modes. Then we build
# rank3 truncation r = 3 U = U2[:,:r] Sig = diag(Sig2)[:r,:r] V = Vh2.conj().T[:,:r] # build A tilde Atil = dot(dot(dot(U.conj().T, Y), V), inv(Sig)) mu,W = eig(Atil)
Each eigenvalue in
Now build the exact DMD modes.
# build DMD modes Phi = dot(dot(dot(Y, V), inv(Sig)), W)
The columns of
This is where the DMD algorithm technically ends. Equipped with the eigendecomposition of
# compute time evolution b = dot(pinv(Phi), X[:,0]) Psi = np.zeros([r, len(t)], dtype='complex') for i,_t in enumerate(t): Psi[:,i] = multiply(power(mu, _t/dt), b)
The three plots above are the time dynamics of the three DMD modes. Notice how all three are oscillating. Furthermore, the second mode appears to grow exponentially, which is confirmed by the eigenvalue plot.
If you wish to construct an approximation to the original data matrix, simply multiply
# compute DMD reconstruction D2 = dot(Phi, Psi) np.allclose(D, D2) # True
Dynamical Systems
I only consider two interpretations of the equation
In this case, the operator
where
In this way, each pair of columns vectors from
Using the DMD, we find the eigendecomposition of
The corresponding function in continuous time
What is truly amazing is that we have just defined explicit functions in time using nothing but data! This is a good example of equationfree modeling.
The second interpretation of
In this case, the operator
After computing the DMD, the function in time is very similar to before (i.e., of the difference equation). If
Examples
Here are a few of examples of how to use the DMD with different types of experimental data. For convenience, I condense the DMD code into a single method and define some helper methods to check linear consistency and confirm my solutions.
def nullspace(A, atol=1e13, rtol=0): # from http://scipycookbook.readthedocs.io/items/RankNullspace.html A = np.atleast_2d(A) u, s, vh = svd(A) tol = max(atol, rtol * s[0]) nnz = (s >= tol).sum() ns = vh[nnz:].conj().T return ns def check_linear_consistency(X, Y, show_warning=True): # tests linear consistency of two matrices (i.e., whenever Xc=0, then Yc=0) A = dot(Y, nullspace(X)) total = A.shape[1] z = np.zeros([total, 1]) fails = 0 for i in range(total): if not np.allclose(z, A[:,i]): fails += 1 if fails > 0 and show_warning: warn('linear consistency check failed {} out of {}'.format(fails, total)) return fails, total def dmd(X, Y, truncate=None): U2,Sig2,Vh2 = svd(X, False) # SVD of input matrix r = len(Sig2) if truncate is None else truncate # rank truncation U = U2[:,:r] Sig = diag(Sig2)[:r,:r] V = Vh2.conj().T[:,:r] Atil = dot(dot(dot(U.conj().T, Y), V), inv(Sig)) # build A tilde mu,W = eig(Atil) Phi = dot(dot(dot(Y, V), inv(Sig)), W) # build DMD modes return mu, Phi def check_dmd_result(X, Y, mu, Phi, show_warning=True): b = np.allclose(Y, dot(dot(dot(Phi, diag(mu)), pinv(Phi)), X)) if not b and show_warning: warn('dmd result does not satisfy Y=AX')
3D Vector Field: Difference Equation
In this first example, we are given a set of trajectories in 3 dimensions produced by an vector field. I used the following system of differential equations to produce the trajectories.
Watch the video below for an animation of the dynamic behavior.
You can generate a dataset with the following code. Basically, it randomly chooses a set of initial points and simulates each of them forward in time. Then it concatenates all of the individual trajectories together to produce
# convenience method for integrating def integrate(rhs, tspan, y0): # empty array to hold solution y0 = np.asarray(y0) # force to numpy array Y = np.empty((len(y0), tspan.size), dtype=y0.dtype) Y[:, 0] = y0 # autodetect complex systems _ode = complex_ode if np.iscomplexobj(y0) else ode # create explicit RungeKutta integrator of order (4)5 r = _ode(rhs).set_integrator('dopri5') r.set_initial_value(Y[:, 0], tspan[0]) # run the integration for i, t in enumerate(tspan): if not r.successful(): break if i == 0: continue # skip the initial position r.integrate(t) Y[:, i] = r.y # return solution return Y # the righthand side of our ODE def saddle_focus(t, y, rho = 1, omega = 1, gamma = 1): return dot(np.array([ [rho, omega, 0], [omega, rho, 0], [0, 0, gamma] ]), y) # generate data dt = 0.01 tspan = np.arange(0, 50, dt) X = np.zeros([3, 0]) Y = np.zeros([3, 0]) z_cutoff = 15 # to truncate the trajectory if it gets too large in z for i in range(20): theta0 = 2*pi*np.random.rand() x0 = 3*cos(theta0) y0 = 3*sin(theta0) z0 = np.random.randn() D = integrate(saddle_focus, tspan, [x0, y0, z0]) D = D[:,np.abs(D[2,:]) < z_cutoff] X = np.concatenate([X, D[:,:1]], axis=1) Y = np.concatenate([Y, D[:,1:]], axis=1) check_linear_consistency(X, Y); # raise a warning if not linearly consistent
Now compute the DMD, just like before.
mu,Phi = dmd(X, Y, 3) check_dmd_result(X, Y, mu, Phi) # raise a warning if Y=AX not satisfied
The DMD computation was successful. Given an initial condition
# simulate an initial condition forward in time b = dot(pinv(Phi), np.array([3.45, 0.2, 0.98]).T) Psi = np.zeros([3, len(tspan)], dtype='complex') for i,_t in enumerate(tspan): Psi[:,i] = multiply(power(mu, _t/dt), b) Dtil = real(dot(Phi, Psi))
Of course, instead of using the DMD, it probably would have been a lot easier just to compute
3D Vector Field: Differential Equation
This example is very, very similar to the previous example. We shall use the same “mystery” vector field to create the dataset. The only difference is that the following code directly samples the output of the system of differential equations instead of simulating trajectories. In other words, we generate the data by randomly generating a set of position vectors,
# randomly generate a set of position vectors X = np.random.rand(3, 5000) * 10  5 # compute the corresponding velocity vectors Y = saddle_focus(0, X) # compute the DMD for the "unknown" data mu,Phi = dmd(X, Y, 3) check_dmd_result(X, Y, mu, Phi) # raise a warning if Y=AX not satisfied
The decomposition is successful. As extra confirmation, when we examine the eigenvalues, we notice that they are the same eigenvalues as the original system that produced our dataset. There are two oscillatingdecaying modes (complex conjugates) and one exponentially growing mode. Watch the video again to witness these exact time dynamics.
Spherical coordinates
This example demonstrates that the spatial dimension of
The mystery dataset is that of a
# define time and space domains q = np.linspace(0, 2 * pi, 40) # theta p = np.linspace(0, pi, 40) # phi t = np.linspace(0, 80, 800) # time dt = t[1]  t[0] Qm,Pm = np.meshgrid(q, p) # each column of the dataset is a flattened thetaphi matrix R = np.zeros([len(q)*len(p), len(t)], dtype='complex') # must be complex for i in range(len(t)): R[:,i] = ( 4 # a stationary mode + 0.4 * cos(10 * Qm) * cos(10 * Pm) * exp(t[i] * 3.2j) + 0.1 * cos(16 * Qm) * cos(22 * Pm) * exp(t[i] * 5.0j) + 0.6 * exp(t[i] * 0.7j) + 0.7 * cos(5 * Qm) * exp(t[i] * 1.5j)  0.3 * cos(9 * Pm) * exp(t[i] * 0.3j) ).flatten() # subtract off the mean (approx. 4) to prevent corruption of modes R = R.mean() # extract X,Y matrices X = R[:,:1] Y = R[:,1:] check_linear_consistency(X, Y)
It is important to note that during the spheroid investigation, I found that the existence of the stationary, nonzero mode (i.e., the constant radius of 4) caused the DMD to fail. One would think that the constant mode would be the easiest mode to identify, but it was not. Instead, the DMD tried to account for the influence of constant mode by injecting it into the other modes, thereby corrupting the results. My conclusion is that, under certain circumstances, the singular value corresponding to the constant mode can shrink so small that the SVDbased DMD misses it completely.
Given that this is a potential problem for some datasets, it is probably a good idea to preprocess your dataset by removing the mean of each signal in the data matrix. While your at it, you could divide each signal by its variance to completely normalize the data.
# compute DMD r = 5 mu,Phi = dmd(X, Y, r) check_dmd_result(X, Y, mu, Phi)
Here are visualizations of the five DMD modes. Each pseudocolor plot represents a different scalar function
And here are the time dynamics of the first three modes.
Limitations
The DMD has several known limitations. First of all, it doesn’t handle translational and rotational invariance particularly well. Secondly, it can fail completely when transient time behavior is present. The follow examples illustrate these problems.
Translational Invariance
The following dataset is very simple. It consists of a single mode (a Gaussian) translating along the spatial domain as the system evolves. Although one would think that the DMD would handle this cleanly, the opposite happens. Instead of picking up a single, welldefined singular value, the SVD picks up many.
# define time and space domains x = np.linspace(10, 10, 50) t = np.linspace(0, 10, 100) dt = t[2]  t[1] Xm,Tm = np.meshgrid(x, t) # create data with a single mode traveling both spatially and in time D = exp(power((XmTm+5)/2, 2)) D = D.T # extract inputoutput matrices X = D[:,:1] Y = D[:,1:] check_linear_consistency(X, Y)
The plot on the left shows the evolution of the system; the plot on the right shows its singular values. It turns out that close to ten DMD modes are needed to correctly approximate the system! Consider the following plots in which the original, true dynamics are compared to the superpositions of an increasing number of modes.
Transient time behavior
In this final example, we look at a dataset containing transient time dynamics. Specifically, the data shows a Gaussian popping in and out of existence. Unfortunately, the DMD cannot accurately decompose this data.
# define time and space domains x = np.linspace(10, 10, 50) t = np.linspace(0, 10, 100) dt = t[2]  t[1] Xm,Tm = np.meshgrid(x, t) # create data with a single transient mode D = exp(power((Xm)/4, 2)) * exp(power((Tm5)/2, 2)) D = D.astype('complex').T # extract inputoutput matrices X = D[:,:1] Y = D[:,1:] check_linear_consistency(X, Y) # do dmd r = 1 mu,Phi = dmd(X, Y, r) check_dmd_result(X, Y, mu, Phi) # time evolution b = dot(pinv(Phi), X[:,0]) Psi = np.zeros([r, len(t)], dtype='complex') for i,_t in enumerate(t): Psi[:,i] = multiply(power(mu, _t/dt), b)
Although the DMD correctly identifies the mode, it fails completely to identify the timebehavior. This is understandable if we consider that the timebehavior of the DMD time series depends on eigenvalues, which are only capable of characterizing combinations of exponential growth (the real part of the eigenvalue) and oscillation (the imaginary part).
The interesting thing about this system is that an ideal decomposition could potentially consist of a superposition of a single mode (as shown in the figure) with various eigenvalues. Imagine the single mode being multiplied by a linear combination of many orthogonal sines and cosines (a Fourier series) that approximates the true time dynamics. Unfortunately, a single application of the SVDbased DMD is incapable of yielding the same DMD mode various times with different eigenvalues.
Furthermore, it is important to note that even if we were able to correctly extract the time behavior as a large set of eigenvalues, the solution’s predictive capabilities would not be reliable without fully understanding the transient behavior itself. Transient behavior, by its very nature, is nonpermanent.
The multiresolution DMD (mrDMD) attempts to alleviate the transient time behavior problem by means of a recursive application of the DMD.
Conclusion
Despite its limitations, the DMD is a very powerful tool for analyzing and predicting dynamical systems. All data scientists from all backgrounds should have a good understanding of the DMD and how to apply it. The purpose of this post is to provide some theory behind the DMD and to provide practical code examples in Python that can be used with realworld data. After studying the formal definition of the DMD, walking through the algorithm stepbystep, and playing with several simple examples of its use – including examples in which it fails – it is my hope that this post has provided you an even clearer understanding of the DMD and how it might be applied in your research or engineering project.
There are many extensions to the DMD. Future work will likely result in posts on some of these extensions, including the multiresolution DMD (mrDMD) and the sparse DMD (sDMD).
References

Kutz, J. Nathan. Datadriven modeling & scientific computation: methods for complex systems & big data. OUP Oxford, 2013. ↩

Tu, Jonathan H., et al. “On dynamic mode decomposition: theory and applications.” arXiv preprint arXiv:1312.0041 (2013). ↩

Wikipedia contributors. “Composition operator.” Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, 28 Dec. 2015. Web. 24 Jul. 2016. ↩

Wikipedia contributors. “Mooreâ€“Penrose pseudoinverse.” Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, 12 May. 2016. Web. 24 Jul. 2016. ↩