Sampling a zero mean GMRF with block-circulant precision
Henrique Ap. Laureano
email: henrique.laureano@kaust.edu.sa
## [1] "Mon Jul 31 10:50:09 2017"First, generating a block-circulant precision matrix
options(width=125)
cm <- function(x) { # cm: circulant matrix
n <- length(x)
suppressWarnings(
matrix(x[matrix(1:n, n+1, n+1, byrow = TRUE)[c(1, n:2), 1:n]], n, n)
)
} ; cm(letters[1:4])## [,1] [,2] [,3] [,4]
## [1,] "a" "b" "c" "d"
## [2,] "d" "a" "b" "c"
## [3,] "c" "d" "a" "b"
## [4,] "b" "c" "d" "a"cols <- function(ms, n) { # ms: matrices (list of)
if (n == 0) return(ms)
c( tail(ms, n), head(ms, -n) )
}
rcols <- function(n, ms) do.call(rbind, cols(ms, n))
bcm <- function(ms) { # bcm: Block-Circulant Matrix
n <- length(ms)
do.call( cbind, lapply(0:(n-1), rcols, ms) )
}
### Three different circulating matrix of dimension 3 x 3
( ms <- list( cm(c(5, 2, 2)), cm(c(2, 0, 0)), cm(1:3) ) )## [[1]]
## [,1] [,2] [,3]
## [1,] 5 2 2
## [2,] 2 5 2
## [3,] 2 2 5
##
## [[2]]
## [,1] [,2] [,3]
## [1,] 2 0 0
## [2,] 0 2 0
## [3,] 0 0 2
##
## [[3]]
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 3 1 2
## [3,] 2 3 1( th <- bcm(ms) )## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 5 2 2 1 2 3 2 0 0
## [2,] 2 5 2 3 1 2 0 2 0
## [3,] 2 2 5 2 3 1 0 0 2
## [4,] 2 0 0 5 2 2 1 2 3
## [5,] 0 2 0 2 5 2 3 1 2
## [6,] 0 0 2 2 2 5 2 3 1
## [7,] 1 2 3 2 0 0 5 2 2
## [8,] 3 1 2 0 2 0 2 5 2
## [9,] 2 3 1 0 0 2 2 2 5chol(th) # Showing that the resulting matrix is SPD## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 2.236068 0.8944272 0.8944272 0.4472136 0.89442719 1.3416408 0.89442719 0.0000000 0.00000000
## [2,] 0.000000 2.0493902 0.5855400 1.2686701 0.09759001 0.3903600 -0.39036003 0.9759001 0.00000000
## [3,] 0.000000 0.0000000 1.9639610 0.4364358 1.09108945 -0.2182179 -0.29095719 -0.2909572 1.01835015
## [4,] 0.000000 0.0000000 0.0000000 1.7320508 0.57735027 0.5773503 0.70565033 0.5132002 1.47545069
## [5,] 0.000000 0.0000000 0.0000000 0.0000000 1.63299316 0.4082483 1.31546671 0.5670115 0.02268046
## [6,] 0.000000 0.0000000 0.0000000 0.0000000 0.00000000 1.5811388 -0.03513642 1.2824793 0.22838672
## [7,] 0.000000 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 1.31656118 0.9367839 0.93678391
## [8,] 0.000000 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.9250780 0.38458297
## [9,] 0.000000 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.84134723Implementation of the algorithm steps:
1: Sample \(\mathbf{z}\), where \({\rm Re}(z_{ij}) \overset{{\rm iid}}{\sim} N(0, 1)\) and \({\rm Im}(z_{ij}) \overset{{\rm iid}}{\sim} N(0, 1)\)
round(z <- matrix(rnorm(nrow(th)) + rnorm(nrow(th)) * 1i
, nrow = nrow(th)
, ncol = ncol(th)), 1)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] -0.5-1.2i -0.5-1.2i -0.5-1.2i -0.5-1.2i -0.5-1.2i -0.5-1.2i -0.5-1.2i -0.5-1.2i -0.5-1.2i
## [2,] 0.1-0.9i 0.1-0.9i 0.1-0.9i 0.1-0.9i 0.1-0.9i 0.1-0.9i 0.1-0.9i 0.1-0.9i 0.1-0.9i
## [3,] 1.1-0.1i 1.1-0.1i 1.1-0.1i 1.1-0.1i 1.1-0.1i 1.1-0.1i 1.1-0.1i 1.1-0.1i 1.1-0.1i
## [4,] 0.2-1.2i 0.2-1.2i 0.2-1.2i 0.2-1.2i 0.2-1.2i 0.2-1.2i 0.2-1.2i 0.2-1.2i 0.2-1.2i
## [5,] -0.2+0.0i -0.2+0.0i -0.2+0.0i -0.2+0.0i -0.2+0.0i -0.2+0.0i -0.2+0.0i -0.2+0.0i -0.2+0.0i
## [6,] -0.2-1.6i -0.2-1.6i -0.2-1.6i -0.2-1.6i -0.2-1.6i -0.2-1.6i -0.2-1.6i -0.2-1.6i -0.2-1.6i
## [7,] -0.7-0.3i -0.7-0.3i -0.7-0.3i -0.7-0.3i -0.7-0.3i -0.7-0.3i -0.7-0.3i -0.7-0.3i -0.7-0.3i
## [8,] -1.7-1.6i -1.7-1.6i -1.7-1.6i -1.7-1.6i -1.7-1.6i -1.7-1.6i -1.7-1.6i -1.7-1.6i -1.7-1.6i
## [9,] 0.1+1.8i 0.1+1.8i 0.1+1.8i 0.1+1.8i 0.1+1.8i 0.1+1.8i 0.1+1.8i 0.1+1.8i 0.1+1.8i2: Compute the (real) eigenvalues, \(\mathbf{\Lambda} = \sqrt{n\cdot N}{\rm DFT2}(\theta)\)
On what
- \(n\) is the number of rows in the matrix \(\theta\);
- \(N\) is the number of columns in the matriz \(\theta\);
- \(\theta\) I assumed it is the block-circulant precision matrix;
- DFT2 (two-dimensional discrete Fourier transform matrix):
- matrix with elements:
\[ \frac{1}{nN} \sum_{{i}'=0}^{n-1} \sum_{{j}'=0}^{N-1} \theta_{{i}'{j}'} \exp \bigg( -2 \pi {\rm i} (\frac{i{i}'}{n} + \frac{j{j}'}{N}) \bigg), \quad {\rm i} = \sqrt{-1}. \]
By steps, for understanding:
For a given \(i\) e \(j\), 1 and 1, for example, we calculated \(\sum_{{i}'=0}^{n-1} \sum_{{j}'=0}^{N-1} \theta_{{i}'{j}'}\) with
# with the following line of code we select all the necessary elements of
# the matrix mx and order them in a way in which the correct
# multiplication is done with the exponential of the next part of the
# equation
c(t( mx[1:(n-1), 1:(N-1)] )) # in sum()For a given \(i\) e \(j\), 1 and 1, for example, we calculated \(\sum_{{i}'=0}^{n-1} \sum_{{j}'=0}^{N-1} \exp \bigg( -2 \pi {\rm i} (\frac{i{i}'}{n} + \frac{j{j}'}{N}) \bigg)\) with
# in sum()
exp(-2 * pi * 1i * apply(expand.grid(i*1:(n-1)/n, j*1:(N-1)/N), 1, sum))
\(\square\)
Ok, now, the complete implementation:
dft2 <- function(mx) {
n <- dim(mx)[1]
N <- dim(mx)[2]
mx.dft2 <- matrix(NA, nrow = n, ncol = N)
for (i in 1:n) {
for (j in 1:N) {
mx.dft2[i, j] <- ( 1/sqrt(n*N) ) *
sum( c(t( mx[1:(n-1), 1:(N-1)] )) *
exp( -2 * pi * 1i *
apply(expand.grid(i*1:(n-1)/n, j*1:(N-1)/N), 1, sum)
) )
} }
return(mx.dft2) } ; n <- dim(th)[1] ; N <- dim(th)[2]
round( lambda <- sqrt(nrow(th) * ncol(th)) * dft2(th), 1 )## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] -3.1- 1.9i 7.1+14.6i -5.9+1.0i -6.5+ 0.6i 5.1- 5.2i -5.9-0.7i -8.2+ 4.3i 36.8-12.6i -19.4+0.1i
## [2,] -19.0- 4.8i -0.3- 3.0i 2.0+2.1i -3.8+12.7i 1.5+ 0.7i 2.0+0.3i 24.3-12.6i 4.8+ 3.3i -11.5+1.2i
## [3,] 0.8- 2.9i -4.3- 5.0i -2.0+0.0i -2.5- 0.4i -2.5- 1.3i 29.5+6.1i -4.3+ 3.3i 0.8+ 1.2i -15.5-0.9i
## [4,] 0.7- 3.6i 7.9- 4.8i -2.1-0.6i -2.6- 1.0i 27.3+18.3i -2.1-2.4i -4.4+ 2.7i -9.3- 7.1i -15.6-1.5i
## [5,] -9.3+ 7.1i -4.4- 2.7i -2.1+2.4i 27.3-18.3i -2.6+ 1.0i -2.1+0.6i 7.9+ 4.8i 0.7+ 3.6i -15.6+1.5i
## [6,] 0.8- 1.2i -4.3- 3.3i 29.5-6.1i -2.5+ 1.3i -2.5+ 0.4i -2.0+0.0i -4.3+ 5.0i 0.8+ 2.9i -15.5+0.9i
## [7,] 4.8- 3.3i 24.3+12.6i 2.0-0.3i 1.5- 0.7i -3.8-12.7i 2.0-2.1i -0.3+ 3.0i -19.0+ 4.8i -11.5-1.2i
## [8,] 36.8+12.6i -8.2- 4.3i -5.9+0.7i 5.1+ 5.2i -6.5- 0.6i -5.9-1.0i 7.1-14.6i -3.1+ 1.9i -19.4-0.1i
## [9,] -12.7- 2.0i -17.8- 4.2i -15.5+0.9i -16.0+ 0.5i -16.0- 0.5i -15.5-0.9i -17.8+ 4.2i -12.7+ 2.0i 124.0+0.0iProblem / doubt
Via spectral decomposition, for example, we should also obtain these eigenvalues, \(\mathbf{\Lambda}\), but this is not what happens:
round(eigen(th)$values, 1)## [1] 17.0+0.0i 5.0+3.5i 5.0-3.5i 3.5+2.6i 3.5-2.6i 2.0+3.5i 2.0-3.5i 3.5+0.9i 3.5-0.9iThe “real” eigenvalues, written in this step of the algorithm also caused me doubt. What does that mean? It’s not just to get the real part of \(\mathbf{\Lambda}\), in this way step 3 fails.
Looking to the implementation, I don’t see anything wrong (wrong implemented), and you?
3: \(\boldsymbol{\upsilon}= {\rm DFT2}((\mathbf{\Lambda}\text{elementwisepower}(-\frac{1}{2})) \odot\mathbf{z})\)
round(upsilon <- dft2( lambda**(-1/2)*z ), 1)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] -0.1+0.2i -0.5+0.0i 0.1-0.5i -0.1-0.2i -0.4+0.4i 0.5-0.6i -0.2-0.1i 0.1+0.9i 0.4-0.2i
## [2,] 0.2-0.7i -0.2+0.2i 0.1-0.1i 0.6+0.0i 0.2+0.9i 0.0+0.3i 0.4+0.1i -0.5-0.1i -0.7-0.7i
## [3,] -0.2+0.1i -0.1-0.3i -0.1-0.1i -0.2+0.1i -0.1-0.6i 0.0+0.0i 0.2+0.2i 0.1+0.0i 0.4+0.5i
## [4,] 0.3+0.0i 0.4+0.0i 0.2+0.8i -0.3-0.1i 0.4+0.1i -0.5-0.2i -0.6-0.5i 0.3-0.6i -0.2+0.4i
## [5,] 0.3+0.3i 0.2-0.2i -0.1+0.4i 0.2+0.2i 0.0+0.1i -0.3-0.3i -0.3-0.2i 0.3-0.2i -0.4-0.2i
## [6,] -0.2+0.4i 0.2+0.1i 0.2+0.0i -0.4+0.2i -0.1-0.5i 0.1+0.2i 0.0+0.3i -0.2-0.7i 0.3+0.0i
## [7,] 0.1+0.2i 0.3+0.1i -0.7-0.2i -0.5-0.8i 0.1-0.5i -0.2+0.2i 0.2-0.4i 0.5+0.2i 0.3+1.1i
## [8,] -0.1+0.4i -0.2-0.6i 0.0-0.6i -0.2+0.4i 0.5+0.3i 0.2-0.2i -0.2+0.5i 0.1-0.2i 0.0+0.0i
## [9,] -0.3-0.9i 0.0+0.7i 0.3+0.2i 0.9+0.2i -0.8-0.1i 0.1+0.4i 0.6+0.0i -0.8+0.6i 0.0-1.0i4: \(\mathbf{x} = {\rm Re}(\boldsymbol{\upsilon})\)
x <- Re(upsilon)5: Return \(\mathbf{x}\)
round(x, 4)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] -0.0666 -0.5033 0.0911 -0.0520 -0.3605 0.5393 -0.2205 0.1243 0.4484
## [2,] 0.1966 -0.1977 0.1042 0.6213 0.2273 -0.0449 0.3624 -0.5292 -0.7399
## [3,] -0.1978 -0.1319 -0.0919 -0.1983 -0.0758 0.0333 0.1641 0.1340 0.3644
## [4,] 0.2904 0.3964 0.1630 -0.3311 0.3637 -0.4703 -0.5530 0.3491 -0.2081
## [5,] 0.3150 0.1734 -0.0723 0.2306 0.0125 -0.2525 -0.3178 0.2763 -0.3653
## [6,] -0.1510 0.2461 0.2498 -0.4086 -0.0673 0.1275 -0.0474 -0.2126 0.2635
## [7,] 0.0736 0.2579 -0.7367 -0.4987 0.1475 -0.2293 0.1848 0.5329 0.2678
## [8,] -0.1384 -0.2011 -0.0268 -0.2220 0.5211 0.1536 -0.1811 0.1004 -0.0057
## [9,] -0.3217 -0.0397 0.3196 0.8588 -0.7685 0.1432 0.6085 -0.7752 -0.0251
\(\blacksquare\)