A simple algorithm for generating positive-semidefinite matrices

Natural distributions on positive semidefinite matrices are Wishart distributions.

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 

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

Example code (Python):

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

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.