How to create an arbitrary covariance matrix

Question

For example, in R, the MASS::mvrnorm() function is useful for generating data to demonstrate various things in statistics. It takes a mandatory Sigma argument which is a symmetric matrix specifying the covariance matrix of the variables. How would I create a symmetric $n\times n$ matrix with arbitrary entries?

Answer 1

Create an $n\times n$ matrix $A$ with arbitrary values

and then use $\Sigma = A^T A$ as your covariance matrix.

For example

n <- 4  
A <- matrix(runif(n^2)*2-1, ncol=n) 
Sigma <- t(A) %*% A

Answer 2

I like to have control over the objects I create, even when they might be arbitrary.

Consider, then, that all possible $n\times n$ covariance matrices $\Sigma$ can be expressed in the form

$$\Sigma= P^\prime\ \text{Diagonal}(\sigma_1,\sigma_2,\ldots, \sigma_n)\ P$$

where $P$ is an orthogonal matrix and $\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_n \ge 0$.

Geometrically this describes a covariance structure with a range of principal components of sizes $\sigma_i$. These components point in the directions of the rows of $P$. See the figures at Making sense of principal component analysis, eigenvectors & eigenvalues for examples with $n=3$. Setting the $\sigma_i$ will set the magnitudes of the covariances and their relative sizes, thereby determining any desired ellipsoidal shape. The rows of $P$ orient the axes of the shape as you prefer.

One algebraic and computing benefit of this approach is that when $\sigma_n \gt 0$, $\Sigma$ is readily inverted (which is a common operation on covariance matrices):

$$\Sigma^{-1} = P^\prime\ \text{Diagonal}(1/\sigma_1, 1/\sigma_2, \ldots, 1/\sigma_n)\ P.$$

Don't care about the directions, but only about the ranges of sizes of the the $\sigma_i$? That's fine: you can easily generate a random orthogonal matrix. Just wrap $n^2$ iid standard Normal values into a square matrix and then orthogonalize it. It will almost surely work (provided $n$ isn't huge). The QR decomposition will do that, as in this code

n <- 5
p <- qr.Q(qr(matrix(rnorm(n^2), n)))

This works because the $n$-variate multinormal distribution so generated is "elliptical": it is invariant under all rotations and reflections (through the origin). Thus, all orthogonal matrices are generated uniformly, as argued at How to generate uniformly distributed points on the surface of the 3-d unit sphere?.

A quick way to obtain $\Sigma$ from $P$ and the $\sigma_i$, once you have specified or created them, uses crossprod and exploits R's re-use of arrays in arithmetic operations, as in this example with $\sigma=(\sigma_1, \ldots, \sigma_5) = (5,4,3,2,1)$:

Sigma <- crossprod(p, p*(5:1))

As a check, the Singular Value decomposition should return both $\sigma$ and $P^\prime$. You may inspect it with the command

svd(Sigma)

The inverse of Sigma of course is obtained merely by changing the multiplication by $\sigma$ into a division:

Tau <- crossprod(p, p/(5:1))

You may verify this by viewing zapsmall(Sigma %*% Tau), which should be the $n\times n$ identity matrix. A generalized inverse (essential for regression calculations) is obtained by replacing any $\sigma_i \ne 0$ by $1/\sigma_i$, exactly as above, but keeping any zeros among the $\sigma_i$ as they were.

Answer 3

You can simulate random positive definite matrices from the Wishart distribution using the function "rWishart" from stats (included in base)

n <- 4
rWishart(1,n,diag(n))

From the documentation for rWishart:

Usage should be:

rWishart(n, df, Sigma)

where,

• n: the number of samples.
• df: the degrees of freedom, i.e. the number of dimensions of the matrix.
• Sigma: a positive definite scaling matrix.

Answer 4

There is a package specifically for that, clusterGeneration (written among other by Harry Joe, a big name in that field).

There are two main functions:

• genPositiveDefMat generate a covariance matrix, 4 different methods
• rcorrmatrix : generate a correlation matrix

Quick example:

library(clusterGeneration)
#> Loading required package: MASS
genPositiveDefMat("unifcorrmat",dim=3)
#> $egvalues
#> [1] 15.408962  5.673916  1.228842
#> 
#> $Sigma
#>          [,1]     [,2]     [,3]
#> [1,] 6.714871 1.643449 6.530493
#> [2,] 1.643449 6.568033 2.312455
#> [3,] 6.530493 2.312455 9.028815
genPositiveDefMat("eigen",dim=3)
#> $egvalues
#> [1] 8.409136 4.076442 2.256715
#> 
#> $Sigma
#>            [,1]       [,2]      [,3]
#> [1,]  2.3217300 -0.1467812 0.5220522
#> [2,] -0.1467812  4.1126757 0.5049819
#> [3,]  0.5220522  0.5049819 8.3078880

^{Created on 2019-10-27 by the reprex package (v0.3.0)}

Finally, note that an alternative approach is to do a first try from scratch, then use Matrix::nearPD() to make your matrix positive-definite.

Answer 5

In my case I still want the eigenvalues to be randomly drawn so, to control the condition number of the covariance matrix does not exceed, say, $\kappa > 1$, one can generate $\sigma_k = 1 + (\kappa - 1) \cdot \mathcal{U}_k$ where $\mathcal{U}_k$ is uniformly sampled in $[0,1]$ for all $k=1,\dots,n$.

Following @whuber (very comprehensive) algorithm, here is a code snippet in Python with all kinds of assertions and data reshaping.

import numpy as np
d = 5
sigmas = np.arange(1, 6).reshape(-1, 1)  # [1, 2, 3, 4, 5]
# or sigmas = (1 + (kappa - 1) * np.random.uniform(size = d)).reshape(-1, 1)
Q, R = np.linalg.qr(np.random.normal(size = d**2).reshape(d, d))

S = Q.T @ (sigmas * Q)
# Check that S is psd
assert np.all(np.linalg.eigvals(S) > 0)
# Check that S is symmetric
assert np.allclose(S.T, S)
S
_, sigmas_, Q_ = np.linalg.svd(S)
# Check SVD retrieves original parameters
assert np.allclose(np.abs(Q), np.abs(Q_[::-1]))  # up to sign
assert np.allclose(sigmas.reshape(-1), sigmas_[::-1])