30
$\begingroup$

Could anybody tell me how to generate random symmetric positive definite matrices using MATLAB?

$\endgroup$

closed as off-topic by user21820, Xander Henderson, Namaste, Simply Beautiful Art, user061703 Jun 22 '18 at 13:34

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question is missing context or other details: Please improve the question by providing additional context, which ideally includes your thoughts on the problem and any attempts you have made to solve it. This information helps others identify where you have difficulties and helps them write answers appropriate to your experience level." – user21820, Xander Henderson, Namaste, Simply Beautiful Art
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 1
    $\begingroup$ Random with which distribution? $\endgroup$ – Lord Soth Apr 11 '13 at 5:46
  • 1
    $\begingroup$ 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. $\endgroup$ – Daryl Apr 11 '13 at 6:08
  • 6
    $\begingroup$ The set of symmetric positive definite matrices is not compact, so such a thing as uniform distribution does not exist. $\endgroup$ – Jyrki Lahtonen Apr 11 '13 at 9:59
  • 4
    $\begingroup$ 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 :-) $\endgroup$ – Jyrki Lahtonen Apr 11 '13 at 10:10
  • 1
    $\begingroup$ 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 $\endgroup$ – srijan Apr 11 '13 at 10:25
30
$\begingroup$

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)

$\endgroup$
  • 2
    $\begingroup$ @ 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. $\endgroup$ – loup blanc Aug 22 '16 at 23:21
  • $\begingroup$ @ 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). $\endgroup$ – loup blanc Aug 22 '16 at 23:31
  • $\begingroup$ @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. $\endgroup$ – Daryl Aug 23 '16 at 8:21
  • 1
    $\begingroup$ @ 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. $\endgroup$ – loup blanc Aug 23 '16 at 10:24
  • 2
    $\begingroup$ @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 ill-conditioned matrices. If your points were relevant, then indeed this answer would be moot. $\endgroup$ – Daryl Aug 23 '16 at 11:22
16
$\begingroup$

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.

$\endgroup$
12
$\begingroup$

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 
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – user40167 Jul 22 '16 at 4:11
  • $\begingroup$ I don't understand how my answer does not produce non-positive 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. $\endgroup$ – Daryl Aug 22 '16 at 21:38
  • 1
    $\begingroup$ @ s_majee , your diagonal matrix $D$ is absolutely useless. Your result $Q^TDQ$ can be written $R^TR$, that is the first Daryl result. $\endgroup$ – loup blanc Aug 23 '16 at 0:01
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Surprised no one else came up with that solution before, it is the clearest one to me. $\endgroup$ – Sheljohn Aug 22 '16 at 17:48
  • $\begingroup$ @ 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 web.math.princeton.edu/mathlab/projects/ranmatrices/yl/… ). Thus you cannot randomly choose (in particular, uniform distribution) the eigenvalues; unfortunately the eigenvalues must also meet supplementary conditions. $\endgroup$ – loup blanc Aug 22 '16 at 23:51
  • $\begingroup$ @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. $\endgroup$ – texasflood Aug 23 '16 at 14:33
  • $\begingroup$ @ 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). $\endgroup$ – loup blanc Aug 23 '16 at 18:15
  • $\begingroup$ @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 $\endgroup$ – texasflood Aug 24 '16 at 9:21
1
$\begingroup$

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}$.

$\endgroup$
  • 1
    $\begingroup$ @ 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$. $\endgroup$ – loup blanc Aug 23 '16 at 0:24
  • $\begingroup$ @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). $\endgroup$ – Matt L. Aug 23 '16 at 14:03
  • 1
    $\begingroup$ @ 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. $\endgroup$ – loup blanc Aug 23 '16 at 18:26

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