Parallel inversion of huge covariance matrices

An extremely common bottleneck encountered in statistical learning algorithms is inversion of huge covariance matrices, examples being in evaluating Gaussian likelihoods for a large number of data points. We propose general parallel algorithms for inverting positive definite matrices, which are nearly rank deficient. Such matrix inversions are needed in Gaussian process computations, among other settings, and remain a bottleneck even with the increasing literature on low rank approximations. We propose a general class of algorithms for parallelizing computations to dramatically speed up computation time by orders of magnitude exploiting multicore architectures. We implement our algorithm on a cloud computing platform, providing pseudo and actual code. The algorithm can be easily implemented on any multicore parallel computing resource. Some illustrations are provided to give a flavor for the gains and what becomes possible in freeing up this bottleneck.


Introduction
We consider a symmetric positive definite matrix K of order n where n is very large; on the order of 10,000s to millions or more. Our interest is in evaluation of K −1 , which cannot be computed sufficiently quickly using current algorithms. Even if we could compute the inverse, substantial numeric instability and inaccuracies would be expected due to propagation of errors arising from finite machine precision arithmetic (Trefethen and Bau III, 1997).
Typically in statistical applications, one needs to evaluate expressions such as K −1 A, where A is some matrix of appropriate order. Instead of directly computing K −1 , one popular approach is to consider the QR decomposition for K = QR and evaluate the expression, R −1 Q T A (Press et al., 2007). The matrix R is lower triangular and therefore R −1 can be evaluated by backward substitution, which requires O(n 2 ) flops instead of the O(n 3 ) flops for usual inversion. QR decomposition is known to be relatively more stable and efficient than other standard algorithms (Cox and Higham, 1997), a close competitor being the Cholesky decomposition.
The problem is that QR decomposition for K requires O(n 3 ) computations, which is of the same order as that of inversion and therefore prohibitively expensive for large n. Even with QR decomposition being relatively stable, for very large n, finite precision numerical errors are large. This is accentuated when the matrix K has a small spectrum, in the sense of fast decaying eigenvalues, as is typically obtained for covariance matrices obtained from smooth covariance kernels and in large least square problems among a host of other areas. However, this small spectrum, while being the bane of full rank algorithms, makes low-rank algorithms highly accurate.
We propose a new class of algorithms combining ideas from recent developments in linear algebra to find approximate low-rank decompositions. These low-rank decompositions provide several orders of magnitude improvement in speed while also improving accuracy. We propose a general blocking method to enable implementation in parallel on distributed computing architectures. We also consider accuracy of these approximations and provide bounds which are obtained with high probability. Our approach amalgamates ideas from three recent but apparently unrelated developments in numerical linear algebra and machine learning, which we briefly outline below.
There has been an increasing literature on approximating a matrix B of order n × n by BΩ where Ω is a matrix with random entries of order n × d with d << n (Halko et al., 2011). Originally motivated by Johnson Lindenstrauss transforms (Johnson and Lindenstrauss, 1984), these results have been used in a host of application areas, including compressive sensing (Donoho, 2006;Candès et al., 2006) and approximate Gaussian process regression (Banerjee et al., 2013). Typically Ω is populated by random entries, such as draws from Gaussian or Rademacher distributions. The product BΩ involves O(n 2 d) flops and can be itself expensive. Recent developments have shown that structured random matrices, such as random Hadamard transforms, may improve the efficiency significantly, bringing down the number of flops required to the order of n log d (Woolfe et al., 2008). In our case, we show that with certain classes of structured random matrices, Ω of order n × d, we have that KΩ has approximately the same column space as that of K (in a sense to be made precise later), and that the product KΩ may be computed efficiently. In addition, such structured random matrices spread the information of the matrix K, so that we obtain bounds on the minimum and maximum eigenvalues of KΩ, implying that approximate decompositions using KΩ are almost entirely devoid of usual inaccuracies arising from numerical conditioning.
Having efficiently obtained a tall skinny matrix KΩ from the square positive definite matrix K, we now consider recent developments which show that QR decompositions of such tall skinny matrices may be done efficiently in parallel (Constantine and Gleich, 2011;Agullo et al., 2010). The key idea is to partition the tall matrix into blocks and efficiently combine QR decomposition of each of these blocks. The main considerations are the choice of the number of blocks and column size d of the projection matrix Ω, which regulate computation time and accuracy. We provide some theoretical pointers while empirically demonstrating gains achieved on modest architectures. Approximate QR decompositions are useful in obtaining other approximate matrix decompositions, such as approximate spectral decompositions or matrix products. Several recent articles have focused on these issues and Halko et al. (2011) provides an excellent review for algorithms which switch between approximate decompositions. For illustration, we focus on the scenario when the principal interest is in evaluating a Gaussian type likelihood and fine tune the algorithms in this context. We finally consider performance versus other competitors.

Notational Preliminaries
We begin with some notations which we shall use through the rest of the article. In general we will represent matrices in the upper case, A, B etc and row or column vectors in the lower case, a, b etc. We shall use Frobenius and spectral matrix norms: A F = i,j a(i, j) 2 and A 2 = max{ Ax : x = 1} respectively, with the notation A ψ to imply the conclusion holds for both norms. We begin with a realvalued positive definite matrix K ∈ R n×n and k(i, j) will denote the (i, j)th element of K. A spectral decomposition of K = U DU T can be partitioned as, Since K is positive definite, the eigenvalues are positive and we assume without loss of generality that the diagonal matrix D of eigenvalues contains them in descending order of magnitude and the eigenvector matrix U is orthonormal. D m is therefore the matrix containing the m largest eigenvalues and U m the corresponding eigenvector matrix. By the Eckart Young theorem (Stewart, 1993), the best rank m approximation K m to K in terms of minimizing K − K m ψ is given by K m = U m D m U T m . We shall focus on the case where the matrix K has fast decaying eigenvalues. Some examples of problems where this happens is in large least squares optimization problems, Ax−b , with K = A T A, or with K being a covariance matrix generated for a Gaussian process at a dense set of locations with a smooth covariance function (Banerjee et al., 2013;Drineas et al., 2011). For some of our results, we will assume that eigenvalues of K decay at an exponential rate, by which we assume that there exists positive constants λ 1 , λ 2 such that d i,i ≤ λ 1 e λ2i ∀i. λ 1 , λ 2 shall be referred to as the proportionality and rate constants of exponential decay respectively.
An important consideration in this article will be the numerical stability of the algorithms and decompositions. One way in which the stability of a matrix decomposition can be measured is by how well the matrix in question is conditioned. The condition number of a positive semidefinite matrix A is given by c(A) = σ l /σ s , where σ l , σ s are the largest and smallest eigenvalues of the matrix respectively. For the positive definite matrix K it would be c(K) = d(1, 1)/d(n, n).

The Main Algorithm
Our fast approximate inversion algorithm uses three different ideas from numerical techniques and randomized linear algebra. The first is to approximate the spectrum of the large matrix K by post-multiplying it with a random matrix Ω, the second is to consider incorporating special structure in the random matrix Ω for faster evaluation of the product and lastly, blocking strategies, to enable computation of decomposition of the product KΩ on multicore architectures.

Column space approximation
As mentioned, we shall be considering positive definite matrices K with very fast decaying spectrums. In section §4, we provide theoretical justification to show that eigenvalues of positive definite matrices, produced as discrete realizations of a wide class of positive definite kernels, decay very fast. For a positive definite matrix K with fast decaying spectrums, it is reasonable to expect that if Ω is a n × r matrix with r n, populated with independent entries, then the product KΩ will capture most of the information in K, or approximate the column space of K. Such approximation schemes have been in focus in the rapidly expanding field of randomized linear algebra (refer to Halko et al. (2011) for a review). Banerjee et al. (2013) consider this method in the context of approximations for Gaussian processes for large data sets. In general, it is difficult to measure the amount of information captured in KΩ from K. One way to do this is to consider the generalized projection of K onto KΩ and then consider differences. Letting Q be an orthonormal basis for the column space of KΩ, we may consider the error, K − QQ T K ψ , where Q T Q = I, by virtue of orthonormality. In general, we may consider an approximation scheme with two objectives in mind.
1. Fix a target rank m and try to minimize the matrix norm error when using Ω having the target rank.
2. Fix a target matrix norm error and try to achieve that error with high probability with the smallest possible Ω.
We shall consider each of objectives in our projection approximation designs. The simplest possible choice for the random matrix Ω is a submatrix of columns from the n × n identity matrix (this submatrix is sometimes called a permutation submatrix, denoted by P ). This amounts to choosing r columns at random from the matrix K. Low rank approximations using a submatrix have been explored in a variety of contexts, including subset of regressors (Quinonero Candela and Rasmussen, 2005) and least squares (Drineas et al., 2011). With the advent of compressive sensing (Donoho, 2006;Candès et al., 2006), it has been shown that instead of restricting Ω to permutation submatrices, more general random matrices often have better performance (Drineas et al., 2011). Common choices of the random matrix Ω, which have been used in the compressive sensing and matrix approximation literature are: 1. elements of Ω are independent and identically distributed Gaussian random variables appropriately scaled, 2. elements of Ω are independent and identically distributed Rademacher variables, 3. elements of Ω are independently sampled from the uniform distribution on the unit sphere.
In general it has been shown that most choices, based on sampling from a centered distribution and then appropriately scaling, work in achieving accurate error bounds with high probability (Candès et al., 2006).

Approximations via structured random matrices
The biggest drawback of the approximation techniques in §3.4.1 is that the matrix product KΩ may itself be prohibitively expensive O(n 2 r) for very large K. Instead we may form the random matrix Ω such that the product KΩ can be evaluated faster, exploiting properties of the structured distribution.
Definition 1. We shall call a random matrix a structured random matrix if it is of the form Ω = cRT P , where • c is an appropriate scaling constant, such that the columns of Ω are orthonormal, • R is an n × n diagonal matrix whose diagonal elements are independent Rademacher entries, • T is an n × n appropriate orthogonal transform, facilitating the fast multiply, • P is an n × r permutation submatrix.
There are different choices of orthogonal transforms we may use, all of which facilitate the fast multiply, examples for T being the discrete Fourier matrix, discrete cosine matrix, discrete Hartley matrix or the Walsh-Hadamard matrix. The Walsh-Hadamard matrix in particular has been in focus in recent literature (Tropp, 2011) and tight bounds on the approximation accuracy have been obtained. The Walsh Hadamard matrix of order n × n is defined as, The Walsh-Hadamard and discrete Fourier transforms have some optimality with respect to preserving matrix coherence (Tropp, 2011;Boutsidis and Gittens, 2012), but a disadvantage is Walsh-Hadamard transforms exist only for positive integers which are powers of 2, while the discrete Fourier transforms may involve complex numbers. The alternatives, Hartley and discrete cosine transforms while being marginally suboptimal do not have these problems (Avron et al., 2010). In the following section, we formalize the theoretical setting and describe probabilistic error bounds for column space approximations using transforms with general orthonormal matrices (therefore valid for any of the aforementioned transforms).

Blocked decompositions for parallelization
Our focus has been on obtaining a low rank rectangular matrixK of order n × r, with r n, as an approximation to K, such that the column space ofK is very close to the column space of K. To measure the approximation error, we used the error in column space approximation via projection, measured as, ||K − QQ T K ψ , where Q is an orthogonal basis for the column space ofK. Such a Q is typically obtained from the QR decomposition ofK, where Q is an n × r matrix with orthonormal columns and R is an upper triangular matrix (lower triangular in some conventions). Obtaining a QR decomposition of an n × r matrix has computational cost O(nr 2 ). With the matrix multiplication for forming KΩ, the eventual computational cost is O(n 2 r) (for n r), which may be too expensive for very large n. To further reduce computational burden, we may consider parallelizing the computations. Matrix multiplication is trivially parallelizable, while most matrix decompositions are not (Choi et al., 1994). State-ofthe-art linear algebra algorithms, for example in ScaLAPACK routines (Blackford et al., 1997), QR, Cholesky or SV D decompositions, involve dense manipulations of the whole matrix, which are not parallelizable or cannot be blocked. In addition to computational cost, storing the full matrixK may be problematic in terms of memory requirements for very large n. However, in our case, we are interested in the QR decomposition of a tall matrix, number of rows being much greater than the number of columns, for which blocking strategies have been recently developed, exploiting modern parallel computing platforms like MapReduce. We elaborate on the blocking strategies next.
For notational clarity assume that n = 8r and consider a partition ofK, In the above each K i is of size 2r × r. Denote the QR factorization of the small matrix K i = Q i R i . Along the lines of the T SQR factorization presented in Constantine and Gleich (2011) we consider the following blocked factorization, Each of the small matrices Q i being orthogonal, the left factor of the decomposition in (1), Q b 1 , is formed by stacking the small matrices Q i as diagonal blocks, in an orthogonal matrix. It is also an orthogonal basis for the column space ofK and in many application areas it suffices to work with Q b 1 . In some other applications, we may require evaluation of the product Q b 1 x, where x is some vector of order r × 1, which can be evaluated in a similar blocked fashion, utilizing the small Q i 's. The matrix R b 1 , formed by stacking the small upper triangular matrices R i , is not upper triangular and hence the decomposition Q b 1 R b 1 is not a QR decomposition ofK. To achieve the QR decomposition ofK consider the following sequence of blocked and finally the QR decomposition of the 2r × r matrix R b 2 , where we can give the final QR decomposition ofK as, R b 3 by virtue of the construction is an upper triangular matrix, whileQ, a product of orthogonal matrices is itself orthogonal. A similar iteration can be carried out for any n, r by appropriately adjusting the block sizes. In some instances, if only (R b 3 ) −1 is required, the matrix product for formingQ need not be evaluated at all. On the other hand, if a product of the formQ T x is required for some some vector x, the entire product maybe evaluated via blocking. We give more details about such strategies in our examples.
There are other ways of constructing a sequence of small blocked QR decompositions leading to a QR decomposition of the tall matrixK. The scheme presented in (4) is the one allowing for maximum assimilation of computations at the same time. Depending on the architecture being used, this may not be the most optimal choice. Below we present another analogous blocked QR scheme, which allows for lesser number of computations to be carried on simultaneously, but possibly uses lesser communication between units of a multicore architecture. Assume, again for notational clarity that n = 8r. The matrixK is now partitioned into 8 blocks, where, Q 1 R 1 is the orthogonal matrix from the QR decomposition of K 1 , Q 2 is from the QR decomposition of R 1 K 2 , etc and I denotes identity matrix of appropriate size. The matrices R 1 , R 2 . . . are intentionally omitted from (5)  . This alternative blocking scheme involves minimum communication between units of the multicore architecture. The optimal schemes will usually be a mixture of the strategies (2) and (5). For more variants of blocking algorithms to get QR decompositions for tall matrices, we refer the readers to Constantine and Gleich (2011); Agullo et al. (2010). We now consider the theoretical gains possible from the blocked algorithm, while considering block size. Consider the ideal scenario, when there is no overhead or communication cost and each unit of the multicore architecture has equal speed. We consider the strategy corresponding to (2). Each of the small QR decompositions involve an 2r × r matrix, which has O(r 3 ) cost. At the first iteration, as in (2), using the same number of units as the number of blocks, say b, we can perform all the computations in O(r 3 ). In fact, it is easy to see that the computational order of O(r 3 ) holds true for each of the subsequent iterations. The total number of iterations needed, following the first strategy would be log 2 (b) + 1, which brings the total computational cost to O(log 2 (b)r 3 ) for the QR ofK. Forming the product KΩ and additional matrix multiplications (which will almost always be needed in applications after the QR decomposition), bring the cost to O(n 2 r) without parallelization, whereas with the blocking algorithm, provided we perform the matrix multiplications with the same blocked structure, we achieve the whole computation in O(log 2 (b)r 3 ). Following the first strategy and from the above discussion, in the ideal scenario, the number of blocks to be used is b = n 2r + 1. The theoretical possible maximum speed-up, assuming no communication cost, is where we measure the speed-up from Amdahl's equation (Mattson et al., 2005), as proportion of speed of the unparallelized computations to the speed of the parallelized computations. In this subsection we derive some motivating results justifying our assertion in §3.1 that positive definite matrices obtained as dense discrete realizations of positive definite kernels will have a very fast decaying spectrum, enabling us to obtain good approximations using much lower rank matrices. The ideas are not entirely new, we adapt abstract results from theory of integral equations and stochastic series expansions to our context. To prepare the background, we start with the well-known Mercer's theorem for positive definite functions.
Definition 2. Let D be a compact metric space and C a function, C : D × D → R + ∪ {0}. C is said to be positive definite if for all n, scalars c 1 , . . . , c n ∈ R and x 1 , . . . , x n ∈ D, we have i j c i c j C(x i , x j ) > 0. It follows trivially from the definition that if K is an n × n matrix with k(i, j) = c(x i , x j ) for x 1 , . . . , x n ∈ D, then K is a positive definite matrix.
With this definition, we give the following version of Mercer's theorem for the positive definite function C(·, ·) (Kühn, 1987), where D is a compact metric space, there exists a real valued scalar sequence, {λ i } ∈ l 1 , λ 1 ≥ λ 2 ≥ . . . ≥ 0 and an orthonormal sequence of functions, {φ i (·)} ∈ L 2 (D), φ i (·) : D → R, such that, ∀x 1 , x 2 ∈ D and this sequence converges in the mean square sense and uniformly in D.
We omit the the proof of this theorem and refer the readers to (Kühn, 1987).Mercer's theorem is a generalization of the spectral decomposition for matrices extended to functions. {λ i }, {φ(·)} are referred to as the eigenvalues and eigenfunctions respectively of C(·, ·). Alternatively, suppose D ⊂ R d and considering Lebesgue measure, we have the following integral equation for the eigenvalues and eigenfunctions, analogous to the definition of matrix eigenvalues and eigenfunctions, ∀x ∈ D, i ∈ N. The eigenvalues and eigenfunctions of C(·, ·) can be approximated using the eigenvalues and eigenvectors of the positive definite matrix K, obtained by evaluating C(·, ·) at any set of n locations x 1 , . . . x n ∈ D. To see this, consider a discrete approximation of the integral equation (6), using the points x 1 , . . . , x n , with equal weights in the quadrature rule, Substituting y = x 1 , . . . , x n in (7) we get a system of linear equations which is equivalent to the matrix eigenproblem, KU i = d(i, i)U i , i = 1, . . . , n, where U i is the i th row of the eigenvector matrix U . This approximation is related to the Galerkin method for integral equations and the Nystrom approximation (Baker and Baker, 1977;Delves and Mohamed, 1988). It also corresponds to finite truncation of the expansion from the Mercer's theorem. In general it can be shown that d(i, i)/n and √ n u(i, j) converge to λ i and φ i (x j ) respectively as n → ∞ (Baker and Baker, 1977).
The accuracy of the finite truncation of the Mercer expansion has been in recent focus and error bounds have been obtained depending on the degree of smoothness of the positive definite function (Todor, 2006;Schwab and Todor, 2006), in the context of stochastic uncertainty quantification. To borrow from these ideas and adapt the abstract results in our case, we begin with definitions quantifying the smoothness of the covariance functions.
Definition 3. Let C(·, ·) be a positive definite function on D × D, where D is a compact metric space. We call C(·, ·) piecewise Sobolev (p, q) smooth, for some p, q ∈ N, if there exists a finite disjoint partition {D j , j ∈ J} of D, with D ⊂ ∪ j∈J D j , such that for any pair (j 1 , j 2 ) ∈ J, C(·, ·) is Sobolev (p, q) on D j2 ×D j2 .
Definition 4. Let C(·, ·) be a positive definite function on D × D, where D is a compact metric space. We call C(·, ·) piecewise analytic smooth, for some p, q ∈ N, if there exists a finite disjoint partition {D j , j ∈ J} of D, with D ⊂ ∪ j∈J D j , such that for any pair (j 1 , j 2 ) ∈ J, C(·, ·) is analytic on D j2 × D j2 .
Commonly used covariance functions such as the squared exponential, C(x, y) = c 1 exp(−c 2 x − y 2 ) and the Matern function, depending on the value of the parameter ν, fall into one of the above categories. In fact the squared exponential function is Sobolev smooth on any compact domain and we can give stronger results for decay of its eigenvalues.
Using the couple of definitions categorizing smoothness, we give the following result, which is analogous to a functional version of the Eckart-Young theorem.
Lemma 2. Let D ⊂ R d . For any m ∈ N, and let F m be an m dimensional closed subspace of positive definite functions on D × D, where by m dimensional we mean that ∃ sequence of orthonormal function {ψ j (·)} m j=1 on D such that{ψ 2 j (·)} m j=1 spans F m , with positive coefficients by virtue of positive definiteness. Let C Fm (·, ·) be the projection of C(·, ·) onto F m , and infimum of the errors in approximating C(·, ·) by m dimensional functions ∈ F m , ie E m = inf Fm { C − C Fm 2 , C m ∈ F m }. Then, (i) If C(·, ·) is piecewise Sobolev (p, q) ∀(p, q) ∃ for each s > 1, a positive constant c s depending only on C(·, ·), s, such that E m ≤ c s j≥m j −s . Henceforth we call this bound E s m .
(ii) If C(·, ·) is piecewise analytic ∃ , positive constants c 1 , c 2 depending only on C(·, ·), such that E m ≤ c s j≥m c 1 exp c2m 1/d . Henceforth we call this bound E a m .
Proof. Let the Mercer expansion for C(x, y) be j≥1 λ j φ j (x)φ j (y). Define a random function, f (x) = j≥1 η j λ j φ j (x), where η j are independent standard normal random variables. Let C Fm (x, y) = m j=1 θ j ψ j (x)ψ j (y) corresponding to the basis {ψ 2 (·)} m j=1 for F m . Similar to f , corresponding to C Fm (·, ·), define a random function, f Fm (x) = j≥1 ζ j θ j ψ j (x), where ζ j 's are independent random normal variables. Defining f − f Fm = E(f − f Fm ) 2 , and applying theorem 2.7 in Schwab and Todor (2006), we get that, E m = j≥m λ j . Then the result follows by applications of Corollary 3.3 and Proposition 3.5 in Todor (2006).
The above lemma quantifies the finite truncation accuracy of the Mercer theorem for smooth kernels. The bounds obtained are optimal and in general cannot be improved. The spaces F m 's can be made more general to encompass all square integrable bivariate functions, but the proof becomes more involved in that case and we omit it for the sake of brevity. In case of the squared exponential kernel, which is smooth of all orders, the finite truncation is even sharper and we give the following corollary for it: Corollary 3. Let C(x, y) = θ 1 exp −θ2 x−y 2 , for x, y ∈ D, compact ⊂ R d and θ 1 , θ 2 > 0. Using all notations as in lemma 2, ∃ positive constant c θ1,θ2 such that E m ≤ c θ1,θ2 j≥m θ j 1/d 2 Γ(j 1/d /2) . Henceforth we shall call this bound E e r .
The proof is exactly similar to the proof of Lemma 2 and an application of Proposition 3.6 in Todor (2006).
We now describe some results relating to the matrix eigenvalues, using the strength of Lemma 2 above. Proof. We begin with the discrete approximate solution of the integral equation, using the Galerkin technique, described in equation (7). Note that this method coincides with the Raleigh-Ritz method with the identity matrix in Baker and Baker (1977), since C(·, ·) is symmetric positive definite. Applying then theorem 3.31 on page 322 of Baker and Baker (1977), we have d(m, m) ≤ nλ m , where λ m is the m th eigenvalue from the Mercer expansion of C(·, ·). The result then follows by a straight-forward application of Lemma 2 above and Corollary 3.3, Proposition 3.5 in Todor (2006). For the squared exponential kernel, it is straightforward application of Corollary 3 above and Proposition 3.6 in Todor (2006).
This theorem shows that for most covariance functions, the positive definite matrices that are generated by their finite realizations have eigenvalues that decay extremely fast, which was our assertion at the start of the section. In the simulated experiments we compute eigenvalues of some such covariance kernels and show that empirical results support our theory.

Approximation accuracy and condition numbers
We now present some results regarding the accuracy of our approximation algorithms. We shall be concerned with the column space approximation error, when we useK = Ω instead of K, measured as K − QQ T K ψ , where Ω is a structured random matrix as in definition 1 and Q is orthonormal basis for the columns ofK. Such a Q can be obtained from the QR decomposition ofK.
Let Ω be an n × r structured random matrix as formulated in definition 1 and Q be the left factor from the QR decomposition of KΩ. Choose r, k such that 4[ √ k + √ 8 ln kn] 2 ln K ≤ r ≤ n. With probability at least (1 − O(1/k)), the following hold true, (i) If C(·, ·) is piecewise Sobolev smooth ∀(p, q), then (a) K − QQ T K 2 ≤ (1 + 7n/r)c s r −s , where c s is a positive constant depending on s for any s > 1; (b) K − QQ T K F ≤ (1 + 7n/r)E s r . (ii) If C(·, ·) is piecewise analytic, then (a) K − QQ T K 2 ≤ n(1 + 7n/r)c 1 exp c2m 1/d , where c 1 , c 2 are positive constants ; (b) K − QQ T K F ≤ n(1 + 7n/r)E a r . (iii) In particular, for a squared exponential kernel, C(x, y) = θ 1 exp −θ2 x−y 2 , then (a) K − QQ T K 2 ≤ (1 + 7n/r)c θ1,θ2 θ m 1/d 2 Γ(m 1/d /2) , where c θ1,θ2 is a positive constant depending on θ 1 , θ 2 ; (b) K − QQ T K F ≤ (1 + 7n/r)E e r , where E s r , E a r , E e r are bounds as in Theorem 2 and Corollary 3. Proof. First note that, where ((K) TK ) + denotes the Moore-Penrose inverse of ((K) TK ). Let Q k R k denote the QR decomposition ofK, so that we have, Let K r be the best rank r approximation to K and let the spectral decomposition of K = U DU T . Then from the Eckart-Young theorem, K − K m 2 = d(r, r) and K − K m 2 = n j=r d(j, j). The result then follows by an application of our Theorem 4 and Theorem 11.2 in Halko et al. (2011).
In addition to providing accurate column space approximations, these orthogonal transforms spread the information of the matrix K and improve its conditioning, while preserving its geometry. By preserving the geometry, we mean preserving the norms of the eigenvectors -we explore more of this aspect empirically in the simulations. Any low rank approximation improves the conditioning, but it has been shown (Boutsidis and Gittens, 2012;Avron et al., 2010) that projections using the orthogonal transforms as above, improve conditioning substantially beyond what is achieved by just using any low rank approximation with high probability in special cases. Halko et al. (2011); Boutsidis and Gittens (2012) obtain bounds on the eigenvalues of AΩ where Ω is a structured random matrix and A is orthogonal. We are in interested in the stability of the numerical system,KR −1 k , as explained in the introduction of this article, where the QR decomposition ofK = Q k R k . We present the following result, without proof, trivially modifying results from Boutsidis and Gittens (2012); Avron et al. (2010), Theorem 6. Fix 0 < < (1/3) and choose 0 < δ < 1, r such that 6 2 [ (r) + 8 ln n/δ] 2 ln 2n/δ ≤ r ≤ n. Then with probability at least 1 − 2δ, we have the condition number of the linear system of interest, This improvement in conditioning number causes huge improvement in learning algorithms over other competing approaches as we demonstrate later.

Simulation Examples
We begin with empirical investigation into the eigenvalue decays of matrices realized from commonly used covariance kernels. We consider an equispaced grid of 100 points in [0, 1], to generate the 100×100 covariance matrices with 100 positive eigenvalues. We start with such a moderate size, 100, to illustrate that even for this small size, eigenvalues decay extremely fast, as was hypothesized and theoretically proven earlier.
This first one we use is the squared exponential kernel, The parameter θ 1 is a scaling parameter and does not change the eigendirections. In our present simulation, we set θ 1 = 1. The parameter θ 2 controls the smoothness of the covariance, as it decays with the distance x − y . We consider a range of θ 2 values, 0.05, 0.5, 1, 1.5, 2, 10, respectively, considering negligible decay with distance and very fast decay with distance. The plot of the eigenvalues is shown in figure 1. The figure shows similar rates of decay across the range of values of the smoothness parameter θ 2 . Depending on the smoothness parameter, the value of the largest eigenvalue changes by several orders, but in all cases, by the 10th largest eigenvalue, the sequence has approximately converged to 0, which indicates that for a low rank approximation, a choice of rank 10 would be sufficient for our purposes. This second one we use is the Matern covariance kernel (Williams and Rasmussen, 1996), where B ν is the modified Bessel function of the second kind. The parameter θ 1 is a scaling parameter and does not change the eigendirections, analogous to θ 1 of the squared exponential kernel. In our present simulation, we set θ 1 = 1, θ 2 = 1. The parameter ν controls the smoothness of the covariance, analogous to the parameter θ 2 of the squared exponential kernel. As a special case, as ν → ∞, we obtain the squared exponential kernel. In typical spatial applications, ν is considered to be within [0.5, 3], as the data are rarely informative about the smoothness levels beyond these. We consider a range of ν values, 0.5, 1, 1.5, 2, 2.5, 3, respectively, considering negligible decay with distance and very fast decay with distance. The plot of the eigenvalues is shown in figure 2. The figure shows similar rates of decay across the range of values of the smoothness parameter ν, exactly as exhibited with the squared exponential kernel. Once again, a low rank approximation of rank 10 would suffice in this case for approximating the 100 × 100 matrix.
Using the same set of simulations, we compute conditioning numbers, in turn for each of the covariance kernels. We compute the conditioning number of the the full 100 × 100 covariance matrix and its best rank m approximations for m = 5, 10, 15, 20, 50 respectively. Results for the squared exponential kernel are tabulated in table 1 and for the Matern covariance kernel in table 2. The full covariance matrix is extremely ill-conditioned in each case, the best rank m approximations improve the conditioning, as is to be expected. Even with the best rank m approximations, the conditioning numbers are still very large indicating very unstable linear systems. In the table we omit the exact digits for the very large numbers and just indicate their orders in terms of powers of 10. In case of the squared exponential kernel, condition numbers decrease from left to right and from top to bottom. This pattern is in general not true for the condition numbers of the Matern kernel, whose condition numbers are several orders smaller than that for the squared exponential but still to large for it to be stable linear systems. In general for the Matern kernel, as ν becomes larger, the kernel becomes smoother, the condition numbers are larger and the eigenvalues decay faster.
We next move to the simulations to consider the effect of blocking and gain in efficiency from parallelization. For large sample sizes, we will get very poor and time consuming inverse estimates, as demonstrated by the large conditioning numbers from the simulations of the previous section. To circumvent this, we start therefore with a known spectral decomposition and apply our approximations, pretending that the decomposition was unknown. We consider the known spectral decomposition as K = EDE T , where E is an orthonormalized matrix with randomly generated iid Gaussian entries. D is the diagonal matrix of eigenvalues in decreasing order or magnitude and for this simulation example we assume exponential decay with scale = 1 and rate 0.01. We generate this covariance matrix K for sample sizes n = 1000, 5000, 10000, 50000 respectively and apply our blocked approximate inversion algorithm. For the random matrix Ω used for the projection, we use three different choices, a matrix with scaled iid standard Gaussian entries, a structured random Hartley transform and a structured random discrete cosine transform. We have a total 64 cores at our disposal and to get a flavor of the gains possible by parallelization, we use 8, 16, 32, 64 cores in turn and measure the gain in efficiency as the ratio of time taken to the time for the algorithm without parallelization.
We report the gain in efficiency in figure 3. The figure reveals the efficiency gain by using the blocked algorithm. It shows that efficiency gain increases with number of cores in the larger sample sizes. Specifically when the sample size is 50, 000, it seems that more than 64 cores could have further increased the efficiency. The gains reported are obviously lower than the theoretical maximum possible gains, but are still substantial and have potential of further increase with larger number of cores in huge sample sizes.

Real data examples
We consider a real data examples with a very large data set, to demonstrate the large data handling power and efficiency gain via our approach. We consider subsets of the actual data available and variables such that the data fit into a Gaussian process framework.
The example we consider is a subset of the Sloan sky digital survey (Stoughton et al., 2007). In the survey, redshift and photometric measurements were made for a very large number of galaxies. The photometric measurements consist of 5 band measurements and are available for most of the galaxies, while the redshift measurements are not available for a large number of galaxies due to spectroscopy constraints. The main interest is therefore in predicting the red-shifts, given the other 5 measurements. For this example we consider a subset comprising of 180, 000 galaxies, whose photometric measurements and red-shift measurements are known. We hold out 30, 000 galaxies and pretend that redshifts are unknown for these galaxies, while using the remaining 150, 000 to train our Gaussian process model.
For the computations we use both a squared exponential kernel and a Matern kernel. The data are scaled and centered before calculations so that we are not concerned with inference regarding the θ 1 parameters for the covariance kernel. For the Matern kernel, we use a discrete uniform prior for the parameter ν based on 100 equispaced points on the grid [0.5, 3]. For the squared exponential kernel, we use a 100 trials, each having a random selection 1, 000 points from the training set to set a band for θ 2 , to estimate range of covariance of the data. It appears that it is sufficient to consider θ 2 in the range [0.05, 1]. We therefore use a grid on 100 uniformly spaced points in [0.05, 1] and analogous to the prior for ν, place a discrete uniform prior on θ 2 of the squared exponential kernel.
For data sets of this size it becomes extremely inefficient to use predictive processes or knot based approaches per se, without blocking or parallelization and knot selection is almost practically infeasible to attempt. We work with the modified predictive processes approach of Finley et al. (2009), without knot selection as competitors to our approach. The modified predictive processes is conceptually equivalent to the fully independent training conditional approach (Quinonero Candela and Rasmussen, 2005) or the Gaussian process approximation with pseudo inputs, with (Snelson and Ghahramani, 2006), with the same priors for the parameters as our approach, so that the only difference is in the way the covariance inversion is done.
For each of the approaches, including ours, one issue is how to select rank of the approximate projections. We use an approximate estimate of the covariance decay to come up with the low rank projection to be used. To calculate this we select 1, 000 data points spread out across the training set, in the following approximate manner, first select the pair of points which are most distant from each other, then select the point which is almost equidistant from the pair selected, and so on, till we have a well spread out selection. Here distance is measured as the Euclidean distance between the 5 variate photometric measurement vectors. Using these 1, 000 data points, we calculate the prior covariance matrix, using the Gaussian covariance kernel and the Matern kernel for each value of the smoothness parameter on its grid of values (Our simulation examples have revealed that the decay of eigenvalues depends on this smoothing parameter). We then estimate the rank such that Frobenius norm error in using the best low rank approximation from the Eckart Young theorem for the full covariance matrix would be no more than 0.001 (bearing in mind that this is an estimate of the best possible error, actual error in different projections will be more). In the computations for our approach, we use a different random projection to the targeted rank for each grid value of the smoothing parameter, while for the knot based approaches, a different random selection of knots is used.
For our approach, we use 3 different choices of the projection matrix Ω, a matrix with scaled iid standard Gaussian entries, a structured random matrix with the Hartley transform and a structure random matrix with the discrete cosine transform. We use blocking as in our main algorithm for each of the three types of projection matrices on all 64 cores of the distributed computing structure. The knot based approach has no obvious method of blocking or for it to be used on a grid computing environment. The MCMC algorithm in each case is run for 10, 000 iterations, with the first 2, 000 discarded as burn-in. Usual diagnostics do not show evidence of non-convergence of the MCMC algorithms and we report the effective sample size for the smoothing parameter as obtained from the CODA package in R Plummer et al. (2006). Time taken by each method is computed and reported. The time calculations are made by recording the exact CPU times for iterations, excluding the time for Frobenius norm accuracy and conditioning number calculations and also excluding the time taken to set up parallel Matlab workers via matlabpool. We measure predictive accuracy, of the predictions in the hold out set, based on the relative mean squared error.
We tabulate the results from the experiment in the table 3. The performance of any of the projection based approaches is substantially better than the knot based PP in terms of any of the parameters reported, predictive accuracy, time taken or the effective sample size of the smoothing parameter. Time taken by the PP in case of this large dataset was several days, as compared to a few hours for the blocked projection approaches. The improvement in predictive accuracy can probably be attributed to much better conditioned stable linear systems being used inherently and hence improved inference. The discrete cosine transform and Hartley transform perform marginally better than the random Gaussian projection, however the structured random transforms show marked improvement in the time taken, possibly due to the faster matrix multiplies.
There is little to choose between the discrete cosine transform or the discrete Hartley transform.

Discussion
In this article we present a new method of blocked approximate inversion for large positive definite matrices, using subsampled random orthogonal transforms, motivated by several recent developments in random linear algebra. The blocking strategy comes from development of parallel QR algorithms for tall matrices and in our case tall matrices are obtained from large positive definite matrices as the first projection step of a low rank approximation. We have also presented a comprehensive set of theoretical results quantifying the decay of eigenvalues of covariance matrices, which in turns helps bound the errors from the low rank projections. The examples show marked improvement in terms of numerical stability, prediction accuracy and efficiency and should be routinely applicable in a wide variety of scenarios.
One interesting future direction is the investigation of several structured random projections taken together. One may break up a target rank into several parts, the first of which is spent on a random projection, and the remaining pieces are learnt through the direction of maximum accuracy. Both theoretical (finding projections in the directions of maximum accuracy, which may be the direction of maximum information descent) and empirical investigation of this issue is currently being investigated. Table 1: Table of condition numbers for the squared exponential kernel for different values of the smoothing parameter θ 2 versus truncation levels, for truncating to the best possible approximation according to the Eckart-Young theorem. The rows represent the levels of truncation and the columns, the values of the smoothing parameter.