+function A = random_positive_definite_matrix(integerN)
+ %
+ % Generate a random, symmetric positive-definite (SPD) matrix.
+ %
+ % Since all (SPD) matrices are diagonalizable and have positive
+ % eigenvalues, it seems likely that we can generate an SPD matrix by
+ % combining random orthogonal matrices with a diagonal matrix of
+ % random eigenvalues.
+ %
+ % I have no proof/evidence that this approach is sound.
+ %
+ % INPUT:
+ %
+ % - ``integerN`` -- The dimension of the resulting matrix.
+ %
+ % OUTPUT:
+ %
+ % - ``A`` -- A symmetric, positive definite matrix.
+ %
+ U = random_orthogonal_matrix(integerN);
+ d = unifrnd(eps, realmax, 1, integerN);
+ D = diag(d);
+ A = U*D*U';
+end