--- /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