Finding latent structures in data is an essential task in diverse fields, from signal processing through [biton2019optoacoustic] fluid dynamics analysis [ohmichi2017preconditioned] to machine learning [kaliroff2019self] . These structures commonly stand in the center of broad applications such as denoising, forecasting and classifying
. These structures commonly stand in the center of broad applications such as denoising, forecasting and classifying[shaham2018spectralnet] to name a few. These structures get different forms in different disciplines. For example, in image processing, the structures can be the probabilistic of repetitive patches in different scales [shaham2019singan], in signal processing they can be a sum of eigenfunctions [burger2016spectral], in fluid dynamics the structures are represented as a sum of modes [schmid2010dynamic] and in machine learning they might be recurrence of words [kuang2017crime]. Despite of this diversity, different techniques from different disciplines might share similar principles.
In this work, we formulate a signal decomposition framework, related to the eigenfunctions of -homogeneous operators, using insights from fluid dynamics and based on dynamic mode decomposition [schmid2010dynamic]. Our analysis shows that DMD produces inaccurate estimations when applied to uniformly sampled observations of nonlinear homogeneous flows. Thus, we suggest an adaptive sampling rate that allows DMD to fully recover the underlying dynamics. Moreover, as our flows are governed by operators with real-valued spectrum, we propose a new method named Symmetric DMD (S-DMD) which guarantees that the output matrix is symmetric.
Recently, a new signal decomposition, related the eigenfunction of the -homogeneous operator, (), was first introduced in [cohen2020Introducing]. This decomposition is based on the analysis of the solution to the PDE
It was found that the solution of this equation has finite support in time. Namely, the solution reaches its steady state in finite time. Moreover, the solution is variable separated when is an eigenfunction of , i.e. admits the nonlinear eigenvalue problem [cohen:hal-01870019]. Following these attributes of the solution, the authors of [cohen:hal-01870019] formulated a nonlinear spectral decomposition framework. This framework reveals the structures that stay unchanged under the flow. It generalized the previous study [gilboa2013spectral][gilboa2014total], focused on zero-homogeneous operators which are variational derivatives of absolutely one-homogeneous functionals.
The basis of this decomposition is the amplitude change rate of a signal under the flow, termed as the decay profile. Knowing the decay profile and the time of convergence to steady state, different modes of the signal can be distinguished. In addition, the decay profile might reveal the scale of the spatial structure. This concept was generalised in [Katzir2017Thesis, gilboa2018nonlinear]. In this work, the authors suggested to find a dictionary of decay profiles that compactly represents the dynamics. Relating different parts of the signal to the same decay profile reveals and connects between latent structures and scale within a signal.
In the fluid dynamics community, finding modes has drawn research attention for decades. A particularly successful approach in this context is based on the relation between Koopman theory [koopman1931hamiltonian, mezic2005spectral] and dynamical systems. In practice, an approximation of the linear yet infinite-dimensional Koopman operator can be computed using the Dynamic Mode Decomposition method [schmid2010dynamic] which combines dimensionality reduction with a linear solve. Several works have shown that DMD is an effective tool for analyzing various aspects of nonlinear flows including their decay profiles [askham2018variable, leroux2016dynamic, gueniat2015dynamic].
However, the baseline DMD algorithm was shown to be sensitive to noisy data [kutz2016dynamic]. Specifically, DMD produces a systematic bias in its spectrum estimates which is not alleviated when more data is present [hemati2017biasing]. The effect of small sensor noise on DMD and
the Koopman expansion has also been explored and characterized [bagheri2013effects]. Several attempts have been made to debias DMD with respect to noise including the use of the forward and backward dynamics [dawson2016characterizing], formulating the problem as a total least squares optimization [hemati2017biasing], and using a variational approach [azencot2019consistent]. The authors in [nonomura2018dynamic, nonomura2019extended] used Kalman filters to cope with the noise, and Williams et al.
used Kalman filters to cope with the noise, and Williams et al.[williams2015data] suggested to extend the basis of the sampled data snapshots, while [li2017extended][azencot2020forecasting]. In this work, we incorporate prior knowledge about the dynamics and their spectrum. In our approach, we utilize prior known information about the spectrum and force an algebraic form of the estimated operator, in some sense this approach is close to [askham2018variable] where exponential model was fitted. We formulate the Symmetric DMD (S-DMD), which is DMD when the matrix has only real eigenvalues.
Main contributions and structure of the paper
In this work, we study the connections between dynamic mode decomposition and nonlinear spectral analysis of -homogeneous flows. Our contributions can be summarized as follows.
First, we initialize a flow with a nonlinear eigenfunction, showing that DMD is not applicable in this setting. Nevertheless, a temporal re-parametrization of the data sampling leads to DMD results that are identical to previous spectral decompositions [cohen2020Introducing].
Second, we suggest a particular re-scaling of the time axis [cohen2019stable] such that a unique DMD solution is obtained. Our re-scaling yields an exponentially decaying solution, even though the operator is nonlinear.
Third, for dynamics with a real-valued spectrum, we improve the baseline DMD which may produce complex values in the noisy setting by proposing the S-DMD variant which yields symmetric operators.
Finally, we show how the nonlinear spectral framework of [cohen2020Introducing] can be approximated with our tools orders of magnitude faster compared to previous work.
Our paper is organized as follows. We briefly describe some of the necessary components and background in Sec. 2. Then, we formulate the problem of representing homogeneous flows with dynamic mode decomposition. In Sec. 3, we describe the adaptation to the sampling rate of the dynamics, and the S-DMD method that forces the eigenvalues to be real. We demonstrate how our method behaves on homogeneous and monotonic flows, and we show its robustness to noise in Sec. 4. We conclude our work and discuss a few future directions in Sec. 5.
Methods in image processing, based on functional analysis, is usually an optimization problem, compound from fidelity, , and regularization, , terms, i.e.
And the solution of is the recovered image. A common fidelity term is norm of the difference between the processed image, , and the given image, . Kuijper chose the regularization term as a norm of the gradient of [kuijper2007p]
known as the -Dirichlet energy. Thus, we can control the level of the smoothness by choosing . In addition, Kuijper suggested the following scale space as a filter
where the initial condition, , is the observed image and is the -Laplacian operator, defined as
This flow has a finite support in time for any initial condition when . Namely, it reaches its steady state in finite time for an arbitrary . Moreover, there exists function that holds
for some real non-positive number (see discussion [cohen2020Introducing]). These function are called eigenfunctions. If the initial condition is an eigenfunction of then the solution of (p-Flow) is variable separated and given by
where . This kind of solution we term as shape preserving flow. And since non-positive and when , one can calculate exactly what is the extinction time by (see Fig. (b)b)
In general, all the discussion above is valid for any -homogeneous operator when , when -homogeneous operator holds
for any . The -Laplacian operator is private case of -homogeneous operator, when .
2.1 The p-framework
The -framework is based on the decay profile of the -flow, eq. (4). This framework includes transform, inverse-transform, filtering, and spectrum. These analysis tools are based on the fractional derivatives and integrals of the solution of (p-Flow). The fractional derivatives and integrals are with respect to and of the order of . The transform is the fractional derivative of the solution of the PDE, eq. (p-Flow). The inverse-transform is defined by fractional integral of the same order of the transform. The filtering is fractional integral of the transform multiplied by . This function attenuates and enhances the transform in different scales. The spectrum is a function of , resulted from the inner product between the transform and the image (the initial condition of (p-Flow)). For more details we refer the reader to [cohen2020Introducing].
2.2 Explicit scheme of the p-flow
Applying the standard explicit scheme of the -flow raise two problems. The first is the step size. Until recently there was not CFL condition that allows to use this scheme. The second problem is when goes to zero. Then, the relation inside the divergence is not defined. These two problems were solved with regularization (adding small enough to the denominator) . In that case, the CFL is where is the dimension of the signal [baravdish2015backward]. However, this solution harm the analytic solution of the -flow. Such that, even for eigenfunction as an initial condition the solution is not decay polynomially.
On the other hand, Cohen et al. suggested to split the solution of these problems [cohen2019stable]. The limit of the fraction when goes to zero they solved with different calculation which avoids the dividing by zero. And as to the step size they put forward a generalized CFL condition. We summarize here their conclusions.
The explicit scheme of the dynamics (p-Flow) is given by
where is the initial condition. Analysing this scheme, they show empirically that this scheme is stable when the step size is constant (any constant). Moreover, they proved that this scheme converges to zero with the adaptive step size policy, given by
2.3 Dynamic Mode Decomposition (Dmd)
DMD is an analysis tool, reveling the main spatial structures in a fluid flow. By means of these structures, termed as modes, we can describe the entire dynamics. In what following, we summarize the main steps in the standard DMD [schmid2010dynamic] . Note, that we denote vectors with bullet font and matrices with capital letters.
. Note, that we denote vectors with bullet font and matrices with capital letters.
Forming data matrices
The data is a sequence of snapshots of the dynamic, denoted by . Each sample from this sequence, , belongs to . Form this sequence, we form two matrices that the th column in one matrix is preceding to the th column to the other, i.e.,
It is presumption, the data can be represented concisely but accurately with much lower dimension. Based on this assumption, Schmid used the Singular Vector Decomposition (SVD) to represent the matrix ; namely, decomposed this matrix as where is a diagonal matrix and the dimensions of are , and , respectively. The superscript denotes the conjugate transpose operator.
In this work, we denote as and as the submatrices, containing the first columns of and , respectively; and is the a submatrix of , containing the first rows and columns of . And the dimensionality reduction of the data is defined as
Then, the dimensionality reduced of the th snapshot is .
Mode, spectrum and coordination calculation
The relation between and can be approximated as,
where , termed as the DMD matrix. Then, the th sample can be expressed as
If the matrix is diagonalizable, i.e. can be writen as where the columns of , , are the right eigenvectors of
, are the right eigenvectors of, is a diagonal matrix containing the eigenvalues and the rows of , , are the left eigenvectors of . And in addition, the relation between and is . Thus, we can expand back to by . These identities and eq. (15) yields
The definitions of the modes, spectrum and coordinates result from eq. (16). The modes, denoted by , are defined as
The eigenvalues of , , is the spectrum of the dynamics. The coordinates are , where . The following algorithm, Algo. 1, summarizes the procedure of DMD.
3 Dmd for homogeneous and decreasing dynamics
In this section, we formulate a mode decomposition for flows that are governed by a homogeneous operator whose associated spectrum is real-valued. We will first show that DMD can not faithfully represent the dynamics of a -homogeneous operator , and we discuss an adaptation to the dynamics which allows for such a decomposition. We will then propose a modification to the baseline DMD approach to account for real-valued spectra even when noisy data samples are provided.
3.1 Dmd for a homogeneous flow
The dynamic mode decomposition method makes two underlying assumptions. First, the nonlinear dynamics can be represented linearly after dimensionality reduction and second, the data are uniformly sampled in time. Therefore, we can formulate the dynamics as defined in Alg. 1, State 5 [schmid2010dynamic]. Moreover, we can decompose a given initial condition as a linear combination of the eigenvectors of . The above structure allows Askham and Kutz [askham2018variable] to interpret DMD as an exponential fitting problem where the data are linear combinations of eigenvectors scaled by powers of the corresponding eigenvalues.
Other advanced algorithms of DMD lean only on one assumption, the linearity of the latent structure, whilst the data is sampled arbitrarily (known or unknown time points) [gueniat2015dynamic, leroux2016dynamic]. The underlying presumption in these articles is that the data is a linear combination of exponential functions. However, to best of our knowledge, no attention was drawn to flow that inherently can not be represented as a summation of exponential functions, for example, dynamics with finite support in time. This is emphasized in dynamics of the (p-Flow), where the decay profile is polynomial (Eq. (4)). Moreover, if the PDE (p-Flow) is initialised with an eigenfunction, trying to find a linear mapping to this nonlinear dynamics is doomed to failure. Namely, the DMD algorithm is ill pose since the Moore-Penrose pseudoinverse operator will be applied on a matrix with high condition number (analytically it should be infinite). This state holds even if after dimensionality reduction to any order. This problem is demonstrated in the following proposition.
Proposition 1 (Dmd an arbitrarily sampled -flow, initialized with an eigenfunction).
If the -flow is initialized with an eigenfunction then there is no linear mapping that relates to , unless is a geometric series.
Let us examine only three snapshots, . According to the linear mapping . On the other hand, according to (9) the constrain leads to . Therefore, to avoid inconsistent constraints the series must be geometric. ∎
Thus, the dynamic (9) can not be represented as a linear system unless the recurrence relating to (eq. (10)) is geometric. This recurrence is geometric under the adaptive step size policy (eq. (8)) [cohen2019stable]. And the solution of the explicit scheme is given by .
Proposition 2 (Dmd of the -flow, sampled adaptively and initialized with eigenfunction).
If the -flow is initialized with an eigenfunction and the time sampling, , is (8), then the eigenvalues are on circle with radius and only one mode is not zero.
Given relation we can rewrite it as
The characteristic polynomial of the matrix is
The eigenvalues of , which are also the eigenvalues of the matrix , are given by
And the corresponding eigenvectors are
As discussed in [rowley2009spectral] an eigenvector of is given by . Thus, these eigenvector are
As a result from this proposition, the linearity of the system is not necessary to the DMD to recover the system precisely. The exponential decay profile is enough in this case. Moreover, we can see that a adaptive step size policy is essential to DMD in dynamics derived from -homogeneous operator. Although the proof is done on -Laplacian operator this is correct to any -homogeneous operator, not only for but for any value of and the poof is trivial.
Generally, in continuous time, not only solutions of linear systems can be expressed as a summation of exponential functions. For example, one-homogeneous operators may induce exponential decay solution. The decay profile depends on the homogeneity order of the operator [cohen:hal-01870019]. If the initial condition, , is an eigenfunction then the decay profile is polynomial when , hyperbolic when and exponential when . However, non linear rescaling of the time variable makes this flow decays exponentially for any value of . For instance, we can rescale the time variable for the -flow by multiplying the -Laplacian operator by
Then, the new -flow, termed as the adaptive -flow, is defined by
One can see that the operator is one-homogeneous. Therefore, if the is an eigenfunction then the solution of (Adap-Flow) decays exponentially. Note, that is generalized Rayleigh quotient, as discussed in [cohen2019stable]. In the following theorem, we list some attributes of the solution of (Adap-Flow).
Theorem 1 (Decay profile and convergence of the adaptive -flow).
Let and be the -Dirichlet energy and the solution of (Adap-Flow). Then,
converges to zero exponentially.
If the initial condition, , is an eigenfunction of the solution is .
For an arbitrary initial condition, , every element from decays exponentially.
Using the chain rule of Brezis[brezis1973ope], we can write
The operator from eq. (Adap-Flow) is a one-homogeneous operator. Moreover, if is an eigenfunction of then the corresponding eigenvalue equals one. Then, the solution is [cohen2020Introducing] .
According to triangle inequality we can say that
Therefore, every element, , decays exponentially.
3.2 Dmd for a symmetric operator
According to the Theo. 1, third branch, one can conclude that every part of the initial condition, , decays exponentially under the adaptive step size policy. This decay is monotonically without any oscillations. In addition, it was discussed in [cohen2019stable] that the adaptive step size flow can yield negative eigenvalues but not complex. Thus, the demand on the DMD matrix to be symmetric and real is called for (a symmetric and real matrix has real eigenvalues).
This demand coincides with the analytic expression of non-linear diffusion in [weickert1998anisotropic, chapter 3.4]. In this monograph, Weickert investigated the following nonlinear ODE
where is a tensor
is a tensorfor continuous and semi-continuous settings. Weickert showed that it is possible to construct a semi-continuous model from this ODE that holds the following
where is symmetry.
By mean of this formulation of , we can embed this demand in the standard DMD. The following algorithm, Algo. 2, summarizes the DMD for real eigenvalues.
3.3 p-Decomposition using S-Dmd
The explicit scheme of the adaptive -flow is defined as
Note, that the is not bold since this is a general case of snapshot. Using Algo. 2, we can approximate this dynamic by
where is the DMD matrix, and are the samples and in coordination of the columns of , respectively. Thus, we can say that the dynamic can be approximated by
where is defined as and is the initial condition. Moreover, it is almost easy to see that the modes, , are the eigenvectors of the matrix with the same eigenvalues. We can write
On the other hand, the matrix approximates the dynamics, described in (20). Compering between these expressions, we can rewrite
We can rephrase this equation as follows
Then, we can say that is an approximation of an eigenfunction of the -Laplacian operator and the corresponding “eigenvalue”, , is
Thus, we can relate this mode to an “extinction” time, according to eq. (5).
To recap, we can approximate the -decomposition, formulated in [cohen2020Introducing]. Using the S-DMD with rescaled time axis, we can match a mode to an “extinction” time. By mean of that, we can define the discrete spectrum, where the horizon axis of the spectrum is the extinction time of the mode and the vertical axis is the coordinate related to the mode. Here, we reformulate the -framework, using the aforementioned approximation.
Definition 1 (Discrete -decomposition).
The discrete -decomposition of an image is the set modes .
As discussed above, the initial condition, , can be approximated as . The error, , depends on the dimensionality, [gavish2014optimal]. Note, that is an orthogonal set. Therefore, if the error is zero then this decomposition holds the Parseval’s identity. Namely, we can say .
Definition 2 (Discrete p-spectrum).
The discrete -spectrum of an image is the set, , where and are the extinction time and the coefficient of the mode , respectively.
Definition 3 (Discrete -filtering).
The discrete -filtering is amplification (or attenuation) of modes, i.e. the filtered image is give by
In this section, we demonstrate the contributions of this article and exhibit the theory we discuss above. First, we examine the robustness to noise of the S-DMD compared to DMD [schmid2010dynamic], tlsDMD [hemati2017biasing] and fbDMD [dawson2016characterizing]. Then, We showcase the results from Theo. 2 by applying the adaptive step size policy on the -flow. Then, we implement the discrete -framework via S-DMD. We compare between the running time of -decomposition, introduced in [cohen2020Introducing], and the discrete -decomposition Def. 2
Symmetric DMD (S-Dmd)
Here, we implement the S-DMD on a discrete stable linear system
The eigenvalues are and the initial condition is a normalized summation of the eigenvectors (namely ). We approximate the eigenvalues of this system based on 8 snapshots in presence of white Gaussian noise. We repeat our experiment times and average of each of the methods. In Fig. 9 we showcase the results and plot the ellipses which enclose the region of of the estimates that are closest to the true eigenvalue for each of the techniques (see [dawson2016characterizing]). One can see that for this kind of systems, DMD is superior on the tlsDMD and the fbDMD. In particularly, a method that takes into account the system and its inverse is doom to fail for every stable system since the inverse system is not stable. Therefore, whilst the fbDMD has good performances when the roots are on the unit cycle (BIBO stability) it fails when the roots are in the unit cycle.
Dmd derived from the adaptive step size policy
Let us recall that the solution of the -flow decay polynomially when it is initialized with an eigenfunction (see Fig. (b)b). And Theorem 1 shows that with a scale time policy the decay is exponential. Moreover, according to Prop. 2 is enough only one nonzero mode to describes accurately the solution of -flow initialised with an eigenfunction.
Here, we demonstrate the discrete decomposition we formulate in Def. 1. This decomposition, containing the triplets mode-time-coefficient , approximated the entire solution of (p-Flow). The approximation is a summation of a finite set of exponential functions. Each function is characterized with a shape (mode), decay rate (“extinction time”), and power (coefficient). Following these definitions we formulate the spectrum in Def. 2.
We bring it here as a reference to the approximated decomposition. The -spectrum, defined in Def. 2, is shown in Fig. (a)a. Then, the Figs. (b)b-(e)e are disjoint parts of the image. Each image represents different part of the spectrum. The image in Fig. (b)b represents the all parts in the spectrum in the segment (the blue part). The image in Fig. (c)c represents the all parts in the spectrum in the segment (the red part). Similarly, the images in Figs. (d)d and (e)e are related to the segments and (the yellow and the purple parts).
The approximation of this decomposition with the same value of is shown in Fig. 21. The discrete spectrum is shown in Fig. (f)f. The and axes indicate the time index and the power (squared coefficient) of the modes. Namely, the spectrum is the set of points . Then, we split the spectrum to four groups. The first part contains all the modes, whose “extinction time” hold (the blue part in the spectrum). This part of the spectrum is shown Fig. (g)g. The second group, shown in Fig. (h)h, contains all the modes with “extinction time” in the segment (the red part in the spectrum). In the same way, the images in Figs. (i)i and (j)j contain the mode with “extinction time” in and (yellow and purple parts in the spectrum), respectively. Generally, one can see that the S-DMD algorithm “fits” the texture and small scale components in the image to small “extinction time”.
Run time Vs. Image size
The most prominent advantage of this method is the running time. The -decomposition requires evaluating the (p-Flow) with uniform step size, which limit to small step size (the step size depends on the accuracy level). Then, for every pixel in the image fractional derivative is calculated. This might involve FFT and IFFT for every point in the image. In Fig. 22, we show the change rate in the computation time versus the size of the image.
The -axis indicates the size of the image (number of pixels) in log scale and -axis indicates the running time for compute the -spectrum in log scale. It seems that the running time increasing exponentially with the image size. However, the running time of the approximated -decomposition is smaller than the accurate one in at least order of magnitude.
The last part of the discrete -decomposition we would like to demonstrate here is the discrete filtering, formulated in Def. 3.
In Fig. 27 we demonstrate filtering with S-DMD. The initial condition (Fig. (a)a) is an eigenfunction, , with additional noise, . The discrete spectrum is brought in Fig. (d)d. One can see that different scales of the signal are represented in different places in the time axes. In Figs. (b)b and (c)c, we show the the blue and the red parts of the spectrum.
We test this algorithm again on a natural image (Fig. (a)a). We corrupt the image with additional noise, , (Fig. (b)b) . Filtering out the modes with the lower “extinction time”, we get the result in Fig. (c)c.
As expected, we loose some small details but the texture in general is preserved.
In this work, we introduce time rescaling to -homogeneous flow. This enables to implement the DMD on these flows. This rescale is essential to DMD and causes this decomposition to be much more accurate and even to recover the entire flow identically. In addition, we formulate the Symmetric DMD (S-DMD) to force the DMD matrix to be symmetric.
These two adaptations, the time rescale and the S-DMD, allow us to formulate the -framework approximation. This framework leans on the solution of (p-Flow), which is -homogeneous and symmetric. The adaptations of the sampling rate and to the DMD spectrum answer the these two attributes of the -flow. Thus, we can approximate the -decomposition, avoiding the numeric fractional derivative, necessary to this decomposition.
List of Symbols
|The system state in the th step in|
|matrix, approximating a linear mapping between and|
|Singular Vector Decomposition (SVD) of|
|dimensional reduced of , respectively|
|matrix,approximating a linear mapping between and|
|column vector, the th right eigenvectors of|
|row column, the th left eigenvectors of|
|are eigenvalues of .|
|are the modes of the dynamics.|
|are the coordinates of the modes respectively.|
|The -Dirichlet energy|
|The Laplacian operator|
|The inverse -transform|
Appendix A Finding a symmetric Dmd matrix
We are looking for linear mapping, , between and when the mapping is symmetric, i.e.
In addition, according to spectral theorem every symmetric real matrix can be diagonalized. Therefore, we can express the matrix as when and are orthogonal and diagonal matrices. Then, we can take squared root of the matrix and rewrite this expression as . Then, we can reformulate the optimization problem as
Note, that is over the complex field and denotes for the transform operator. Embedding the constrain in the optimization expression, we get
Using and the derivatives
we get that the minimizer, , holds
Substituting with , we get that the minimizer, , of (28) holds the following Sylvester equation