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.

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:

$

xt \sim \mathcal{N}(A x{t-1}, Q) \ yt \sim \mathcal{N}(C xt, R) $

Where:

  • \[x_t\]

    is the hidden state at time $t$

  • \[y_t\]

    is the observed data 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

This can equivalently be written in equation form:

$

xt = A x{t-1} + \epsilont \ yt = C xt + \etat $

Where:

  • \[\epsilon_t \sim \mathcal{N}(0, Q)\]

    is the process noise

  • \[\eta_t \sim \mathcal{N}(0, R)\]

    is the observation noise

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 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. Thus our generative model is given by:

\[x_t \sim \mathcal{N}(A x_{t-1}, Q) \\ y_t \sim \text{Poisson}(\text{exp}(Cx_t + b)) \]

Where:

  • \[b\]

    is a bias term

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 priori (MAP) estimate of the latent state path is:

\[\underset{x}{\text{argmax}} \left[ \log p(x_0) + \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_0) + \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. But, 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.

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. Notably, in the case of Gaussian Observations (which is 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 diomension (d x T) x (d x T), this is still computationally efficient, as the Markov structure of the model, renders the Hessian block-tridiagonal, and thus the inversion is not intractable.

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.

StateSpaceDynamics.GaussianObservationModel — Type
GaussianObservationModel{T<:Real} <: AbstractObservationModel

Represents the observation model of a Linear Dynamical System with Gaussian noise.

Fields

  • C::Matrix{T}: Observation matrix
  • R::Matrix{T}: Observation noise covariance
StateSpaceDynamics.GaussianObservationModel — Method
GaussianObservationModel(; C, R, obs_dim, latent_dim)

Construct a GaussianObservationModel with the given parameters or random initializations.

Arguments

  • C::Matrix{T}=Matrix{T}(undef, 0, 0): Observation matrix
  • R::Matrix{T}=Matrix{T}(undef, 0, 0): Observation noise covariance
  • obs_dim::Int: Dimension of the observations (required if C or R is not provided.)
  • latent_dim::Int: Dimension of the latent state (required if C is not provided.)
StateSpaceDynamics.GaussianStateModel — Type
GaussianStateModel{T<:Real} <: AbstractStateModel

Represents the state model of a Linear Dynamical System with Gaussian noise.

Fields

  • A::Matrix{T}: Transition matrix
  • Q::Matrix{T}: Process noise covariance
  • x0::Vector{T}: Initial state
  • P0::Matrix{T}: Initial state covariance
StateSpaceDynamics.GaussianStateModel — Method
GaussianStateModel(; A, Q, x0, P0, latent_dim)

Construct a GaussianStateModel with the given parameters or random initializations.

Arguments

  • A::Matrix{T}=Matrix{T}(undef, 0, 0): Transition matrix
  • Q::Matrix{T}=Matrix{T}(undef, 0, 0): Process noise covariance
  • x0::Vector{T}=Vector{T}(undef, 0): Initial state
  • P0::Matrix{T}=Matrix{T}(undef, 0, 0): Initial state covariance
  • latent_dim::Int: Dimension of the latent state (required if any matrix is not provided.)
StateSpaceDynamics.LinearDynamicalSystem — Type
LinearDynamicalSystem{S<:AbstractStateModel, O<:AbstractObservationModel}

Represents a unified Linear Dynamical System with customizable state and observation models.

Fields

  • state_model::S: The state model (e.g., GaussianStateModel)
  • obs_model::O: The observation model (e.g., GaussianObservationModel or PoissonObservationModel)
  • latent_dim::Int: Dimension of the latent state
  • obs_dim::Int: Dimension of the observations
  • fit_bool::Vector{Bool}: Vector indicating which parameters to fit during optimization
StateSpaceDynamics.PoissonObservationModel — Type
PoissonObservationModel{T<:Real} <: AbstractObservationModel

Represents the observation model of a Linear Dynamical System with Poisson observations.

Fields

  • C::Matrix{T}: Observation matrix
  • log_d::Vector{T}: Mean firing rate vector (log space)
StateSpaceDynamics.PoissonObservationModel — Method
PoissonObservationModel(; C, log_d, obs_dim, latent_dim)

Construct a PoissonObservationModel with the given parameters or random initializations.

Arguments

  • C::Matrix{T}=Matrix{T}(undef, 0, 0): Observation matrix
  • log_d::Vector{T}=Vector{T}(undef, 0): Mean firing rate vector (log space)
  • obs_dim::Int: Dimension of the observations (required if any matrix is not provided.)
  • latent_dim::Int: Dimension of the latent state (required if C is not provided.)
StateSpaceDynamics.GaussianLDS — Method
GaussianLDS(; A, C, Q, R, x0, P0, fit_bool, obs_dim, latent_dim)

Construct a Linear Dynamical System with Gaussian state and observation models.

Arguments

  • A::Matrix{T}=Matrix{T}(undef, 0, 0): Transition matrix
  • C::Matrix{T}=Matrix{T}(undef, 0, 0): Observation matrix
  • Q::Matrix{T}=Matrix{T}(undef, 0, 0): Process noise covariance
  • R::Matrix{T}=Matrix{T}(undef, 0, 0): Observation noise covariance
  • x0::Vector{T}=Vector{T}(undef, 0): Initial state
  • P0::Matrix{T}=Matrix{T}(undef, 0, 0): Initial state covariance
  • fit_bool::Vector{Bool}=fill(true, 6): Vector indicating which parameters to fit during optimization
  • obs_dim::Int: Dimension of the observations (required if C or R is not provided.)
  • latent_dim::Int: Dimension of the latent state (required if A, Q, x0, P0, or C is not provided.)
StateSpaceDynamics.Gradient — Method
Gradient(lds::LinearDynamicalSystem{S,O}, y::AbstractMatrix{T}, x::AbstractMatrix{T}) where {T<:Real, S<:GaussianStateModel{T}, O<:GaussianObservationModel{T}}

Compute the gradient of the log-likelihood with respect to the latent states for a linear dynamical system.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System.
  • y::AbstractMatrix{T}: The observed data.
  • x::AbstractMatrix{T}: The latent states.
  • w::Vector{Float64}: coeffcients to weight the data.

Returns

  • grad::Matrix{T}: Gradient of the log-likelihood with respect to the latent states.
StateSpaceDynamics.Gradient — Method
Gradient(lds::LinearDynamicalSystem{S,O}, y::Matrix{T}, x::Matrix{T}) where {T<:Real, S<:GaussianStateModel{T}, O<:PoissonObservationModel{T}}

Calculate the gradient of the log-likelihood of a Poisson Linear Dynamical System model for a single trial.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System model.
  • y::Matrix{T}: The observed data. Dimensions: (obsdim, Tsteps)
  • x::Matrix{T}: The latent state variables. Dimensions: (latentdim, Tsteps)
  • w::Vector{T}: Weights for each observation in the log-likelihood calculation. Not currently used.

Returns

  • grad::Matrix{T}: The gradient of the log-likelihood. Dimensions: (latentdim, Tsteps)

Note

The gradient is computed with respect to the latent states x. Each row of the returned gradient corresponds to the gradient for a single time step.

StateSpaceDynamics.Hessian — Method
Hessian(lds::LinearDynamicalSystem{S,O}, y::AbstractMatrix{T}, x::AbstractMatrix{T}) where {T<:Real, S<:GaussianStateModel{T}, O<:GaussianObservationModel{T}}

Construct the Hessian matrix of the log-likelihood of the LDS model given a set of observations.

This function is used for the direct optimization of the log-likelihood as advocated by Paninski et al. (2009). The block tridiagonal structure of the Hessian is exploited to reduce the number of parameters that need to be computed, and to reduce the memory requirements. Together with the gradient, this allows for Kalman Smoothing to be performed by simply solving a linear system of equations:

̂xₙ₊₁ = ̂xₙ - H \ ∇

where ̂xₙ is the current smoothed state estimate, H is the Hessian matrix, and ∇ is the gradient of the log-likelihood.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System.
  • y::AbstractMatrix{T}: Matrix of observations.
  • x::AbstractMatrix{T}: Matrix of latent states.
  • w::Vector{Float64}: coeffcients to weight the data.

Returns

  • H::Matrix{T}: Hessian matrix of the log-likelihood.
  • H_diag::Vector{Matrix{T}}: Main diagonal blocks of the Hessian.
  • H_super::Vector{Matrix{T}}: Super-diagonal blocks of the Hessian.
  • H_sub::Vector{Matrix{T}}: Sub-diagonal blocks of the Hessian.

Note

  • x is not used in this function, but is required to match the function signature of other Hessian calculations e.g., in PoissonLDS.
StateSpaceDynamics.Hessian — Method
Hessian(lds::LinearDynamicalSystem{S,O}, y::AbstractMatrix{T}, x::AbstractMatrix{T}) where {T<:Real, S<:GaussianStateModel{T}, O<:PoissonObservationModel{T}}

Calculate the Hessian matrix of the log-likelihood for a Poisson Linear Dynamical System.

This function computes the Hessian matrix, which represents the second-order partial derivatives of the log-likelihood with respect to the latent states.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System with Poisson observations.
  • y::AbstractMatrix{T}: The observed data. Dimensions: (obsdim, Tsteps)
  • x::AbstractMatrix{T}: The current estimate of latent states. Dimensions: (latentdim, Tsteps)
  • w::Vector{T}: Weights for each observation in the log-likelihood calculation. Not currently used.

Returns

  • H::Matrix{T}: The full Hessian matrix.
  • H_diag::Vector{Matrix{T}}: The main diagonal blocks of the Hessian.
  • H_super::Vector{Matrix{T}}: The super-diagonal blocks of the Hessian.
  • H_sub::Vector{Matrix{T}}: The sub-diagonal blocks of the Hessian.
StateSpaceDynamics.PoissonLDS — Method
PoissonLDS(; A, C, Q, log_d, x0, P0, refractory_period, fit_bool, obs_dim, latent_dim)

Construct a Linear Dynamical System with Gaussian state and Poisson observation models.

Arguments

  • A::Matrix{T}=Matrix{T}(undef, 0, 0): Transition matrix
  • C::Matrix{T}=Matrix{T}(undef, 0, 0): Observation matrix
  • Q::Matrix{T}=Matrix{T}(undef, 0, 0): Process noise covariance
  • log_d::Vector{T}=Vector{T}(undef, 0): Mean firing rate vector (log space)
  • x0::Vector{T}=Vector{T}(undef, 0): Initial state
  • P0::Matrix{T}=Matrix{T}(undef, 0, 0): Initial state covariance
  • refractory_period::Int=1: Refractory period
  • fit_bool::Vector{Bool}=fill(true, 7): Vector indicating which parameters to fit during optimization
  • obs_dim::Int: Dimension of the observations (required if C, D, or log_d is not provided.)
  • latent_dim::Int: Dimension of the latent state (required if A, Q, x0, P0, or C is not provided.)
StateSpaceDynamics.Q_function — Function
Q(A, Q, H, R, P0, x0, E_z, E_zz, E_zz_prev, y)

Calculate the complete Q-function for the EM algorithm in a Linear Dynamical System.

Arguments

  • A::Matrix{<:Real}: The state transition matrix.
  • Q::AbstractMatrix{<:Real}: The process noise covariance matrix (or its Cholesky factor).
  • H::Matrix{<:Real}: The observation matrix.
  • R::AbstractMatrix{<:Real}: The observation noise covariance matrix (or its Cholesky factor).
  • P0::AbstractMatrix{<:Real}: The initial state covariance matrix (or its Cholesky factor).
  • x0::Vector{<:Real}: The initial state mean.
  • E_z::Matrix{<:Real}: The expected latent states, size (state_dim, T).
  • E_zz::Array{<:Real, 3}: The expected value of zt * zt', size (statedim, statedim, T).
  • E_zz_prev::Array{<:Real, 3}: The expected value of zt * z{t-1}', size (statedim, statedim, T).
  • y::Matrix{<:Real}: The observed data, size (obs_dim, T).

Returns

  • Q_val::Float64: The complete Q-function value.
StateSpaceDynamics.Q_function — Method
Q_function(A::Matrix{T}, Q::Matrix{T}, C::Matrix{T}, log_d::Vector{T}, x0::Vector{T}, P0::Matrix{T}, E_z::Matrix{T}, E_zz::Array{T, 3}, E_zz_prev::Array{T, 3}, P_smooth::Array{T, 3}, y::Matrix{T})

Calculate the Q-function for the Linear Dynamical System.

Arguments

  • A::Matrix{T}: The transition matrix.
  • Q::Matrix{T}: The process noise covariance matrix.
  • C::Matrix{T}: The observation matrix.
  • log_d::Vector{T}: The mean firing rate vector in log space.
  • x0::Vector{T}: The initial state mean.
  • P0::Matrix{T}: The initial state covariance matrix.
  • E_z::Matrix{T}: The expected latent states.
  • E_zz::Array{T, 3}: The expected latent states x the latent states.
  • E_zz_prev::Array{T, 3}: The expected latent states x the previous latent states.
  • P_smooth::Array{T, 3}: The smoothed state covariances.
  • y::Matrix{T}: The observed data.

Returns

  • Float64: The Q-function for the Linear Dynamical System.
StateSpaceDynamics.Q_obs — Function
Q_obs(H, R, E_z, E_zz, y)

Calculate the observation component of the Q-function for the EM algorithm in a Linear Dynamical System.

Arguments

  • H::Matrix{<:Real}: The observation matrix.
  • R::AbstractMatrix{<:Real}: The observation noise covariance matrix (or its Cholesky factor).
  • E_z::Matrix{<:Real}: The expected latent states, size (state_dim, T).
  • E_zz::Array{<:Real, 3}: The expected value of zt * zt', size (statedim, statedim, T).
  • y::Matrix{<:Real}: The observed data, size (obs_dim, T).

Returns

  • Q_val::Float64: The observation component of the Q-function.
StateSpaceDynamics.Q_obs — Method
Q_obs(H, R, E_z, E_zz, y)

Calculate the a single time step observation component of the Q-function for the EM algorithm in a Linear Dynamical System before the R^-1 is accounted for.

Arguments

  • H::Matrix{<:Real}: The observation matrix.
  • R::AbstractMatrix{<:Real}: The observation noise covariance matrix (or its Cholesky factor).
  • E_z::Vector{<:Real}: The expected latent states at time t, size (state_dim).
  • E_zz::Matrix{<:Real}: The expected value of zt * zt' at time t, size (statedim, statedim).
  • y::Vector{<:Real}: The observed data at time t, size (obs_dim).

Returns

  • q::Float64: The observation component at time t of the Q-function prior to R^-1.
StateSpaceDynamics.Q_observation_model — Method
Q_observation_model(C::Matrix{<:Real}, D::Matrix{<:Real}, log_d::Vector{<:Real}, E_z::Array{<:Real}, E_zz::Array{<:Real}, y::Array{<:Real})

Calculate the Q-function for the observation model.

Arguments

  • C::Matrix{<:Real}: The observation matrix.
  • log_d::Vector{<:Real}: The mean firing rate vector in log space.
  • E_z::Array{<:Real}: The expected latent states.
  • E_zz::Array{<:Real}: The expected latent states x the latent states.
  • y::Array{<:Real}: The observed data.

Returns

  • Float64: The Q-function for the observation model.
StateSpaceDynamics.Q_state — Method
Q_state(A, Q, P0, x0, E_z, E_zz, E_zz_prev)

Calculate the state component of the Q-function for the EM algorithm in a Linear Dynamical System.

Arguments

  • A::Matrix{<:Real}: The state transition matrix.
  • Q::AbstractMatrix{<:Real}: The process noise covariance matrix (or its Cholesky factor).
  • P0::AbstractMatrix{<:Real}: The initial state covariance matrix (or its Cholesky factor).
  • x0::Vector{<:Real}: The initial state mean.
  • E_z::Matrix{<:Real}: The expected latent states, size (state_dim, T).
  • E_zz::Array{<:Real, 3}: The expected value of zt * zt', size (statedim, statedim, T).
  • E_zz_prev::Array{<:Real, 3}: The expected value of zt * z{t-1}', size (statedim, statedim, T).

Returns

  • Q_val::Float64: The state component of the Q-function.
StateSpaceDynamics.Q_state — Method
Q_state(A::Matrix{T}, Q::Matrix{T}, P0::Matrix{T}, x0::Vector{T}, E_z::Array{T, 3}, E_zz::Array{T, 4}, E_zz_prev::Array{T, 4}) where T<:Real

Calculates the Q-function for the state model over multiple trials.

Arguments

  • A::Matrix{T}: The transition matrix.
  • Q::Matrix{T}: The process noise covariance matrix.
  • P0::Matrix{T}: The initial state covariance matrix.
  • x0::Vector{T}: The initial state mean.
  • E_z::Array{T, 3}: The expected latent states.
  • E_zz::Array{T, 4}: The expected latent states x the latent states.
  • E_zz_prev::Array{T, 4}: The expected latent states x the previous latent states.

Returns

  • Float64: The Q-function for the state model.
StateSpaceDynamics.calculate_elbo — Method
calculate_elbo(lds::LinearDynamicalSystem{S,O}, E_z::Array{T,3}, E_zz::Array{T,4}, E_zz_prev::Array{T,4}, p_smooth::Array{T,4}, y::Array{T,3}) where {T<:Real, S<:GaussianStateModel{T}, O<:AbstractObservationModel{T}}

Calculate the Evidence Lower Bound (ELBO) for a Linear Dynamical System.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System struct
  • E_z::Array{T,3}: Expected latent states, size (statedim, statedim, Tsteps, ntrials)
  • E_zz::Array{T,4}: Expected zt * zt', size (statedim, statedim, Tsteps, ntrials, state_dim)
  • E_zz_prev::Array{T,4}: Expected zt * z{t-1}', size (statedim, statedim, Tsteps, ntrials, state_dim)
  • p_smooth::Array{T,4}: Smoothed state covariances, size (statedim, statedim, Tsteps, ntrials, state_dim)
  • y::Array{T,3}: Observed data, size (obsdim, Tsteps, n_trials)

Returns

  • elbo::T: The Evidence Lower Bound (ELBO) for the LDS.

Note

  • For a GaussianLDS the ELBO is equivalent to the total marginal likelihood
StateSpaceDynamics.calculate_elbo — Method
calculate_elbo(plds::LinearDynamicalSystem{S,O}, E_z::Array{T, 3}, E_zz::Array{T, 4}, 
               E_zz_prev::Array{T, 4}, P_smooth::Array{T, 4}, y::Array{T, 3}) where 
               {T<:Real, S<:GaussianStateModel{T}, O<:PoissonObservationModel{T}}

Calculate the Evidence Lower Bound (ELBO) for a Poisson Linear Dynamical System (PLDS).

Arguments

  • plds::LinearDynamicalSystem{S,O}: The PLDS model.
  • E_z::Array{T, 3}: Expected values of latent states. Dimensions: (statedim, tsteps, n_trials).
  • E_zz::Array{T, 4}: Expected values of latent state outer products. Dimensions: (statedim, statedim, tsteps, ntrials).
  • E_zz_prev::Array{T, 4}: Expected values of latent state outer products with previous time step. Dimensions: (state dimension, state dimension, tsteps-1, ntrials).
  • P_smooth::Array{T, 4}: Smoothed covariance matrices. Dimensions: (state dimension, state dimension, tsteps, ntrials).
  • y::Array{T, 3}: Observed data. Dimensions: (obsdim, tsteps, n_trials).

Returns

  • elbo::Float64: The calculated Evidence Lower Bound.

Description

This function computes the ELBO for a PLDS model, which consists of two main components:

  1. The expected complete log-likelihood (ECLL), calculated using the Q_function.
  2. The entropy of the variational distribution, calculated using gaussian entropy.

The ELBO is then computed as: ELBO = ECLL - Entropy.

Note

Ensure that the dimensions of input arrays match the expected dimensions as described in the arguments section.

StateSpaceDynamics.estep — Method
estep(lds::LinearDynamicalSystem{S,O}, y::Array{T,3}) where {T<:Real, S<:GaussianStateModel{T}, O<:AbstractObservationModel{T}}

Perform the E-step of the EM algorithm for a Linear Dynamical System, treating all input as multi-trial.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System struct
  • y::Array{T,3}: Observed data, size (obsdim, Tsteps, ntrials) Note: For single-trial data, use y[1:1, :, :] to create a 3D array with ntrials = 1

Returns

  • E_z::Array{T,3}: Expected latent states, size (statedim, statedim, Tsteps, ntrials)
  • E_zz::Array{T,4}: Expected zt * zt', size (statedim, statedim, Tsteps, ntrials, state_dim)
  • E_zz_prev::Array{T,4}: Expected zt * z{t-1}', size (statedim, statedim, Tsteps, ntrials, state_dim)
  • x_smooth::Array{T,3}: Smoothed state estimates, size (statedim, statedim, Tsteps, ntrials)
  • p_smooth::Array{T,4}: Smoothed state covariances, size (statedim, statedim, Tsteps, ntrials, state_dim)
  • ml::T: Total marginal likelihood (log-likelihood) of the data across all trials

Note

  • This function first smooths the data using the smooth function, then computes sufficient statistics.
  • It treats all input as multi-trial, with single-trial being a special case where n_trials = 1.
StateSpaceDynamics.fit! — Method
fit!(lds::LinearDynamicalSystem{S,O}, y::Matrix{T}; 
     max_iter::Int=1000, 
     tol::Real=1e-12, 
     ) where {T<:Real, S<:GaussianStateModel{T}, O<:GaussianObservationModel{T}}

Fit a Linear Dynamical System using the Expectation-Maximization (EM) algorithm with Kalman smoothing.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System to be fitted.
  • y::Matrix{T}: Observed data, size (obsdim, Tsteps).

Keyword Arguments

  • max_iter::Int=1000: Maximum number of EM iterations.
  • tol::Real=1e-12: Convergence tolerance for log-likelihood change.

Returns

  • mls::Vector{T}: Vector of log-likelihood values for each iteration.
StateSpaceDynamics.gradient_observation_model! — Method
gradient_observation_model!(grad::AbstractVector{T}, C::AbstractMatrix{T}, log_d::AbstractVector{T}, E_z::AbstractArray{T}, P_smooth::AbstractArray{T}, y::Array{T}) where T<:Real

Compute the gradient of the Q-function with respect to the observation model parameters (C and log_d) for a Poisson Linear Dynamical System.

Arguments

  • grad::AbstractVector{T}: Pre-allocated vector to store the computed gradient.
  • C::AbstractMatrix{T}: The observation matrix. Dimensions: (obsdim, latentdim)
  • log_d::AbstractVector{T}: The log of the baseline firing rates. Dimensions: (obs_dim,)
  • E_z::AbstractArray{T}: The expected latent states. Dimensions: (latentdim, tsteps, n_trials)
  • P_smooth::AbstractArray{T}: The smoothed state covariances. Dimensions: (latentdim, latentdim, tsteps, ntrials)
  • y::Array{T}: The observed data. Dimensions: (obsdim, tsteps, N-trials)

Note

This function modifies grad in-place. The gradient is computed for the negative Q-function, as we're minimizing -Q in optimization routines.

StateSpaceDynamics.initialize_FilterSmooth — Method
initialize_FilterSmooth(model::LinearDynamicalSystem, num_obs::Int) -> FilterSmooth{T}

Initialize a FilterSmooth object for a given linear dynamical system model and number of observations.

Arguments

  • model::LinearDynamicalSystem: The linear dynamical system model containing system parameters, including the latent dimensionality (latent_dim).

  • num_obs::Int: The number of observations (time steps) for which to initialize the smoothing filters.

Returns

  • FilterSmooth{T}: A FilterSmooth instance with all fields initialized to zero arrays. The dimensions of the arrays are determined by the number of states (latent_dim) from the model and the specified number of observations (num_obs).

Example

```julia

Assume model is an instance of LinearDynamicalSystem with latent_dim = 10

numobservations = 100 filtersmooth = initializeFilterSmooth(model, numobservations)

filter_smooth now contains zero-initialized arrays for smoothing operations

StateSpaceDynamics.loglikelihood — Method
loglikelihood(x::Array{T, 3}, lds::LinearDynamicalSystem{S,O}, y::Array{T, 3}) where {T<:Real, S<:GaussianStateModel{T}, O<:PoissonObservationModel{T}}

Calculate the complete-data log-likelihood of a Poisson Linear Dynamical System model for multiple trials.

Arguments

  • x::Array{T, 3}: The latent state variables. Dimensions: (latentdim, TSteps, n_trials)
  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System model.
  • y::Array{T, 3}: The observed data. Dimensions: (obsdim, Tsteps, n_trials)

Returns

  • ll::T: The log-likelihood value.

Examples

lds = PoissonLDS(obs_dim=4, latent_dim=3)
x, y = sample(lds, 100, 10)  # 10 trials, 100 time steps each
ll = loglikelihood(x, lds, y)
StateSpaceDynamics.loglikelihood — Method
loglikelihood(x::Matrix{T}, lds::LinearDynamicalSystem{S,O}, y::Matrix{T}) where {T<:Real, S<:GaussianStateModel, O<:PoissonObservationModel}

Calculate the complete-data log-likelihood of a Poisson Linear Dynamical System model for a single trial.

Arguments

  • x::Matrix{T}: The latent state variables. Dimensions: (latentdim, Tsteps)
  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System model.
  • y::Matrix{T}: The observed data. Dimensions: (obsdim, Tsteps)
  • w::Vector{T}: Weights for each observation in the log-likelihood calculation. Not currently used.

Returns

  • ll::T: The log-likelihood value.

Examples

lds = PoissonLDS(obs_dim=4, latent_dim=3)
x, y = sample(lds, 100, 1)  # 1 trial, 100 time steps
ll = loglikelihood(x, lds, y)
StateSpaceDynamics.loglikelihood — Method
loglikelihood(x::AbstractMatrix{T}, lds::LinearDynamicalSystem{S,O}, y::AbstractMatrix{T}) where {T<:Real, S<:GaussianStateModel{T}, O<:GaussianObservationModel{T}}

Calculate the complete-data log-likelihood of a linear dynamical system (LDS) given the observed data.

Arguments

  • x::AbstractMatrix{T}: The state sequence of the LDS.
  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System.
  • y::AbstractMatrix{T}: The observed data.
  • w::Vector{Float64}: coeffcients to weight the data.

Returns

  • ll::T: The complete-data log-likelihood of the LDS.
StateSpaceDynamics.mstep! — Method
mstep!(lds::LinearDynamicalSystem{S,O}, E_z::Array{T,3}, E_zz::Array{T,4}, E_zz_prev::Array{T,4}, p_smooth::Array{T, 4}, y::Array{T,3}) where {T<:Real, S<:GaussianStateModel{T}, O<:GaussianObservationModel{T}}

Perform the M-step of the EM algorithm for a Linear Dynamical System with multi-trial data.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System struct
  • E_z::Array{T,3}: Expected latent states, size (statedim, statedim, Tsteps, ntrials)
  • E_zz::Array{T,4}: Expected zt * zt', size (statedim, statedim, Tsteps, ntrials)
  • E_zz_prev::Array{T,4}: Expected zt * z{t-1}', size (statedim, statedim, Tsteps, ntrials)
  • p_smooth::Array{T,4}: Smoothed state covariances, size (statedim, statedim, Tsteps, ntrials) (not used)
  • y::Array{T,3}: Observed data, size (obsdim, Tsteps, n_trials)

Note

  • This function modifies lds in-place by updating all model parameters.
  • Updates are performed only for parameters where the corresponding fit_bool is true.
  • All update functions now handle multi-trial data.
  • P_smooth is required but not used in the M-step so that the function signature matches the PoissonLDS version.
StateSpaceDynamics.mstep! — Method
mstep!(plds::LinearDynamicalSystem{S,O}, E_z::Array{T,3}, E_zz::Array{T,4}, E_zz_Prev{T,4}, p_smooth{T,4}, y::Array{T,3}) where {T<:Real, S<:GaussianStateModel{T}, O<:PoissonObservationModel{T}}

Perform the M-step of the EM algorithm for a Poisson Linear Dynamical System with multi-trial data.

Arguments

  • plds::LinearDynamicalSystem{S,O}: The Poisson Linear Dynamical System struct.
  • E_z::Array{T,3}: Expected latent states, size (statedim, statedim, Tsteps, ntrials)
  • E_zz::Array{T,4}: Expected zt * zt', size (statedim, statedim, Tsteps, ntrials)
  • E_zz_prev::Array{T,4}: Expected zt * z{t-1}', size (statedim, statedim, Tsteps, ntrials)
  • p_smooth::Array{T,4}: Smoothed state covariances, size (statedim, statedim, Tsteps, ntrials)
  • y::Array{T,3}: Observed data, size (obsdim, Tsteps, n_trials)

Note

  • This function modifies plds in-place by updating all model parameters.
StateSpaceDynamics.obsparams — Method
obsparams(lds::LinearDynamicalSystem{S,O}) where {S<:AbstractStateModel,O<:AbstractObservationModel}

Extract the observation parameters from a Linear Dynamical System.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System.

Returns

  • params::Vector{Vector{Real}}: Vector of observation parameters.
StateSpaceDynamics.sample — Method
sample(lds::LinearDynamicalSystem{S,O}, T_steps::Int, n_trials::Int) where {T<:Real, S<:GaussianStateModel{T}, O<:GaussianObservationModel{T}}

Sample from a Linear Dynamical System (LDS) model for multiple trials.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System model.
  • T_steps::Int: The number of time steps to sample for each trial.
  • n_trials::Int: The number of trials to sample.=

Returns

  • x::Array{T, 3}: The latent state variables. Dimensions: (latentdim, TSteps, n_trials)
  • y::Array{T, 3}: The observed data. Dimensions: (obsdim, Tsteps, n_trials)

Examples

lds = GaussianLDS(obs_dim=4, latent_dim=3)
x, y = sample(lds, 10, 100)  # 10 trials, 100 time steps each
StateSpaceDynamics.sample — Method
sample(lds::LinearDynamicalSystem{S,O}, T_steps::Int, n_trials::Int) where {T<:Real, S<:GaussianStateModel{T}, O<:PoissonObservationModel{T}}

Sample from a Poisson Linear Dynamical System (LDS) model for multiple trials.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System model.
  • T_steps::Int: The number of time steps to sample for each trial.
  • n_trials::Int: The number of trials to sample.

Returns

  • x::Array{T, 3}: The latent state variables. Dimensions: (latentdim, TSteps, n_trials)
  • y::Array{Int, 3}: The observed data. Dimensions: (obsdim, Tsteps, n_trials)

Examples

lds = LinearDynamicalSystem(obs_dim=4, latent_dim=3)
x, y = sample(lds, 100, 10)  # 10 trials, 100 time steps each
StateSpaceDynamics.smooth — Method
smooth(lds::LinearDynamicalSystem{S,O}, y::Array{T,3}) where {T<:Real, S<:GaussianStateModel{T}, O<:GaussianObservationModel{T}}

This function performs direct smoothing for a linear dynamical system (LDS) given the system parameters and the observed data for multiple trials.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The LDS object representing the system parameters.
  • y::Array{T,3}: The observed data array with dimensions (obsdim, tiemsteps, n_trials).

Returns

  • x::Array{T,3}: The optimal state estimates with dimensions (ntrials, timesteps, latent_dim).
  • p_smooth::Array{T,4}: The posterior covariance matrices with dimensions (latentdim, latentdim, timesteps, ntrials).
  • inverse_offdiag::Array{T,4}: The inverse off-diagonal matrices with dimensions (latentdim, latentdim, timesteps, ntrials).

Example

lds = GaussianLDS(obs_dim=4, latent_dim=3)
y = randn(5, 100, 4)  # 5 trials, 100 time steps, 4 observed dimension
x, p_smooth, inverse_offdiag = smooth(lds, y)
StateSpaceDynamics.smooth — Method
smooth(lds::LinearDynamicalSystem{S,O}, y::Matrix{T}) where {T<:Real, S<:GaussianStateModel{T}, O<:GaussianObservationModel{T}}

This function performs direct smoothing for a linear dynamical system (LDS) given the system parameters and the observed data.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The LDS object representing the system parameters.
  • y::Matrix{T}: The observed data matrix.
  • w::Vector{Float64}: coeffcients to weight the data.

Returns

  • x::Matrix{T}: The optimal state estimate.
  • p_smooth::Array{T, 3}: The posterior covariance matrix.
  • inverse_offdiag::Array{T, 3}: The inverse off-diagonal matrix.
  • Q_val::T: The Q-function value.

Example

lds = GaussianLDS(obs_dim=4, latent_dim=3)
y = randn(100, 4)  # 100 time steps, 4 observed dimensions
x, p_smooth, inverse_offdiag, Q_val = DirectSmoother(lds, y)
StateSpaceDynamics.stateparams — Method
stateparams(lds::LinearDynamicalSystem{S,O}) where {S<:AbstractStateModel,O<:AbstractObservationModel}

Extract the state parameters from a Linear Dynamical System.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System.

Returns

  • params::Vector{Vector{Real}}: Vector of state parameters.
StateSpaceDynamics.sufficient_statistics — Method
sufficient_statistics(x_smooth::Array{T,3}, p_smooth::Array{T,4}, p_smooth_t1::Array{T,4}) where T <: Real

Compute sufficient statistics for the EM algorithm in a Linear Dynamical System.

Arguments

  • x_smooth::Array{T,3}: Smoothed state estimates, size (statedim, statedim, Tsteps, ntrials)
  • p_smooth::Array{T,4}: Smoothed state covariances, size (statedim, statedim, Tsteps, ntrials, state_dim)
  • p_smooth_t1::Array{T,4}: Lag-one covariance smoother, size (statedim, statedim, Tsteps, ntrials, state_dim)

Returns

  • E_z::Array{T,3}: Expected latent states, size (statedim, statedim, Tsteps, ntrials)
  • E_zz::Array{T,4}: Expected zt * zt', size (statedim, statedim, Tsteps, ntrials, state_dim)
  • E_zz_prev::Array{T,4}: Expected zt * z{t-1}', size (statedim, statedim, Tsteps, ntrials, state_dim)

Note

  • The function computes the expected values for all trials.
  • For single-trial data, use inputs with n_trials = 1.
StateSpaceDynamics.update_C! — Method
update_C!(lds::LinearDynamicalSystem{S,O}, E_z::Array{T,3}, E_zz::Array{T,4}, y::Array{T,3}) where {T<:Real, S<:GaussianStateModel{T}, O<:GaussianObservationModel{T}}

Update the observation matrix C of the Linear Dynamical System.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System struct
  • E_z::Array{T,3}: Expected latent states, size (statedim, statedim, Tsteps, ntrials)
  • E_zz::Array{T,4}: Expected zt * zt', size (statedim, statedim, Tsteps, ntrials)
  • y::Array{T,3}: Observed data, size (obsdim, Tsteps, n_trials)

Note

  • This function modifies lds in-place.
  • The update is only performed if lds.fit_bool[5] is true.
  • The result is averaged across all trials.
StateSpaceDynamics.update_Q! — Method
update_Q!(lds::LinearDynamicalSystem{S,O}, E_zz::Array{T, 4}, E_zz_prev::Array{T, 4}) where {T<:Real, S<:GaussianStateModel{T}, O<:AbstractObservationModel{T}}

Update the process noise covariance matrix Q of the Linear Dynamical System.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System struct
  • E_zz::Array{T, 4}: Expected zt * zt', size (statedim, statedim, Tsteps, ntrials)
  • E_zz_prev::Array{T, 4}: Expected zt * z{t-1}', size (statedim, statedim, Tsteps, ntrials)

Note

  • This function modifies lds in-place.
  • The update is only performed if lds.fit_bool[4] is true.
  • The result is averaged across all trials.
StateSpaceDynamics.update_R! — Method
update_R!(lds::LinearDynamicalSystem{S,O}, E_z::Array{T,3}, E_zz::Array{T,4}, y::Array{T,3}) where {T<:Real, S<:GaussianStateModel{T}, O<:GaussianObservationModel{T}}

Update the observation noise covariance matrix R of the Linear Dynamical System.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System struct
  • E_z::Array{T,3}: Expected latent states, size (statedim, statedim, Tsteps, ntrials)
  • E_zz::Array{T,4}: Expected zt * zt', size (statedim, statedim, Tsteps, ntrials)
  • y::Array{T,3}: Observed data, size (obsdim, Tsteps, n_trials)

Note

  • This function modifies lds in-place.
  • The update is only performed if lds.fit_bool[6] is true.
  • The result is averaged across all trials.
StateSpaceDynamics.update_initial_state_covariance! — Method
update_initial_state_covariance!(lds::LinearDynamicalSystem{S,O}, E_z::Array{T,3}, E_zz::Array{T,4}) where {T<:Real, S<:GaussianStateModel{T}, O<:AbstractObservationModel{T}}

Update the initial state covariance of the Linear Dynamical System using the average across all trials.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System struct
  • E_z::Array{T,3}: Expected latent states, size (statedim, statedim, Tsteps, ntrials)
  • E_zz::Array{T,4}: Expected zt * zt', size (statedim, statedim, Tsteps, ntrials, state_dim)

Note

  • This function modifies lds in-place.
  • The update is only performed if lds.fit_bool[2] is true.
  • The initial state covariance is computed as the average of the first time step across all trials.
StateSpaceDynamics.update_initial_state_mean! — Method
update_initial_state_mean!(lds::LinearDynamicalSystem{S,O}, E_z::Array{T,3}) where {T<:Real, S<:GaussianStateModel{T}, O<:AbstractObservationModel{T}}

Update the initial state mean of the Linear Dynamical System using the average across all trials.

Arguments

  • lds::LinearDynamicalSystem{S,O}: The Linear Dynamical System struct
  • E_z::Array{T,3}: Expected latent states, size (statedim, statedim, Tsteps, ntrials)

Note

  • This function modifies lds in-place.
  • The update is only performed if lds.fit_bool[1] is true.
  • The initial state mean is computed as the average of the first time step across all trials.
StateSpaceDynamics.update_observation_model! — Method
update_observation_model!(plds::LinearDynamicalSystem{S,O}, E_z::Array{T, 3}, P_smooth::Array{T, 4}, y::Array{T, 3}) where {T<:Real, S<:GaussianStateModel{T}, O<:PoissonObservationModel{T}}

Update the observation model parameters of a Poisson Linear Dynamical System using gradient-based optimization.

Arguments

  • plds::LinearDynamicalSystem{S,O}: The Poisson Linear Dynamical System model.
  • E_z::Array{T, 3}: The expected latent states. Dimensions: (latentdim, TSteps, n_trials)
  • P_smooth::Array{T, 4}: The smoothed state covariances. Dimensions: (latentdim, TSteps, ntrials, latentdim)
  • y::Array{T, 3}: The observed data. Dimensions: (obsdim, Tsteps, n_trials)

Note

This function modifies plds in-place by updating the observation model parameters (C and logd). The optimization is performed only if `plds.fitbool[5]` is true.