Hamiltonian Monte Carlo Inference of Marginalized Linear Mixed-Effects Models

Authors: jinlin  lai
daniel  sheldon
justin  domke
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 15:31:21

Abstract: Bayesian reasoning in linear mixed-effects models (LMMs) is challenging and often requires advanced sampling techniques like Markov chain Monte Carlo (MCMC).
A common approach is to write the model in a probabilistic programming language and then sample via Hamiltonian Monte Carlo (HMC).
However, there are many ways a user can transform a model that make inference more or less efficient.
In particular, marginalizing some variables can greatly improve inference but is difficult for users to do manually.
We develop an algorithm to easily marginalize random effects in LMMs.
A naive approach introduces cubic time operations within an inference algorithm like HMC, but we reduce the running time to linear using fast linear algebra techniques.
We show that marginalization is always beneficial when applicable and highlight improvements in various models, especially ones from cognitive sciences.

Main paper: 
 introduction: 
 Bayesian hierarchical models account for complicated relationships in data by introducing hierarchical structures [23] . Among hierarchical models, linear mixed effects models (LMMs) are widely used in various scientific disciplines, including ecology [31] , medicine [7] , psychology [41] , neuroscience [77] and cognitive science [47] . Solving LMMs involves inferring latent variables, such as fixed and random effects, based on the observed data. Fixed effects are shared by all observations, while random effects vary across different groups within the data. LMMs are often implemented using probabilistic programming languages (PPLs), which isolate inference from modeling: users write a program representing the model and the PPL automatically executes a suitable inference algorithm. Variants of Hamiltonian Monte Carlo (HMC) [15] are dominant in many PPLs today and are widely used for LMMs. For example, BRMS [8] is an influential R package that allows users to write regression-style formulas that are automatically translated to Stan programs [9] representing an LMM, and then Stan's HMC implementation is called to generate posterior samples.We develop techniques that allow users to easily transform their models to analytically marginalize random effect variables from LMMs to improve the efficiency of HMC. Marginalization has several benefits. First, there are often pathologies in LMMs that hinder efficient HMC sampling. A notable one is the "funnel" shape created by correlation between variance parameters and parameters for fixed or random effects [45] . Marginalization [35] and other program transformations [26] have been shown to be useful in addressing such pathologies. Second, marginalization reduces the number H of latent variables for HMC. The complexity of HMC is about O(H 5/4 ) [11, 46] , so it is desirable to run HMC on a subset of variables if marginalization can be done efficiently. Our methods enable marginalization of random effects in LMMs with a linear Gaussian structure, which includes models with normal and log-normal likelihoods as well as other likelihoods for continuous data based on transforming a normal distribution. Note that our methods are not limited to HMC, and could be applied to many inference algorithms.There are several challenges to efficient marginalization. The automatic marginalization algorithm of [35] can be applied to LMMs but is limited to scalar random variables, so it requires users to construct the LMM as a graphical model with separate variables for each effect and observation. Another alternative is to model the relationships between effects and observations with a design matrix and marginalize effects using properties of multivariate normal distributions. We call this the "vectorized approach" since it can leverage vectorization to accelerate computations. Unfortunately, vectorized marginalization leads to a dense covariance matrix over the observations and thus cubic time for evaluating the log-density within HMC, when the log-density of the original could be evaluated in linear time. Our main technical contribution is to accelerate vectorized marginalization for LMMs using fast linear algebra: we show that marginalization for a single random effect can be achieved with linear time complexity and can significantly accelerate HMC compared to both the original model and non-vectorized marginalization.We implement vectorized marginalization for LMMs in NumPyro [5, 54] via simple classes users can use to express their models. We evaluate our approach on a variety of real LMMs from past scientific investigations, including nine models and datasets from cognitive sciences, and find that marginalization is always beneficial. Our findings suggest that practitioners should marginalize group-level effects whenever applicable in Bayesian hierarchical inference. 
background: 
 To motivate our problem, we present an example model. In [72] , a set of experiments were run to examine the relationship between human pupil and attention load. A total of N = 2228 measurements of pupil sizes from M = 20 subjects were taken under different attention load levels. Specifically, in the ith measurement, the pupil size y i ∈ R + of subject g i ∈ {1, 2, ..., k} under attention load c i ∈ {0, 1, 2, 3, 4, 5} was recorded. Pupil size can be assumed to have linear relationship y i ≈ θ 0 + θ 1 c i with respect to the attention load c i , where both the slope θ 1 and intercept θ 0 split into fixed and random effects:y i = α + u gi,1 + c i (β + u gi,2 ) + ϵ, ϵ ∼ N (0, σ 2 ),where α, β are variables for fixed effects and u •,• are variables for subject-specific random effects. Bayesian hierarchical modeling assigns priors to each unknown variable: α ∼ N (1000, 500 2 ), β ∼ N (0, 100), σ ∼ N + (0, 1000), T ∼ N + (0, diag(1000 2 , 1000 2 )),L u ∼ LKJCholesky(2, 1), [u j,1 , u j,2 ] ∼ N (0, TL u L T u T), j = 1, 2, ..., k.A half-normal distribution (N + ) and an LKJ distribution (LKJCholesky) [36] are used as a prior on the covariance matrix. Inference for the unknown parameters determining the relationship between pupil size and attention load can be performed by writing a probabilistic program and running HMC. For example, in NumPyro, the regression model for all measurements may be implemented as below. The code above uses advanced indexing and vectorization techniques in numpy, where u,g,c,y are all vectors or matrices. We further observe that, conditioned on α, β, σ, T, L u , the distribution of all u j and all y i form a multivariate normal distribution. Theoretically it is possible to analytically integrate u out from the model to improve inference efficiency. But it is not straightforward for users to transform the probabilistic program to do so, and, as we will see, if done in the most obvious way, may not make the model more efficient for HMC. To be more clear about how marginalization can be implemented, we rearrange the model into a canonical form that focuses on the random effects. All observations are collected into the vector y = [y 1 , ..., y N ] T and random effects into the vector uu 1 u 2 y 1 y 2 y 3= [u 1,1 , u 1,2 , ..., u k,1 , u k,2 ] T .Then, we can writeu ∼ N (µ, Σ u ), y ∼ N (Au + b, Σ y ),where µ, Σ u , A, b, Σ y are functions of α, β, σ, T, L u , g i , c i . Note that y i only depends on the entry u gi of u. The corresponding graphical model has a tree structure, as demonstrated in Figure 1 . This tree structure has several benefits: first, matrix multiplications like Au and A T y can be done efficiently; second, we will see that it leads to a block-diagonal structure that facilitates efficient inversion in a key matrix that appears later.For more general LMMs with more than one class of random effects we generalize the canonical form asEQUATIONwhere p(Θ) is the distribution for global variables (including fixed effects), p(u i |Θ) is the distribution for random effects and p(y|Θ, u 1 , ..., u L ) is the distribution for observations. Notationally this generalization further adds an index to each random effect to specify its class. A user might specify the model directly in this canonical form, or in another syntax (e.g., the formula syntax of BRMS) that is compiled to this form. Each pair (u i , A i ) specifies a class of random effects for a particular classification of the observations (e.g., by subject, age, gender, etc.). Each classification contains multiple groups and different classifications are distinct from one another. Each observation belongs to one group for each classification. The vectoru i = [u T i,1 , u T i,2 , ..., u T i,ki ]T contains random effects for the ith classification (e.g., subject, age, or gender), consisting of k i groups (e.g., one subject, age, or gender), with u i,j containing the random effects (e.g., slope and intercept) for the jth group. We denote the number of observations as dim(y) = N , and the number of random effects per group as dim(u i,j ) = d. Any covariates-such as c i in the pupil size example-are considered constants and not represented in the notation. In LMMs, the number d is related to the number of covariates and is usually small. The total number of random effects for u i is denoted as dim(u i ) = M i = k i d. The matrix A i therefore has size N × M i , and encodes the group structure for u i by mapping random effects (together with covariates) to observations. Each row of A i encodes the assignment of an observation to one group, so it has at most d nonzero elements. Therefore, the complexity of computing A i u i is O(N d), as A has at most N d nonzero elements. Henceforth, we omit the dependence on Θ for µ, Σ u , A, b, Σ y for simplicity.Marginalizing u i It is possible to analytically marginalize variables in this model: since the mean of y is linear in each u i and all of these variables are normally distributed, the joint distribution of (y, u 1 , . . . , u L ) is also multivariate normal. We will focus for most of the paper on marginalizing the random effects u i for a single i in order to leverage the tree structure mentioned earlier, but return in Section 4 to the idea of marginalizing many effects. Locally, u i and y form the conditional distribution p(u i , y|Θ, u -i ) = p(u i |Θ)p(y|Θ, u -i , u i ). Marginalized MCMC rewrites this conditional distribution as p(u i , y|Θ, u -i ) = p(y|Θ, u -i )p(u i |Θ, y, u -i ), which reverses the dependence between u i and y [35] . During sampling, u i is marginalized from the HMC procedure by using p(y|Θ, u -i ) as the likelihood function and p(Θ, u -i ) as the distribution of latent variables. After HMC sampling, u i is recovered through ancestral sampling from p(u i |Θ, y, u -i ) given posterior samples of (Θ, u -i ). The reversal requires analytical forms of p(y|Θ, u -i ) and p(u i |Θ, y, u -i ), which can be obtained via standard marginalization and conditioning operations on multivariate normal distributions [e.g., 6]EQUATIONwhereM = Σ ui A T i (A i Σ ui A T i + Σ y ) -1 .Marginalization introduces the benefit of sampling in a lower dimensional space, but the cost depends on the complexity of evaluating the log-density functions of these two distributions in order to run HMC. Table 1 : Time complexities of different HMC approaches for the submodel involved in marginalization. Initialization is done once before the HMC loop. The log density is computed within each step of the leapfrog integrator. Recovery is performed for each sample from HMC. N is the number of observations, M is the dimension for one class of random effects, D is the dimension for all classes of random effects, L is the number of classes, d is the dimension for an effect of a group in a class.Approach Initialization Log density Recoveryp(u i , y|Θ, u -i ) HMC - O(M d 2 + N Ld) - Naive marginalization - O(M 3 + N 3 ) O(M 3 + N 3 ) Marginalize with lemmas - O(M d 2 + N Ld + N d 2 ) O(M d 2 + N Ld + N d 2 ) p(v, y|Θ) HMC - O(Dd 2 + N Ld) - Naive marginalization - O(D 3 + N 3 ) O(D 3 + N 3 ) Marginalize with assumptions O(D 3 + N L 2 d 2 ) O(D 2 + N Ld) O(D 2 + N Ld)In practice, the original model usually has structure that makes evaluating its density very efficient, which is lost by naive marginalization. For example, the observations in y are usually conditionally independent, making Σ y diagonal; also, Σ-i ) = -1 2 det (E) -1 2 z T E -1 z + C, where z = y -j̸ =i A j u j -A i µ i -b.E is not diagonal and without using additional structure will trigger O(N 3 ) operations within each step of the leapfrog integrator within HMC. For the recovery distribution p(u i |Θ, y, u -i ), E will be inverted when calculating M. Also, a Cholesky decomposition for the covariance (I -MA i )Σ ui should be computed for sampling, which takes O(M 3 i ) time. These cubic time operations are prohibitively expensive for large datasets. We summarize the complexities of different approaches in Table 1 . In Section 3, we discuss how to marginalize one group of random effects with lemmas from linear algebra. In Section 4, we discuss how to marginalize all random effects with additional assumptions. 
marginalization with fast linear algebra: 
 We now show how to speed up calculations with the marginalized model using fast linear algebra methods. In particular, we use the matrix inversion lemma and matrix determinant lemma together with special structure in the relevant matrices. In this section, we sometimes omit the subscript i such as for A i and Σ ui for simplicity. The steps in log density evaluation and recovery are summarized in Algorithm 1, and in Algorithm 2 in the appendix, with comments about their implementation and cost. We mainly use sparsity and tree-structure in A to make operations faster. As an overview, computing z takes O(N Ld) time for L sparse matrix multiplications of time O(N d) each. Also, evaluating As and A T t both take O(N d) for any s ∈ R M and any t ∈ R N . With tree-structure, we will see thatA T Σ -1y A is block-diagonal and can be computed efficiently.The two main bottlenecks when evaluating log p(y|Θ, u -i ) are computing det(E) and z T E -1 z.With the matrix determinant lemma [32] , we have thatEQUATIONBy the matrix inversion lemma or the Woodbury formula [53] we have thatEQUATIONAlgorithm 1 Evaluating log p(y|Θ, u -i ). Each A i is an N × M i sparse matrix with N d elements and tree structure.Σ y is N × N diagonal. Σ u is M × M block-diagonal with block size d. 1: z = y -j̸ =i A j u j -A i µ i -b ▷ Sparse matrix multiplication in O(N Ld) time 2: F = Σ -1 u + A T Σ -1 y A ▷ Block diagonal computation in O((M + N )d 2 ) time 3: x = A T Σ -1 y z ▷ Sparse matrix multiplication in O(N d) time 4: a = log det(F) + log det(Σ u ) + log det(Σ y ) ▷ Determinants in O(M d 2 ) time 5: b = z T Σ -1 y z -x T F -1 x ▷ Quadratic form in O(N + M d) time 6: return -1 2 (a + b) + CBy using the facts that Σ u is block-diagonal, Σ y is diagonal, and A has N d nonzero elements, the quantities det(Σ u ), det(Σ y ), z T Σ -1 y z, andA T Σ -1 y z can each be calculated in O(M d 2 + N d) time.Equations ( 3) and ( 4) contain the expressionsF -1 or det(F) for the M × M matrix F := Σ -1 u + A T Σ -1 y A, which both require O(M 3) time when done naively. The following theorem shows that these quantities can be computed inO((M + N )d 2 ) for LMMs. Theorem 1. If Σ y is diagonal, Σ u is block-diagonal with blocks of size d × d, then F = Σ -1 u + A T Σ -1 y A is also block-diagonal with d × d blocks and computing A T Σ -1 y A takes O(N d 2 ).Proof. The proof uses the tree-structure in A. For details, see Appendix B.1.Therefore, it is O((M + N )d 2 )to compute det(F) and F -1 . Combined with other parts in the formulas, the overall complexity isO(M d 2 + N Ld + N d 2 ).In LMMs, d is usually small, so the complexity with marginalization can be viewed as the same as the complexity without marginalization.Different from evaluating log p(y|Θ, u -i ), ancestral sampling from p(u i |Θ, y, u -i ) is only performed once for each posterior sample. When sampling from p(u i |Θ, y, u -i ), computing M directly is also costly. With the matrix inversion lemma, we haveEQUATIONWith this expression, the mean variable µ + Mz, then is evaluated in O((M + N )d 2 ), by computing Σ -1 u +A T Σ -1 y A in the same way as Line 2 of Algorithm 1. For the covariance variable (I-MA)Σ u , we have from substituting Equation 5that(I -MA)Σ u = Σ u -Σ u A T Σ -1 y AΣ u + Σ u A T Σ -1 y A(Σ -1 u + A T Σ -1 y A) -1 A T Σ -1 y AΣ u = Σ u (Σ -1 u -A T Σ -1 y A + A T Σ -1 y A(Σ -1 u + A T Σ -1 y A) -1 A T Σ -1 y A)Σ u . Note that Σ -1 u , A T Σ -1 y A, and F = Σ -1 u + A T Σ -1y A are all block diagonal, so the result of adding and multiplying them is also block diagonal. For a block diagonal matrix with k blocks of size d × d, the time complexity for a Cholesky decomposition isO(kd 3 ) = O(M d 2 ). Combined with the complexity of computing z, the recovery step takes O(M d 2 + N Ld + N d 2 ) time. 
marginalizing multiple effects with additional assumptions: 
 We have shown that it is efficient to marginalize one class of random effects. With additional practical assumptions, it is possible to marginalize all classes of random effects for efficient HMC inference. Instead of separating different classes of random effects, LMMs can also be written asv ∼ N (µ, Σ v ), y ∼ N (Bv + b, Σ y ), where B = [A 1 , ..., A L ] and v = [u T 1 , ..., u T L ] T . We define that D = L i=1 M i .The matrix inversion and determinant lemmas can still be applied to marginalize v out, but the combined matrix B does not have the special structure of A i we exploited in Section 3. More specifically, the computation of det(F) and the evaluation ofF -1 for F = Σ -1 v + B T Σ -1y B both become non-trivial. We introduce additional assumptions to show that they can be solved faster in some special cases. For the general case, see the discussion section.The assumption we make is that Σ v = τ v I and Σ y = τ y I, where τ v , τ y are scalars that either belong to Θ or are fixed non-random parameters. This means that all effects share the same variance and all observations share the same noise scale. These assumptions are not as restrictive as it may appear. If the underlying distribution is u i ∼ N (µ i , σ 2 i I) where σ i is a fixed parameter, it is possible to reparameterize this distribution asu ′ i ∼ N (0, I), A ′ i = σ i A i , b ′ = b + Bµ i , and use u ′ i , A ′ i , b ′ in place of u i , A i , b.Then Σ v becomes a scaled identity matrix. Also, in many models, the noise scale for different observations is the same, making Σ y a scaled identity matrix as well.In practice, if the assumptions are satisfied, marginalization can be done in O(D 2 + N d) time with O(D 3 + N L 2 d 2 ) preprocessing. Details are provided in Appendix B.3. 
related work: 
 While many works aim to improve HMC directly [71, 30, 58, 73] , a number of other works focus on model transformation. Non-centered parameterization [49] is a widely used trick among MCMC users to alleviate slow sampling in difficult posterior distributions. However, there is no general way to know whether a non-centered parameterization will be beneficial [76] . Variationally inferred parameterization [26] proposes to learn a model parameterization from a specified family that will lead to effective sampling. In Parno and Marzouk [52] and Hoffman et al. [33] , preconditioners for HMC are learned to transform the model to be approximately isotropic Gaussians. Marginalization differs from reparameterization in that it reduces the problem dimension as well as potentially alleviating difficult characteristics such as funnels, so it has two mechanisms to improve MCMC efficiency. The Laplace approximation (LA) is one way to approximately marginalize variables in MCMC [59, 40, 65] , but it may be difficult to quantify the error or recover the marginalized variables.Marginalization, or Rao-Blackwellization, has been an important topic in Bayesian inference and probabilistic programming. In Gibbs sampling, marginalization is usually called collapsing [37] . Collapsed Gibbs sampling has been developed for latent Dirichlet allocation [56] and LMMs [50] . We explore marginalization in the context of HMC, which induces different considerations. Methods with HMC do not have to make the conditional distributions of the marginalized model tractable. Marginalization is also related to symbolic inference in probabilistic programming. Hakaru [44] and PSI [21, 22] are systems for performing exact Bayesian inference by symbolically marginalizing all latent variables. To marginalize discrete variables, Gorinova et al. [27] propose an information flow type system. Another line of related work is delayed sampling [43, 3] , which automates marginalization of variables within Rao-Blackwellized particle filters [42] . Lai et al. [35] developed an automatic system for marginalizing variables in HMC, but is limited to scalar variables so cannot leverage vectorization and forces users to write models with univariate distributions. Linear algebra tricks have been widely utilized in various machine learning algorithms, such as ridge regression [68] , Gaussian processes [61] and Kalman filters [60] . Recently, frameworks [62, 20, 57] have been proposed to ease the implementation of fast linear algebras in machine learning algorithms. Marginalization in Bayesian models may be an interesting application of those frameworks.Fast and scalable inference for LMMs has been studied in the context of maximum likelihood estimation [19] , variational EM [24] , Gibbs sampling [51] and numerical integration [28] . We are the first to consider speeding up the inference of LMMs with HMC. There is also a recent trend in integrating random effects into deep neural networks for correlated data [67] or personalization [66, 64, 74] with parameters estimated by maximum likelihood.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: We provide the data and code in the supplemental material. Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experiments: 
 We conduct experiments on LMMs from various disciplines using the default no-U-turn sampler (NUTS) [34] from NumPyro [5, 54] , which has an adaptive step size with dual averaging, adaptive and diagonal mass matrix, target acceptance probability of 0.8, and maximum tree depth of 10. For the ETH instructor evaluation model, we set the maximum tree depth to 12 to overcome difficulties performing inference without marginalization in preliminary experiments. For all models, we use weakly informative priors unless specified. In general, our conclusion is insensitive to the choice of hyperparameters and priors. For all experiments, we collect 10,000 warm up samples for tuning, and 100,000 samples for evaluation, and evaluate performance via effective sample size (ESS) and running time. Cross-effects models are a type of LMM that have more than one class of random effects (i.e. L > 1). Usually each observation belongs to one subject group (e.g. individuals, animals) and one item group (e.g. questions, objects). The correlation among latent variables can create severely challenging geometry that slows down the sampling of HMC. With our idea, it is possible to marginalize one or more group of effects from the model, reducing the dimension of latent space for faster sampling and better geometry. ETH instructor evaluations An example cross-effects model describes university lecture evaluations by students at ETH [4] . The dataset records N = 73421 ratings, where each rating y n comes from student s n for professor p n teaching a course from department d n , with t n indicating whether the professor is teaching outside their own department. There are a total of M 1 = 2972 students, M 2 = 1128 professors and M 3 = 14 departments. We use a version of the model from the document of Tensorflow probability [13] . The model is Likelihood : y n ∼ N (u 1,sn + u 2,pn + u 3,dn + α + βt n , σ 2 ), Prior : u 1,i ∼ N (0, 1), u 2,j ∼ N (0, 1), u 3,k ∼ N (0, 1), α ∼ N (0, 5), β ∼ N (0, 1), σ ∼ N + (0, 1), where1 ≤ i ≤ M 1 , 1 ≤ j ≤ M 2 and 1 ≤ k ≤ M 3 .Given the dataset, we wish to learn about the latent variables u 1 , u 2 , u 3 , α, β and σ. HMC is the most direct way to sample those variables, but the dimension and complicated relations make it inefficient. Marginalization can be applied to one of the effects, u 1 , u 2 or u 3 . We report the running time of sampling from the model with and without marginalization in Table 2 . We found that marginalizing any group of random effects improves the sampling speed of HMC. However, the improvements are not necessarily predicted by the dimension of marginalized variable: HMC is faster when marginalizing u 3 than when marginalizing u 1 even though u 1 has 200-times higher dimension than u 3 . In Figure 2 , the ESS for each variable is reported. Without marginalization, sampling u 2 and u 3 are both difficult compared to sampling u 1 , and HMC becomes more efficient when marginalizing either of these variables, so we conjecture that u 2 and u 3 are responsible for the difficulty for sampling in the original model. In this model, all random effects are independent and have the same variance, so Σ u is a scaled identity matrix and we can marginalize all random effects efficiently. This approach is observed to be the most efficient in our experiments, despite having quadratic complexity in D. Overall, marginalization never hurts ESS, and runs faster. We expect that any marginalization strategy works better than HMC in the original model, a finding which will be consistent across experiments. Additional results of this experiment, including trace plots and R diagnosis, are included in Figure 6 and Table 5 in the Appendix.To tackle bad geometry in statistical models, another model transformation is non-centered parameterization, or reparameterization [49] . Reparameterization converts the distribution of z ∼ N (µ, σ 2 ) into ϵ ∼ N (0, 1) and z = ϵσ + µ. Reparameterization is especially useful for funnel shapes in hierarchical models. We note that when applicable, marginalization is able to solve a broader class of problems. We compare marginalization and reparameteriation on the grouse ticks model. The dataset [4] contains observations y of the the number of ticks on the heads of red grouse chicks in the field. Each observation y k comes from brood b k in location l k during year e k at altitude a k , where year and altitude give fixed effects, and there are random effects u 1 and u 2 corresponding to brood and location. There are N = 403 observations, M 1 = 118 broods and M 2 = 63 locations. We define the hierarchical model as follows:Likelihood :EQUATIONt ) Prior : µ 1 ∼ N (0, 1), σ 1 ∼ HalfCauchy(5), µ 2 ∼ N (0, 1), σ 2 ∼ HalfCauchy (5) ,EQUATION2 ), σ t ∼ HalfCauchy (5) , where i = 1, ..., M 1 , j = 1, ..., M 2 , k = 1, ..., N and each y k is observed. The correlation between σ and u creates the funnel shape that makes vanilla HMC inefficient. Nevertheless, it is possible to apply either marginalization or reparameterization to each random effect. In Figure 3 , we plot the distributions of samples for variable pairs (σ 1 , u 1,1 ) and (σ 2 , u 2,1 ) with different combinations of marginalization and reparameterization. There is a difficult correlation between σ 2 and u 2 . After applying marginalization or reparameterization to u 2 , HMC manages to explore the funnel region (at low values of σ 1 ). However, we find that only samplers that marginalize u 2 report zero divergent transitions after warm-up. Such behavior is consistent with different random seeds. See Table 6 in the Appendix. Also, the distribution of divergent samples is related to specific parameters when reparameterizing u 2 , implying that reparameterization introduces pathologies that create challenges for HMC inference. In addition, we find that reparameterization does not improve the running time of HMC, while marginalizing u 2 speeds up sampling by about 20%.In theory, marginalization with LMMs can be done by constructing a graphical model for scalar random variables and performing automatic marginalization as in [35] . But it is more efficient to marginalize in a vectorized way. We demonstrate the benefits from vectorization in Table 3 . Both marginalization strategies are performed on two hierarchical linear regression models, the electric company model [23] and the pulmonary fibrosis model [63] . We find that vectorized marginalization is much more efficient for sampling from the two models. dillonE1 [12] dutch [17] eeg [48] english [69] gg05 [29] mandarin [75] mandarin2 [70] pupil [72] stroop [ Figure 4 : Experimental results for the 9 cognitive science datasets with and without marginalization. Each experiment is performed 5 times with different random seeds. Marginalization usually improves sampling speed measured by iterations per second (iter/s) and sample efficiency measured by ESS per iteration (ESS/iter).Hierarchical Bayesian inference with LMMs has wide applications in cognitive science [47] . We highlight the effectiveness of marginalization with 9 datasets from cognitive science (Table 4 ). They cover various settings, with one or two random effects, normal or log-normal likelihoods, on CPU or GPU. Experiments that are slow on CPU are performed on GPU. Each dataset corresponds to an LMM where both the intercept and the coefficient include random effects. Details of all the models can be found in Appendix D. Results are summarized in Figure 4 . Marginalization usually improves the sampling speed of HMC and consistently improves efficiency measured by ESS per iteration.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results? Answer: [Yes] Justification: The hyperparameters for all models are included in the appendix. Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.The code is available at https://github.com/lll6924/hamiltonian_lme.git.38th Conference on Neural Information Processing Systems (NeurIPS2024).Each A T i Aj is O(N d 2 ), as a corollary of Theorem 1. 
discussion: 
 There are several promising directions for future work.Marginalization is related to Rao-Blackwellization. This paper focuses on marginalization, which improves the speed of obtaining samples from the remaining variables by improving mixing times, reducing the cost per iteration, or both. Combining marginalization with Rao-Blackwellization is an interesting avenue for future work. More formally, if one is interested in some expectation E (Θ,u)∼p(Θ,u|y) [f (Θ, u)] in an LMM, there is a Monte Carlo estimatorE 1 = 1 N N i=1 f (Θ i , u i ),where (Θ i , u i ) ∼ p(Θ, u|y) and N is the sample size. Marginalization is a trick to improve the efficiency of the posterior sampling, so that we can achieve the same estimation variance with smaller N or less runtime . At the same time, we also have access to a conditional distribution that is useful for Rao-Blackwellization. If the effects variable u can be marginalized we have both an approximate posterior for p(Θ|y) and an analytical conditional distribution p(u|Θ, y). With Rao-Blackwellization we have thatE (Θ,u)∼p(Θ,u|y) [f (Θ, u)] = E Θ∼p(Θ|y) [E u∼p(u|Θ,y) [f (Θ, u)]].In such case, another Monte Carlo estimator can be constructed:E 2 = 1 N N i=1 E u∼p(u|Θ,y) f (Θ i , u) ,where Θ i ∼ p(Θ|y). For some functions, such as those that are polynomial in u, the inner expectation can be computed exactly using properties of Gaussians. In other cases, the inner expectation can be estimated cheaply via Monte Carlo using exact samples from p(u|Θ i , y).In Section 4, we proposed to marginalize multiple classes of random effects by assuming a scaled identity covariance matrix. To marginalize multiple effects in general models, a possibility is to compute z T E -1 z and estimate det(E) and the corresponding gradients with conjugate gradient (CG) solvers [14, 20] . However, this approach uses stochastic estimators for the determinant and gradients, which introduce bias into the HMC dynamics. These biases can be corrected through pseudo-marginalization [2] , but it is unclear how significantly the extra stochasticity will affect the sampling. Another possible way to marginalize multiple effects for LMMs is to introduced the balanced levels assumption [50] . We leave these ideas for future exploration.In this work, we only consider normal or log-normal likelihoods, but our method can be easily generalized to other deterministic transformation of normal likelihood. This implies that marginalization can benefit regression with most continuous predictors given proper link functions. Another potential future direction is to marginalize classification models with probit regressions [1] . Marginalization will turn probit models into multivariate probit models as AΣ u A T + Σ y is a dense covariance matrix, which may require a simulation-based method [10] or variational Bayes [39] . It will be interesting to see how ideas from multivariate probit regression could be fit into an HMC pipeline. In a broader context, marginalization is related to data augmentation techniques that "create" conjugacy for non-normal likelihoods or non-normal effects. Those techniques were developed for Gibbs sampling, e.g. [18, 55] , but may also be useful for HMC.We have developed a tool to speed up the HMC inference for LMMs. In our implementation, the marginalized likelihood p(y|Θ, u -i ) is defined as a special type of parametric distribution available to the user, and the recovery distribution p(u i |Θ, u -i , y) is a function called after sampling. In our experiments, marginalization never hurt sampling efficiency measured by ESS/s, and usually helped. Thus, it would be desirable to always marginalize one group of random effects when the model is an LMM. Future work could aim to automatically apply such transformations to user-specified LMMs.There are two possible high-level approaches. The first is to perform marginalization starting with a model described using a high-level abstraction such as an R formula. Then, when compiling the high-level model description into a concrete model (e.g., a probabilistic program), we can marginalize one or more of the effects using our methods. The second is to perform marginalization starting with a user-written probabilistic program representing an LMM. In this case, some compilation or program tracing technique will be needed to convert the user's program to a model representation suitable for manipulation. For example, Lai et al. [35] used program tracing to construct a graphical model representation that could be programmatically analyzed and transformed. To apply this methodology to LMMs, a special parser would also be needed to match the models to LMMs.We summarize the important symbols used in the paper. We first review the tree structure of the matrix A. A is an N × M matrix where every block of d columns corresponds to the effects for one group (e.g., an individual subject, age, school, or gender). For example, if N = 3, k = 2 and d = 2, one possible graphical model is as below. Each u j ∈ R 2 . If the coefficients are all 1s, thenu 1 u 2 y 1 y 2 y 3A = 1 1 0 0 1 1 0 0 0 0 1 1 .To generalize, if for y i , the grouping variable is g i , then in the ith row of A, only A i,j:k can be nonzero for j = (g i -1)d + 1 and k = g i d. We consider three representations of the matrix A. By rows,A =    r 1 r 2 ... r N    , by columns, A = (c 1 c 2 ... c kd ) ,and by blocks of d columns,A = A :,1:d A :,d+1:2d ... A :,(k-1)d+1:kd = (C 1 C 2 ... C k )where each C i (i = 1, 2, ..., k) is N × d. Now we restate and prove Theorem 1.Theorem 1. If Σ y is diagonal, Σ u is block-diagonal with blocks of size d × d, then F = Σ -1 u + A T Σ -1 y A is also block-diagonal with d × d blocks and computing A T Σ -1 y A takes O(N d 2 ) time.Proof. The theorem has two parts: (a) the property of F, and (b) the computation of F. We address them with the three representations of A.(a) F is block-diagonal. Because Σ u is block-diagonal, Σ -1u is also block-diagonal with the same sizes. Also, Σ y is diagonal, so the block-diagonality of A T Σ -1 y A is the same as A T A. We consider the column representation of A, thenA T A =    C T 1 C T 2 ... C T k    (C 1 C 2 ... C k ) =    C T 1 C 1 C T 1 C 2 ... C T 1 C k C T 2 C 1 C T 2 C 2 ... C T 2 C k ... C T k C 1 C T k C 2 ... C T k C k    . For 1 ≤ i ≤ k, C T i C i is d × d. For 1 ≤ i < j ≤ k, C T i C j =     c T (i-1)d+1 c T (i-1)d+2 ... c T id     c (j-1)d+1 c (j-1)d+2 ... c jd .The following lemma shows that C T i C j = 0. Lemma 1. For any 1 ≤ i < j ≤ k and 1 ≤ s, t ≤ d, it holds that c T (i-1)d+s c (j-1)d+t = 0.Proof. The lemma can be proved by contradiction. Suppose c T (i-1)d+s c (j-1)d+t ̸ = 0. Then there exists an index n such that c (i-1)d+s [n] ̸ = 0 and c (j-1)d+t [n] ̸ = 0. This means that in the nth row of A, both A n,(i-1)d+s and A n,(j-1)d+t are non-zero. This contradicts with the tree-structure where only one group of d elements can be non-zero in a row.A T A is block-diagonal, thus Σ -1 u + A T Σ -1 y A is also block-diagonal and each block is d × d. (b) The computation of A T Σ -1 y A is O(N d 2 ). Since Σ y is diagonal, A ′ = Σ -1y A has the same pattern of zeros and nonzeros as A. We consider the row representations such thatA ′ =    r ′ 1 r ′ 2 ... r ′ N    . Then A T A ′ = r T 1 r T 2 ... r N    r ′ 1 r ′ 2 ... r ′ N    = N i=1 r T i r ′ i .note that each of r i and r ′ i has d non-zero elements. So computingA T A ′ is O(N d 2 ).B.2 Pseudocode for recovery after marginalizing one group of random effectsAlgorithm 2 Sampling from p(u i |Θ, y, u -i ) z = y -j̸ =i A j u j -A i µ i -b ▷ Sparse matrix multiplication in O(N Ld) time G = A T Σ -1 y A ▷ Block diagonal computation in O(N d 2 ) time F = Σ -1 u + G ▷ Block diagonal computation in O(M d 2 ) time µ = µ i + Σ u (I -GF -1 )A T Σ -1 y z ▷ Sparse matrix multiplication in O((M + N )d) time L = Σ u Cholesky(Σ -1 u -G + GF -1 G) ▷ Cholesky of block diagonal matrix in O(M d 2 ) time return u i ∼ Normal(µ, LL T ).With the assumptions of scaled identity covariance matrices, all effects can be marginalized with a preprocessing of the eigendecomposition of B T B.(a) Preprocessing before HMC. We computeB T B = QΛQ T .In LMMs, the computation of B T B is O(N L 2 d 2 ) 2 , and the eigendecomposition of it is O(D 3 ). So the overall complexity for preprocessing is O(D 3 + N L 2 d 2 ). Compared with the HMC sampling loop that takes thousands of steps and visits the model hundreds of times each step, the cost of preprocessing is not expensive. In our attempt to marginalize all random effects for the instructor evaluation model in Section 6.1, this step takes less than 10 seconds.(b) Marginalized likelihood during HMC. During HMC sampling, the log density log p(y|Θ) would be calculated, which islog p(y|Θ) = - 1 2 det (E) - 1 2 z T E -1 z + Cwhere z = y -Bµb and E = BΣ v B T + Σ y . The computation of z takes O(N Ld) time. With the two lemmas, we havedet(E) = det(Σ -1 v + B T Σ -1 y B) det(Σ v ) det(Σ y ), z T E -1 z = z T Σ -1 y z -z T Σ -1 y B(Σ -1 v + B T Σ -1 y B) -1 B T Σ -1 y z.A shared matrix in the formulas isF = Σ -1 v + B T Σ -1 y B.ThenF = Σ -1 v + B T Σ -1 y B = Q 1 τ v I + 1 τ y Λ Q T .With the trick, evaluating det(F) reduced to O(D) time as det(Q) = 1. Also z T E -1 z becomesz T E -1 z = z T Σ -1 y z -z T Σ -1 y BQ 1 τ v I + 1 τ y Λ -1 Q T B T Σ -1 y z.Note that B T Σ -1 y z ∈ R D can be computed in O(N Ld) time, but its multiplication with Q T takes O(D 2 ) time. Given that 1 τv I + 1 τy Λ is diagonal, the complexity of evaluating log p(y|Θ) once is then O(D 2 + N Ld).(c) Ancestral sampling after HMC. In the recovery step, we perform ancestral sampling from p(v|Θ, y). To efficiently generate samples, we give the following theorem. Theorem 2. If Σ v = τ v I, Σ y = τ y I and B T B = QΛQ T , then v|Θ, y ∼ N (µ v|Θ,y , Σ v|Θ,y ), whereµ v|Θ,y = µ + τ v τ y B T - 1 τ y QΛ 1 τ v + Λ τ y -1 Q T B T z, Σ v|Θ,y = Q τ v I - τ v 2 τ y Λ + τ v 2 τ y 2 1 τ v + Λ τ y -1 Λ 2 Q T .In Theorem 2, from z, we can apply matrix multiplications from right to left to get µ v|Θ,y . The whole computation takes O(D 2 + N Ld). To generate normal samples a Cholseky factorization forΣ v|Θ,y is required. But τ v I -τv 2 τy Λ + τv 2 τy 2 1 τv + Λ τy -1Λ 2 is diagonal, so it can be obtained in O(D 2 ) time as well. Now we prove Theorem 2.Proof. µ v|Θ,y and Σ v|Θ,y can both be derived algebraically. • The answer NA means that the paper does not include experiments.µ v|Θ,y = µ + Mz = µ + Σ v B T (BΣ v B T + Σ y ) -1 z = µ + Σ v B T (Σ -1 y -Σ -1 y B(Σ -1 v + B T Σ -1 y B) -1 B T Σ -1 y )z = µ + τ v τ y (B T - 1 τ y B T B(Σ -1 v + B T Σ -1 y B) -1 B T )z = µ + τ v τ y (B T - 1 τ y QΛQ T Q 1 τ v I + 1 τ y Λ -1 Q T B T )z = µ + τ v τ y (B T - 1 τ y QΛ 1 τ v I + 1 τ y Λ -1 Q T B T )z. Σ v|Θ,y = (I -MB)Σ v = (I -Σ v B T (BΣ v B T + Σ y ) -1 B)Σ v = τ v I -τ v 2 B T (BΣ v B T + Σ y ) -1 B = τ v I -τ v 2 B T (Σ -1 y -Σ -1 y B(Σ -1 v + B T Σ -1 y B) -1 B T Σ -1 y )B = τ v I -τ v 2 B T Σ -1 y B + τ v 2 B T Σ -1 y B(Σ -1 v + B T Σ -1 y B) -1 B T Σ -1 y B = τ v I - τ v 2 τ y B T B + τ v 2 τ y 2 B T B(Σ -1 v + B T Σ -1 y B) -1 B T B = τ v I - τ v 2 τ y QΛQ T + τ v 2 τ y 2 QΛQ T Q 1 τ v I + 1 τ y Λ -1 Q T QΛQ T = τ v I - τ v 2 τ y QΛQ T + τ v 2 τ y 2 QΛ 1 τ v I + 1 τ y Λ -1 ΛQ T = Q τ v I - τ v 2 τ y Λ + τ v 2 τ y 2 1 τ v + Λ τ y -1 Λ 2 Q T .• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 
 
 Appendix:  


Decision: This manuscript considers Bayesian inference in a popular class of models, linear mixed-effects models. The authors discuss strategies to improve the performance of Hamiltonian Monte Carlo sampling in this class of models, including how to use ideas from numerical linear algebra to accelerate the running time dramatically.

During the initial review period, the reviewers generally praised the work, identifying only some minor concerns, which were adequately addressed during the author-reviewer discussion period. Following the discussion period, there was consensus in favor of accepting the work.

Reviews: review: paper_summary: paper_summary: This paper considers the problem of Bayesian inference in linear mixed-effects models. This is a popular class of models used throughout many areas of application. Inference in such models depends on being able to determine values of latent variables (fixed and observed effects) given the observed data. This is typically done by writing a model of interest in a probabilistic programming language and then running an automated inference algorithm on it. For example, the Stan package runs a variant of Hamiltonian Monte Carlo on the model. In general, this can be a complex problem from a sampling point of view. The main contribution of this paper is to speed up vectorized random effects marginalization, which simplifies the sampling problem and leads to faster sampling. The authors note that marginalization tends to always be beneficial, and illustrate their proposed methodology on a well-chosen set of real-life examples.
strength_weakness: strength_weakness: The proposed methods are limited to Gaussian models, which the authors note. I would have appreciated some more commentary in this direction. The marginalization approach seems sufficiently interesting to be used widely, whenever possible.
questions: questions: What sort of prior structures allow for this to be applied?

Are there methods that can be developed that perform approximate marginalization in some way?

How hard is it to automate marginalization in practice?
limitations: limitations: Yes.
paper_summary: paper_summary: The paper studies marginalization in linear mixed-effects models with a Gaussian likelihood.
Marginalizing out variables is a common trick to improve the convergence of MCMC.
As pointed out, for the linear mixed-effects models with a Gaussian likelihood, naively doing so will actually increase the computation as there are a few linear algebra operations that is costly.
The paper introduces two tricks using the matrix inversion lemma to speed them up.
Some experiments are performed to demonstrate the effectiveness and to compare against the reparameterization alternatives.
strength_weakness: strength_weakness: ## Models studied are narrowed
Only models with Gaussian or log-Gaussian likelihood is studied.

## Experiments are not enough
1. In the experiments there is only naive baselines. Other alternative automatic marginalization methods are discussed in the related work section but not compared against.
2. Only 1 dataset (section 6.1) and 1 toy target (section 6.2) are studied.
questions: questions: Is the method specific to HMC?
limitations: limitations: Limitations are discussed in section 7 however I do think those limitations make the paper a bit weak (not enough contributions) for acceptance.
paper_summary: paper_summary: The paper addresses the challenge of marginalizing latent variables in linear mixed-effects models (LMMs) to improve Hamiltonian Monte Carlo (HMC) sampling. The authors propose an algorithm for automatic marginalization of random effects in LMMs, leveraging fast linear algebra techniques. The experiments evaluate the performance of the proposed marginalisation on LMMs from several disciplines by using the No-U-Turn sampler in NumPyro and assessing performance via effective sample size (ESS), number of divergent transitions, and running time.
strength_weakness: strength_weakness: My two main concerns about the paper are regarding the experiments and about the significance of the results.

### Experimental Details: 

Overall, I find that the experiments lack sufficient details and metrics:
* There is not enough information about hyperparameters used for each method and model, e.g. the step size for NUTS. Does the approach rely solely on NumPyro's adaptation routine? If yes, what parameters were used for that?
* The evaluation can use additional diagnostics and metrics to provide more evidence that marginalization improves inference: for example trace plots, and visualisation of divergences can be very useful. Effective sample size alone is not enough for properly assessing the quality of inference.
* I don’t think the number of divergences is a reliable metric, without additional context, such as their location in the posterior distribution or the step size used. In the case with section 6.2, the divergences are not neccesarily indicative of any problems. In my opinion, the 12 divergences observed for R1, R2 are just 0.012% of all drawn samples, which is negligible. Where are those divergences located in the posterior plot? Are they grouped in one place? 
* The choice of priors for various models seems to me quite broad and lacks justification. 

### Significance of Results

The claim that marginalization is always beneficial is not novel, as it follows from the Rao-Blackwell theorem. The paper mentions Rao-Blackwellization, but doesn’t go into details.

Edit: the rebuttal clarified on all of these points.
questions: questions: 1. Can you provide more detailed convergence diagnostics for your experiments, such as trace plots, for example?
2. Are divergences grouped in specific regions of the posterior distribution?
3. What is the impact of varying HMC parameters (e.g., step size, or adaptation routine parameters) on the number of divergences and the effective sample size?
4. Could you elaborate on the choice of priors in your experiments? Why were these particular priors selected?

Edit: the rebuttal has answered all my questions.
limitations: limitations: The last section of the paper addresses a few limitations of the approach, including addressing distributions beyon normal. However, the paper can benefit of further discussion on one, in my opinion, fundament limitation of this approch, which is that it only considers cases where analytical marginalization is possible. How many practical cases are covered by this constraint?
paper_summary: paper_summary: This work considers Bayesian inference via MCMC methods (specifically: Hamiltonian Monte Carlo) in linear mixed-effects models (LMMs). Due to linearity and Gaussianity, the random effects can be analytically integrated out so that the MCMC algorithm needs to only target a distribution on a smaller space which does not include the random effects. This typically leads to more efficient MCMC algorithms. However, as the authors show, a naive marginalisation is costly, incurring a complexity that is cubed in the number of observations and the size of the vector of random effects. The main contribution is to show that, by exploiting sparsity, the random effects can be integrated out in a way that incurs (typically) the same computational complexity as evaluating the likelihood without any marginalization. The authors demonstrate that marginalization in this way improves the performance of the HMC algorithm compared with "no marginalization" both in terms of the effective sample size and effective sample size per second.
strength_weakness: strength_weakness: **Title**

I don't think the title (and the emphasis on HMC in the abstract and introduction) is entirely adequate. The wording "marginalized Hamiltonian Monte Carlo" implies to me that some component specific to HMC (e.g. the momentum variable) is analytically integrated out. However, the marginalization here is completely orthogonal to HMC. Indeed, it could be exploited within any other MCMC algorithm or even within other Monte Carlo schemes like sequential Monte Carlo samplers). 


**Presentation**

Overall, I think the presentation is good. Just a few minor comments:

1. In Table 1 and Figures 2, 3 and 4, referring to the "no marginalization" strategy as "HMC" is confusing since all methods use HMC. This is also inconsistent with, e.g., Table 2. Better replace "HMC" by "no marginalization".
2. A number of symbols are not or not clearly defined, e.g., on Page 2, define $\mathbf{u}_j$, $\mathbf{u}$ (is only defined much later), and $k$ (should $k$ be $M$?). More generally, the fact that symbols without the subscript $i$ denote the collection over all $i$ should be mentioned (unless I missed it).
3. From the context, it is clear what the authors mean by "recovery" and "ancestral sampling" but unless these terms are somehow standard in the literature, I think they should be more clearly explained to avoid confusion.
4. Perhaps the distinction between "classifications" and "groups" (around Lines 88 to 92) could be made more clear, especially because the examples in parentheses say "subject", "gender", and "age" in both cases which is a bit confusing.
5. I couldn't quite understand Section 6.3. What does "marginalization in a vectorized way" mean?
6. Check the bibliography; it has many mistakes, e.g., missing capital letters in proper nouns and journal names.


**Other**
* L241: document -> documentation
questions: questions: 1. Can this be extended to multivariate observations?
2. As the authors mention, the "no marginalization" implementation seems to be slower than the proposed implementation that performs (some) marginalization. This seems surprising. Do the authors have an explanation for this?
limitations: limitations: I don't foresee any negative societal impact.

Comments: title: Rebuttal Response
comment: Thanks for the clarification.
Two of my main concerns are resolved and I've raised my score accordingly.
Dear authors, thank you for your response. I appreciate the extensive rebuttal and the additional results you presented in the attached pdf. Thank you, in particular, for the plots visualising where divergences happen - I think these really help highlight that there is indeed a problem with the geometry that your method seems to address. 

This has addressed all major concerns I had and I will increase my score accordingly. Looking forward to seeing the updated version of the paper when available.
Thank you for the clarifications!
Thank you for these clarifications.