Random projections for Bayesian regression
Leo N. Geppert1 , Katja Ickstadt1 ,
Alexander Munteanu2 , Jens Quedenfeld2 , and Christian Sohler2

arXiv:1504.06122v2 [stat.CO] 30 Nov 2015

1

Department of Statistics
{geppert, ickstadt}@statistik.uni-dortmund.de
2
Department of Computer Science
{alexander.munteanu, jens.quedenfeld, christian.sohler}@tu-dortmund.de
Technische Universität Dortmund
44221 Dortmund, Germany
November 30, 2015

Abstract
This article deals with random projections applied as a data reduction technique for Bayesian
regression analysis. We show sufficient conditions under which the entire d-dimensional distribution is approximately preserved under random projections by reducing the number of data
points from n to k ∈ O(poly(d/ε)) in the case n  d. Under mild assumptions, we prove that
evaluating a Gaussian likelihood function based on the projected data instead of the original
data yields a (1 + O(ε))-approximation in terms of the `2 Wasserstein distance. Our main result
shows that the posterior distribution of Bayesian linear regression is approximated up to a small
error depending on only an ε-fraction of its defining parameters. This holds when using arbitrary
Gaussian priors or the degenerate case of uniform distributions over Rd for β. Our empirical
evaluations involve different simulated settings of Bayesian linear regression. Our experiments
underline that the proposed method is able to recover the regression model up to small error
while considerably reducing the total running time.

1

Introduction

In this paper we consider linear regression. Using a linear map Π ∈ Rk×n whose choice is still to
be defined, we transform the original data set [X, Y ] ∈ Rn×(d+1) into a sketch, i.e., a substitute
data set, [ΠX, ΠY ] ∈ Rk×(d+1) that is considerably smaller. Therefore, the likelihood function can
be evaluated faster than on the original data. Moreover, we will show that the likelihood is very
similar to the original one. In the context of Bayesian regression we have the likelihood L(β|X, Y )
and additional prior information ppre (β) given in form of a prior distribution over the parameters
β ∈ Rd which we would like to estimate. Our main result is to show that the resulting posterior
distribution
ppost (β|X, Y ) ∝ f (Y |β, X) · ppre (β)
1

(1)

will also be well approximated within a small error. Please note that Bayes’ Theorem (1) contains
the probability density function f (Y |β, X). From now on, we will concentrate on its interpretation
as likelihood function L(β|X, Y ) as a function of the unknown parameter vector β as given in (2)
ppost (β|X, Y ) ∝ L(β|X, Y ) · ppre (β).

(2)

The main idea of our approach is given in the following scheme:
Π

[X, Y ] −−−−→ [ΠX, ΠY ]
↓
↓
ppost (β|X, Y ) ≈ε ppost (β|ΠX, ΠY ).
More specifically, we can reduce the number of observations from the number of input points
n to a target dimension k ∈ O(poly(d/ε)), which, in particular, is independent of n. Thus, the
running time of all subsequent calculations does not further depend on n. For instance, a Markov
Chain Monte Carlo (MCMC) sampling algorithm may be used to obtain samples from the unknown
distribution. Using the reduced data set will speed up the computations considerably. The samples
remain sufficiently accurate to resemble the original distribution and also to make statistical predictions that are nearly indistinguishable from the predictions that would have been made based on the
original sample. Note, that mathematically it is possible to achieve a similar reduction by setting
Π = X T without incurring any error. From the resulting matrix [X T X, X T Y ] it would be possible
to compute or evaluate the exact likelihood respectively posterior in some analytically tractable
cases, which is the standard textbook approach for classical linear regression and Bayesian linear regression with Gaussian prior and independent normal error, (cf. Bishop 2006; Hastie et al. 2009).
However, this approach is not numerically stable leading to ill-conditioned covariance matrices,
which is known from the literature (cf. Golub 1965; Lawson and Hanson 1995; Golub and van
Loan 2013) and evident from our experiments. For that reason, the exact calculation is not an
option in general. Instead, approximation algorithms are an alternative, where the error can be
controlled.
Using no reduction technique is also not an option, since then even likelihood evaluations depend
at least linearly on n. This does not pose a problem for small data sets. For larger n this is also
still possible, but employing reduction techniques can already be beneficial to reduce the running
time. For data sets that do not fit into the working memory, intelligent solutions are needed to
avoid frequent swapping to slower secondary memory. However, for really massive data or infinite
data streams, already the much simpler problem of computing the exact ordinary least squares
estimator in one pass over the data requires at least Ω(n) space (Clarkson and Woodruff 2009).
This makes the task impossible on finite memory machines when n grows large enough.
There are different computational models to deal with massive data sets in streaming and
distributed environments. We focus on the streaming model, formally introduced in Muthukrishnan
(2005). A data stream algorithm is given an input stream of items, like numerical values, points
in Rd or edges of a graph at a high rate. As the items arrive one by one it maintains some sort
of summary of the data that is observed so far. This can be a subsample or a linear sketch as
described above. The linearity allows for flexible dynamic updates of the sketch as we will discuss
later. At any time, the memory used by the algorithm is restricted to be sublinear, usually at most
polylogarithmic in the number of items. For geometrical problems the dependence on the dimension
2

d is often restricted to be at most a small polynomial. Also, the algorithm is allowed to make only
one single pass over the data. With these restrictions in mind, it is clear that the sketching matrix
Π ∈ Rk×n cannot be explicitly stored in memory. In fact it has to fulfill the following criteria to
allow for a streaming algorithm:
1. Π[X, Y ] approximates [X, Y ] well in the above sense.
2. Π can be stored succinctly.
3. We can efficiently generate the entries of Π.
We will see that the structured randomized constructions of Π can be provably achieved using
random bits of only limited independence. This means that the entries of Π need not be fully independent. However, if we choose a small number of entries, they behave as if they were independent.
In particular the entries can be stored implicitly, meeting the independence requirements by using
hash functions. These can be evaluated very efficiently and make the memory dependency on n
only logarithmic.
Although the techniques presented in the following are employed to make the computations
possible in a streaming setting, the results are of interest also in the non-streaming setting whenever
large data sets can be reduced to meet time and memory constraints.

2

Background and Related Work

Dimensionality reduction techniques like principal component analysis (PCA) (Jolliffe 2002) and
random projections have been widely used in Statistics and Computer Science. However, their
focus is usually on reducing the number of variables. Our method aims to reduce the number of
observations while keeping the algebraic structure of the data. This leads to a speed-up in the
subsequent (frequentist or Bayesian) regression analysis, because the running times of commonly
used algorithms heavily depend on n. Basic techniques based on PCA include principle component
regression and partial least squares (Hastie et al. 2009). More recent results using PCA stem from
the theory of core-sets for the k-means clustering problem and address the problem of computing
a small set of data that approximates the original point set with respect to the given objective
up to little, say (1 ± ε), error (Feldman et al. 2013). The concept of core-sets is related to the
early works of Madigan et al. (2002) and DuMouchel et al. (1999) on data-squashing that seeks for
data compression based on the likelihood of the observations. One of the more recent contributions
is Quiroz et al. (2015). They suggest inclusion probabilities proportional to the observations’
contribution to the likelihood, which is approximated by a Gaussian Process or a thin-plate spline
approximation. Data-squashing can lead to a considerable reduction in the necessary number of
observations, however there is a lack of approximation guarantees. These references show that in
the advent of massive data sets, besides the effort in reducing dimensionality, there is also need to
reduce the number of observations without incurring loss of too much statistical information.
Random projections have been studied in the context of low-rank approximation
(Cohen et al. 2015), least squares regression (Sarlós 2006; Clarkson and Woodruff 2009), Gaussian process regression (Banerjee et al. 2013), clustering problems (Boutsidis et al. 2010; Kerber
and Raghvendra 2014; Cohen et al. 2015), classification tasks (Paul et al. 2014) and compressed
3

sensing (Candès et al. 2006; Donoho 2006). Random projections are used similarly to our work,
to approximate a collection of subspaces, consisting only of sparse vectors (Baraniuk et al. 2007).
Also, Bayesian inference has been proposed for efficient computation in compressed sensing (Ji and
Carin 2007).
Recently there has been a series of works dealing with the statistical aspects of randomized linear
algebra algorithms. In Raskutti and Mahoney (2015) and Ma et al. (2014), the statistical properties
of subsampling approaches based on the statistical leverage scores of the data are investigated in
detail. Deviating from the worst case algorithmic perspective, it was shown in Ma et al. (2014) that
on average the leverage scores behave quite uniformly if the data is generated following a standard
linear regression model. In Yang et al. (2015), several sketching and subsampling methods are
used for fast preconditioning before solving the ordinary least squares (OLS) estimators on the
subsampled data using state of the art OLS solvers. Moreover, they give parallel and distributed
algorithms and extensive empirical evaluations on large scale data for this task. Our work, while
aware of providing worst case guarantees, continues the discussion of statistical properties to the
Bayesian setting.
Bayesian regression analysis for large scale data sets has been considered before. Guhaniyogi and
Dunson (2014) proposed reducing the number of variables via random projections as a preprocessing
step in the large d, small n scenario. They show that under several assumptions the approximation
converges to the desired posterior distribution, which is not possible in general, since it was shown
in Boutsidis and Magdon-Ismail (2014) that such dimensionality reduction oblivious to the target
variable causes additive error in the worst case.
For the large n, small d case, Balakrishnan and Madigan (2006) proposed a one-pass algorithm
that reads the data block-wise and performs a certain number of MCMC steps. When the next
block is read, the algorithm keeps or replaces some of the data points based on weights that keep
track of the importance of the data. The selection rule is justified empirically but lacks theoretical
guarantees. Theoretical support is only given in the univariate case based on central limit theorems
for Sequential Monte Carlo methods.
When allowing more passes over the data, TSQR (Demmel et al. 2012) is a QR decomposition
which works especially well in the large n, small d case and can easily be parallelized in the
MapReduce setting, as studied by Constantine and Gleich (2011) and Benson et al. (2013). TSQR
provides a numerically stable decomposition, which can be used as a preprocessing step prior to
MCMC inference which can be conducted with high accuracy. This method, however, depending on
the computational setting, has a number of limitations. It only works when the data is given rowby-row and the expensive decomposition has to be carried out a linear number of times, resulting
in a total lower bound on the running time of Ω(nd2 ) (cf. Demmel et al. 2012). The method is
restricted to `2 regression. Our method of random projections is capable of going beyond these
limitations. In particular, it can be extended to `p regression and more generally to robust M estimators (Clarkson and Woodruff 2015). This flexibility comes at the price of a loss in accuracy,
however, this loss is controllable by bounding parameters and does not lead to invalid inference.
Another approach by Bardenet et al. (2014) tries to subsample the data to approximate the
acceptance rule in each iteration of an MCMC sampler. The decision is shown to be similar to
the original with high probability in each step. The number of samples is highly dependent on

4

the variance of the logarithm of likelihood ratios. The method may be useful for interesting and
intractable cases when the variance can be bounded.
While frequentist linear regression can be solved straightforwardly by computing the projection
of the target variable to the subspace spanned by the data, Bayesian regression is typically computationally more demanding. In some cases, it is possible to calculate the posterior distribution
analytically, but in general this is not true. For that reason, an approximation of the posterior
distribution is needed. MCMC methods are one possibility and standard in Bayesian analysis.
They are reliable, but can take considerable time before they converge and sample from the desired
posterior distribution. Moreover, the running time grows with the number of observations in the
data set.
The main bottleneck of a lot of Bayesian analysis methods including MCMC is the repeated
evaluation of the likelihood. The running time of each evaluation grows linearly with the number
of observations in the data set. There are several approaches trying to reduce the computational
effort for Bayesian regression analysis by employing different algorithms that may perform more efficient in certain settings. Approximate Bayesian Computing (ABC) and Integrated Nested Laplace
Approximations (INLA) both fall into this category. The main idea behind ABC is to avoid the
exact evaluations by approximating the likelihood function using simulations (Csillery et al. 2010).
INLA (Rue et al. 2009; Martins et al. 2013) on the other hand is an approximation of the posterior
distribution that is applicable to models that fall into the class of so-called latent Gaussian models.
Both methods can lead to a considerable speed-up compared to standard MCMC algorithms.
Note however, that the speed-up is achieved by changing the algorithm which is used to conduct
the analysis. This is different in our approach, which reduces the number of observations in the data
set while approximately retaining its statistical properties. The running times of many algorithms
including MCMC algorithms highly depend on the number of observations, which means that our
proposed method also results in a speed-up of the analysis. In this article, we focus on MCMC
methods for the analysis, but in principle, as our method provably approximates the posterior, all
algorithms that assess the posterior can be employed. We did not use ABC since it is only suitable
for summary statistics of very low dimension (d < 10) (Beaumont et al. 2002; Csillery et al. 2010).
However, we have tried INLA on a small scale, achieving comparable results as with MCMC,
making the running time of the analysis independent of n. Likewise one could consider calculating
the exact formulae for analytically tractable cases of the posterior. However, we concentrate on
MCMC methods because of their general applicability and reliability.
New directions in Bayesian data analysis in the context of Big Data are surveyed in Welling
et al. (2014). Our work directly suits the criteria that is proposed in that reference for the large n
case in streaming as well as distributed computation environments.

3
3.1

Preliminaries
General notation

For the sake of a brief presentation, we introduce some notation. We denote by [n] = {1, . . . , n} the
set of all positive integers up to n. For a probability measure λ, let Eλ [X] be the expected value
of X with respect to λ. We skip the subscript in E [X] if the probability measure is clear from the
context. For a matrix M ∈ Rn×d we let M = U ΣV T denote its singular value decomposition (SVD),
5

where U ∈ Rn×d and V ∈ Rd×d are unitary matrices spanning the columnspace and rowspace of
M , respectively. Σ ∈ Rd×d is a diagonal matrix, whose elements σ1 ≥ . . . ≥ σd are the singular
values of M . We denote by σmax = σ1 the largest and by σmin = σd the smallest singular value
of M and write σi (M ) to make clear the σi belong to M . The trace of M T M equals the sum of

P
the squared singular values, i.e., tr M T M = di=1 σi2 (M ). We assume w.l.o.g. all matrices to
have full rank and stress that all our proofs carry out similarly to our presentation if the matrices
are of lower rank. One might even use knowledge about lower rank to reduce the space and time
complexities to bounds that only depend on the rank rather than on the number of variables.

3.2

Bayesian regression

A linear regression model is given in the following equation:
Y = Xβ + ξ.
Y ∈ Rn is a random variable containing the values of the response, where n is the number of
observations in the data set. X ∈ Rn×d is a matrix containing the values of the d independent variables. We denote by ξ ∼ N (0, ς 2 In ) an n-dimensional random vector that models the unobservable
error term. The dependent variable Y then follows a Gaussian distribution, Y ∼ N (Xβ, ς 2 In ). The
corresponding probability density function is


1
− 12
−n
2
2
f (y|Xβ, Σ) = (2π) |Σ| exp − 2 kXβ − yk2 ,
2ς
where Σ = ς 2 In .
In a Bayesian setting, β ∈ Rd is the unknown parameter vector which is assumed to follow an
unknown distribution p(β|X, Y ) called the posterior distribution. Prior knowledge about β can be
modeled using the prior distribution p(β). The posterior distribution is a compromise between the
prior distribution and the observed data.
In general, the posterior distribution cannot be calculated analytically. In this paper, we determine the posterior distribution employing Markov Chain Monte Carlo methods, even though the
posterior is explicitly known. Regardless of the computational problems related to these explicit
formulae (cf. Section 1), we focus on MCMC, because this work forms the basis for further research
on more complex models where analytical solutions are not obtainable. Possible extensions are hierarchical models and mixtures of normal distributions (cf. Section 7). Furthermore, we follow this
strategy to rule out possible interaction effects between sketching and MCMC that might occur
even in these basic cases. However, our empirical evaluation indicates that there are none.

3.3

Norms and metrics

Before going into details about random projections and subspace embeddings, let us first define
the matrix and vector norms used in this paper as well as the metric that we are going to use for
quantifying the distance between distributions.
Definition 1 (spectral norm). The spectral or operator norm of a matrix A ∈ Rn×d is defined as
kAk2 =

kAxk2
,
x∈Rd \{0} kxk2
sup

6

1
P
2 ) 2 denotes the Euclidean vector norm for any y ∈ Rm .
y
where kyk2 = ( m
i=1 i

A useful fact that is straightforward from Definition 1 is that the spectral norm of a matrix M
equals its largest singular value, i.e., we have kM k2 = σmax (M ) (cf. Horn and Johnson 1990).
In order to quantify the distance between probability measures and in particular between the
original posterior and its approximated counterpart we will need some further definitions. For this
sake, given two probability measures γ, ν over Rd , let Λ(γ, ν) denote the set of all joint probability
measures on Rd × Rd with marginals γ and ν, respectively.
Definition 2 (Wasserstein distance, cf. Villani (2009)). Given two probability measures γ, ν on
Rd the `2 Wasserstein distance between γ and ν is defined as

1
Z
2

1
2
W2 (γ, ν) =
inf
kx − yk2 dλ(x, y)
= inf Eλ kx − yk22 2
λ∈Λ(γ,ν) Rd ×Rd

λ∈Λ(γ,ν)

From the definition of the Wasserstein distance we can derive a measure of how much points
drawn from a given distribution will spread from the origin. The Wasserstein weight can be thought
of as a norm of a probability measure.
Definition 3 (Wasserstein weight). We define the `2 Wasserstein weight of a probability measure
γ as
Z
W2 (γ) = W2 (γ, δ) =
Rd

kxk22

1

2

dγ


1
= Eγ kxk22 2

where δ denotes the Dirac delta function.

3.4

Random projections and ε-subspace embeddings

The following definition of so called ε-subspace embeddings will be central to our work. Such an
embedding can be used to reduce the size of a given data matrix while preserving the algebraic
structure of its spanned subspace up to (1 ± ε) distortion. Before we summarize several methods
to construct a subspace embedding for a given input matrix, we give a formal definition. Here and
in the rest of the paper we assume 0 < ε ≤ 1/2.
Definition 4 (ε-subspace embedding). Given a matrix U ∈ Rn×d with orthonormal columns, an
integer k ≤ n and an approximation parameter 0 < ε ≤ 1/2, an ε-subspace embedding for U is a
map Π : Rn → Rk such that
(1 − ε)kU xk22 ≤ kΠU xk22 ≤ (1 + ε)kU xk22

(3)

holds for all x ∈ Rd , or, equivalently
kU T ΠT ΠU − Id k2 ≤ ε.

(4)

Inequality (3) is mainly used in this paper, but Inequality (4) is more instructive in the sense
that it makes clear that the embedded subspace is close to the identity, not involving much scale
or rotation.
Note that an ε-subspace embedding Π for the columnspace of a matrix M preserves its squared
singular values up to (1 ± ε) distortion, which in particular means that it also preserves its rank.
We prove this claim for completeness.
7

Observation 1. Let Π be an ε-subspace embedding for the columnspace of M ∈ Rn×d . Then
(1 − ε) σi2 (M ) ≤ σi2 (ΠM ) ≤ (1 + ε) σi2 (M )
and
(1 − 2ε) σi−2 (M ) ≤ σi−2 (ΠM ) ≤ (1 + 2ε) σi−2 (M ).
Proof. For the first claim, we make use of a min-max representation of the singular values that is
known as the Courant-Fischer theorem (cf. Horn and Johnson 1990). In the following derivation
we choose x∗ to be the maximizer of (5) and S ∗ the minimizer of (6).
σi2 (ΠM ) =
≤

min

max

S∈R(i−1)×d Sx=0,kxk2 =1
max
kΠM xk22
S ∗ x=0,kxk2 =1

kΠM xk22
(5)

= kΠM x∗ k22
≤ (1 + ε) kM x∗ k22
≤ (1 + ε)
= (1 + ε)
= (1 +

max

S ∗ x=0,kxk2 =1

min

kM xk22
max

S∈R(i−1)×d Sx=0,kxk2 =1
ε) σi2 (M ).

kM xk22

(6)

The lower bound can be derived analogously using the lower bound of (3).
Now we use the first claim to prove the second. To this end, we bound the difference
1
1
− 2
2
σi (M ) σi (ΠM )

=
=

|σi2 (M ) − σi2 (ΠM )|
εσi2 (M )
≤
σi2 (M )σi2 (ΠM )
(1 − ε) σi4 (M )
ε
σ −2 (M ) ≤ 2ε σi−2 (M ).
(1 − ε) i

There are several ways to construct an ε-subspace embedding. One of the more recent methods
is using a so called graph-sparsifier, which was initially introduced for the efficient construction of
sparse sub-graphs with good expansion properties (Batson et al. 2012). The work of Boutsidis
et al. (2013) adapted the technique to work for ordinary least-squares regression. While the initial
construction was deterministic, they also gave alternative constructions combining the deterministic decision rules with non-uniform random sampling techniques. Another approach is subspace
preserving sampling of rows from the data matrix. This technique was introduced by Drineas et al.
(2006) for `2 regression and generalized to more general subspace sampling for the p-norm (Dasgupta et al. 2009). The sampling is done proportional to the so called statistical leverage scores.
These techniques have recently been analyzed and extended in a statistical setting as opposed to the
algorithmic worst case setting (Raskutti and Mahoney 2015; Ma et al. 2014). All the aforementioned methods are in principle applicable whenever it is possible to read the input multiple times.
For instance, one needs two passes over the data to perform the subspace sampling procedures, one
for preprocessing the input matrix and another for computing the probabilities and for the actual

8

sampling. This way one might reach a stronger reduction or better statistical properties since their
(possibly random) construction depends on the input itself and therefore uses more information.
In principle our approximation results are independent of the actual method used to calculate
the embedding as long as the property given in Definition 4 is fulfilled. However, when the number
of observations grows really massive or we deal with an infinite stream, then the data can only be
read once, given time and space constraints. In order to use ε-subspace embeddings in a single-pass
streaming algorithm, we consider the approach of so called oblivious subspace embeddings in this
paper. These can be viewed as distributions over appropriately structured k × n matrices from
which we can draw a realization Π independent of the input matrix. It is then guaranteed that
for any fixed matrix U as in Definition 4 and failure probability 0 < α ≤ 1/2, Π is an ε-subspace
embedding with probability at least 1 − α. The results of our work are always conditioned on the
event that the map Π is an ε-subspace embedding omitting to further mention the error probability
α. The reader should keep in mind that there is the aforementioned possibility of failure during
the phase of sketching the data.
Instead of a subspace embedding, one might consider Π = X T a suitable choice for a sketching
matrix leading to a sketch of size d × (d + 1) from which the exact likelihood and posterior can
be characterized in some analytically tractable cases. However, it is well known that using the
resulting matrix X T [X, Y ] may result in numerical instabilities leading to bad conditioning of the
covariance matrix X T X. This effect can occur in the presence of collinearities or slight variations
from orthogonality independent of the size of the data, and may lead to highly inaccurate frequentist
estimators, (cf. Lawson and Hanson 1995). In a Bayesian setting, as we consider in this work, the
instabilities result in extremely large variances of the MCMC sample, leading to simulations that
do not converge. We will observe this behavior in Section 5.
We therefore consider the following approaches for obtaining oblivious ε-subspace embeddings:
1. The Rademacher Matrix (RAD): Π is obtained by choosing each entry independently
from {−1, 1} with equal probability. The matrix is then rescaled by √1k . This method has been
shown by Sarlós (2006) to form an ε-subspace embedding with probability at least 1 − α when
choosing essentially k = O( d log(d/α)
). This was later improved to k = O( d+log(1/α)
) in Clarkson
ε2
ε2
and Woodruff (2009), which was recently shown to be optimal by Nelson and Nguyên (2014). While
this method yields the best reduction among the different constructions that we consider in the
present work, the RAD embedding has the disadvantage that we need Θ(ndk) time to apply it to
an n × d matrix when streaming the input in general. If the input is given row by row or at least
block by block, our implementation applies a fast matrix multiplication algorithm to each block.
We remark that it is provably sufficient that the {−1, 1}-entries in each row of the RAD matrix
are basically four wise independent, i.e., when considering up to four entries of the same row, these
behave as if they were fully independent. Such random numbers can be generated using a hashing
scheme that generates BCH codes using a seed of size O(log n). This has first been noticed by
Alon et al. (1999). In our implementation we have used the four wise independent BCH scheme as
described in Rusu and Dobra (2007).
2. The Subsampled Randomized Hadamard Transform (SRHT) (originally from Ailon
and Liberty (2009)) is an embedding that is chosen to be Π = RHm D where D is an m×m diagonal
matrix where each entry is independently chosen from {−1, 1} with equal probability. The value of
m is assumed to be a power of two. It is convenient to choose the smallest such number that is not

9

smaller than n. Hm is the Hadamard-matrix of order m and R is a k×m row sampling matrix. That
is, each row of R contains exactly one 1-entry and is 0 everywhere else. The index of the 1-entry
is chosen uniformly from [m] i.i.d. for every row. The matrix is then rescaled by √1k . Since m is
often larger than n, the input data must be padded with 0-entries to compute the product ΠX. Of
course, it is not necessary to do that explicitly since all multiplications by zero can be omitted. The
target dimension needed to form an ε-subspace embedding with probability at least
1 − α using this
√ √
( d+ log n)2 log(d/α)
family of matrices was shown by Boutsidis and Gittens (2013) to be k = O(
),
ε2
which improved upon previous results from Drineas et al. (2011). Using this method, we have a
small dependency on n, which is negligible whenever n = O(exp(d)). This is often true in practice
when d is reasonably large. Compared to the RAD method, the dependency on the dimension d
is worse by essentially a factor of O(log d). It is known that k = Ω(d log d) is necessary due to the
sampling based approach. This was shown by reduction from the coupon collectors problem, i.e.,
solving one problem can be reduced to solving the other. See Halko et al. (2011) for details. The
benefit that we get is that due to the inductive structure of the Hadamard matrix, the embedding
can be applied in O(nd log k) time, which is considerably faster. It has been noticed in the original
paper (Ailon and Liberty 2009) that the construction is closely related to four wise independent
BCH codes. To our knowledge, there is no explicit proof that it is sufficient to use random bits of
little independence. However, we use again the four wise BCH scheme for the implicit construction
of the matrix D and the linear congruency generator from the standard library of C++ 11 for the
uniform subsampling matrix R. We will see in the empirical evaluation that this works well in
practice.
3. The Clarkson Woodruff (CW) sketch (Clarkson and Woodruff 2013) is the most recent
construction that we consider in this article. In this case the embedding is obtained as Π = ΦD.
The n × n matrix D is constructed in the same way as the diagonal matrix in the SRHT case.
Given a random map h : [n] → [k] such that for every i ∈ [n] its image is chosen to be h(i) = t ∈ [k]
with probability k1 , again Φ is a binary matrix whose 1-entries can be defined by Φh(i),i = 1. All
other entries are 0. This is obviously the fastest embedding, due to its sparse construction. It
can be applied to any matrix X ∈ Rn×d in O(nnz(X)) = O(nd) time, where nnz(X) denotes the
number of non-zero entries in X. This is referred to as input sparsity time and is clearly optimal up
to small constants, since this is the time needed to actually read the input from a data stream or
external memory, which dominates the sketching phase. However, its disadvantage is that the target
dimension is k = Ω(d2 ) (Nelson and Nguyen 2013b). Roughly spoken, this is necessary due to the
need to obliviously and perfectly hash d of the standard basis vectors spanning Rn . Improved upper
2
bounds over the original ones of Clarkson and Woodruff (2013) show that k = O( εd2 α ) is sufficient
to draw an ε-subspace embedding from this distribution of matrices with probability at least 1 − α
(Nelson and Nguyen 2013a). This reference also shows that it is sufficient to use only four wise
independent random bits to generate the diagonal matrix D. Again, in our implementation we
use the four wise independent BCH scheme from Rusu and Dobra (2007). Moreover, Φ can be
constructed using only pairwise independent entries. This can be achieved very efficiently using the
fast universal hashing scheme introduced by Dietzfelbinger et al. (1997) which we have employed
in our implementation. The space requirement is only O(log n) for a hash function from this class.
For a really fast implementation using bit-wise operations the actual size parameters of the sketch
are chosen to be the smallest powers of two that are larger than the required sizes n and k.

10

sketching method
RAD
SRHT
CW

target dimension


O d+log(1/α)
2
ε


d·log(d/α)
O
2
ε
 2
d
O ε2 α

running time
O (ndk)
O (nd log k)
O (nnz(X)) = O(nd)

Table 1: Comparison of the three considered ε-subspace embeddings; nnz(X) denotes the number of non-zero entries
in X, α denotes the failure probability.

Table 1 summarizes the above discussion, in particular the trade-off behavior between time and
space complexity of the presented sketching methods. While in general one is interested in the
fastest possible application time, memory constraints might make it impossible to apply the CW
sketch due to its quadratic dependency on d. Taking it the other way, for a fixed sketching size,
CW will give the weakest approximation guarantee (cf. Yang et al. 2015). For really large d, even
the O(d log d) factor of SRHT might be too large so that we have to use the slowest RAD sketching
method.

3.5

Extension to the streaming model

The presented reduction techniques are of interest whenever we deal with medium to large sized
data for reducing time and space requirements. However, when the data grows massive, we need to
put more importance on the computational requirements. We therefore want to briefly discuss and
give references to some of these technical details. For example, while the dimensions of the resulting
sketches do not depend on n, this is not true for the embedding matrices Π ∈ Rk×n . However,
due to the structured constructions that we have surveyed above, we stress that the sketching
matrices can be stored implicitly by using the different hash functions of limited independence.
The hash functions used in our implementations are the four wise independent BCH scheme used
in the seminal work of Alon et al. (1999) and the universal hashing scheme by Dietzfelbinger et al.
(1997). These can be evaluated very efficiently using bit-wise operations and can be stored using a
seed whose size is only O(log n). Note that even this small dependency on n is only needed in the
sketching phase. After the sketch has been computed, the space requirements will be independent
of n. A survey and evaluation of alternative hashing schemes can be found in Rusu and Dobra
(2007).
The linearity of the embeddings allows for efficient application in sequential streaming and in distributed environments, see e.g. Clarkson and Woodruff (2009); Woodruff and Zhang (2013); Kannan
et al. (2014). The sketches can be updated in the most flexible dynamic setting, which is commonly
referred to as the turnstile model (Muthukrishnan 2005). In this model, think of an initial matrix
of all zero values. The stream consists of updates of the form (i, j, u) meaning that the entry Xij will
be updated to Xij + u. A single entry can be defined by one single update or by a sequence of not
necessarily consecutive updates. For example a stream S = {. . . , (i, j, +5), . . . , (i, j, −3), . . .} will
result in Xij = 2. Even deletions are possible in this setting by using negative updates. Clearly this
also allows for additive updates of rows or columns, each consisting of consecutive single updates
to all the entries in the same row or column. At first sight this model might seem very technical

11

and unnatural. But the usual form of storing data in a table is not appropriate or performant for
massive data sets. The data is rather stored as a sequence of (key, value) pairs in arbitrary order.
For dealing with such unstructured data, the design of algorithms working in the turnstile model
is of high importance.
For distributed computations, note that the embedding matrices can be communicated efficiently to every machine in a computing cluster of l machines. This is due to the small implicit
P
representation by hash functions. Now, suppose the data is given as X = li=1 X (i) where X (i)
is stored on the machine with index i ∈ [l]. Note that by the above data representation in form
of updates, X (i) can consist of rows, columns or single entries of X. Again, multiple updates to
the same entry are possible and may be distributed to different machines. Every machine i ∈ [l]
can compute a small sketch on its own share X (i) of the data and efficiently communicate it to
one dedicated central server. A sketch of the entire data set can be obtained by summing up the
P
single sketches since ΠX = li=1 ΠX (i) . More details can be found in Kannan et al. (2014). Recent
implementations of similar distributed and parallel approaches for OLS are given by Yang et al.
(2015).
The above discussions make clear that our methods suit the criteria that need to be satisfied
when dealing with Big Data (cf. Welling et al. 2014). Namely, the number of data items that need
to be accessed at a time is only a small subset of the whole data set, particularly independent of
n. The algorithms should be amenable to distributed computing environments like MapReduce.

4

Theory

In this section we introduce and develop the theoretical foundations of our approach and will
combine them with existing results on ordinary least squares regression to bound the Wasserstein
distance between the original likelihood function and its counterpart that is defined only on the
considerably smaller sketch. Empirical evaluations supporting and complementing our theoretical
results will be conducted in the subsequent section.

4.1

Embedding the likelihood

The following observation is standard (cf. Givens and Shortt (1984); Kannan and Vempala (2009))
and will be helpful in bounding the `2 Wasserstein distance of two Gaussian measures. It allows us
to derive such a bound by inspecting their means and their covariances separately.
Observation 2. Let Z1 , Z2 ∈ Rd be random variables with finite first moments m1 , m2 < ∞ and
let Z1m = Z1 − m1 , respectively, Z2m = Z2 − m2 be their mean-centered counterparts. Then it holds
that




E kZ1 − Z2 k22 = km1 − m2 k22 + E kZ1m − Z2m k22 .
Proof.




E kZ1 − Z2 k22 = E kZ1m − Z2m + m1 − m2 k22


= E kZ1m − Z2m k22 + km1 − m2 k22 + 2 (m1 − m2 )T E [Z1m − Z2m ]
|
{z
}


=0
= E kZ1m − Z2m k22 + km1 − m2 k22

12

In our first lemma we show that using an ε-subspace embedding Π for the columnspace of
[X, Y ], we can approximate the least squares regression problem up to a factor of 1 + ε. That is,
we can find a solution ν by projecting ΠY into the columnspace of ΠX such that kXν − Y k2 ≤
(1 + ε) minβ∈Rd kXβ − Y k2 . Similar proofs can be found in Clarkson et al. (2013) and Boutsidis
et al. (2013). We repeat the result here for completeness.
Lemma 5. Given X ∈ Rn×d , Y ∈ Rn , let Π be an (ε/3)-subspace embedding for the columnspace of
[X, Y ]. Now let γ = argminβ∈Rd kXβ − Y k22 and similarly define ν = argminβ∈Rd kΠ(Xβ − Y )k22 .
Then
kXν − Y k22 ≤ (1 + ε) kXγ − Y k22 .
Proof. Let [X, Y ] = U ΣV T denote the SVD of [X, Y ]. Now define η1 = ΣV T [γ T , −1]T and η2 =
ΣV T [ν T , −1]T . Using this notation we can rewrite U η1 = Xγ − Y and similarly U η2 = Xν − Y .
We have that
(1 − ε/3) kU η2 k22 ≤ kΠU η2 k22 ≤ kΠU η1 k22 ≤ (1 + ε/3) kU η1 k22 .
The first and the last inequality are direct applications of the subspace embedding property (3),
whereas the middle inequality follows from the optimality of ν in the embedded subspace. Now,
by rearranging and resubstituting terms, this yields


1 + ε/3
2
kXν − Y k2 ≤
kXγ − Y k22 ≤ (1 + ε) kXγ − Y k22 .
1 − ε/3

√
√
One can even show that a distortion of order ε, i.e., an O( ε)-subspace embedding is already
enough to get the result. This was shown by using a more complicated proof taking the geometry
of the least squares solution into account and using the property that the solution is obtained
by an orthogonal projection onto the columnspace spanned by the data matrix (cf. Clarkson and
Woodruff 2009). Putting it the other way around, by using an (ε/3)-subspace embedding as in
Lemma 5, we even have
kXν − Y k22 ≤ (1 + ε2 ) kXγ − Y k22 .
(7)
In the following, we investigate the distributions proportional to the likelihood functions p ∝
L(β|X, Y ) and p0 ∝ L(β|ΠX, ΠY ) and bound their Wasserstein distance.
We begin our contribution with a bound on the distance of their means γ and ν, respectively.
We generalize upon previous results for the specific embedding methods to arbitrary ε-subspace
embeddings.
Lemma 6. Given X ∈ Rn×d , Y ∈ Rn , let Π be an (ε/3)-subspace embedding for the columnspace of
[X, Y ]. Now let γ = argminβ∈Rd kXβ − Y k22 and similarly define ν = argminβ∈Rd kΠ(Xβ − Y )k22 .
Then
ε2
kXγ − Y k22 .
kγ − νk22 ≤ 2
σmin (X)

13

Proof. Let X = U ΣV T denote the SVD of X. Let η = V T (γ − ν). First note that γ and ν are both
contained in the columnspace of V (cf. Sarlós 2006) which means that V T is a proper rotation
with respect to γ − ν. Thus,
kX(γ − ν)k22 = kU ΣV T (γ − ν)k22 = kΣV T (γ − ν)k22
X
X
2
=
σi2 ηi2 ≥
σmin
ηi2
2
2
= σmin
kV T (γ − ν)k22 = σmin
kγ − νk22 .

Consequently, it remains to bound kX(γ − ν)k22 . This can be done by using the fact that the
minimizer γ is obtained by projecting Y orthogonally onto the columnspace of X. Therefore, we
have X T (Xγ − Y ) = 0 (cf. Clarkson and Woodruff 2009). Furthermore, by Equation (7) it holds
that kXν − Y k22 ≤ (1 + ε2 )kXγ − Y k22 . Now by plugging this into the Pythagorean theorem and
rearranging we get that
kX(γ − ν)k22 = kXν − Y k22 − kXγ − Y k22 ≤ ε2 kXγ − Y k22 .
Putting all together this yields the proposition
kγ − νk22 ≤

ε2
1
2
2
kX(γ
−
ν)k
≤
2
2 (X)
2 (X) kXγ − Y k2 .
σmin
σmin

Now that we have derived a bound on the distance of the different means, recall that by
Observation 2, we can assume w.l.o.g. γ = ν = 0 when we consider the variances. Namely,


it remains to derive a bound on inf E kZ1m − Z2m k22 , i.e., the least expected squared Euclidean
distance of two points drawn from a joint distribution whose marginals are the mean-centered
original distribution and its embedded counterpart. Of course we can bound this quantity by
explicitly defining a properly chosen joint distribution and bounding the expected squared distance
for its particular choice. This is the idea that yields our next lemma.
Lemma 7. Let p ∝ L(β|X, Y ) and p0 ∝ L(β|ΠX, ΠY ). Let Z1m , Z2m be the mean-centered versions
of the random variables Z1 ∼ p and Z2 ∼ p0 that are distributed according to p and p0 respectively.
Then we have



inf Eλ kZ1m − Z2m k22 ≤ ε2 tr (X T X)−1 .
λ∈Λ(p,p0 )

Proof. Our plan is to design a joint distribution that deterministically maps points from one distribution to another in such a way that we can bound the distance of every pair of points. This
can be done by utilizing the Dirac delta function δ(·), which is a degenerate probability density
function that concentrates all probability mass at zero and has zero density otherwise. Given a
bijection g : Rd → Rd we can define such a joint distribution λ ∈ Λ(p, p0 ) through its conditional
distributions λ(x | y) = δ(x − g(y)) for every y ∈ Rd . It therefore remains to define g.
According to Observation 1, when applying the embedding Π, the columnspace of a matrix is
expanded or contracted, respectively, by a factor of at most (1 ± ε). We will make use of this
fact in the following way. Let X = U ΣV T and ΠX = Ũ Σ̃Ṽ T denote the SVDs of X and ΠX,
respectively. Now, to define the x-y-pairs that will be mapped to each other by g, we consider
14

vectors x, x0 , y, y 0 ∈ Rd where x0 and y 0 are contained in the columnspaces of V and Ṽ , respectively.
To obtain the bijection g, let the vectors have the following properties for arbitrary, but fixed radius
c ≥ 0:
1. kx0 k2 = ky 0 k2 = c
2. x = ΣV T x0
3. y = Σ̃Ṽ T y 0
4. ∃τ > 0 : x = τ y.
Observe that by the first property x0 and y 0 lie on a d-dimensional sphere with radius c centered at
0. Therefore, there exists a rotation matrix R ∈ Rd×d such that y 0 = Rx0 . Note that such a map is
bijective by definition. The second item defines a map of such spheres to ellipsoids (also centered
at 0) given by ΣV T . Recall that x0 was chosen from the columnspace of V . Thus, this map can be
seen as bijection between the d-dimensional vector space and a d-dimensional subspace contained
in n-dimensional space. The third property is defined analogously. The fourth property urges that
x and y both lie on a ray emanating from 0. Note that any such ray intersects each ellipsoid exactly
once.
Our bijection can be defined accordingly as
g : Rd → Rd
x

7→ Σ̃Ṽ T RV Σ−1 x

by composing the map ΣV T , defined in the second item, with the rotation R and finally with Σ̃Ṽ T
from the third property. The map is bijective since it is obtained as the composition of bijections.
Now, in order to bound the distance kZ1m − Z2m k22 for any realization of (Z1m , Z2m ) according to
their joint distribution defined above, we can derive a bound on the parameter τ . Substituting the
second and third properties into the fourth, we get that
ΣV T x0 = τ Σ̃Ṽ T y 0
which can be rearranged to
y 0T y 0 τ

= (y 0T Ṽ )Σ̃−1 Σ(V T x0 )
X
σi
=
(y 0T Ṽ )i (V T x0 )i
σ̃i
X
σi
≤
(y 0T Ṽ )i (V T x0 )i √
σi 1 − ε
X
0T
≤ (1 + ε)
(y Ṽ )i (V T x0 )i
≤ (1 + ε) c2 .

√
The first inequality follows from σ̃i ≥ 1 − ε σi and the second from the assumption ε ≤ 1/2. This
eventually means that τ ≤ (1 + ε) since y 0T y 0 = c2 by the first property.
√
A lower bound of τ ≥ (1 − ε) can be derived analogously by using σ̃i ≤ 1 + ε σi .

15

Now we can conclude our proof. It follows that






Eλ0 kZ1m − Z2m k22 ≤ Eλ kZ1m − Z2m k22 ≤ Eλ kεZ1m k22
λ0 ∈Λ(p,p0 )



= ε2 Eλ kZ1m k22 = ε2 tr (X T X)−1 .
inf

The last equality holds since the expected squared norm of the mean-centered random variable is
just the trace of its covariance matrix.
Combining the above results we get the following lemma.
Lemma 8. Let Π be an (ε/3)-subspace embedding for the columnspace of X. Let p ∝ L(β|X, Y )
and p0 ∝ L(β|ΠX, ΠY ). Then
W22 (p, p0 ) ≤


ε2
2
2
T
−1
kXγ
−
Y
k
+
ε
tr
(X
X)
.
2
2 (X)
σmin

Proof. The lemma follows from Definition 2, Observation 2, Lemma 6 and Lemma 7.
Under mild assumptions we can argue that this leads to a (1 + O(ε))-approximation of the
likelihood with respect to the Wasserstein weight (see Definition 3).
Corollary 9. Let Π be an (ε/3)-subspace embedding for the columnspace of X. Let p ∝ L(β|X, Y )
and similarly let p0 ∝ L(β|ΠX, ΠY ). Let κ(X) = σmax (X)/σmin (X) be the condition number of X.
Assume that for some ρ ∈ (0, 1] we have kXγk2 ≥ ρkY k2 . Then


κ(X)
0
W2 (p ) ≤ 1 +
ε W2 (p).
ρ
Proof. By definition, the squared `2 Wasserstein weight of p equals its second moment. Since p is
a Gaussian measure with mean γ and covariance matrix (X T X)−1 , we thus have
W22 (p) = kγk22 + tr (X T X)−1



and similarly

W22 (p0 ) = kνk22 + tr (X T ΠT ΠX)−1 .
Since Π is an (ε/3)-subspace embedding for the columnspace of X we know from Observation
1, that all the squared singular values of X are approximated up to less than (1 ± ε) error and so
are their inverses. Therefore, we have that


tr (X T ΠT ΠX)−1 ≤ (1 + ε) tr (X T X)−1 .

(8)

It remains to bound kνk22 . To this end we use the assumption that for some ρ ∈ (0, 1] we have
kXγk2 ≥ ρkY k2 . By X T (Xγ − Y ) = 0 and applying the Pythagorean Theorem this means that


1
kXγk22
2
2
2
2
kXγ − Y k2 = kY k2 − kXγk2 ≤ kXγk2
−
1
≤
.
ρ2
ρ2

16

Now we can apply the triangle inequality, Lemma 6, Inequality (9) and Definition 1 to get
kνk2 ≤ kγk2 + kν − γk2
ε
≤ kγk2 +
kXγ − Y k2
σmin (X)
ε
≤ kγk2 +
kXγk2
ρσmin (X)
ε
≤ kγk2 +
kXk2 kγk2
ρσmin (X)
ε
= kγk2 + κ(X)kγk2
ρ


κ(X)
=
1+
ε kγk2 .
ρ
Combining this with Inequality (8), the claim follows since
(1 +

κ(X) 2
ρ ε)

κ(X)
ρ

≥ 1 and therefore (1 + ε) ≤

and finally taking square roots on both sides.

We stress that the assumption that there exists some constant ρ ∈ (0, 1] such that kXγk2 ≥
ρkY k2 is very natural and mild in the setting of linear regression since it means that at least a
constant fraction of the dependent variable Y can be explained within the columnspace of the data
X (cf. Drineas et al. 2006). If this is not true, then a linear model is not appropriate at all for the
given data.

4.2

Bayesian regression

So far we have shown that using subspace embeddings to compress a given data set for regression
yields a good approximation to the likelihood. Note that in a Bayesian regression setting Lemma
8 already implies a similar approximation error for the posterior distribution if the priors for
β are chosen to be uniform distributions over R. This is an improper, non-informative choice,
ppre (β) = 1Rd . From this, it follows that
ppost (β|X, Y ) ∝ L(β|X, Y ) · 1Rd
= L(β|X, Y ).
The remaining term is just the Gaussian likelihood which is proper. For regression models, especially on data sets with large n, this covers a considerable amount of the cases of interest (cf.
Gelman et al. 2014). We will extend this to arbitrary Gaussian priors ppre (β) leading to our main
result: an approximation guarantee for Gaussian Bayesian regression in its most general form.
To this end, let m be the mean of the prior distribution and let S be derived from its covariance
matrix by Σ = ς 2 (S T S)−1 . Now, the posterior distribution is given by
ppost (β|X, Y ) ∝ L(β|X, Y ) · ppre (β)




1
1
1
1
2
2
=
· exp − 2 kXβ − Y k2 ·
− 2 kS(β − m)k2 .
d
1 · exp
2ς
2ς
(2πς 2 )n/2
(2π) 2 |Σ| 2
Thus, we know that up to some constants that are independent of β, the exponent of the posterior
can be described by
kXβ − Y k22 + kS(β − m)k22
(9)
17

which contains all the information to define the mean and covariance structure of the posterior
distribution. Now let
" #
"
#
X
Y
Z=
and z =
.
S
Sm
With these definitions we can rewrite Equation (9) above as kZβ − zk22 . This, in turn, can be
treated as a (frequentist) regression problem to which we can apply Lemma 8. We just have to
use a subspace embedding for the columnspace of [Z, z] instead of only embedding [X, Y ]. We will
see that it is not necessary to do this explicitly. More precisely, embedding only the data matrix
is sufficient to have a subspace embedding for the entire columnspace defined by the data and the
prior information and, therefore, to have a proper approximation of the posterior distribution. This
is formalized in the following lemma.
Lemma 10. Let M = [M1T , M2T ]T ∈ R(n1 +n2 )×d be an arbitrary matrix. Suppose Π is an ε-subspace
embedding for the columnspace of M1 . Let In2 ∈ R(n2 ×n2 ) be the identity matrix. Then
"
#
Π 0
P =
∈ R(k+n2 )×(n1 +n2 )
0 In2
is an ε-subspace embedding for the columnspace of M .
Proof. Fix an arbitrary x ∈ Rd . We have
|kP M xk22 − kM xk22 | = |kΠM1 xk22 + kM2 xk22 − kM1 xk22 − kM2 xk22 |
= |kΠM1 xk22 − kM1 xk22 |
≤ εkM1 xk22
≤ ε(kM1 xk22 + kM2 xk22 )
= εkM xk22
which concludes the proof by singular value decomposition M = U ΣV T and surjectivity of the
linear map ΣV T .
This lemma finally enables us to prove our main theoretical result.
Theorem 11. Let Π be an (ε/3)-subspace embedding for the columnspace of X. Let ppre (β) be an
arbitrary Gaussian distribution with mean m and covariance matrix Σ = ς 2 (S T S)−1 . Let
" #
"
#
X
Y
Z=
and z =
.
S
Sm
Let µ = argminβ∈Rd kZβ − zk2 be the posterior mean. Let p ∝ L(β|X, Y ) · ppre (β) and p0 ∝
L(β|ΠX, ΠY ) · ppre (β). Then
W22 (p, p0 ) ≤


ε2
kZµ − zk22 + ε2 tr (Z T Z)−1 .
2
σmin (Z)

18

Proof. From our previous reasoning we know that approximating the posterior distribution can
be reduced to approximating a likelihood function that is defined in terms of the data as well as
the parameters of the prior distribution. This has been shown by rewriting Equation (9) above as
kZβ − zk22 . For that reason, we can apply Lemma 8 to get the desired result if we are given an
(ε/3)-subspace embedding for the columnspace of Z. Using Lemma 10 we know that for this, it
is sufficient to use an (ε/3)-subspace embedding for the columnspace of [X, Y ] independent of the
covariance and mean that define the prior distribution.
Similar to Corollary 9 we have the following result concerning the posterior distribution.
Corollary 12. Let Π be an (ε/3)-subspace embedding for the columnspace of X. Let ppre (β) be an
arbitrary Gaussian distribution with mean m and covariance matrix Σ = ς 2 (S T S)−1 . Let
" #
"
#
X
Y
Z=
and z =
.
S
Sm
Let µ = argminβ∈Rd kZβ − zk2 be the posterior mean. Let p ∝ L(β|X, Y ) · ppre (β) and p0 ∝
L(β|ΠX, ΠY ) · ppre (β). Let κ(Z) be the condition number of Z. Assume that for some ρ ∈ (0, 1]
we have kZµk2 ≥ ρkzk2 . Then we have


κ(Z)
0
W2 (p ) ≤ 1 +
ε W2 (p).
ρ
Both Theorem 11 and Corollary 12 show that the sketch preserves the expected value and the
covariance structure of the posterior distribution very well. Note that for normal distributions,
these parameters fully characterize the distribution as they are sufficient statistics. Therefore, one
can see the corresponding parameters based on the sketched data set as very accurate approximate
sufficient statistics for the posterior distribution.

5

Simulation Study

To validate the proposed method empirically, we conduct a simulation study. For this, we employ
MCMC methods to obtain the posterior distributions for the parameters of the Bayesian regressions. Please note that the sketching techniques can also be combined with other methods. We
concentrate on MCMC methods, however, as they are very reliable, widely-used and allow for easy
checking of convergence. The different sketching methods were implemented as described in Section
3.4 and technically more detailed in the given references. All codes were written in the C++ 11
programming language and compiled using GCC 4.7. For fast matrix multiplications we employed
the LAPACK 3.5.0 library where applicable. Our R-package RaProR uses the above implementations. The package is available on its project website (Geppert et al. 2015). All simulations
were done using R, version 3.1.2 (R Core Team 2014) and the R-package rstan, version 2.3 (Stan
Development Team 2013). The simulations were conducted on an Intel Xeon E5430 quad-core
CPU running at 2.66 GHz using 16 GB DDR2 memory on a Debian GNU Linux 7.8 distribution.
The hard drive used was a Seagate Momentus 7200.4 G-Force 500 GB, SATA 3Gb/s HDD with
7200 rpm and 16 MB cache.

19

5.1

Data generation

For the simulation study, we create a set of data sets. We vary the number of observations n,
the number of variables d, and the standard deviation of the error term ς. The variation of n
is introduced to monitor whether the running time of the analyses based on sketches is indeed
independent of n and also to see how well the proposed method deals with growing data sets. We
choose values of n ∈ {50 000, 100 000, 500 000, 1 000 000}. The size of the sketches depends mainly
on the number of variables in the data set. For this reason, we conduct simulations with two values
of d, d = 50 and d = 100. The reason for choosing rather small values of n and d is that our aim
is to compare the results of the MCMC on the sketched data sets to the results on the respective
original data set. The sketching methods presented here can handle larger values of d and arbitrary
values of n, but employing MCMC on the original data set then becomes unfeasible. The standard
deviation of the error term is important, because the goodness of the approximation also depends
on the goodness of the model (cf. Lemma 8 and Theorem 11). Here, we choose ς ∈ {1, 2, 5, 10},
thus ranging from very well-fitting models to models with quite high error variance.
The generated true values of β follow a zero-inflated Poisson distribution, where the expected
value of the Poisson distribution is 3 and the probability of a component exhibiting an excess
zero is equal to 0.5. This means that the components of β have no influence, i.e. are 0, with
probability 0.5 and follow a Poi(3) distribution with probability 0.5. All components that follow a
Poi(3) distribution are multiplied with (−1) with probability 0.5. The data set X is obtained in
two steps. At first, a d-dimensional vector that represents the column means is drawn randomly
from a N (0, 25) distribution. In a second step, the actual values of X are drawn from a normal
distribution with the column mean as expected value and variance of four. The variance in the
columns of X is thus lower than the error variance for two of our choices and the same or less for
the other two choices. Y is then generated by multiplying X with β and adding the error term, in
accordance with the model.

5.2

Regression model

We employ a standard Bayesian linear regression model (cf. Section 3.2)
Y ∼ N (Xβ, ς 2 In )
with independent uniform priors over R for all components of β, which are improper, non-informative
prior distributions. For ς, the uniform prior is limited to the positive part of R. We choose an improper uniform prior rather than an inverse gamma prior with small values for the hyperparameters
as Gelman (2006) indicates that such priors can have a skewing effect on the posterior distribution.
When using improper prior distributions, it is necessary to ensure that the posterior distributions
are proper. For our choice of uniform distributions, this does not pose a problem, since the uniform
priors are represented by indicators, which are constant over R, ppre (β) = 1Rd , and for that reason
do not influence the likelihood. Thus, it follows that
ppost (β|X, Y ) ∝ L(β|X, Y ) · 1Rd
= L(β|X, Y ).

20

2000

●

1000

● ●

run−time in minutes

500

●

● ●

200

●
● ●

100

●

50

●

●

●

20
●

10
5

●
●

2

●
●

50

100

200

500

1000

2000

5000

10000

20000

50000

number of observations

Figure 1: Running times for simulated data sets with varying number of observations and d = 50 variables

The remaining term is just the Gaussian likelihood which is proper with respect to both, β and ς
(cf. Gelman et al. 2014). Although closed-form expressions are known for this model, we employ
MCMC for the reasons motivated in Sections 2 and 3.2.
Our theoretical guarantees comprise the posterior distributions of β. Our model is thus more
general than the theoretical results. As an alternative, ς can be fixed to an estimated value obtained
√
using kΠXβ − ΠY k2 / n. The results we present in the following are all based on β.

5.3

Preliminary simulations

In a first step, we consider the running times necessary to carry out Bayesian regression for (nonembedded) data sets with an increasing number of observations n, employing the No-U-TurnSampler (Hoffman and Gelman 2014), which is implemented in the R-package rstan. We use
standard settings and four chains, which are sampled in parallel. The resulting running times are
plotted in Figure 1. The running time depends at least linearly on n, with occasional jumps that
are probably due to swaps to external memory, which slows down the computation by a large factor.
There is an outlier for the case of n = 50, which takes a lot more time to compute than n = 100
and even more than n = 200. Since we have d = 50 variables, the total number of parameters (52)
is higher than the number of observations for this case only. While the linear dependency on the
number of observations does not pose a problem for small to medium-sized data sets, for big data
settings MCMC methods become unfeasible. This underlines the usefulness of embedded data sets.
Before we consider our embedding methods, we pick up on the idea of using Π = X T as sketching
matrix. If the data set is given row by row, the sketch ΠX = X T X can be computed efficiently
P T
using the tensor product of the row vectors X T X =
x x , which is a sum of d × d matrices
Pi Ti
with rank one. Analoguously, we have ΠY = X T Y =
xi Yi . This results in a sketch of size
d × (d + 1), which is the smallest possible sketch for problems of full rank and has no error at all in
the sense that the exact likelihood respectively posterior distribution can be calculated from these
matrices. We have argued before that this method is not numerically stable in general. As we
focus on MCMC in this work, we show how this effect influences the run of the MCMC sampler.
We tried analyzing Bayesian linear regression models based on X T [X, Y ], using some of the data
21

5.0e+09
−5.0e+09
−1.5e+10

β41

0

2000

4000

6000

8000

10000

Iteration

Figure 2: Traceplot of MCMC sample (with four chains) for one parameter based on data set obtained using Π = X T
as sketching matrix. Original data set contains n = 10 000 observations and d = 50 variables

sets described in Section 5.1 and also a smaller data set with n = 10 000, which was generated in
the same way. We have found that the models do not converge in practice. Increasing the number
of iterations does not seem to be a remedy as the variance of the MCMC sample grows with more
iterations. Figure 2 shows an exemplary traceplot for one parameter, consisting of four chains. The
range of the sample is enormous. The variation is high for all of the chains, they exhibit standard
deviations in the order of 109 . When reducing the number of iterations from the 10 000 in Figure
2 to, say, 5 000, the range of the MCMC sample decreases markedly. However, there is no sign of
convergence with minimum and maximum around −4 · 109 and 4 · 109 , respectively.
We can deal with both of these issues using subspace embeddings as we can underline in our
next experiments. We conduct a series of simulations that aims at comparing our proposed method
to the standard method on the original data. To obtain the subspace embeddings, we employ the
three approaches described in Section 3. As approximation parameter, we choose ε = 0.1 and
ε = 0.2 for all three methods. We do not recommend using values of ε > 0.2. Table 2 contains
the number of observations of the sketches depending on the number of variables and the value of
d
the approximation parameters. The sizes for RAD and SRHT are both set to k = d d log
e to be
ε2
comparable. They differ by one due to rounding errors. For the CW sketch we used k equal to the
d2
smallest power of two larger than 20ε
2 . Please note that the CW embeddings generally result in a
higher number of observations due to the quadratic dependency on d. However, the opposite is true
for 50 variables. This is due to the constants that we used. The constants are set to 1 in the case
of RAD and SRHT. For the RAD method this was empirically evaluated by Venkatasubramanian
and Wang (2011). For CW the constant may be much smaller as indicated by a lower bound in
1
Nelson and Nguyen (2013b). Preliminary experiments on small scale led to our choice of 20
.

5.4

Comparison of posterior means

To evaluate the results, we first compare the posterior means of the models based on the embedded
data sets with the posterior means of the model based on the original data set. Table 3 contains
an overview of the sum of squared distances between the embedded data sets’ posterior means and

22

d

ε

RAD

SRHT

CW

50
50

0.1
0.2

20 546
5 136

20 547
5 137

16 384
4 096

100
100

0.1
0.2

47 174
11 793

47 175
11 794

65 536
16 384

Table 2: Number of observations of the sketches for different values of d and ε

n

sketch

ε

ς=1

ς=2

ς=5

ς = 10

5 · 104
5 · 104
5 · 104
5 · 104
5 · 104
5 · 104
1 · 105
1 · 105
1 · 105
1 · 105
1 · 105
1 · 105

RAD
RAD
SRHT
SRHT
CW
CW
RAD
RAD
SRHT
SRHT
CW
CW

0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2

0.052
0.014
0.001
0.004
0.025
0.016

0.025
0.781
0.009
0.077
0.004
0.040

0.021
0.892
0.021
0.093
0.021
0.156
0.836
0.061
0.025
0.171
0.056
2.624

0.834
1.512
0.165
0.757
0.195
0.915
0.958
0.777
0.964
0.617
3.844
2.937

Table 3: Sum of squared distances between posterior mean values of the original model and models based on the
respective sketches

those of the original model. Geometrically, this is the squared Euclidean distance of the posterior
mean vectors. As indicated by Theorem 11, the sum of squared distances grows with the standard
deviation of the error term. There does not seem to be a systematic difference in performance
between the different sketching methods. With larger ε, we usually, but not necessarily observe an
increase in the distance. Please note that some values are missing, because the original models did
not converge within reasonable time bounds.
In addition to the comparison to the original models’ mean, we also compare the posterior
means to the true means. Table 4 contains the sum of the squared distances between the true
mean for d = 50 and varying values of ς. The general picture looks very similar to the results in
Table 3. The original model often exhibits the smallest sum of squared distances, but sometimes
models based on embedded data sets are closer to the true mean. Again, there does not seem to be
a systematic difference between the sketching methods. The squared distances do not seem to be
influenced by the value of n, with some squared distances even exhibiting smaller values for larger
n.

23

n

sketch

5 · 104
5 · 104
5 · 104
5 · 104
5 · 104
5 · 104
5 · 104
1 · 105
1 · 105
1 · 105
1 · 105
1 · 105
1 · 105
1 · 105
5 · 105
5 · 105
5 · 105
5 · 105
5 · 105
5 · 105
1 · 106
1 · 106
1 · 106
1 · 106
1 · 106
1 · 106

none
RAD
RAD
SRHT
SRHT
CW
CW
none
RAD
RAD
SRHT
SRHT
CW
CW
RAD
RAD
SRHT
SRHT
CW
CW
RAD
RAD
SRHT
SRHT
CW
CW

ε

ς=1

ς=2

ς=5

ς = 10

0.1
0.2
0.1
0.2
0.1
0.2

0.000
0.048
0.012
0.001
0.005
0.022
0.014

0.003
0.016
0.710
0.018
0.059
0.011
0.056

0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2

0.007
0.033
0.030
0.007
0.004
0.011
0.009
0.045
0.009
0.016
0.027
0.050
0.001
0.080
0.002
0.002
0.002
0.003

0.031
0.009
0.136
0.125
0.232
0.072
0.223
0.322
0.147
0.033
0.097
0.009
0.016
0.011
0.010
0.183
0.289
0.047

0.065
0.124
0.506
0.032
0.046
0.046
0.089
0.065
1.354
0.117
0.040
0.387
0.022
3.484
0.563
1.729
0.418
0.085
1.305
0.135
0.126
0.072
0.599
2.029
1.202
0.100

4.614
1.718
10.845
3.372
8.721
6.474
1.870
0.035
0.679
0.579
0.696
0.453
4.496
3.473
12.920
0.658
0.059
2.978
0.153
3.579
3.967
1.357
0.288
4.329
4.445
0.395

Table 4: Sum of squared distances between true mean values and posterior means of models based on the respective
sketches

5.5

Comparison of fitted values

After this comparison on the level of parameters – whose number is not changed by sketching – we
will compare the models on the level of observations, of which the sketches contain merely a fraction
of the number of observations in the original data set. We multiply X with the posterior mean
vector of β, where this posterior mean can be based on the original data set or on the respective
sketches. In a frequentist sense, these are fitted values Ŷ , but all X values are taken from the
original data set, not necessarily from the data set the model is based on. This is done to see how
close the approximation is on the level of Y values for both ε = 0.1 and ε = 0.2. Figure 3 is a
scatterplot which contains smoothed densities. The fitted values based on the original model are on
the x-axis while the fitted values based on the CW sketch (with ε = 0.1) are on the y-axis. Darker
shades of black stand for more observations. Even though the fitted values are based on one of the

24

50
0
−50
−100

based on embedded model using CW (ε = 0.1)

−150
−200
−200

−150

−100

−50

0

50

based on model using original data set

Figure 3: Comparison of fitted values based on the original data set with n = 50 000, d = 50, ς = 10 and a CW
sketch with ε = 0.1. Darker shades of black stand for more observations
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

4
2

●
●
●
●
●
●
●

●
●
●
●
●
●

●
●
●
●
●
●

●
●
●
●
●

●
●
●
●
●
●

●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●

0

CW (ε = 0.1)

SRHT (ε = 0.1)

RAD (ε = 0.1)

−4

●
●
●
●
●
●
●
●
●
●
●

CW (ε = 0.2)

●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●

SRHT (ε = 0.2)

●
●
●
●

RAD (ε = 0.2)

−2

Figure 4: Difference of fitted values according to models based on the respective sketching methods and fitted values
according to model based on original data set with n = 50 000, d = 50, ς = 10

data sets with the highest standard deviation of the error (n = 50 000, ς = 10), all values are close or
reasonably close to the bisecting line. This means that the fitted values obtained by the two models
do not differ by much. To get a better overview, Figure 4 depicts the distances between the fitted
values as boxplots. Here, all three sketching methods with both ε = 0.1 and ε = 0.2 are included.
All six sets of distances are centered around zero. The effect of the approximation parameter ε is
evident from the boxplot, the variation is larger for ε = 0.2 regardless of the sketching method.
When fixing ε, all three sketching methods exhibit similar results, although the RAD sketch seems
to introduce slightly more variation into the differences than the other two sketching methods, thus
prediction for the data used in learning the model is highly accurate.
We have also generated additional data that has not been used in learning the model and
employed the posterior mean to predict y-values for these data. The results are very similar to
those described above and presented in Figures 3 and 4. Sketching, again, introduces a little more

25

β33

●

●
●

●
●
●

●
●
●

●
●
●
●
●
●

CW (ε = 0.2)

●

●
●

SRHT (ε = 0.2)

●
●
●

RAD (ε = 0.2)

SRHT (ε = 0.2)

RAD (ε = 0.2)

CW (ε = 0.1)

SRHT (ε = 0.1)

RAD (ε = 0.1)

original model

●
●

●
●
●

CW (ε = 0.1)

●
●
●
●
●
●
●

●
●

●
●
●

●
●

original model

●
●

●
●
●
●

−0.1

●
●

CW (ε = 0.2)

●
●

0.0

●
●
●
●

SRHT (ε = 0.1)

●
●

●
●
●
●
●

●
●
●

RAD (ε = 0.1)

●

0.1

●
●

●
●
●

●
●
●
●
●

0.2

●
●

−4.15

−4.05

−3.95

−3.85

β11

Figure 5: Boxplots of MCMC sample for two parameters based on data set with n = 50 000, d = 50, ς = 5 and
respective sketches

variation, depending on ε. More formal treatment of prediction accuracy based on similar sketching
techniques for the OLS solution is given in Raskutti and Mahoney (2015).

5.6

Comparison of posterior distributions

As we have conducted Bayesian regression the model consists not only of a mean value, but of
a whole posterior distribution for each parameter. Figure 5 contains two exemplary boxplots
of MCMC samples representing marginal posterior distributions. The original data set contains
n = 50 000 observations, d = 50 variables and has an error standard deviation of ς = 5. The
medians of the MCMC samples based on the original data set are well-represented by the MCMC
samples based on sketches. Even though the median based on an embedded data set might be
higher or lower for certain parameters, we did not find any systematic biases. The embedding
introduces additional variation, which depends on the value of the approximation parameter ε, but
does not seem to be influenced by the choice of sketching method.
In regression, a common task is the identification of important variables by means of variable
selection. In a Bayesian setting, this can be done via credible intervals, among other approaches.
Our results indicate that the identification of important variables is quite accurately possible based
on the resulting approximate models. However, one has to take the additional variation into
account. Exemplarily, when using 95% credible intervals as criterion, one should not compare
the endpoints of the credible interval to a fixed value µ. Instead, take the extra variation in the
posterior distribution and also possible small shifts of the mean and median into account. We
therefore recommend using smaller values of ε when aiming at variable selection (see Figure 5).

5.7

Comparison of running times

One of our aims is to make Bayesian regression feasible on very large data sets. After ensuring
that the results are close to those obtained by the analysis on the original data set, we will now
contemplate the running time required. The total running time is composed of the time required
for the analysis and the time required for the necessary preliminary steps: reading the data set into
26

n

sketch

5 · 104
5 · 104
5 · 104
5 · 104
5 · 104
5 · 104
5 · 104
1 · 105
1 · 105
1 · 105
1 · 105
1 · 105
1 · 105
1 · 105
5 · 105
5 · 105
5 · 105
5 · 105
5 · 105
5 · 105
5 · 105
1 · 106
1 · 106
1 · 106
1 · 106
1 · 106
1 · 106
1 · 106

none
RAD
RAD
SRHT
SRHT
CW
CW
none
RAD
RAD
SRHT
SRHT
CW
CW
none
RAD
RAD
SRHT
SRHT
CW
CW
none
RAD
RAD
SRHT
SRHT
CW
CW

ε
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2
0.1
0.2

ς=1
0.32
1.60
0.40
0.03
0.02
0.01
0.01
0.69
3.27
0.76
0.05
0.04
0.02
0.02
5.49
16.88
3.73
0.20
0.19
0.09
0.09
18.23
51.77
7.92
0.41
0.39
0.19
0.21

Preprocessing
ς = 2 ς = 5 ς = 10
0.41
1.68
0.39
0.02
0.02
0.01
0.01
0.83
3.41
0.84
0.06
0.04
0.03
0.02
5.16
15.96
4.00
0.21
0.18
0.09
0.08
12.88
147.42
8.46
0.49
1.44
0.27
0.19

0.43
1.68
0.42
0.03
0.02
0.01
0.01
1.02
3.41
0.84
0.05
0.05
0.02
0.02
5.92
16.10
4.03
0.20
0.19
0.09
0.09
12.59
33.75
8.38
0.61
0.68
0.38
0.45

0.44
1.73
0.43
0.03
0.02
0.01
0.01
1.05
3.27
0.80
0.05
0.04
0.02
0.02
5.71
16.36
3.85
0.21
0.19
0.10
0.09
14.09
34.71
8.21
0.62
0.68
0.46
0.39

ς=1

Analysis
ς=2
ς=5

1095.55
315.12
23.17
317.34
29.04
375.22
26.92

749.10
213.42
26.00
278.39
30.65
293.89
25.77

278.87
21.44
284.96
23.72
257.50
21.94

ς = 10

260.80
23.21
282.20
26.82
278.20
26.29

616.47
156.79
17.48
181.94
23.26
164.56
20.57
2035.81
167.24
17.52
128.60
21.52
186.95
21.22

498.68
154.73
21.81
166.41
26.00
171.82
22.94
1617.28
182.92
23.65
196.48
22.70
198.78
23.45

279.81
27.37
310.20
31.37
335.74
26.03

313.33
27.19
308.01
25.62
300.32
25.23

165.85
17.22
190.78
22.26
189.33
24.39

198.40
19.58
190.18
24.76
166.92
22.86

209.19
21.27
341.12
26.61
281.72
28.58

279.03
19.93
264.99
31.32
232.40
19.50

215.78
22.87
294.04
19.69
175.49
22.05

145.64
23.43
154.77
23.69
144.02
9.72

Table 5: Running times for data sets with d = 50. Columns 4 to 7 (“Preprocessing”) contain the running times of
the sketching methods in minutes, for the original data set, the values represent the time required to read the data
set into memory, which is a prerequisite for every sketching method. The four columns to the right (“Analysis”)
contain the running times of the Bayesian linear regression analysis in minutes

27

2000
1500
1000
500

running times in minutes

0
n
sketch

5 ⋅ 104

1 ⋅ 105

original

1 ⋅ 105

5 ⋅ 104
RAD

5 ⋅ 104

1 ⋅ 105

SRHT

1 ⋅ 105

5 ⋅ 104
CW

Figure 6: Total running times in minutes for data sets with n ∈ {50 000, 100 000}, d = 50, ς = 5 and approximation
parameter ε = 0.1. For the sketched data sets, the total running time consists of the time for reading, sketching and
analyzing the data set. For the original data set, the sketching time is 0 since this step is not applied

memory and calculating the sketch. Table 5 contains the running times required for the Bayesian
regression analysis (four rightmost columns headed “Analysis”) and the running times for the
preliminary steps (columns headed “Preprocessing”). For the original data sets, the preliminary
steps consist of the time required for reading the data set into memory, for all other cases, Table 5
gives the running times required for constructing the sketch. The total running time of an analysis
on a sketch is obtained by summing up the reading time, the sketching time, and the time required
for the Bayesian linear regression analysis.
Table 5 suggests that both the reading times and the sketching times are stable for data sets
of the same size, with the possible exception of an outlier for n = 1 000 000 and ς = 2. Both the
reading times and the sketching times grow with the size of the data set. The sketching times grow
for smaller values of ε. For RAD sketches with ε = 0.1, the sketching takes longer than the reading
of the data set, for all other combinations the opposite is true. CW sketches require the shortest
amount of running time of the three sketching methods.
Although only a few of the original data sets could be analyzed, the “Analysis” values in Table 5
indicate that the running times for the Bayesian analysis increase with the number of observations
in the data set. The running times for the sketched data sets on the other hand show no systematic
increase for growing values of n. There is some variation in the running times, but this seems to be
more random chance than trend. Larger values of ε lead to shorter running times in the analysis,
which indicates that the trade-off between time and goodness of the approximation is present both
in calculating and analyzing the sketch. The running time of the analyses is similar for all three
sketching methods. However, the running times do seem to depend on the value of ς. For higher
standard deviations of the error term, the required running time tends to become less.
Figure 6 exemplarily shows the total running times in minutes for data sets with d = 50 and
ς = 5. This comprises reading, sketching (if applicable) and analyzing the data set. For the
sketches, ε = 0.1 was used. To illustrate the potential improvement with increasing n, Figure 6
contains the total running times for n = 50 000 and n = 100 000 side by side. For the original data

28

sets, we observe a large increase in running time as n doubles, whereas the running times on the
sketched data sets hardly change.
All results presented so far are based on data sets with d = 50 variables. We have conducted
Bayesian analyses on data sets with d = 100 with the same parameter values for n and ς. The
findings from these simulations are similar to those for d = 50. One exception is that the analysis of
CW sketches now takes longer than the analysis of the respective original data sets for n = 50 000.
This changes for larger data sets. One strength of the CW method, however, is its efficiency
when dealing with streaming data when the number of variables is not too large. Recall that its
dependency on d is quadratic, which means that the sketch sizes are already very large for medium
dimensional problems with d ≈ 500 or d ≈ 1 000. Setting the target dimension to a lower value is
not a remedy, because this results in a weaker approximation guarantee as also noticed by Yang
et al. (2015). In that case one should consider using one of the slower and denser sketches.

5.8

Streaming example and concluding remarks

The data sets considered so far were chosen to be small enough to allow for Bayesian analysis on
the full data set in a convenient time. This is by far not what might be called Big Data. To
show that our approach is suitable for analyzing really big amounts of data, we have generated
and at the same time sketched a data set. The generation of the data followed the same rules as
the other simulated data sets, but with ς = 0.1. The data set’s original size is 109 × 100 double
precision values. This corresponds to about 2 TB in CSV format or at least 750 GB in a binary
representation. We have used the CW sketching method resulting in a sketch of size 65 536 × 100,
which requires only around 140 MB of space in CSV format and fits into RAM. Clearly, we cannot
compare the results to those on the original data set, but calculating the sum of squared distances
between the true mean values of β and the posterior mean of the model adds up to 3.741 · 10−6 .
The Bayesian regression analysis took 2781 minutes.
In some cases the algebraic structure might strongly depend on a few observations or variables.
In such situations it is important to identify these or to retain their contribution in the reduced data
set. So far, our model assumptions did not suffer from such ill-behaved situations, but now we assess
the performance of our method in this case. We construct data sets where an important part of the
target variable falls into a subspace that is spanned by a small constant number of observations.
Uniform random subsampling will pick these only with probability O( n1 ). Oblivious subsampling
techniques in general will have trouble identifying the important observations. In contrast, oblivious
subspace embedding techniques preserve these effects with high probability. This effect is observed
in practice even when comparing one sketch against the best of 1000 repetitions of uniform random
subsampling.
In conclusion, the simulation study indicates that our proposed method works well for simulated
data sets, which are generally well-suited for conducting Bayesian linear regression. But even with
a high variance of the error term (and thus a relatively bad model fit), our proposal leads to results
similar to those one would obtain on the original data set. The running time of the analysis with
the proposed sketches is largely independent of n, giving advantages for very large n. Since the
embeddings can be read in sequentially, it is not necessary to load the whole data set into the
memory at once, which reduces the required memory.

29

Variable

Description

Remark

cnt
season
yr
hour
holiday
weekday
weathersit
atemp
hum
windspeed

number of rental bikes used
season of the year
year (2011 or 2012)
hour (0 to 23)
public holiday
day of the week
weather (“clear” to “rain”)
apparent temperature
humidity
windspeed

Y -variable
factor (4 levels)
factor (2 levels)
factor (24 levels)
factor (2 levels)
factor (7 levels)
factor (3 levels)
standardized
standardized
standardized

Table 6: Variables from the bike sharing data set used in the model

For CW embeddings, reading the data in and calculating the sketch only takes marginally longer
than only reading the data in. In our experiments, we found that reading in and sketching takes
around 1.01 to 1.04 times longer. This factor is typically higher for small data sets and lower for
larger data sets (cf. Table 5). However, when the number of variables is large it may be favorable
to use SRHT, whose sketching time is only slightly larger but has considerably smaller embedding
dimension.

6

Real Data Example

As a real data example, we consider the bike sharing data set taken from Fanaee-T and Gama
(2014), which is available in the UCI Machine Learning Repository (Lichman 2013). This is
only meant as an exemplary application of the methods to a real data scenario and should not
be mistaken for a complete statistical analysis of the data set. The bike sharing data set contains
the number of rental bike users per hour over two years as well as additional information about
the day and the weather. See Table 6 for an overview of the variables we use in the model. The
data set contains some additional variables we do not employ, because they are highly correlated
with the variables present in the model. We also made a change to the factor levels of the variable
weathersit. In the original data set, this variable contains 4 levels. The fourth level only is present
3 times out of total of n = 17 379 hours in the data set. To avoid any problems with such an
underrepresented level, we combine levels 3 and 4 to obtain a factor with 3 levels. The original
levels 3 and 4 stand for “light rain” and “heavy rain”. The new level 3 can easily be interpreted as
the presence of rain. For a more detailed description of all variables, please refer to the data set’s
web page on the UCI Machine Learning Repository. 1
The variable cnt contains the number of rental bikes used per hour and is thus a count variable.
However, there are around 850 distinct values, which makes analyzing cnt as a continuous variable
reasonable. When analyzing such variables with a linear regression model, transforming them using
1

http://archive.ics.uci.edu/ml/datasets/Bike+Sharing+Dataset

30

d

ε

RAD

SRHT

CW

40
40

0.15
0.20

6 767
3 807

6 767
3 807

8 192
4 096

Table 7: Number of observations of the sketches for the bike sharing example. Different values of ε are used for
RAD and SRHT sketches; the target dimension of CW sketches is chosen to be the power of two closest to the size
of the RAD and SRHT sketches

the square-root is a common procedure. After the transformation, the values of cnt show some
bi-modality, but fit reasonably well to the assumption of a normal distribution.
We use the transformed variable cnt as Y -variable and all other variables in Table 6 as Xvariables. To handle the factor variables, we use the R-function model.matrix to create a design
matrix, which is then passed on to RaProR and rstan. The resulting design matrix contains
n = 17 379 observations and d = 39 variables plus the intercept. Again, we calculate sketches for
all three methods and with two different settings of ε. Because of the size of the data set relative
to the number of variables, we choose ε = 0.15 and ε = 0.2 for the RAD and SRHT sketches. For
CW, we choose values of k that are closest to the target dimension of the other sketches. This
results in the values given in Table 7.
The Bayesian model based on the original data set suggests that all mentioned variables are
important for the modeling of the number of bikes used per hour. Figure A.1 in the appendix
gives an overview of the posterior distributions for the model based on the original data set. As
one might expect, the weather has a strong influence. More bikes are rented when the apparent
temperature is high and, to a lesser extent, when the humidity and the wind speed are low. In
clear or partly cloudy weather, the number of rented bikes is highest, but the negative influence
of heavier clouds is comparatively small. Rainy weather, however, reduces the number of bike
users more substantially. In addition to that, fall seems to be the most popular seasons for bike
sharing. Spring and summer also have positive effects in comparison to winter, but the effect sizes
are smaller. This might seem surprising at first, especially as the number of rental bike users is
highest in summer. This might partly be an effect of the apparent temperature, which is generally
higher in summer.
There is also a distinct hourly effect. During night time, especially between midnight and 5am,
the number of rented bikes is greatly reduced. On the other hand, between 7am and 9pm, a lot of
bikes are used, with two peaks at 8am and 5pm/6pm. This might indicate that the service is used
by people transiting to and from work. Holidays – which only includes days that would otherwise
be a working day, so only Monday to Friday – have a negative influence on the number of rental
bike users. When taking the days of the week into account, Friday and Saturday have the highest
positive effect while Sunday seems to be the least popular day. All of these effects based on days
have a small influence compared to the variables mentioned before. Lastly, the variable yr also has
a positive influence which indicates a positive trend for this bike sharing service.
All conclusions can similarly be derived from the models based on the sketches. Following our
approach in Section 5, we first compare the resulting posterior mean values of β. Table 8 shows the
sum of squared distances between posterior mean values of the original model and models based
on the three sketching methods, using two values of ε each. There is a general increase in the sum
31

ε

RAD

SRHT

CW

0.15
0.2

1.790
6.511

2.349
2.732

0.907
1.657

Table 8: Sum of squared distances between posterior mean values of the original model and models based on the
respective sketches for the bike sharing data set
1.5
1.0

●
●
●
●
●
●
●
●
●
●

0.5

●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●

0.0
●
●
●
●
●
●

●
●
●
●
●
●
●
●

−1.0

RAD (ε = 0.2)

CW (ε = 0.15)

SRHT (ε = 0.15)

RAD (ε = 0.15)

●
●

●
●
●
●
●
●

●
●
●
●
●
●

CW (ε = 0.2)

●
●
●
●

SRHT (ε = 0.2)

−0.5

Figure 7: Difference of fitted values according to models based on the respective sketching methods and fitted values
according to model based on the original bike sharing data set

of squared distances for ε = 0.2 compared to ε = 0.15, but the amount differs depending on the
sketching method. This should not be over-interpreted, however. As these values represent only
one realization of a random subspace embedding, there is no evidence for systematic differences.
To see the effects of the differences in the posterior means of β on the level of the y-variable,
which is the number of rental bikes used, we compare the fitted values as in Section 5.5. Again,
we multiply the original design matrix X with the posterior means of β, where the posterior mean
values are obtained from the model on the original design matrix and the models on the respective
sketches. Figure 7 contains the six resulting boxplots. As in the simulation study, the differences
of the fitted values are centered around zero with only small deviations. Further analysis indicates
that the higher deviations occur when the number of bikes used is high, which means that the
majority of differences in the fitted values are small relative to the observed value for that data
point.
As a last step of our short analysis of the bike sharing data set, we will concentrate on the
posterior distributions of the factors that take the weather situation into account. This factor has
three levels in our data set. The first level stands for clear weather, which also includes partly
cloudy hours, the second level stands for heavier clouds or mist while the third level models rainy
weather, which includes light rain, heavy rain, thunderstorms, and snow. The different levels occur
with different relative frequencies: around 66% of the observed hours fall into level 1, 26% into
level 2 and 8% into level 3.
Figure 8 shows boxplots of the MCMC samples based on the original design matrix and the
sketches. The values represent the marginal posterior distributions of the two dummy variables
32

rain vs. clear

●
●
●

●
●

●
●
●
●
●

●
●
●
●

●

●
●
●
●
●

RAD (ε = 0.15)

original model

CW (ε = 0.2)

SRHT (ε = 0.2)

RAD (ε = 0.2)

CW (ε = 0.15)

SRHT (ε = 0.15)

RAD (ε = 0.15)

●
●
●
●

original model

●
●
●
●

●
●
●
●
●

●
●
●
●

●
●
●

●
●
●

●
●
●
●

CW (ε = 0.2)

●
●
●

●
●
●
●
●

SRHT (ε = 0.2)

●
●

●
●
●
●

●
●

RAD (ε = 0.2)

●
●

●
●
●
●

●
●
●
●

CW (ε = 0.15)

●
●
●

●
●
●
●
●

●
●

SRHT (ε = 0.15)

●
●
●

−2.5

●
●
●
●

●
●
●
●

−3.0

−0.5

0.0

0.5

−2.0

heavy clouds vs. clear

Figure 8: Boxplots of MCMC sample for the two weather situation parameters based on the original data set and
all sketches

associated with the variable weathersit. The scales of the boxplots are chosen such that one unit
is of the same length on both y-axes. This allows for easy comparison of the variation in the
two posterior distributions. Again, we can see that the embedding introduces additional variation,
which depends on the value of the approximation parameter ε, but does not seem to be influenced
by the choice of the sketching method. In addition, the posterior distributions for the factor “heavy
clouds” show less variation compared to the posterior distributions for “rain”. This is as one might
expect as the number of occurrences is more than three times higher for “heavy clouds” and a larger
number of observations tends to reduce the uncertainty. Nonetheless, it is interesting to observe
that the variation introduced by the embedding seems to be a factor of the variation present in the
original model.
While the effect of the factor “rain” is undoubtedly negative according to the original model
and all sketches used here, the effect of factor “heavy clouds” is close to zero. In the original model,
“heavy clouds” would also be seen as an influential factor when using the 95% credible interval as a
criterion. The conclusion is the same for all sketching methods when using ε = 0.15. However, when
using ε = 0.2 and CW, “heavy clouds” would be seen as not influential. This stresses again that the
endpoints of credible intervals based on sketches exhibit some additional variation and inference
based on them may change, depending on the variation in the original model and the choice of ε.
If variable selection is a focus of the regression analysis, we recommend choosing reasonably small
ε (cf. Figure 5 and its discussion in Section 5.6).
This example underlines that our method also works well on real world applications when the
original data follows the model assumptions reasonably well.

7

Conclusion

Our paper deals with random projections as a data reduction technique for Bayesian regression. We
have shown how projections can be applied to compress the columnspace of a given data matrix
with only little distortion. The size of the reduced data set is independent of the number n of
observations in the original data set. Therefore, subsequent computations can operate within time
33

and space bounds that are also independent of n, regardless of which algorithm is actually used.
While our focus was on MCMC and the No-U-Turn-Sampler in particular, we tried INLA as well and
observed a considerable reduction in running time while achieving very similar results. However,
our proposed reduction method is not limited to these approaches, making it highly flexible.
The presented embedding techniques allow for fast application to the data set and do not need
the embedding matrices to be stored explicitly. Thus, only very little memory is needed while
sketching the data. Furthermore, we have surveyed their useful properties when the computations
are performed in sequential streaming as well as in distributed environments. These scenarios are
highly desirable when dealing with Big Data (cf. Welling et al. 2014).
We consider the situation where the likelihood is modeled using standard linear regression with a
Gaussian error term. We show that the likelihood is approximated within small error. Furthermore,
if an arbitrary Gaussian distribution is used as prior distribution, the desired posterior distribution
is also well approximated within small error. This includes the case of a uniform prior distribution
over Rd , an improper, non-informative choice that is widely used (cf. Gelman et al. 2014). We
also show the results to be (1 + O(ε)) approximations of the distributions of interest in the context
of Bayesian linear regression. As the structure of both mean and variance is preserved up to small
errors, approximate sufficient statistics for the posterior distributions are well-recovered. This gives
the user all the information that is needed for Bayesian regression analysis.
In our simulation experiments, we found that the approximation works well for simulated data
sets. All three sketching methods we considered lead to results that are very similar to Bayesian
regression on the full data set and the true underlying values. The running time for the MCMC
analysis based on the sketches is independent of the number of observations n. The calculation of
the embedding does depend on n, but requires little more time than the necessity of only reading
the data. This is especially true when using the CW method. But for larger dimensions a CW
embedding can be too large. In such a case, the denser SRHT construction also performs very well
and is preferable because of its lower dependency on d. RAD has even lower dependency on d but
takes considerably more time to calculate.
We have applied the methods to a bike sharing data set by Fanaee-T and Gama (2014). The
approximation also works well on this real data example, giving very similar results to the original
Bayesian regression, while adding little additional variation to the posterior distributions.
For future research, we would like to generalize our results to other classes of distributions
for the likelihood and to more general priors. As a first step, we have used hierarchical models
involving normal, Gamma and exponential distributions as hyperpriors. For normal and Gamma
distributions, the results seem promising, whereas using exponential distributions seems more challenging. The recent results on frequentist `p regression of Woodruff and Zhang (2013) might give
rise to efficient streaming algorithms also in the Bayesian regression setting. Another interesting
direction would be to consider Gaussian mixtures, since they allow to approximate any continuous
distribution.
In real-world applications one might exhibit the domain-specific structure to further reduce the
time and space bounds when these indicate that the data itself is of low rank or allows for sparse
solution vectors.

34

Acknowledgements We would like two thank the anonymous referees for their valuable comments and suggestions which helped to improve this manuscript. This work has been supported
by Deutsche Forschungsgemeinschaft (DFG) within the Collaborative Research Center SFB 876
“Providing Information by Resource-Constrained Analysis”, project C4.

References
Ailon, N. and Liberty, E. (2009). Fast dimension reduction using Rademacher series on dual BCH
codes. Discrete Comput. Geom., 42(4):615–630.
Alon, N., Matias, Y., and Szegedy, M. (1999). The space complexity of approximating the frequency
moments. J. Comput. Syst. Sci., 58(1):137–147.
Balakrishnan, S. and Madigan, D. (2006). A one-pass sequential Monte Carlo method for Bayesian
analysis of massive datasets. Bayesian Anal., 1(2):345–361.
Banerjee, A., Dunson, D. B., and Tokdar, S. T. (2013). Efficient Gaussian process regression for
large datasets. Biometrika, 100(1):75–89.
Baraniuk, R., Davenport, M., Devore, R., and Wakin, M. (2007). A simple proof of the restricted
isometry property for random matrices. Constr. Approx., 28(3).
Bardenet, R., Doucet, A., and Holmes, C. C. (2014). Towards scaling up Markov Chain Monte
Carlo: an adaptive subsampling approach. In Proc. of ICML, pages 405–413.
Batson, J. D., Spielman, D. A., and Srivastava, N. (2012). Twice-Ramanujan sparsifiers. SIAM J.
Comput., 41(6):1704–1721.
Beaumont, M. A., Zhang, W., and Balding, D. J. (2002). Approximate Bayesian Computation in
population genetics. Genetics, 162(4):2025–2035.
Benson, A. R., Gleich, D. F., and Demmel, J. (2013). Direct QR factorizations for tall-and-skinny
matrices in MapReduce architectures. In Proc. of IEEE Int. Conf. on Big Data, pages 264–272.
Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer.
Boutsidis, C., Drineas, P., and Magdon-Ismail, M. (2013). Near-optimal coresets for least-squares
regression. IEEE Trans. Inf. Theory, 59(10):6880–6892.
Boutsidis, C. and Gittens, A. (2013). Improved matrix algorithms via the Subsampled Randomized
Hadamard Transform. SIAM J. Matrix Anal. Appl., 34(3):1301–1340.
Boutsidis, C. and Magdon-Ismail, M. (2014). Faster SVD-truncated regularized least-squares. In
IEEE Symp. Inf. Theory, pages 1321–1325.
Boutsidis, C., Zouzias, A., and Drineas, P. (2010). Random projections for k-means clustering. In
Proc. of NIPS, pages 298–306.

35

Candès, E. J., Romberg, J. K., and Tao, T. (2006). Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Trans. Inf. Theory,
52(2):489–509.
Clarkson, K. L., Drineas, P., Magdon-Ismail, M., Mahoney, M. W., Meng, X., and Woodruff, D. P.
(2013). The fast Cauchy transform and faster robust linear regression. In Proc. of SODA, pages
466–477.
Clarkson, K. L. and Woodruff, D. P. (2009). Numerical linear algebra in the streaming model. In
Proc. of STOC, pages 205–214.
Clarkson, K. L. and Woodruff, D. P. (2013). Low rank approximation and regression in input
sparsity time. In Proc. of STOC, pages 81–90.
Clarkson, K. L. and Woodruff, D. P. (2015). Sketching for M -estimators: A unified approach to
robust regression. In Proc. of SODA, pages 921–939.
Cohen, M. B., Elder, S., Musco, C., Musco, C., and Persu, M. (2015). Dimensionality reduction
for k-means clustering and low rank approximation. In Proc. of STOC.
Constantine, P. G. and Gleich, D. F. (2011). Tall and skinny QR factorizations in MapReduce
architectures. In Proc. of Int. Workshop on MapReduce and Appl., pages 43–50. ACM.
Csillery, K., Blum, M., Gaggiotti, O., and Francois, O. (2010). Approximate Bayesian Computation
(ABC) in practice. Trends Ecol. Evol., 25(7):410–418.
Dasgupta, A., Drineas, P., Harb, B., Kumar, R., and Mahoney, M. W. (2009). Sampling algorithms
and coresets for `p regression. SIAM J. Comput., 38(5):2060–2078.
Demmel, J., Grigori, L., Hoemmen, M., and Langou, J. (2012). Communication-optimal parallel
and sequential QR and LU factorizations. SIAM J. Sci. Comput., 34(1).
Dietzfelbinger, M., Hagerup, T., Katajainen, J., and Penttonen, M. (1997). A reliable randomized
algorithm for the closest-pair problem. J. Algorithms, 25(1):19–51.
Donoho, D. L. (2006). Compressed sensing. IEEE Trans. Inf. Theory, 52(2):1289–1306.
Drineas, P., Mahoney, M. W., and Muthukrishnan, S. (2006). Sampling algorithms for `2 regression
and applications. In Proc. of SODA, pages 1127–1136.
Drineas, P., Mahoney, M. W., Muthukrishnan, S., and Sarlós, T. (2011). Faster least squares
approximation. Numer. Math., 117(2):219–249.
DuMouchel, W., Volinsky, C., Johnson, T., Cortes, C., and Pregibon, D. (1999). Squashing flat
files flatter. In Proc. of KDD, pages 6–15.
Fanaee-T, H. and Gama, J. (2014). Event labeling combining ensemble detectors and background
knowledge. Prog. in AI, 2(2-3):113–127.
Feldman, D., Schmidt, M., and Sohler, C. (2013). Turning big data into tiny data: Constant-size
coresets for k -means, PCA and projective clustering. In Proc. of SODA, pages 1434–1453.
36

Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models (comment
on article by Browne and Draper). Bayesian Anal., 1(3):515–534.
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2014).
Bayesian Data Analysis. Texts in Statistical Science. Chapman & Hall/CRC, 3rd edition.
Geppert, L. N., Ickstadt, K., Munteanu, A., Quedenfeld, J., and Sohler, C. (2015). RaProR:
Random Projections for Bayesian linear Regression, R-package, Version 1.0.
Givens, C. R. and Shortt, R. M. (1984). A class of Wasserstein metrics for probability distributions.
Michigan Math. J., 31(2):231–240.
Golub, G. H. (1965). Numerical methods for solving linear least squares problems. Numerische
Mathematik, 7(3):206–216.
Golub, G. H. and van Loan, C. F. (2013). Matrix computations. Johns Hopkins University Press,
4th edition.
Guhaniyogi, R. and Dunson, D. B. (2014). Bayesian compressed regression. J. Amer. Stat. Assoc.,
published online.
Halko, N., Martinsson, P.-G., and Tropp, J. A. (2011). Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Rev., 53(2):217–
288.
Hastie, T., Tibshirani, R., and Friedman, J. (2009). The Elements of Statistical Learning: Data
mining, Inference, and Prediction. Springer, 2nd edition.
Hoffman, M. D. and Gelman, A. (2014). The No-U-Turn sampler: Adaptively setting path lengths
in Hamiltonian Monte Carlo. J. Mach. Learn. Res., 15:1593–1623.
Horn, R. and Johnson, C. (1990). Matrix Analysis. Cambridge University Press.
Ji, S. and Carin, L. (2007). Bayesian compressive sensing and projection optimization. In Proc. of
ICML, pages 377–384.
Jolliffe, I. (2002). Principal component analysis. Springer, 2nd edition.
Kannan, R. and Vempala, S. (2009). Spectral algorithms. Found. Trends Theor. Comput. Sci.,
4(3-4):157–288.
Kannan, R., Vempala, S., and Woodruff, D. P. (2014). Principal component analysis and higher
correlations for distributed data. In Proc. of COLT, pages 1040–1057.
Kerber, M. and Raghvendra, S. (2014). Approximation and streaming algorithms for projective
clustering via random projections. CoRR, abs/1407.2063.
Lawson, C. L. and Hanson, R. J. (1995). Solving least squares problems. Classics in applied
mathematics. SIAM, Philadelphia (Pa.).
Lichman, M. (2013). UCI machine learning repository.
37

Ma, P., Mahoney, M. W., and Yu, B. (2014). A statistical perspective on algorithmic leveraging.
In Proc. of ICML, pages 91–99.
Madigan, D., Raghavan, N., DuMouchel, W., Nason, M., Posse, C., and Ridgeway, G. (2002).
Likelihood-based data squashing: A modeling approach to instance construction. Data Min.
Knowl. Discov., 6(2):173–190.
Martins, T., Simpson, D., Lindgren, F., and Rue, H. (2013). Bayesian computing with INLA: New
features. Comput. Stat. Data Anal., 67:68–83.
Muthukrishnan, S. (2005). Data streams: Algorithms and applications. Found. Trends Theor.
Comput. Sci., 1(2).
Nelson, J. and Nguyen, H. L. (2013a). OSNAP: Faster numerical linear algebra algorithms via
sparser subspace embeddings. In Proc. of FOCS, pages 117–126.
Nelson, J. and Nguyen, H. L. (2013b). Sparsity lower bounds for dimensionality reducing maps. In
Proc. of STOC, pages 101–110.
Nelson, J. and Nguyên, H. L. (2014). Lower bounds for oblivious subspace embeddings. In Proc.
of ICALP, Part I, pages 883–894.
Paul, S., Boutsidis, C., Magdon-Ismail, M., and Drineas, P. (2014). Random projections for linear
support vector machines. ACM Trans. Knowl. Discov. Data, 8(4):22.
Quiroz, M., Villani, M., and Kohn, R. (2015). Speeding up MCMC by efficient data subsampling.
stat.ME, abs/1404.4178.
R Core Team (2014). R: A Language and Environment for Statistical Computing. R Foundation
for Statistical Computing, Vienna, Austria.
Raskutti, G. and Mahoney, M. (2015). Statistical and algorithmic perspectives on randomized
sketching for ordinary least-squares. In Proc. of ICML, pages 617–625.
Rue, H., Martino, S., and Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian
models by using integrated nested Laplace approximations (with discussion). J. R. Statist. Soc.
B, 71:319–392.
Rusu, F. and Dobra, A. (2007). Pseudo-random number generation for sketch-based estimations.
ACM Trans. Database Syst., 32(2):11.
Sarlós, T. (2006). Improved approximation algorithms for large matrices via random projections.
In Proc. of FOCS, pages 143–152.
Stan Development Team (2013). Stan: A C++ Library for Probability and Sampling, Version 2.3.
Venkatasubramanian, S. and Wang, Q. (2011). The Johnson-Lindenstrauss transform: An empirical
study. In Proc. of ALENEX, pages 164–173.
Villani, C. (2009). Optimal transport: Old and new. Grundlehren der mathematischen Wissenschaften. Springer, Berlin.
38

Welling, M., Teh, Y. W., Andrieu, C., Kominiarczuk, J., Meeds, T., Shahbaba, B., and Vollmer, S.
(2014). Bayesian inference & Big Data: A snapshot from a workshop. ISBA Bulletin, 21(4):8–11.
Woodruff, D. P. and Zhang, Q. (2013). Subspace embeddings and `p -regression using exponential
random variables. In Proc. of COLT, pages 546–567.
Yang, J., Meng, X., and Mahoney, M. W. (2015). Implementing randomized matrix algorithms in
parallel and distributed environments. CoRR, abs/1502.03032.

39

Appendix

●

●

●

●

Spring

●
●

●
●

Summer

●
●
●
●

Fall

●
●

Year2

●
●

●

●
●
●

●

Hour5

●
●
●

●
●

Hour3

●
●

●

Hour2

●

●
●

Hour1

Hour4

10

5

0

−5
Intercept

●
●

●
●

Hour6

●

●
●

Hour7

●
●

●

●

●

Hour11

●

●

Hour9
Hour10

●
●

●

Hour8

●

●

Hour15

●
●

●

Hour14

●
●

●
●

Hour13

●

●
●

Hour12

●

●

Hour16

●
●

Hour23

●
●

●
●

Hour22

●
●

●

Hour21

●
●
●

●

Hour20

●
●

●
●

Hour19

●
●

●
●

Hour18

●
●

●
●

Hour17

●

●

●

●

●

●
●

Thursday

●
●

Wednesday

●

Monday
Tuesday

●
●
●

●

Holiday

●

●
●
●

Heavy clouds

●
●

●
●

Friday
Saturday
●
●

●
●

Rain

●
●

●
●

Humidity

−5

●
●
●

●
●

Wind speed

●
●
●
●
●

●
●
●

App. Temp.

0

5

10

Figure A.1: Boxplots of the MCMC samples for all β parameters of the bike sharing data set based on the original
model

40

