From: Michael Orlitzky Date: Fri, 22 Mar 2013 04:57:45 +0000 (-0400) Subject: Add the random_orthogonal_matrix() function and its tests. X-Git-Url: http://gitweb.michael.orlitzky.com/?p=octave.git;a=commitdiff_plain;h=f7394188abd4c5ca3a8e6a74ac16b97a4966df5e Add the random_orthogonal_matrix() function and its tests. --- 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