From 0fa5265cad135de24f0cc87195103396fdff3b91 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Fri, 22 Mar 2013 01:48:19 -0400
Subject: [PATCH 1/1] Add random_positive_definite_matrix() and its tests.

---
 random_positive_definite_matrix.m             | 24 +++++++++++++++++++
 tests/random_positive_definite_matrix_tests.m |  9 +++++++
 2 files changed, 33 insertions(+)
 create mode 100644 random_positive_definite_matrix.m
 create mode 100644 tests/random_positive_definite_matrix_tests.m

diff --git a/random_positive_definite_matrix.m b/random_positive_definite_matrix.m
new file mode 100644
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
diff --git a/tests/random_positive_definite_matrix_tests.m b/tests/random_positive_definite_matrix_tests.m
new file mode 100644
index 0000000..01afb4c
--- /dev/null
+++ b/tests/random_positive_definite_matrix_tests.m
@@ -0,0 +1,9 @@
+for n = [ 5, 10, 25, 50, 100 ]
+  A = random_positive_definite_matrix(n);
+  [R, P] = chol(A);
+
+  # chol() will set P=0 if A was symmetric positive-definite.
+  expected = 0;
+  msg = sprintf('random_positive_definite_matrix(%d) has a chol()', n);
+  unit_test_equals(msg, expected, P);
+end
-- 
2.43.2