Fast Bayesian inference of block Nearest Neighbor Gaussian models for large data

This paper presents the development of a spatial block-Nearest Neighbor Gaussian process (blockNNGP) for location-referenced large spatial data. The key idea behind this approach is to divide the spatial domain into several blocks which are dependent under some constraints. The cross-blocks capture the large-scale spatial dependence, while each block captures the small-scale spatial dependence. The resulting blockNNGP enjoys Markov properties reflected on its sparse precision matrix. It is embedded as a prior within the class of latent Gaussian models, thus fast Bayesian inference is obtained using the integrated nested Laplace approximation. The performance of the blockNNGP is illustrated on simulated examples, a comparison of our approach with other methods for analyzing large spatial data and applications with Gaussian and non-Gaussian real data.

The author would like to thank the referee for the suggestions that we believe drastically improved the context of the paper. Also, Zaida C. Quiroz would like to thank the Pontificia Universidad Católica del Perú for partial financial support through the project [DGI-2019-740]. Marcos O. Prates would like to acknowledge partial funding support from Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) grants 436948/2018-4 and 307547/2018-4 and Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG) grant APQ-01837-22.

Authors and Affiliations

1. Department of Sciences, Pontificia Universidad Católica del Perú, Lima, Peru

   Zaida C. Quiroz

2. Department of Statistics, Universidade Federal de Minas Gerais, Belo Horizonte, Brazil

   Marcos O. Prates

3. Department of Statistics, University of Connecticut, Storrs, USA

   Dipak K. Dey

4. Computer, Electrical and Mathematical Science and Engineering, King Abdullah University of Science and Technology, Thuwal, Saudi Arabia

   H.åvard Rue

Corresponding author

Correspondence to Zaida C. Quiroz.

Conflict of interest

The authors declare that they have no conflict of interest.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Apendix: Proof of proposition 2

Matrix Analysis Background

Theorem A1: A matrix \({\varvec{B}} \in {\varvec{\Re }}^{m\times n}\) is full column rank if and only if \({\varvec{B^{T}B}}\) is invertible

Theorem A2: The determinant of an \(n\times n\) matrix \({\varvec{B}}\) is 0 if and only if the matrix \({\varvec{B}}\) is not invertible.

Theorem A3: Let \({\varvec{T_n}}\) be a triangular matrix (either upper or lower) of order n. Let \(|{\varvec{T_n}}|\) be the determinant of \({\varvec{T_n}}\). Then \(|{\varvec{T_n}}|\) is equal to the product of all the diagonal elements of \({\varvec{T_n}},\) that is, \(|{\varvec{T_n}}| = \prod _{k=1}^{n}(a_{kk}).\)

Proposition A1: If \({\varvec{B}}\) is positive definite (p.d.), then if \({\varvec{S}}\) has full column rank, then \({\varvec{S^{T}BS}}\) is positive definite.

Corollary A1: If \({\varvec{B}}\) is positive definite, then \({\varvec{B^{-1}}}\) is positive definite.

Proof

Without loss of generality, assume that the data were reordered by blocks. From known properties of Gaussian distributions,

where \({\varvec{B_{b_k}}} = {\varvec{C_{b_k, N(b_k)}}} {\varvec{C^{-1}_{N(b_k)}}}\) and \({\varvec{F_{b_k}}} = {\varvec{C_{b_k}}} - {\varvec{C_{b_k, N(b_k)}}} {\varvec{C^{-1}_{N(b_k)}}}{\varvec{C_{ N(b_k),b_k}}}.\) Hence,

Let \({\varvec{w_{b_k}}} - {\varvec{B_{b_k}}} {\varvec{w_{N(b_k)}}} = {\varvec{B^\star _{b_k}}} {\varvec{w_S}}\), and j be the j-th observation of block \({\varvec{b_k}}\), then \(\forall k = 1, \dots , M\), \(i = 1, \dots , n\) and \( j = 1, \dots , n_{bk}\):

and

where \({\varvec{B^{\star }_{b_k} (j)}}= (B^{\star }_{b_k} (j,1),B^{\star }_{b_k} (j,2),\dots , B^{\star }_{b_k} (j,n))\). From these definitions, \( {\varvec{B^{\star }_{b_k}}}\) is a matrix with i-th column full of zeros if \(s_i\notin {\varvec{s_{b_k}}}\) or \( s_i\notin {\varvec{N(s_{b_k})}}\). Since the data were reordered by blocks and the neighbor blocks are from the past, \({\varvec{B^{\star }_{b_k}}}\) has also the next form:

where \({\varvec{A_k}}\) is a \(n_{bk}\times n_{bk}\) matrix and \({\varvec{R_k}}\) is a \(n_{bk}\times \sum _{r=1}^{k-1} n_{br}\) matrix with at least one column with none null-element if \(nb \ne 0\).

Then,

Let \(\sum _{k=1}^M ({\varvec{B_{b_k}}}^\star )^T {\varvec{F^{-1}_{b_k}}} {\varvec{B_{b_k}}}^{\star } = ({\varvec{B^{\star }_s}})^T {\varvec{F_s ^{-1}}} {\varvec{B^{\star }_s}} \), where \({\varvec{B_{s}}} = \left[ \begin{array}{ccccc} {\varvec{B_{b1}^{\star }}} &{} \ldots &{} \ldots &{} \dots &{} {\varvec{B_{bM}^{\star }}} \\ \end{array}\right] \) and \({\varvec{F_s ^{-1}}} = \text{ diag }( {\varvec{F_{b_k} ^{-1}}})\). \( {\varvec{F_s ^{-1}}} \) is a block diagonal matrix And given that \( {\varvec{B^{\star }_{b_k}}}\) is a matrix with i-th column full of zeros for \(i> \sum _{r=1}^{k} n_{br}\), then \({\varvec{B_s}}\) is a block matrix and lower triangular.

Then, \(\widetilde{p}({\varvec{w_S}}) \propto \frac{1}{ \prod _{k=1}^M |{\varvec{F_{b_k}}}|^{1/2}} \exp \left\{ -\frac{1}{2} {\varvec{w_S^T}} ( {\varvec{B_s^T F_s ^{-1}B_s}}) {\varvec{w_S}} \right\} \) and \( {\varvec{{\widetilde{C}}_s^{-1}}} = {\varvec{{\widetilde{Q}}_s}}= {\varvec{B_s^T F_s ^{-1}B_s}}.\)

We also need to prove that \({\varvec{\widetilde{Q}_s}}\) is positive definite. From properties of the Normal distribution, the covariance of the conditional distribution of \({\varvec{w_{b_k}}}|{\varvec{w_{N(b_k)}}}\) is also p.d. (by Schur complement conditions), then

is p.d. Moreover, \({\varvec{F_s}} = diag({\varvec{F_{bk}}})\) and a block diagonal matrix is p.d. iff each diagonal block is positive definite, so given that \({\varvec{F_{bk}}}\) is p.d. and \({\varvec{F_s}}\) is block diagonal with blocks \({\varvec{F_{bk}}}\) p.d then \({\varvec{F_s}}\) is p.d. By Corollary A1, \({\varvec{F_s}}\) is p.d. then \({\varvec{F^{-1}_s}}\) is p.d. By Theorem A1, \({\varvec{B_s}}\) has full column rank if and only iff \({\varvec{R_s}} = {\varvec{B^T_s B_s}}\) is invertible. By Theorem A2, the inverse of \({\varvec{R_s}}\) exists if and only if \(|{\varvec{R_s}}|\ne 0\). Using the well-known matrix theorems, we can prove the following: \(|{\varvec{R_s}}|= |{\varvec{B^T_s B_s}}| = |{\varvec{B^T_s}}| |{\varvec{B_s}}| \ne 0\) if \(|{\varvec{B^T_s}}| = |{\varvec{B_s}}| \ne 0\). Given that \({\varvec{B_s}}\) is a lower triangular matrix, by Theorem A3, \(|{\varvec{B_s}}| = \prod _{k=1}^{n}(b_{kk})\). And, \(b_{kk} =1, \forall k\), then \(|{\varvec{B_s}}| \ne 0.\) So, \({\varvec{R_s}}\) is invertible and \({\varvec{B_s}}\) has full column rank. By Proposition A1, given that \({\varvec{B_s}}\) has full column rank, and \({\varvec{F_s^{-1}}}\) is p.d. then \({\varvec{\widetilde{Q}_s}} = {\varvec{\widetilde{C}_s^{-1}}} = {\varvec{B^T_s F_s^{-1} B_s}}\) is p.d.

Also by corollary A1, \({\varvec{\widetilde{C}_s^{-1}}}\) is p.d. then \({\varvec{\widetilde{C}_s}}\) is p.d. Finally, since \({\widetilde{p}}({\varvec{w_S}}) \propto \frac{1}{ \prod _{k=1}^M |{\varvec{F_{b_k}}}|^{1/2}} \exp \left\{ -\frac{1}{2} {\varvec{w_S}}^T ( {\varvec{{\widetilde{C}}_s}}^{-1}) {\varvec{w_S}} \right\} \), \({\varvec{{\widetilde{C}}_s}}^{-1} = {\varvec{B_s^T F_s ^{-1}B_s}}\), and \({\varvec{\widetilde{C}_s}}\) is p.d., then \({\widetilde{p}}({\varvec{w_S}})\) is a pdf of a multivariate normal distribution. \(\square \)

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

Reprints and permissions