42
$\begingroup$

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?

$\endgroup$
1

5 Answers 5

41
$\begingroup$

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
$\endgroup$
7
  • $\begingroup$ Likewise, Sigma <- A + t(A). $\endgroup$
    – rsl
    May 31, 2016 at 8:45
  • 16
    $\begingroup$ @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 $\endgroup$
    – Henry
    May 31, 2016 at 10:30
  • $\begingroup$ Yes, I noticed that R returns error in the event my suggested way produced unsuitable matrix. $\endgroup$
    – rsl
    May 31, 2016 at 18:32
  • 8
    $\begingroup$ Note that if you prefer a correlation matrix for better interpretability, there is the ?cov2cor function, which can be applied subsequently. $\endgroup$ May 31, 2016 at 22:58
  • 1
    $\begingroup$ @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$ $\endgroup$
    – Henry
    Jan 11, 2018 at 9:26
50
$\begingroup$

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.

$\endgroup$
5
  • 1
    $\begingroup$ It might help to demonstrate how to use the rows of $P$ to orient the axes as preferred. $\endgroup$ May 31, 2016 at 23:33
  • 1
    $\begingroup$ Might be worth mentioning that the singular values in svd(Sigma) will be reordered -- that confused me for a minute. $\endgroup$
    – FrankD
    Mar 28, 2018 at 10:45
  • 1
    $\begingroup$ @whuber could you edit please the line of code to parametrize n: crossprod(p, p*(n:1)). In order to improve reproducibility, good answer! $\endgroup$ Feb 2, 2020 at 15:54
  • $\begingroup$ @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. $\endgroup$
    – whuber
    Feb 2, 2020 at 16:37
  • $\begingroup$ @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. $\endgroup$
    – whuber
    Mar 13, 2023 at 15:31
11
$\begingroup$

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.
$\endgroup$
4
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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])
$\endgroup$

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.