{'fraction_non_alphanumeric': 0.026813085289547726, 'fraction_numerical': 0.00770587494225843, 'mean_word_length': 4.753442860594347, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Parallelized integrated nested Laplace approximations for fast Bayesian inference

There is a growing demand for performing larger-scale Bayesian inference tasks, arising from greater data availability and higher-dimensional model parameter spaces. In this work we present parallelization strategies for the methodology of integrated nested Laplace approximations (INLA), a popular framework for performing approximate Bayesian inference on the class of Latent Gaussian models. Our approach makes use of nested OpenMP parallelism, a parallel line search procedure using robust regression in INLA's optimization phase and the state-of-the-art sparse linear solver PARDISO. We leverage mutually independent function evaluations in the algorithm as well as advanced sparse linear algebra techniques. This way we can flexibly utilize the power of today's multi-core architectures. We demonstrate the performance of our new parallelization scheme on a number of different real-world applications. The introduction of parallelism leads to speedups of a factor 10 and more for all larger models. Our work is already integrated in the current version of the open-source R-INLA package, making its improved performance conveniently available to all users.


Introduction
The methodology of integrated nested Laplace approximations (INLA) has become a widely spread framework for performing complex Bayesian inference tasks, see e.g. [4,38,46,46]. INLA is applicable to a wide subclass of additive Bayesian hierarchical models. It employs a deterministic approximation scheme that relies on Gaussian Markov random fileds (GMRFs) in the latent parameter space which enables the usage of efficient sparse linear algebra techniques [45]. A user-friendly implementation of INLA is available in the form of an R-package under the same name, referred to as R-INLA. It can be downloaded and installed as described in www.r-inla.org with the complete source code available on GitHub 1 . Since its inception there have been many papers exploring advancements of theoretical concepts of the INLA methodology, also leading to a constant evolution of its implementation. Their agglomeration forms an impressive repertoire for fast, versatile and reliable approximate Bayesian inference, see e.g. [31,34], including a wide variety of applications, see e.g. [1,5,7,10,13,21,23,24,28,29,32,33,36,40,41,49,50]. Algorithmic concepts and improvements concerning performance matters of the continuously growing software library have been discussed much less, even though they are a key component to INLA's success. In this work we want to invite the reader to look at the INLA methodology with us from a more algorithmic point of view. We will present a much more performant implementation of R-INLA making use of OpenMP [15] parallelism and the state-of-the-art sparse linear solver PARDISO [9].
The introduction of OpenMP allows for the simultaneous execution of multiple tasks within the algorithm operating on shared memory. We have parallelized the computationally intensive operations within INLA whenever the dependency structure allowed it, which also opened the door for us to add more robust approximation strategies. Internally, the PARDISO library is also utilizing OpenMP which leads to a double layer of parallelism within INLA. A large quantity of the algorithm's overall runtime is spent in an optimization routine to determine the posterior mode of the model's hyperparameters. The mode is found iteratively using a quasi-Newton method that requires gradient information of the functional being maximized. Approximating the gradient in each iteration is a computationally expensive task, which we were, however, able to parallelize. Additionally, INLA performs a line search procedure in every iteration of the optimization to find a suitable parameter configuration for the next step. Instead of sequentially looking for a good candidate we have parallelized this process, again, using OpenMP while also including a robust regression scheme for more stable results. This can allow us to save time in each iteration and can also reduce the total number of iterations required. To compute the posterior marginal distributions of the latent parameters of the model partial matrix inversions are necessary. For all computations INLA generally works with the sparse precision matrices instead of their dense counterparts, the corresponding covariance matrices. However, to obtain the marginal variances, inversions of often high-dimensional precision matrices are necessary. A time-consuming task which we have adapted to be executed in parallel through the usage of OpenMP. The collection of these strategies forms the first or top layer parallelism.
The original implementation of INLA employs the sparse linear algebra library TAUCS [54] to perform the required numerically intensive core operations. Foremost, these include Cholesky decompositions of matrices with recurring sparsity patterns and the partial matrix inversions. TAUCS is a well-designed and efficient library which operates, however, sequentially and whose support was discontinued in the early 2000's. Additionally, TAUCS does not include a partial inversion routine which was therefore implemented by the authors of R-INLA, see [44,45], in a sequential manner. While this was done with much thought, effort and consideration, providing reliable results, improved parallel implementations have emerged since [27,58]. They were developed by experts in the field of high-performance computing and sparse numerical solvers. The plan of integrating PARDISO into INLA was first described by Nierkerk et al. in [56]. The authors discuss the need for the usage of faster numerical solvers within INLA to support evolving statistical models of higher complexity and continuously growing availability of data. They also show the efficiency and scalability with which PARDISO performs the beforementioned numerically intensive core operations like Cholesky decompositions and partial matrix inversions. Especially for the latter task there are not many performant libraries available which make PARDISO a particularly great fit [56]. Hence, the second layer of parallelism is formed by PARDISO as well as other parallelized linear algebra operations like matrixmatrix or matrix-vector products. We will present performance results for three different real-world applications, quantifying the improvement compared to previous implementations in terms of speedup and scalability through the various parallel OpenMP strategies. The case studies include a complex joint survival model containing 50 hyperparameters [47], a brain activation model using fMRI data with very high-dimensional spatial domains [52] and a smaller scale model describing spatial variation in Leukemia survival data [31]. All corresponding code and data is publicly available.
The rest of the paper is organized as follows. We will begin with a brief introduction of the INLA methodology, presenting the class of applicable models and underlying statistical concepts, with a particular focus on the computations they entail. Additionally we will provide an overview of the fundamental numerical methods and sparse linear algebra operations involved in the implementation. Then we will introduce the new parallelization schemes, before demonstrating our improvements on the various applications.

Latent Gaussian Models
The INLA methodology is applicable to the class of latent Gaussian models (LGMs). They comprise a subclass of hierarchical Bayesian additive models among which there are many of the frequently used statistical models, like regression, mixed or spatio-temporal models as well as many others, see e.g. [4,34,45,46] for details. Each observation y i is assumed to belong to a distribution from the exponential family and is associated with the additive linear predictor η i through a link function, where η i is defined as The linear predictor has fixed effects β with covariates Z i and nonlinear random effects u i with covariates w i . The observations are assumed to be conditionally independent given the parameters, such that where θ are the hyperparamters of the model. The latent parameter space of the model has traditionally been defined as x = (β 0 , β, u, η), but can now also be formulated without the linear predictor as x = (β 0 , β, u), which was recently presented in [57]. The concepts discussed in this work are equally applicable to the two formulations. In both cases we assume that x forms a Gaussian Markov random field (GMRF) with zero mean and sparse precision matrix Q(θ), i.e. x ∼ N (0, Q(θ)).
One of the key components to INLA's computational efficiency is this sparse GMRF structure.
It allows for a very high-dimensional parameter space without performance degradation due to the employed sparse linear algebra operations [43]. In addition one adopts a prior π(θ) for the hyperparameters θ. Thus, forming a three-stage model consisting of the hyperparamter distribution, the latent field and the likelihood. It is also possible account for constraints imposed on the latent parameter space of the model [43]. For simplicity, we will not discuss this case in further detail, as it is similar to the unconstrained framework when it comes to parallelization strategies. For a further discussion of the model specifications or Bayesian hierarchical models see [11,46].
(3) The term π G (x|θ, y) describes a Gaussian approximation of π(x|θ, y) evaluated at the mode of x for a given θ. Through an iterative optimization algorithm the hyperparameter configuration θ * for which θ * = arg maxπ(θ|y) is found. Once the maximum is determined, the area around the mode is explored and a set of selected evaluation points {θ k } K k=1 is chosen, at whichπ(θ k |y) is evaluated to approximate the posterior marginals for each hyperparameter θ j using numerical integration. There is no global information available onπ(θ|y), instead it can only be evaluated for fixed values of θ. Each function evaluation is, however, a computationally expensive operation, mainly for two reasons. Let us for simplicity first consider the case where the likelihood is Gaussian. To evaluate the different terms of Eq. (3) for a given θ, we require the normalizing constant of the prior of the latent effects π(x|θ). As mentioned previously, x forms a GMRF, and hence the normalizing constant is computed from the determinant of its precision matrix. The computational complexity of computing the determinant of a dense matrix is in general O(n 3 ), where n denotes the size of the matrix [51]. For sparse matrices this cost can be lowered through employing specialized solution methods but it nevertheless continues to pose an expensive computational task for high-dimensional latent parameter spaces [12]. To evaluate π G (x|θ, y) for a given θ, one requires its normalizing constant as well as its mode which is dependent on θ. This requires the computation of the determinant and solving a sparse linear system of size n × n, respectively, which are again computationally expensive tasks. Fortunately, the sparsity pattern of the precision matrix of π G is maintained when conditioning on the data y which is an essential ingredient to the scalability of the INLA methodology [45]. Determining the remaining terms of Eq. (3) is comparatively cheap from a computational perspective. When considering the case, where the likelihood is non-Gaussian, we have to perform an inner or nested optimization, to obtain a good approximation π G (x|θ, y) [45]. This requires the repeated evaluation of its normalizing constant and mode during each iteration, thus further increasing the computational cost.

Approximating π(x i |y)
The posterior marginals π(x i |y) can be written as π(x i |y) = π(θ|y)π(x i |θ, y)dθ, for which an approximation to the first factor of the right-hand side and numerical integration points {θ k } K k=1 are already known from the previous step. It remains to estimate π(x i |θ, y). If the likelihood is normally distributed, the Gaussian approximation π G (x|θ, y) from Eq. (3) is exact. The posterior marginals of π(x i |θ, y) can then be approximated by extracting the latent parameter configuration x at the mode of π(θ|y) from which one can directly infer the marginal means x i |x −i , θ, y. Here, x −i denotes the vector x without the i-th entry. The marginal variances are computed from π G , again, evaluated at the modal configuration. This, however, requires at least a partial inversion of the precision matrix of π G , a computationally intensive operation especially for high-dimensional latent parameters spaces. In the non-Gaussian case, the approximation procedure is more involved. An update for the correction of the conditional mean of x i |x −i , θ, y for each i needs to be performed. Additionally the marginal variances are computed at all integration points {θ k } K k=1 which entails K partial matrix inversions. This allows INLA to more accurately account for the overall shape and potential skewness in the distributions. Using this information an approximation to the marginal distributions π(x i |θ, y) are constructed, which are then used in Eq. (4) to compute the posterior marginals using numerical integration. For a detailed overview we refer to [45]. Generally all computations are performed in log-scale as this is favorable for dependence structures and in terms of numerical robustness.

Sparse Linear Algebra
In this section we want to provide a quick overview of the underlying sparse linear algebra techniques leveraged by the R-INLA package. A sparse matrix is simply defined as a matrix where the majority of its elements are equal to zero. For sparse matrix algebra it is customary to only save the non-zero entries of the matrix in so-called compressed matrix formats, that store the non-zero values and their corresponding position instead of large quantities of zero entries. This allows to massively reduce the required memory and computations but hence, also demands for solvers that are especially tailored to sparse problems. These arise, however, very commonly and are therefore extensively researched, for introductory purposes see [12,48].
During the first stage of the INLA algorithm, see Sec. 2.2, the log determinant of different precision matrices Q (to obtain the corresponding normalizing constant) and the solution to linear systems of equations of the form Qx = b (to obtain the mode of the associated normal distribution) are required. Since the arising precision matrices are symmetric positive-(semi)definite, it is most efficient to perform a Cholesky decomposition, where the original matrix Q factors into Q = LL T , with L being a lower-triangular matrix [2]. The log determinant of Q can easily be computed from the diagonal entries of L as n i=1 2 log(L ii ) and the system Qx = b can quickly be solved for x using forward-backward substitution. That means first solving Ly = b for y and then L T x = y for x. These operations can be executed efficiently due to the lower-triangular structure of L.

Matrix Reordering
Additionally it is in most cases advantageous to apply a symmetric permutation to the matrix Q (and hence b, respectively) before computing the Cholesky decomposition. Even if the matrix Q is very sparse, L can have a large fill-in, meaning that there are many entries which are non-zero in L but equal to zero in Q. In the elimination process that is used to compute L, similar to classical Gaussian elimination, one can see how these non-zeros arise [12]. As an illustrative example, we consider a symmetric positive-definite arrowhead matrix Q 1 whose sparsity pattern can be seen in Fig. 1. Its corresponding Cholesky factor L 1 is computed recursively, starting from the first diagonal entry. To compute the off-diagonal entry (L 1 ) 32 , we use that ( are. More generally one can observe a dependency structure on previous columns of the i-th and jth row. By applying a symmetric permutation to Q 1 , we can obtain Q 2 . If we now compute L 2 we can see that (L 2 ) 32 = 0, since (Q 2 ) 32 , (L 2 ) 21 and (L 2 ) 32 are all equal to zero. Fig. 1: Sparsity patterns of the four-by-four symmetric positive-definite arrowhead matrices Q 1 and Q 2 , where each x stands for a non-zero entry and Q 2 is a symmetric permutation of Q 1 . The sparsity patterns of their Cholesky factors are represented by L 1 and L 2 , respectively. The matrix L 1 is dense i.e. exhibits a large fill-in while L 2 preserves the sparsity pattern of Q 2 .
In general, fill-in can often be drastically reduced by finding a suitable matrix reordering of Q which then lowers the overall memory and computation requirements of L.
For each of the different arising precision matrices in INLA the sparsity pattern never changes throughout the optimization phase but only the numerical values of the non-zero entries. Hence, it is sufficient to only compute suitable reorderings once throughout the entire algorithm. Finding a favorable permutation is, however, a challenging task and has hence been, an active area of research [8,19].
A simple greedy strategy is the minimum degree ordering, where columns are successively permuted to minimize the number of non-zero off-diagonal entries of the current pivot element [18]. Another popular class of reordering strategies called nested dissection employs graph partitioning techniques [8,17]. A symmetric matrix can be associated with an undirected graph, where each diagonal matrix entry represents a node and each nonzero off-diagonal entry an edge between the corresponding nodes. The graph is partitioned into two roughly equal-sized independent subgraphs G A , G B which are only connected through a separating set G S . The partition is chosen such that the size of the separating set is minimized. For the corresponding matrix this means that the associated submatrices Q A and Q B can be written as block matrices only connected through Q S . Hence, all fill-in that can occur is within Q A , Q B and Q S , see Fig. 2. For large graphs the computational cost is often too high to compute an optimal separating set. Therefore, multilevel strategies are used to recursively coarsen a large graph by merging connected vertices until it is reduced to a tractable size. For the coarse graph an edge cut can be defined which is then recursively refined while being propagated back to the original large graph. From the final edge cut we can deduce a separating set G S [22]. The same strategy of finding roughly balanced minimal separating sets can be independently reapplied to Q A and Q B , leading to the nested dissection approach. A popular state-of-the-art library that computes such heuristic matrix reorderings is called METIS [22] and is used by both, TAUCS as well as PAR-DISO. While the previously employed TAUCS library is limited to using fill-in reducing permutations, PARDISO can additionally perform much We can observe that Q was permuted such that it is separated in two independent submatrices Q A and Q B , which are only connected through entries that are now placed in the two rows, which we refer to as Q S . This way there is no fill-in between Q A and Q B but only within the before-mentioned submatrices. Thus, reducing the number of non-zeros in the Cholesky factorization of Q perm . The strategy can be applied recursively to each of the components Q A and Q B .
of the Cholesky factorization in parallel, making use of the independent submatrices arising from the nested dissection reordering. This will be explained in detail in Sec. 3.3.

Partial Inversion
In the second stage of the INLA algorithm, see Sec. 2.3, an estimate for the variances of the posterior marginals of the latent parameters, i.e. π(x i |y) for all i, is obtained. To do so INLA uses information from the Gaussian approximations π G (x|θ k , y) at the various integration points {θ k } K k=1 . All information about the marginal variances only exists, however, in form of the precision matrices. As the variances Σ ii correspond to (Q −1 ) ii , they can only be obtained through matrix inversions. If the full precision matrices had to be inverted this would become very time and memory consuming if not infeasible for high-dimensional latent parameter spaces as matrix inversions have a complexity of O(n 3 ) for a matrix of size n × n. Fortunately there is an alternative method for computing the marginal variances that is much more efficient. There are two different approaches to derive this recursive strategy. One is using the conditional independence properties of GMRFs and their associated graphical structure, and was first described by Rue and Martino in [44]. The other one is based on a particular way of writing matrix identities without a further interpretation but describing the same recursion and was developed by Takahashi [53] almost 50 years ago. We will begin by considering the former. The solution x to the problem L T x = z where z ∼ N (0, I) is a sample from a GMRF with zero mean and precision matrix Q = LL T [43]. Since L T is an upper triangular matrix, we can use a backward solve and recursively compute for i = n − 1, ..., 1. We exclude all terms L ki from the summation directly that equate to zero. This way it becomes clear that the more zeros we have in L, the less computations are necessary. If we use the known expected value and variance of z as well as Eq. (5), we can deduce the covariance matrix Σ of x from first principles, obtaining where δ ij is one if i = j and zero otherwise. The entries can only be computed recursively starting at i = n, traversing the matrix from the bottom right to the top left.
Equivalently we can derive Eq. (6) without using statistical properties. Instead we consider the slightly altered decomposition Q = LL T = V DV T , i.e. L = V D 1/2 , where D is a diagonal matrix and V a lower triangular matrix with ones on the diagonal. The following matrix identity was proposed by Takahashi in [53], As Σ is a symmetric matrix, it is enough to compute its upper triangular part. The term since V has a unit diagonal. Writing out Eq. (7) as sums we obtain which we can compute recursively starting from Σ nn and where again, δ ij is equal to one for i = j and zero otherwise. Eq. (6) and (8) are equal since L = V D 1/2 . It is possible to use these recursions to compute all entries of Σ, however, in this case they do not give us a computational advantage over traditional inversion algorithms. In both cases we have a complexity of O(n 3 ). If we are instead only interested in particular entries of Σ, e.g. only the diagonal, and if additionally L (and then likewise V ) are sparse, a tremendous amount of computational cost can be saved, as only the entries Σ ij for which L ij is non-zero need to be computed. Hence, the less non-zeros we have in the factor L, the smaller is the computational need. This highlights the importance of finding suitable permutations as described in the previous section. PARDISO and INLA's previous selected inversion routine both employ such a partial inversion strategy. While the latter is sequential, PAR-DISO is using parallelized computations for a shorter time to solution, exploiting the particular matrix structure given through the permutation, see Sec. 3.3 for details. The marginal variances Σ ii that we obtain are then repermuted to correspond to the original ordering of the parameters x. Afterwards, they can be used as described in Sec. 2.3 to obtain estimates of the marginal posterior distributions of the latent parameters.

New Parallelization Scheme
In the following we present the novel parallelization strategies that have been added to the current R-INLA implementation. There are four significant updates in the algorithm that we want to discuss in this work. The first two comprise parallelized versions of the gradient computations and the posterior marginal variance estimations. Additionally, we introduced a parallel line search strategy within the optimization routine of the algorithm. The parallelizations are realized through the usage of OpenMP, an application programming interface that supports multiprocessing for shared-memory architectures and allows to efficiently execute multiple tasks at once. The last update applies to all stages of the algorithm and consists of the incorporation of the PARDISO library into R-INLA to handle all required major sparse linear algebra operations.

Parallelization of Function Evaluations
The logarithmic scale is introduced to enhance numeric stability. The sign switch turns f into a minimization problem to find the mode θ * of the posterior. To solve the optimization problem INLA uses a BFGS-algorithm [37] which resembles a Newton method, without requiring knowledge about the second order derivative. Information on the gradient is, however, needed. We estimate it numerically using a finite difference scheme [26]. Each directional derivative ∂f /∂θ i in the gradient is approximated using either a first order forward or central difference. As an example we show the central difference approximation along the coordinate axes The vector ǫ i is of the same dimension as θ and contains only zeros except for the i-th component which contains a small value ǫ > 0.
Often the directional derivatives are computed according to this canonical basis, however, noncanonical bases can just as well be used. INLA makes use of knowledge from previous iterations to choose directional derivatives exhibiting more robust numerical properties and hence faster overall convergence, see [16] for details. Independently of the choice of basis, the directional derivatives are computed for each component of θ and each time entail one or two function evaluations of f . We can, however, see from Eq. (10) that all function evaluations are independent from each other.
In the new parallelization scheme we are using this to our advantage by computing the function values, i.e. f (θ + ǫ 1 ), f (θ − ǫ 1 ), f (θ + ǫ 2 ), ... simultaneously in each iteration. For a central difference scheme we need two times as many function evaluations as there are hyperparameters, which can now all computed in parallel instead of sequentially, while also computing f (θ) itself. So, if e.g. dim(θ) = 3 we can have a theoretical speedup of 7 during the gradient computation. This means the introduced parallelism allows the gradient to be computed at almost no further cost in addition to the already necessary iterative evaluations of f (θ) for fixed values of θ. There are other methods to obtain gradient information, in particular automatic differentiation [6]. In the ideal case it offers a highly accurate solution at runtimes comparable to a single function evaluation [6], similar to the parallelized finite difference approximation. It requires, however, an adaptation of the basic linear algebra operations and then allows for parallelism within these operations. This is not trivially done for non-standard cases [55] and would require major changes in the PARDISO software library. Additionally, if the likelihood is non-Gaussian we have the inner optimization routine within each function evaluation, for which we are not aware of existing automatic differentiation methods.
After the mode ofπ(θ|y) is found, we can also perform the function evaluations of the integration points {θ k } K k=1 in parallel. Additionally, we have to estimate the posterior marginal variances, computing the partial inverses of the precision matrices at the integration points, see Sec. 2.3. We also implemented this to be executed in parallel, as these computations are also independent from each other.

Parallel Line Search
During the l-th iteration of the optimization routine we approximate the gradient ∇f θ l at θ l . If a full Newton scheme was employed, one would additionally have to compute the Hessian of f . When using a BFGS algorithm, an approximation to the Hessian B l is formed in each iteration, utilizing information from the previous gradients. Each next iterate is computed as where 0 < α l < 1 is the step size in iteration l. Determining a suitable value for α l is crucial to the convergence of the overall algorithm and is done through a line search procedure, see e.g. [37] for an overview. The most common approach is to use a mixture of artificially chosen upper and lower bounds, as well as the value from the previous iteration, to propose the next step size α l . Then, a check is performed to see if a certain condition is met to either accept or reject the suggested step size. If it is rejected, the step size is reduced, and the check is performed again, if it is accepted, θ l is updated as described in Eq. (11) and one proceeds to the next iteration. For this check, the function f is evaluated at the potential new candidate which, as discussed previously, is an expensive operation to perform, but necessary in order to determine the validity of the new step. Hence, every time a step is rejected another sequentially computed function evaluation is required. To be more efficient we introduce a parallel line search strategy to make better use of the available resources. Instead of computing different values for α l sequentially, we define a search interval I = (θ l , θ l − γ l p l ), where p l := [B l (f θ l )] −1 ∇f θ l and γ l > 0 is an upper bound to α l . Hence, I contains all possible solutions of Eq. (11) for 0 < α l ≤ γ l . Depending on the number of available cores we define a number of points θ li on I for which we evaluate all f (θ li ) in parallel. The easiest option would be to now choose the candidate θ li that minimizes f as the next iterate. However, it has shown to be advantageous to fit a second order polynomial q through the newly evaluated points using robust regression [3,42]. This way slight inaccuracies in the function evaluations (that can e.g. arise from more complicated choices of likelihood and require an inner optimization loop) are counter balanced. We additionally add two more evaluation points in positive p l direction close to θ l to stabilize the polynomial fitting process close to the global optimum. This does not increase the overall runtime as they can be evaluated in parallel with the other θ li . Robust regression differs from regular regression in the sense that each pair (θ li , f (θ li )) gets assigned a weight w θ l i , which will make it more or less influential on the overall fitting process.
There are a number of commonly used weighting functions w. We have chosen to use the so-called bisquare weighting as this has numerically shown to be most suitable. After finding the second order polynomial q its minimum is determined on I and chosen as the next iterate θ l+1 . In Fig. 3 we can see an illustration of the new strategy. The advantages are the mitigation of inherently sequential function evaluations in the reject-accept check, as well as an improved step size choice. This can lower the number of required iterations until convergence, and hence also the overall runtime, see Sec. 4 for numerical results.

Parallelization within PARDISO
The PARDISO library offers routines to efficiently provide solutions to various problems commonly arising in the field of sparse linear algebra [14]. In R-INLA the PARDISO implementation for Cholesky decomposition, subsequent forward-backward substitution and partial inversion are invoked. PARDISO is a state-of-the-art direct solver that has been in active development for many years. To achieve excellent performance it relies on numerous strategies ranging from thorough algorithm selection to detailed hardware considerations.
In this work we want to solely focus on the employed parallelism and refer the interested reader to [9] for a comprehensive overview.
As described in Sec. 2.4.1, a symmetric permutation found through nested dissection is applied to Q before computing its Cholesky decomposition L. This reordering technique does not only reduce the arising fill-in but also creates independent block matrices Q A , Q B , see Fig. 2, that are only connected through a small submatrix Q S . Hence, it is possible to factorize Q A and Q B in parallel, before factorizing Q S . The reordering generated through nested dissection has a recursive pattern, and thus Q A and Q B themselves exhibit the same structure. Both of them, again, contain two independent blocks and a connecting submatrix which corresponds to the separating set of the associated graph. The reordered matrix can therefore be factorized in parallel, starting from the set of smallest block matrices. If sufficient compute resources are a available this introduction of parallelism drastically reduces the computational time and especially exhibits much better scaling. While it is not possible to give precise bounds on timings without specifying more graph theoretical concepts and notation, one can see that the recursive subdivisions induce a logarithmic growth. For details see [39].
A similar principle can be applied to determine the partial inverse of Q, for which we have seen in Sec. 2.4.2 that only the non-zero entries of L are required. We can compute the same subgraphs in parallel, as in the Cholesky decomposition. This time, however, we first compute the values of the submatrix connected to the separating set, and referred to as Q S in Fig. 2. Once we have recursively computed all necessary inverse elements belonging to the indices of Q S , we can determine the inverse elements of Σ with indices belonging to Q A and Q B in parallel. This is possible since Q A and Q B are only connected through Q S and these entries have already been computed at this point. Since for larger matrices, each submatrix Q A , Q B was again permuted following the same strategy, we can recursively traverse the matrix, in the opposite direction as in the Cholesky decomposition, computing the required entries of Σ in many parallel regions [9].
PARDISO uses OpenMP to carry out the simultaneous computations. Thus, we obtain a nested OpenMP structure, as we have multiple matrix factorizations or inversions carried out in parallel, see Sec. 3.1 and 3.2. Additionally to the parallelism within PARDISO, OpenMP is used on the second level to parallelize matrix operations e.g. during the computation of matrix-matrix or matrix-vector products and while assembling matrices.

Results & Benchmarks
In this section we present performance results of the newly introduced parallelization strategies for various applications. They are based on previously published models that make use of the R-INLA package. We will not go into the details of the individual application as all the details can be found in the original publications, and instead focus on the computational aspects involved. All numerical experiments for Case Study I & II were performed on a single node machine with 755 GB of main memory and 26 dual-socket Intel(R) Xeon(R) Gold 6230R CPU @ 2.10GHz, totaling 52 cores. The large number of cores allows us to demonstrate the full performance gains that can be obtained through parallelism. All numerical experiments for Case Study III were performed on an Apple M1 Mac mini with 16 GB of memory, 4 performance and 4 efficiency cores. We chose this machine to illustrate that we also obtain performance gains through parallelism on regular desktop computers, laptops and notebooks, although of course to a smaller extent than compared to a larger computer architecture.
Case Study I : Joint survival modeling of randomized clinical trial. [47], with R package INLAjoint 2 , a joint survival model consisting of multivariate longitudinal markers paired with competing risks of events. The different submodels are linked to each other through shared or correlated random effects. The authors consider a randomized placebo controlled trial for the treatment of the rare autoimmune disease primary biliary cholangitis (PBC) which affects the liver. Their model has a particularly large number of hyperparamters, d(θ) = 50, due to the 7 different correlated likelihoods.   [52] that functional brain responses can be reliably estimated using cortical surface-based spatial Bayesian generalized linear models (GLMs). Functional magnetic resonance imaging (fMRI) data from individual subjects is used to identify areas of significant activation during a task or stimulus. The authors use the stochastic partial differential equations approach [25,30,31] for manifolds to define a spatial GMRF field on the surface of the brain employing geodesic distances. This allows to flexibly encode spatial dependency structures in the model while maintaining sparsity in the precision matrices. Their approach defines a GLM that fits the framework presented in 2.1. The likelihood is assumed to follow a Gaussian distribution and we therefore do not require an inner optimization loop for every function evaluation. More details can be found in [35,52] which also includes information about the corresponding R package BayesfMRI 3 . For this work we fit a single subject, single repetition set up using fMRI data for 4 tasks. The dimension of the hyperparameter space d(θ) = 9, two for the parametrization of the spatial field of each task and a noise term for the observations. The images are preprocessed for standardization and the surface of the brain is discretized using a triangular mesh. We use two different spatial discretizations that are determined by the resampling size, which subsequently influence the dimension of the latent parameter space and the number of considered observations.  As a final example we consider a study on spatial variation in Leukemia survival data from [20,31]. The additive linear predictor of the hazard function includes 5 fixed effects and a spatial GMRF random field, which induces a three-dimensional hyperparamter space. The fitted model provides

Leveraging Level 1 Parallelism
The level 1 parallelism enables simultaneous function evaluations as described in Sec. 3.1 which are relevant for the gradient computation, the final Hessian approximation during the optimization phase as well as the partial inversion scheme at the mode. The speedup that can theoretically be obtained through parallelization depends on the number of hyperparamters. If a forward difference scheme is employed a maximum speedup of dim(θ)+1 is possible during the gradient computation, for a central difference scheme it is 2 · dim(θ). If not specified otherwise R-INLA employs a mix of forward and central difference approximations. The theoretical speedup during the Hessian approximation and the partial inversion scheme depends on the number of integration points which in turn depends on the chosen integration scheme as well as the dimension of θ. While these are the computationally most costly operations, there are other nonparallelizable steps required. Additionally, setting up a parallel OpenMP environment always creates overhead. Hence, the theoretical speedup can almost never be attained. For relatively small applications the observed speedup through parallelism is usually not as significant as for larger problems, because the sequential parts and the induced overhead make up a larger part of the total time. On the other hand, larger applications using a very large number of threads might not exhibit further speed up after a certain thread count as the memory bandwidth can become the limiting factor. We will be able to observe both effects when looking at performance results of parallelizing level 1. In Fig. 4  When employing more than 20 threads the parallel efficiency starts to reduce. Case Study II has less hyperparameters, hence the maximum attainable speedup on level 1 is lower, nevertheless providing significant improvements for up to 10 threads with a speedup of almost 6 compared to the single-threaded version. We also analyse the effects of level 1 parallelism for Case Study III. This model is of much smaller dimension than the first two and hence a larger part of the overall runtime is dedicated to sequential operations like setting up or the initializations of the parallel regions. Nevertheless we observe a speedup of a factor of 2 when using 4 instead of 1 thread.

Leveraging Parallel Line Search
The accuracy of robust regression clearly depends on the number of available regressors, which in     turn is dependent on the number of available threads. In our experience making use of the parallel line search implementation only becomes advantageous if 8 threads or more are in use on level 1. Then there are sufficiently many evaluation points in the search interval I to adequately represent the original function. The parallel linesearch becomes especially relevant for models with non-Gaussian likelihoods as in this case the inaccuracies in the function evaluations increase due to the arising inner iteration. In Fig. 5, we show the total runtime of Case Study I using varying numbers of threads on level 1 with and without enabling the parallel line search. Since this model has a large number of hyperparameters the computation of the posterior marginals of the latent parameters using the simplified Laplace approximation strategy would be the dominating the overall compute time. In this section we, therefore, used the computationally cheaper empirical Bayes' approximation, to be able to put the emphasis on the optimization phase of the algorithm, where the parallel line search feature is relevant. We can see that it lowers the overall time to solution without requiring more threads. It additionally exhibits a more stable behavior when it comes to reduction in runtime over an increasing thread count.

Leveraging Level 2 Parallelism
The level 2 parallelism occurs within each call of the PARDISO library as well as other matrix operations like matrix-matrix multiplications, following the concepts described in Sec. 3.3. The parallelization on level 2 only starts to make a noticeable difference for larger latent parameter spaces. The single-threaded version of the PAR-DISO library is already very efficient for small to moderate sized problems. The full effects of the parallelization start becoming visible for the latent parameter spaces of dimensions in the range of 10 4 and larger, especially for models with a three-dimensional structure.
We consider both examples of Case Study II in the performance analysis of Fig. 6. For a fixed number of threads on level 1, we observe a continuous speedup for increasing numbers of threads on level 2 up until 8 threads for the medium-sized Case Study II, leading to a maximum speedup of almost 3 over the single-threaded version. In turn the larger-sized Case Study II continues to exhibit a speedup also beyond 8 threads, showing that larger latent parameter space benefit more from level 2 parallelism, leading to a speedup of almost 5 at 16 threads over the single-threaded version.

Leveraging combined Parallelism
In this section we want to present results for combined parallelization strategies and discuss how to choose a favorable set up. The ideal approach clearly depends on the available infrastructure and the problem at hand. In almost all cases the number of available cores will be limited. Often this limit will easily be reached as the total number of threads is equal to (# threads L1) · (# threads L2). This poses the question of how best utilize the available resources. We have found that hyperthreading often does not significantly increase the performance (or is even counter productive) and therefore only consider thread counts that are within the number of physical cores. In the previous two sections we have seen that smaller problems benefit less from larger thread counts as the overhead of thread initialization can overshadow the gain. We can see from both Fig. 4 & 6 that the parallel efficiency is highest for smaller thread counts on both levels. This implies that its often favorable to distribute the available number of threads among both levels. Additionally the structure of the model should be taking into account. If the number of hyperparamters is large it is advantageous to focus more on level 1 parallelism whereas a very high-dimensional latent parameter space will benefit more from level 2 parallelism. This behavior can be observed empirically in Fig. 7, where we present the runtime of the INLA algorithm in seconds for different thread configurations. It becomes clear that the introduction of parallelism leads to a tremendous reduction in runtime. For Case Study I, we observe speedups of a factor 20 and more over the single-threaded version. For the medium-sized Case Study II we observe a speedup of a factor up to 15. While for the large Case Study II we observe a speedup of more than 10 times compared to the singlethreaded version. For sufficiently large multi-core architectures, the newly introduced updates open the door to previously unfeasible modeling scales, and drastically reduce waiting times for users.

Discussion
In this paper we presented novel parallelization strategies for the INLA methodology. While INLA has always been using optimized sparse linear algebra operations to be computationally efficient, we have taken this a step further by introducing parallelism using OpenMP. This allows for the simultaneous execution of computationally expensive tasks which are distributed over a specified number of cores. Our approach contains two layers of parallelism. On the first layer we have parallelized independent function evaluations during INLA's optimization phase to determine and explore the posterior mode of the hyperparameters. Each evaluation requires the Cholesky factorization, a computationally expensive task, of at least one sparse precision matrix whose size is tied to the potentially very high-dimensional latent parameter space of the model. We have additionally introduced a parallel line search routine using robust regression to parallelize and stabilize the search for each next step size, thus being able to lower the overall number of required iterations. We were also able to parallelize large parts of the computation of the posterior marginal variances of the latent parameters which require a computationally expensive partial matrix inversion. The state-of-the-art sparse solver PARDISO was included in INLA to handle the most expensive linear algebra operations and makes up the second layer of parallelism. PARDISO makes use of intelligent matrix reordering techniques that allow large parts of the otherwise sequential computations to be parallelized. Thus, PARDISO internally employs OpenMP for simultaneous operations within a Cholesky factorization, the solution to a linear system of equations or a partial matrix inversion. We empirically demonstrated the performance and scalability of the parallelization strategies on three case studies. Each of them uses   a different type of model, hence posing different computational challenges. For all larger applications we observed speeds up of a factor of 10-25 times compared to the single-threaded version. Hence, lowering runtimes for some models from almost seven hours to less than 30 minutes or from more than two hours to just over ten minutes. For smaller applications we still observe speedups but to a smaller extend. Given that larger multi-core computer architectures are becoming more and more accessible, these computational advancements allow users to perform inference using more complex models at higher resolution in shorter runtimes with INLA.