randrot

 Generates uniformly random rotation matrices.

 function Q = randrot(n, N)

 Q is an n-by-n-by-N array such that each slice Q(:, :, i) is a random
 orthogonal matrix of size n of determinant +1 (i.e., a matrix in SO(n)),
 sampled from the Haar measure (uniform distribution).

 By default, N = 1.

 Complexity: N times O(n^3).
 Theory in Diaconis and Shahshahani 1987 for the uniformity on O(n);
 With details in Mezzadri 2007,
 "How to generate random matrices from the classical compact groups."

 To ensure matrices in SO(n), we permute the two first columns when
 the determinant is -1.

 See also: randskew qr_unique randunitary

0001 function Q = randrot(n, N)
0002 % Generates uniformly random rotation matrices.
0003 %
0004 % function Q = randrot(n, N)
0005 %
0006 % Q is an n-by-n-by-N array such that each slice Q(:, :, i) is a random
0007 % orthogonal matrix of size n of determinant +1 (i.e., a matrix in SO(n)),
0008 % sampled from the Haar measure (uniform distribution).
0009 %
0010 % By default, N = 1.
0011 %
0012 % Complexity: N times O(n^3).
0013 % Theory in Diaconis and Shahshahani 1987 for the uniformity on O(n);
0014 % With details in Mezzadri 2007,
0015 % "How to generate random matrices from the classical compact groups."
0016 %
0017 % To ensure matrices in SO(n), we permute the two first columns when
0018 % the determinant is -1.
0019 %
0020 % See also: randskew qr_unique randunitary
0021 
0022 % This file is part of Manopt: www.manopt.org.
0023 % Original author: Nicolas Boumal, Sept. 25, 2012.
0024 % Contributors:
0025 % Change log:
0026 %   June 18, 2019 (NB)
0027 %       Now generating all initial random matrices in one shot (which
0028 %       should be more efficient) and calling qr_unique.
0029 
0030 
0031     if nargin < 2
0032         N = 1;
0033     end
0034     
0035     if n == 1
0036         Q = ones(1, 1, N);
0037         return;
0038     end
0039     
0040     % Generated as such, Q is uniformly distributed over O(n): the group
0041     % of orthogonal matrices; see Mezzadri 2007.
0042     Q = qr_unique(randn(n, n, N));
0043     
0044     for k = 1 : N
0045         
0046         % If a slice of Q is in O(n) but not in SO(n), we permute its two
0047         % first columns to negate its determinant. This ensures the new
0048         % slice is in SO(n), uniformly distributed.
0049         if det(Q(:, :, k)) < 0
0050             Q(:, [1 2], k) = Q(:, [2 1], k);
0051         end
0052         
0053     end
0054 
0055 end