Title: LINEARIZING VISUAL PROCESSES WITH DEEP GEN-ERATIVE MODELS

Abstract: This work studies the problem of modeling non-linear visual processes by leveraging deep generative architectures for learning linear, Gaussian models of observed sequences. We propose a joint learning framework, combining a multivariate autoregressive model and deep convolutional generative networks. After justification of theoretical assumptions of linearization, we propose an architecture that allows Variational Autoencoders and Generative Adversarial Networks to simultaneously learn the non-linear observation as well as the linear state-transition model from a sequence of observed frames. Finally, we demonstrate our approach on conceptual toy examples and dynamic textures.

1. INTRODUCTION:

While classification of image and video with Convolutional Neural Networks (CNN) is becoming an established practice, unsupervised learning and generative modeling remain to be challenging problems in deep learning. A generative model of a visual process enables the possibility of generating sequences of video frames such that the appearance as well as the dynamics approximately resemble the original training process without copying it. This procedure is typically referred to as video generation (Vondrick & Torralba (2017) ; De Souza et al. (2017) ) or video synthesis (Liu et al. (2017b) ). More technically, this means that in addition to a suitable probability model for the individual frames, a probabilistic description for the frame-to-frame transition is also necessary. Analysis and reproduction of visual processes simplifies considerably, if this transition can be described as a multivariate autoregressive (MAR) model, i.e., as a combination of linear transformations and Gaussian noise. For instance, linear transformations are easily invertible and by means of spectral analysis, it can be studied how such a process behaves in the long term.

Realistically, most frame transitions in real-world visual processes unlikely are linear functions. Nevertheless, unsupervised learning has come up with many approaches to fit MAR models to realworld processes, for instance by using linear low-rank approximations, as proposed by Doretto et al. (2003) , or sparse approximations of the frames, as proposed by Wei et al. (2017) , or applying the kernel trick to them (Chan & Vasconcelos (2007) ).

The success of Generative Adversarial Networks (GAN) introduced by Goodfellow et al. (2014) and Variational Autoencoders (VAE) introduced by Kingma & Welling (2013) has led to an increased interest in deep generative learning and it seems natural to apply such techniques to sequential processes. We approach this idea from the perspective of linearization, in order to keep the model as simple as possible. In an analogous way as physicists transforming non-linear differential equations into linear ones by means of an appropriate change of variables, our approach is to learn latent representations of visual processes, such that the latent state-to-state transition can be described by an MAR model. To this end, we jointly learn a non-linear observation and a linear state transition function by introducing a dynamic layer that can be used in conjunction with deep generative architectures such as GANs and VAEs.

NOTATION:

In this paper, letters that appear both in roman and italicized type refer to random variables and realizations thereof, respectively. If not stated otherwise, expected values of more than one random variables assume statistical independence.

2. RELATED WORK:

The deep predictive coding network (Chalasani & Principe (2013)) , can be considered one of the earliest approaches of combining MAR models with deep learning. More recently, such a combination was studied by Liu et al. (2018) . The work by Goroshin et al. (2015) deals with linearizing transformations under uncertainty via neural networks. It resembles our work in that we also focus on representation learning rather than on a particular type of process. However, these works do not employ VAEs or GANs.

By contrast, Watter et al. (2015) combine Linear Dynamic Systems (LDSs) with VAEs. However, the focus of their work is on control rather than on synthesis. Furthermore, their model is locally linear and the transition distribution is modeled in the variational bound, whereas we model it as a separate layer. This also is the main difference to the work of Krishnan et al. (2015) , in which VAEs are used as Kalman Filters. The work of Johnson et al. (2016) combines VAEs with linear dynamic models for forecasting images in video sequences. Mathematically, it is a well-thought out approach, but proposes a training objective that is considerably more complex than the one proposed in this work. GANs in combination with MAR models have been studied in experiments of Han et al. (2017) , but the MAR model was learned separately from the GAN.

Theoretical groundwork regarding learned visual transformations has been done by Cohen & Welling (2014; 2016) ; Hinton et al. (2011) and Memisevic (2013) . In particular, the dynamic layer proposed in this work bears resemblance to techniques employed in image-to-image translation (Liu et al. (2017a) ; Zhu et al. (2017) ).

The synthesis of video dynamics by means of neural networks has been discussed, among others by Vondrick et al. (2016) ; Xue et al. (2016) ; Srivastava et al. (2015) and Xie et al. (2017) . We would like to emphasize the difference between video synthesis which is discussed in our work and the prediction of video frames as studied, for instance, by Mathieu et al. (2015) . While the former refers to the problem of finding a probabilistic model for the spatial and temporal behavior of a visual process, the latter refers to finding a deterministic mapping from a set of previous frames to one or several future frames to come. As a consequence, a video synthesis model needs to take care of long-term behavior. Additionally, the probabilistic nature of video synthesis makes it considerably harder to evaluate the generated frames, since a unique ground truth cannot be provided, ruling out classical quantitative quality measures such as mean squared error or structured similarity.

Finally, the core contribution of this work is a combination of neural networks with Markov processes. This has been the subject of many works in the recent past. For a broad overview of results in this field, the reader is referred to Chapter 20 of the book Deep Learning by Goodfellow et al. (2016) .

3. VISUAL PROCESSES:



3.1. DYNAMIC SYSTEMS:

The Dynamic Texture (DT) model by Doretto et al. (2003) has popularized LDSs in the modeling of visual processes. Typically, an LDS is of the following form

EQUATION

where h t ∈ R n is the low dimensional state space variable at time t, A ∈ R n×n the state transition matrix, y t ∈ R d the observation at time t and C ∈ R d×n the observation matrix. The vector ȳ ∈ R d represents a constant offset in the observation space. The input term v t is modeled as i.i.d. zero-mean Gaussian noise, and is independent of h t .

Real-world visual processes are often highly non-linear and non-Gaussian, hence are considerably harder to work with. In the following, we define a more general dynamic model that is applicable to a broad class of visual processes. To keep the problem tractable, we make the assumption of two properties that were already implicitly assumed in the classical DT model.

Assumption 1. The visual processes of interest are stationary first-order Markov processes.

The stationarity property guarantees that synthesis of new sequences will not diverge and the Markov property facilitates the statistical inference. However, we believe that the method presented in this work can be generalized to Markov processes of higher orders. We refer to Appendix A.2 for details.

Consider the following generic dynamic equation of observations y in an n-dimensional manifold M ⊂ R d . y t+1 = ψ vt (ϕ(y t )).

(2) Again, y t ∈ R d denotes the observation at time t. The self map ϕ : M → M models the predicted frame transition and describes the deterministic part of the dynamics. In this work, it is assumed to be differentiable. Furthermore, the function ψ vt : M → M describes a displacement by v t in the in the tangent space T yt M of M at y t , followed by a retraction onto M. The displacement is sampled from i.i.d. zero-mean Gaussian noise. Eq. ( 1) is a special case of Eq. ( 2), where ϕ(y) = CAC + (y -ȳ) + ȳ and ψ v (y) = y + Cv.

The model Eq. ( 2) describes a much broader class of visual processes than Eq. ( 1). However, unlike model Eq. ( 2), the state transition in the linear model Eq. ( 1) enables straightforward prediction, generation, and analysis of observations. It is thus of great interest to find a model that linearizes real-world visual processes, so that in some latent state space representation, the state transition admits the MAR model of the first line of Eq. ( 1). Specifically, this work focuses on the following non-linear dynamic system model, i.e., a linear state transition and a non-linear observation mapping

EQUATION

where n < d and C : R n → M ⊂ R d is differentiable. Transforming a non-linear model to this form makes video synthesis as easy as sampling from autoregressive noise.

It is worth noticing that the model in Eq. ( 3) is not unique with respect to changes of basis in the state space (Doretto et al. (2003) ; Afsari & Vidal (2013) ). However, if C is implemented via a neural network, we can ensure that it accounts for a possible change of basis, e.g. by adding a linear layer to the input of the network. We can thus make the following assumption on the latent samples h t without loss of generality. Assumption 2. The latent samples h t abide a standard normal distribution, i.e., h t ∼ N ( • ; 0, I). Remark 1. If the state transition matrix A is given, and the process is stationary, i.e., p(h t ) = p(h t+1 ) for all t, then Assumption 2 essentially identifies the process noise model. Namely, we have

EQUATION

In order to make sure that the latent states h t remain standard Gaussian in sequential synthesis scenarios, i.e., E ht+1 h t+1 h t+1 = I, we just need to ensure that the process noise is zero-mean and has the covariance matrix I -AA . ♦

Before we propose an algorithm to jointly learn C and A by means of deep generative architecture, we further characterize the problem in order to investigate its feasibility and motivate our approach. For a model Eq. ( 3) to substitute Eq. ( 2), it needs to transform the transition ϕ to a multiplication by a matrix and the perturbation ψ vt to a superposition by zero-mean Gaussian noise. The latter is generally easy to achieve, if the perturbation is sufficiently small and C is a diffeomorphism. In the remainder of the section, we thus focus on linearization of ϕ.

3.2. LOCAL LINEARIZATION:

A common linearization method in control theory is to approximate the dynamic system equation by a first-order Taylor polynomial around an equilibrium point (Perko (2013) ). Clearly, for this to be possible, such a point needs to exist. We thus propose the following assumption on the transition function ϕ of Eq. ( 2).

Assumption 3. The differentiable self map ϕ : M → M has at least one fixed point y * .

However, such an approach might fail for real-world processes due to the curse of dimensionality.

A remedy is to employ an representation function Γ that maps the system observations to a lowerdimensional latent space prior to performing the linearization. We call such a function a local linearizer. Definition 1. Let M ⊂ R d be an n-dimensional manifold, ϕ : R d → R d be differentiable on ϕ(M) = M, and Γ : R d → R n be a diffeomorphism on M. The map Γ is said to be a local linearizer at y * ∈ M of ϕ , if there exists a matrix Φ ∈ R n×n such that the following equality holds true with y ∈ R d .

EQUATION

Generally speaking, local linearization is made possible by moving the fixed point to the origin of a new coordinate system. More precisely, the following proposition holds. Proposition 1. Let M ⊂ R d be an n-dimensional manifold. Furthermore, let ϕ : R d → R d be differentiable on ϕ(M) = M, and Γ : R d → R n be a diffeomorphism on M. If y * ∈ M is a fixed point of ϕ, then the following map,

EQUATION

locally linearizes ϕ, if the rows and columns of the Jacobi matrix J ϕ of ϕ at y * lie in the row space of the Jacobi matrix J Γ of Γ at y * .

Proof. See Appendix A.1

The Jacobi matrix property is only needed for consistency with Definition 1, where the limit y → y * is approached from any possible direction in R d . It can be ignored, if we restrict the analysis entirely to the manifold M. Proposition 1 states that if ϕ is a differentiable self map on M with a fixed point y * ∈ M, and a diffeomorphism from M to R n exists, then there is a representation in which ϕ can be approximated by a linear function for points on M that are not too far away from y * .

Note that even if the Jacobi matrix J ϕ can be sufficiently well estimated from the data, the linear approximation directly in the observation space R d will likely be worse than via reparametrization of a local linearizer Γ that functions as a chart of M. The reason is that predictions directly via J ϕ can not be assumed to remain on M. Figuratively speaking, the aim of the local linearizer is to bend the space spanned by the rows and columns of J ϕ to match the shape of M.

3.3. SEPARATE LEARNING:

Local linearizability is a useful concept to analyze the feasibility of the problem at hand. However, it does not take into account that two local linearizers can have significantly different linearization properties on a global scale. Moreover, it does not provide any instructions on how to find an appropriate representation Γ and matrix Φ. Traditionally, this task is approached by learning Γ and Φ separately while considering the sampled observations of the process globally, as will be described in the following.

To characterize the problem of global linearization, we need to introduce a measure. It is sensible to consider the expected squared distance between the result of a transformation ϕ and its linear approximation. Since it is not possible to have an analytic expression for the distribution on the data manifold M, the latent space is often considered as a heuristic. Because it is a common model assumption for deep generative models and in accordance with Assumption 2, we can assume standard Gaussian distribution for the latent space.

Let Γ : M → R n denote a data representation mapping and ϕ : M → M the transition function to be linearized. Consider the following expression

EQUATION

The denominator 2n is a normalization factor. If Γ(y) has standard Gaussian distribution, then q(Γ, ϕ, Φ) < 1 denotes that the linear prediction with Φ is smaller than the expected distance of two independently drawn samples. A common method of linearizing ϕ, is to first choose a representation Γ that is assumed to have good linearization properties, and then inferring Φ by minimizing Eq. (7) (Doretto et al. (2003) ; Chan & Vasconcelos (2007) ; Han et al. (2017) ). In that case, the solution is given by

EQUATION

if Γ(y) has standard Gaussian distribution.

However, the drawback of this approach is the difficulty to find an appropriate model for Γ that can be assumed to linearize the transformation ϕ. Moreover, separate learning implicitly assumes that a small linearization error in the embedded space will carry over to a small error in the observation space, but such an assumption is hard to justify given the high dimension of the problem. We therefore propose to approach the problem by learning the representation and the linear transition jointly, by approximating the data distribution directly via deep generative models.

4. JOINT LEARNING VIA DEEP GENERATIVE MODELS:



4.1. MARKOV ASSUMPTIONS:

We now turn to the problem of finding a model Eq. ( 3) that approximates Eq. ( 2). In the following, we will provide the preliminaries to lay out a training procedure to infer C and A from observed sequences of visual processes. A sequence of length N of a visual process Y is viewed as a realization of the random variable y N = [y 1 , . . . , y N ]. The according sequence in the latent space is a realization of the random variable

EQUATION

In order to model Y by Eq. (3), we need to make sure that the probability distributions of y N and ỹN coincide for any N ∈ N. By taking into account Assumption 1, this is equivalent to demanding that the joint probability of two succeeding frames coincide.

The joint probability distribution of h t and h t+1 is zero-mean Gaussian and, more specifically,

EQUATION

holds due to Assumption 2. To summarize, we are looking for the function C and a matrix A, such that the random variable

EQUATION

has the same probability distribution as y 2 .

4.2. DEEP GENERATIVE MODELS:

Estimation of probability distributions for high-dimensional data is still an ongoing research topic in deep learning. The problem is typically framed as an approximation of a function f θ , parameterized by a vector θ in a finite-dimensional euclidean space. The purpose of f θ is to map realizations of low-dimensional, standard Gaussian noise to samples that abide the data distribution. Typically, the function f θ is realized by a neural network with trainable weights θ.

Most of the widely employed deep generative models, including the GAN, and the VAE, adopt this approach. What these models differ in, is merely the way the network parameters θ are trained. For this work, we employ the Wasserstein GAN and the VAE to evaluate the proposed method. Due to space constraints, we will not review these techniques here and refer the reader to Arjovsky et al. (2017) for the Wasserstein GAN and to Doersch (2016) for the VAE. For the following sections, it is sufficient to assume that these techniques are capable of approximating a function f θ that maps from a low-dimensional standard Gaussian to a high-dimensional data distribution.

θ = (η, A, B) f θ (•) N (x; 0, I) x B• A• C η (•) C η (•) h 1 h 2 x 1 x 2 C η (•)

s We now focus on how to learn a matrix A and a function C from a finite sequence y N ∈ R d×N of a visual process such that it can be described by Eq. ( 3). To this end, recall that, as discussed in Section 4.1, we only need to consider the joint probability of two succeeding observed frames. The first step is thus to generate a training set of frame pairs from y N as {s 1 , . . . , s N -1 }, where s t = y t y t+1 (12) denotes a vectorized pair of frames. We treat the samples as realizations of the random variable s and are looking for an architecture to learn the function C and a matrix A such that the probability distribution of the random variable

EQUATION

where h 1 , h 2 have the joint distribution described by Eq. ( 10), coincides with the probability distribution of s, i.e.

EQUATION

)

It turns out that we can use a deep generative architecture to accomplish this task. Remember, that a function f θ learned by a deep generative model is capable of transforming samples x of standard Gaussian noise to samples from a high dimensional data distribution. Let θ = (A, B, η) contain the matrices A, B ∈ R n×n such that the constraint AA + BB = I n (15) is fulfilled, and the parameter η that defines the model for C η = C. Consider the following definition for f θ .

EQUATION

)

Assume that we can train f θ , as defined in Eq. ( 16), to map x ∼ N (x; 0, I 2n ) to s = f θ (x) with the probability distribution s ∼ p(s).

(17) Then, A and C = C η indeed fulfill the condition in Eq. ( 14) for h 1 , h 2 with joint probability distribution Eq. ( 10). In order to implement f θ by a neural network, we propose the architecture depicted in Fig. 1 . The first layer is linear and will be referred to as the dynamic layer in the following. The dynamic layer implements the multiplication with the matrix

EQUATION

)

The output of the dynamic layer is split into an upper half h 1 and a lower half h 2 and both halves are fed to the subnetwork that implements the observer function C η . The weights of the dynamic layer contain the matrices A, B. Thus, they can be trained along with η, by back-propagation of the stochastic gradient. However, we need to make sure that the stationarity constraint in Eq. ( 15) is not violated. This can be done by adding a regularizer to the loss function L of the Deep Generative Network. We thus substitute L(θ) by

EQUATION

where λ > 0 is a regularizer constant. Here, L refers to the variational bound for the VAE or to the Generator loss for the (Wasserstein) GAN. 2015)) as the foundation for implementing the observer function C η . We use it in combination with an affine layer at the input. This layer serves three purposes. Firstly, it accounts for changes of bases in the latent space in order to conform with Assumption 2. Secondly, it includes a bias so that the fixed point of the transition function can be transformed to the origin according to Proposition 1. Finally, it further reduces the latent dimension to make the search for the transition matrix A feasible.

We train C η and A simultaneously by means of the dynamic layer introduced in Section 4.3. To this end, we integrate the construction in Fig. 1 as the decoder of a VAE (without batch normalization) and the generator of a Wasserstein GAN (with batch normalization), respectively. As the critic (discriminator) network of the Wasserstein GAN, we use the original discriminator of the DCGAN but double the number of output channels in order to match the dimensions of the training frame pairs. As the encoder of the VAE, we use the discriminator of the DCGAN, but adapt the number of output channels to be 2n, where n is the latent dimension of the model. The latent dimension is set to n = 10 for all experiments.

Synthesis is performed by sampling from the MAR model described by

EQUATION

and mapping the latent states to the observation space by means of C η . The initial latent state h 0 is sampled from Gaussian white noise for the Wasserstein GAN experiments. For the VAE experiments, it is estimated by applying the encoder to an observation frame pair from the training sequence. PyTorch code for reproducing the experiments will be made available online upon publication, along with the entire set of results.

5.2. VARIATIONAL AUTOENCODER:

An isotropic standard Gaussian model N (s; f θ (x), σ 2 ) is assumed for the conditional data distribution. The Adam optimizer with a step size of 2.5 * 10 -4 is used to train the architecture. The regularizer constant is set to λ = 100. We carry out two series of experiments.

In the first series of experiments, we use the MNIST dataset to generate repeating sequences of hand written digits. We choose σ = 4 and use sequences of length 10000 for training the network for 25 epochs. Most of the number sequences we test can be well reproduced. Fig. 2 depicts the synthesis results for the sequences 0123401234... and 4567845678.... Occasional jumps occur, due to the stochastic nature of the model. An observation we make is that higher values of σ improve the probability of synthesizing the correct sequence, but decrease the variability of digit shapes.

To investigate the linearization of geometrical transformations, we employ Category 0 of the Small NORB dataset that contains pictures of miniature animals under 6 different lighting conditions, 9 elevational and 18 azimuthal poses. We train our architecture for 100 epochs with σ = √ 5 to synthesize sequences depicting azimuthal rotations of 20 • per frame. Generally, the synthesis provides good results. For each one of the five animals, the synthesis yields clearly recognizable rotation movements. However, the algorithm tends to confuse pairs of opposite poses, as can be observed in Fig. 3a , where the poses are flipped after the fourth frame in each row. With regard to Section 3.2, two explanations for this can be identified. On the one hand, it is possible that the assumptions of Proposition 1 are not fulfilled, and no diffeomorphism, i.e. a bijective mapping exists from the image manifold to the embedded space. In that case, the best we can do is finding a non-injective local diffeomorphism that assigns more than one point on the image manifold to a point in the embedded space. On the other hand, it is thinkable, that an actual diffeomorphism exists, but that the euclidean distance implicitly minimized by the VAE is a bad estimation of the geodesic distance on the manifold. Fig. 3b depicts the learned fixed point of the transition function. It can be thought of as a rotationally invariant structure under limited exposure of light from above.

5.3. WASSERSTEIN GENERATIVE ADVERSARIAL NETWORK:

We employ the Wasserstein GAN for dynamic texture synthesis experiments and set the regularizer to λ = 1. We firstly test our method on the cropped UCLA-50 dataset of grayscale DTs. The dataset contains 50 classes of 4 sequences of length 75 each. We use one class at a time for the experiments. The architecture is trained for 2000 epochs via RMSPROP with step size 2.5 * 10 -4 . Fig. 4 depicts the results for the classes candle and fountain-c-far.

Finally, we evaluate the method on three RGB sequences. We run the optimization via RMSPROP with step size 5 * 10 -5 for 900 epochs. Fig. 5 depicts the result for the firepot, springwater and waterfountain sequences taken from Xie et al. (2017) . Despite the low complexity of our model, we believe that the difference in quality compared to the state of the art results reported in Xie et al. ( 2017) is negligible.

6. CONCLUSION:

This work presents an approach to learn embedded MAR models from image sequences. We motivate the feasibility of this approach by introducing the concept of local linearizability and propose a joint learning procedure that employs deep generative models in combination with an additional linear component, the dynamic layer. We report first positive results on low-resolution visual processes, where a first-order Markov property can be assumed, and hope to shed some light on the nature of linearization. A possible future research direction is improving the theoretical understanding of linearizing representations and their applicability outside of stationary visual processes. (

EQUATION

This requirement is fulfilled if the Jacobi matrices of ϕ and φ coincide, i.e.

EQUATION

The Jacobian of φ at y * is given, according to he chain rule, by

EQUATION

where J Γ -1 ∈ R d×n is the Jacobian matrix of Γ -1 at Γ(y * ). A matrix Φ ∈ R n×n can be always found, such that Eq. ( 23) is fulfilled, if the columns and rows of J ϕ lie in the column space of J Γ -1 and the row space of J Γ , respectively. The column space of J Γ -1 coincides with the row space of J Γ , due to the identity J Γ J Γ -1 = I n . The statement of the proposition follows.

A.2 GENERALIZATION TO HIGHER-ORDER MARKOV PROCESSES

The proposed method is designed to linearize first-order Markov processes only. This is in line with our empirical evaluation, in which the method was not capable of synthesizing sequences with persisting motions of non-constant velocities, e.g. the videos of running animals in Xie et al. (2017) . However, the presented approach can be theoretically generalized to Markov processes of order m ≥ 1. We present a procedure that could be employed for this purpose. Further investigation is necessary to evaluate the feasibility and practical applicability of such an approach.

The general procedure is as follows.

1. Set up an appropriate block matrix model for the joint probability of m + 1 succeeding observations. 2. Build a dynamic layer that performs a multiplication with a lower-triangular (m + 1) × (m + 1) block matrix F .

3. Introduce regularizers to preserve the block Toeplitz structure of the covariance matrix, 4. Compute the MAR parameters from the resulting covariance matrix.

We illustrate the procedure for the case m = 2.

1. Due to the stationarity assumption, the covariance matrix for three succeeding observations h t , h t+1 , h t+2 has the form 

Σ 2 =   I n Cov(

2. The matrix describing the dynamic layer has the form

EQUATION

The outputs of the dynamic layer will thus have the distribution

EQUATION

We thus need to achieve

EQUATION

)

3. This can be ensured by using the regularizer

EQUATION

with λ 1 , λ 2 , λ 3 > 0. 4. A second-order MAR model has the form h t+2 = A 0 h t + A 1 h t+1 + Bv t , v t ∼ N (v t , 0, I n ).

(30)

:



Figures and Tables:

Figure 1: Generative network with dynamic layer. The dashed arrow indicates the trainable parameters

Figure 2: Synthesis result for sequences of MNIST digits with VAE

Figure 3: Synthesis result for NORB sequences (a) and learned fixed point (b)

Figure 5: Synthesis result for RGB sequences. Frames are rescaled to match the original aspect ratio

h t , h t+1 ) Cov(h t , h t+2 ) Cov(h t , h t+1 ) I n Cov(h t , h t+1 ) Cov(h t , h t+2 )Cov(h t , h t+1 ) I

Table content: <table/>

Peer review (iclr2019):
 This paper proposes a new deep generative model for sequences, particularly image sequences and video. By utilizing a linear structure in part of the model, this work aims to offer more efficient sampling, while also aiming for a simpler training procedure. Unfortunately there are some important mistakes in the mathematical justification of the method provided in Section 3. One is the claim that "Eq. (1) is a special case of Eq. (2) ..." and later "The model Eq. (2) describes a much broader class of visual processes than Eq. (1)." A counterexample is to take n=2, d=1, and A = [[ cos(.pi/4) -sin(.pi/4) ] [ sin(.pi/4) cos(.pi/4) ]] C = [[ 1 0 ]] .bar{y} = 0 Since C^+ is [[ 1 ] [ 0 ]], the provided formula does not recover the periodic dynamics of Eq. (1) and is instead just a mean-reverting random walk. To make the counterexample more extreme, we can take n=1000 and d=1 with an A matrix that exhibits 500 distinct frequencies (with eigenvalues in complex conjugate pairs) and with C = [[ 1 1 ... 1 ]], then clearly in Eq. (1) the dynamical behavior is a sum of 500 oscillators, which can*t be recovered by Eq. (2) (in which y_{t+1} depends only on y_{t}). More precisely, with d=1 Eq. (1) can represent a stationary process with any rational power spectral density (up to the right degree, depending on n) but Eq. (2), being first-order Markov on the observation sequence (stated as Assumption 1) can*t. (It*s incorrectly implied near the bottom of p. 2 that processes following the dynamics of Eq. (1) are first-order Markov on the y sequence, but they*re not; that*s the purpose of having a latent state.) There are some other issues in Section 3 that are not mistakes but are unclear. It*s inaccurate to describe C as a diffeomorphism (just before Sec 3.2) and .Gamma : R^d .to R^n as a diffeomorphism, since diffeomorphisms must be invertible and here n < d, though it*s clear from context that you mean their appropriate restrictions are diffeomorphisms. Remark 1 states the process noise must have covariance I - AA* without stating that, in addition to choosing a basis for the latent space in which states have a marginal standard Gaussian distribution, it must be chosen so that I - AA* is positive definite. The initial state distribution is not taken to be the stationary distribution of the latent state dynamics, yet a basis can only be chosen to satisfy Assumption 2 in stationarity. Some of the mathematics in Section 3.2 is unnecessary; .Gamma essentially serves as a coordinate chart, and notions of fixed points on manifolds (as in Proposition 1) is textbook stuff. In the description of the learning procedure, no overall objective is specified, and the treatment of the constraint in Eq. (15) is weak (a squared Frobenius penalty is simply added into the objective). One consequence of these correctness and clarity issues is that it*s hard to understand how to relate this work to other ideas. In addition to having a general related works section, this paper should contrast the proposed model and fitting procedure with others. For example, how does this compare to fitting a nonlinear AR(1) model? Or an autoregressive model for p(y_{t+1} | y_t, ..., y_1) parameterized by an RNN? Or any of the various latent variable generative models for videos based on VAE architectures? The other main issue with the current paper is a lack of baselines. As far as I can tell, there are no baselines of comparison in the experiments section. One clear baseline would be a nonlinear AR(1) model parameterizing p(h_{t+1} | h_t) directly, perhaps as a Gaussian with mean and variance that depend on h_t through a neural network. But there are lots of video models at this point, of both the autoregressive and latent variable flavor, and this paper must include some comparisons. Moreover, the modeling tasks considered here seem nonstandard. To make up for the clarity and correctness issues with the mathematical justification of the approach in the paper as currently written, the results would have to be especially impressive compared to strong baselines, but the experiments section does not give enough information to make that evaluation. Overall, there are some interesting ideas here but the paper needs more work to polish the mathematics and to produce convincing experiments. Some minor suggestions: * It*s confusing to write "linear, Gaussian models of observed sequences" in the abstract, since the models are not (completely) linear and the resulting distribution on observations is not Gaussian. * Since the model in Eq. (3) is one of those considered in the SVAE paper (namely in Sec 2.2 of the SVAE paper), and probably also in earlier works, it*s worth citing some related work around Eq. (3). * It might be interesting to draw a comparison to subspace identification methods like 4SID, as described in "Subspace Identification for Linear Systems" by van Overschee and de Moor, or spectral methods for HMMs, as in "A Spectral Algorithm for learning Hidden Markov Models" by Hsu, Kakade, and Zhang. * It might be helpful to rename "MAR" to "VAR" for "vector autoregressive", since "MA" often stands for "moving average" in the context of VARMA models. * Writing "due to the curse of dimensionality" in Sec 3.2 is vague; what do you mean there?