On the Stochastic Stability of Deep Markov Models
Year: 2022

Abstract:
Deep Markov models (DMM) are generative models that are scalable and expressive generalization of Markov models for representation, learning, and inference problems. However, the fundamental stochastic stability guarantees of such models have not been thoroughly investigated. In this paper, we provide sufficient conditions of DMM’s stochastic stability as defined in the context of dynamical systems and propose a stability analysis method based on the contraction of probabilistic maps modeled by deep neural networks. We make connections between the spectral properties of neural network’s weights and different types of used activation functions on the stability and overall dynamic behavior of DMMs with Gaussian distributions. Based on the theory, we propose a few practical methods for designing constrained DMMs with guaranteed stability. We empirically substantiate our theoretical results via intuitive numerical experiments using the proposed stability constraints.

heading: 1 Introduction
text: Modeling, analysis, and control of dynamical systems are of utmost importance for various physical and engineered systems such as fluid dynamics, oscillators, power grids, transportation networks, and autonomous driving, to name just a few. The systems are generally subjected to uncertainties arising from a plethora of factors such as exogenous noise, plant-model mismatch, and unmodeled system dynamics, which have led researchers to model the dynamics in stochastic frameworks. One of the most commonly used probabilistic frameworks to model dynamical system is the Hidden Markov Models (HMMs) [Rabiner and Juang, 1986, Eddy, 1996] which in their vanilla form have been extensively investigated for representation, learning, and inference problems [Ghahramani and Jordan, 1997, Cappé et al., 2006, Beal et al., 2002]. One of its variants, the Gaussian state space models, have been used in the systems and control community for decades [Beckers and Hirche, 2016, Eleftheriadis et al., 2017].
It has been shown that expressivity of Markov models to emulate complex dynamics and sequential behaviors is greatly improved by parametrizing such models using deep neural networks, giving rise to Deep Markov Models (DMMs) [Krishnan et al., 2017]. Main research activities have been focused on the inference of these models. In particular, works such as Awiszus and Rosenhahn [2018], Liu et al. [2019], Mustafa et al. [2019], Qu et al. [2019] proposed parametrizing the probability distributions using deep neural networks for modeling various complex dynamical systems. Despite the rising popularity of DMMs, many of their theoretical properties such as robustness to perturbations, and stability of the generated trajectories remain open research questions. Many natural systems exhibit complex yet stable dynamical behavior that is described by converging trajectories towards an attractor set [Brayton and Tong, 1979]. Additionally, safety-critical systems such as an autonomous
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
driving call for formal verification methods to ensure safe operation. Thus the ability to guarantee stability in data-driven models could lead to improved generalization or act as a safety certificate in real-world applications.
In this paper, we propose a new analytical method to assess stochastic stability of DMMs. More specifically, we utilize spectral analysis of deep neural networks (DNNs) modeling DMM’s distributions. This allows us to make connections between the stability of deterministic DNNs and the stochastic stability of deep Markov models. As a main theoretical contribution we provide sufficient conditions for stochastic stability of DMMs. Based on the proposed theory we introduce several practical methods for design of constrained DMM with stability guarantees. In summary, the main contributions of this paper include:
1. Stability analysis method for deep Markov models: we base the analysis on the operator norms of deep neural networks modeling mean and variance of the DMM’s distributions. 2. Sufficient conditions for stochastic stability of deep Markov models: we show the sufficiency of the operator norm-based contraction conditions for DMM’s deep neural networks by leveraging Banach fixed point theorem. 3. Stability constrained deep Markov models: we introduce a set of methods for the design of provably stable deep Markov models. 4. Numerical case studies: we analyze connections between the design parameters of neural networks, stochastic stability, and operator norms of deep Markov models.

heading: 2 Related Work
text: Deep Markov models (DMMs) have been used as a scalable and expressive generalization of Hidden Markov Models (HMM) for learning probabilistic generative models of complex high-dimensional dynamical systems from sequential data [Rezende et al., 2014, Krishnan et al., 2017, Fraccaro et al., 2016]. DMMs have been successfully applied to speech recognition problems [Li et al., 2013, Prasetio et al.], control problems [Shashua and Mannor, 2017], human pose forecasting [Toyer et al., 2017], fault detection [Wang et al., 2018], climate data forecasting [Che et al., 2018], molecular dynamics [Wu et al., 2018], or as internal models in model-based deep RL applied to automatic trading [Ferreira, 2020]. Several modifications of DMMs have been proposed to handle incomplete data [Tan et al., 2019], or multi-rate time series [Che et al., 2018] multivariate time series [Montanez et al., 2015], training DMMs in unsupervised settings [Tran et al., 2016], or architectures inspired by Kalman Filters [Krishnan et al., 2015, Shashua and Mannor, 2017, Becker et al., 2019]. However, works focusing on formal analysis to ensure stability guarantees for DMMs are missing.
Stability notions and analysis for stochastic dynamic systems have been studied in the automatic control literature in various forms, depending on the representation of the system dynamics. Some classical results on stochastic stability for analysis and control can be found in McLane [1971], Willems and Willems [1976]. These results are presented mainly for stochastic differential equation (SDE) models. Khasminskii [2012] discusses different notions of stochastic stability, among them mean-square based stability notions have gained interest in works such as Lu and Skelton [2002], Farmer et al. [2009], Elia et al. [2013], Nandanoori et al. [2018], Wu et al. [2019]. We resort to such mean-square based stability notions when analyzing probabilistic state transition models parametrized by deep neural networks. We show that when the stochastic transitions are modeled by DNNs, the probabilistic stability requirements can be translated to deterministic stability notions of nonlinear discrete-time dynamics [Khalil, 2002].
In recent years, deep neural networks have been extensively studied from the viewpoint of dynamical systems [Chen et al., 2018, Raissi et al., 2017, Ciccone et al., 2018], allowing for the application of stability analysis methods to DNNs. For instance, Manek and Kolter [2019] proposed neural Lyapunov function to stabilize learned neural dynamics of autonomous system. Haber and Ruthotto [2017] interpret residual connections in neural networks as Euler discretization of ODEs and provide stability guarantees of ResNet. Goel and Klivans [2017] makes connections between eigenvalue decay and learnability of neural networks. Engelken et al. [2020], Vogt et al. [2020] studies the Lyapunov spectrum to of the input-output Jacobian of recurrent neural networks (RNNs) to assess RNN’s stability. In this work we leverage advances in the spectral analysis of deep neural networks and apply them to derive stochastic stability guarantees for DMMs.
Besides analytical methods, many authors introduced provably stable neural architectures [Haber et al., 2019, Greydanus et al., 2019, Cranmer et al., 2020] or stability constraints [John et al., 2017]. Another popular strategy is to employ stabilizing regularizations. This can be achieved by minimizing eigenvalues of squared weights [Ludwig et al., 2014], using symplectic weights [Haber and Ruthotto, 2017], orthogonal parametrizations [Mhammedi et al., 2017], Perron-Frobenius theorem [Tuor et al., 2020], Gershgorin dics theorem [Lechner et al., 2020], or via singular value decomposition (SVD) [Zhang et al., 2018]. In this paper we leverage different weight factorization to empirically validate the proposed theoretical guarantees on stochastic stability of DMMs.

heading: 3 Methodology
text: This section presents stochastic stability analysis method for deep Markov models (DMM). First, we demonstrate the equivalence of DNNs with pointwise affine (PWA) functions. Next, we recall the definition of DMM with transition probabilities modeled by deep neural networks (DNNs). We introduce definitions of stochastic stability and show how can we leverage deterministic stability analysis in the probabilistic context. Finally, we will leverage the equivalence of DNNs with PWA maps to pose sufficient stability stochastic conditions for DMMs based on contraction of PWA maps.

heading: 3.1 Equivalence of Deep Neural Networks with Pointwise Affine Maps
text: Let us consider deep neural network (DNN) ψ : Rn → Rm parametrized by θψ = {Aψ0 , . . .A ψ L ,b0, . . . ,bL} with hidden layers 1 ≤ l ≤ L with bias given as follows:
ψθψ (x) = A ψ Lh ψ L + bL (1a)
hψl = v(A ψ l−1h ψ l−1 + bl−1) (1b)
with hψ0 = x, and v : Rnz → Rnz representing element-wise application of an activation function to vector elements such that v(z) := [v(z1) . . .v(znz )] ᵀ.
Lemma 1. For a multi-layer feedforward neural networkψθψ (3.2) with arbitrary activation function v, there exists an equivalent pointwise affine map (PWA) parametrized by x which satisfies:
ψθψ (x) = Aψ(x)x + bψ(x). (2)
Where Aψ(x) is a parameter varying matrix given as:
Aψ(x)x = A ψ LΛ ψ zL−1A ψ L−1 . . .Λ ψ z0A ψ 0 x (3)
And bψ(x) is a parameter varying vector, both parametrized by input vector x given by following recurrent formula:
bψ(x) = b ψ L (4)
bψl := A ψ i Λ ψ zl−1 bψl−1 + A ψ i σl−1(0) + bl, l ∈ N L 1 (5)
with bψ0 = b0, and i representing index of the network layer. Here Λ ψ zl represents parameter varying diagonal matrix of activation patterns defined as:
σ(z) =  σ(z1)−σ(0) z1 . . .
σ(zn)−σ(0) zn
 z + σ(0)... σ(0)  = Λψz z + σ(0) (6) Proof. First lets observe the following:
σ(z) = σ(z1)... σ(zn)  =  z1(σ(z1)−σ(0)+σ(0)) z1 ...
zn(σ(zn)−σ(0)+σ(0)) zn
 =  σ(z1)−σ(0) z1 . . .
σ(zn)−σ(0) zn
 z + σ(0)... σ(0)  (7)
Remember σ(0) − σ(0) = 0, and zizi = 1 are identity elements of addition and multiplication, respectively. Thus (7) demonstrates the equivalence σ(z) = Λψz z +σ(0) as given in (6). Then if we let zl = A ψ l xl + bl, we can represent a neural network layer in a parameter varying affine form:
σl(A ψ l xl + bl) = Λ ψ zl (Aψl xl + bl) + σ(0) = Λ ψ zl Aψl xl + Λ ψ zl bl + σl(0) (8)
Now for simplicity of exposition lets assume only activations with trivial null space, i.e. σ(0) = 0. Thus σ(z) = Λψz z. By composition, a DNN ψθψ (x) can now be formulated as a parameter-varying affine map Aψ(x)x + bψ(x), parametrized by input x
ψθψ (x) := Aψ(x)x + bψ(x) =
AψLΛ ψ zL−1(A ψ L−1(. . .Λ ψ z1(A ψ 1 Λ ψ z0(A ψ 0 x + b0) + b1) . . .) + bL−1)x + bL
Aψ(x)x = A ψ LΛ ψ zL−1A ψ L−1 . . .Λ ψ z0A ψ 0 x
bψ(x) = A ψ L . . .A ψ 2 Λ ψ z1A ψ 1 Λ ψ z0b0 + A ψ L . . .A ψ 2 Λ ψ z1A ψ 1 σ0(0)+ AψL . . .A ψ 2 Λ ψ z1b1 + A ψ L . . .A ψ 2 σ1(0) + . . .+ A ψ LΛ ψ zL−1bL−1 + A ψ LσL−1(0) + bL
(9)
Hence, each input feature vector x generates a unique affine map Aψ(x)x + bψ(x) of the DNN ψθψ (x). Thus proving the equivalence of DNN map (3.2) with the form (2). The case with σ(0) 6= 0 can be derived following the same algebraic operations as as above.

heading: 3.2 Deep Markov Models
text: We consider a dynamical system with latent state variables xt ∈ Rn, and the observed variables yt ∈ Rm. The transition from xt to the next time step xt+1, and the outputs yt are modeled by probabilistic transitions. Over a horizon of T time steps with a step size ∆t, we assume the Markov property to embed structural independence conditions in the dynamic state evolution, i.e.,
xt+1 ⊥ x0:t−1 | xt, (10)
Thus having latent state dynamics characterized by one-time-step conditional distribution P (xt+1|xt). The joint distribution over the latent states and the observations is given by,
P (x0:T ,y0:T ) = P (x0)P (y0|x0) T−1∏ t=0 P (xt+1|xt)P (yt|xt). (11)
More explicitly, we consider the following probabilistic transition and the emission maps:
xt+1 ∼ N (Kα(xt,∆t), Lβ(xt,∆t)) (Transition) (12a) yt ∼M(Fκ(xt)) (Emission) (12b)
with the initial condition xt=0 = x0. Here, N and M denote the probability distributions. For the transition mapping, N denotes a Gaussian distribution with the mean vector Kα(xt,∆t), and covariances Lβ(xt,∆t). The distributionM can be arbitrary with its distribution characterized by the map Fκ(xt). In this paper, we are interested in the stability characterizations of the latent state dynamics given by transition map (12a), thereby assuming full state observability.
We are interested in expressing the dynamics of complex systems using (12), therefore it is suitable to expand the expressivity of our model by parametrizing the conditional distribution P (xt+1|xt) by deep neural networks (DNNs) given as,
Kα(xt,∆t) = fθf (xt), (13) vec(Lβ(xt,∆t)) = gθg(xt), (14)
where f : Rn → Rn and g : Rn → Rn2 are two deep neural networks parametrized by θf and θg, respectively. And vec(·) denotes standard vectorization operation. Therefore, the probabilistic transition dynamics (12a) can be characterized by analysing the stability and boundedness of deep neural networks f and g.

heading: 3.3 Stability of Deep Markov Models
text: To this end, we bring forth a few stability notions in the context of stochastic state transitions. In stochastic dynamics and control literature [Khasminskii, 2012], various different notions of stability for stochastic state transitions, such as mean-square stability, almost-sure stability and stability via convergence in probability, have been discussed. In this article, since we are interested in the latent state trajectories of the dynamic systems, we consider the mean-square stability as defined in [Willems and Willems, 1976, Nandanoori et al., 2018]. The dynamic system is said to be mean-square stable if the first and the second moment converge over time. Definition 2. The stochastic process xt ∈ Rn is mean-square stable (MSS) if and only if there exists µ ∈ Rn,Σ ∈ Rn×n, such that limt→∞ E(xt) = µ, and limt→∞ E(xtxTt ) = Σ.
The MSS condition from Definition 2 requires the dynamics (13) to have stable equilibrium x̄e = µ, where x̄ denotes the mean state vector. The definition also requires the covariance to converge. To express the dynamic behavior of complex system, the second moment convergence criterion can be relaxed by only requiring it to be norm bounded in order to ensure stochastic stability given as,
‖gθg(xt)‖p < K, K > 0, ∀t. (15) Here ‖ · ‖p denotes any appropriate vector norm, e.g., L2−norm. The bound on the covariances will depend on the extent of stochasticity that the dynamic system encounters in an uncertain environment. However, in the later parts, we consider the convergence scenario for g as in Definition 2, rather than merely boundedness, in prescribing the sufficient conditions for the MSS-type stability.
Let us first analyze the mean dynamics characterized by (13), and its equilibrium x̄e = µ which satisfies the stationarity condition fθf (x̄e) = x̄e. We have the mean state vector x̄t evolving under the following dynamics:
x̄t+1 = fθf (x̄t). (16) (16) allows us to analyze the asymptotic stability of the DMM mean dynamics. Now in the main result of this section we leverage the fact that the dynamic characteristics of the deep neural networks fθf , gθg around a point x̄t can be evaluated by obtaining their exact pointwise affine forms (PWA) (2). Based on this equivalence we formulate Theorem 3 and Corollary 5 as follows with sufficient conditions for the stability of deep Markov models. Theorem 3. The deep Markov model (12) which is parametrized by deep neural networks (13)-(14) remains globally stable in the mean-square sense if the following holds. The mean neural network fθf (x) is a contractive map for any x in the domain of fθf (x). The variance network gθg(x) is bounded for any x in the domain of gθg(x). Or more formally:
‖Af (x)‖p < 1 (17a)
||Ag(x)||p + ||bg(x)||p ‖x‖p < 1, (17b)
∀x ∈ Domain(fθf (x),gθg(x)). (17c)
Proof. First we prove the sufficiency of the contraction condition of the mean dynamics (17a). We base the proof on the equivalence of multi-layer neural networks with pointwise affine maps (2). An affine map is a contraction if the 2-norm of its linear part is bounded below one, i.e. ||A||2 < 1. Thus it follows that the condition (17a) and equivalence (2) imply a contractive mean neural network fθf (x). The sufficiency of the contraction condition on mean square stable (MSS) equilibrium in the sense of Definition 2 follows directly from the Banach fixed-point theorem, which states that every contractive map converges towards single point equilibirum. Hence condition (17a) implies convergent mean transition dynamics:
µ = fθf (µ) = limt→∞ fθf (x̄t) (18)
Now we show the sufficiency of (17b) to guarantee the boundedness of the covariance matrix elements (15) by bounding the p-norm of the covariance neural network ‖gθg(x)‖p. Please note that using the form (2) gives us ‖gθg(x)‖p = ‖Ag(x)x + bg(x)‖p yielding following inequalities:
‖gθg(x)‖p ≤ ‖Ag(x)x‖p + ‖bg(x)‖p, (19a) ||gθg(x)||p ‖x‖p ≤ ‖Ag(x)‖p + ||bg(x)||p ‖x‖p . (19b)
We show that (19b) gives in fact a local Lipschitz constant of the variance network gθg(x). We exploit the point-wise affine nature of the neural network’s form (2) and the fact that the norm of a linear operator A is equivalent to its minimal Lipschitz constant KAmin = ||A||p [Huster et al., 2018]. Thus we can compute the local Lipschitz constants of a neural network gθg(x) as:
Kg(x) = ||Ag(x)||p + ||bg(x)||p ‖x‖p . (20)
Applying the upper bound (17b) on the local Lipschitz constant (20) guarantees the contraction of the variance neural network fθg(x) towards a fixed steady state Σ.
Remark 4. To guarantee stochastic stability, the condition (17b) can be relaxed as given in (15) to bounded second moment (15) with maxxKg(x) < K, where K > 0. Corollary 5. The deep Markov model (12) which is parametrized by deep neural networks (13)-(14) remains globally stable in the mean-square sense if the following holds: All weights Afi of the mean network are fθf contractive maps. All activation scaling matrices Λ f zi of the mean network are non-expanding. Norms of all weights Agj and activation scaling matrices Λ g zj of the variance network gθg are upper bounded by 1. Or more formally:
‖Afi‖p < 1, ||Λfzi ||p ≤ 1 i ∈ N Lf 1 , (21a)
‖Agj ‖p < 1, ||Λ g zj ||p ≤ 1, j ∈ N Lg 1 , (21b)
∀x ∈ Domain(fθf (x),gθg(x)). (21c)
Proof. First we show the sufficiency ‖Ai‖p < 1 of contractive weights and non-expanding activation scaling matrices ‖Λzi‖ ≤ 1 to guarantee the contractivity of arbitrary deep neural networks. Assuming general non-square weights Ai ∈ Rni×mi we use the submultiplicativity of the induced p-norms to upper bound the norm of a products of m matrices given as:
‖A1 . . .Am‖p ≤ ‖A1‖p . . . ‖Am‖p (22) Now by applying (22) to the linear parts (3) of the mean neural network fθf in the pointwise affine form (2) with ‖Afi‖p < 1, ∀i ∈ NL0 , ‖Λfzj‖p ≤ 1, ∀j ∈ N L 1 , it yields ‖Af (x)‖p < 1 over the entire domain of fθf (x), thus with p = 2 satisfying the contraction condition ||Af (x)||2 < 1 for affine maps. The submultiplicativity (22) naturally applies also to the variance network gθg(x) thus implying the contraction towards a fixed point given the conditions (21).
Remark 6. To guarantee stochastic stability, we can relax the upper bound of the second moment as K = ∏L i c
AcΛ, where cA > 0, and cΛ > 0 represent the relaxed upper bounds of the operator norms in condition (21b). Thus satisfying the relaxed boundedness condition on the variance via (15).
Assuming the contraction conditions (17) or (21) hold and the mean neural network fθf has zero bias, then the DMM’s mean network fθf (x) is equivalent with stable parameter varying linear map (3) with equilibrium in the origin, i.e. x̄ = 0. In the case with non-zero bias in fθf , the corresponding PWA map (2) has non-zero equilibrium x̄ 6= 0. Both conditions (17) or (21) are sufficient for a convergence of a DMM (12) to a stable equilibrium x̄. However, they do not provide bounds of the admissible values of the equilibrium x̄. The corresponding equilibrium bounds are provided in the supplementary material.

heading: 3.4 Design of Stable Deep Markov Models
text: In this section, we provide a set of practical design methods for provably stable DMM (12). Based on the Corollary 5 the use of contractive activation functions together with contractive weights for both mean and variance network will guarantee the stability of DMM by design. In particular, the conditions (21) on bounded norm of transitions’ activation scaling matrices ||Λfzi || ≤ 1, ||Λ g zj || ≤ 1 implies Lipschitz continuous activation functions with constant K ≤ 1. Conveniently, this condition is satisfied for many popular activation functions such as ReLU, LeakyReLU, or tanh. The contractivity conditions (21) on weight matrices ‖Afi‖p < 1, ‖A g j ‖p < 1, respectively, can be enforced by employing various matrix factorizations proposed in the deep neural network litertature. Examples include singular value decomposition (SVD) [Zhang et al., 2018], Perron-Frobenius (PF) [Tuor et al., 2020], and Gershgorin discs (GD) [Lechner et al., 2020] factorizations given below.
PF weights: This factorization applies Perron-Frobenius theorem for constraining the dominant eigenvalue of the square nonnegative matrices. Based on this theorem, we can construct the weight matrix A with bounded eigenvalues as follows:
M = λmax − (λmax − λmin)g(M′) (23a)
Ai,j = exp(A′ij)∑nx k=1 exp(A′ik) Mi,j (23b)
here M represents the damping factor parameterized by the matrix M′ ∈ Rnx×nx , while A′ ∈ Rnx×nx represents the second parameter matrix encoding the stable weights A. The lower and upper bound of the dominant eigenvalue are given by λmin and λmax, respectively.
SVD weights: Inspired by singular value decomposition (SVD), this method decomposes a possibly non-square weight matrix A = UΣV into two unitary matrices U and V, and a diagonal matrix Σ with singular values on its diagonal. The orthogonality of U and V is enforced via penalties:
Lreg = ||I−UU>||2 + ||I−U>U||2 + ||I−VV>||2 + ||I−V>V||2 (24) An alternative approach to penalties introduced in Zhang et al. [2018] is to use Householder reflectors to represent unitary matrices U and V. The constraints λmin and λmax on the singular values λ can be implemented by clamping and scaling given as:
Σ = diag(λmax − (λmax − λmin) · σ(λ)) (25)
GD weights: This method supporting square matrices leverages the Gershgorin dics theorem [Varga, 2004]. It says that all eigenvalues λi of the weight A can be bounded in the complex plane with center λ and radius r given by the formula:
A = diag ( r s1 , ..., rsn ) M + diag (λ, ..., λ) (26)
Where M ∈ Rn×n with mi,j ∼ U(0, 1), except mi,i = 0 is learnable parameter matrix. While diagonal matrices diag ( r s1 , ..., rsn ) , and diag (λ, ..., λ) represent radii and centers of the bounded
eigenvalues, where sj = ∑ i 6=jmi,j .
Parametric stability constraints: The disadvantage of enforcing the global stability conditions as given via Corollary 5 is their negative effect on the expressivity of the DMM, resulting in dynamics with a single point or line attractors. This will effectively prevent the DMM from expressing more complex attractors such limit cycles or chaotic attractors. As a more expressive alternative we introduce the use of parameter varying bounds in the conditions (17), and (21), such as:
p(x) < ‖Af (x)‖p < p(x) (27)
Where p(x) : Rnx → R, and p(x) : Rnx → R are scalar valued functions parametrizing lower and upper bounds of operator norm of the DMM’s mean transition dynamics. Similar parametric constraints can be applied to the variance bounds in (17), or weight norm constraints in (21). This approach allows us to control the contractivity of the DMMs depending on the position in the state space. This allows us to increase the expressivity of the DMM, e.g., by partitioning the state space into constrained and unconstrained regions resulting in DMM with hybrid or switching dynamics. In particular, we could divide the state space to outer contractive regions (where conditions (17) hold) and inner relaxed regions allowing for more complex trajectories to emerge. This parametrization will effectively generate non-empty attractor set in which it is possible to learn arbitrary attractor shape. The proposed state space partitioning method is inspired by the Bendinxon-Dulac criteria on periodic solutions of differential equations [McCluskey and Muldowney, 1998].

heading: 4 Numerical Case Studies
text: In this section we empirically validate the conditions given in Theorem 3 and Corollary 5 by investigating the dynamics of DMM’s transition maps (12) whose mean fθf (x) and variance gθg(x) are parametrized by neural networks with different spectral distributions of their weights and activation scaling matrices (6). We apply spectral analysis to the PWA forms (2) of neural networks modeling the mean and variance maps to obtain the corresponding spectra of DMMs. We performed the experiments using the probabilistic programming language Pyro [Bingham et al., 2019].

heading: 4.1 Design of the Experiments
text: Since the stability of DMM (12) depends on the transition dynamcis, in all of the case studies we consider a fully observable model with identity as an emission map. We parametrize the mean and variances of the transition map fθf (x) (13) and gθg(x) (14) by feedforward neural networks. Given the mean neural network fθf (x) we generate a set of different transition dynamcis by changing activation functions v(x) ∈ {ReLU,Tanh,Sigmoid,SELU,Softplus}, layer depth L ∈ {1, 2, 4, 8}, and presence of bias b ∈ {True,False}. For the variance network gθg(x) we use ReLU activations. For both, mean and variance networks we initialize their weights with desired spectral properties via design methods described in Section 3.4. In particular we use SVD, PF, and GD factorizations to bound the weight’s singular values in a prescribed range. We generate three categories of weights Ai: stable with operator norm strictly below one ||Ai||p < 1, marginally stable with norm close one ||Ai||p ≈ 1, and unstable with norm larger than one ||Ai||p > 1.

heading: 4.2 Stability Analysis of Deep Markov Models
text: In order to provide intuitive visualisations of the dynamics in the phase space, in this section, we focus on two dimensional system. Fig. 1 visualizes the phase portraits of randomly generated DMM’s probabilistic transition maps of the mean fθf (x) and variance fθg(x) networks with constrained operator norms enforced using PF weights from Section 3.4. Figures in the first row demonstrate that DMMs with asymptotically stable mean transition dynamics ||Af ||p < 1 with bounded variances ‖Ag‖p < K, K > 0 generate stable single point attractors. Hence they validate the sufficient conditions of Theorem 3. Figures in the second row display dynamics of DMM with marginally stable mean ||Af ||p ≈ 1. Due to the non-dissipativeness of the mean transition dynamics, the trajectories converge to a line attractor only if the variance is a converging map ‖Ag‖p < 1, thus having a dissipative second moment. In case of marginally stable variance, ‖Ag‖p ≈ 1 the energy conserving nature of the mean and variance together generate random walk type trajectories along the direction of the mean’s line attractor. For the cases with unstable variance ‖Ag‖p > 1, the overall dynamics behaves close to a Brownian motion with a degree of randomness, which is positively
correlated with the variance network’s operator norm. Figures in the third row show diverging dynamics of DMM with unstable mean ||Af ||p > 1. With converging variance ‖Ag‖p < 1 the diverging stochastic trajectories stay close to the mean direction. While, for both marginal and unstable variances the stochastic trajectories diverge in all directions.
4.3 Effect of Biases and Depth on the Stability of Deep Markov Models
In Fig. 2, we experiment with biases and depths of mean fθf (x) and variance gθg(x) networks.
Effect of biases: In the left panels of Fig. 2, we demonstrate the dynamics of DMM with tanh activations and SVD factorized weights resulting in asymptotically stable transition maps, thus ||Af (x)|| < 1, and ||Ag(x)|| < 1. Fig. 2(a) shows the scenario without any bias whereas Fig. 2(d) shows the scenario where both fθf (x),gθg(x) have bias terms. It demonstrates that the general contractive nature of the stable behavior, as given via conditions (17), does not change with addition of biases. Instead, the biases shift the region of attraction by generating non-zero equilibrium points. This shift is correlated with absolute value of the aggregate bias term of the PWA form (2). For the norm bounds on the equilibria of stable DMM see supplementary material.
Effect of depth: The right panels of Fig. 2 demonstrate the dynamics of DMMs with increasing number of layers using ReLU activations and SVD weights close to marginal stability 0.99 < ||Ai(x)|| < 1. It can be seen that with increase in the number of layers, the convergence of trajectories toward origin becomes less uncertain. In this case, a larger number of mildly contractive layers results in stabilizing behavior, demonstrating the effect of the norm submultiplicativity on the operator norm of the mean transition dynamics (22). However, one needs to be careful about the delicate balance between stability, and ability to efficiently train the parameters of DMM with gradient-based optimization. With increasing depth, the norm product of the contractive layers (22) will eventually result in a network with a very small operator norm thus causing the vanishing gradient problem. Analogously, the exploding gradient problem will occur for DMM parametrized with very deep neural networks with non-contractive layers, i.e. ||Ai(x)|| > 1 The use of parametric stability constraints (27) as a function of depth could be an efficient strategy for avoiding the vanishing and exploding gradients by keeping the overall dynamics norm bounded.

heading: 4.4 Deep Markov Models with Parametrized Stability Constraints
text: In Fig. 3 we demonstrate the use of parametrized stability constraints (27) in the design of stable DMMs without compromising the expressivity as it is in the case of restrictive single point attractors
enforced via (17) and (21). In particular, we design two DMMs with randomly generated weights with three phase space regions with different mean transition dynamics, (i) an inner expanding region ||Af (x)|| > 1, xi ∈ R1, (ii) a middle marginal region ||Af (xi)|| ≈ 1, xi ∈ R2, and (iii) an outer contractive region ||Af (xi)|| < 1, xi ∈ R3. Where the regions are given asR1 = {x|0 ≤ ||x||2 < 20},R2 = {x|20 ≤ ||x||2 < 40}, andR3 = {x|40 ≤ ||x||2}, respectively. The variance dynamics in both cases is kept being contractive ||Ag(xi)|| < 1. From Fig. 3 it is apparent that the overall dynamics of the DMMs with parametrized constraints (27) is able to generate stochastic periodic behavior while remaining bounded within prescribed region of attraction, thus providing high degree of expressivity while being provable stable. On the left (Fig. 3 (a) and (c)) we show phase plots, and on the right (Fig. 3 (b) and (d)) corresponding time series trajectories. As a potential extension, we envision learning the constraints bounds p(x), and p(x) in (27) using penalty methods.

heading: 5 Conclusion
text: In this paper, we introduce a new stability analysis method for deep Markov models (DMMs). As the main result, we provide sufficient conditions for the stochastic stability and introduce a set of practical methods for designing provably stable DMMs. In particular, we discuss the use of contractive weight matrices factorizations and stability conditions for the activation functions. Furthermore, we propose using novel parametric stability constraints allowing expression of more complex stochastic dynamics while remaining contractive towards non-empty region of attraction. The proposed theory is supported by numerical experiments, with design guidelines considering weight factorizations, choices of activation functions, network depth, or use of the bias terms for guaranteed stability. In future work, we aim to derive the stability guarantees for a broader family of probability distributions modeled by normalizing flows. We also aim to expand the theory to partially observable Markov decision processed (POMDP) to derive closed-loop stability guarantees in the context of deep reinforcement learning.
Acknowledgments and Disclosure of Funding
We acknowledge our colleagues Aaron Tuor, Mia Skomski, Soumya Vasisht, and Draguna Vrabie for their contributions to related topics that served as a base for developing the method presented in this paper. We would like to thank Craig Baker and David Rolnick for fruitful discussions that helped improve the technical quality of the presented ideas. Also, we want to thank our anonymous reviewers for their constructive feedback and suggestions.
This research was supported by the U.S. Department of Energy, through the Office of Advanced Scientific Computing Research’s “Data-Driven Decision Control for Complex Systems (DnC2S)” project. Pacific Northwest National Laboratory is operated by Battelle Memorial Institute for the U.S. Department of Energy under Contract No. DE-AC05-76RL01830. Oak Ridge National Laboratory is operated by UT-Battelle LLC for the U.S. Department of Energy under contract number DE-AC0500OR22725.

heading: 6 Limitations and Broader Impact
text: Stability is the major concern of many safety-critical systems. The significance of the presented work lies in the proposed stability conditions, analysis, and design methods allowing the generation of provably stable DMMs. Furthermore, the presented methods for the design of stable DMMs could be of large significance in the context of stochastic control systems with stability guarantees. Thus, the authors believe that besides academic relevance, the methods presented in this paper have the potential for practical impact in many real-world applications such as unmanned autonomous vehicles, robotics, or process control applications.
The authors are aware of the limitations of the presented numerical case studies focusing on smallscale DMMs with two-dimensional state space. This choice was made for the sake of the visualizations of the state space trajectories allowing us to provide intuitive examples of the presented theoretical results. As part of the future work, the authors plan to use the proposed stability analysis and design methods for DMMs on real-world datasets with higher dimensional state space.
The presented work falls into the basic research category. As such, authors are not aware of any potential direct negative societal impact of the proposed work. On the contrary, authors of this paper believe that the presented theory is a minor contribution towards general knowledge, which accumulation has been historically proved to inherently benefit all humanity.

Review summaries:
The paper discusses stable deep Markov models (DMM), which can be used to model stochastic dynamical systems, restricting to the case where state transition probabilities are given by normal distributions with mean and covariance matrices parameterized by two deep neural nets. One of the main contributions of this paper is to study sufficient conditions under which such DMMs are stable (in the mean-square sense), the second is to provide a list of practical (re)parameterization of the DNNs that satisfy the sufficient conditions.  The paper concludes with a series of small scale numerical experiments aimed at showing the possible behaviours of the class of proposed models. 

This paper analyzes deep Markov models from a theoretical perspective. More specifically, it derives conditions on the construction/parameters of models which allow for a specific notion of stochastic stability. This is accomplished by first proving that a multilayer perceptron can be equivalently represented as a pointwise affine map (PWA). It is then proven that achieving specific properties on the norms of these PWAs leads to global stability, and that bounding the norms of the layers in neural networks leads to achieving these desired properties within the corresponding PWAs. There is a discussion on how one can achieve these properties when training a neural network, followed by empirical analysis demonstrating the connections between the properties of DMMs and their stability. 

This paper analyses stochastic stability of Deep Markov Models (DMMs) based on the contraction of probabilistic maps modeled by Deep Neural Networks (DNNs). To this end, DMMs are modeled using latent space representations of a dynamical system, and Gaussian distribution is used in modeling the transition and emission dynamics through DNN representation. With this modeling assumption, the authors perform stochastic stability analysis of DMM focused on bounding the norm of the model parameters of DNN. They state sufficient conditions for achieving stability and outline techniques for doing so. Finally, their claim is substantiated through small-scale experimentation of DMMs.