What is a Linear Dynamical System?

A Linear Dynamical System (LDS) is a mathematical model used to describe how a system evolves over time. These systems are a subset of state-space models, where the hidden state dynamics are continuous. What makes these models linear is that the latent dynamics evolve according to a linear function of the previous state. The observations, however, can be related to the hidden state through a nonlinear link function.

At its core, an LDS defines:

• A state transition function: how the internal state evolves from one time step to the next.
• An observation function: how the internal state generates the observed data.

LinearDynamicalSystem{T<:Real, S<:AbstractStateModel{T}, O<:AbstractObservationModel{T}}

The Gaussian Linear Dynamical System

The Gaussian Linear Dynamical System — typically just referred to as an LDS — is a specific type of linear dynamical system where both the state transition and observation functions are linear, and all noise is Gaussian.

The generative model is given by:

\[\begin{aligned} x_t &\sim \mathcal{N}(A x_{t-1} + B u_t + b, Q) \\ y_t &\sim \mathcal{N}(C x_t + D v_t + d, R) \end{aligned}\]

Where:

• $x_t$ is the hidden state at time $t$
• $y_t$ is the observed data at time $t$
• $u_t$ is the dynamics input at time $t$
• $v_t$ is the observation input at time $t$
• $A$ is the state transition matrix
• $C$ is the observation matrix
• $Q$ is the process noise covariance
• $R$ is the observation noise covariance
• $B$ and $D$ are input matrices
• $b$ and $d$ are bias terms

This can equivalently be written in equation form:

\[\begin{aligned} x_t &= A x_{t-1} + b + \epsilon_t \\ y_t &= C x_t + d + \eta_t \end{aligned}\]

Where:

• $\epsilon_t \sim \mathcal{N}(0, Q)$ is the process noise
• $\eta_t \sim \mathcal{N}(0, R)$ is the observation noise

GaussianStateModel{T<:Real, M<:AbstractMatrix{T}, V<:AbstractVector{T}}

GaussianObservationModel{T<:Real, M<:AbstractMatrix{T}, V<:AbstractVector{T}}

The Poisson Linear Dynamical System

The Poisson Linear Dynamical System is a variant of the LDS where the observations are modeled as counts. This is useful in fields like neuroscience where we are often interested in modeling spike count data. To relate the spiking data to the Gaussian latent variable, we use a nonlinear link function, specifically the exponential function.

The generative model is given by:

\[\begin{aligned} x_t &\sim \mathcal{N}(A x_{t-1} + b, Q) \\ y_t &\sim \text{Poisson}(\exp(C x_t + d)) \end{aligned}\]

Where $d$ is a bias term.

PoissonObservationModel{
    T<:Real,
    M<:AbstractMatrix{T},
    V<:AbstractVector{T}
} <: AbstractObservationModel{T}

Sampling from Linear Dynamical Systems

You can generate synthetic data from fitted LDS models. Pass a scalar tsteps::Integer for a single trial, or a vector of per-trial lengths to sample a multi-trial dataset (trial lengths may differ):

Random.rand([rng,] lds, tsteps::Integer; latent_inputs=nothing, obs_inputs=nothing)
Random.rand([rng,] lds, tsteps_per_trial::AbstractVector{<:Integer};
            latent_inputs=nothing, obs_inputs=nothing)

Inference in Linear Dynamical Systems

In StateSpaceDynamics.jl, we directly maximize the complete-data log-likelihood function with respect to the latent states given the data and the parameters of the model. In other words, the maximum a posteriori (MAP) estimate of the latent state path is:

\[\underset{x}{\text{argmax}} \left[ \log p(x_1) + \sum_{t=2}^T \log p(x_t \mid x_{t-1}) + \sum_{t=1}^T \log p(y_t \mid x_t) \right]\]

This MAP estimation approach has the same computational complexity as traditional Kalman filtering and smoothing — $\mathcal{O}(T)$ — but is significantly more flexible. Notably, it can handle nonlinear observations and non-Gaussian noise while still yielding exact MAP estimates, unlike approximate techniques such as the Extended Kalman Filter (EKF) or Unscented Kalman Filter (UKF).

Newton's Method for Latent State Optimization

To find the MAP trajectory, we iteratively optimize the latent states using Newton's method. The update equation at each iteration is:

\[x^{(i+1)} = x^{(i)} - \left[ \nabla^2 \mathcal{L}(x^{(i)}) \right]^{-1} \nabla \mathcal{L}(x^{(i)})\]

Where:

• $\mathcal{L}(x)$ is the complete-data log-likelihood:

\[\mathcal{L}(x) = \log p(x_1) + \sum_{t=2}^T \log p(x_t \mid x_{t-1}) + \sum_{t=1}^T \log p(y_t \mid x_t)\]

• $\nabla \mathcal{L}(x)$ is the gradient of the full log-likelihood with respect to all latent states
• $\nabla^2 \mathcal{L}(x)$ is the Hessian of the full log-likelihood

This update is performed over the entire latent state sequence $x_{1:T}$, and repeated until convergence.

For Gaussian models, $\mathcal{L}(x)$ is quadratic and Newton's method converges in a single step — recovering the exact Kalman smoother solution. For non-Gaussian models, the Hessian is not constant and the optimization is more complex. However, the MAP estimate can still be computed efficiently using the same approach as the optimization problem is still convex.

smooth(lds, y::AbstractMatrix)

Laplace Approximation of Posterior for Non-Conjugate Observation Models

In the case of non-Gaussian observations, we can use a Laplace approximation to compute the posterior distribution of the latent states. For Gaussian observations (which are conjugate with the Gaussian state model), the posterior is also Gaussian and is the exact posterior. However, for non-Gaussian observations, we can approximate the posterior using a Gaussian distribution centered at the MAP estimate of the latent states. This approximation is given by:

\[p(x \mid y) \approx \mathcal{N}(x^{*}, -\left[ \nabla^2 \mathcal{L}(x^{*}) \right]^{-1})\]

Where:

• $x^{*}$ is the MAP estimate of the latent states
• $\nabla^2 \mathcal{L}(x^{*})$ is the Hessian of the log-likelihood at the MAP estimate

Despite the requirement of inverting a Hessian of dimension $(d \times T) \times (d \times T)$, this is still computationally efficient, as the Markov structure of the model renders the Hessian block-tridiagonal, and thus the inversion is tractable.

Learning in Linear Dynamical Systems

Given the latent structure of state-space models, we must rely on either the Expectation-Maximization (EM) or Variational Inference (VI) approaches to learn the parameters of the model. StateSpaceDynamics.jl supports both EM and VI. For LDS models, we can use Laplace EM, where we approximate the posterior of the latent state path using the Laplace approximation as outlined above. Using these approximate posteriors (or exact ones in the Gaussian case), we can apply closed-form updates for the model parameters.

fit!(lds, y; max_iter=100, tol=1e-6, progress=true,
     latent_inputs=nothing, obs_inputs=nothing)

Inverse-Wishart Priors on Covariances (MAP)

To encourage well-conditioned covariance estimates, StateSpaceDynamics.jl supports Inverse-Wishart (IW) priors on the covariance matrices of a Linear Dynamical System. These can be used to impose shrinkage toward a target scale, yielding maximum a posteriori (MAP) estimates instead of pure MLEs.

You can attach an IW prior to:

• The process noise covariance Q_prior
• The initial state covariance P0_prior
• The observation noise covariance R_prior (Gaussian models only)

IWPrior{T<:Real, M<:AbstractMatrix}

Definition

For a covariance matrix $\Sigma \in \mathbb{R}^{d\times d}$,

\[\Sigma \sim \text{IW}(\Psi, \nu), \qquad p(\Sigma)\propto |\Sigma|^{-(\nu+d+1)/2}\exp\!\Big[-\tfrac12\operatorname{tr}(\Psi\,\Sigma^{-1})\Big].\]

• $\Psi$ is the scale matrix (positive-definite).
• $\nu$ is the degrees of freedom. A proper mode exists when $\nu > d+1$.

Given $n$ effective samples and sample statistic $S$, the posterior is

\[\Sigma\mid\text{data}\sim\text{IW}(\Psi+S,\;\nu+n)\]

and the MAP (mode) update used internally is

\[\widehat{\Sigma}_{\text{MAP}} =\frac{\Psi+S}{\nu+n+d+1}.\]

Example: Adding IW Priors to an LDS

using LinearAlgebra, Random, StateSpaceDynamics

rng = MersenneTwister(42)
D, P = 3, 4

# Proper SPD scale matrices (avoid UniformScaling)
ΨQ  = diagm(0 => fill(0.01, D))
ΨP0 = diagm(0 => fill(0.01, D))
ΨR  = diagm(0 => fill(0.01, P))

Qprior  = IWPrior(Ψ = ΨQ,  ν = D + 3.0)
P0prior = IWPrior(Ψ = ΨP0, ν = D + 3.0)
Rprior  = IWPrior(Ψ = ΨR,  ν = P + 3.0)

# Gaussian LDS
A = 0.9I + 0.05randn(rng, D, D)
Q = Matrix(I, D, D) .* 0.3
b = zeros(D); x0 = zeros(D); P0 = Matrix(I, D, D) .* 0.8

C = randn(rng, P, D)
R = Matrix(I, P, P) .* 0.25
d = 0.1 .* randn(rng, P)

gsm = GaussianStateModel(A=A, Q=Q, b=b, x0=x0, P0=P0,
                         Q_prior=Qprior, P0_prior=P0prior)
gom = GaussianObservationModel(C=C, R=R, d=d, R_prior=Rprior)
lds = LinearDynamicalSystem(gsm, gom)

X, Y = rand(rng, lds, fill(100, 5))   # 5 trials of 100 timesteps each
fit!(lds, Y; max_iter=20, progress=false)

Effect on Learning

• When a prior is present, the M-step uses the MAP mode formula above.
• The reported ELBO includes IW log-prior terms (up to constants), so convergence still tracks the true MAP objective.
• With small datasets or high-dimensional states, the priors keep $Q$, $P_0$, and $R$ positive-definite and well-conditioned.

Choosing $\Psi$ and $\nu$

Tip: In Julia, I is a UniformScaling (not a matrix). Use Matrix(I, d, d) or diagm(0 => fill(..., d)) to build $Ψ$.

Matrix-Normal Priors on Regression Coefficients (MAP)

The M-step updates for the dynamics and emission parameters are (multivariate) linear regressions: the dynamics update fits the stacked coefficient matrix $[A\; b]$ and the emission update fits $[C\; d]$. StateSpaceDynamics.jl supports a matrix-normal (MN) prior on these coefficient matrices, turning each regression into a MAP (ridge-style) update:

\[W_{\text{MAP}} = (XY^\top + M_0 \Lambda)(XX + \Lambda)^{-1}\]

where $M_0$ is the prior mean and $\Lambda$ the column precision. Setting $M_0 = 0$ recovers ordinary ridge regression; $\Lambda = 0$ recovers the unpenalized least-squares update.

You can attach an MN prior to:

• The stacked dynamics matrix via AB_prior on a GaussianStateModel
• The stacked emission matrix via CD_prior on a GaussianObservationModel or PoissonObservationModel

The MN prior is kept separate from the covariance priors on purpose: pairing an AB_prior/CD_prior with the matching Q_prior/R_prior (IWPrior) recovers a full matrix-normal-inverse-Wishart (MNIW) prior, but each half can also be used on its own — e.g. the Poisson observation model has no noise covariance, so only the MN half applies there.

MNPrior{T<:Real, M<:AbstractMatrix}