Gaussian-process factor analysis (GPFA) is a dimensionality reduction method [Yu2009] that extracts smooth, low-dimensional latent trajectories from noisy, high-dimensional time series data. GPFA applies factor analysis (FA) to observed data to simultaneously reduce their dimensionality and smooth the resulting low-dimensional trajectories by fitting a Gaussian process (GP) model to them. See the below :ref:`notes` for its main differences to the popular PCA dimensionality redunction method.
This documentation describes BlockInvGPFA, a Python implementation of the GPFA model that incorporates several performance and usability improvements over existing implementations (e.g., from the Elephant toolbox [Denker2018]).
The essential ingredients of the generative model underlying GPFA are observations \boldsymbol{X} and latent variables \boldsymbol{Z}. For each high-dimensional time series in a dataset, the x_\textrm{dim}-dimensional time series is split into T time bins of size \delta t, leading to the x_\textrm{dim} \times T matrix \boldsymbol{X} that contains these data. The observations in each time bin are assumed to be stochastically generated by a z_\textrm{dim}-dimensional latent variable (z_\textrm{dim} < x_\textrm{dim}) through some noisy linear mapping. Across time, these latent variables are collected in the z_\textrm{dim} \times T matrix \boldsymbol{Z}. A dataset usually contains a set of observations \boldsymbol{X}'s that each have their associated latent variabes \boldsymbol{Z}, and that all have the same observation and latent dimensions, x_\textrm{dim} and z_\textrm{dim}, but might differ in their number of time bins T.
GPFA's generative model is defined by the emission and the latent time-course models:
Emission model, p(\boldsymbol{X} | \boldsymbol{Z}): GPFA's emission model assumes that, in each time bin t, the observations \boldsymbol{x}_{:,t} are generated independently through a stochastic linear mapping from the corresponding latent variables \boldsymbol{z}_{:,t},
\boldsymbol{x}_{:,t} | \boldsymbol{z}_{:,t} \sim \mathcal{N} \left( \boldsymbol{C} \boldsymbol{z}_{:,t} + \boldsymbol{d}, \boldsymbol{R} \right) ,where \boldsymbol{x}_{:,t} and \boldsymbol{z}_{:,t} are the t th column (i.e., the t th time bin) of \boldsymbol{X} and \boldsymbol{Z}, respectively, \mathcal{N}\left( \boldsymbol{\mu}, \boldsymbol{\Sigma}\right) denotes the multivariate Gaussian with mean \boldsymbol{\mu} and covariance \boldsymbol{\Sigma}, and the loading matrix \boldsymbol{C} (x_\textrm{dim} \times z_\textrm{dim}), bias vector \boldsymbol{d} (x_\textrm{dim}), and diagonal (as in FA) emission covariance matrix \boldsymbol{R} (x_\textrm{dim} \times x_\textrm{dim}) are model parameters.
Latent time-course model, p(\boldsymbol{Z}): GPFA assumes that the latent variables evolve in each latent dimension i (or factor) independently according to a Gaussian process,
\boldsymbol{z}_{i,:} \sim \mathcal{GP} \left( \boldsymbol{0}, \boldsymbol{K} \right) ,where \boldsymbol{z}_{i,:} is the i th row (i.e., i th latent dimension/factor) of \boldsymbol{Z}, and \mathcal{GP} \left( \boldsymbol{0}, \boldsymbol{K} \right) denotes a zero-mean Gaussian process with covariance kernel \boldsymbol{K}.
Overall, the model has emission model parameters \boldsymbol{C}, \boldsymbol{d}, and \boldsymbol{R} and time-course model parameters that are the parameters of the chosen GP kernel \boldsymbol{K}. When GPFA is fit to a set of observations, \boldsymbol{X}_1, \boldsymbol{X}_2, \dots, GPFA adjusts these parameters and associated latent variables to best explain these observations.
GPFA differs from the popular PCA in two significant ways:
- It performs dimensionality reduction by factor analysis, which supports different residual variances in different data dimensions, whereas PCA assumed this residual variance to be the same across all dimensions.
- PCA does not assume any temperal dependencies between different latent variables. GPFA, in contrast assumes latent trajecories to be smooth across time, and models this assumption by a Gaussian process.
This example illustrates the application of BlockInvGPFA to data analysis by encompassing synthetic data generation, BlockInvGPFA model fitting, and latent variable extraction. It emphasizes a step-by-step explanation of utilizing BlockInvGPFA for data analysis, highlighting key functionalities of the library.
import numpy as np
from blockinvgpfa import BlockInvGPFA
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, RBF, WhiteKernel
# Simulation parameters
rng_seeds = [0, 42]
z_dim = 2
x_dim = 10
tau_f = 0.6
sigma_n = 0.001
sigma_f = 1 - sigma_n
bin_size = 0.05 # [s]
num_obs = len(rng_seeds)
T_per_obs = 400
kernel = ConstantKernel(sigma_f, constant_value_bounds='fixed') * RBF(length_scale=tau_f) + \
ConstantKernel(sigma_n, constant_value_bounds='fixed') * WhiteKernel(noise_level=1, noise_level_bounds='fixed')
tsdt = np.arange(0, T_per_obs) * bin_size
np.random.seed(2)
C = np.random.uniform(0, 2, (x_dim, z_dim)) # loading matrix
sqrtR = np.random.uniform(0, 0.5, x_dim)
# Generate latent and observation sequences in line with GPFA model
X = []
Z = []
for n in range(num_obs):
np.random.seed(rng_seeds[n])
# Sample latent sequence according to GPFA latent time-course model
gp_model = GaussianProcessRegressor(kernel=kernel)
z = gp_model.sample_y(tsdt[:, np.newaxis], n_samples=z_dim, random_state=rng_seeds[n]).T
Z.append(z)
# Sample observations according to GPFA emission model
x = C @ z + np.random.normal(0, sqrtR[:, np.newaxis] , (x_dim, T_per_obs))
X.append(x)
blockinv_gpfa = BlockInvGPFA(bin_size=bin_size, z_dim=z_dim, em_tol=1e-3, verbose=True)
blockinv_gpfa.fit(X) # this will take a while
results = blockinv_gpfa.predict(returned_data=['Z_mu_orth', 'Z_mu'])
print(blockinv_gpfa.variance_explained())Initializing parameters using factor analysis...
Fitting GP parameters by EM...
Fitting has converged after 135 EM iterations.
(0.9834339959622203, array([0.80069798, 0.18273602]))Note that, due to differences in the linear algebra packages used on different machine architectures and operating systems, the above results might differ across machines.
Let us now visualize the true orthogonalized latents against the inferred orthogonalized latents.
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rcParams.update({'font.size': 20})
mpl.rcParams.update({'lines.linewidth': 3})
# orthogonalize true latents using inferred C
true_Z = Z
true_Z_orth = [blockinv_gpfa.OrthTrans_ @ Zn for Zn in true_Z]
# Extract Z_orth and pZ_mu from results
Z_orth = results[0]['Z_mu_orth']
Z_mu = results[0]['Z_mu']
fig, axs = plt.subplots(2, 1, figsize=(13, 13))
# plot first (orthonormalized) latent from first sequence (flipping the sign
# of the inferred latent, as this is arbitrary and happens to be wrongly inferred)
axs[0].plot(tsdt, true_Z_orth[0][0,:], label='true $Z_{1,:}$ (orth)', c='blue')
axs[0].plot(tsdt, -Z_orth[0][0,:], label='inf $Z_{1,:}$ (orth)', c='orange')
axs[0].legend()
axs[0].set_title('Comparing inferred to true latents')
axs[0].set_ylabel('$Z_{1,:}$ from first sequence')
# plot first latent from second sequence
axs[1].plot(tsdt, true_Z_orth[1][0,:], label='true $Z_{1,:}$ (orth)', c='blue')
axs[1].plot(tsdt, -Z_orth[1][0,:], label='inf $Z_{1,:}$ (orth)', c='orange')
axs[1].legend()
axs[1].set_xlabel('time $t$')
axs[1].set_ylabel('$Z_{1,:}$ from second sequence')
plt.tight_layout()
plt.show()
You can find a more complex example that applies BlockInvGPFA to neural data used by [Even2019] (available at [Even2022]) in an :doc:`neural_example`. To run this example yourself, download the :download:`Jupyter notebook containing the code <neural_example.ipynb>`.
The code was adjusted and extended from the Elephant ephys analysis toolkit which is a direct Python translation of Byron Yu's MATLAB implementation.
| [Denker2018] | Denker, M., Yegenoglu, A., and Gr"un, S. "Collaborative HPC-enabled workflows on the HBP Collaboratory using the Elephant framework." In Neuroinformatics 2018, 2018, p. P19. |
.. toctree::
:maxdepth: 2
:hidden:
blockinv_gpfa_class
preprocessing_class
neural_example
installation
contribute
acknowledgments
authors
citation