Matrixfree conditional simulations of GMRF Somak Dutta Joint Work - PowerPoint PPT Presentation

Matrix–free conditional simulations of GMRF Somak Dutta Joint Work with Debashis Mondal. University of Chicago June 24, 2014

Data on a regular grid. Image of an dummy array of plot D 1 . Black = missing observations.

Mixed linear model. y = T τ + Fx + ǫ . Array dimension = r × c . (VERY LARGE). y = n × 1 response vector. τ = m × 1 vector of fixed effects. T = n × m known design matrix. x = rc × 1 vector of underlying spatial random field. F = known sparse incidence matrix - Fx gives back the values of the spatial random field on n observed plots. ǫ ∼ N ( 0 , λ − 1 y I n ) : nugget effects.

Intrinsic auto-regression model for x . y = T τ + Fx + ǫ . ◮ x is Gaussian with sparse singular precision matrix W , � � � � ( x i , j − x i − 1 , j ) 2 + λ 01 ( x i , j − x i , j − 1 ) 2 . x T Wx = λ 10 ◮ W has analytically known spectral decomposition W = P ( λ 01 D 01 + λ 10 D 10 ) P T . ◮ P correspond to the two dimensional discrete cosine transformation.

Conditional simulation. Interested in sampling from: x | y ∼ N ( λ y A − 1 F T ( y − T τ ) , A − 1 ) , A = λ y F T F + W . ◮ Step 1: First draw z ∼ N ( 0 , I ). ◮ Traditional way: Compute Cholesky factor L such that LL T = A . And let x = L − 1 z . ◮ Costs: memory = O (( rc ) 1 . 5 ) , #FLOPs = O (( rc ) 2 ) . ◮ We will create algorithm that has costs: memory = O ( rc ) , #FLOPs = O ( rc log rc )

An “exact” method ◮ x | y ∼ N ( A − 1 ( y − T τ ) , A − 1 ) , A = λ y F T F + W ◮ Analytically known spectral decomposition: W = PDP T . ◮ Square root of A : 1 then SS T = A PD 1 / 2 ] y F T S = [ λ 2 Simulation algorithm: ◮ Strike 1: First draw z ∼ N ( 0 , I ). ◮ Strike 2: Sample with A as the covariance matrix b = Sz + λ y ( y − T τ ) ∼ N ( λ y ( y − T τ ) , A ) ◮ Strike 3: Solve x = A − 1 b ∼ N ( A − 1 ( y − T τ ) , A − 1 )

Lanczos algorithm and Incomplete Cholesky Preconditioner To solve: Ax = b using Lanczos algorithm (Dutta and Mondal, 2012). ◮ Condition number of A → ∞ . ◮ L = incomplete Cholesky factorization (lower triangular): LL T ≈ A ⇒ L − 1 AL − T ≈ I . ◮ solve L − 1 AL − T x 1 = L − 1 b , then x = L − T x 1 . ◮ Geometric convergence of Lanczos algo in O (log rc ) iterations.

Arsenic concentration in Bangladesh (Dutta and Mondal, 2013). Arsenic conc. (in ppb) ● ● ● 0 − 0.5 ● ● ● ● ● ● ● ● 0.5 − 10 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 10 − 50 ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 50 − 150 ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 150 − 1660 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 1027 1020 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● 606 ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● 463 ● ● ● ● ● ●● ● ● 418 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 6 ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● log−arsenic ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 4 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 2 ● ● ● ● ● ● ● ● ● ● ● ● ● 0 ● ● ● Embed the data in a 500 x 300 array.