🚀 go-pugleaf

>I'm looking for methods of creating correlated random variables.
>I am especially interested in "fast" methods.  I already know that
>given a covariance matrix C_{ij}, and a set of independent random variables,
>u_i, I can create a set of random variables y = Au, such that
>the covariance matrix of the y_i is C_{ij}.  The trick is to find
>"Cholesky factorization" of C, i.e a lower-triangular matrix A, such that
> A A^{T} = C.  This is an N^3 problem, of great interest to linear
>algebraists, as it is closely related to solving symmetric positive
>definite systems of equations. 

>The problem is that Cholesky factorization requires O(N^3) time, and
>I have N of order 10^7.  However, my covariance matrix is not
>completely arbitrary.  What I am really trying to do is create
>representative set of samples of f(t) where f is a stochastic process
>with a smooth, well-defined correlation function (call it g).  The covariance
>matrix has the form:

>C_{ij} = g(t_i - t_j)

Hi, John,

My E-mail has been bounced, so I am posting this.

You can generate a set of correlated random numbers based upon a known
auto-correlation function through Fourier transform.  Since you know the
correlation function, this method may be well suited for your purpose.

The method is based upon the fact that the auto-correlation of 
a (stationary series) and the power spectrum are a Fourier transform
pair. So if you take your correlation function and apply Fourier tansform,
you get the power spectrum and, hence, the amplitude spectrum. Then you
can add a random phase term (e.g., random phase distributed uniformly
between 0 and 2*pi). Inverse transforming the made-up Fourier spectrum
yields the correlated random numbers you want. Since you use the same
correlation function, for every set of random phases you put in you
get a different set of correlated random number. Yet, they all have
the same covariance. The Fourier transform can be done by FFT, so the
procedure is really fast. (in fact, I used the method you proposed for some of
my work before,  but its slow speed prompted me to look into the Fourier 
transform approach.) 

This method can also be easily extended to generate correlated random 
numbers in multidimensions.One drawback of the approach is that it's only 
good for generating stationary random numbers.

------------------------------------------------------------
  Dr. Yaoguo Li    Department of Geophysics and Astronomy
  li@geop.ubc.ca   University of British Columbia, Canada
------------------------------------------------------------