Generating Correlated Random Vectors

Developing Models for Kalman Filters

For purposes of studying our noise models, we want to be able to exercise our state transition equations with noise sources that have known covariance properties. In this installment, we will discuss a way to simulate the correlated noise sources.

The Choleski factorization

Variance matrices have a symmetry property, Q = QT. This is because each covariance matrix term Qi,j results from averaging pairs of terms xi · xj, which are always equal to the commuted product terms xj · xi. The main diagonal terms of a covariance matrix result from products of one random term with itself, yielding a sum of squared values, so it is never negative. Excluding pathological cases, all of the main diagonal terms of the covariance matrix will be positive and will dominate the off-diagonal terms, which could be of either sign.

What is a little less obvious is that covariance matrices (and other matrices like them) can be factored into "square root" matrix forms. That is, for any such (symmetric, real-valued, positive definite) matrix Q, it is possible to find a matrix factor q such that

  Q    =    q  qT   

In fact, an infinite number of such factorizations are possible. To show this, suppose that you already know of one factorization Q = q qT . We can note that any orthonormal matrix U will have the property that its transpose is its own inverse. Consequently:

 ( q  U )  ( q  U ) T   =   ( q  U )  ( UT qT )    =    q  qT 

So if q is a factor matrix, ( q U ) will be a similar sort of "square root" matrix also.

But there is one particular factorization that is computationally very useful: the one for which the q matrix is lower triangular. This one is easy to determine, and easy to use.


Finding triangular factors

Let's step through a factorization of a symmetric 3x3 covariance matrix X. Once you see the pattern, you can extend this to any order matrix. If you are a purist, you are welcome to construct a formal proof by integer induction.

We start with the factored form that we want to find. We know the terms of the X covariance matrix. When we finish, we want to know the terms in factor matrix x, with its corresponding transpose. As shown in the following illustration, the grayed-out terms are all equal to zero and are of no consequence.

We first observe that the X1,1 term results from the dot product of the first row of factor matrix x with the first column of its transpose, which turns out to be the same vector terms but with the terms in a different orientation. Since only the first term of the vectors is non-zero, this dot product is relatively trivial. We can determine immediately that the factor matrix must satisfy the condition

X 1,1   =   x1,1 * x1,1     →     x1,1 = √ X1,1

We can color highlight in red the terms that have been processed and that have known values.

If we now consider work downwards and consider the next two terms in the first column of the X matrix, these terms must be formed by dot products between the other two rows and of the factor matrix and the known first column of the transposed matrix (which happens to have only one non-zero term). These dot products also involve only one multiplication each, and the following conditions must hold.

X1,2   =   x1,2 * x1,1

X1,3   =   x1,3 * x1,1

Knowing the value of x1,1 it quickly follows that

x1,2     =     X1,2 / x1,1

x1,3     =     X1,3 / x1,1

So far, we have the following processed terms.

But when we take advantage of the fact that the second factor matrix will be the transpose of the first one, some additional terms are known and need not be calculated again.

Now we can start working down the second column of the X covariance matrix. The X2,2 term results from the product of the second row of the factor matrix with its transpose, so the following conditions must be satisfied.

X2,2   =   x1,2 * x1,2 * x2,2 * x2,2

We already know the value of x1,2 , so we quickly determine that

x2,2    =     √ ( X2,2 - x1,2 * x1,2)

Working the way down the rest of the second column, we must satisfy

X2,3   =   x1,2 * x1,3 + x2,2 * x2,3

Knowing the values of x1,2 and x1,3, we can calculate

x2,3   =   ( X2,3 - x1,2 * x1,3) / x2,2)

Again taking advantage of the symmetry, and reviewing how many terms remain, we find that only one unknown term remains.

For this term, we can determine that it must satisfy the condition

X3,3 = x1,3 * x1,3 + x2,3 * x2,3 + x3,3 * x3,3

x3,3    =     √ ( X3,3 - x1,3 * x1,3 - x2,3 * x2,3)

The factorization is now done. To verify, calculate the product of the two triangular factors and check that the result reproduces the original covariance matrix.

Once you are completely comfortable with the Choleski decomposition, you can forget about all of the internal details, and simply let the Octave / Matlab chol command to do all of these matrix calculations for you.


A quick example

Define a real-valued, symmetric initial matrix.

Qmat =
   1.0000000   0.8100000   1.0000000   0.0100000
   0.8100000   0.7186000   0.7475000   0.0706000
   1.0000000   0.7475000   1.0664062  -0.0018750
   0.0100000   0.0706000  -0.0018750   1.7187000

Apply the chol command to produce the upper and lower symmetric triangular factors.

f2 = chol(Qmat)
f1 = f2'

f1 =
   1.00000   0.00000   0.00000   0.00000
   0.81000   0.25000   0.00000   0.00000
   1.00000  -0.25000   0.06250   0.00000
   0.01000   0.25000   0.81000   1.00000

f2 =
   1.00000   0.81000   1.00000   0.01000
   0.00000   0.25000  -0.25000   0.25000
   0.00000   0.00000   0.06250   0.81000
   0.00000   0.00000   0.00000   1.00000

Multiply the two factors together again and compare to the original symmetric matrix to see if they match within statistical tolerances.

restored = f1 * f2

restored =
   1.0000000   0.8100000   1.0000000   0.0100000
   0.8100000   0.7186000   0.7475000   0.0706000
   1.0000000   0.7475000   1.0664062  -0.0018750
   0.0100000   0.0706000  -0.0018750   1.7187000


Yes, but what good is this?

For right now, here is the particular application we care about. Suppose that you start with three statistically independent streams of random values with zero mean and variance 1.0, collected to form a column vector x. As you know from past articles, when a huge number of the rank-one matrix products x xT are averaged, the results will converge to a covariance matrix. Since each random term is independent from the others, this matrix converges to an identity matrix.

So now, suppose that the desired covariance matrix is Q. Compute its Choleski factorization q qT.

      Q   =   q qT

Now take the same set of random vectors x as before, and transform them by calculating the matrix-vector product q x. Observe that the rank-one product terms using these vectors will look like the following.

      (q x) (q x)T  =  q x xT qT  =  (q) (x xT) (qT)

Now average a huge number of these products to obtain the expected value. The outer q terms are just constant multipliers, so this converges to

      q E( x xT) qT  =

      q I qT  =

      q qT  = 


And there you go! You can now take a set of independent, uncorrelated random variables and combine them in such a way that you get a new set of random variables showing the noise correlation properties you want.


A demonstration

Let's build some random vectors and verify that this works. Generate vectors of independent, normally-distributed random vectors with variance 1.0. Apply the Choleski matrix to each vector, and then calculate the estimated covariance matrix using rank-one updates. Then let's see if the covariance estimate is the one we wanted.

% Script to verify correlation in a generated random sequence.

% Define desired correlation matrix and its factors
Qmat = [
   0.1000000   0.0800000   0.0000000   0.0100000;
   0.0800000   0.7186000   0.0275000   0.0000000;
   0.0000000   0.0275000   0.0964062  -0.0018750;
   0.0100000   0.0000000  -0.0018750   0.2564062 ]
 qfact2 = chol(Qmat);
 qfact1 = qfact2'
 estm = zeros(4,4);
 for term=1:250000
   vect= qfact1 * randn(4,1);
   estm = estm + vect*vect';
 estm = (1/25000)*estm
estm =
   1.0005e+00   8.1153e-01   2.5921e-04   1.0314e-01
   8.1153e-01   7.1906e+00   2.6385e-01   5.8694e-04
   2.5921e-04   2.6385e-01   9.6338e-01  -1.7193e-02
   1.0314e-01   5.8694e-04  -1.7193e-02   2.5680e+00

Next time, we will apply this to the study of covariance propagation in state transition equations.



Contents of the "Developing Models for Kalman Filters" section of the website, including this page, are licensed under a Creative Commons Attribution-ShareAlike 4.0 International License unless otherwise specifically noted. For complete information about the terms of this license, see The license allows usage and derivative works for individual or commercial purposes under certain restrictions. For permissions beyond the scope of this license, see the contact information page for this site.