Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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

share|improve this question
    
Random with which distribution? –  Lord Soth Apr 11 '13 at 5:46
1  
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
2  
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
3  
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
1  
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

4 Answers 4

up vote 8 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);

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 are known and stored in rc, then the command

A = sprandsym(n,density,rc);

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

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

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|improve this answer
    
Dear I need matlab code. I am not using any statistical software. –  srijan Apr 11 '13 at 8:24
3  
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 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|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

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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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