I'm trying to generate many draws (i.e., realizations) of a Gaussian process $e_i(t)$, $1\leq t \leq T$ with mean 0 and covariance function $\gamma(s,t)=\exp(-|t-s|)$.

Is there an efficient way to do this that wouldn't involve computing the square root of a $T \times T$ covariance matrix? Alternatively can anyone recommend an R package to do this?

  • 2
    $\begingroup$ It's a stationary process (looks close to a simple version of an OU process). Is it uniformly sampled? $\endgroup$
    – cardinal
    May 24, 2012 at 13:53
  • $\begingroup$ The R package mvtnorm has rmvnorm(n, mean, sigma) where sigma is the covariance matrix; you'd have to construct the covariance matrix for your sampled / selected $t$s yourself, though. $\endgroup$
    – jbowman
    May 24, 2012 at 14:06
  • 3
    $\begingroup$ @jb Presumably $T$ is huge, otherwise the OP wouldn't be asking to avoid the matrix decomposition (which is implicit in rmvnorm). $\endgroup$
    – whuber
    May 24, 2012 at 14:31
  • 1
    $\begingroup$ @cardinal I agree, this is a Ornstein-Uhlenbeck Gaussian process. (It would be great if the "Ornstein Uhlenbeck" keyword could be edited into the question and/or title. It would get this question the more traffic it deserves) $\endgroup$
    – redmoskito
    Dec 5, 2014 at 23:12

2 Answers 2


Yes. There is a very efficient (linear time) algorithm, and the intuition for it comes directly from the uniformly-sampled case.

Suppose we have a partition of $[0,T]$ such that $0=t_0 < t_1 < t_2 < \cdots < t_n = T$.

Uniformly sampled case

In this case we have $t_i = i \Delta$ where $\Delta = T/n$. Let $X_i := X(t_i)$ denote the value of the discretely sampled process at time $t_i$.

It is easy to see that the $X_i$ form an AR(1) process with correlation $\rho = \exp(-\Delta)$. Hence, we can generate a sample path $\{X_t\}$ for the partition as follows $$ X_{i+1} = \rho X_i + \sqrt{1-\rho^2} Z_{i+1} \>, $$ where $Z_i$ are iid $\mathcal N(0,1)$ and $X_0 = Z_0$.

General case

We might then imagine that it could be possible to do this for a general partition. In particular, let $\Delta_i = t_{i+1} - t_i$ and $\rho_i = \exp(-\Delta_i)$. We have that $$ \gamma(t_i,t_{i+1}) = \rho_i \>, $$ and so we might guess that $$ X_{i+1} = \rho_i X_i + \sqrt{1-\rho_i^2} Z_{i+1} \>. $$

Indeed, $\mathbb E X_{i+1} X_i = \rho_i$ and so we at least have the correlation with the neighboring term correct.

The result now follows by telescoping via the tower property of conditional expectation. Namely, $$ \newcommand{\e}{\mathbb E} \e X_i X_{i-\ell} = \e( \e(X_i X_{i-\ell} \mid X_{i-1} )) = \rho_{i-1} \mathbb E X_{i-1} X_{i-\ell} = \cdots = \prod_{k=1}^\ell \rho_{i-k} \>, $$ and the product telescopes in the following way $$ \prod_{k=1}^\ell \rho_{i-k} = \exp\Big(-\sum_{k=1}^\ell \Delta_{i-k}\Big) = \exp(t_{i-\ell} - t_i) = \gamma(t_{i-\ell},t_i) \>. $$

This proves the result. Hence the process can be generated on an arbitrary partition from a sequence of iid $\mathcal N(0,1)$ random variables in $O(n)$ time where $n$ is the size of the partition.

NB: This is an exact sampling technique in that it provides a sampled version of the desired process with the exactly correct finite-dimensional distributions. This is in contrast to Euler (and other) discretization schemes for more general SDEs, which incur a bias due to the approximation via discretization.

  • $\begingroup$ Just a few more remarks. (1) To get a good idea of what the continuous time process looks like, $n$ and $T$ must be chosen so that $\Delta$ is small, say less than $0.1$. (2) The inverse covariance (precision) matrix for the timeseries vector is tri-diagonal, as is its Cholesky root. $\endgroup$
    – Yves
    May 24, 2012 at 17:55
  • $\begingroup$ @Yves: Thanks for your comments. To be clear, the procedure I've outlined gives an exact realization of the continuous-time process sampled on the corresponding partition; in particular, there is no discretization error like there is in typical Euler-scheme approximation to more general SDEs. The inverse Cholesky, as shown by the construction in the answer has nonzero terms only on the diagonal and lower off-diagonal, so it's a little simpler than tridiagonal. $\endgroup$
    – cardinal
    May 24, 2012 at 18:13
  • $\begingroup$ Great answer! Does this generalize to the general OU process with arbitrary scale, $\gamma(t_i, t_j) = \exp(\alpha \; |t_i - t_j|)$? It seems like it might. $\endgroup$
    – redmoskito
    Dec 5, 2014 at 23:14

Calculate the decomposed covariance matrix by incomplete Cholesky decomposition or any other matrix decomposition technique. Decomposed matrix should be TxM, where M is only a fraction of T.


  • 2
    $\begingroup$ Can you give an explicit form of the Cholesky decomposition here? I think that the answer by cardinal achieves just that, if you think about it, by expressing $X_i$ as a function of the history. $\endgroup$
    – StasK
    May 25, 2012 at 2:54
  • 2
    $\begingroup$ The algorithm is a little too long to summarize. You can find an excellent description here: Kernel ICA, page 20. Note that this algorithm is incomplete, meaning it doesn't calculate the entire decomposition but rather an approximation (hence it is much faster). I have published code for this algorithm in the KMBOX toolbox, you can download it here: km_kernel_icd. $\endgroup$
    – Steven
    Jun 6, 2012 at 20:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.