Could anybody tell me how to generate random symmetric positive definite matrices using MATLAB?
 327
 4
 23
 11,854
 9
 64
 108

1Random with which distribution? – Lord Soth Apr 11 '13 at 05:46

@LordSoth Uniform distribution – Srijan Apr 11 '13 at 05:47

1I can give an algorithm that will generate a "random" symmetric, positive definite matrix, but the entries are by no means uniformly distributed, if they follow a standard distribution at all. – Daryl Apr 11 '13 at 06:08

@Daryl That may also work for me. I would be thankful to you dear. – Srijan Apr 11 '13 at 06:12

8The set of symmetric positive definite matrices is not compact, so such a thing as uniform distribution does not exist. – Jyrki Lahtonen Apr 11 '13 at 09:59

@JyrkiLahtonen Thanks for info. I am unaware of this fact. I have to use random symmetric positive definite matrix for my work. – Srijan Apr 11 '13 at 10:06

6Whoever tells you to do that should then also specify the distribution. If not directly, then via a description of the random process that you are expected to study. I suspect that Wishart would be good one (see Johnny's answer). But really your task has not been fully specified, so your responsibility might be to go to your boss, and ask for more information  informing him/her about the danger of "garbage in/ garbage out" simulations to be done otherwise :) – Jyrki Lahtonen Apr 11 '13 at 10:10

1I agree with you . My task is to compute weighted moore penrose inverse $A^{+}_{M,N}$ for randomly generated matrices, where $M$ and $N$ are given symmetric positive definite matrix. I have just figured out that for a gien matrix $A$, $AA'$ is a symmetric positive definite matrix. Perhaps this may work – Srijan Apr 11 '13 at 10:25
5 Answers
The algorithm I described in the comments is elaborated below. I will use $\tt{MATLAB}$ notation.
function A = generateSPDmatrix(n)
% Generate a dense n x n symmetric, positive definite matrix
A = rand(n,n); % generate a random n x n matrix
% construct a symmetric matrix using either
A = 0.5*(A+A'); OR
A = A*A';
% The first is significantly faster: O(n^2) compared to O(n^3)
% since A(i,j) < 1 by construction and a symmetric diagonally dominant matrix
% is symmetric positive definite, which can be ensured by adding nI
A = A + n*eye(n);
end
Several changes are able to be used in the case of a sparse matrix.
function A = generatesparseSPDmatrix(n,density)
% Generate a sparse n x n symmetric, positive definite matrix with
% approximately density*n*n non zeros
A = sprandsym(n,density); % generate a random n x n matrix
% since A(i,j) < 1 by construction and a symmetric diagonally dominant matrix
% is symmetric positive definite, which can be ensured by adding nI
A = A + n*speye(n);
end
In fact, if the desired eigenvalues of the random matrix are known and stored in the vector rc
, then the command
A = sprandsym(n,density,rc);
will construct the desired matrix. (Source: MATLAB sprandsym website)
 5,343
 3
 24
 39

4@ Daryl , your solution $AA^T$ works. Yet your solution $(A+A^T)/2+nI_n$ is diagonally dominant; thus it is not random amongst the symmetric $>0$ matrices. – Aug 22 '16 at 23:21

@ Daryl , if "rand(n,n)" randomly give $A_{i,j}\in(1,1)$, then you obtain a very special matrix; it is better to choose the $(a_{i,j})$ i.i.d. and following a normal law (for example). – Aug 22 '16 at 23:31

1@loupblanc The OP did not define what they meant by random matrix, so there is no "correct" distribution for the matrix entries and both are correct solutions to the problem. The matrix $(A+A^T)/2 + nI$ is a random matrix, but the entries will probably not follow any standard distribution. Also, rand gives entries in $[0,1)$, and in conjunction with my first statement, I don't see the relevance of your second comment. – Daryl Aug 23 '16 at 08:21

1@ Daryl , a random matrix $A=[a_{i,j}]$ is (in general) so that the $(a_{i,j})$ are i.i.d. and follow a standard probability law generally $N(0,1)$ or an uniform law in $[a,a]$ The properties of the eigenvalues of a symm. $A>0$ are similar in both cases (see Tao's papers). In particular, your random entries have $1/2$ as mean, that is original. Your matrix $(A+A^T)/2+nI$ is diagonally dominant and moreover has only positive entries; if you don't see why such a matrix is not random amongst the symm $>0$ ones, then I can do nothing for you. – Aug 23 '16 at 10:24

@ Daryl , Even more serious: the distribution and estimates of the eigenvalues of your matrix (when $n$ is large) do not satisfy the standard results obtained by researchers as Maestro Tao (cf. my comments below to the Matt L.'s and texasflood's posts). Finally, if you do numerical experiments using this type of matrix, then the results are biased and are not admissible. – Aug 23 '16 at 10:25

@ Daryl , the last but not the least: if your $(A_{i,j})$ are randomly chosen in $[0,1)$, then $AA^T$ has only positive entries and consequently, is not random amongst the symm $>0$ matrices. – Aug 23 '16 at 10:30

3@loupblanc They are all good points, but are *irrelevant* to this question as the OP did not define that the entries were to have a particular distribution. Furthermore, if no distribution is specified, one cannot expect any "standard results" to hold.This approach is standard in numerical analysis for generating SPD matrices to test algorithms; in which case the only bias is due to well conditioned matrices and a more nuanced matrix generation strategy is required to apply the algorithm with illconditioned matrices. If your points were relevant, then indeed this answer would be moot. – Daryl Aug 23 '16 at 11:22

@ Daryl , I understand your point of view. Yet, for example, a diagonally dominant matrix is (in general ) very stable; then there exist algorithms that work well, when using such a matrix, and badly, when using a not diagonally dominant matrix. On the other hand, the OP is far from being a specialist of our problem of course, he's making progress in almost three years . It is our responsibility to give him the right assumptions. – Aug 23 '16 at 23:34

@Daryl. Your first approach does not always produce PSD matrix in MATLAB (verified in MATLAB). Link for my [code](https://mega.nz/#!pUsglDKZ!0HZ41YDPJpbVTT7I_ZjjOoVVku8mK9NIh8RsU0MZThM). Run the file in MATLAb a couple of times and sigma matrix is sometimes not PSD as raised by chol factorization error. – CKM Oct 21 '16 at 09:37

@chandresh I cannot get the link to work to show me the code; it keeps redirecting me to ask if I want to see the code. It should be clear the resulting matrix is symmetric. The first approach will produce a diagonally dominant matrix, which is guaranteed to be positive definite. If your code does not produce a diagonally dominant matrix, it can't be a replication of my above code (unless MATLAB has changed its rand implementation). If it is diagonally dominant, the issue lies within the chol method, which may complain due to rounding error. Please post an A matrix for which this is not SPD. – Daryl Oct 21 '16 at 10:56

@Daryl. I checked and the link is working fine for me. Try to download the code through your browser. I'm sorry its `positive definite` error not PSD. A =[ 2.5767 4.2095; 4.2095 5.5344]. Try this [link](https://drive.google.com/file/d/0BzaBVVWfPCqIM0RkbTBjd1dUVWM/view?usp=sharing) – CKM Oct 21 '16 at 11:52
A usual way in Bayesian statistics is to sample from a probability measure on real symmetric positivedefinite matrices such as Wishart (or InverseWishart).
I don't use Matlab but a quick check on Google gives this command (available in the Statistics toolbox):
W = wishrnd(Sigma,df)
where Sigma
is some userfixed positive definite matrix such as the identity and df
are degrees of freedom.
 2,163
 1
 18
 25
While Daryl's answer is great, it gives symmetric positive definite matrices with very high probability , but that probability is not 1. This method gives a random matrix being symmetric positive definite matrix with probability 1.
My answer relies on the fact that a positive definite matrix has positive eigenvalues.
The matrix symmetric positive definite matrix A can be written as , A = Q'DQ , where Q is a random matrix and D is a diagonal matrix with positive diagonal elements. The elements of Q and D can be randomly chosen to make a random A.
The matlab code below does exactly that
function A = random_cov(n)
Q = randn(n,n);
eigen_mean = 2;
% can be made anything, even zero
% used to shift the mode of the distribution
A = Q' * diag(abs(eigen_mean+randn(n,1))) * Q;
return

I don't see why Q' * Q isn't symmetric positive definite with probability 1. Q is invertible with probability 1, which is sufficient t ensure it. – Jul 22 '16 at 04:11

I don't understand how my answer does not produce nonpositive definite matrices. The matrix returned from either of my functions is absolutely diagonally dominant, which is a sufficient condition for a matrix to be positive definite. – Daryl Aug 22 '16 at 21:38

1@ s_majee , your diagonal matrix $D$ is absolutely useless. Your result $Q^TDQ$ can be written $R^TR$, that is the first Daryl result. – Aug 23 '16 at 00:01
For the case where you want a complex matrix (which not all previous answers address), you can do
Q = orth(randn(n) + randn(n)*i);
D = diag(abs(randn(n, 1)) + 0.3);
A = Q*D*Q';
If you want a real matrix, do
Q = orth(randn(n));
D = diag(abs(randn(n, 1)) + 0.3);
A = Q*D*Q';
If you want a semi positive definite matrix, remove the 0.3. One may also change the 0.3 to any other appropriate positive number depending on how positive definite they want the matrix to be guaranteed to be.
This also works in Octave.
 545
 4
 13

Surprised no one else came up with that solution before, it is the clearest one to me. – Jonathan H Aug 22 '16 at 17:48

1@ texasflood , it's a bad idea (I assume that $Q$ is orthogonal). If $A$ is symmetric $>0$ (where the $A_{i,j}$ are i.i.d. Gaussian distributed, then the distribution of the spacing between adjacent eigenvalues is very special (see http://web.math.princeton.edu/mathlab/projects/ranmatrices/yl/randmtx.PDF ). Thus you cannot randomly choose (in particular, uniform distribution) the eigenvalues; unfortunately the eigenvalues must also meet supplementary conditions. – Aug 22 '16 at 23:51

1@loupblanc Hmmm I see part of your point, I'll have to think about that. What I don't get is how do you ensure that $A$ is Hermitian and $> 0$ if $A_{i,j}$ are i.i.d. and Gaussian? That statement makes no sense. Also you say that the eigenvalues are being chosen from a uniform distribution. They are not, they are normally distributed. – texasflood Aug 23 '16 at 14:33

1@ texasflood , of course, the symmetric $>0$ matrix is $AA^T$ (perhaps you did not realise...). More seriously, I did not see the "n" in "randn" (line 2 in your procedure); here, each eigenvalue of your $A$ is in $[0.3,1.3]$ and their distributions are absolutely not the distributions associated to a random symm. $>0$ matrix (see the Tao's paper or my comment of Matt L.'s post). – Aug 23 '16 at 18:15

@loupblanc Ah OK I see. I guess different people will want different distributions depending on their use case, and some might not even care. I wonder if there is a way to distribute the eigenvalues so that the matrices generated in my method have the same distribution as a p.d. $AA^T$ with normally distributed elements – texasflood Aug 24 '16 at 09:21

@texasflood Thank you very much for your answer. I am sorry too for my late response. Nice way to generate complex matrix. – Srijan Oct 27 '16 at 10:39
The following is not computationally efficient but very simple. You could fill a matrix $\bf A$ with random values, computed for some desired distribution. Then you define a new matrix $\bf B = \bf{A} + \bf{A}^T$ in order to get a symmetric matrix. Then you use matlab to compute the eigenvalues of this matrix. If $\mathbf{B}$ doesn't happen to be positive definite, construct a new matrix matrix by
$$\bf{C} = \bf{B} + (\lambda_{min} + \delta)\bf{I}$$
where $\lambda_{min}$ is the absolute value of the smallest eigenvalue of $\bf{B}$ and $\delta$ is some small positive constant which defines the smallest eigenvalue of the your final matrix $\bf{C}$.
 10,199
 15
 25

1@ Matt L. , that does not work. Indeed let $A$ be symmetric $>0$ (where the $A_{i,j}$ are i.i.d. Gaussian distributed) and let $\lambda_1\geq\cdots\geq \lambda_n$ be the spectrum of $A$. Then, when $n$ is large, $\lambda_1^2\approx 4n$ and $E(\lambda_n)^2=O(1/n)$. Then it seems (to me) very difficult to choose a correct $\delta$. – Aug 23 '16 at 00:24

@loupblanc: There is no "correct" $\delta$. Just choose any positive value for it, and that will be the minimum eigenvalue of the new matrix $\mathbf{C}$ (assuming that $\mathbf{B}$ was not already positive definite). – Matt L. Aug 23 '16 at 14:03

1@ Matt L. , just not. About the choice of $\lambda_n$ (for example), the selection of $\delta$ is not at all arbitrary. More important: the eigenvalues of a random symm. $>0$ matrix are not the translated values of the eigenvalues of a random symmetric matrix $A+A^T$. The distributions are not of the same type. – Aug 23 '16 at 18:26

Along similar lines, you can take $L$, the lower triangular matrix of A, and $LL^T$ will be positive definite. – Evan Zamir Jan 11 '21 at 21:04