# How to create an arbitrary covariance matrix

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?

• Jun 2, 2016 at 13:50

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

• Likewise, Sigma <- A + t(A).
– rsl
May 31, 2016 at 8:45
• @MoazzemHossen: Your suggestion will produce a symmetric matrix, but it may not always be positive semidefinite (e.g. your suggestion could produce a matrix with negative eigenvalues) and so it may not be suitable as a covariance matrix May 31, 2016 at 10:30
• Yes, I noticed that R returns error in the event my suggested way produced unsuitable matrix.
– rsl
May 31, 2016 at 18:32
• Note that if you prefer a correlation matrix for better interpretability, there is the ?cov2cor function, which can be applied subsequently. May 31, 2016 at 22:58
• @B11b: You need your covariance matrix to be positive semi-definite. That would put some limits on the covariance values, not totally obvious ones when $n \gt 2$ Jan 11, 2018 at 9:26

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.

• It might help to demonstrate how to use the rows of $P$ to orient the axes as preferred. May 31, 2016 at 23:33
• Might be worth mentioning that the singular values in svd(Sigma) will be reordered -- that confused me for a minute. Mar 28, 2018 at 10:45
• @whuber could you edit please the line of code to parametrize n: crossprod(p, p*(n:1)). In order to improve reproducibility, good answer! Feb 2, 2020 at 15:54
• @CristóbalAlcázar Thank you for your comment. Because that's a specific example and is not intended for general purposes, I think it's best to keep the values hard-coded as they are.
– whuber
Feb 2, 2020 at 16:37
• @FrankD Many implementations of SVD do not specify that the singular values they return will be ordered, but in every case I have seen -- including the svd function in R -- they are ordered by decreasing value. Since in this post I initially specified an ordered sequence of singular values, the results returned by svd ought to be identical to the input, even up to order.
– whuber
Mar 13 at 15:31

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.

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)
genPositiveDefMat("unifcorrmat",dim=3)
#> $$egvalues #>  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 #>  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.

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])