---
title: "Circulant Embedding"
author: "Luke Tierney"
output:
  html_document
---

A Gaussian autocorrelation vector is created by

```{r}
f1 <- function(n, b) exp(- b * (0 : (n - 1)) ^ 2)
```

A circulant version is created by

```{r}
f2 <- function(n, b) exp(- b * pmin(0 : (n - 1), n : 1) ^ 2)
```

```{r}
plot(f1(16, 0.5), type = "l")
lines(f2(16, 0.5), col = "red")
```

The covariance matrix corresponding to `f1` can be created using the
`toeplitz` function:

```{r}
m1 <- toeplitz(f1(5, 0.5))
m1
```

This is identical to the upper left corner of the matrix for `f2(10, 0)`

```{r}
m2 <- toeplitz(f2(10, 0.5))
identical(m1, m2[1 : 5, 1 : 5])
```

A normal vector with covariance matrix corresponding to `f1(n, b)` can be
generated with

```{r}
gen1 <- function(n, b) {
    a <- f1(n, b)
    z <- rnorm(n)
    chol(toeplitz(a)) %*% z
}
```

An alternative is to use the larger circulant covariance matrix
corresponding for `f2(2 * n, b)`, keeping only the first half:

```{r}
gen2 <- function(n, b) {
    N <- 2 * n
    a <- f2(N, b)
    z <- rnorm(N)
    x <- Re(fft(sqrt(fft(a)) * fft(z), inverse = TRUE)) / N
    x[1 : n]
}
```

For larger sample sizes the circulant approach is much faster.

```{r}
n <- 2048
system.time(gen1(n, 0.5))
system.time(gen2(n, 0.5))
```

If `b` is constant then the Cholesky factor and the `sqrt(fft(a))`
vector can be pre-computed and re-used. This reduces but does not
eliminate the advantage of the circulant approach.

Usually the circulant embedding approach only extends the region of interest by 10% or 20%,
which is often enough for a good approximation.

Circulant embedding is particularly effective for spatial and volume data.