From f7394188abd4c5ca3a8e6a74ac16b97a4966df5e Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Fri, 22 Mar 2013 00:57:45 -0400
Subject: [PATCH] Add the random_orthogonal_matrix() function and its tests.

---
 random_orthogonal_matrix.m             | 25 +++++++++++++++++++++++++
 tests/random_orthogonal_matrix_tests.m |  8 ++++++++
 2 files changed, 33 insertions(+)
 create mode 100644 random_orthogonal_matrix.m
 create mode 100644 tests/random_orthogonal_matrix_tests.m

diff --git a/random_orthogonal_matrix.m b/random_orthogonal_matrix.m
new file mode 100644
index 0000000..8c7cea5
--- /dev/null
+++ b/random_orthogonal_matrix.m
@@ -0,0 +1,25 @@
+function U = random_orthogonal_matrix(integerN)
+  %
+  % Generate a random orthogonal matrix.
+  %
+  % INPUT:
+  %
+  %   - ``integerN`` -- The dimension of the resulting square matrix.
+  %
+  % OUTPUT:
+  %
+  %   - ``U`` -- An orthogonal matrix of dimension integerN.
+  %
+  % REFERENCES:
+  %
+  %   1. G.W. Stewart, Efficient Generation of Random Orthogonal Matrices,
+  %      SIAM J. Numer. Analysis, 1980, pp. 403--409, Section 3.
+  %
+
+  % We begin by computing a random matrix A.
+  A = rand(integerN);
+
+  % The Q-R decomposition of A will give us an orthogonal factor of A.
+  % See the reference for why this doesn't suck.
+  [U, R] = qr(A);
+end
diff --git a/tests/random_orthogonal_matrix_tests.m b/tests/random_orthogonal_matrix_tests.m
new file mode 100644
index 0000000..f17e639
--- /dev/null
+++ b/tests/random_orthogonal_matrix_tests.m
@@ -0,0 +1,8 @@
+for n = [ 1, 10, 25, 50, 100 ]
+  U = random_orthogonal_matrix(n);
+  actual = norm(U);
+  expected = 1;
+
+  msg = sprintf('random_orthogonal_matrix(%d) has norm 1', n);
+  unit_test_equals(msg, expected, actual);
+end
-- 
2.43.2