Geodesic Monte Carlo on Embedded Manifolds

Low-dimensional encoding of decisions in parietal cortex reflects long-term training history

Neurons in parietal cortex exhibit task-related activity during decision-making tasks. However, it remains unclear how long-term training to perform different tasks over months or even years shapes neural computations and representations. We examine lateral intraparietal area (LIP) responses during a visual motion delayed-match-to-category task. We consider two pairs of male macaque monkeys with different training histories: one trained only on the categorization task, and another first trained to perform fine motion-direction discrimination (i.e., pretrained). We introduce a novel analytical approach-generalized multilinear models-to quantify low-dimensional, task-relevant components in population activity. During the categorization task, we found stronger cosine-like motion-direction tuning in the pretrained monkeys than in the category-only monkeys, and that the pretrained monkeys' performance depended more heavily on fine discrimination between sample and test stimuli. These results suggest that sensory representations in LIP depend on the sequence of tasks that the animals have learned, underscoring the importance of considering training history in studies with complex behavioral tasks.

Principled, practical, flexible, fast: a new approach to phylogenetic factor analysis

Biological phenotypes are products of complex evolutionary processes in which selective forces influence multiple biological trait measurements in unknown ways. Phylogenetic comparative methods seek to disentangle these relationships across the evolutionary history of a group of organisms. Unfortunately, most existing methods fail to accommodate high-dimensional data with dozens or even thousands of observations per taxon. Phylogenetic factor analysis offers a solution to the challenge of dimensionality. However, scientists seeking to employ this modeling framework confront numerous modeling and implementation decisions, the details of which pose computational and replicability challenges.We develop new inference techniques that increase both the computational efficiency and modeling flexibility of phylogenetic factor analysis. To facilitate adoption of these new methods, we present a practical analysis plan that guides researchers through the web of complex modeling decisions. We codify this analysis plan in an automated pipeline that distills the potentially overwhelming array of decisions into a small handful of (typically binary) choices.We demonstrate the utility of these methods and analysis plan in four real-world problems of varying scales. Specifically, we study floral phenotype and pollination in columbines, domestication in industrial yeast, life history in mammals, and brain morphology in New World monkeys.General and impactful community employment of these methods requires a data scientific analysis plan that balances flexibility, speed and ease of use, while minimizing model and algorithm tuning. Even in the presence of non-trivial phylogenetic model constraints, we show that one may analytically address latent factor uncertainty in a way that (a) aids model flexibility, (b) accelerates computation (by as much as 500-fold) and (c) decreases required tuning. These efforts coalesce to create an accessible Bayesian approach to high-dimensional phylogenetic comparative methods on large trees.

Nonparametric Bayesian Regression and Classification on Manifolds, With Applications to 3D Cochlear Shapes

Advanced shape analysis studies such as regression and classification need to be performed on curved manifolds, where often, there is a lack of standard statistical formulations. To overcome these limitations, we introduce a novel machine-learning method on the shape space of curves that avoids direct inference on infinite-dimensional spaces and instead performs Bayesian inference with spherical Gaussian processes decomposition. As an application, we study the shape of the cochlear spiral-shaped cavity within the petrous part of the temporal bone. This problem is particularly challenging due to the relationship between shape and gender, especially in children. Experimental results for both synthetic and real data show improved performance compared to state-of-the-art methods.

Evaluating the Implicit Midpoint Integrator for Riemannian Manifold Hamiltonian Monte Carlo

Riemannian manifold Hamiltonian Monte Carlo is traditionally carried out using the generalized leapfrog integrator. However, this integrator is not the only choice and other integrators yielding valid Markov chain transition operators may be considered. In this work, we examine the implicit midpoint integrator as an alternative to the generalized leapfrog integrator. We discuss advantages and disadvantages of the implicit midpoint integrator for Hamiltonian Monte Carlo, its theoretical properties, and an empirical assessment of the critical attributes of such an integrator for Hamiltonian Monte Carlo: energy conservation, volume preservation, and reversibility. Empirically, we find that while leapfrog iterations are faster, the implicit midpoint integrator has better energy conservation, leading to higher acceptance rates, as well as better conservation of volume and better reversibility, arguably yielding a more accurate sampling procedure.

Flexible Bayesian Dynamic Modeling of Correlation and Covariance Matrices

Modeling correlation (and covariance) matrices can be challenging due to the positive-definiteness constraint and potential high-dimensionality. Our approach is to decompose the covariance matrix into the correlation and variance matrices and propose a novel Bayesian framework based on modeling the correlations as products of unit vectors. By specifying a wide range of distributions on a sphere (e.g. the squared-Dirichlet distribution), the proposed approach induces flexible prior distributions for covariance matrices (that go beyond the commonly used inverse-Wishart prior). For modeling real-life spatio-temporal processes with complex dependence structures, we extend our method to dynamic cases and introduce unit-vector Gaussian process priors in order to capture the evolution of correlation among components of a multivariate time series. To handle the intractability of the resulting posterior, we introduce the adaptive Δ-Spherical Hamiltonian Monte Carlo. We demonstrate the validity and flexibility of our proposed framework in a simulation study of periodic processes and an analysis of rat's local field potential activity in a complex sequence memory task.

Simulating sticky particles: A Monte Carlo method to sample a stratification

Many problems in materials science and biology involve particles interacting with strong, short-ranged bonds that can break and form on experimental timescales. Treating such bonds as constraints can significantly speed up sampling their equilibrium distribution, and there are several methods to sample probability distributions subject to fixed constraints. We introduce a Monte Carlo method to handle the case when constraints can break and form. More generally, the method samples a probability distribution on a stratification: a collection of manifolds of different dimensions, where the lower-dimensional manifolds lie on the boundaries of the higher-dimensional manifolds. We show several applications of the method in polymer physics, self-assembly of colloids, and volume calculation in high dimensions.

Bayesian adaptation of chaos representations using variational inference and sampling on geodesics

A novel approach is presented for constructing polynomial chaos representations of scalar quantities of interest (QoI) that extends previously developed methods for adaptation in Homogeneous Chaos spaces. In this work, we develop a Bayesian formulation of the problem that characterizes the posterior distributions of the series coefficients and the adaptation rotation matrix acting on the Gaussian input variables. The adaptation matrix is thus construed as a new parameter of the map from input to QoI, estimated through Bayesian inference. For the computation of the coefficients' posterior distribution, we use a variational inference approach that approximates the posterior with a member of the same exponential family as the prior, such that it minimizes a Kullback-Leibler criterion. On the other hand, the posterior distribution of the rotation matrix is explored by employing a Geodesic Monte Carlo sampling approach, consisting of a variation of the Hamiltonian Monte Carlo algorithm for embedded manifolds, in our case, the Stiefel manifold of orthonormal matrices. The performance of our method is demonstrated through a series of numerical examples, including the problem of multiphase flow in heterogeneous porous media.

Intrinsic stochastic differential equations as jets

We explain how Itô stochastic differential equations (SDEs) on manifolds may be defined using 2-jets of smooth functions. We show how this relationship can be interpreted in terms of a convergent numerical scheme. We also show how jets can be used to derive graphical representations of Itô SDEs, and we show how jets can be used to derive the differential operators associated with SDEs in a coordinate-free manner. We relate jets to vector flows, giving a geometric interpretation of the Itô–Stratonovich transformation. We show how percentiles can be used to give an alternative coordinate-free interpretation of the coefficients of one-dimensional SDEs. We relate this to the jet approach. This allows us to interpret the coefficients of SDEs in terms of ‘fan diagrams’. In particular, the median of an SDE solution is associated with the drift of the SDE in Stratonovich form for small times.

Bayesian inference of elastic properties with resonant ultrasound spectroscopy

Bayesian modeling and Hamiltonian Monte Carlo (HMC) are utilized to formulate a robust algorithm capable of simultaneously estimating anisotropic elastic properties and crystallographic orientation of a specimen from a list of measured resonance frequencies collected via Resonance Ultrasound Spectroscopy (RUS). Unlike typical optimization procedures which yield point estimates of the unknown parameters, computing a Bayesian posterior yields probability distributions for the unknown parameters, and HMC is an efficient way to compute this posterior. The algorithms described are demonstrated on RUS data collected from two parallelepiped specimens of structural metal alloys. First, the elastic constants for a specimen of fine-grain polycrystalline Ti-6Al-4 V with random crystallographic texture and isotropic elastic symmetry are estimated. Second, the elastic constants and crystallographic orientation for a single crystal Ni-based superalloy CMSX-4 specimen are accurately determined, using only measurements of the specimen geometry, mass, and resonance frequencies. The unique contributions of this paper are as follows: the application of HMC for sampling the Bayesian posterior of a probabilistic RUS model, and the procedure for simultaneous estimation of elastic constants and lattice-specimen misorientation. Compared to previous approaches these algorithms demonstrate superior convergence behavior, particularly when the initial parameterization is unknown, and enable substantially simplified experimental procedures.

Geodesic Lagrangian Monte Carlo over the space of positive definite matrices: with application to Bayesian spectral density estimation

We present geodesic Lagrangian Monte Carlo, an extension of Hamiltonian Monte Carlo for sampling from posterior distributions defined on general Riemannian manifolds. We apply this new algorithm to Bayesian inference on symmetric or Hermitian positive definite matrices. To do so, we exploit the Riemannian structure induced by Cartan's canonical metric. The geodesics that correspond to this metric are available in closed-form and-within the context of Lagrangian Monte Carlo-provide a principled way to travel around the space of positive definite matrices. Our method improves Bayesian inference on such matrices by allowing for a broad range of priors, so we are not limited to conjugate priors only. In the context of spectral density estimation, we use the (non-conjugate) complex reference prior as an example modeling option made available by the algorithm. Results based on simulated and real-world multivariate time series are presented in this context, and future directions are outlined.

A Bayesian supervised dual-dimensionality reduction model for simultaneous decoding of LFP and spike train signals

Neuroscientists are increasingly collecting multimodal data during experiments and observational studies. Different data modalities-such as EEG, fMRI, LFP, and spike trains-offer different views of the complex systems contributing to neural phenomena. Here, we focus on joint modeling of LFP and spike train data, and present a novel Bayesian method for neural decoding to infer behavioral and experimental conditions. This model performs supervised dual-dimensionality reduction: it learns low-dimensional representations of two different sources of information that not only explain variation in the input data itself, but also predict extra-neuronal outcomes. Despite being one probabilistic unit, the model consists of multiple modules: exponential PCA and wavelet PCA are used for dimensionality reduction in the spike train and LFP modules, respectively; these modules simultaneously interface with a Bayesian binary regression module. We demonstrate how this model may be used for prediction, parametric inference, and identification of influential predictors. In prediction, the hierarchical model outperforms other models trained on LFP alone, spike train alone, and combined LFP and spike train data. We compare two methods for modeling the loading matrix and find them to perform similarly. Finally, model parameters and their posterior distributions yield scientific insights.

Gradient-based MCMC samplers for dynamic causal modelling

In this technical note, we derive two MCMC (Markov chain Monte Carlo) samplers for dynamic causal models (DCMs). Specifically, we use (a) Hamiltonian MCMC (HMC-E) where sampling is simulated using Hamilton's equation of motion and (b) Langevin Monte Carlo algorithm (LMC-R and LMC-E) that simulates the Langevin diffusion of samples using gradients either on a Euclidean (E) or on a Riemannian (R) manifold. While LMC-R requires minimal tuning, the implementation of HMC-E is heavily dependent on its tuning parameters. These parameters are therefore optimised by learning a Gaussian process model of the time-normalised sample correlation matrix. This allows one to formulate an objective function that balances tuning parameter exploration and exploitation, furnishing an intervention-free inference scheme. Using neural mass models (NMMs)-a class of biophysically motivated DCMs-we find that HMC-E is statistically more efficient than LMC-R (with a Riemannian metric); yet both gradient-based samplers are far superior to the random walk Metropolis algorithm, which proves inadequate to steer away from dynamical instability.

Spherical Hamiltonian Monte Carlo for Constrained Target Distributions

Statistical models with constrained probability distributions are abundant in machine learning. Some examples include regression models with norm constraints (e.g., Lasso), probit models, many copula models, and Latent Dirichlet Allocation (LDA) models. Bayesian inference involving probability distributions confined to constrained domains could be quite challenging for commonly used sampling algorithms. For such problems, we propose a novel Markov Chain Monte Carlo (MCMC) method that provides a general and computationally efficient framework for handling boundary conditions. Our method first maps the D-dimensional constrained domain of parameters to the unit ball [Formula: see text], then augments it to a D-dimensional sphere SD such that the original boundary corresponds to the equator of SD . This way, our method handles the constraints implicitly by moving freely on the sphere generating proposals that remain within boundaries when mapped back to the original space. To improve the computational efficiency of our algorithm, we divide the dynamics into several parts such that the resulting split dynamics has a partial analytical solution as a geodesic flow on the sphere. We apply our method to several examples including truncated Gaussian, Bayesian Lasso, Bayesian bridge regression, and a copula model for identifying synchrony among multiple neurons. Our results show that the proposed method can provide a natural and efficient framework for handling several types of constraints on target distributions.