skcuda.linalg.dmd

skcuda.linalg.dmd(a_gpu, k=None, modes='exact', return_amplitudes=False, return_vandermonde=False, handle=None)

Dynamic Mode Decomposition.

Dynamic Mode Decomposition (DMD) is a data processing algorithm which allows to decompose a matrix a in space and time. The matrix a is decomposed as a = FBV, where the columns of F contain the dynamic modes. The modes are ordered corresponding to the amplitudes stored in the diagonal matrix B. V is a Vandermonde matrix describing the temporal evolution.

Parameters:

• a_gpu (pycuda.gpuarray.GPUArray) – Real/complex input matrix a with dimensions (m, n).
• k (int, optional) – If k < (n-1) low-rank Dynamic Mode Decomposition is computed.
• modes ({‘standard’, ‘exact’}) –

  ‘standard’ : uses the standard definition to compute the dynamic modes,

  F = U * W.

  ’exact’ : computes the exact dynamic modes, F = Y * V * (S**-1) * W.

• return_amplitudes (bool {True, False}) – True: return amplitudes in addition to dynamic modes.
• return_vandermonde (bool {True, False}) – True: return Vandermonde matrix in addition to dynamic modes and amplitudes.
• handle (int) – CUBLAS context. If no context is specified, the default handle from skcuda.misc._global_cublas_handle is used.

Returns:

• f_gpu (pycuda.gpuarray.GPUArray) – Matrix containing the dynamic modes of shape (m, n-1) or (m, k).
• b_gpu (pycuda.gpuarray.GPUArray) – 1-D array containing the amplitudes of length min(n-1, k).
• v_gpu (pycuda.gpuarray.GPUArray) – Vandermonde matrix of shape (n-1, n-1) or (k, n-1).

Notes

Double precision is only supported if the standard version of the CULA Dense toolkit is installed.

This function destroys the contents of the input matrix.

Arrays are assumed to be stored in column-major order, i.e., order=’F’.

References

M. R. Jovanovic, P. J. Schmid, and J. W. Nichols. “Low-rank and sparse dynamic mode decomposition.” Center for Turbulence Research Annual Research Briefs (2012): 139-152.

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

>>> #Numpy
>>> import numpy as np
>>> #Plot libs
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.mplot3d import Axes3D
>>> from matplotlib import cm
>>> #GPU DMD libs
>>> import pycuda.gpuarray as gpuarray
>>> import pycuda.autoinit
>>> from skcuda import linalg, rlinalg
>>> linalg.init()

>>> # Define time and space discretizations
>>> x=np.linspace( -15, 15, 200)
>>> t=np.linspace(0, 8*np.pi , 80)
>>> dt=t[2]-t[1]
>>> X, T = np.meshgrid(x,t)
>>> # Create two patio-temporal patterns
>>> F1 = 0.5* np.cos(X)*(1.+0.* T)
>>> F2 = ( (1./np.cosh(X)) * np.tanh(X)) *(2.*np.exp(1j*2.8*T))
>>> # Add both signals
>>> F = (F1+F2)

>>> #Plot dataset
>>> fig = plt.figure()
>>> ax = fig.add_subplot(231, projection='3d')
>>> ax = fig.gca(projection='3d')
>>> surf = ax.plot_surface(X, T, F, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=True)
>>> ax.set_zlim(-1, 1)
>>> plt.title('F')
>>> ax = fig.add_subplot(232, projection='3d')
>>> ax = fig.gca(projection='3d')
>>> surf = ax.plot_surface(X, T, F1, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
>>> ax.set_zlim(-1, 1)
>>> plt.title('F1')
>>> ax = fig.add_subplot(233, projection='3d')
>>> ax = fig.gca(projection='3d')
>>> surf = ax.plot_surface(X, T, F2, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
>>> ax.set_zlim(-1, 1)
>>> plt.title('F2')

>>> #Dynamic Mode Decomposition
>>> F_gpu = np.array(F.T, np.complex64, order='F')
>>> F_gpu = gpuarray.to_gpu(F_gpu)
>>> Fmodes_gpu, b_gpu, V_gpu, omega_gpu = linalg.dmd(F_gpu, k=2, modes='exact', return_amplitudes=True, return_vandermonde=True)
>>> omega = omega_gpu.get()
>>> plt.scatter(omega.real, omega.imag, marker='o', c='r')

>>> #Recover original signal
>>> F1tilde = np.dot(Fmodes_gpu[:,0:1].get() , np.dot(b_gpu[0].get(), V_gpu[0:1,:].get() ) )
>>> F2tilde = np.dot(Fmodes_gpu[:,1:2].get() , np.dot(b_gpu[1].get(), V_gpu[1:2,:].get() ) )

>>> #Plot DMD modes
>>> #Mode 0
>>> ax = fig.add_subplot(235, projection='3d')
>>> ax = fig.gca(projection='3d')
>>> surf = ax.plot_surface(X[0:F1tilde.shape[1],:], T[0:F1tilde.shape[1],:], F1tilde.T, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
>>> ax.set_zlim(-1, 1)
>>> plt.title('F1_tilde')
>>> #Mode 1
>>> ax = fig.add_subplot(236, projection='3d')
>>> ax = fig.gca(projection='3d')
>>> surf = ax.plot_surface(X[0:F2tilde.shape[1],:], T[0:F2tilde.shape[1],:], F2tilde.T, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
>>> ax.set_zlim(-1, 1)
>>> plt.title('F2_tilde')
>>> plt.show()