{'fraction_non_alphanumeric': 0.03938611304261448, 'fraction_numerical': 0.017979537101420483, 'mean_word_length': 4.338016967126193, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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}}
A Fast Linear Regression via SVD and Marginalization

We describe a numerical scheme for evaluating the posterior moments of Bayesian linear regression models with partial pooling of the coefficients. The principal analytical tool of the evaluation is a change of basis from coefficient space to the space of singular vectors of the matrix of predictors. After this change of basis and an analytical integration, we reduce the problem of finding moments of a density over k + m dimensions, to finding moments of an m-dimensional density, where k is the number of coefficients and k + m is the dimension of the posterior. Moments can then be computed using, for example, MCMC, the trapezoid rule, or adaptive Gaussian quadrature. An evaluation of the SVD of the matrix of predictors is the dominant computational cost and is performed once during the precomputation stage. We demonstrate numerical results of the algorithm. The scheme described in this paper generalizes naturally to multilevel and multi-group hierarchical regression models where normal-normal parameters appear.


Introduction
Linear regression is a ubiquitous tool for statistical modeling in a range of applications including social sciences, epidemiology, biochemistry, and environmental sciences ( [Gelman et al., 2013, Gelman and Hill, 2007, Greenland, 2000, Merlo et al., 2005, Bardini et al., 2017).
A common bottleneck for applied statistical modeling workflow is the computational cost of model evaluation. Since posterior distributions in statistical models are often high dimensional and computationally intractable, various techniques have been used to approximate posterior moments. Standard approaches often involve a variety of techniques including Markov chain Monte Carlo (MCMC) or using a suitable approximation of the posterior.
In this paper, we describe an approach for reducing the computational costs for a particular class of regression models -those that contain parameters θ ∈ R k such that θ has a normal prior and normal likelihood. These models represent only a subset of regression models that appear in applications. We focus our attention in this paper on normal-normal models because they have well known analytical properties and are more computationally tractable than the vast majority of multilevel models. A broader class of models, including logistic regression, contain distributions that are less amenable to the techniques of this paper and will require other analytical and computational tools. Mathematically, marginalization of normal-normal parameters is well-known and has been applied to the posterior by, for example, [Lindley and Smith, 1972]. Our contribution is to provide a stable, accurate, and fast algorithm for marginalization. The primary numerical tool used in the algorithm is the singular value decomposition (SVD) of the data matrix. As a mathematical and statistical tool, SVD has been known since at least 1936 (see [Eckart and Young, 1936]). Use of the SVD as a practical and efficient numerical algorithm only started gaining popularity much later, with the first widely used scheme introduced in [Golub and Kahan, 1965]. Due in large part to advances in computing power, use of the SVD as a tool in applied mathematics, statistics, and data science has been gaining significant popularity in recent years, however efficient evaluation of SVDs and related matrix decompositions is still an active area of research (see [Hastie et al. 2015], [Halko et al., 2011], [Shamir et al., 2016]).
Similar schemes to ours are used in the software packages lme4 ( [Bates et al., 2015]) and INLA ( [Rue et al., 2017]). There are several differences between the problems they address and their computational techniques, and those that we shall discuss here. While lme4 finds maximum likelihood and restricted maximum likelihood estimates, our goal is to find posterior moments. The software package INLA uses Laplace approximation on the posterior for a general choice of likelihood functions, whereas our algorithm is focused on fast and accurate solutions for only a particular class of densities: those with normal-normal parameters.
The approach presented in this paper analytically marginalizes the normal-normal parameters of a model using a change of variables. After marginalization, posterior moments can be computed using standard techniques on the lower-dimensional density. In particular, for a model that contains k + m total variables, k of which are normalnormal, our scheme converts the problem of evaluating expectations of a density in k + m dimensions to finding expectations of an m-dimensional density. After marginalization, we evaluate the m-dimensional posterior density in O(k) operations. Without the change of variables, standard evaluation of marginal densities that relies on determinant evaluation requires at least O(k 3 ) operations.
We illustrate our scheme on the problem of evaluating the marginal expectations of the unnormalized density where γ > 0 is a constant, σ 1 , σ 2 > 0, and β ∈ R k . We assume that X is a fixed n × k matrix, y ∈ R n is fixed, and the normalizing constant of (1) is unknown. For fixed n, k ∈ N, the algorithm is nearly identical when X is an n × k matrix to when X is a k × n matrix. In the case where k n, see [Kwon et al., 2011] for a similar approach. Using the standard notation of Bayesian models, density q is the unnormalized posterior of the model In Appendix A, we include Stan code that can be used to sample from density (1). Statistical model (2) is a standard model of Bayesian statistics and appears when seeking to model an outcome, y, as a linear combination of related predictors, the columns of X. In [Gelman and Hill, 2007], these models are described in detail and are used in the estimation of the distribution of radon levels in houses in Minnesota.
Density (1) is also closely related to posterior densities that appear in genomewide association studies (GWAS; see [Zhu and Stephens, 2017], [Meuwissen, et al., 2001], [Azevedo et al., 2015]) which can be used to identify genomic regions containing genes linked with a specific trait, such as height. Using the notation of (1), each row of matrix X corresponds to a person, each column of X represents a genomic location, entries of X indicate genotypes, and y corresponds to the trait. Due to technical advances in genome sequencing over the last ten years, it is now feasible to collect large amounts of sequencing data. GWAS models can contain data on up to millions of people and often between hundreds and thousands of genome locations (see [Linner et al., 2019]). As a result, efficient computational tools are required for model evaluation.
The number of operations required by the scheme of this paper scales like O(nk 2 ) with a small constant. The key analytical tool is a change of variables of β such that the terms, in (1) are converted to a diagonal quadratic form in R k . After that change of variables, expectations over q are analytically converted from integrals over R k+2 to integrals over R 2 . The remaining 2-dimensional integrals can be computed to high accuracy using classical numerical techniques including, for example, adaptive Gaussian quadrature or even the 2-dimensional trapezoid rule. The schemes used to evaluate the expectations of (1) generalize naturally to evaluation of expectations of multilevel and multigroup posterior distributions including, for example, the two-group posterior of the form, where X is a n × k matrix, y ∈ R n , k 1 and k 2 are non-negative integers satisfying k 1 + k 2 = k, and the vector t ∈ R m .
For models where m is large, MCMC can be used to evaluate the m-dimensional expectations, with, for example, Stan [Carpenter et al., 2017]. The m-dimensional distribution has two qualities that make it preferable to its high-dimensional counterpart. First, it requires only O(k) operations to evaluate the integrand, and second, the geometry of the m dimensional marginal distribution will allow for more efficient sampling.
The structure of this paper is as follows. In the following section we describe the analytic integration that transforms (1) from a k + 2-dimensional problem to a 2-dimensional problem. Section 3 includes formulas that will allow for the evaluation of posterior moments using the 2-dimensional density. In Sections 4 and 5 we provide formulas for evaluating covariances of (1). In Section 6, we discuss the numerical results of the implementation of the algorithm. Conclusions and generalizations of the algorithm of this paper are presented in Section 7. Appendix A provides Stan code that can be used to sample from (1), and Appendix B includes proofs of the formulas of this paper.

Analytic Integration of β
In this section, we describe how we analytically marginalize the normal-normal parameter β of density (1). We include proofs of all formulas in Appendix B.
We start by in marginalizing β using a change of variables that converts the quadratic forms in (1) into diagonal quadratic forms. The resulting integral in the new variable, z, is Gaussian, and the coefficients of z i and z 2 i are available analytically. The change of variables is given by the right orthogonal matrix of the singular value decomposition (SVD) of X. That is, we set where the SVD of X is We define λ i to be the i th element of the diagonal of D. The elements of diagonal need not be sorted. After this change of variables, we obtain the following identity for the last two terms of (1). A proof can be found in Lemma B.5 in Appendix B.
Formula 2.1. where and After performing the change of variables z = V t β and using (7), we now have an expression for density (1) in a form that allows us to use the well-known properties of a Gaussian with diagonal covariance. The following identity uses these properties and provides a formula for analytically reducing expectations of (1) from integrals over k + 2 dimensions to integrals over 2 dimensions. After the formula is applied, we have a new density,q, over only 2 dimensions. See Theorem B.6 in Appendix B for a proof.
where a 2,i is defined in (8), a 1,i is defined in (9), a 0 is defined in (10), and γ is a constant.
Remark 2.1. Certain Bayesian models might contain correlated priors on β that will result in posteriors such as (28) of Section 4. For such models, we perform the change of variables that uses the fact that two diagonal forms over β can be simultaneously diagonalized.
We include in Figure 1 a plot of the density of q as a function of σ 1 and β 1 for fixed σ 2 and randomly chosen X and y. Figure 2 shows a plot of q as a function of σ 2 and β for fixed σ 1 . Figure 3 provides an illustration ofq, obtained after the change of variables and marginalization described in this section.

Evaluation of Posterior Means
Now that we have reduced the k + 2-dimensional density q to the 2-dimensional densityq, it remains to recover the posterior moments of q usingq. We first observe that moments of σ 1 and σ 2 with respect to q are equivalent to moments of σ 1 and σ 2 overq. That is, and E q (σ 2 ) = Eq(σ 2 ).
As for moments of β, we use (13) and standard properties of Gaussians to obtain the following formula.
where q is defined in (1),q is defined in (13), a 2,i is defined in (8), and a 1,i is defined in (9).
As an immediate consequence of (16), we are able to evaluate the posterior expectation of z as an expectation of a 2-dimensional density:  (1)) with respect to σ 1 and β 1 , where γ = 8, n = 100, k = 10, and data were randomly generated.  (1)) with respect to σ 2 and β 1 , for the same parameters as Figure 1.
We then transform those expectations back to expectations over the desired basis, β using the matrix V computed in (6). Specifically, using linearity of expectation and (17) (13)) using the same q as Figure 1, where n = 100, k = 10, and data were randomly generated.

Covariance of β
In addition to facilitating the rapid evaluation of posterior means, the change of variables described in Section 2 is also useful for the evaluation of higher moments. Equation (7) shows that after the change of variables from β to z, the resulting density is a Gaussian in z with a diagonal covariance matrix. Additionally, for each z i , using equation (7) and standard properties of Gaussians, we have the following identity.
where µ z i is the expectation of z i ,q is defined in (13), and a 2,i is defined in (8).
The second moments of the posterior of β are obtained as a linear transformation of the posterior variances of z. In particular, denoting the expectation of β by µ β and the expectation of z by µ z , we have We observe that due to the independence of all z i , is diagonal and we can therefore evaluate the k × k posterior covariance matrix of β by evaluating var(z i ) for i = 1, ..., k and then applying two orthogonal matrices. Specifically, combining Formula 4.1 and (20), we obtain 5 Variance of σ 1 and σ 2 Higher moments of σ 1 and σ 2 with respect to q can be evaluated directly as higher moments of σ 1 and σ 2 with respect toq. That is, for all j ∈ {2, 3, ..., }, we have and E q ((σ 2 − µ σ 2 ) j ) = Eq((σ 2 − µ σ 2 ) j ).

Numerical Experiments
Algorithm 1 was implemented in Fortran. We used the GFortran compiler on a 2.6 GHz 6-Core Intel Core i7 MacBook Pro. All examples were run in double precision arithmetic.
The matrix X and vector y were randomly generated as follows. Each entry of X was generated with an independent Gaussian with mean 0 and variance 1. The vector y was created by first randomly generating a vector β ∈ R k , each entry of which is an independent Gaussian with mean 0 and variance 1. The vector y was set to the value of Xβ + where ∈ R n contains standard normal iid entries. We generated y this way in order to ensure that the E(β i ) were not all small in magnitude. We set γ of (1) to 8.
In Table 2 and Figure 5, we compare the performance of Algorithm 1 to two alternative schemes for computing posterior expectations -one in which we analytically marginalize via equation (12) and then integrate the 2-dimensional density via MCMC using Stan. In the other, we use Stan's MCMC sampling on the full k + 2 dimensional posterior. When using MCMC with Stan, we took 10,000 posterior draws. In Table 2 and Figure 5 we denote Algorithm 1 by "SVD-Trap". The algorithm that uses Stan on the marginal 2dimensional density is labeled "SVD-MCMC", and "MCMC" corresponds to the algorithm that uses only MCMC sampling in Stan.
In the appendix, we include Stan code to sample from the marginal densityq of (13).
Remark 6.1. In the numerical integration stage of algorithm 1, we use the trapezoid rule with 200 nodes in each direction. Because the integrand is smooth and vanishes near the boundary, convergence of the integral is super-algebraic when using the trapezoid rule (see [Stoer and Bulirsch, 1992]). A rectangular grid with 200 points in each direction is satisfactory for obtaining approximately double precision accuracy. In problems with large numbers of non-normal-normal parameters, MCMC algorithms such as Hamiltonian Monte Carlo or other methods can be used.
In Tables 1 and 2, n and k represent the size of the n × k random matrix X. The column labeled "max error" provides the maximum absolute error of the expectations of σ 1 , σ 2 , and β i for i ∈ {1, 2, . . . , n}. The true solution was evaluated using trapezoid rule with 500 nodes in each direction in extended precision.
In Table 1, "Precompute time (s)" denotes the time in seconds of all computations until numerical integration. These times are dominated by the cost of SVD (36). The total time of the numerical integration in addition to the matrix-vector product (18) is given in "integrate time (s)." The final column of Table 1, "total time (s)", provides the total time of precomputation and integration.

Generalizations and Conclusions
In this paper, we present a numerical scheme for the evaluation of the expectations of a particular class of distributions that appear in Bayesian statistics; posterior dsitributions of linear regression problems with normal-normal parameters.
The scheme presented generalizes naturally to several classes of distributions that appear frequently in Bayesian statistics. We list several examples of posteriors whose expectations can be evaluated using this method.  (1) Table 2: Scaling of computation times for evaluation of expectations of q (see (1)) using three different algorithms: i) SVD-Trap: Algorithm 1 of this paper, ii) SVD-MCMC: marginalization with MCMC integration ofq using Stan, and iii) MCMC: full MCMC integration of q using Stan. 1. The choice of priors for σ 1 , and σ 2 in this document were log normal and half-normal. This choice did not substantially impact the algorithm and can be generalized. Adaptive Gaussian quadrature can be used for the numerical integration step of the algorithm for a more general choice of prior on σ 1 and σ 2 .
2. Multilevel regression problems with more than two levels.
3. Regression problems with multiple groups such as the two-group model with posterior where X is a n × k matrix, y ∈ R n , and k 1 and k 2 are non-negative integers satisfying 3. Regression problems with correlated priors on β: For regression problems with large numbers of non-normal-normal parameters, marginal expectations can be computed using, for example, MCMC in Stan. For such problems, the algorithm of this paper would convert an MCMC evaluation from k + m dimensions to m dimensions, where k is the number of normal-normal parameters.

Acknowledgements
The authors are grateful to Ben Bales and Mitzi Morris for useful discussions.

A Code
The following Stan code allows for sampling from the distribution corresponding to the probability density function proportional to (1). y~normal(X*beta, sigma2); beta~normal(0, sigma1); sigma1~lognormal(0, 0.25); sigma2~normal(0, 1); } The following Stan program samples from the marginal densityq (see (13)). The data input yty corresponds to y t y of (10), lam is the vector of singular values of X, and w is the vector w in equation (11). We include R code for computing yty, lam, and w after the following Stan code.

B Proofs
In this appendix, we include proofs of the formulas provided in this paper. For increased readability, this appendix is self-contained.

B.1 Mathematical Preliminaries and Notation
In this section, we introduce notation and elementary mathematical identities that will be used throughout the remainder of this section.
We will be denoting by 1 the vector of ones We denote the i th component of a vector v by v i . The following two well-known identities give the normalizing constant and expectation of a Gaussian distribution.

B.2 Analytic Integration of β
We denote the singular value decomposition of X by where U is an orthogonal n × k matrix, V is an orthogonal k × k matrix, and D is a k × k diagonal matrix. We define z ∈ R k by the formula The following lemma, which will be used in the proof of Lemma B.5, gives an expression for the second to last term of the exponent in (1) after a change of variables.
Lemma B.3. For all β ∈ R k , and y ∈ R n , z is defined in (37), and λ i is the i th entry on the diagonal of D (see (36)).
The following lemma provides an equation for the last term of the exponent in (1). The identity will be used in Lemma B.5.
where β ∈ R k , z is defined in (37), and V is defined in (36).
The following formula combines Lemma B.3 and Lemma B.4 to convert the final two terms of (1) into a Gaussian in k dimensions. where and a 0 = − y t y 2σ 2 2 where z is defined in (37), w is defined in (39) and V is defined in (36). Xβ − y 2 − 1 2σ 2 1 β 2 = a 0 + k i=1 a 1,i z i − a 2,i z 2 i .
We obtain equation (45) by completing the square in equation (49).
The following theorem is the principal analytical apparatus of this note. It provides a formula for the k-dimensional integrals that appear in (29), (30), and (31).
The following theorem provides a formula for the expectation of z (see (37)). We use this formula, in combination with an orthogonal transformation, to obtain the expectation of β.