Self-supervised contrastive learning performs non-linear system identification

## Abstract

Abstract Self-supervised learning (SSL) approaches have brought tremendous success across many tasks and domains. It has been argued that these successes can be attributed to a link between SSL and identifiable representation learning: Temporal structure and auxiliary variables ensure that latent representations are related to the true underlying generative factors of the data. Here, we deepen this connection and show that SSL can perform system identification in latent space. We propose DynCL , a framework to uncover linear, switching linear and non-linear dynamics under a non-linear observation model, give theoretical guarantees and validate them empirically.
 Code: github.com/dynamical-inference/dyncl

## 1 Introduction

The identification and modeling of dynamics from observational data is a long-standing problem in machine learning, engineering and science.
A discrete-time dynamical system with latent variables ${\bm{x}}$, observable variables ${\bm{y}}$, control signal ${\bm{u}}$ and noise ${\bm{\varepsilon}},{\bm{\nu}}$ can take the form

$$ $\displaystyle{\bm{x}}_{t+1}$ $\displaystyle={\bm{f}}({\bm{x}}_{t})+{\bm{u}}_{t}+{\bm{\varepsilon}}_{t}$ (1) $\displaystyle{\bm{y}}_{t}$ $\displaystyle={\bm{g}}({\bm{x}}_{t})+{\bm{\nu}}_{t}.$ $$

and we aim to infer the functions ${\bm{f}}$ and ${\bm{g}}$ from a time-series of observations and, when available, control signals.
Numerous algorithms have been developed to tackle special cases of this problem formulation, ranging from classical system identification methods [, ] to recent generative models [, , ]. Yet, it remains an open challenge to improve the generality, interpretability and efficiency of these inference techniques, especially when ${\bm{f}}$ and ${\bm{g}}$ are non-linear functions.

Contrastive learning (CL) and next-token prediction tasks have become important backbones of modern machine learning systems for learning from sequential data, proving highly effective for building meaningful latent representations [, , , , , , ]. An emerging view is a connection between these algorithms and learning of “world models” [, ]. Yet, non-linear system identification in such sequence-learning algorithms is poorly theoretically studied.

In this work, we revisit and extend contrastive learning in the context of system identification. We uncover several surprising facts about its out-of-the-box effectiveness in identifying dynamics and unveil common design choices in SSL systems used in practice.
Our theoretical study extends identifiability results [, , , , ] for CL towards dynamical systems.
While our theory makes several predictions about capabilities of standard CL, it also highlights shortcomings.
To overcome these and enable interpretable dynamics inference across a range of data generating processes, we propose a general framework for linear and non-linear system identification with CL (Figure 1).

Figure: Figure 1: (a), the assumed data-generating process: ${\bm{y}}$ represents the observable input variables, ${\bm{x}}$ denotes the latent variables, and ${\bm{u}}$ is the control input. (b), general formulation of our method: It consists of an encoder ${\bm{h}}$ that is shared across the reference, positive, and negative samples. The method also includes a dynamics model $\hat{\bm{f}}$. Additionally, a (possibly latent) variable ${\bm{z}}$ can be used to parameterize the dynamics model (see section 4).

#### Background.

An influential motivation of our work is Contrastive Predictive Coding [CPC; ]. CPC can be recovered as a special case of our framework when using an RNN dynamics model. Related works have emerged across different modalities: wav2vec [], TCN [] and CPCv2 [].
In the field of system identification, notable approaches include the Extended Kalman Filter (EKF) [] and NARMAX []. Additionally, several works have also explored generative models for general dynamics [] and switching dynamics, e.g. rSLDS []. In the Nonlinear ICA literature, identifiable algorithms for time-series data, such as Time Contrastive Learning [TCL; ] for non-stationary processes and Permutation Contrastive Learning [PCL; ] for stationary data have been proposed, with recent advances like SNICA [] for more generally structured data-generating processes.

In contrast to previous work, we focus on bridging time-series representation learning through contrastive learning with the identification of dynamical systems, both theoretically and empirically. Moreover, by not relying on an explicit data-generating model, our framework offers greater flexibility. We extend and discuss the connections to related work in more detail in Appendix C.

#### Contributions.

We extend the existing theory on contrastive learning for time series learning and make adaptations to common inference frameworks. We introduce our CL variant in section 2, and give an identifiability result for both the latent space and the dynamics model in section 3. These theoretical results are later empirically validated. We then propose a practical way to parameterize switching linear dynamics in section 4 and demonstrate that this formulation corroborates our theory for both switching linear system dynamics and non-linear dynamics in sections 5-6.

## 2 Contrastive learning for time-series

In contrastive learning, we aim to model similarities between pairs of data points (Figure 1b). Our full model $\psi$ is specified by the log-likelihood

$$ $\log p_{\psi}({\bm{y}}|{\bm{y}}^{+},N)=\psi({\bm{y}},{\bm{y}}^{+})-\log\!\! \sum_{{\bm{y}}^{-}\in N\cup\{{\bm{y}}^{+}\}}\!\!\exp(\psi({\bm{y}},{\bm{y}}^{- })).$ (2) $$

where ${\bm{y}}$ is often called the reference or anchor sample, ${\bm{y}}^{+}$ is a positive sample, and ${\bm{y}}^{-}\in N$ are negative examples.
The model $\psi$ itself is parameterized as a composition of an encoder, a dynamics model, and a similarity function and will be defined further below. We fit the model by minimizing the negative log-likelihood on the time series,

$$ $\min_{\psi}\mathcal{L}[\psi]=\min_{\psi}\ \mathbb{E}_{t,t_{1},\dots,t_{M}\sim U (1,T)}[-\log p_{\psi}({\bm{y}}_{t+1}|{\bm{y}}_{t},\{{\bm{y}}_{t_{m}}\}_{m=1}^{ M})]$ (3) $$

where positive examples are just adjacent points in the time-series, and negatives are sampled uniformly across the dataset. $U(1,T)$ denotes a uniform distribution across the discrete timesteps.

To attain favourable properties for identifying the latent dynamics, we carefully design the hypothesis class for $\psi$. The motivation for this particular design will become clear later. To define the full model, a composition of several functions is necessary. Recall from Eq. 1 that the dynamics model is given as ${\bm{f}}$ and the mixing function is ${\bm{g}}$. Correspondingly, our model is composed of the encoder ${\bm{h}}:\mathbb{R}^{D}\mapsto\mathbb{R}^{d}$ (de-mixing), the dynamics model $\hat{\bm{f}}:\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$, the similarity function $\phi:\mathbb{R}^{d}\times\mathbb{R}^{d}\mapsto\mathbb{R}$ and a correction term $\alpha:\mathbb{R}^{d}\mapsto\mathbb{R}$. We define their composition as(^1^11 Note that we can equivalently write $\phi(\tilde{\bm{h}}({\bm{x}})),\tilde{\bm{h}}^{\prime}({\bm{x}}^{\prime}))$ using two asymmetric encoder functions, see additional results in Appendix D.)

$$ $\psi({\bm{y}},{\bm{y}}^{\prime}):=\phi(\hat{\bm{f}}({\bm{h}}({\bm{y}})),{\bm{h }}({\bm{y}}^{\prime}))-\alpha({\bm{y}}^{\prime}),$ (4) $$

and call the resulting algorithm DynCL.
Intuitively, we obtain two observed samples $({\bm{y}},{\bm{y}}^{\prime})$ which are first mapped to the latent space, $({\bm{h}}({\bm{y}}),{\bm{h}}({\bm{y}}^{\prime}))$. Then, the dynamics model is applied, and the resulting points are compared through the similarity function $\phi$. A key insight is that the similarity function $\phi$ will be informed by the form of control signal ${\bm{u}}_{t}$. In the simplest form, the control can be chosen as isotropic Gaussian noise, which results in a negative squared Euclidean norm for $\phi$.
Note the additional term $\alpha({\bm{y}}^{\prime})$ is a correction applied to account for non-uniform marginal distributions. It can be parameterized as a kernel density estimate (KDE) with $\log\hat{q}({\bm{h}}({\bm{y}}^{\prime}))\approx\log q({\bm{x}}^{\prime})$ around the datapoints. While on very special cases, the KDE makes a difference in empirical performance (Appendix B, Fig. 9) and is required for our theory, we found that on the time-series datasets considered, empirically it was possible to drop this term without loss in performance (i.e., $\alpha=0$).

## 3 Structural identifiability of non-linear latent dynamics

We now study the aforementioned model theoretically.
The key components of our theory along with our notion of linear identifiability [, ] are visualized in Figure 2.
We are interested in two properties. First, linear identifiability of the latent space: The composition of mixing function ${\bm{g}}$ and model encoder ${\bm{h}}$ should recover the ground-truth latents up to a linear transform. Second, identifiability of the (non-linear) dynamics model: We would like to relate the estimated dynamics ${\hat{{\bm{f}}}}$ to the underlying ground-truth dynamics ${\bm{f}}$. This property is also called *structural identifiability* []. Our model operates on a subclass of Eq. 1 with the following properties:

Data-generating process.
A discrete-time dynamical system is defined as

$$ $\displaystyle{\bm{x}}_{t+1}={\bm{f}}({\bm{x}}_{t})+{\bm{u}}_{t},\qquad{\bm{y}} _{t}={\bm{g}}({\bm{x}}_{t}),$ (5) $$

where ${\bm{x}}_{t}\in\mathbb{R}^{d}$ are latent variables, ${\bm{f}}:\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$ is a bijective dynamics model, ${\bm{u}}_{t}\in\mathbb{R}^{d}$ a control or noise signal, and ${\bm{g}}:\mathbb{R}^{d}\mapsto\mathbb{R}^{D}$ is a non-linear injective mapping from latents to observables ${\bm{y}}_{t}\in\mathbb{R}^{D}$, $d\leq D$.
We sample a total number of $T$ timesteps.

We proceed by stating our main result:

###### Theorem 1 (Contrastive estimation of non-linear dynamics) .

Assume that

- •
A time-series dataset $\{{\bm{y}}_{t}\}_{t=1}^{T}$ is generated according to the ground-truth dynamical system in Eq. 5 with a bijective dynamics model ${\bm{f}}$ and an injective mixing function ${\bm{g}}$.
- •
The system input follows an iid normal distribution, $p({\bm{u}}_{t})=\mathcal{N}({\bm{u}}_{t}|0,{\bm{\Sigma}}_{u})$.
- •
The model $\psi$ is composed of an encoder ${\bm{h}}$, a dynamics model ${\hat{{\bm{f}}}}$, a correction term $\alpha$, and the similarity metric $\phi({\bm{u}},{\bm{v}})=-\|{\bm{u}}-{\bm{v}}\|^{2}$ and attains the global minimizer of Eq. 3.

Then, in the limit of $T\rightarrow\infty$ for any point ${\bm{x}}$ in the support of the data marginal distribution:

- (a)
The composition of mixing and de-mixing ${\bm{h}}({\bm{g}}({\bm{x}}))={\bm{L}}{\bm{x}}+{\bm{b}}$ is a bijective affine transform, and ${\bm{L}}={\bm{Q}}{\bm{\Sigma}}_{u}^{-1/2}$ with unknown orthogonal transform ${\bm{Q}}\in\mathbb{R}^{d\times d}$ and offset ${\bm{b}}\in\mathbb{R}^{d}$.
- (b)
The estimated dynamics $\hat{\bm{f}}$ are bijective and identify the true dynamics ${\bm{f}}$ up to the relation ${\hat{{\bm{f}}}}({\bm{x}})={\bm{L}}{\bm{f}}({\bm{L}}^{-1}({\bm{x}}-{\bm{b}}))+
{\bm{b}}$.

###### Proof.

See Appendix A for the full proof, and see Figure 2 for a graphical intuition of both results.
∎

Figure: Figure 2: Graphical intuition behind Theorem 1. Left: the ground truth latent space is mapped to observables through the injective mixing function ${\bm{g}}$. Our model maps back into the latent space. The composition of mixing and de-mixing by the model is an affine transform. Right: dynamics in the ground-truth space are mapped to the latent space. By observing variations introduced by the control signal ${\bm{u}}$, our model is able to infer the ground-truth dynamics up to an affine transform.
Refer to caption: x1.png

With this main result in place, we can make statements for several systems of interest; specifically linear dynamics in latent space:

###### Corollary 1 .

Contrastive learning without dynamics model (${\hat{{\bm{f}}}}\!=\!1$) cannot identify latent dynamics.

This is because due to the indeterminacy in the theorem, we would require that ${\bm{L}}^{\top}{\bm{L}}={\bm{\Sigma}}_{u}^{-1}={\bm{A}}$ for the case of an identity dynamics model, which reduces the system to the case of Brownian motion. This setting corresponds to purely symmetric contrastive learning systems.
We can fix this case by either decoupling the two encoders (Appendix D), or take a more structured approach and parameterize a dynamics model with a dynamics matrix:

###### Corollary 2 .

For a linear dynamical system ${\bm{f}}({\bm{x}})={\bm{A}}{\bm{x}}$ and dynamics model ${\hat{{\bm{f}}}}({\bm{x}})=\hat{\bm{A}}{\bm{x}}$, we identify the latents up to ${\bm{h}}({\bm{g}}({\bm{x}}))={\bm{L}}{\bm{x}}+{\bm{b}}$ and dynamics with $\hat{\bm{A}}={\bm{L}}{\bm{A}}{\bm{L}}^{-1}$.

This means that simultaneously fitting the system dynamics and encoding model allows us to recover the system matrix up to an indeterminacy.

###### Theorem 1 (Contrastive estimation of non-linear dynamics) .

###### Proof.

###### Corollary 1 .

###### Corollary 2 .

## 4 ∇ ∇ \nabla ∇ -SLDS: Towards non-linear dynamics estimation

Figure: Figure 3: The core components of the $\nabla$-SLDS model is the differentiable parameterization of the switching process.
Refer to caption: x2.png

#### Piecewise linear approximation of dynamics.

Our theoretical results suggest that contrastive learning allows the fitting of non-linear bijective dynamics. This is a compelling result, but in practice it requires the use of a powerful, yet easy to parameterize dynamics model. One option is to use an RNN [, ] or Transformer [] model to perform this link across timescales.
An alternative option is to linearize the system, which we propose in the following.

We propose a new forward model for differentiable switching linear dynamics ($\nabla$-SLDS) in latent space. The estimation is outlined in Figure 3.
This model allows fast estimation of switching dynamics and can be easily integrated into the DynCL algorithm.
The dynamics model has a trainable bank ${\bm{\mathsfit{W}}}=[{\bm{W}}_{1},\dots,{\bm{W}}_{K}]$ of possible dynamics matrices. $K$ is a hyperparameter. The dynamics depend on a latent variable $k_{t}$ and are defined as

$$ $\displaystyle\hat{\bm{f}}({\bm{x}}_{t};{\bm{\mathsfit{W}}},k_{t})={\bm{W}}_{k_ {t}}{\bm{x}}_{t},\quad k_{t}=\mathrm{argmin}_{k}\|{\bm{W}}_{k}{\bm{x}}_{t}-{ \bm{x}}_{t+1}\|^{2}.$ (6) $$

During training, we approximate the argmin using the Gumbel-Softmax trick [] without hard sampling,

$$ $\displaystyle\hat{\bm{f}}({\bm{x}}_{t};{\bm{\mathsfit{W}}},{\bm{z}}_{t})=(\sum _{k=1}^{K}z_{t,k}{\bm{W}}_{k}){\bm{x}}_{t},\ z_{t,k}=\frac{\exp(\lambda_{k}/ \tau)}{\sum_{j}\exp(\lambda_{j}/\tau)},\ \lambda_{k}=\frac{1}{\|{\bm{W}}_{k}{ \bm{x}}_{t}-{\bm{x}}_{t+1}\|^{2}}+g_{k}.$ (7) $$

Note that the dynamics model $\hat{\bm{f}}({\bm{x}}_{t};{\bm{\mathsfit{W}}},{\bm{z}}_{t})$
depends on an additional latent variable ${\bm{z}}_{t}=[z_{t,1},\dots,z_{t,K}]^{\top}$ which contains probabilities to parametrize the dynamics. During inference, we can relate $k_{t}=\operatorname*{arg\,max}_{k}z_{t,k}$ to obtain hard assignments.
The variables $g_{k}$ are samples from the Gumbel distribution [] and we use a temperature $\tau$ to control the smoothness of the resulting probabilities. During pilot experiments, we found that the reciprocal parameterization of the logits outperforms other choices for computing an argmin, like flipping the sign.

#### From linear switching to non-linear dynamics.

Non-linear system dynamics of the general form in Eq. 5 can be approximated using our switching model.
Consider a continuous-time non-linear dynamical system with latent dynamics
$\dot{\bm{x}}={\bm{f}}({\bm{x}})$,
which we approximate around reference points $\{{\tilde{\bm{x}}}_{k}\}_{k=1}^{K}$ using a first-order tailor expansion,

$$ ${\bm{f}}({\bm{x}})\approx\tilde{\bm{f}}({\bm{x}})={\bm{f}}({\tilde{\bm{x}}}_{k })+{\bm{J}}_{\bm{f}}({\tilde{\bm{x}}}_{k})({\bm{x}}-{\tilde{\bm{x}}}_{k})$ (8) $$

where we denote the Jacobian matrix of ${\bm{f}}$ with ${\bm{J}}_{\bm{f}}$.
We evaluate the equation at each point $t$ using the best reference point ${\tilde{\bm{x}}}_{k}$. We obtain system matrices and control signals

$$ $\displaystyle{\bm{A}}_{k}={\bm{J}}_{\bm{f}}({\tilde{\bm{x}}}_{k}),\quad{\bm{w} }_{k}={\bm{f}}({\tilde{\bm{x}}}_{k})-{\bm{J}}_{\bm{f}}({\tilde{\bm{x}}}_{k}){ \tilde{\bm{x}}}_{k}$ (9) $$

which can be modeled with the $\nabla$-SLDS model ${\hat{{\bm{f}}}}({\bm{x}}_{t};k_{t})$:

$$ $\displaystyle{\bm{x}}_{t+1}=({\bm{A}}_{k_{t}}{\bm{x}}_{t}+{\bm{w}}_{k_{t}})+{ \bm{u}}_{t}=:{\hat{{\bm{f}}}}({\bm{x}}_{t};k_{t})+{\bm{u}}_{t}.$ (10) $$

While a theoretical guarantee for this general case is beyond the scope of this work, we give an empirical evaluation on Lorenz attractor dynamics below.
Note, as the number of “basis points” of $\nabla$-SLDS approaches the number of timesteps, we could trivially approach perfect estimation capability of the latents as we store the exact value of ${\bm{f}}$ at every point. However, this comes at the expense of having less points to estimate each individual dynamics matrix. Empirically, we used 100–200 matrices for datasets of 1M samples.

## 5 Experiments

To verify our theory, we implement a benchmark dataset for studying the effects of various model choices. We generate data with 1M samples per time-series, either as a single sequence or across multiple trials. Our experiments rigorously evaluate different variants of contrastive learning algorithms within our DynCL framework.

#### Data generation.

Our datasets are generated following Eq. 5 by simulating latent variables ${\bm{x}}$ that evolve according to a dynamical system. These latent variables are then passed through a nonlinear mixing function ${\bm{g}}$ to produce the observable data ${\bm{y}}$. The mixing function ${\bm{g}}$ consists of a nonlinear injective component which is parameterized by a randomly initialized 4-layer MLP [], and a linear map to a 50-dimensional space. The final mixing function is defined as their composition. We ensure the injectivity of the resulting function by monitoring the condition number of each matrix layer, following previous work [, ].

#### LDS.

We simulate 1M datapoints in 3d space following ${\bm{f}}(x_{t})={\bm{A}}x_{t}$ with system noise standard deviation $\sigma_{u}=0.01$ and choose ${\bm{A}}$ to be an orthogonal matrix to ensure stable dynamics with all eigenvalues equal to 1. We do so by taking the product of multiple rotation matrices, one for each possible plane to rotate around with rotation angles being randomly chosen to be -5° or 5°.

#### SLDS.

We simulate switching linear dynamical systems with ${\bm{f}}({\bm{x}}_{t};k_{t})={\bm{A}}_{k_{t}}{\bm{x}}_{t}$ and system noise standard deviation $\sigma_{u}=0.0001$. We choose ${\bm{A}}_{k}$ to be an orthogonal matrix ensuring that all eigenvalues are 1, which guarantees system stability. Specifically, we set ${\bm{A}}_{k}$ to be a rotation matrix with varying rotation angles (5°, 10°, 20°). The latent dimensionality is 6. The number of samples is 1M. We use 1000 trials, and each trial consists of 1000 samples. The noise ${\bm{u}}$ is sampled from a Normal distribution with varying standard deviations $\sigma$. We use $k={0,1,\dots,K}$ distinct modes following a mode sequence $i_{t}$. The mode sequence $i_{t}$ follows a Markov chain with symmetric transition matrix and uniform prior: $i_{0}\sim\text{Cat}(\pi)$, where $\pi_{j}=\frac{1}{K}$ for all $j$. At each timestep, $i_{t+1}\sim\text{Cat}(\Pi_{i_{t}})$, where $\Pi$ is a transition matrix with uniform off-diagonal probabilities. The transition probability is $10^{-4}$.
Example data is visualized in Figure 4 and Appendix E.

#### Non-linear dynamics.

We simulate 1M points of a Lorenz system, with equations

$$ $\displaystyle{\bm{f}}({\bm{x}}_{t})={\bm{x}}_{t}+dt[\sigma(x_{2,t}-x_{1,t}),x_ {1,t}((\rho-x_{3,t})-x_{2,t}),(x_{1,t}x_{2,t}-\beta x_{3,t})]^{\top}$ (11) $$

and parameters $\sigma=10,\,\beta=\frac{8}{3},\,\rho=28$ and system noise standard deviation $\sigma_{u}=0.001$. The observable data, ${\bm{y}}$. We then apply our non-linear mixing function as for other datasets.
The time step $dt$ is varied during experiments.

#### Model estimation.

For the feature encoder ${\bm{h}}$, baseline and our model we use an MLP with three layers followed by GELU activations []. Each layer has 512 units. We train on batches with 2048 samples each (reference and positive) and use 32,768 negative samples. We use the Adam optimizer [] with learning rates $3\times 10^{-4}$ for LDS data, $10^{-3}$ for SLDS data, and $10^{-4}$ for Lorenz system data.
Our baseline model is standard self-supervised contrastive learning with the InfoNCE loss, akin to the CEBRA-time model (with symmetric encoders, i.e., without a dynamics model) used by Schneider et al. []. For DynCL, we add an LDS or $\nabla$-SLDS dynamics model for fitting.

**Table 1: Comprehensive categorization of ground-truth dynamical processes ${\bm{f}}$ and model configurations ${\hat{{\bm{f}}}}$. Note, the mixing function ${\bm{g}}$ is always assumed non-linear, and ${\bm{h}}$ is always parameterized as a neural network. For all metrics, we report mean and standard deviation across 3 datasets (5 for Lorenz) and 3 experiment repeats.**
| data |  | model | results |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| ${\bm{f}}$ | $p({\bm{u}})$ |  | $\hat{\bm{f}}$ | theory | %$R^{2}$ $\uparrow$ | LDS$\downarrow$ |
| identity | Normal |  | identity | ✓ | 99.56 $\pm$ 0.21 | 0.00 $\pm$ 0.00 |
| identity | Normal |  | LDS | ✓ | 99.31 $\pm$ 0.43 | 0.04 $\pm$ 0.01 |
| LDS (low $\Delta t$) | Normal (large $\sigma$) |  | identity | – | 89.23 $\pm$ 4.46 | 8.53 $\pm$ 0.05 |
| LDS | Normal |  | identity | ✗ | 73.56 $\pm$ 24.45 | 21.24 $\pm$ 0.31 |
| LDS | Normal |  | LDS | ✓ | 99.03 $\pm$ 0.41 | 0.38 $\pm$ 0.34 |
| LDS | Normal |  | GT | ✓ | 99.46 $\pm$ 0.39 | 0.17 $\pm$ 0.06 |
|  |  |  |  |  | %$R^{2}$ $\uparrow$ | %dyn$R^{2}$ $\uparrow$ |
| SLDS | Normal |  | identity | ✗ | 76.80 $\pm$ 7.40 | 85.47 $\pm$ 8.07 |
| SLDS | Normal |  | $\nabla$-SLDS | (✓) | 99.52 $\pm$ 0.05 | 99.93 $\pm$ 0.01 |
| SLDS | Normal |  | GT | (✓) | 99.20 $\pm$ 0.10 | 99.97 $\pm$ 0.00 |
| Lorenz (small $\Delta t$) | Normal (large $\sigma$) |  | identity | – | 99.74 $\pm$ 0.36 | 99.94 $\pm$ 0.07 |
| Lorenz (small $\Delta t$) | Normal (large $\sigma$) |  | LDS | – | 98.31 $\pm$ 2.55 | 97.21 $\pm$ 5.90 |
| Lorenz (small $\Delta t$) | Normal (large $\sigma$) |  | $\nabla$-SLDS | – | 94.14 $\pm$ 4.34 | 94.20 $\pm$ 6.57 |
| Lorenz | Normal |  | identity | ✗ | 40.99 $\pm$ 8.58 | 27.02 $\pm$ 8.72 |
| Lorenz | Normal |  | LDS | ✗ | 81.20 $\pm$ 16.93 | 80.30 $\pm$ 14.13 |
| Lorenz | Normal |  | $\nabla$-SLDS | (✓) | 94.08 $\pm$ 2.75 | 93.91 $\pm$ 5.32 |

#### Evaluation metrics.

Our metrics are informed by the result in Theorem 1 and measure empirical identifiability up to affine transformation of the latent space and its underlying linear or non-linear dynamics.
To account for the affine indeterminacy, we explicitly estimate ${\bm{L}}$, ${\bm{b}}$ for ${\bm{x}}={\bm{L}}\hat{\bm{x}}+{\bm{b}}$ which allows us to transform recovered latents $\hat{\bm{x}}$ into the space of ground truth latents ${\bm{x}}$.
In those cases, where the inverse transform $\hat{\bm{x}}={\bm{L}}^{-1}({\bm{x}}-{\bm{b}})$ is required, for the purpose of numerical stability we estimate it from data rather than computing an explicit inverse of ${\bm{L}}$.
This results in estimates for ${\bm{L}}_{1}$,${\bm{b}}_{1}$ and ${\bm{L}}_{2}$,${\bm{b}}_{2}$, which we fit via linear regression:

$$ $\displaystyle\min_{{\bm{L}}_{1},{\bm{b}}_{1}}\|\hat{{\bm{x}}}-({\bm{L}}_{1}{ \bm{x}}+{\bm{b}}_{1})\|_{2}^{2}\qquad\text{ and }\qquad\min_{{\bm{L}}_{2},{\bm {b}}_{2}}\|{\bm{x}}-({\bm{L}}_{2}\hat{{\bm{x}}}+{\bm{b}}_{2})\|_{2}^{2}.$ (12) $$

To evaluate the identifiability of the representation, we measure the $R^{2}$ score between the true latents and the optimally aligned recovered latents (and with “r2_score” denote the scikit-learn implementation with default parameters which averages over the latent dimensions with a uniform weight):

$$ $R^{2}({\bm{x}},\hat{\bm{x}})=\text{r2_score}({\bm{x}},{\bm{L}}_{2}\hat{\bm{x} }+{\bm{b}}_{2}).$ (13) $$

We propose two metrics as direct measures of identifiability for the recovered linear dynamics $\hat{\bm{f}}=\hat{\bm{A}}$. First, the LDS error, which is suitable only for linear dynamics models, computes the norm of the difference between the true dynamics matrix ${\bm{A}}$ and the estimated dynamics matrix $\hat{\bm{A}}$ by accounting for the linear transformation between the true and recovered latent spaces. The LDS error (related to the metric for Dynamical Similarity Analysis, []) is then computed as (cf. Corollary 2):

$$ $\text{LDS}({\bm{A}},\hat{\bm{A}})=\|{\bm{A}}-{\bm{L}}_{1}\hat{\bm{A}}{\bm{L}}_ {2}\|_{F}\approx\|{\bm{A}}-{\bm{L}}^{-1}\hat{\bm{A}}{\bm{L}}\|_{F}$ (14) $$

As a second, more general identifiability metric for the recovered dynamics $\hat{\bm{f}}$, we introduce dyn$R^{2}$, which builds on Theorem 1 to evaluate the identifiability of non-linear dynamics. This metric computes the $R^{2}$ between the predicted dynamics $\hat{\bm{f}}$ and the true dynamics ${\bm{f}}$, corrected for the linear transformation between the two latent spaces. Specifically, motivated by Theorem 1(b), we compute:

$$ $\text{dyn}R^{2}({\bm{f}},\hat{\bm{f}})=\text{r2_score}(\hat{\bm{f}}(\hat{\bm{ x}}),{\bm{L}}_{1}{\bm{f}}({\bm{L}}_{2}\hat{\bm{x}}+{\bm{b}}_{2})+{\bm{b}}_{1})$ (15) $$

Finally, when evaluating switching linear dynamics, we compute the accuracy for assigning the correct modes at any point in time. To compute the cluster accuracy in the case of SLDS ground truth dynamics, we leverage the Hungarian algorithm to match the estimated latent variables modeling mode switches to the ground truth modes, and then proceed to compute the accuracy.

All metrics are estimated on the dataset the model is fit on. See Appendix F for additional discussion on estimating metrics on independently sampled dynamics.
See Appendix G for additional discussion on the Dyn$R^{2}$ metric and additional variants.

#### Implementation.

Experiments were carried out on a compute cluster with A100 cards. On each card, we ran $\sim$3 experiments simultaneously. The experiments run for this paper comprised about 120 days of A100 compute time. We will open source our benchmark suite for identifiable dynamics learning upon publication of the paper.

## 6 Results

Figure: Figure 4: Switching linear dynamics: (a) example ground-truth dynamics in latent space for four matrices ${\bm{A}}_{k}$. (b) $R^{2}$ metric for different noise levels as we increase the angles used for data generation. We compare a baseline (no dynamics) to $\nabla$-SLDS and a model fitted with ground-truth dynamics. (c) cluster accuracies for models shown in (b).
Refer to caption: x3.png

Figure: Figure 5: Contrastive learning of 3D non-linear dynamics following a Lorenz attractor model. (a), left to right: ground truth dynamics for 10k samples with $dt=0.0005$ and $\sigma=0.1$, estimation results for baseline (identity dynamics), DynCL with $\nabla$-SLDS, estimated mode sequence. (b), empirical identifiability ($R^{2}$) between baseline (BAS) and $\nabla$-SLDS for varying numbers of discrete states $K$. (c, d), same layout but for $dt=0.01$ and $\sigma=0.001$.
Refer to caption: x4.png

### 6.1 Verification of the theory for linear dynamics

A variety of settings are relevant and interesting for system identification. We summarize our results in Table 1. In the main paper, we consider the Euclidean case which is most relevant for practical dynamical systems and more general. We provide additional results corroborating our theory with empirical results for the case of a vMF conditional and dot-product similarity in Appendix D.

#### Suitable dynamics models enable identification of latents and dynamics.

For all considered classes of models, we show in Table 1 that DynCL effectively identifies the correct dynamics. For linear dynamics (LDS), DynCL reaches an $R^{2}$ of 99.0%, close to the oracle performance (99.5%). Most importantly, the LDS error of our method (0.38) is substantially closer to the oracle (0.17) compared to the baseline model (21.24).
In the case of switching linear dynamics (SLDS), DynCL also shows strong performance, both in terms of latent $R^{2}$ (99.5%) and dynamics $R^{2}$ (99.9%) outperforming the respective baselines (76.8% $R^{2}$ and 85.5% dynamics $R^{2}$).
For non-linear dynamics, the baseline model fails entirely (41.0%/27.0%), while $\nabla$-SLDS dynamics can be fitted with 94.1% $R^{2}$ for latents and $93.9$% dynamics $R^{2}$. We also clearly see the strength of our piecewise-linear approximation, as the LDS dynamics models only reaches $81.2\%$ latent identifiability and $80.3\%$ dynamics $R^{2}$.

#### Learning noisy dynamics does not require a dynamics model.

If the variance of the distribution for ${\bm{u}}_{t}$ dominates the changes actually introduced by the dynamics, we find that the baseline model is also able to identify the latent space underlying the system. Intuitively, the change introduced by the dynamical system is then negligible compared to the noise. In Table 1 (“large $\sigma$”), we show that recovery is possible for cases with small angles, both in the linear and non-linear case. While in some cases, this learning setup might be applicable in practice, it seems generally unrealistic to be able to perturb the system beyond the actual dynamics.
As we scale the dynamics to larger values (Figure 4, panel b and c), the estimation scheme breaks again.
However, this property offers an explanation for the success of existing contrastive estimation algorithms like CEBRA-time [] which successfully estimate dynamics in absence of a dynamics model.

#### Symmetric encoders cannot identify non-trivial dynamics.

In the more general case where the dynamics dominates the system behavior, the baseline cannot identify linear dynamics (or more complicated systems). In the general LDS and SLDS cases, the baseline fails to identify the ground truth dynamics (Table 1) as predicted by Corollary 1 (rows marked with ✗).
For identity dynamics, the baseline is able to identify the latents (99.56% $R^{2}$) but breaks as soon as linear dynamics are introduced (73.56% $R^{2}$).

### 6.2 Approximation of non-linear dynamics

Next, we study in more details how the DynCL can identify piecewise linear or non-linear latent dynamics using the $\nabla$-SLDS dynamics model.

#### Identification of switching dynamics

Switching dynamics are depicted in Fig. 4a for four different modes of the 10 degrees dataset.
DynCL obtains high $R^{2}$ for various choices of dynamics (Fig. 4b) and additionally identifies the correct mode sequence (Fig. 4c) for all noise levels and variants of the underlying dynamics. As we increase the rotation angle used to generate the matrices, the gap between baseline and our model increases substantially.

#### Non-linear dynamics

Figure 5 depicts the Lorenz system as an example of a non-linear dynamical system for different choices of algorithms. The ground truth dynamics vary in the ratio between $dt$/$\sigma$ and we show the full range in panels b/c.
When the noise dominates the dynamics (panel a), the baseline is able to estimate also the nonlinear dynamics accurately, with 99.7%.
However, as we move to lower noise cases (panel b), performance reduces to 41.0%.
Our switching dynamics model is able to estimate the system with high $R^{2}$ in both cases (94.14% and 94.08%).
However, note that in this non-linear case, we are primarily succeeding at estimating the latent space, the estimated dynamics model did not meaningfully outperform an identity model (Appendix ).

### 6.3 Ablation studies

For practitioners leveraging contrastive learning for statistical analysis, it is important to know the trade-offs in empirical performance in relation to various parameters. In real-world experiments, the most important factors are the size of the dataset, the trial-structure of the dataset, the latent dimensionality we can expect to recover, and the degree of non-linearity between latents and observables. We consider these factors of influence: As a reference, we use the SLDS system with a 6D latent space, 1M samples (1k trials $\times$ 1k samples), $L=4$ mixing layers, 10 degrees for the rotation matrices, 65,536 negative samples per batch; batch size 2,048, and learning rate $10^{-3}$.

#### Impact of dataset size (Fig. 6 a).

We keep the number of trials fixed to 1k.
As we vary the sample size per trial, $R^{2}$ degrades for smaller dataset, and for the given setting we need at least 100 points per trial to attain identifiability empirically.
We outperform the baseline model in all cases.

#### Impact of trials (Fig. 6 b).

We next simulate a fixed number of 1M datapoints, which we split into trials of varying length. We consider 1k, 10k, 100k, and 1M as trial lengths. Performance is stable for the different settings, even for cases with small trial length (and less observed switching points).
DynCL consistently outperforms the baseline algorithm and attains stable performance close to the theoretical maximum given by the ground-truth dynamics.

Figure: Figure 6: Variations and ablations for the SLDS. We compare the $\nabla$-SLDS model to the ground-truth switching dynamics (oracle) and a standard CL model without dynamics (baseline). All variations are with respect to the setting with 1M timesteps (1k trials $\times$ 1k samples), $L=4$ mixing layers, $d=6$ latent dimensionality, 5 modes, and $p=0.0001$ switching probability.
Refer to caption: x5.png

#### Impact of non-linear mixing (Fig. 6 c).

All main experiments have been conducted with $L=4$ mixing layers in the mixing function ${\bm{g}}$.
Performance of DynCL stays at the theoretical maximum as we increase the number of mixing layers. As we move beyond four layers, both oracle performance in $R^{2}$ and our model declines, hinting that either (1) more data or (2) a larger model is required to recover the dynamics successfully in these cases.

#### Impact of dimensionality (Fig. 6 d).

Increasing latent dimensionality does not meaningfully impact performance of our model. We found that for higher dimensions, it is crucial to use a large number of negative examples (65k) for successful training.

#### Number of modes for switching linear dynamics fitting (Fig. 6 d).

Increasing the number of modes in the dataset leads to more successful fitting of the $R^{2}$ for the baseline model, but to a decline in accuracy. This might be due to the increased variance: While this helps the model to identify the latent space (dynamics appear more like noise), it still fails to identify the underlying dynamics model, unlike DynCL which attains high $R^{2}$ and cluster accuracy throughout.

#### Robustness to changes in switching probability (Fig. 6 e).

Figure: Figure 7: Impact of modes for non-linear dynamics in the Lorenz system for different noise levels, averaged over all $dt$.
Refer to caption: x6.png

Finally, we vary the switching probability. Higher switching probability causes shorter modes, which are harder to fix by the $\nabla$-SLDS dynamics model. Our model obtains high empirical identifiability throughout the experiment, but the accuracy metric begins to decline when $p=0.1$ and $p=0.2$.

#### Number of modes for non-linear dynamics fitting (Fig. 7 ).

We study the effect of increasing the number of matrices in the parameter bank ${\bm{\mathsfit{W}}}$ in the $\nabla$-SLDS model. The figure depicts the impact of increasing the number of modes for DynCL on the non-linear Lorenz dataset. We observe that increasing modes to 200 improves performance, but eventually converges to a stable maximum for all noise levels.

## 7 Discussion

The DynCL framework is versatile and allows to study the performance of contrastive learning in conjunction with different dynamics models. By exploring different special cases (identity, linear, switching linear), our study can be regarded as a categorization of different forms of contrastive learning and makes predictions about their behavior in practice.

In comparison to contrastive predictive coding [CPC; ] or wav2vec [], DynCL generalizes the concept of training contrastive learning models with (explicit) dynamics models. CPC uses an RNN encoder followed by linear projection, while wav2vec leverages CNNs dynamics models and affine projections. Theorem 1 applies to both these models, and offers an explanation for their successful empirical performance.

For applications in scientific data analysis, CEBRA [] uses supervised or self-supervised contrastive learning, either with symmetric encoders or asymmetric encoder functions. While our results show that such an algorithm is able to identify dynamics for a sufficient amount of system noise, adding dynamics models is required as the system dynamics dominate. Hence, the DynCL approach with LDS or $\nabla$-SLDS dynamics generalises the self-supervised mode of CEBRA and makes it applicable for a broader class of problems.

Finally, there is a connection to the joint embedding predictive architecture [JEPA; , ]. The model setup of DynCL (in particular with the $\nabla$-SLDS dynamics model) can be regarded as a special case of JEPA, but with symmetric encoders to leverage distillation of the system dynamics into the predictor (the dynamics model). In contrast to JEPA, the use of symmetric encoders again requires use of a contrastive loss for avoiding collapse and, more importantly, serves as the foundation for our theoretical result.

A limitation of the present study is its exclusive focus on simulated data which clearly corroborates our theory but does not yet demonstrate real-world applicability. However, our simulated data bears the signatures of real-world datasets (multi-trial structures, varying degrees of dimensionality, number of modes, and different forms of dynamics). A challenge is the availability of real-world benchmark datasets for dynamics identification. We believe that rigorous evaluation of different estimation methods on such datasets will continue to show the promise of contrastive learning for dynamics identification. Integrating recent benchmarks like DynaDojo [] with realistic mixing functions (${\bm{g}}$) offers a promising direction for evaluating latent dynamics models.

## 8 Conclusion

We proposed a first identifiable, end-to-end, non-generative inference algorithm for latent switching dynamics along with an empirically successful parameterization of non-linear dynamics. Our results point towards the empirical effectiveness of contrastive learning across time-series, and back these empirical successes by theory. We show empirical identifiability with limited data for linear, switching linear and non-linear dynamics. Our results add to the understanding of SSL’s empirical success, will guide the design of future contrastive learning algorithms and most importantly, make SSL amenable for computational statistics and data analysis.

### Author Contributions

RGL and TS: Methodology, Software, Investigation, Writing–Editing.
StS: Conceptualization, Methodology, Formal Analysis, Writing–Original Draft and Writing–Editing.

### Acknowledgments

We thank Luisa Eck and Stephen Jiang for discussions on the theory, and Lilly May for input on paper figures.

## Appendix A Proof of the main result

We re-state the Theorem from the main paper, and provide a full proof below:

Theorem 1 (Contrastive estimation of non-linear dynamics).
Assume that

- •
A time-series dataset $\{{\bm{y}}_{t}\}_{t=1}^{T}$ is generated according to the ground-truth dynamical system in Eq. 5 with a bijective dynamics model ${\bm{f}}$ and an injective mixing function ${\bm{g}}$.
- •
The system input follows an iid normal distribution, $p({\bm{u}}_{t})=\mathcal{N}({\bm{u}}_{t}|0,{\bm{\Sigma}}_{u})$.
- •
The model $\psi$ is composed of an encoder ${\bm{h}}$, a dynamics model ${\hat{{\bm{f}}}}$, a correction term $\alpha$, and the similarity metric $\phi({\bm{u}},{\bm{v}})=-\|{\bm{u}}-{\bm{v}}\|^{2}$ and attains the global minimizer of Eq. 3.

Then, in the limit of $T\rightarrow\infty$ for any point ${\bm{x}}$ in the support of the data marginal distribution:

- (a)
The composition of mixing and de-mixing ${\bm{h}}({\bm{g}}({\bm{x}}))={\bm{L}}{\bm{x}}+{\bm{b}}$ is a bijective affine transform, and ${\bm{L}}={\bm{Q}}{\bm{\Sigma}}_{u}^{-1/2}$ with unknown orthogonal transform ${\bm{Q}}\in\mathbb{R}^{d\times d}$ and offset ${\bm{b}}\in\mathbb{R}^{d}$.
- (b)
The estimated dynamics $\hat{\bm{f}}$ are bijective and identify the true dynamics ${\bm{f}}$ up to the relation ${\hat{{\bm{f}}}}({\bm{x}})={\bm{L}}{\bm{f}}({\bm{L}}^{-1}({\bm{x}}-{\bm{b}}))+
{\bm{b}}$.

###### Proof.

Our proof proceeds in three steps:
First, we leverage existing theory [, ] to arrive at the minimizer of the contrastive loss, and relate the limited sample loss function to the asymptotic case.
Second, we derive the statement about achieving successful demixing in Theorem 1(a).
Finally, we derive the statement in Theorem 1(a) about structural identifiability of the dynamics model.

###### Proof.

#### Step 1: Minimizer of the InfoNCE loss.

By the assumption ${\bm{u}}_{t}$ is normally distributed, we obtain the positive sample conditional distribution

$$ $\displaystyle p({\bm{x}}_{t+1}|{\bm{x}}_{t})$ $\displaystyle=\mathcal{N}({\bm{x}}_{t+1}|{\bm{f}}({\bm{x}}_{t}),{\bm{\Sigma}}_ {u}).$ (16) The negative sample distribution $q({\bm{x}}_{t})$ is obtained by sampling $t$ uniformly from all time-steps and can hence be written as a Gaussian mixture along the dynamics imposed by ${\bm{f}}$ , $\displaystyle q({\bm{x}})$ $\displaystyle=\frac{1}{T}\sum_{t=1}^{T}p_{u}({\bm{x}}-{\bm{x}}_{t})$ (17) $\displaystyle=\frac{1}{T}\sum_{t=1}^{T}\mathcal{N}({\bm{x}}-{\bm{x}}_{t}|{\bm{ f}}({\bm{x}}_{t-1}),{\bm{\Sigma}}_{u}).$ (18) $$

We use these definitions of $p$ and $q$ to study the asymptotic case of our loss function. For $T\rightarrow\infty$, due to Wang & Isola [] we can rewrite the limit of our loss function (Eq. 3) as

$$ $\displaystyle\mathcal{L}[\psi]=$ $\displaystyle\lim_{T\rightarrow\infty}\mathbb{E}_{t,N}[-\log p_{\psi}({\bm{y}} _{t+1}|{\bm{y}}_{t},N)]-\log T$ (19) $\displaystyle=$ $\displaystyle\int q({\bm{y}})\left[-\int p({\bm{y}}^{\prime}|{\bm{y}})\psi({ \bm{y}},{\bm{y}}^{\prime})d{\bm{y}}^{\prime}+\log\int q({\bm{y}}^{\prime})\exp [\psi({\bm{y}},{\bm{y}}^{\prime})]d{\bm{y}}^{\prime}\right].$ $$

It was shown [Proposition 1, ] that this loss function is convex in $\psi$ with the unique minimizer

$$ $\displaystyle\psi({\bm{g}}({\bm{x}}),{\bm{g}}({\bm{x}}^{\prime}))$ $\displaystyle=\log\frac{p({\bm{x}}^{\prime}|{\bm{x}})}{q({\bm{x}}^{\prime})}+c ({\bm{x}}),$ (20) where $c:\mathbb{R}^{d}\mapsto\mathbb{R}$ is an arbitrary scalar-valued function. Note that we also expressed ${\bm{y}}={\bm{g}}({\bm{x}}),{\bm{y}}^{\prime}={\bm{g}}({\bm{x}}^{\prime})$ to continue the proof in terms of the relation between original and estimated latents. We insert the definition of the model on the left hand side. Let us also denote ${\bm{h}}\circ{\bm{g}}=:{\bm{r}}$ , and the definition of the ground-truth generating process on the right hand side to obtain $\displaystyle\phi(\hat{\bm{f}}({\bm{r}}({\bm{x}})),{\bm{r}}({\bm{x}}^{\prime}) )-\alpha({\bm{x}}^{\prime})$ $\displaystyle=\log p({\bm{x}}^{\prime}|{\bm{f}}({\bm{x}}))-\log q({\bm{x}}^{ \prime})+c({\bm{x}}).$ (21) Inserting the potential 2 2 2 This is feasible in practice by parameterizing $\alpha({\bm{x}})$ as a kernel density estimate, but empirically often not required. See Appendix B for additional technical details. as $\alpha({\bm{x}})=\log q({\bm{x}})$ simplifies the equation to $\displaystyle\phi(\hat{\bm{f}}({\bm{r}}({\bm{x}})),{\bm{r}}({\bm{x}}^{\prime}))$ $\displaystyle=\log p({\bm{x}}^{\prime}|{\bm{f}}({\bm{x}}))+c({\bm{x}}).$ (22) From here onwards, we will use that $\phi$ is the negative squared Euclidean norm and correspondingly, the positive conditional is a normal distribution with concentration ${\bm{\Lambda}}={\bm{\Sigma}}_{u}^{-1}$ , $\displaystyle-\|\hat{\bm{f}}({\bm{r}}({\bm{x}}))-{\bm{r}}({\bm{x}}^{\prime})\| _{2}^{2}$ $\displaystyle=-({\bm{f}}({\bm{x}})-{\bm{x}}^{\prime})^{\top}{\bm{\Lambda}}({ \bm{f}}({\bm{x}})-{\bm{x}}^{\prime})+c^{\prime}({\bm{x}})$ (23) $$

where we pulled the normalization constant of $p$ into the function $c^{\prime}$ for brevity.

#### Step 2: Properties of the feature encoder.

Starting from the last equation, we compute the derivative with respect to ${\bm{x}}$ and ${\bm{x}}^{\prime}$ on both sides and obtain

$$ $\displaystyle{\bm{J}}^{\top}_{\bm{r}}({\bm{x}}^{\prime}){\bm{J}}_{\hat{{\bm{f} }}}({\bm{r}}({\bm{x}})){\bm{J}}_{\bm{r}}({\bm{x}})$ $\displaystyle={\bm{\Lambda}}{\bm{J}}_{\bm{f}}({\bm{x}}).$ (24) Because this equation holds for any ${\bm{x}}^{\prime}\in\operatorname{supp}q$ independently of ${\bm{x}}$ , we can conclude that the Jacobian matrix of ${\bm{r}}$ needs to be constant. From there it follows that ${\bm{r}}$ is affine. Let us write ${\bm{r}}({\bm{x}})={\bm{L}}{\bm{x}}+{\bm{b}}$ . Then, the Jacobian matrix ${\bm{J}}_{\bm{r}}={\bm{L}}$ . Inserting this yields $\displaystyle{\bm{L}}^{\top}{\bm{J}}_{\hat{{\bm{f}}}}({\bm{L}}{\bm{x}}+{\bm{b} }){\bm{L}}$ $\displaystyle={\bm{\Lambda}}{\bm{J}}_{\bm{f}}({\bm{x}}).$ (25) $$

We next establish that ${\bm{L}}$ has full rank: because the dynamics function ${\bm{f}}$ is bijective by assumption, ${\bm{J}}_{\bm{f}}({\bm{x}})$ has full rank $d$. ${\bm{\Lambda}}$ has full rank by assumption about the distribution $p({\bm{u}})$. All matrices on the LHS are square and need to have full rank as well for any point ${\bm{x}}$. Hence, we can conclude that ${\bm{L}}$ has full rank, and likewise ${\bm{J}}_{{\hat{{\bm{f}}}}}$. From there, we conclude that ${\bm{r}}$ and ${\hat{{\bm{f}}}}$ are bijective.

Next, we derive additional constraints on the matrix ${\bm{L}}$.
We insert the solution for ${\bm{r}}$ obtained so far in Eq. 23,

$$ $\displaystyle-\|\hat{\bm{f}}({\bm{L}}{\bm{x}}+{\bm{b}})-{\bm{L}}{\bm{x}}^{ \prime}-{\bm{b}}\|^{2}$ $\displaystyle=-({\bm{f}}({\bm{x}})-{\bm{x}}^{\prime})^{\top}{\bm{\Lambda}}({ \bm{f}}({\bm{x}})-{\bm{x}}^{\prime})+c^{\prime}({\bm{x}})$ (26) and take the derivative twice with respect to ${\bm{x}}^{\prime}$ , to obtain $\displaystyle{\bm{L}}^{\top}{\bm{L}}$ $\displaystyle={\bm{\Lambda}}\quad\Leftrightarrow\quad{\bm{L}}^{\top}={\bm{ \Lambda}}{\bm{L}}^{-1}.$ (27) $$

Without loss of generality, we introduce ${\bm{Q}}\in\mathbb{R}^{d\times d}$ to write ${\bm{L}}$ in terms of ${\bm{\Lambda}}$ as ${\bm{L}}={\bm{Q}}{\bm{\Lambda}}^{1/2}$. Inserting into the previous equation lets us conclude

$$ $\displaystyle{\bm{L}}^{\top}{\bm{L}}$ $\displaystyle={\bm{\Lambda}}^{1/2}{\bm{Q}}^{\top}{\bm{Q}}{\bm{\Lambda}}^{1/2}= {\bm{\Lambda}}$ (28) $\displaystyle{\bm{Q}}^{\top}{\bm{Q}}$ $\displaystyle={\bm{\Lambda}}^{-1/2}{\bm{\Lambda}}{\bm{\Lambda}}^{-1/2}={\bm{I}},$ (29) $$

from which follows that ${\bm{Q}}$ is an orthogonal matrix. Hence, ${\bm{L}}$ is a composition of an orthogonal transform and ${\bm{\Sigma}}_{u}^{-1/2}$, concluding the first part of the proof for statement (a).

#### Step 3: Dynamics.

To derive part (b), we start at the condition Eq. 23 again to determine the value of $c^{\prime}({\bm{x}})$. We can consider two special cases:

$$ $\displaystyle{\bm{f}}({\bm{x}})={\bm{x}}^{\prime}:\quad c^{\prime}({\bm{x}})$ $\displaystyle=-\|\hat{\bm{f}}({\bm{r}}({\bm{x}}))-{\bm{r}}({\bm{x}}^{\prime}) \|^{2}\leq 0$ (30) $\displaystyle\hat{\bm{f}}({\bm{r}}({\bm{x}}))={\bm{r}}({\bm{x}}^{\prime}): \quad c^{\prime}({\bm{x}})$ $\displaystyle=({\bm{f}}({\bm{x}})-{\bm{x}}^{\prime})^{\top}{\bm{\Lambda}}({\bm {f}}({\bm{x}})-{\bm{x}}^{\prime})\geq 0$ (31) $$

When combining both conditions for points where ${\bm{f}}({\bm{x}})={\bm{x}}^{\prime}$ and $\hat{\bm{f}}({\bm{r}}({\bm{x}}))={\bm{r}}({\bm{x}}^{\prime})$ the only admissible solution is $c^{\prime}({\bm{x}})=0$ for points with $\hat{\bm{f}}({\bm{r}}({\bm{x}}))={\bm{r}}({\bm{f}}({\bm{x}}))$, i.e. $\hat{\bm{f}}({\bm{x}})={\bm{r}}({\bm{f}}({\bm{r}}^{-1}({\bm{x}})))$, hinting at the final solution.
However, we have not shown yet that this solution is unique.

To show uniqueness, without loss of generality, we use the ansatz

$$ $\displaystyle{\bm{f}}({\bm{x}})$ $\displaystyle={\bm{A}}_{1}{\hat{{\bm{f}}}}({\bm{L}}{\bm{x}}+{\bm{b}})+{\bm{d}} _{1}+{\bm{\varepsilon}}({\bm{x}}),$ (32) and computing the derivative with respect to ${\bm{x}}$ yields $\displaystyle{\bm{J}}_{\bm{f}}({\bm{x}})$ $\displaystyle={\bm{A}}_{1}{\bm{J}}_{\hat{{\bm{f}}}}({\bm{L}}{\bm{x}}+{\bm{b}}) {\bm{L}}+{\bm{J}}_{\bm{\varepsilon}}({\bm{x}}).$ (33) We insert this into Eq. 25 and obtain $\displaystyle{\bm{L}}^{\top}{\bm{J}}_{\hat{{\bm{f}}}}({\bm{L}}{\bm{x}}+{\bm{b} }){\bm{L}}$ $\displaystyle={\bm{\Lambda}}{\bm{A}}_{1}{\bm{J}}_{\hat{{\bm{f}}}}({\bm{L}}{\bm {x}}+{\bm{b}}){\bm{L}}+{\bm{\Lambda}}{\bm{J}}_{\bm{\varepsilon}}({\bm{x}})$ (34) $\displaystyle({\bm{L}}^{\top}-{\bm{\Lambda}}{\bm{A}}_{1}){\bm{J}}_{\hat{{\bm{f }}}}({\bm{L}}{\bm{x}}+{\bm{b}}){\bm{L}}$ $\displaystyle={\bm{\Lambda}}{\bm{J}}_{\bm{\varepsilon}}({\bm{x}})$ (35) $\displaystyle({\bm{L}}^{\top}-{\bm{\Lambda}}{\bm{A}}_{1})$ $\displaystyle={\bm{\Lambda}}{\bm{J}}_{\bm{\varepsilon}}({\bm{x}}){\bm{L}}^{-1} {\bm{J}}^{-1}_{\hat{{\bm{f}}}}({\bm{L}}{\bm{x}}+{\bm{b}})$ (36) The left hand side is a constant, hence the same needs to hold true for the right hand side. Without loss of generality, let us introduce an arbitrary matrix ${\bm{A}}_{2}$ we set as this constant, $\displaystyle{\bm{J}}_{\bm{\varepsilon}}({\bm{x}}){\bm{L}}^{-1}{\bm{J}}^{-1}_{ \hat{{\bm{f}}}}({\bm{L}}{\bm{x}}+{\bm{b}})$ $\displaystyle={\bm{A}}_{2}$ (38) $\displaystyle{\bm{J}}_{\bm{\varepsilon}}({\bm{x}})$ $\displaystyle={\bm{A}}_{2}{\bm{J}}_{\hat{{\bm{f}}}}({\bm{L}}{\bm{x}}+{\bm{b}}) {\bm{L}}$ (39) which only admits the solution $\displaystyle{\bm{\varepsilon}}({\bm{x}})$ $\displaystyle={\bm{A}}_{2}{\hat{{\bm{f}}}}({\bm{L}}{\bm{x}}+{\bm{b}})+{\bm{d}} _{2},$ (40) where we introduced an additional integration constant ${\bm{d}}_{2}$ . Inserting this into the ansatz in Eq. 32 gives $\displaystyle{\bm{f}}({\bm{x}})$ $\displaystyle={\bm{A}}_{1}{\hat{{\bm{f}}}}({\bm{L}}{\bm{x}}+{\bm{b}})+{\bm{d}} _{1}+{\bm{\varepsilon}}({\bm{x}}),$ (41) $\displaystyle{\bm{f}}({\bm{x}})$ $\displaystyle=({\bm{A}}_{1}+{\bm{A}}_{2}){\hat{{\bm{f}}}}({\bm{L}}{\bm{x}}+{ \bm{b}})+({\bm{d}}_{1}+{\bm{d}}_{2}).$ (42) Using the shorthand ${\bm{A}}={\bm{A}}_{1}+{\bm{A}}_{2}$ , ${\bm{d}}={\bm{d}}_{1}+{\bm{d}}_{2}$ we can repeat the steps in Eqs. 33 – LABEL:eq:ansatz-v1-last-step to arrive at the condition $\displaystyle({\bm{L}}^{\top}-{\bm{\Lambda}}{\bm{A}}){\bm{J}}_{\hat{{\bm{f}}}} ({\bm{L}}{\bm{x}}+{\bm{b}}){\bm{L}}$ $\displaystyle=0.$ (43) Since all matrices have full rank, the only valid solution is ${\bm{A}}={\bm{\Lambda}}^{-1}{\bm{L}}^{\top}={\bm{L}}^{-1}$ . Inserting back into the ansatz yields the refined solution $\displaystyle{\bm{f}}({\bm{x}})$ $\displaystyle={\bm{L}}^{-1}{\hat{{\bm{f}}}}({\bm{L}}{\bm{x}}+{\bm{b}})+{\bm{d}},$ (44) and for brevity, we let ${\bm{\xi}}={\hat{{\bm{f}}}}({\bm{r}}({\bm{x}}))={\hat{{\bm{f}}}}({\bm{L}}{\bm{ x}}+{\bm{b}})$ : $\displaystyle{\bm{f}}({\bm{x}})$ $\displaystyle={\bm{L}}^{-1}{\bm{\xi}}+{\bm{d}}.$ (45) $$

We then insert the current solution into Eq. 23 and input ${\bm{r}}$ which gives

$$ $\displaystyle\|{\bm{\xi}}-{\bm{L}}{\bm{x}}^{\prime}-{\bm{b}}\|^{2}$ $\displaystyle=({\bm{L}}^{-1}{\bm{\xi}}+{\bm{d}}-{\bm{x}}^{\prime})^{\top}{\bm{ \Lambda}}({\bm{L}}^{-1}{\bm{\xi}}+{\bm{d}}-{\bm{x}}^{\prime})+c^{\prime}({\bm{ x}})$ (46) $\displaystyle=({\bm{L}}^{-1}{\bm{\xi}}+{\bm{d}}-{\bm{x}}^{\prime})^{\top}{\bm{ L}}^{\top}{\bm{L}}({\bm{L}}^{-1}{\bm{\xi}}+{\bm{d}}-{\bm{x}}^{\prime})+c^{ \prime}({\bm{x}})$ (47) $\displaystyle=\|{\bm{\xi}}+{\bm{L}}{\bm{d}}-{\bm{L}}{\bm{x}}^{\prime}\|^{2}+c^ {\prime}({\bm{x}})$ (48) $\displaystyle c^{\prime}({\bm{x}})$ $\displaystyle=\|{\bm{\xi}}-{\bm{L}}{\bm{x}}^{\prime}-{\bm{b}}\|^{2}-\|{\bm{\xi }}-{\bm{L}}{\bm{x}}^{\prime}+{\bm{L}}{\bm{d}}\|^{2}$ (49) Let us denote ${\bm{v}}={\bm{\xi}}-{\bm{L}}{\bm{x}}^{\prime}$ and note that ${\bm{v}}$ and ${\bm{x}}$ remain independent variables. We then get $\displaystyle c^{\prime}({\bm{x}})$ $\displaystyle=\|{\bm{v}}-{\bm{b}}\|^{2}-\|{\bm{v}}+{\bm{L}}{\bm{d}}\|^{2}$ (50) $\displaystyle=-2{\bm{v}}^{\top}({\bm{b}}+{\bm{L}}{\bm{d}})+\|{\bm{b}}\|^{2}-\| {\bm{L}}{\bm{d}}\|^{2}$ (51) Because ${\bm{v}}$ and ${\bm{x}}$ vary independently and the equation is true for any pair of these points, both sides of the equation need to be independent of their respective variables. This is true only if ${\bm{b}}=-{\bm{L}}{\bm{d}}$ . Hence, it follows that $\displaystyle c^{\prime}({\bm{x}})$ $\displaystyle=0\quad\text{and}\quad{\bm{d}}=-{\bm{L}}^{-1}{\bm{b}}.$ (52) $$

Inserting ${\bm{d}}$ into Eq. 44 gives the final solution,

$$ $\displaystyle{\bm{f}}({\bm{x}})$ $\displaystyle={\bm{L}}^{-1}{\hat{{\bm{f}}}}({\bm{L}}{\bm{x}}+{\bm{b}})-{\bm{L} }^{-1}{\bm{b}}.$ (53) Solving for ${\hat{{\bm{f}}}}$ gives us $\displaystyle{\hat{{\bm{f}}}}({\bm{L}}{\bm{x}}+{\bm{b}})$ $\displaystyle={\bm{L}}{\bm{f}}({\bm{x}})+{\bm{b}}$ (54) $\displaystyle{\hat{{\bm{f}}}}({\bm{x}})$ $\displaystyle={\bm{L}}{\bm{f}}({\bm{L}}^{-1}({\bm{x}}-{\bm{b}}))+{\bm{b}}=({ \bm{r}}\circ{\bm{f}}\circ{\bm{r}}^{-1})({\bm{x}})$ (55) $$

which concludes the proof.
∎

## Appendix B Kernel density estimate correction

Theorem 1 requires to include a “potential function” $\alpha$ into our model. In this section, we discuss how this function can be approximated by a kernel density estimate (KDE) in practice. The KDE intuitively corrects for the case of non-uniform marginal distributions.
Correcting with $\alpha$ overcomes the limitation of requiring a uniform marginal distribution discussed before []. While other solutions have been discussed, such as training a separate MLP [], the KDE solution discussed below is conceptually simpler and non-parametric.

Figure: Figure 8: Introducing KDE into the loss allows to compensate for non-uniform marginal distribution. We show performance in terms of $R^{2}$ across datasets with increasingly non-uniform marginal. We replicate the data-generating process and experimental setup performed by Zimmermann et al. [, Figure 2].
Refer to caption: x7.png

For the models considered in the main paper, we considered representation learning in Euclidean space, while Appendix D contains some additional experiments for the very common case of training embeddings on the hypersphere. For both cases, we can parameterize appropriate KDEs.

For the Euclidean case, we use the KDE based on the squared Euclidean norm,

$$ $\hat{q}({\bm{x}})=\frac{1}{{\epsilon}M}\sum_{i=1}^{M}\exp\left(-\frac{\|{\bm{x }}-{\bm{x}}_{i}\|^{2}}{{\epsilon}}\right),\quad{\bm{x}}_{i}\sim q({\bm{x}}).$ (56) $$

We note that in the limit ${\epsilon}\rightarrow 0$, $M\rightarrow\infty$, this estimate converges to the correct distribution, $\hat{q}({\bm{x}})\rightarrow q({\bm{x}})$. This is also the case used in Theorem 1.
However, this estimate depends on the ground truth latents ${\bm{x}}_{i}$, which are not accessible during training. Hence, we need to find an expression that depends on the observable data. We leverage the feature encoder ${\bm{h}}$ to express the estimator as

$$ $\hat{q}_{\bm{h}}({\bm{y}})=\frac{1}{{\epsilon}M}\sum_{i=1}^{M}\exp\left(-\frac {\|{\bm{h}}({\bm{y}})-{\bm{h}}({\bm{y}}_{i})\|^{2}}{{\epsilon}}\right),\quad{ \bm{y}}_{i}\sim q({\bm{y}}).$ (57) $$

We can express this estimator in terms of the final solution, ${\bm{r}}({\bm{x}})={\bm{h}}({\bm{g}}({\bm{x}}))={\bm{Q}}{\bm{\Sigma}}_{u}^{-1/
2}{\bm{x}}+{\bm{b}}$ in the theorem. If we express the solution in terms of the ground truth latents again, the orthogonal matrix ${\bm{Q}}$ vanishes and we obtain

$$ $\hat{q}_{\bm{h}}({\bm{x}})=\frac{1}{{\epsilon}M}\sum_{i=1}^{M}\exp\left(-\frac {\|{\bm{\Sigma}}^{-1/2}_{u}({\bm{x}}-{\bm{x}}_{i})\|^{2}}{{\epsilon}}\right). \quad{\bm{y}}_{i}\sim q({\bm{y}})$ (58) $$

This corresponds to a KDE using a Mahalanobis distance with covariance matrix ${\bm{\Sigma}}_{u}$, which is a valid KDE of $q$.

We can derive a similar argument when computing embeddings and dynamics on the hypersphere. When, a von Mises-Fisher distribution is suitable to express the KDE, and we obtain

$$ $\hat{q}({\bm{x}})=\frac{C_{p}(\kappa)}{M}\sum_{i=1}^{M}\exp(\kappa{\bm{x}}^{ \top}{\bm{x}}_{i}),\quad{\bm{x}}_{i}\sim q({\bm{x}})$ (59) $$

where $C_{p}(\kappa)$ is the normalization constant of the von Mises-Fisher distribution.
This again approaches the correct data distribution for $\hat{q}\rightarrow q$ as $M,\kappa\rightarrow\infty$. Following the same arguments above, but using
${\bm{r}}({\bm{x}})={\bm{h}}({\bm{g}}({\bm{x}}))={\bm{Q}}{\bm{x}}$ as the indeterminacy on the hypersphere, we can express this in terms of the ground truth latents,

$$ $\displaystyle\hat{q}_{\bm{h}}({\bm{x}})$ $\displaystyle=\frac{C_{p}(\kappa)}{N}\sum_{i=1}^{M}\exp\left(\kappa{\bm{r}}({ \bm{x}})^{\top}{\bm{r}}({\bm{x}}_{i})\right),$ (60) $\displaystyle=\frac{C_{p}(\kappa)}{N}\sum_{i=1}^{M}\exp\left(\kappa{\bm{x}}^{ \top}{\bm{x}}_{i}\right),$ (61) $$

which is again a valid KDE.

It is interesting to consider the effect of the KDE on the loss function.
Inserting $\psi\leftarrow\psi-\log\hat{q}$ into the loss function yields

$$ $\displaystyle-\log p_{\psi}({\bm{x}}|{\bm{x}}^{+},N)$ $\displaystyle=-(\psi({\bm{x}}_{i},{\bm{x}}^{+}_{i})-\log\hat{q}({\bm{x}}_{i}^{ +}))+\log\sum_{i=1}^{N}e^{\psi({\bm{x}}_{i},{\bm{x}}_{j}^{-})-\log\hat{q}({\bm {x}}_{j}^{-})},$ (62) $\displaystyle=-\psi({\bm{x}}_{i},{\bm{x}}^{+}_{i})+\log\hat{q}({\bm{x}}_{i}^{+ })+\log\sum_{i=1}^{N}\frac{1}{\hat{q}({\bm{x}}_{j}^{-})}e^{\psi({\bm{x}}_{i},{ \bm{x}}_{j}^{-})},$ (63) $\displaystyle=-\psi({\bm{x}}_{i},{\bm{x}}^{+}_{i})+\log\sum_{i=1}^{N}\frac{ \hat{q}_{h}({\bm{x}}_{i}^{+})}{\hat{q}_{h}({\bm{x}}_{j}^{-})}e^{\psi({\bm{x}}_ {i},{\bm{x}}_{j}^{-})},$ (64) $\displaystyle=-\psi({\bm{x}}_{i},{\bm{x}}^{+}_{i})+\log\sum_{i=1}^{N}w_{h}({ \bm{x}}_{i}^{+},{\bm{x}}_{j}^{-})e^{\psi({\bm{x}}_{i},{\bm{x}}_{j}^{-})}$ (65) $$

with the importance weights $w_{h}({\bm{x}}_{i}^{+},{\bm{x}}_{j}^{-})=\frac{\hat{q}_{h}({\bm{x}}_{i}^{+})}{
\hat{q}_{h}({\bm{x}}_{j}^{-})}$. Intuitively, this means that the negative examples are re-weighted according to the density ratio between the current positive and each negative sample.

#### Empirical motivation

Figure 8 shows preliminary results on applying this KDE correction to contrastive learning models. We followed the setting from Zimmermann et al. [] and re-produced the experiment reported in Fig. 2 in their paper. We use 3D latents, a 4-layer MLP as non-linear mixing function with a final projection layer to 50D observed data. The reference, positive and negative distributions are all vMFs parameterized according to $\kappa$ (x-axis) in the case of the reference and negative distribution and $\kappa_{p}$ for the positive distribution.

The grey curve shows the decline in empirical identifiability ($R^{2}$) as the uniformity assumption is violated by an increasing concentration $\kappa$ (x-axis). Applying a KDE correction to the data resulted in substantially improved performance (red lines).

However, when testing the method directly on the dynamical systems considered in the paper, we did not found a substantial improvement in performance. One hypothesis for this is that the distribution of points on the data manifold (not necessarily the whole $\mathbb{R}^{d}$ is already sufficiently uniform.
Hence, while the theory requires inclusion of the KDE term (and it did not degrade results), we suggest to drop this computationally expensive term when applying the method on real-world datasets that are approximately uniform.

## Appendix C Additional Related Work

#### Contrastive learning.

An influential and conceptual motivation for our work is Contrastive Predictive Coding (CPC) [] which uses the InfoNCE loss with an additional non-linear projection head implemented as an RNN to aggregate information from multiple timesteps. Then, an affine projection is used for multiple forward prediction steps. However, here the “dynamics model” is not explicitly parameterized, limiting its interpretability.
Similar frameworks have been successfully applied across various domains, including audio, vision, and language, giving rise to applications such as wav2vec [], time contrastive networks for video [TCN; ] or CPCv2 [].

#### Non-Contrastive learning.

data2vec [] and JEPA [] learns a representation by trying to predict missing information in latent space, using an MSE loss. JEPA uses asymmetric encoders, and on top a predictor model in latent space parameterized by a neural net.

#### System identification.

In system identification, a problem closely related to the one addressed in this work is known as ”nonlinear system identification. Widely used algorithms for this problem include Extended Kalman Filter (EKM) [] and Nonlinear Autoregressive Moving Average with Exogenous inputs (NARMAX) []. EKF is based on linearizing ${\bm{g}}$ and ${\bm{f}}$ using a first-order Taylor-series approximation and then apply the Kalman Filter (KF) to the linearized functions. NARMAX, on the other hand, typically employs a power-form polynomial representation to model the non-linearities.
In neuroscience, practical (generative algorithms) include systems modeling linear dynamics [fLDS; ] or non-linear dynamics modelled by RNNs [LFADS; ]. Hurwitz et al. [] provide a detailed summary of additional algorithms.

#### Nonlinear ICA.

The field of Nonlinear ICA has also extensively studied a similar problem. For example, Time Contrastive Learning (TCL) [] uses a contrastive loss to predict the segment-ID of multivariate time-series which was shown to perform Non-linear ICA. Permutation Contrastive Learning (PCL) [] permutes the time series and aims to distinguish positive and negatives pair.

#### Switching Linear Dynamical Systems.

Several foundational papers have addressed this problem [, , ], , leading to a variety of extensions and variants. For example, Recurrent SLDSs [, ] address state-dependent switching by changing the switch transition distribution to $p(y_{t}|y_{t-1},x_{t-1})$, allowing for more flexible dependencies on previous states. Another extension, Explicit duration SLDS introduces additional latent variables to model the distribution of switch durations explicitly []. Some approaches relax the assumption of linear dynamics, such as in the case of SNLDS and RSSSM, where the dynamics model is assumed to be nonlinear [, ]. In the context of Nonlinear Independent Component Analysis (ICA), recent extensions include structured data generating processes (e.g., SNICA []) which were shown to be useful for the inference of switching dynamics. In this vein, [] proposed additional identifiability theory for the switching case. Other approaches, based on Neural Ordinary Differential Equations (Neural ODEs) [, ], or methods aimed at discovering switching dynamics within recurrent neural networks (RNNs) [], also present interesting avenues for modeling switching dynamics.

#### Deep state-space models.

Recently, (deep) state-space models (SSMs) such as S4 or Mamba [, ] have emerged as a promising architecture. These models are particularly well-suited for capturing long-range dependencies, making them an attractive choice for sequence modeling tasks.

#### Symbolic Regression.

An alternative approach to modeling dynamical systems is the use of symbolic regression, which aims to directly infer explicit symbolic mathematical expressions governing the underlying dynamical laws. Examples include Sparse Identification of Nonlinear Dynamics [SINdy,[]], as well as more recent transformer-based models [, ], which have demonstrated promise in discovering interpretable representations of dynamical systems.

## Appendix D Von Mises–Fisher (vMF) conditional distributions

In the main paper, we have shown experimental results that verify Theorem 1 in the case of Normal distributed positive conditional distribution and using the Euclidean distance. This approach has allowed for modeling latents and their dynamics in Euclidean space, which we argue is the most practical setting to apply DynCL in. However, self-supervised learning methods and especially contrastive learning have commonly been applied to produce representations on the hypersphere and using the dot-product distance [, , , , ].

Here we validate empirically that Theorem 1 equally holds under the assumption of vMF conditional distributions and using the dot-product distance $\phi({\bm{x}},{\bm{y}})={\bm{x}}^{\top}{\bm{y}}$ as part of the loss. We run experiments as in Table 1 for the case where the true dynamics model ${\bm{f}}$ is a linear dynamical system.
Additionally, we vary the setting similar to Figure 4 to show increasing $\Delta t$ (angular velocity).

We compare:

- •
DynCL (ours) – with linear dynamics:
${\hat{{\bm{f}}}}({\bm{x}})=\hat{\bm{A}}{\bm{x}}$.
- •
GTD – the ground-truth dynamics model (LDS)
${\hat{{\bm{f}}}}({\bm{x}})={\bm{A}}{\bm{x}}$.
- •
No dynamics – the baseline setting we use throughout the paper
${\hat{{\bm{f}}}}({\bm{x}})={\bm{x}}$.
- •
Asymmetric – a variation on the baseline setting that uses asymmetric encoders (one for reference, one for positive or negative) which would be a possible fix of Corollary 1. We can obtain this setting by skipping the explicit dynamics modeling, and defining two encoders ${\bm{h}}_{1},{\bm{h}}_{2}$ which relate as follows:
${\bm{h}}\circ{\bm{f}}:={\bm{h}}_{1}$,
${\bm{h}}:={\bm{h}}_{2}$.

Figure: Figure 9: Our findings from Table 1 hold equally under vMF noise distribution when using LDS ground truth dynamics. We show empirical identifiability of the latents in terms of $R^{2}$ under varying a) angles of the rotation dynamics i.e. angular velocity $\Delta t$ (x-axis) and b) the magnitude of the dynamics noise $\sigma$ (panels, left: low noise, right: high noise)
Refer to caption: x8.png

Similar to our results for the Euclidean case (Table 1), in Figure 9 we show results that experimentally verify Theorem 1 for latent dynamics on hypersphere and using vMF as conditional distribution. Both for low (left panel) and high (right panel) variance of the conditional distribution we can see that DynCL effectively identifies the ground truth latents on par with the oracle (GTD) model performance. On the other hand, the baseline, standard time contrastive learning without dynamics, can not identify the ground truth latents with underlying linear dynamics as predicted by Corollary 1. This prediction is only violated in the case where the variance of the noise distribution is high enough, such that the noise dominates the changes introduced by the actual dynamics. This is the case for dynamics with rotations up to 4 degrees for $\sigma=0.01$ and angles up to 10 degrees for $\sigma=0.1$.

## Appendix E Additional plots for SLDS

Figure: Figure 10: Visualizations of 6D linear dynamical systems at $\sigma=0.0001$ (left) and $\sigma=0.01$ for 10 degree rotations. These systems are used in our SLDS experiments.

## Appendix F Generalization - Train- vs Testset

In the main paper, all metrics are evaluated using the full training dataset of the respective experiment. We argue that this is sufficient for showing the efficacy of our model and verifying claims from the theory in section 3 because a) in self-supervised learning, the model learns generalizable representations through pretext tasks, making overfitting less of a concern; b) the metrics we are interested in are about uncovering the true underlying latent representation and dynamics of the available training data, not of new data; and c) most importantly, we ensured that the training dataset is large enough to approximate the full data distribution.

Nonetheless, here we show a series of control experiments to re-evaluate models on new and unseen data. We do so by following the same data generating process of the given experiment (same dynamics model and mixing function) and sample completely new trials (10% of the number of trials of the training dataset). Every new trial starts at a random new starting point, with a randomly sampled new mode sequence and regenerated with different seeds for the dynamics noise.

First, we re-evaluated every experiment of Figure 6 on the test dataset generated as described above and show these results in Figure 11. Comparing those results to the train dataset version of Figure 6 shows that there is almost no difference in the performance (with regard to identifiability and systems identification). The difference are so small, that visually comparing the results almost becomes impractical, so we additionally provide the exact numbers of the first panel (variations on the number of samples per trial) in Table 2.

Finally, to qualitatively and quantitatively show the difference between the train and test datasets we provide a) depictions of the ground truth latents of 5 random test trials and their closest possible matching trial from the training set in Figure 12 and b) a distribution of the distances (in terms of $R^{2}$ between the data from the test and train trials) between all test trials of one of a random test set and their closest trial from the training dataset in Figure 13.

Figure: Figure 11: Same as Figure 6 (Variations and ablations for SLDS), but re-evaluated on a newly generated test data with different starting points.
Refer to caption: x10.png

**Table 2: A detailed view on the #samples panel from Figure 6 and 11 showing the difference between train and test set evaluation.**
|  | DynCL (ours) | CL w/o dynamics | CL w/ ground truth dynamics |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
|  | $R^{2}$ (train) | $R^{2}$ (test) | $R^{2}$ (train) | $R^{2}$ (test) | $R^{2}$ (train) | $R^{2}$ (test) |
| # samples |  |  |  |  |  |  |
| 10 | 0.991 $\pm$ 0.00137 | 0.989 $\pm$ 0.00172 | 0.923 $\pm$ 0.03703 | 0.917 $\pm$ 0.04000 | 0.954 $\pm$ 0.00397 | 0.952 $\pm$ 0.00442 |
| 100 | 0.995 $\pm$ 0.00106 | 0.994 $\pm$ 0.00108 | 0.786 $\pm$ 0.06794 | 0.791 $\pm$ 0.06931 | 0.990 $\pm$ 0.00107 | 0.990 $\pm$ 0.00099 |
| 1000 | 0.995 $\pm$ 0.00074 | 0.995 $\pm$ 0.00078 | 0.765 $\pm$ 0.06671 | 0.770 $\pm$ 0.06973 | 0.990 $\pm$ 0.00507 | 0.990 $\pm$ 0.00511 |
| 10000 | 0.996 $\pm$ 0.00046 | 0.996 $\pm$ 0.00044 | 0.694 $\pm$ 0.06937 | 0.801 $\pm$ 0.10568 | 0.991 $\pm$ 0.00509 | 0.991 $\pm$ 0.00425 |

Figure: Figure 12: Ground truth latents from five random trials of the testsets for Figure 11 and their closest match within the corresponding trainset. The closest match is evaluated by computing the $R^{2}$-Score between a given trial from the testset and every possible trial.

Figure: Figure 13: Histogram of all $R^{2}$-Scores between every trial from the testset and its closest possible match from the trainset as shown in Figure 12.

## Appendix G Variations and additional baselines for the Dyn R 2 superscript 𝑅 2 R^{2} italic_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT Metric

#### Method.

As an addition to Table 1, we analyse the Dyn$R^{2}$ in more detail. In Table 3 we show variants for the metric. Firstly, we modify the number of forward prediction steps,

$$ $\displaystyle{\bm{f}}^{n}({\bm{x}}):=(\underbrace{{\bm{f}}\circ\cdots\circ{\bm {f}}}_{n\text{ times}})({\bm{x}})$ (66) $$

and respectively for ${\hat{{\bm{f}}}}^{n}$ in relation to ${\hat{{\bm{f}}}}$.
We then consider two variants of Eq. 15. Firstly, we perform multiple forward predictions ($n>1$) and compare the resulting embeddings:

$$ $\text{r2_score}(\hat{\bm{f}}^{n}(\hat{\bm{x}}),{\bm{L}}_{1}{\bm{f}}^{n}({\bm{ L}}_{2}\hat{\bm{x}}+{\bm{b}}_{2})+{\bm{b}}_{1}).$ (67) $$

A rationale for this metric is that the prediction task becomes increasingly difficult with an increasing number of time steps, and errors accumulate faster.

Secondly, as an additional control, we replace ${\hat{{\bm{f}}}}$ with the identity, and compute

$$ $\text{r2_score}(\hat{\bm{x}},{\bm{L}}_{1}{\bm{f}}^{n}({\bm{L}}_{2}\hat{\bm{x} }+{\bm{b}}_{2})+{\bm{b}}_{1}).$ (68) $$

This metric can be regarded as a naiive baseline/control for comparing performance of the dynamics model. If the dyn$R^{2}$ is not significantly larger than this value, we cannot conclude to have obtained meaningful dynamics.

For the lower part of Table 1, we report the resulting metrics in Table 3, setting the number of forward steps $n$ to $1$ or $10$, and using either the original metric (Eq. 66), or the control (Eq. 67).

#### Results.

For the SLDS system, we can corroborate our results further: the baseline model obtains a dyn$R^{2}$ of around 85% for single step prediction, both for the original and control metric. Our $\nabla$-SLDS model and the ground truth dynamical model obtain over 99.9% well above the level of the control metric which remains at around 95%. The high value of the control metric is due to the small change introduced by a single timestep, and should be considered when using and interpreting the metric. If more steps are performed, the performance of the $\nabla$-SLDS model drops to about 95.5% vs. chance level for the control metric, again highlighting the high performance of our model, but also the room for improvement, as the oracle model stays at above 99% as expected.

For the Lorenz system, we do not see a significant difference between original dyn$R^{2}$metric and dyn$R^{2}$ control for any of the considered algorithms. Yet, as noted in the main paper, $\nabla$-SLDS is the only dynamics model that gets a high $R^{2}$ of 94.08%, vs. the lower 81.20% for a single LDS model, or 40.99% for the baseline model.
In other words, while DynCL with the $\nabla$-SLDS dynamics model falls short of identifying the true underlying dynamics for this non-linear chaotic system, without DynCL we wouldn’t even identify the latents. We leave optimizing the parameterization of the dynamics model to identify non-linear chaotic systems for future work.

**Table 3: Extended metrics for dynamics models including additional variations on the $dynR^{2}$ metric where Control is replacing $\hat{\bm{f}}$ with the identity and 10 Steps is applying $\hat{\bm{f}}$ (and ${\bm{f}}$) 10 times, i.e., predicting 10 steps forward instead of only one step as is done in the Original version.**
| data |  | model | % dyn$R^{2}$, n=1 step | % dyn$R^{2}$, n=10 steps |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| ${\bm{f}}$ | $p({\bm{u}})$ |  | $\hat{\bm{f}}$ | theory | Original | Control | Original | Control |
| SLDS | Normal |  | identity | ✗ | 85.47 $\pm$ 8.07 | 84.54 $\pm$ 7.31 | 2.78 $\pm$ 9.34 | -58.62 $\pm$ 7.22 |
| SLDS | Normal |  | $\nabla$-SLDS | (✓) | 99.93 $\pm$ 0.01 | 95.15 $\pm$ 0.68 | 95.53 $\pm$ 0.47 | -124.65 $\pm$ 6.97 |
| SLDS | Normal |  | GT | (✓) | 99.97 $\pm$ 0.00 | 94.94 $\pm$ 0.68 | 99.36 $\pm$ 0.18 | -129.32 $\pm$ 6.95 |
| Lorenz | Normal |  | identity | ✗ | 27.02 $\pm$ 8.72 | 27.17 $\pm$ 8.74 | 22.87 $\pm$ 7.13 | 24.85 $\pm$ 7.14 |
| Lorenz | Normal |  | LDS | ✗ | 80.30 $\pm$ 14.13 | 82.98 $\pm$ 12.64 | -13.07 $\pm$ 41.03 | 42.08 $\pm$ 26.14 |
| Lorenz | Normal |  | $\nabla$-SLDS | (✓) | 93.91 $\pm$ 5.32 | 93.70 $\pm$ 5.11 | 34.48 $\pm$ 6.47 | 55.75 $\pm$ 6.01 |

## Peer Reviews

confidence: 2
contribution: 3
initial_rating: 6
presentation: 2
questions: None
soundness: 3
strengths: Applying contrastive learning to recover latent dynamics is itself a relatively new approach and the paper is well organised. The proofs use standard jacobian analysis tools and is easy to follow.
summary: The main contribution of this paper is to provide identifiability theory for contrastive learning on time-series data with non-linear mixing, in the same spirit as time-contrastive learning (Hyvärinen & Morioka, 2016) for non-linear ICA. However, the authors discard the independence assumptions typically made in non-linear ICA with respect to the latent variables, and instead define a dynamical system as the data generating process. The proof operates under the assumption that the mapping from latent states to observables is injective but not necessarily linear, which is exploited to show that the composition of mixing and de-mixing by the model is an affine transform. As such, the estimated dynamics via contrastive estimation identify the true dynamics up to affine transformation in the latent space. There are experiments that corroborate the validity of this approach.
weaknesses: The paper needs refinement, with minor typos and inadequately captioned figures. While studying the identifiability of time-series contrastive learning might be novel, all the technical tools require carefully controlled assumptions and specific behavior of Jacobians under contrastive loss minimization, which typically do not hold in practice. Nonetheless, such assumptions are common in the literature on the identifiability of dynamical systems from observed time series.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: See the weaknesses section.
soundness: 3
strengths: 1. The connection between contrastive learning and dynamic system identification is novel and could lead to simple encoder-only implementations favored in practice.
2. The synthetic experiments, especially the ablation studies, are extensive and investigate many aspects of the theoretical results.
summary: The paper explores the use of self-supervised contrastive learning for dynamic system identification. It connects self-supervised learning (SSL) with identifiable representation learning, showing that SSL can identify system dynamics in latent space. The authors propose a model to uncover linear, switching linear, and non-linear dynamics under a non-linear observation model, providing theoretical guarantees and empirical validation.
weaknesses: 1. The contributions (comparison with previous work), theoretical techniques, and novelty are very lightly discussed. I would appreciate a detailed discussion comparing this work's conditions with those in recent literature on temporal causal representation learning (e.g., [1] and its follow-up work). This would aid the contextualization of this work and make the contribution more transparent. 

2. This work strikes me as theoretically oriented. There is a lack of discussion on the theoretical conditions, limitations, and implications. Just to name a few: 

(a) the theorem requires the length for each time series to be infinite, which appears to be quite a strong assumption and not necessitated in recent work (e.g., [1], admittedly I am not an expert and would appreciate correction if I have misunderstood the statements). Why is this necessary, and would it be possible to show results for finite-length series?

(b) How does the control distribution influence the identification results? what conditions should it meet (does it have to be a time-independent Gaussian distribution)?

(c) Line 175: what is the significance of $A$ -- is it previously defined? Could you comment on lines 175-177?

3. Some minor issues:

(a) Line 145: what is (...)?

(b) Line 165: shouldn't $L^{\top} L$ be a matrix?


[1] Temporally Disentangled Representation Learning. Yao et al. NeurIPS 2022.

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 3
questions: My main questions can be summarized as follows:

1. What is the marginal utility of this method rather than various other latent nonlinear dynamics estimation methods, e.g. (Watter et al., 2015) (which in fact also imposes locally linear latent dynamics), which do not make strong assumptions on identifiability?

2. How does enforcing the identifiability/invertibility conditions in this paper affect the method's performance in partially observed settings? This could be as simple as the LQG setting detailed above. Does this strong inductive bias translate to large errors when it is not satisfied (which is typically the case when only provided with observations of a ground-truth belief/latent state)?

3. As detailed above, what is the role of the control input? What is the effective difference of the proposed setting and an autonomous (latent) dynamical system $x_{t+1} = f(x_t) + \epsilon_t$?

4. From a practical perspective, how are the reference points for computing first-order linear approximation in the switching case chosen? Also, do these need to be recomputed per iteration, since the parameterization of $\hat{\mathbf{f}}$ changes per iteration? 

**References:**

Watter et al., "Embed to Control: A Locally Linear Latent Dynamics Model for Control from Raw Images", 2015.
soundness: 3
strengths: Automatic identification of latent variables or dynamics is of critical importance in modern machine / reinforcement learning. The method the authors propose follows a line of self-supervised methods in contrastive learning. In comparison to its closest relative in time-contrastive learning (Hyvarinen and Morioka, 2016), the proposed method is seemingly more well-fit for fitting non-linear time-series data by fitting a latent time-series, rather than predicting a categorical label as in the aforementioned paper.

Since a main inductive bias built into the base method is that the latent dynamics are linear, the proposed method of iterative Jacobian linearizations is a sensible adaptation, and seems to benefit performance significantly.

Numerically, the proposed method appears to make contrastive methods more robustly performant.
summary: In this paper, the authors propose a system-identification scheme for non-linear observations of non-linear time series data. In particular, they propose a modified contrastive learning set-up that posits linear latent dynamics. Compared to prior works in (time) contrastive learning, this directly enforces a notion of sequential temporal consistency, and seems to provide some benefit in system identification settings. Some supporting theory is provided, demonstrating that if the underlying dynamics are linear and invertible, then the proposed method asymptotically recovers the true dynamics up to affine ambiguity. For general non-linear systems, a (soft) switched-linear system heuristic is proposed, where Jacobian linearizations are applied at user-provided reference points.
weaknesses: In my opinion, the paper leaves quite a few critical questions unanswered, and in general suffers from a lack of polish. In its current state, I cannot recommend the paper for acceptance. The main weaknesses in my eyes are the following:

The paper claims to perform latent nonlinear system identification. This is a key desideratum in various fields such as reinforcement learning and continuous control, and thus has a rich history and literature. However, the assumptions in this paper--and notably how these inductive biases propagate to the algorithm design--severely restrict the applicability of the method without further evidence. Notably, a design assumption in this paper is that the observer function (i.e."mixing function") is invertible. This is a very strong assumption in the context of *non-linear system identification*, where even the foundational theory of linear system identification does not presume: in the Linear-Quadratic Gaussian (LQG) model, where the underlying state evolves linearly $x_{t+1} = Ax_t + Bu_t + w_t$, and observations are a linear function of state $y_t = Cx_t + v_t$ (ignoring the control input term for simplicity), the classical set-up has $d_y < d_x$, such that the observations are per-timestep a low-dimensional measurement of the underlying state. This immediately rules out the mixing function $g(x) = Cx$ being invertible, and this is precisely the motivation for notions such as observability/detectability. Partial observability presents the key challenge in non-linear sysID or reinforcement learning. In particular, it is well-known in controls and RL that ignoring partial observability and imposing a Markovian model (which this paper does implicitly by enforcing the state estimate as a function solely of the current observation) can lead to very undesirable outcomes. In the contrastive learning literature, partial observability is usually not a central issue, often because it is irrelevant for the motivating application (e.g. in computer vision), but one must address this problem for time-series data. In fact, the cited Time-Contrastive Learning method (Hyvarinen and Tomioka, 2016), despite making the same assumption in theory, actually propose a method that is more amenable to partial observability, since they predict categorical labels to *chunks* of observed data.

Regarding the polish of the paper, there are various typos and lacking definitions that make the paper hard to parse at times. The minor ones that I have caught are listed below. A particularly confusing point is the role of the control input $u_t$. The paper presents the control input as entering the latent dynamics directly. However, it is typically the case that the control input enters the state through a (possibly state-dependent) actuation matrix $\mathbf B(x_t) u_t$. In any case, how the control input enters the dynamics should be dependent on the parameterization of the dynamics, e.g. the affine ambiguity in $\mathbf L$ in the paper, which is not reflected in the authors' method as far as I can tell. Furthermore, it is unclear if the control input is available to the learner (which is usually the case in sysID), or if it is playing the role of stochastic noise, which eq (9) seems to suggest compared to eq (1). In either case, what role is the control input playing here: in the authors' set-up, there is no need to learn the actuation matrix, and the experiments involve learning a low-noise, nearly deterministic Lorenz system, which rules out some persistency of excitation effect (Tsiamis and Pappas, 2019).

**Minor comments/typos:**

Figure 1: x -> $x$

Page 3: "linear identifiability (...)", missing eqref?

Theorem 1: "bijective dynamics model $\mathbf f$", should probably mathematically define what that means.

Theorem 1: $\lambda$ is not defined in the main paper, only in the appendix.

Corollary 1: "$\hat{\mathbf f} :=1$", seems to be bad notation.

Beginning of Sec 4: "non-lineary" -> "non-linearly"

Equation (7): where is $g_k$ defined? Possible hash collision with mixing function notation.

Table 1: should probably introduce acronym "LDS" = Linear Dynamical System somewhere

Table 1: What does LDS$\downarrow$ mean?

Table 1: What do $\mathbf L$, $\mathbf L'$, $(\checkmark)$, $\mathbf I$ in the theory column indicate?

Between (11) and (12): "tailor" -> "Taylor"

Before Sec 5: "matices" -> "matrices"

Implementation paragraph: possibly missing number of A100 cards?

Eq (23): where is $p_u$ defined?

Eq (25): what does $p_{D}(y)$ denote precisely?

After Eq (50): "which is probably still fine because $\exp(-\|\mathbf L \cdot \|^2)$ is a valid kernel function (?)". This probably needs to be formalized/reworded.

**References:**

Anastasios Tsiamis and George J. Pappas, "Finite Sample Analysis of Stochastic System Identification", 2019.

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 2
questions: - It is confusing why it requires 120 GPU days on A100 for these synthetic data. Methods like TCL require only dozens of minutes for one synthetic data. What is the detailed model structure with approximate parameter size?
soundness: 3
strengths: - The theorem is both interesting and novel, demonstrating that the learned latent variables can be identifiable up to linear transformations even in the absence of independent noise, provided that some other assumptions are met.
- This theorem provides valuable insights into the mechanisms underlying contrastive self-supervised learning methods.
- A lot of ablation studies and visualization experiments are conducted, which makes the paper more convincing.
summary: The paper "Self-Supervised Contrastive Learning Performs Non-Linear System Identification" explores contrastive learning (CL) for identifying non-linear temporal representations. It presents proof of the identifiability of latent variables up to a linear transformation, removing the requirement for independent noise. The proposed model, DynCL, is validated using synthetic data to support the theoretical findings. Additionally, the authors introduce a model called delta-SLDS, designed to capture switching between linear and non-linear dynamic systems.
weaknesses: - Notations are not clear. For example, in Eq (3), the meanings of $x, x', x''$ should be mentioned in advance.
- In theorem (1) and its proof, the assumption is not aligned with Equation (1-2), where all noise disappears. Further discussion is required.
- The difference in theorem part should be compared with previous works on CL more detailed, since it is a work focusing on theorem. Making the difference more clear will make it more readable.
- In experiment parts
    - Baselines like TCL should be compared
    - By identifiability, some metrics like MCC should be compared even though they are not component-wise identifiable.
- Lots of typos:
    - line 133 (supp figure), line 145 (...), seems unfinished part
    - line 250: tailor -> talyor
    - line 637: Theorem 2 -> Theorem A1
    - footnote of page 13: broken reference  
    - line 659, 665: missing reference

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. You use the Gumbel-softmax to approximate the argmin. In my experience, the proper turning of the schedule of its temperature parameter is quite challenging. Is it something you noticed as well?
soundness: 3
strengths: 1. Overall this is an interesting paper, that gives an insight on why CL techniques are effective for system identification
2. The introduced CL variant is well presented and theoretically grounded, and could inspire further theoretical research and models
3. DynCL can effectively identify latent states and system dynamics in the experiments on simulated data
4. The authors present a good selection of ablation studies to demonstrate the impact of the different modeling/parametric choices
summary: This paper focuses on contrastive learning (CL) methods for dynamical systems.
The authors show that under certain assumptions CL performs system identification and can therefore uncover the latent dynamics of the data.
The theoretical findings are applied to switching linear dynamics and non-linear dynamics, and are demonstrated from an empirical point of view using simulated data.
weaknesses: MORE REALISTIC EXPERIMENTS

As stated by the authors in the limitations of this work, the focus of this paper is only on simulated data. While I understand their point of view, and I also agree that the theoretical contribution/simulated experiments are also valuable by themselves, I fear that the impact of this work in this current state will be more limited than it could be if there was a better demonstration of real world applicability.

You could for example try to apply your method to some of the datasets used in "Discovering State Variables Hidden in Experimental Data" (https://www.cs.columbia.edu/~bchen/neural-state-variables/).

If the above is too challenging to achieve, you should at least try to discuss more in detail what each of you theoretical assumptions means in practice, and what you expect to happen if they are not met in real-world experiments. For example, the fact that $p(u_t)$ is a normal distribution seems quite strict in many applications. 

BASELINE

The baseline you use in your experiment seems quite weak, as it does not even use a dynamics model. Have you tried other approaches, for example models doing next-token prediction tasks?




CLARITY

There are several missing definitions/clarifications in the paper that make it a bit harder to follow:
1. N in (3) is not defined
2. "Supp figure" in line 133 is unclear
3. Not sure what "(…)" in line 145 means
4. The name $\nabla$-SLDS is never formally defined
5. In table 1 you have a column called "theory" with different options. What do these option represent exactly?
6. The abbreviation "GT", which I assume stands for ground truth is used in many places but never defined
7. What is $\pi$ in equation (8)?
8. The vMF abbreviation is never defined
9. The DynCL results from Table 1 should be discussed more in depth.

Meta-review: The paper studies representation learning for system identification of the time-series data. Namely, under the assumption that the ground-true time-evolution happens in the latent space, the authors examine the question of whether one can identify the system, i.e. learn the latent representation and the dynamics in the latent space from the time-series observations.

Building on the ideas of contrastive learning, the authors propose a model and prove a theoretical result for this model that in the limit of infinite time-series, the dynamical system can be identified up to a linear transformation (notably, assuming that the observation is obtained from the latent by an injective map). They apply their theory to identify the switching linear dynamics (linear dynamics in which parameters change in time by switching between a finite number of values). Furthermore, the authors argue that the non-linear dynamics can be approximated within the family of switching linear dynamics by the corresponding linearization of the dynamics around the time-series points in the latent space. Finally, the authors empirically validate the proposed framework for linear, switching linear, and Lorenz dynamical systems generating synthetic data in 3 and 6 dimensions.

Despite the reviewers' unanimous low confidence and the much room for improvement left in the presentation, I'm leaning toward accepting the paper. The main reason for this is that different reviewers found different parts of the work to be appealing. This signals that the paper might be relevant beyond a single community of researchers.

Decision: Accept (Poster)