Kernel Embedding based Variational Approach for Low-dimensional Approximation of Dynamical Systems

It has been shown that complex nonlinear processes can be accurately described by linear models in many science and engineering fields, including wireless communications [19, 56], molecular dynamics [54, 4, 55], fluid dynamics [28, 43] and control theory [1], where the linear models can be expressed by a unified formula

and the expectation operator can be removed for deterministic systems. In such models, the state variable $\mathbf{x}$ is tranformed into a feature space by the transformation $\mathbf{f}(\mathbf{x})=\left(f_{1}(\mathbf{x}),\ldots,f_{m}(\mathbf{x})\right)^{\top}$, and the dynamics with lag time $\tau$ is characterized by a linear time-invariant system in the feature space. Then, all the dynamical properties of the system can be quantitatively analyzed after estimating the transition matrix $\mathbf{K}$ from data via linear regression. A special case of the linear models is Markov state models [39, 35] for conformational dynamics, which is equivalent to the well-known Ulam’s method [8, 14]. In a Markov state model, the feature mapping $\mathbf{f}$ consists of indicator functions of subsets of state and $\mathbf{K}=[K_{ij}]$ represents the transition probability from subset $i$ to subset $j$. Besides Markov state models and the Ulam’s method, a large number of similar modeling methods, e.g., dynamic mode decomposition [38, 2, 48, 21], time-lagged independent component analysis (TICA) [30, 34, 41], extended dynamic mode decomposition (EDMD) [50, 16, 40], Markov transition models [52], variational approach of conformation dynamics (VAC) [31, 32, 33], variational Koopman models [54] and their variants based on kernel ebmeddings [15, 16] and tensors [14, 33], are proposed by using different feature mappings.

From the perspective of operator theory, all the models in the form of (1) can be interpreted as algebraic representations of transfer operators of systems, including Frobenius-Perron (FP) operators and Koopman operators, and some of them are universal approximators for nonlinear dynamical systems under the assumption that the dimension of feature space is large enough [20] or the infinite-dimensional kernel mappings are utilized [45]. However, due to the limitation of computational cost and requirements of dynamical analysis, the low-dimensional approximation of transfer operators is still a critical and challenging problem in applications [17].

One common way to solve this problem is to identify the dominant dynamical structures, e.g., metastable states [9, 37], cycles [6] and coherent sets [11], and achieve the corresponding low-dimensional representations via spectral clustering. But this strategy assumes that an accurate high-dimensional model is known a priori, which is often violated especially for large-scale systems.

Another way for deterministic systems is to seek the feature mapping $\mathbf{f}$ by minimizing the regression error of (1) under the constraint that the state variable can also be accurately reconstructed from $\mathbf{f}$ [23, 24]. Notice that the constraint is necessary, otherwise a trivial but uninformative model with $\mathbf{f}(\mathbf{x})\equiv 1$ and $\mathbf{K}=1$ could be found. Some similar methods are developed for stochastic systems by considering (1) as a conditional generative model, where the parameters of $\mathbf{f}$ can be trained based on the likelihood or the other statistical criteria [51, 25]. However, these methods are applicable only if $\mathbf{f}$ are non-negative functions and usually involves the intractable probability density estimation.

In recent years, the variational approach has led to great progress for low-dimensional dynamical modeling, which was first proposed for time-reversible processes [31, 32, 27, 54] and extended to non-reversible processes in [53]. In contrast with the other methods, this approach provides a general and unified framework for data-driven model choice, reduction and optimization of dynamical systems based on the presented variational scores related to approximation errors of transfer operators. It can be easily integrated with deep learning to effectively analyze high-dimensional time series in an end-to-end manner [26, 3]. The existing variational principle based methods suffer from two drawbacks: First, it is necessary to assume that the transfer operator is Hilbert-Schimdt (HS) or compact as an operator between two weighted $\mathcal{L}^{2}$ spaces so that the maximum values of variational scores exist. But there is no easy way to test the assumption especially when we do not have strong prior knowledge regarding the system. Specifically, it can be proved that the assumption does not hold for most deterministic systems. Second, even for stochastic systems which satisfies the assumption, the common variational scores are possibly sensitive to small modeling variations, which could affect the effectiveness of the variational approach.

In this work, we introduce a kernel embedding based variational approach for dynamical systems (KVAD) using the theory of kernel embedding of functions [44, 46, 47, 45], where the modeling error is measured by using the distance between kernel embeddings of transition densities. The kernel based variational score in KVAD provides a robust and smooth quantification of differences between transfer operators, and is proved to be bounded for general dynamical systems, including deterministic and stochastic systems. Hence, it can effectively overcome the difficulties of existing variational methods, and expands significantly the range of applications. Like the previous variational scores, the kernel based score can also be consistently estimated from trajectory data without solving any intermediate statistical problem. Therefore, we develop a data-driven KVAD algorithm by considering $\mathbf{f}$ as a linear superposition of a given set of basis functions. Furthermore, we establish a relationship between KVAD, diffusion maps[5] and the singular components of transfer operators. Finally, the effectiveness the proposed algorithm is demonstrated by numerical experiments.

For a Markovian dynamical system in the state space $\mathbb{M}\subset\mathbb{R}^{D}$ , its dynamics can be completely characterized by the transition density

and the time evolution of the system state distribution can be formulated as

Here $\mathbf{x}_{t}$ denotes the state of the system at time $t$ and $p_{t}$ is the probability density of $\mathbf{x}_{t}$. The transfer operator $\mathcal{P}_{\tau}$ is called the Perron-Frobenius (PF) operator³³3Another commonly used transfer operator for Markovian dynamics is Koopman operator [49], which describes the evolution of observables instead of probability densities, and is the dual of the PF operator. In this paper, we focus only on the PF operator for convenience of analysis., which is a linear but usually infinite-dimensional operator. Notice that the deterministic dynamics in the form of $\mathbf{x}_{t+\tau}=\Theta_{\tau}(\mathbf{x}_{t})$ is a specific case of the Markovian dynamics, where $p_{\tau}(\mathbf{x},\cdot)=\delta_{\Theta_{\tau}(\mathbf{x})}(\cdot)$ is a Dirac function centered at $\Theta_{\tau}(\mathbf{x})$, and the corresponding PF operator is given by

By further assuming that the conditional distribution of $\mathbf{x}_{t+\tau}$ for given $\mathbf{x}_{t}=\mathbf{x}$ can always be represented by a linear combination of $m$ density basis functions $\mathbf{q}=\left(q_{1},\ldots,q_{m}\right)^{\top}:\mathbb{M}\to\mathbb{R}^{m}$, we obtain a finite-dimensional approximation of the transition density:

The feature mapping $\mathbf{f}(\mathbf{x})=\left(f_{1}(\mathbf{x}),\ldots,f_{m}(\mathbf{x})\right)^{\top}$ are real-valued observables of the state $\mathbf{x}_{t}=\mathbf{x}$, and provide a sufficient statistic for predicting the future states. Based on this approximation, the time evolution equation (3) of the state distribution can then be transformed into a linear evolution model of the feature functions $\mathbf{f}$ in the form of (1) with the transition matrix

and many dynamical properties of the Markov system, including spectral components, coherent sets and the stationary distribution, can be efficiently from the linear model.

It is shown in [20] that Eq. (4) provides a universal approximator of Markovian dynamics if the set of basis function is rich enough. But in this paper, we focus on a more practically problem: Given a small $m$, find $\mathbf{f}$ and $\mathbf{q}$ with $\mathrm{dim}(\mathbf{f})=\mathrm{dim}(\mathbf{q})=m$ such that the modeling error of (4) is minimized.

In this section, we develop a new variational principle for Markovian dynamics based on the kernel embedding of trainstion densities.

Assuming that $\kappa:\mathbb{M}\times\mathbb{M}\to\mathbb{R}$ is a universal kernel and bounded by $B$, $\mathcal{E}p_{\tau}(\mathbf{x},\cdot)$ is also bounded in $\mathbb{H}$ with

Motivated by this conclusion, we propose a new measure for approximation errors of PF operators

by replacing the $\mathcal{L}_{\rho_{1}^{-1}}^{2}$ norm with $\left\|\cdot\right\|_{\mathcal{E}}$ in (6), which is finite for both deterministic and stochastic systems if $\left\|\hat{p}_{\tau}(\mathbf{x},\cdot)\right\|_{\mathcal{E}}<\infty$. In contrast with Eq. (6), the new measure provides a more general way to quantify modeling errors of dynamics. Furthermore, from an application point of view, Eq. (7) provides a more smooth and effective representation of modeling errors of the conditional distributions.

The following proposition shows that Eq. (7) can also be derived from the HS norm of operator error by treating $\mathcal{P}_{\tau}$ as a mapping from $\mathcal{L}_{\rho_{0}^{-1}}^{2}$ to $\mathcal{L}_{\mathcal{E}}^{2}=\left\{q|\left\|q\right\|_{\mathcal{E}}^{2}<\infty\right\}$.

As a result of this proposition, we can find optimal $\mathbf{f}$ and $\mathbf{q}$ by maximizing $\mathcal{R}_{\mathcal{E}}$.

In this section, we derive a data-driven algorithm to estimate the optimal low-dimensional linear models based on the variational principle stated in Proposition 3.

In what follows, we demonstrate the benefits of the KVAD method for studies of nonlinear dynamical systems by two examples, and compare the results from KVAD with VAMP and kernel EDMD, where the basis functions in VAMP and kernel functions in kernel EDMD are the same as those in KVAD. For kernel EDMD, the low-dimensional linear model is achieved by leading eigenvalues and eigenfunctions as in [18], which characterizes invariant subspaces of systems.

In both examples, the feature mapping $\mathbf{f}$ is represented by basis functions

with all components of $\boldsymbol{\theta}_{i}$ and $b_{i}$ randomly drawn from $[-1,1]$, which are widely used in shallow neural networks [12, 7]. The kernel $\kappa$ is selected as the Gaussian kernel with $\sigma=1.5$ for the oscillator and $10$ for the Lorenz system.

Fig. 2 shows estimates of singular values of $\tilde{\mathcal{P}}_{\tau}:\mathcal{L}_{\rho_{0}^{-1}}^{2}\to\mathcal{L}_{\mathcal{E}}^{2}$ (KVAD), singular values of $\mathcal{P}_{\tau}:\mathcal{L}_{\rho_{0}^{-1}}^{2}\to\mathcal{L}_{\rho_{1}^{-1}}^{2}$ (VAMP) and absolute values of eigenvalues of $\mathcal{P}_{\tau}$ (kernel EDMD) with different noise parameters, where singular values must be nonnegative real numbers but eigenvalues could be complex or negative. We see that the singular values and eigenvalues given by VAMP and kernel EDMD decay very slowly. Hence, it is difficult to extract an accurate model from the estimation results of VAMP and kernel EDMD for a small $m$. Especially for VAMP, a large number of singular values are close to $1$ as analyzed in Proposition 1 when the systems are deterministic with $\xi=0$. In contrast, the singular values utilized in KVAD rapidly converges to zero, which implies one can effectively extract the essential part of dynamics from a small number of feature mappings.

The first two singular components of $\tilde{\mathcal{P}}_{\tau}$ for $\xi=0$ obtained by KVAD are shown in Fig. 3 (see Section 4.3).⁵⁵5The $q_{i}$ is approximated by multiple delta functions and hard to visualize. It can be observed that $f_{1},f_{2}$ of the oscillator characterize transitions between left-right and up-down areas separately. Those of the Lorenz systems are related to the two attractor lobes and the transition areas.

It is interesting to note that the singular values of $\tilde{\mathcal{P}}_{\tau}$ given by KVAD are slightly influenced by $\xi$ as illustrated in Fig. 2. Our numerical experience also show that the right singular functions remain almost unchanged for different $\xi$ (see Fig. 5 in Appendix I). This phenomenon can be partially explained by the fact that the variational score $\mathcal{R}_{\epsilon}$ estimated by (17) is not sensitive to small perturbations of $\mathbf{Y}$ if the bandwidth of the kernel is large. More thorough investigations on this phenomenon will be performed in future.

In order to quantitively evaluate the performance of the three methods, we define the following trajectory reconstruction error:

where $\mathbf{x}_{t}$ is the true trajectory data and $\mathbb{E}_{\mathrm{model}}[\mathbf{x}_{t}|\mathbf{x}_{0}]$ is the conditional mean value of $\mathbf{x}_{t}$ obtained by the model. The average error over multiple replicate simulations is minimized if and only if $\mathbb{E}_{\mathrm{model}}[\mathbf{x}_{l\tau}|\mathbf{x}_{0}]$ equals to the exact conditional mean value of $\mathbf{x}_{l\tau}$ for all $l$. For all the three methods,

where $\mathbf{G}$ is the least square solution to the regression problem $\mathbf{x}_{t+\tau}\approx\mathbf{G}^{\top}\mathbf{f}(\mathbf{x}_{t})$ [49]. Fig. 4 summarizes of reconstruction errors of the two systems obtained with different choices of the model dimension $m$ and noise parameter $\xi$, and the superiority of our KVAD is clearly shown.

In this paper, we combine the kernel embedding technique with the variational principle for transfer operators. This provides a powerful and flexible tool for low-dimensional approximation of dynamical system, and effectively addresses the shortcomings and limitations of the existing variational approach. In the proposed KVAD framework, a bounded and well defined distance measure of transfer operators is developed based on kernel embedding of transition densities, and the corresponding variational optimization approach can be applied to a broader range of dynamical systems than the existing variational approaches.

Our future work includes the convergence analysis of KVAD and optimization of kernel functions. From the algorithmic point of view, the main remaining question is how to efficiently perform KVAD learning from big data with deep neural networks. It will be also interesting to apply KVAD to multi-ensemble Markov models [55] for data analysis of enhanced sampling.