Function gpfa_e_step 

fn gpfa_e_step(
    y: &[f64],
    c: &[f64],
    d: &[f64],
    r_diag: &[f64],
    k_all: &[Vec<f64>],
    n_neurons: usize,
    n_bins: usize,
    n_latents: usize,
) -> (Vec<f64>, Vec<f64>)

E-step: joint Gaussian posterior p(x|y) over all latents and time points.

The posterior precision M (see gpfa_precision) is Cholesky-factored once; the same factor yields the posterior mean (M⁻¹ applied to AᵀR⁻¹y) and the posterior covariance (M⁻¹). Working on the n_state-dimensional state avoids the dense (n_neurons·n_bins)² solve of the naive marginal form.