Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Could anybody tell me How to generate random symmetric positive definite matrix using Matlab? Thank you very much for the help and suggestions.

share|cite|improve this question
Random with which distribution? – Lord Soth Apr 11 '13 at 5:46
I 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 6:08
The 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 9:59
Whoever 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
I 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
up vote 15 down vote accepted

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 = 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);


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);


In fact, if the desired eigenvalues of the random are known and stored in rc, then the command

A = sprandsym(n,density,rc);

will construct the desired matrix. (Source: MATLAB sprandsym website)

share|cite|improve this answer
Daryl Thank you very much daryl. It worked for me. Many thanks. – srijan Apr 11 '13 at 9:28
@Daryl, Thanks a lot! – nullgeppetto Nov 20 '13 at 9:39
@Daryl, nice and clear answer! – nullgeppetto Aug 14 '14 at 8:38

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;

share|cite|improve this answer
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. – Kris Joanidis Jul 22 at 4:11

A usual way in Bayesian statistics is to sample from a probability measure on real symmetric positive-definite matrices such as Wishart (or Inverse-Wishart).

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 user-fixed positive definite matrix such as the identity and df are degrees of freedom.

share|cite|improve this answer
Dear I need matlab code. I am not using any statistical software. – srijan Apr 11 '13 at 8:24
That is Matlab code! But you may or may not have Statistics toolbox. Try running 'wishrnd(eye(10),10)'. If it works, good for you. – johnny Apr 11 '13 at 8:28
Thank you very much. This code is working for me. I have one last question can we normalize the entries of the matrix to make them lie on some particular interval say $[0, 1]$ – srijan Apr 11 '13 at 9:05

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. Finally, you construct your 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}$.

share|cite|improve this answer
Thanks for your reply. Can't I just use matrix B? I am not clear why to use matrix C. thanks – srijan Apr 11 '13 at 8:24
Matrix $\bf{B}$ is only symmetric, but not necessarily positive definite. – Matt L. Apr 11 '13 at 8:32
How to find $\lambda_{min}$ in case of random matrices? Eigenvalues will keep on varying every time I run the program. – srijan Apr 11 '13 at 8:43
Of course, every time you compute a new random matrix, you also need to compute the new minimum eigenvalue. – Matt L. Apr 11 '13 at 8:44

protected by Zev Chonoles Mar 4 at 5:56

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site (the association bonus does not count).

Would you like to answer one of these unanswered questions instead?

Not the answer you're looking for? Browse other questions tagged or ask your own question.