49

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.

Scottie T
  • 11,729
  • 10
  • 45
  • 59
BHS
  • 991
  • 3
  • 12
  • 26
  • Some additional ideas can be found [here](https://math.stackexchange.com/questions/69702/algorithm-for-generating-positive-semidefinite-matrices) – HackerBoss Nov 20 '21 at 00:41

6 Answers6

85
  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
ead
  • 32,758
  • 6
  • 90
  • 153
vartec
  • 131,205
  • 36
  • 218
  • 244
  • 5
    From Wikipedia: "For any matrix A, the matrix A*A is positive semidefinite," Link: http://en.wikipedia.org/wiki/Positive-definite_matrix#Negative-definite.2C_semidefinite_and_indefinite_matrices – chillysapien Mar 06 '09 at 16:35
  • This didn't work for me unless I did `np.dot(A,A.T)` as suggested by Erogol – houbysoft Feb 27 '15 at 22:37
  • 7
    What can be said about the distribution of matrices generated this way? – Edward Doolittle May 21 '15 at 20:10
  • Wikipedia lists a few random matrix options. Some of these [Gaussian Ensemble matrices](https://en.wikipedia.org/wiki/Random_matrix) might work better, depending on your (or the OP's) problem. – hobs Sep 20 '15 at 19:56
  • 1
    You don't actually need to import `linalg` in this example. You do however need to import `numpy`. – André C. Andersen Mar 05 '17 at 15:35
  • @EdwardDoolittle If using `randn` rather than `rand`, I believe the entries in the result would be Chi-Squared distributed with the degrees of freedom depending on the matrix size – HackerBoss Nov 20 '21 at 02:23
  • If any significant number of entries in the random matrix `A` are zero, this may generate a matrix `B` which is not full-rank. A bulletproof solution would be to then add small constants along the diagonal, e.g. `B = B + np.eye(B.shape[0]) * 0.01`. Chances of this occurring with the current code are insignificant, but this comes into play more heavily if, instead of `A` being generated as a random dense matrix, it is generated as a sparse or block-sparse matrix. – Greg Kramida May 30 '23 at 13:05
24

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 = ATA (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 = QSQT (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 - 2vvT); 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.

Jason S
  • 184,598
  • 164
  • 608
  • 970
  • 1
    The method based on the Cholesky decomposition only works for _dense_ matrices. Granted, a random Householder matrix will in general be dense as well, so the same thing could be said for the SVD-based method. – ocramz Jan 20 '17 at 05:19
10

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
chillysapien
  • 2,256
  • 1
  • 26
  • 42
8

Natural distributions on positive semidefinite matrices are Wishart distributions.

plasmacel
  • 8,183
  • 7
  • 53
  • 101
Alexandre C.
  • 55,948
  • 11
  • 128
  • 197
4

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.

A. Awotunde
  • 49
  • 2
  • 1
    If A has full rank, AA' is still semidefinite positive. If A has m rows and n columns, then AA' has rank *at most* m. – Alexandre C. Sep 21 '11 at 22:45
  • Well, your statement isn't true Alex. The eignevalues of A'A are always the same as those of AA'. So if A is a square matrix of full rank, then both A'A and AA' are both square symmetric and of full rank. Simply put: If A has full rank, then AA' CANNOT be semidefinite. It must be positive-definite. – A. Awotunde Sep 23 '11 at 04:47
  • 3
    A positive definite matrix is in particular semidefinite positive. – Alexandre C. Sep 23 '11 at 08:48
  • 1
    That's true, but it's of interest to generate p.s.d. matrices that are not p.d., as well as those that are p.d., if one wants a comprehensive test set. – Evgeni Sergeev Nov 19 '13 at 06:10
  • EvgeniSergeev: your statement: "A'*A will give a positive semidefite [sic.] matrix iff and only if A is rank-deficient." and @AlexandreC's statement: "A positive definite matrix is a particular positive semidefinite matrix" cannot both be True. Your random rectangular matrix product recipe does create some positive semidefinite matrices that aren't positive definite, but 50% of the time it produces matrices that aren't even positive semidefinite, at least with [my implementation of your algorithm](https://gist.github.com/hobson/ad451869dd11439cc8ff). – hobs Sep 20 '15 at 02:20
2

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.

Allin Cottrell
  • 41
  • 3