#!/usr/bin/env python """ PyCUDA-based randomized linear algebra functions. """ from __future__ import absolute_import, division from pprint import pprint from string import Template from pycuda.tools import context_dependent_memoize from pycuda.compiler import SourceModule from pycuda.reduction import ReductionKernel from pycuda import curandom from pycuda import cumath import pycuda.gpuarray as gpuarray import pycuda.driver as drv import pycuda.elementwise as el import pycuda.tools as tools import numpy as np from . import cublas from . import misc from . import linalg rand = curandom.MRG32k3aRandomNumberGenerator() import sys if sys.version_info < (3,): range = xrange class LinAlgError(Exception): """Randomized Linear Algebra Error.""" pass try: from . import cula _has_cula = True except (ImportError, OSError): _has_cula = False from .misc import init, add_matvec, div_matvec, mult_matvec from .linalg import hermitian, transpose # Get installation location of C headers: from . import install_headers def rsvd(a_gpu, k=None, p=0, q=0, method="standard", handle=None): """ Randomized Singular Value Decomposition. Randomized algorithm for computing the approximate low-rank singular value decomposition of a rectangular (m, n) matrix `a` with target rank `k << n`. The input matrix a is factored as `a = U * diag(s) * Vt`. The right singluar vectors are the columns of the real or complex unitary matrix `U`. The left singular vectors are the columns of the real or complex unitary matrix `V`. The singular values `s` are non-negative and real numbers. The paramter `p` is a oversampling parameter to improve the approximation. A value between 2 and 10 is recommended. The paramter `q` specifies the number of normlized power iterations (subspace iterations) to reduce the approximation error. This is recommended if the the singular values decay slowly and in practice 1 or 2 iterations achive good results. However, computing power iterations is increasing the computational time. If k > (n/1.5), partial SVD or trancated SVD might be faster. Parameters ---------- a_gpu : pycuda.gpuarray.GPUArray Real/complex input matrix `a` with dimensions `(m, n)`. k : int `k` is the target rank of the low-rank decomposition, k << min(m,n). p : int `p` sets the oversampling parameter (default k=0). q : int `q` sets the number of power iterations (default=0). method : `{'standard', 'fast'}` 'standard' : Standard algorithm as described in [1, 2] 'fast' : Version II algorithm as described in [2] handle : int CUBLAS context. If no context is specified, the default handle from `skcuda.misc._global_cublas_handle` is used. Returns ------- u_gpu : pycuda.gpuarray Right singular values, array of shape `(m, k)`. s_gpu : pycuda.gpuarray Singular values, 1-d array of length `k`. vt_gpu : pycuda.gpuarray Left singular values, array of shape `(k, n)`. 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'. Input matrix of shape `(m, n)`, where `n>m` is not supported yet. References ---------- N. Halko, P. Martinsson, and J. Tropp. "Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions" (2009). (available at `arXiv `_). S. Voronin and P.Martinsson. "RSVDPACK: Subroutines for computing partial singular value decompositions via randomized sampling on single core, multi core, and GPU architectures" (2015). (available at `arXiv `_). Examples -------- >>> import pycuda.gpuarray as gpuarray >>> import pycuda.autoinit >>> import numpy as np >>> from skcuda import linalg, rlinalg >>> linalg.init() >>> rlinalg.init() >>> #Randomized SVD decomposition of the square matrix `a` with single precision. >>> #Note: There is no gain to use rsvd if k > int(n/1.5) >>> a = np.array(np.random.randn(5, 5), np.float32, order='F') >>> a_gpu = gpuarray.to_gpu(a) >>> U, s, Vt = rlinalg.rsvd(a_gpu, k=5, method='standard') >>> np.allclose(a, np.dot(U.get(), np.dot(np.diag(s.get()), Vt.get())), 1e-4) True >>> #Low-rank SVD decomposition with target rank k=2 >>> a = np.array(np.random.randn(5, 5), np.float32, order='F') >>> a_gpu = gpuarray.to_gpu(a) >>> U, s, Vt = rlinalg.rsvd(a_gpu, k=2, method='standard') """ #************************************************************************* #*** Author: N. Benjamin Erichson *** #*** *** #*** License: BSD 3 clause *** #************************************************************************* if not _has_cula: raise NotImplementedError('CULA not installed') if handle is None: handle = misc._global_cublas_handle alloc = misc._global_cublas_allocator # The free version of CULA only supports single precision floating data_type = a_gpu.dtype.type real_type = np.float32 if data_type == np.complex64: cula_func_gesvd = cula.culaDeviceCgesvd cublas_func_gemm = cublas.cublasCgemm copy_func = cublas.cublasCcopy alpha = np.complex64(1.0) beta = np.complex64(0.0) TRANS_type = 'C' isreal = False elif data_type == np.float32: cula_func_gesvd = cula.culaDeviceSgesvd cublas_func_gemm = cublas.cublasSgemm copy_func = cublas.cublasScopy alpha = np.float32(1.0) beta = np.float32(0.0) TRANS_type = 'T' isreal = True else: if cula._libcula_toolkit == 'standard': if data_type == np.complex128: cula_func_gesvd = cula.culaDeviceZgesvd cublas_func_gemm = cublas.cublasZgemm copy_func = cublas.cublasZcopy alpha = np.complex128(1.0) beta = np.complex128(0.0) TRANS_type = 'C' isreal = False elif data_type == np.float64: cula_func_gesvd = cula.culaDeviceDgesvd cublas_func_gemm = cublas.cublasDgemm copy_func = cublas.cublasDcopy alpha = np.float64(1.0) beta = np.float64(0.0) TRANS_type = 'T' isreal = True else: raise ValueError('unsupported type') real_type = np.float64 else: raise ValueError('double precision not supported') #CUDA assumes that arrays are stored in column-major order m, n = np.array(a_gpu.shape, int) if n>m : raise ValueError('input matrix of shape (m,n), where n>m is not supported') #Set k if k == None : raise ValueError('k must be provided') if k > n or k < 1: raise ValueError('k must be 0 < k <= n') kt = k k = k + p if k > n: k=n #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Generate a random sampling matrix O #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if isreal==False: Oimag_gpu = gpuarray.empty((n,k), real_type, order="F", allocator=alloc) Oreal_gpu = gpuarray.empty((n,k), real_type, order="F", allocator=alloc) O_gpu = gpuarray.empty((n,k), data_type, order="F", allocator=alloc) rand.fill_uniform(Oimag_gpu) rand.fill_uniform(Oreal_gpu) O_gpu = Oreal_gpu + 1j * Oimag_gpu O_gpu = O_gpu.T * 2 - 1 #Scale to [-1,1] else: O_gpu = gpuarray.empty((n,k), real_type, order="F", allocator=alloc) rand.fill_uniform(O_gpu) #Draw random samples from a ~ Uniform(-1,1) distribution O_gpu = O_gpu * 2 - 1 #Scale to [-1,1] #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Build sample matrix Y : Y = A * O #Note: Y should approximate the range of A #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Allocate Y Y_gpu = gpuarray.zeros((m,k), data_type, order="F", allocator=alloc) #Dot product Y = A * O cublas_func_gemm(handle, 'n', 'n', m, k, n, alpha, a_gpu.gpudata, m, O_gpu.gpudata, n, beta, Y_gpu.gpudata, m ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Orthogonalize Y using economic QR decomposition: Y=QR #If q > 0 perfrom q subspace iterations #Note: economic QR just returns Q, and destroys Y_gpu #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if q > 0: Z_gpu = gpuarray.empty((n,k), data_type, order="F", allocator=alloc) for i in np.arange(1, q+1 ): if( (2*i-2)%q == 0 ): Y_gpu = linalg.qr(Y_gpu, 'economic', lib='cula') cublas_func_gemm(handle, TRANS_type, 'n', n, k, m, alpha, a_gpu.gpudata, m, Y_gpu.gpudata, m, beta, Z_gpu.gpudata, n ) if( (2*i-1)%q == 0 ): Z_gpu = linalg.qr(Z_gpu, 'economic', lib='cula') cublas_func_gemm(handle, 'n', 'n', m, k, n, alpha, a_gpu.gpudata, m, Z_gpu.gpudata, n, beta, Y_gpu.gpudata, m ) #End for #End if Q_gpu = linalg.qr(Y_gpu, 'economic', lib='cula') #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Project the data matrix a into a lower dimensional subspace #B = Q.T * A #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Allocate B B_gpu = gpuarray.empty((k,n), data_type, order="F", allocator=alloc) cublas_func_gemm(handle, TRANS_type, 'n', k, n, m, alpha, Q_gpu.gpudata, m, a_gpu.gpudata, m, beta, B_gpu.gpudata, k ) if method == 'standard': #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Singular Value Decomposition #Note: B = U" * S * Vt #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #gesvd(jobu, jobvt, m, n, int(a), lda, int(s), int(u), ldu, int(vt), ldvt) #Allocate s, U, Vt for economic SVD #Note: singular values are always real s_gpu = gpuarray.empty(k, real_type, order="F", allocator=alloc) U_gpu = gpuarray.empty((k,k), data_type, order="F", allocator=alloc) Vt_gpu = gpuarray.empty((k,n), data_type, order="F", allocator=alloc) #Economic SVD cula_func_gesvd('S', 'S', k, n, int(B_gpu.gpudata), k, int(s_gpu.gpudata), int(U_gpu.gpudata), k, int(Vt_gpu.gpudata), k) #Compute right singular vectors as U = Q * U" cublas_func_gemm(handle, 'n', 'n', m, k, k, alpha, Q_gpu.gpudata, m, U_gpu.gpudata, k, beta, Q_gpu.gpudata, m ) U_gpu = Q_gpu #Set pointer # Free internal CULA memory: cula.culaFreeBuffers() #Return return U_gpu[ : , 0:kt ], s_gpu[ 0:kt ], Vt_gpu[ 0:kt , : ] elif method == 'fast': #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Orthogonalize B.T using reduced QR decomposition: B.T = Q" * R" #Note: reduced QR returns Q and R, and destroys B_gpu #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if isreal==True: B_gpu = transpose(B_gpu) #transpose B else: B_gpu = hermitian(B_gpu) #transpose B Qstar_gpu, Rstar_gpu = linalg.qr(B_gpu, 'reduced', lib='cula') #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Singular Value Decomposition of R" #Note: R" = U" * S" * Vt" #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #gesvd(jobu, jobvt, m, n, int(a), lda, int(s), int(u), ldu, int(vt), ldvt) #Allocate s, U, Vt for economic SVD #Note: singular values are always real s_gpu = gpuarray.empty(k, real_type, order="F", allocator=alloc) Ustar_gpu = gpuarray.empty((k,k), data_type, order="F", allocator=alloc) Vtstar_gpu = gpuarray.empty((k,k), data_type, order="F", allocator=alloc) #Economic SVD cula_func_gesvd('A', 'A', k, k, int(Rstar_gpu.gpudata), k, int(s_gpu.gpudata), int(Ustar_gpu.gpudata), k, int(Vtstar_gpu.gpudata), k) #Compute right singular vectors as U = Q * Vt.T" cublas_func_gemm(handle, 'n', TRANS_type, m, k, k, alpha, Q_gpu.gpudata, m, Vtstar_gpu.gpudata, k, beta, Q_gpu.gpudata, m ) U_gpu = Q_gpu #Set pointer #Compute left singular vectors as Vt = U".T * Q".T Vt_gpu = gpuarray.empty((k,n), data_type, order="F", allocator=alloc) cublas_func_gemm(handle, TRANS_type, TRANS_type, k, n, k, alpha, Ustar_gpu.gpudata, k, Qstar_gpu.gpudata, n, beta, Vt_gpu.gpudata, k ) # Free internal CULA memory: cula.culaFreeBuffers() #Return return U_gpu[ : , 0:kt ], s_gpu[ 0:kt ], Vt_gpu[ 0:kt , : ] #End if def rdmd(a_gpu, k=None, p=5, q=1, modes='exact', method_rsvd='standard', return_amplitudes=False, return_vandermonde=False, handle=None): """ Randomized 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. p : int `p` sets the oversampling parameter for rSVD (default k=5). q : int `q` sets the number of power iterations for rSVD (default=1). 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`. method_rsvd : `{'standard', 'fast'}` 'standard' : (default) Standard algorithm as described in [1, 2] 'fast' : Version II algorithm as described in [2] 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 ---------- N. B. Erichson and C. Donovan. "Randomized Low-Rank Dynamic Mode Decomposition for Motion Detection" Under Review. N. Halko, P. Martinsson, and J. Tropp. "Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions" (2009). (available at `arXiv `_). J. H. Tu, et al. "On dynamic mode decomposition: theory and applications." arXiv preprint arXiv:1312.0041 (2013). Examples -------- >>> #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() >>> rlinalg.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 = rlinalg.rdmd(F_gpu, k=2, p=0, q=1, 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() """ #************************************************************************* #*** Author: N. Benjamin Erichson *** #*** <2015> *** #*** License: BSD 3 clause *** #************************************************************************* if not _has_cula: raise NotImplementedError('CULA not installed') if handle is None: handle = misc._global_cublas_handle alloc = misc._global_cublas_allocator # The free version of CULA only supports single precision floating data_type = a_gpu.dtype.type real_type = np.float32 if data_type == np.complex64: cula_func_gesvd = cula.culaDeviceCgesvd cublas_func_gemm = cublas.cublasCgemm cublas_func_dgmm = cublas.cublasCdgmm cula_func_gels = cula.culaDeviceCgels copy_func = cublas.cublasCcopy transpose_func = cublas.cublasCgeam alpha = np.complex64(1.0) beta = np.complex64(0.0) TRANS_type = 'C' isreal = False elif data_type == np.float32: cula_func_gesvd = cula.culaDeviceSgesvd cublas_func_gemm = cublas.cublasSgemm cublas_func_dgmm = cublas.cublasSdgmm cula_func_gels = cula.culaDeviceSgels copy_func = cublas.cublasScopy transpose_func = cublas.cublasSgeam alpha = np.float32(1.0) beta = np.float32(0.0) TRANS_type = 'T' isreal = True else: if cula._libcula_toolkit == 'standard': if data_type == np.complex128: cula_func_gesvd = cula.culaDeviceZgesvd cublas_func_gemm = cublas.cublasZgemm cublas_func_dgmm = cublas.cublasZdgmm cula_func_gels = cula.culaDeviceZgels copy_func = cublas.cublasZcopy transpose_func = cublas.cublasZgeam alpha = np.complex128(1.0) beta = np.complex128(0.0) TRANS_type = 'C' isreal = False elif data_type == np.float64: cula_func_gesvd = cula.culaDeviceDgesvd cublas_func_gemm = cublas.cublasDgemm cublas_func_dgmm = cublas.cublasDdgmm cula_func_gels = cula.culaDeviceDgels copy_func = cublas.cublasDcopy transpose_func = cublas.cublasDgeam alpha = np.float64(1.0) beta = np.float64(0.0) TRANS_type = 'T' isreal = True else: raise ValueError('unsupported type') real_type = np.float64 else: raise ValueError('double precision not supported') #CUDA assumes that arrays are stored in column-major order m, n = np.array(a_gpu.shape, int) nx = n-1 #Set k if k == None : k = nx if k > nx or k < 1: raise ValueError('k is not valid') #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Split data into lef and right snapshot sequence #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Note: we need a copy of X_gpu, because SVD destroys X_gpu #While Y_gpu is just a pointer X_gpu = gpuarray.empty((m, n), data_type, order="F", allocator=alloc) copy_func(handle, X_gpu.size, int(a_gpu.gpudata), 1, int(X_gpu.gpudata), 1) X_gpu = X_gpu[:, :nx] Y_gpu = a_gpu[:, 1:] #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Randomized Singular Value Decomposition #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ U_gpu, s_gpu, Vh_gpu = rsvd(X_gpu, k=k, p=p, q=q, method=method_rsvd, handle=handle) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Solve the LS problem to find estimate for M using the pseudo-inverse #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #real: M = U.T * Y * Vt.T * S**-1 #complex: M = U.H * Y * Vt.H * S**-1 #Let G = Y * Vt.H * S**-1, hence M = M * G #Allocate G and M G_gpu = gpuarray.empty((m,k), data_type, order="F", allocator=alloc) M_gpu = gpuarray.empty((k,k), data_type, order="F", allocator=alloc) #i) s = s **-1 (inverse) if data_type == np.complex64 or data_type == np.complex128: s_gpu = 1/s_gpu s_gpu = s_gpu + 1j * gpuarray.zeros_like(s_gpu) else: s_gpu = 1.0/s_gpu #ii) real/complex: scale Vs = Vt* x diag(s**-1) Vs_gpu = gpuarray.empty((nx,k), data_type, order="F", allocator=alloc) lda = max(1, Vh_gpu.strides[1] // Vh_gpu.dtype.itemsize) ldb = max(1, Vs_gpu.strides[1] // Vs_gpu.dtype.itemsize) transpose_func(handle, TRANS_type, TRANS_type, nx, k, alpha, int(Vh_gpu.gpudata), lda, beta, int(Vh_gpu.gpudata), lda, int(Vs_gpu.gpudata), ldb) cublas_func_dgmm(handle, 'r', nx, k, int(Vs_gpu.gpudata), nx, int(s_gpu.gpudata), 1 , int(Vs_gpu.gpudata), nx) #iii) real: G = Y * Vs , complex: G = Y x Vs cublas_func_gemm(handle, 'n', 'n', m, k, nx, alpha, int(Y_gpu.gpudata), m, int(Vs_gpu.gpudata), nx, beta, int(G_gpu.gpudata), m ) #iv) real/complex: M = U* x G cublas_func_gemm(handle, TRANS_type, 'n', k, k, m, alpha, int(U_gpu.gpudata), m, int(G_gpu.gpudata), m, beta, int(M_gpu.gpudata), k ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Eigen Decomposition #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Note: If a_gpu is real the imag part is omitted Vr_gpu, w_gpu = linalg.eig(M_gpu, 'N', 'V', 'F', lib='cula') omega = cumath.log(w_gpu) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Compute DMD Modes #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ F_gpu = gpuarray.empty((m,k), data_type, order="F", allocator=alloc) modes = modes.lower() if modes == 'exact': #Compute (exact) DMD modes: F = Y * V * S**-1 * W = G * W cublas_func_gemm(handle, 'n', 'n', m, k, k, alpha, G_gpu.gpudata, m, Vr_gpu.gpudata, k, beta, G_gpu.gpudata, m ) F_gpu_temp = G_gpu elif modes == 'standard': #Compute (standard) DMD modes: F = U * W cublas_func_gemm(handle, 'n', 'n', m, k, k, alpha, U_gpu.gpudata, m, Vr_gpu.gpudata, k, beta, U_gpu.gpudata, m ) F_gpu_temp = U_gpu else: raise ValueError('Type of modes is not supported, choose "exact" or "standard".') #Copy is required, because gels destroys input copy_func(handle, F_gpu_temp.size, int(F_gpu_temp.gpudata), 1, int(F_gpu.gpudata), 1) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Compute amplitueds b using least-squares: Fb=x1 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if return_amplitudes==True: #x1_gpu = a_gpu[:,0].copy() x1_gpu = gpuarray.empty(m, data_type, order="F", allocator=alloc) copy_func(handle, x1_gpu.size, int(a_gpu[:,0].gpudata), 1, int(x1_gpu.gpudata), 1) cula_func_gels( 'N', m, k, int(1) , F_gpu_temp.gpudata, m, x1_gpu.gpudata, m) b_gpu = x1_gpu #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Compute Vandermonde matrix (CPU) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if return_vandermonde==True: V_gpu = linalg.vander(w_gpu, n=nx) # Free internal CULA memory: cula.culaFreeBuffers() #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Return #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if return_amplitudes==True and return_vandermonde==True: return F_gpu, b_gpu[:k], V_gpu, omega elif return_amplitudes==True and return_vandermonde==False: return F_gpu, b_gpu[:k], omega elif return_amplitudes==False and return_vandermonde==True: return F_gpu, V_gpu, omega else: return F_gpu, omega def cdmd(a_gpu, k=None, c=None, modes='exact', return_amplitudes=False, return_vandermonde=False, handle=None): """ Compressed 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. c : int `p` sets the number of measurements sensors. modes : `{'exact'}` '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 ---------- S. L. Brunton, et al. "Compressed sampling and dynamic mode decomposition." arXiv preprint arXiv:1312.5186 (2013). J. H. Tu, et al. "On dynamic mode decomposition: theory and applications." arXiv preprint arXiv:1312.0041 (2013). Examples -------- >>> #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() >>> rlinalg.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 = rlinalg.cdmd(F_gpu, k=2, c=20, 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() """ #************************************************************************* #*** Author: N. Benjamin Erichson *** #*** <2015> *** #*** License: BSD 3 clause *** #************************************************************************* if not _has_cula: raise NotImplementedError('CULA not installed') if handle is None: handle = misc._global_cublas_handle alloc = misc._global_cublas_allocator # The free version of CULA only supports single precision floating data_type = a_gpu.dtype.type real_type = np.float32 if data_type == np.complex64: cula_func_gesvd = cula.culaDeviceCgesvd cublas_func_gemm = cublas.cublasCgemm cublas_func_dgmm = cublas.cublasCdgmm cula_func_gels = cula.culaDeviceCgels copy_func = cublas.cublasCcopy transpose_func = cublas.cublasCgeam alpha = np.complex64(1.0) beta = np.complex64(0.0) TRANS_type = 'C' isreal = False elif data_type == np.float32: cula_func_gesvd = cula.culaDeviceSgesvd cublas_func_gemm = cublas.cublasSgemm cublas_func_dgmm = cublas.cublasSdgmm cula_func_gels = cula.culaDeviceSgels copy_func = cublas.cublasScopy transpose_func = cublas.cublasSgeam alpha = np.float32(1.0) beta = np.float32(0.0) TRANS_type = 'T' isreal = True else: if cula._libcula_toolkit == 'standard': if data_type == np.complex128: cula_func_gesvd = cula.culaDeviceZgesvd cublas_func_gemm = cublas.cublasZgemm cublas_func_dgmm = cublas.cublasZdgmm cula_func_gels = cula.culaDeviceZgels copy_func = cublas.cublasZcopy transpose_func = cublas.cublasZgeam alpha = np.complex128(1.0) beta = np.complex128(0.0) TRANS_type = 'C' isreal = False elif data_type == np.float64: cula_func_gesvd = cula.culaDeviceDgesvd cublas_func_gemm = cublas.cublasDgemm cublas_func_dgmm = cublas.cublasDdgmm cula_func_gels = cula.culaDeviceDgels copy_func = cublas.cublasDcopy transpose_func = cublas.cublasDgeam alpha = np.float64(1.0) beta = np.float64(0.0) TRANS_type = 'T' isreal = True else: raise ValueError('unsupported type') real_type = np.float64 else: raise ValueError('double precision not supported') #CUDA assumes that arrays are stored in column-major order m, n = np.array(a_gpu.shape, int) nx = n-1 #Set k if k == None : k = nx if k > nx or k < 1: raise ValueError('k is not valid') #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Compress #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if c==None: Ac_gpu = A c=m else: #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Generate a random sensing matrix S #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if isreal==False: Simag_gpu = gpuarray.empty((m,c), real_type, order="F", allocator=alloc) Sreal_gpu = gpuarray.empty((m,c), real_type, order="F", allocator=alloc) S_gpu = gpuarray.empty((c,m), data_type, order="F", allocator=alloc) rand.fill_uniform(Simag_gpu) rand.fill_uniform(Sreal_gpu) S_gpu = Sreal_gpu + 1j * Simag_gpu S_gpu = S_gpu.T * 2 -1 #Scale to [-1,1] else: S_gpu = gpuarray.empty((c,m), real_type, order="F", allocator=alloc) rand.fill_uniform(S_gpu) #Draw random samples from a ~ Uniform(-1,1) distribution S_gpu = S_gpu * 2 - 1 #Scale to [-1,1] #Allocate Ac Ac_gpu = gpuarray.empty((c,n), data_type, order="F", allocator=alloc) #Compress input matrix cublas_func_gemm(handle, 'n', 'n', c, n, m, alpha, int(S_gpu.gpudata), c, int(a_gpu.gpudata), m, beta, int(Ac_gpu.gpudata), c ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Split data into lef and right snapshot sequence #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Note: we need a copy of X_gpu, because SVD destroys X_gpu #While Y_gpu is just a pointer X_gpu = gpuarray.empty((c, n), data_type, order="F", allocator=alloc) copy_func(handle, X_gpu.size, int(Ac_gpu.gpudata), 1, int(X_gpu.gpudata), 1) X_gpu = X_gpu[:, :nx] Y_gpu = Ac_gpu[:, 1:] Yorig_gpu = a_gpu[:, 1:] #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Singular Value Decomposition #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Allocate s, U, Vt for economic SVD #Note: singular values are always real min_s = min(nx,c) s_gpu = gpuarray.zeros(min_s, real_type, order="F", allocator=alloc) U_gpu = gpuarray.zeros((c,min_s), data_type, order="F", allocator=alloc) Vh_gpu = gpuarray.zeros((min_s,nx), data_type, order="F", allocator=alloc) #Economic SVD cula_func_gesvd('S', 'S', c, nx, int(X_gpu.gpudata), c, int(s_gpu.gpudata), int(U_gpu.gpudata), c, int(Vh_gpu.gpudata), min_s) #Low-rank DMD: trancate SVD if k < nx if k != nx: s_gpu = s_gpu[:k] U_gpu = U_gpu[: , :k] Vh_gpu = Vh_gpu[:k , : ] #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Solve the LS problem to find estimate for M using the pseudo-inverse #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #real: M = U.T * Y * Vt.T * S**-1 #complex: M = U.H * Y * Vt.H * S**-1 #Let G = Y * Vt.H * S**-1, hence M = M * G #Allocate G and M G_gpu = gpuarray.zeros((c,k), data_type, order="F", allocator=alloc) M_gpu = gpuarray.zeros((k,k), data_type, order="F", allocator=alloc) #i) s = s **-1 (inverse) if data_type == np.complex64 or data_type == np.complex128: s_gpu = 1/s_gpu s_gpu = s_gpu + 1j * gpuarray.zeros_like(s_gpu) else: s_gpu = 1/s_gpu #ii) real/complex: scale Vs = Vt* x diag(s**-1) Vs_gpu = gpuarray.zeros((nx,k), data_type, order="F", allocator=alloc) lda = max(1, Vh_gpu.strides[1] // Vh_gpu.dtype.itemsize) ldb = max(1, Vs_gpu.strides[1] // Vs_gpu.dtype.itemsize) transpose_func(handle, TRANS_type, TRANS_type, nx, k, 1.0, int(Vh_gpu.gpudata), lda, 0.0, int(Vh_gpu.gpudata), lda, int(Vs_gpu.gpudata), ldb) #End Transpose cublas_func_dgmm(handle, 'r', nx, k, int(Vs_gpu.gpudata), nx, int(s_gpu.gpudata), 1 , int(Vs_gpu.gpudata), nx) #iii) real: G = Y * Vs , complex: G = Y x Vs cublas_func_gemm(handle, 'n', 'n', c, k, nx, alpha, int(Y_gpu.gpudata), c, int(Vs_gpu.gpudata), nx, beta, int(G_gpu.gpudata), c ) #iv) real/complex: M = U* x G cublas_func_gemm(handle, TRANS_type, 'n', k, k, c, alpha, int(U_gpu.gpudata), c, int(G_gpu.gpudata), c, beta, int(M_gpu.gpudata), k ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Eigen Decomposition #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Note: If a_gpu is real the imag part is omitted Vr_gpu, w_gpu = linalg.eig(M_gpu, 'N', 'V', 'F', lib='cula') omega = cumath.log(w_gpu) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Compute DMD Modes #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ F_gpu = gpuarray.empty((m,k), data_type, order="F", allocator=alloc) modes = modes.lower() if modes == 'exact': #Compute (exact) DMD modes: F = Y * V * S**-1 * W = G * W cublas_func_gemm(handle, 'n' , 'n', nx, k, k, alpha, int(Vs_gpu.gpudata), nx, int(Vr_gpu.gpudata), k, beta, int(Vs_gpu.gpudata), nx ) cublas_func_gemm(handle, 'n', 'n', m, k, nx, alpha, Yorig_gpu.gpudata, m, Vs_gpu.gpudata, nx, beta, F_gpu.gpudata, m ) else: raise ValueError('Type of modes is not supported, choose "exact" or "standard".') #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Compute amplitueds b using least-squares: Fb=x1 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if return_amplitudes==True: F_gpu_temp = gpuarray.empty((m,k), data_type, order="F", allocator=alloc) #Copy is required, because gels destroys input copy_func(handle, F_gpu.size, int(F_gpu.gpudata), 1, int(F_gpu_temp.gpudata), 1) #x1_gpu = a_gpu[:,0].copy() x1_gpu = gpuarray.empty(m, data_type, order="F", allocator=alloc) copy_func(handle, x1_gpu.size, int(a_gpu[:,0].gpudata), 1, int(x1_gpu.gpudata), 1) cula_func_gels( 'N', m, k, int(1) , F_gpu_temp.gpudata, m, x1_gpu.gpudata, m) b_gpu = x1_gpu #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Compute Vandermonde matrix (CPU) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if return_vandermonde==True: V_gpu = linalg.vander(w_gpu, n=nx) # Free internal CULA memory: cula.culaFreeBuffers() #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Return #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if return_amplitudes==True and return_vandermonde==True: return F_gpu, b_gpu[:k], V_gpu, omega elif return_amplitudes==True and return_vandermonde==False: return F_gpu, b_gpu[:k], omega elif return_amplitudes==False and return_vandermonde==True: return F_gpu, V_gpu, omega else: return F_gpu, omega if __name__ == "__main__": import doctest doctest.testmod()