Add random_positive_definite_matrix() and its tests.
[octave.git] / random_positive_definite_matrix.m
1 function A = random_positive_definite_matrix(integerN)
2 %
3 % Generate a random, symmetric positive-definite (SPD) matrix.
4 %
5 % Since all (SPD) matrices are diagonalizable and have positive
6 % eigenvalues, it seems likely that we can generate an SPD matrix by
7 % combining random orthogonal matrices with a diagonal matrix of
8 % random eigenvalues.
9 %
10 % I have no proof/evidence that this approach is sound.
11 %
12 % INPUT:
13 %
14 % - ``integerN`` -- The dimension of the resulting matrix.
15 %
16 % OUTPUT:
17 %
18 % - ``A`` -- A symmetric, positive definite matrix.
19 %
20 U = random_orthogonal_matrix(integerN);
21 d = unifrnd(eps, realmax, 1, integerN);
22 D = diag(d);
23 A = U*D*U';
24 end