How to generate random symmetric positive definite matrices using MATLAB? [closed]
Could anybody tell me how to generate random symmetric positive definite matrices using MATLAB?
Solution 1:
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)
Solution 2:
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.
Solution 3:
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
Solution 4:
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.