--- /dev/null
+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
--- /dev/null
+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