Lecture 21 Computational Methods for GPs Colin Rundel 04/10/2017 - - PowerPoint PPT Presentation

Lecture 21

Computational Methods for GPs

Colin Rundel 04/10/2017

1

GPs and Computational Complexity

2

The problem with GPs

Unless you are lucky (or clever), Gaussian process models are difficult to scale to large problems. For a Gaussian process 𝐳 ∼ 𝒪(𝝂, 𝚻): Want to sample 𝐳?

𝝂 + Chol(𝚻) × 𝐚 with 𝑎𝑗 ∼ 𝒪(0, 1) 𝒫 (𝑜3)

Evaluate the (log) likelihood?

−1 2 log |Σ| − 1 2(𝐲 − 𝝂)′ 𝚻−1 (𝐲 − 𝝂) − 𝑜 2 log 2𝜌 𝒫 (𝑜3)

Update covariance parameter?

{Σ}𝑗𝑘 = 𝜏2 exp(−{𝑒}𝑗𝑘𝜚) + 𝜏2

𝑜 1𝑗=𝑘

𝒫 (𝑜2)

3

The problem with GPs

Unless you are lucky (or clever), Gaussian process models are difficult to scale to large problems. For a Gaussian process 𝐳 ∼ 𝒪(𝝂, 𝚻): Want to sample 𝐳?

𝝂 + Chol(𝚻) × 𝐚 with 𝑎𝑗 ∼ 𝒪(0, 1) 𝒫 (𝑜3)

Evaluate the (log) likelihood?

−1 2 log |Σ| − 1 2(𝐲 − 𝝂)′ 𝚻−1 (𝐲 − 𝝂) − 𝑜 2 log 2𝜌 𝒫 (𝑜3)

Update covariance parameter?

{Σ}𝑗𝑘 = 𝜏2 exp(−{𝑒}𝑗𝑘𝜚) + 𝜏2

𝑜 1𝑗=𝑘

𝒫 (𝑜2)

3

The problem with GPs

Unless you are lucky (or clever), Gaussian process models are difficult to scale to large problems. For a Gaussian process 𝐳 ∼ 𝒪(𝝂, 𝚻): Want to sample 𝐳?

𝝂 + Chol(𝚻) × 𝐚 with 𝑎𝑗 ∼ 𝒪(0, 1) 𝒫 (𝑜3)

Evaluate the (log) likelihood?

−1 2 log |Σ| − 1 2(𝐲 − 𝝂)′ 𝚻−1 (𝐲 − 𝝂) − 𝑜 2 log 2𝜌 𝒫 (𝑜3)

Update covariance parameter?

{Σ}𝑗𝑘 = 𝜏2 exp(−{𝑒}𝑗𝑘𝜚) + 𝜏2

𝑜 1𝑗=𝑘

𝒫 (𝑜2)

3

The problem with GPs

Unless you are lucky (or clever), Gaussian process models are difficult to scale to large problems. For a Gaussian process 𝐳 ∼ 𝒪(𝝂, 𝚻): Want to sample 𝐳?

𝝂 + Chol(𝚻) × 𝐚 with 𝑎𝑗 ∼ 𝒪(0, 1) 𝒫 (𝑜3)

Evaluate the (log) likelihood?

−1 2 log |Σ| − 1 2(𝐲 − 𝝂)′ 𝚻−1 (𝐲 − 𝝂) − 𝑜 2 log 2𝜌 𝒫 (𝑜3)

Update covariance parameter?

{Σ}𝑗𝑘 = 𝜏2 exp(−{𝑒}𝑗𝑘𝜚) + 𝜏2

𝑜 1𝑗=𝑘

𝒫 (𝑜2)

3

A simple guide to computational complexity

𝒫 (𝑜) - Linear complexity

• Go for it

𝒫 (𝑜2) - Quadratic complexity - Pray 𝒫 (𝑜3) - Cubic complexity - Give up

4

A simple guide to computational complexity

𝒫 (𝑜) - Linear complexity - Go for it 𝒫 (𝑜2) - Quadratic complexity - Pray 𝒫 (𝑜3) - Cubic complexity - Give up

4

A simple guide to computational complexity

𝒫 (𝑜) - Linear complexity - Go for it 𝒫 (𝑜2) - Quadratic complexity

• Pray

𝒫 (𝑜3) - Cubic complexity - Give up

4

A simple guide to computational complexity

𝒫 (𝑜) - Linear complexity - Go for it 𝒫 (𝑜2) - Quadratic complexity - Pray 𝒫 (𝑜3) - Cubic complexity - Give up

4

A simple guide to computational complexity

𝒫 (𝑜) - Linear complexity - Go for it 𝒫 (𝑜2) - Quadratic complexity - Pray 𝒫 (𝑜3) - Cubic complexity

• Give up

4

A simple guide to computational complexity

𝒫 (𝑜) - Linear complexity - Go for it 𝒫 (𝑜2) - Quadratic complexity - Pray 𝒫 (𝑜3) - Cubic complexity - Give up

4

How bad is the problem?

10 20 30 2500 5000 7500 10000

n time (secs) method

chol inv LU inv QR inv

5

Practice - Migratory Model Prediction

After fitting the GP need to sample from the posterior predictive distribution at ∼ 3000 locations

𝐳𝑞 ∼ 𝒪 (𝜈𝑞 + Σ𝑞𝑝Σ−1

𝑝 (𝑧𝑝 − 𝜈𝑝), Σ𝑞 − Σ𝑞𝑝Σ−1 𝑝 Σ𝑝𝑞)

Step CPU (secs)

• 1. Calc. Σ𝑞, Σ𝑞𝑝, Σ𝑞

1.080

• 2. Calc. chol(Σ𝑞 − Σ𝑞𝑝Σ−1

𝑝 Σ𝑝𝑞)

0.467

• 3. Calc. 𝜈𝑞|𝑝 + chol(Σ𝑞|𝑝) × 𝑎

0.049

• 4. Calc. Allele Prob

0.129 Total 1.732

Total run time for 1000 posterior predictive draws:

• CPU (28.9 min)

6

Practice - Migratory Model Prediction

After fitting the GP need to sample from the posterior predictive distribution at ∼ 3000 locations

𝐳𝑞 ∼ 𝒪 (𝜈𝑞 + Σ𝑞𝑝Σ−1

𝑝 (𝑧𝑝 − 𝜈𝑝), Σ𝑞 − Σ𝑞𝑝Σ−1 𝑝 Σ𝑝𝑞)

Step CPU (secs)

• 1. Calc. Σ𝑞, Σ𝑞𝑝, Σ𝑞

1.080

• 2. Calc. chol(Σ𝑞 − Σ𝑞𝑝Σ−1

𝑝 Σ𝑝𝑞)

0.467

• 3. Calc. 𝜈𝑞|𝑝 + chol(Σ𝑞|𝑝) × 𝑎

0.049

• 4. Calc. Allele Prob

0.129 Total 1.732

Total run time for 1000 posterior predictive draws:

• CPU (28.9 min)

6

A bigger hammer?

Step CPU (secs) CPU+GPU (secs)

• Rel. Perf
• 1. Calc. Σ𝑞, Σ𝑞𝑝, Σ𝑞

1.080 0.046 23.0

• 2. Calc. chol(Σ𝑞 − Σ𝑞𝑝Σ−1

𝑝 Σ𝑝𝑞)

0.467 0.208 2.3

• 3. Calc. 𝜈𝑞|𝑝 + chol(Σ𝑞|𝑝) × 𝑎

0.049 0.052 0.9

• 4. Calc. Allele Prob

0.129 0.127 1.0 Total 1.732 0.465 3.7

Total run time for 1000 posterior predictive draws:

• CPU (28.9 min)
• CPU+GPU (7.8 min)

7

Cholesky CPU vs GPU (P100)

10 20 30 2500 5000 7500 10000

n time (secs) method

chol inv LU inv QR inv

comp

cpu gpu

8

0.1 10.0 2500 5000 7500 10000

n time (secs) method

chol inv LU inv QR inv

comp

cpu gpu

9

Relative Performance

1 10 2500 5000 7500 10000

n Relative performance method

chol inv LU inv QR inv

10

Aside (1) - Matrix Multiplication

0.0 2.5 5.0 7.5 2500 5000 7500 10000

n time (sec) comp

cpu gpu

Matrix Multiplication

11

20 25 30 35 40 45 2500 5000 7500 10000

n time (sec)

Matrix Multiplication − Relative Performance

12

Aside (2) - Memory Limitations

A general covariance is a dense 𝑜 × 𝑜 matrix, meaning it will require 𝑜2× 64-bits to store.

5 10 15 20 10000 20000 30000 40000 50000

n Cov Martrix Size (GB) 13

Other big hammers

bigGP is an R package written by Chris Paciorek (UC Berkeley), et al.

• Specialized distributed implementation of linear algebra operation for GPs
• Designed to run on large super computer clusters
• Uses both shared and distributed memory
• Able to fit models on the order of 𝑜 = 65k (32 GB Cov. matrix)

0.01 0.1 1 10 100 1000 2048 8192 32768 131072 Execution time, seconds (log scale) Matrix dimension, n (log scale) Cholesky Decomposition 6 cores 60 cores 816 cores 12480 cores 49920 cores

14

More scalable solutions?

• Spectral domain / basis functions
• Covariance tapering
• GMRF approximations
• Low-rank approximations
• Nearest-neighbor models

15

Low Rank Approximations

16

Low rank approximations in general

Lets look at the example of the singular value decomposition of a matrix,

𝑁

𝑜×𝑛 = 𝑉 𝑜×𝑜 diag(𝑇) 𝑜×𝑛

𝑊 𝑢

𝑛×𝑛

where 𝑉 are called the left singular vectors, 𝑊 the right singular vectors, and 𝑇 the singular values. Usually the singular values and vectors are

• rdered such that the singular values are in descending order.

The Eckart–Young theorem states that we can construct an approximatation

• f 𝑁 with rank 𝑙 by setting

̃ 𝑇 to contain only the 𝑙 largest singular values

and all other values set to zero.

̃ 𝑁

𝑜×𝑛 = 𝑉 𝑜×𝑜 diag( ̃

𝑇)

𝑜×𝑛

𝑊 𝑢

𝑛×𝑛

= ̃ 𝑉

𝑜×𝑙

diag( ̃

𝑇)

𝑙×𝑙

̃ 𝑊 𝑢

𝑙×𝑛

17

Low rank approximations in general

Lets look at the example of the singular value decomposition of a matrix,

𝑁

𝑜×𝑛 = 𝑉 𝑜×𝑜 diag(𝑇) 𝑜×𝑛

𝑊 𝑢

𝑛×𝑛

where 𝑉 are called the left singular vectors, 𝑊 the right singular vectors, and 𝑇 the singular values. Usually the singular values and vectors are

• rdered such that the singular values are in descending order.

The Eckart–Young theorem states that we can construct an approximatation

• f 𝑁 with rank 𝑙 by setting

̃ 𝑇 to contain only the 𝑙 largest singular values

and all other values set to zero.

̃ 𝑁

𝑜×𝑛 = 𝑉 𝑜×𝑜 diag( ̃

𝑇)

𝑜×𝑛

𝑊 𝑢

𝑛×𝑛

= ̃ 𝑉

𝑜×𝑙

diag( ̃

𝑇)

𝑙×𝑙

̃ 𝑊 𝑢

𝑙×𝑛

17

Example

𝑁 = ⎛ ⎜ ⎜ ⎜ ⎝ 1.000 0.500 0.333 0.250 0.500 0.333 0.250 0.200 0.333 0.250 0.200 0.167 0.250 0.200 0.167 0.143 ⎞ ⎟ ⎟ ⎟ ⎠ = 𝑉 diag(𝑇) 𝑊 𝑢 𝑉 = 𝑊 = ⎛ ⎜ ⎜ ⎜ ⎝ −0.79 0.58 −0.18 −0.03 −0.45 −0.37 0.74 0.33 −0.32 −0.51 −0.10 −0.79 −0.25 −0.51 −0.64 0.51 ⎞ ⎟ ⎟ ⎟ ⎠ 𝑇 = (1.50 0.17 0.01 0.00) Rank 2 approximation: ̃ 𝑁 = ⎛ ⎜ ⎜ ⎜ ⎝ −0.79 0.58 −0.45 −0.37 −0.32 −0.51 −0.25 −0.51 ⎞ ⎟ ⎟ ⎟ ⎠ (1.50 0.00 0.00 0.17) (−0.79 −0.45 −0.32 −0.25 0.58 −0.37 −0.51 −0.51) = ⎛ ⎜ ⎜ ⎜ ⎝ 1.000 0.501 0.333 0.249 0.501 0.330 0.251 0.203 0.333 0.251 0.200 0.166 0.249 0.203 0.166 0.140 ⎞ ⎟ ⎟ ⎟ ⎠

18

Example

𝑁 = ⎛ ⎜ ⎜ ⎜ ⎝ 1.000 0.500 0.333 0.250 0.500 0.333 0.250 0.200 0.333 0.250 0.200 0.167 0.250 0.200 0.167 0.143 ⎞ ⎟ ⎟ ⎟ ⎠ = 𝑉 diag(𝑇) 𝑊 𝑢 𝑉 = 𝑊 = ⎛ ⎜ ⎜ ⎜ ⎝ −0.79 0.58 −0.18 −0.03 −0.45 −0.37 0.74 0.33 −0.32 −0.51 −0.10 −0.79 −0.25 −0.51 −0.64 0.51 ⎞ ⎟ ⎟ ⎟ ⎠ 𝑇 = (1.50 0.17 0.01 0.00) Rank 2 approximation: ̃ 𝑁 = ⎛ ⎜ ⎜ ⎜ ⎝ −0.79 0.58 −0.45 −0.37 −0.32 −0.51 −0.25 −0.51 ⎞ ⎟ ⎟ ⎟ ⎠ (1.50 0.00 0.00 0.17) (−0.79 −0.45 −0.32 −0.25 0.58 −0.37 −0.51 −0.51) = ⎛ ⎜ ⎜ ⎜ ⎝ 1.000 0.501 0.333 0.249 0.501 0.330 0.251 0.203 0.333 0.251 0.200 0.166 0.249 0.203 0.166 0.140 ⎞ ⎟ ⎟ ⎟ ⎠

18

Approximation Error

We can measure the error of the approximation using the Frobenius norm,

‖𝑁 − ̃ 𝑁‖𝐺 = (

𝑛

∑

𝑗=1 𝑜

∑

𝑘=1

(𝑁𝑗𝑘 − ̃ 𝑁𝑗𝑘)2)

1/2

𝑁 − ̃ 𝑁 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 0.00022 −0.00090 0.00012 0.00077 −0.00090 0.00372 −0.00053 −0.00317 0.00012 −0.00053 0.00013 0.00039 0.00077 −0.00317 0.00039 0.00277 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ ‖𝑁 − ̃ 𝑁‖𝐺 = 0.00674

19

Approximation Error

We can measure the error of the approximation using the Frobenius norm,

‖𝑁 − ̃ 𝑁‖𝐺 = (

𝑛

∑

𝑗=1 𝑜

∑

𝑘=1

(𝑁𝑗𝑘 − ̃ 𝑁𝑗𝑘)2)

1/2

𝑁 − ̃ 𝑁 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 0.00022 −0.00090 0.00012 0.00077 −0.00090 0.00372 −0.00053 −0.00317 0.00012 −0.00053 0.00013 0.00039 0.00077 −0.00317 0.00039 0.00277 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ ‖𝑁 − ̃ 𝑁‖𝐺 = 0.00674

19

Cov Mat - Strong dependence (large eff. range):

10 20 30 40 50 2 4 6 8 10 14

SVD

Singular Values 10 20 30 40 50 5 10 15

Low Rank SVD

Rank Error (Frob. norm)

20

Cov Mat - Weak dependence (short eff. range):

10 20 30 40 50 1 2 3 4

SVD

Singular Values 10 20 30 40 50 2 4 6 8 10

Low Rank SVD

Rank Error (Frob. norm)

21

How does this help? (Sherman-Morrison-Woodbury)

There is an immensely useful linear algebra identity, the Sherman-Morrison-Woodbury formula, for the inverse (and determinant) of a decomposed matrix,

̃ 𝑁

𝑜×𝑛 −1

= ( 𝐵

𝑜×𝑛 + 𝑉 𝑜×𝑙 𝑇 𝑙×𝑙 𝑊 𝑢 𝑙×𝑛) −1

= 𝐵−1 − 𝐵−1𝑉 (𝑇−1 + 𝑊 𝑢𝐵−1𝑉)

−1 𝑊 𝑢𝐵−1.

How does this help?

• Imagine that 𝐵 = diag(𝐵), then it is trivial to find 𝐵−1.
• 𝑇−1 is 𝑙 × 𝑙 which is hopefully small, or even better 𝑇 = diag(𝑇).
• (𝑇−1 + 𝑊 𝑢𝐵−1𝑉) is 𝑙 × 𝑙 which is also hopefully small.

22

How does this help? (Sherman-Morrison-Woodbury)

There is an immensely useful linear algebra identity, the Sherman-Morrison-Woodbury formula, for the inverse (and determinant) of a decomposed matrix,

̃ 𝑁

𝑜×𝑛 −1

= ( 𝐵

𝑜×𝑛 + 𝑉 𝑜×𝑙 𝑇 𝑙×𝑙 𝑊 𝑢 𝑙×𝑛) −1

= 𝐵−1 − 𝐵−1𝑉 (𝑇−1 + 𝑊 𝑢𝐵−1𝑉)

−1 𝑊 𝑢𝐵−1.

How does this help?

• Imagine that 𝐵 = diag(𝐵), then it is trivial to find 𝐵−1.
• 𝑇−1 is 𝑙 × 𝑙 which is hopefully small, or even better 𝑇 = diag(𝑇).
• (𝑇−1 + 𝑊 𝑢𝐵−1𝑉) is 𝑙 × 𝑙 which is also hopefully small.

22

Aside - Determinant

Remember for any MVN distribution when evaluating the likelihood

−1 2 log |Σ| − 1 2(𝐲 − 𝝂)′𝚻−1(𝐲 − 𝝂) − 𝑜 2 log 2𝜌

we need the inverse of Σ as well as its determinant.

• For a full rank Cholesky decomposition we get the determinant for

“free’’.

|𝑁| = |𝑀𝑀𝑢| =

𝑜

∏

𝑗=1

(diag(𝑀)𝑗)

2

• For a low rank approximation the Sherman-Morrison-Woodbury

Determinant lemma gives us,

det( ̃ 𝑁) = det(𝐵 + 𝑉𝑇𝑊 𝑢) = det(𝑇−1 + 𝑊 𝑢𝐵−1𝑉) det(𝑇) det(𝐵)

23

Aside - Determinant

Remember for any MVN distribution when evaluating the likelihood

−1 2 log |Σ| − 1 2(𝐲 − 𝝂)′𝚻−1(𝐲 − 𝝂) − 𝑜 2 log 2𝜌

we need the inverse of Σ as well as its determinant.

• For a full rank Cholesky decomposition we get the determinant for

“free’’.

|𝑁| = |𝑀𝑀𝑢| =

𝑜

∏

𝑗=1

(diag(𝑀)𝑗)

2

• For a low rank approximation the Sherman-Morrison-Woodbury

Determinant lemma gives us,

det( ̃ 𝑁) = det(𝐵 + 𝑉𝑇𝑊 𝑢) = det(𝑇−1 + 𝑊 𝑢𝐵−1𝑉) det(𝑇) det(𝐵)

23

Aside - Determinant

Remember for any MVN distribution when evaluating the likelihood

−1 2 log |Σ| − 1 2(𝐲 − 𝝂)′𝚻−1(𝐲 − 𝝂) − 𝑜 2 log 2𝜌

we need the inverse of Σ as well as its determinant.

• For a full rank Cholesky decomposition we get the determinant for

“free’’.

|𝑁| = |𝑀𝑀𝑢| =

𝑜

∏

𝑗=1

(diag(𝑀)𝑗)

2

• For a low rank approximation the Sherman-Morrison-Woodbury

Determinant lemma gives us,

det( ̃ 𝑁) = det(𝐵 + 𝑉𝑇𝑊 𝑢) = det(𝑇−1 + 𝑊 𝑢𝐵−1𝑉) det(𝑇) det(𝐵)

23

Low rank approximations for GPs

For a standard spatial random effects model,

𝑧(𝐭) = 𝑦(𝐭) 𝜸 + 𝑥(𝐭) + 𝜗, 𝜗 ∼ 𝑂(0, 𝜐2𝐽) 𝑥(𝐭) ∼ 𝒪(0, 𝚻(𝐭)), 𝚻(𝐭, 𝐭′) = 𝜏 𝜍(𝐭, 𝐭′|𝜄)

if we can replace 𝚻(𝐭) with a low rank approximation of the form

• 𝚻(𝐭) ≈ 𝐕 𝐓 𝐖𝑢 where
• 𝐕 and 𝐖 are 𝑜 × 𝑙,
• 𝐓 is 𝑙 × 𝑙, and
• 𝐵 = 𝜐2𝐽 or a similar diagonal matrix

24

Predictive Processes

25

Gaussian Predictive Processes

For a rank 𝑙 approximation,

• Pick 𝑙 knot locations 𝐭⋆
• Calculate knot covariance, 𝚻(𝐭⋆), and knot cross-covariance, 𝚻(𝐭, 𝐭⋆)
• Approximate full covariance using

𝚻(𝐭) ≈ 𝚻(𝐭, 𝐭⋆)

𝑜×𝑙

𝚻(𝐭⋆)−1

𝑙×𝑙

𝚻(𝐭⋆, 𝐭)

𝑙×𝑜

.

• PPs systematically underestimates variance (𝜏2) and inflate 𝜐2, Modified

predictive processs corrects this using

𝚻(𝐭) ≈𝚻(𝐭, 𝐭⋆) 𝚻(𝐭⋆)−1 𝚻(𝐭⋆, 𝐭) + diag(𝚻(𝐭) − 𝚻(𝐭, 𝐭⋆) 𝚻(𝐭⋆)−1 𝚻(𝐭⋆, 𝐭)).

Banerjee, Gelfand, Finley, Sang (2008) Finley, Sang, Banerjee, Gelfand (2008)

26

Gaussian Predictive Processes

For a rank 𝑙 approximation,

• Pick 𝑙 knot locations 𝐭⋆
• Calculate knot covariance, 𝚻(𝐭⋆), and knot cross-covariance, 𝚻(𝐭, 𝐭⋆)
• Approximate full covariance using

𝚻(𝐭) ≈ 𝚻(𝐭, 𝐭⋆)

𝑜×𝑙

𝚻(𝐭⋆)−1

𝑙×𝑙

𝚻(𝐭⋆, 𝐭)

𝑙×𝑜

.

• PPs systematically underestimates variance (𝜏2) and inflate 𝜐2, Modified

predictive processs corrects this using

𝚻(𝐭) ≈𝚻(𝐭, 𝐭⋆) 𝚻(𝐭⋆)−1 𝚻(𝐭⋆, 𝐭) + diag(𝚻(𝐭) − 𝚻(𝐭, 𝐭⋆) 𝚻(𝐭⋆)−1 𝚻(𝐭⋆, 𝐭)).

Banerjee, Gelfand, Finley, Sang (2008) Finley, Sang, Banerjee, Gelfand (2008)

26

Gaussian Predictive Processes

For a rank 𝑙 approximation,

• Pick 𝑙 knot locations 𝐭⋆
• Calculate knot covariance, 𝚻(𝐭⋆), and knot cross-covariance, 𝚻(𝐭, 𝐭⋆)
• Approximate full covariance using

𝚻(𝐭) ≈ 𝚻(𝐭, 𝐭⋆)

𝑜×𝑙

𝚻(𝐭⋆)−1

𝑙×𝑙

𝚻(𝐭⋆, 𝐭)

𝑙×𝑜

.

• PPs systematically underestimates variance (𝜏2) and inflate 𝜐2, Modified

predictive processs corrects this using

𝚻(𝐭) ≈𝚻(𝐭, 𝐭⋆) 𝚻(𝐭⋆)−1 𝚻(𝐭⋆, 𝐭) + diag(𝚻(𝐭) − 𝚻(𝐭, 𝐭⋆) 𝚻(𝐭⋆)−1 𝚻(𝐭⋆, 𝐭)).

Banerjee, Gelfand, Finley, Sang (2008) Finley, Sang, Banerjee, Gelfand (2008)

26

Gaussian Predictive Processes

For a rank 𝑙 approximation,

• Pick 𝑙 knot locations 𝐭⋆
• Calculate knot covariance, 𝚻(𝐭⋆), and knot cross-covariance, 𝚻(𝐭, 𝐭⋆)
• Approximate full covariance using

𝚻(𝐭) ≈ 𝚻(𝐭, 𝐭⋆)

𝑜×𝑙

𝚻(𝐭⋆)−1

𝑙×𝑙

𝚻(𝐭⋆, 𝐭)

𝑙×𝑜

.

• PPs systematically underestimates variance (𝜏2) and inflate 𝜐2, Modified

predictive processs corrects this using

𝚻(𝐭) ≈𝚻(𝐭, 𝐭⋆) 𝚻(𝐭⋆)−1 𝚻(𝐭⋆, 𝐭) + diag(𝚻(𝐭) − 𝚻(𝐭, 𝐭⋆) 𝚻(𝐭⋆)−1 𝚻(𝐭⋆, 𝐭)).

Banerjee, Gelfand, Finley, Sang (2008) Finley, Sang, Banerjee, Gelfand (2008)

26

Example

Below we have a surface generate from a squared exponential Gaussian Process where

{Σ}𝑗𝑘 = 𝜏2 exp (−(𝜚 𝑒)2) + 𝜐2𝐽 𝜏2 = 1 𝜚 = 9 𝜐2 = 0.1

True Surface

−3 −2 −1 1 2 3

Observed Data

−3 −2 −1 1 2

27

Predictive Process Model Results

True Field PP − 5 x 5 knots PP − 10 x 10 knots PP − 15 x 15 knots Full GP

• Mod. PP − 5 x 5 knots
• Mod. PP − 10 x 10 knots
• Mod. PP − 15 x 15 knots

28

Performance

30 40 50 50 100 150 200

time error knots

− 25 100 225

model

Full GP PP

• Mod. PP

29

Parameter Estimates

phi sigma.sq tau.sq 4 5 6 7 8 9 1.2 1.5 1.8 2.1 0.1 0.2 0.3 0.4 0.5

Parameter Value model

Full GP PP

• Mod. PP

True

knots

− 25 100 225

30

Random Projections

31

Low Rank Approximations via Random Projections

• 1. Starting with an matrix

𝐁

𝑛×𝑜.

• 2. Draw a Gaussian random matrix

𝛁

𝑜×𝑙+𝑞.

• 3. Form 𝐙 = 𝐁 𝛁 and compute its QR factorization 𝐙 = 𝐑 𝐒
• 4. Form 𝐂 = 𝐑′ 𝐁.
• 5. Compute the SVD of 𝐂 =

̂ 𝐕 𝐓 𝐖′.

• 6. Form the matrix 𝐕 = 𝐑

̂ 𝐕.

• 7. Form

̃ 𝐁 = 𝐕𝐓𝐖′

Resulting approximation has a bounded expected error,

𝐹|𝐁 − 𝐕𝐓𝐖′‖𝐺 ≤ [1 + 4√𝑙 + 𝑞 𝑞 − 1 √min(𝑛, 𝑜)] 𝜏𝑙+1.

Halko, Martinsson, Tropp (2011)

32

Low Rank Approximations via Random Projections

• 1. Starting with an matrix

𝐁

𝑛×𝑜.

• 2. Draw a Gaussian random matrix

𝛁

𝑜×𝑙+𝑞.

• 3. Form 𝐙 = 𝐁 𝛁 and compute its QR factorization 𝐙 = 𝐑 𝐒
• 4. Form 𝐂 = 𝐑′ 𝐁.
• 5. Compute the SVD of 𝐂 =

̂ 𝐕 𝐓 𝐖′.

• 6. Form the matrix 𝐕 = 𝐑

̂ 𝐕.

• 7. Form

̃ 𝐁 = 𝐕𝐓𝐖′

Resulting approximation has a bounded expected error,

𝐹|𝐁 − 𝐕𝐓𝐖′‖𝐺 ≤ [1 + 4√𝑙 + 𝑞 𝑞 − 1 √min(𝑛, 𝑜)] 𝜏𝑙+1.

Halko, Martinsson, Tropp (2011)

32

Low Rank Approximations via Random Projections

• 1. Starting with an matrix

𝐁

𝑛×𝑜.

• 2. Draw a Gaussian random matrix

𝛁

𝑜×𝑙+𝑞.

• 3. Form 𝐙 = 𝐁 𝛁 and compute its QR factorization 𝐙 = 𝐑 𝐒
• 4. Form 𝐂 = 𝐑′ 𝐁.
• 5. Compute the SVD of 𝐂 =

̂ 𝐕 𝐓 𝐖′.

• 6. Form the matrix 𝐕 = 𝐑

̂ 𝐕.

• 7. Form

̃ 𝐁 = 𝐕𝐓𝐖′

Resulting approximation has a bounded expected error,

𝐹|𝐁 − 𝐕𝐓𝐖′‖𝐺 ≤ [1 + 4√𝑙 + 𝑞 𝑞 − 1 √min(𝑛, 𝑜)] 𝜏𝑙+1.

Halko, Martinsson, Tropp (2011)

32

Low Rank Approximations via Random Projections

• 1. Starting with an matrix

𝐁

𝑛×𝑜.

• 2. Draw a Gaussian random matrix

𝛁

𝑜×𝑙+𝑞.

• 3. Form 𝐙 = 𝐁 𝛁 and compute its QR factorization 𝐙 = 𝐑 𝐒
• 4. Form 𝐂 = 𝐑′ 𝐁.
• 5. Compute the SVD of 𝐂 =

̂ 𝐕 𝐓 𝐖′.

• 6. Form the matrix 𝐕 = 𝐑

̂ 𝐕.

• 7. Form

̃ 𝐁 = 𝐕𝐓𝐖′

Resulting approximation has a bounded expected error,

𝐹|𝐁 − 𝐕𝐓𝐖′‖𝐺 ≤ [1 + 4√𝑙 + 𝑞 𝑞 − 1 √min(𝑛, 𝑜)] 𝜏𝑙+1.

Halko, Martinsson, Tropp (2011)

32

Low Rank Approximations via Random Projections

• 1. Starting with an matrix

𝐁

𝑛×𝑜.

• 2. Draw a Gaussian random matrix

𝛁

𝑜×𝑙+𝑞.

• 3. Form 𝐙 = 𝐁 𝛁 and compute its QR factorization 𝐙 = 𝐑 𝐒
• 4. Form 𝐂 = 𝐑′ 𝐁.
• 5. Compute the SVD of 𝐂 =

̂ 𝐕 𝐓 𝐖′.

• 6. Form the matrix 𝐕 = 𝐑

̂ 𝐕.

• 7. Form

̃ 𝐁 = 𝐕𝐓𝐖′

Resulting approximation has a bounded expected error,

𝐹|𝐁 − 𝐕𝐓𝐖′‖𝐺 ≤ [1 + 4√𝑙 + 𝑞 𝑞 − 1 √min(𝑛, 𝑜)] 𝜏𝑙+1.

Halko, Martinsson, Tropp (2011)

32

Low Rank Approximations via Random Projections

• 1. Starting with an matrix

𝐁

𝑛×𝑜.

• 2. Draw a Gaussian random matrix

𝛁

𝑜×𝑙+𝑞.

• 3. Form 𝐙 = 𝐁 𝛁 and compute its QR factorization 𝐙 = 𝐑 𝐒
• 4. Form 𝐂 = 𝐑′ 𝐁.
• 5. Compute the SVD of 𝐂 =

̂ 𝐕 𝐓 𝐖′.

• 6. Form the matrix 𝐕 = 𝐑

̂ 𝐕.

• 7. Form

̃ 𝐁 = 𝐕𝐓𝐖′

Resulting approximation has a bounded expected error,

𝐹|𝐁 − 𝐕𝐓𝐖′‖𝐺 ≤ [1 + 4√𝑙 + 𝑞 𝑞 − 1 √min(𝑛, 𝑜)] 𝜏𝑙+1.

Halko, Martinsson, Tropp (2011)

32

Low Rank Approximations via Random Projections

• 1. Starting with an matrix

𝐁

𝑛×𝑜.

• 2. Draw a Gaussian random matrix

𝛁

𝑜×𝑙+𝑞.

• 3. Form 𝐙 = 𝐁 𝛁 and compute its QR factorization 𝐙 = 𝐑 𝐒
• 4. Form 𝐂 = 𝐑′ 𝐁.
• 5. Compute the SVD of 𝐂 =

̂ 𝐕 𝐓 𝐖′.

• 6. Form the matrix 𝐕 = 𝐑

̂ 𝐕.

• 7. Form

̃ 𝐁 = 𝐕𝐓𝐖′

Resulting approximation has a bounded expected error,

𝐹|𝐁 − 𝐕𝐓𝐖′‖𝐺 ≤ [1 + 4√𝑙 + 𝑞 𝑞 − 1 √min(𝑛, 𝑜)] 𝜏𝑙+1.

Halko, Martinsson, Tropp (2011)

32

Random Matrix Low Rank Approximations and GPs

PreceThe pding algorithm can be modified slightly to take advantage of the positive definite structure of a covariance matrix.

• 1. Starting with an 𝑜 × 𝑜 covariance matrix 𝐁.
• 2. Draw Gaussian random matrix

𝛁

𝑜×𝑙+𝑞.

• 3. Form 𝐙 = 𝐁 𝛁 and compute its QR factorization 𝐙 = 𝐑 𝐒
• 4. Form the 𝐂 = 𝐑′ 𝐁 𝐑.
• 5. Compute the eigen decomposition of 𝐂 =

̂ 𝐕 𝐓 ̂ 𝐕′.

• 6. Form the matrix 𝐕 = 𝐑

̂ 𝐕.

Once again we have a bound on the error,

𝐹‖𝐁 − 𝐕𝐓𝐕′‖𝐺 ≲ 𝑑 ⋅ 𝜏𝑙+1.

Halko, Martinsson, Tropp (2011), Banerjee, Dunson, Tokdar (2012)

33

Random Matrix Low Rank Approximations and GPs

PreceThe pding algorithm can be modified slightly to take advantage of the positive definite structure of a covariance matrix.

• 1. Starting with an 𝑜 × 𝑜 covariance matrix 𝐁.
• 2. Draw Gaussian random matrix

𝛁

𝑜×𝑙+𝑞.

• 3. Form 𝐙 = 𝐁 𝛁 and compute its QR factorization 𝐙 = 𝐑 𝐒
• 4. Form the 𝐂 = 𝐑′ 𝐁 𝐑.
• 5. Compute the eigen decomposition of 𝐂 =

̂ 𝐕 𝐓 ̂ 𝐕′.

• 6. Form the matrix 𝐕 = 𝐑

̂ 𝐕.

Once again we have a bound on the error,

𝐹‖𝐁 − 𝐕𝐓𝐕′‖𝐺 ≲ 𝑑 ⋅ 𝜏𝑙+1.

Halko, Martinsson, Tropp (2011), Banerjee, Dunson, Tokdar (2012)

33

Random Matrix Low Rank Approximations and GPs

PreceThe pding algorithm can be modified slightly to take advantage of the positive definite structure of a covariance matrix.

• 1. Starting with an 𝑜 × 𝑜 covariance matrix 𝐁.
• 2. Draw Gaussian random matrix

𝛁

𝑜×𝑙+𝑞.

• 3. Form 𝐙 = 𝐁 𝛁 and compute its QR factorization 𝐙 = 𝐑 𝐒
• 4. Form the 𝐂 = 𝐑′ 𝐁 𝐑.
• 5. Compute the eigen decomposition of 𝐂 =

̂ 𝐕 𝐓 ̂ 𝐕′.

• 6. Form the matrix 𝐕 = 𝐑

̂ 𝐕.

Once again we have a bound on the error,

𝐹‖𝐁 − 𝐕𝐓𝐕′‖𝐺 ≲ 𝑑 ⋅ 𝜏𝑙+1.

Halko, Martinsson, Tropp (2011), Banerjee, Dunson, Tokdar (2012)

33

Random Matrix Low Rank Approximations and GPs

PreceThe pding algorithm can be modified slightly to take advantage of the positive definite structure of a covariance matrix.

• 1. Starting with an 𝑜 × 𝑜 covariance matrix 𝐁.
• 2. Draw Gaussian random matrix

𝛁

𝑜×𝑙+𝑞.

• 3. Form 𝐙 = 𝐁 𝛁 and compute its QR factorization 𝐙 = 𝐑 𝐒
• 4. Form the 𝐂 = 𝐑′ 𝐁 𝐑.
• 5. Compute the eigen decomposition of 𝐂 =

̂ 𝐕 𝐓 ̂ 𝐕′.

• 6. Form the matrix 𝐕 = 𝐑

̂ 𝐕.

Once again we have a bound on the error,

𝐹‖𝐁 − 𝐕𝐓𝐕′‖𝐺 ≲ 𝑑 ⋅ 𝜏𝑙+1.

Halko, Martinsson, Tropp (2011), Banerjee, Dunson, Tokdar (2012)

33

Random Matrix Low Rank Approximations and GPs

PreceThe pding algorithm can be modified slightly to take advantage of the positive definite structure of a covariance matrix.

• 1. Starting with an 𝑜 × 𝑜 covariance matrix 𝐁.
• 2. Draw Gaussian random matrix

𝛁

𝑜×𝑙+𝑞.

• 3. Form 𝐙 = 𝐁 𝛁 and compute its QR factorization 𝐙 = 𝐑 𝐒
• 4. Form the 𝐂 = 𝐑′ 𝐁 𝐑.
• 5. Compute the eigen decomposition of 𝐂 =

̂ 𝐕 𝐓 ̂ 𝐕′.

• 6. Form the matrix 𝐕 = 𝐑

̂ 𝐕.

Once again we have a bound on the error,

𝐹‖𝐁 − 𝐕𝐓𝐕′‖𝐺 ≲ 𝑑 ⋅ 𝜏𝑙+1.

Halko, Martinsson, Tropp (2011), Banerjee, Dunson, Tokdar (2012)

33

Random Matrix Low Rank Approximations and GPs

PreceThe pding algorithm can be modified slightly to take advantage of the positive definite structure of a covariance matrix.

• 1. Starting with an 𝑜 × 𝑜 covariance matrix 𝐁.
• 2. Draw Gaussian random matrix

𝛁

𝑜×𝑙+𝑞.

• 3. Form 𝐙 = 𝐁 𝛁 and compute its QR factorization 𝐙 = 𝐑 𝐒
• 4. Form the 𝐂 = 𝐑′ 𝐁 𝐑.
• 5. Compute the eigen decomposition of 𝐂 =

̂ 𝐕 𝐓 ̂ 𝐕′.

• 6. Form the matrix 𝐕 = 𝐑

̂ 𝐕.

Once again we have a bound on the error,

𝐹‖𝐁 − 𝐕𝐓𝐕′‖𝐺 ≲ 𝑑 ⋅ 𝜏𝑙+1.

Halko, Martinsson, Tropp (2011), Banerjee, Dunson, Tokdar (2012)

33

Random Matrix Low Rank Approximations and GPs

PreceThe pding algorithm can be modified slightly to take advantage of the positive definite structure of a covariance matrix.

• 1. Starting with an 𝑜 × 𝑜 covariance matrix 𝐁.
• 2. Draw Gaussian random matrix

𝛁

𝑜×𝑙+𝑞.

• 3. Form 𝐙 = 𝐁 𝛁 and compute its QR factorization 𝐙 = 𝐑 𝐒
• 4. Form the 𝐂 = 𝐑′ 𝐁 𝐑.
• 5. Compute the eigen decomposition of 𝐂 =

̂ 𝐕 𝐓 ̂ 𝐕′.

• 6. Form the matrix 𝐕 = 𝐑

̂ 𝐕.

Once again we have a bound on the error,

𝐹‖𝐁 − 𝐕𝐓𝐕′‖𝐺 ≲ 𝑑 ⋅ 𝜏𝑙+1.

Halko, Martinsson, Tropp (2011), Banerjee, Dunson, Tokdar (2012)

33

Low Rank Approximations and GPUs

Both predictive process and random matrix low rank approximations are good candidates for acceleration using GPUs.

• Both use Sherman-Woodbury-Morrison to calculate the inverse

(involves matrix multiplication, addition, and a small matrix inverse).

• Predictive processes involves several covariance matrix calculations

(knots and cross-covariance) and a small matrix inverse.

• Random matrix low rank approximations involves a large matrix

multiplication (𝐁 𝛁) and several small matrix decompositions (QR, eigen).

34

Comparison (𝑜 = 15, 000, 𝑙 = {100, … , 4900})

5 10 15 10 20 30 10 20 30 10 20 30

time time error error method

lr1 lr1 mod pp pp mod

method

lr1 lr1 mod pp pp mod

Strong Dependence Weak Dependence

35

• Rand. Projection LR Depositions for Prediction

This approach can also be used for prediction, if we want to sample

𝐳 ∼ 𝒪(0, 𝚻) Σ ≈ 𝐕𝐓𝐕𝑢 = (𝐕𝐓1/2𝐕𝑢)(𝐕𝐓1/2𝐕𝑢)𝑢

then

𝑧pred = (𝐕 𝐓1/2 𝐕𝑢) × 𝐚 where 𝑎𝑗 ∼ 𝒪(0, 1)

because 𝐕𝑢 𝐕 = 𝐽 since 𝐕 is an orthogonal matrix.

Dehdari, Deutsch (2012)

36

𝑜 = 1000, 𝑞 = 10000

37