A simple algorithm for generating positive-semidefinite matrices

Question

I want to generate positive random semi-definite matrices. I am looking for an algorithm or more preferably an simple implementation of the algorithm in C, matlab, java or any language.

Answer 1

1. generate random matrix
2. multiply it by it's own transposition
3. you have obtained a positive semi-definite matrix.

Example code (Python):

from scipy import random, linalg
matrixSize = 10 
A = random.rand(matrixSize,matrixSize)
B = numpy.dot(A,A.transpose())
print 'random positive semi-define matrix for today is', B

Answer 2

You need to be clear on your definition of "random". What are your constraints on the resulting matrix? Do you want the coefficients to be uniformly or normally distributed? Do you want the eigenvalues to have a particular distribution? (etc.)

There are a number of ways to generate positive semidefinite matrices M, including:

1. Given an arbitrary matrix A, compute M = A^TA (constructing a Cholesky decomposition)
2. Given an arbitrary diagonal matrix S with nonnegative diagonal entries, and an orthonormal matrix Q of the same size, compute M = QSQ^T (constructing a singular value decomposition)

For numerical reasons I'd probably choose the second approach by generating the diagonal matrix with desired properties, then generating Q as the composition of a number of Householder reflections (generate a random vector v, scale to unit length, H = I - 2vv^T); I suspect you'd want to use K * N where N is the size of the matrix M, and K is a number between 1.5-3 (I'm guessing on this) that ensures that it has enough degrees of freedom.

You could also generate an orthonormal matrix Q using Givens rotations: pick 2 distinct values from 1 to N and generate a Givens rotation about that pair of axes, with an angle uniformly distributed from 0 to 2 * pi. Then take K * N of these (same reasoning as above paragraph) and their composition yields Q.

edit: I'd guess (not sure) that if you have coefficients that are independently-generated and normally distributed, then the matrix as a whole would be "normally distributed" (whatever that means). It's true for vectors, at least. (N independently-generated Gaussian random variables, one for each component, gives you a Gaussian random vector) This isn't true for uniformly-distributed components.

Answer 3

If you can generate a random matrix in your chosen language, then by using the property that a matrix multiplied by its transpose is positive semi-definte, you can generate a random positive semi-definite matix

In Matlab it would be as simple as

% Generate a random 3x3 matrix
    A = rand(3,3) 
% Multiply by its tranpose
    PosSemDef = A'*A 

Answer 4

Natural distributions on positive semidefinite matrices are Wishart distributions.

Answer 5

A'*A will give a positive semidefite matrix iff and only if A is rank-deficient. So the answers stated above and that copied from wikipedia are not generally true. To compute a positive semidefinite matrix simply take any rectangular m by n matrix (m < n) and multiply it by its transpose. I.e. if B is an m by n matrix, with m < n, then B'*B is a semidefinite matrix. I hope this helps.

Answer 6

To clarify a little (I hope). Let A be a random matrix (for example, populated by random normal variates), m x n with m >= n. Then if A is of full column rank, A'A will be positive definite. If A is of rank < n then A'A will be positive semidefinite (but not positive definite).

A random normal matrix with m >= n will almost surely be of full rank; to generate a rank-deficient matrix one can append one or more columns that are linear combinations of other columns.