Causal Optimal Coupling for Gaussian Input-Output Distributional Data

1. Problem Statement

The paper addresses the problem of system identification from distributional data where the underlying pairing between input and output observations is unknown.

• Context: Given two independent stochastic processes, an input $U = \{U_t\}_{t=1}^T \sim \mu$ and an output $Y = \{Y_t\}_{t=1}^T \sim \nu$, with known marginal distributions but unknown joint law.
• Objective: Identify an optimal coupling $\pi$ (a joint probability distribution) between $\mu$ and $\nu$ that represents the input-output relationship.
• Constraints:
  1. Marginal Constraints: The coupling must have marginals $\mu$ and $\nu$.
  2. Causality Constraint: The coupling must respect the temporal structure of the system. Specifically, future inputs must be independent of past outputs given past inputs:
     $$(Y_1, \dots, Y_t) \perp\!\!\!\perp (U_{t+1}, \dots, U_T) \mid (U_1, \dots, U_t)$$
     This ensures the model is non-anticipative.
• Optimization Goal: Find the coupling $\pi$ closest to a reference prior model $\gamma$ (encoding structural assumptions) in terms of Kullback-Leibler (KL) divergence, subject to the above constraints. This is formulated as an Entropic Causal Optimal Transport (ECOT) problem:
  $$\min_{\pi \in \Pi_c(\mu, \nu)} D_{KL}(\pi, \gamma)$$
  where $\Pi_c(\mu, \nu)$ is the set of causal couplings.

2. Methodology

The authors propose a solution based on Schrödinger Bridge theory and Sinkhorn iterations, specialized for the Gaussian setting.

A. Reference Model and Cost Function

The reference coupling $\gamma$ is induced by a linear Gaussian causal input-output model:
$$Y_t = h_t^\top U_{1:t} + f_t^\top Y_{1:t-1} + b_t + w_t$$
where $w_t \sim \mathcal{N}(0, \epsilon_t^2)$.

• This choice transforms the KL divergence minimization into a problem with a general time-dependent quadratic cost function:
  $$c(u_{1:T}, y_{1:T}) = \sum_{t=1}^T \frac{1}{2\epsilon_t^2} (y_t - h_t^\top u_{1:t} - f_t^\top y_{1:t-1} - b_t)^2$$
• This cost captures temporal dependencies and interactions, generalizing the standard squared loss ($\sum (y_t - u_t)^2$).

B. Sinkhorn Iterations in the Gaussian Setting

The authors adapt the abstract Sinkhorn algorithm (alternating projections) to the specific case where marginals and the prior are Gaussian.

• Invariance Property: The set of Gaussian distributions is invariant under the projection steps required by the Sinkhorn algorithm. This allows the iterations to be performed entirely via updates to Gaussian parameters (means and covariances) rather than density functions.
• Algorithm Structure:
  1. Initialization: Start with the reference coupling $\pi^{(0)} = \gamma$.
  2. Odd Iteration (Causal Projection): Project onto the set of couplings with fixed input marginal $\mu$ and satisfying the causality constraint ($\Pi_c(\mu, *)$).
     • The solution remains Gaussian.
     • The authors derive Theorem 1, providing explicit formulas for the conditional parameters of the updated coupling. Crucially, the conditional distribution of $U_t$ given pasts is fixed to match $\mu$, while the conditional distribution of $Y_t$ is updated as a convex combination of the prior and a correction term derived from the future constraints.
  3. Even Iteration (Marginal Projection): Project onto the set of couplings with fixed output marginal $\nu$ ($\Pi(*, \nu)$), ignoring the causality constraint for this step.
     • The solution remains Gaussian.
     • Theorem 2 shows that in this step, the conditional distribution of $U$ given $Y$ reverts to the prior $\gamma$, while the marginal of $Y$ is forced to match $\nu$.
• Convergence: The algorithm alternates between these two steps until convergence, yielding a joint Gaussian coupling $\pi^*$.

C. Model Recovery

Once the optimal coupling $\pi^*$ is found, the authors show how to recover the parameters of the optimal linear Gaussian input-output model ($h_t, f_t, b_t, \epsilon_t$) by conditioning $Y_t$ on the history $U_{1:t}, Y_{1:t-1}$ using standard Gaussian conditional formulas (covariance and mean extraction).

3. Key Contributions

1. Formulation of ECOT for System Identification: The paper formalizes the identification of dynamical systems from distributional data as an Entropic Causal Optimal Transport problem, bridging optimal transport theory with control/system identification.
2. Tractable Gaussian Characterization: Unlike previous works that rely on abstract iterations or restrict costs to separable forms, this paper provides a fully explicit, closed-form characterization of Sinkhorn iterations for:
   • General quadratic costs (including time-lagged inputs and outputs).
   • One-sided causality constraints.
   • Gaussian marginals.
3. Algorithmic Implementation: The derivation of specific update rules for the mean and covariance matrices makes the algorithm computationally feasible and implementable, avoiding the need for discretization of the state space.
4. Distinction from Prior Work:
   • Extends Veraguas et al. (2017) which was limited to Markovian inputs and separable costs.
   • Makes the abstract Sinkhorn iterations of Eckstein and Pammer (2023) concrete and tractable for the Gaussian case.
   • Unlike Gunasingam and Wong (2025) or Acciaio et al. (2025), it imposes only one-sided causality (input $\to$ output), which is more appropriate for standard system identification, rather than bicausal constraints.

4. Results

• Numerical Example: The authors tested the algorithm on a discretized time interval ($T=128$) with stationary Gaussian marginals.
  • Setup: Input and output processes were defined with specific covariance kernels. The reference model was a linear Gaussian state-space model.
  • Comparison: The causal solution was compared against a non-causal (standard entropic OT) solution.
• Observations:
  • Causal Coupling: The conditional cross-covariance $Cov(U_t, Y_s | U_{1:\tau})$ was zero for $s < \tau$ (past outputs do not depend on future inputs), strictly adhering to the non-anticipative constraint.
  • Non-Causal Coupling: The non-causal solution exhibited non-zero dependencies on past outputs for future inputs (anticipatory correlations) and displayed sign changes/anti-symmetric patterns, violating physical causality.
  • Convergence: The alternating projection scheme converged rapidly (tolerance $10^{-6}$).

5. Significance and Future Directions

• Significance: This work provides a principled mathematical framework for learning dynamical systems directly from distributional data (e.g., ensembles of trajectories) rather than point-wise paired data. It ensures the learned model respects the fundamental physical law of causality.
• Applications: The framework is applicable to system identification in control theory, signal processing, and any domain where only statistical properties of inputs and outputs are available, but the pairing is lost.
• Future Work: The authors plan to extend the framework to non-linear reference models and non-Gaussian marginals, aiming to apply these methods to generative modeling of complex time-series data.

In summary, the paper successfully bridges the gap between theoretical optimal transport and practical system identification by providing a computationally efficient, causality-preserving algorithm for Gaussian systems.