Add random_positive_definite_matrix() and its tests.

[octave.git] / random_positive_definite_matrix.m

diff --git a/random_positive_definite_matrix.m b/random_positive_definite_matrix.m
new file mode 100644 (file)
index 0000000..fb3be2c
--- /dev/null
+++ b/random_positive_definite_matrix.m
@@ -0,0 +1,24 @@
+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