Computing Haar measure of matrices sampled from SO(n)

Question

I am looking to sample uniform matrices from SO(n).

I know that uniform matrices can be sampled from O(n) by taking the QR decomposition of Gaussian random square matrices and adjusting the sign of the diagonal values (e.g. as found here and in references therein). The uniformity of these matrices can then be tested by assessing the angle and spread of the eigenvalues of the matrices obtained (full code here):

[e.g.

import numpy as np
# to calculate Haar measure for a given matrix `input`
B,_ = np.linalg.eig(input)
evals = np.sort(np.angle(B))
seps = B.shape[0] / (2*np.pi) * np.diff(evals,axis=0) 

]

I note that the distribution of eigenvalues is approximately uniform in the interval $(-\pi, \pi]$ - except with many matrices having eigenvalues of $\pm1$ - is this expected?

Is it possible to uniformly sample from SO(n) by taking matrices sampled from O(n) and forcing the determinant to be +1? (e.g. sample a matrix from O(n); then, if the determinant is -1, simply multiply the first column by -1).

If so, what is the equivalent Haar measure for this group? Using the same computation as above changes the distribution of the eigenvalues significantly:

Given that eigenvalues of a should occur in $\lfloor \frac{n}{2} \rfloor$ complex conjugate pairs (with potentially one eigenvalue of +1, if n is odd), one possible alternative is to define the Haar measure using the argument of the eigenvalues in the interval $[0, \frac{\pi}{\lfloor \frac{n}{2} \rfloor})$. This does indeed create a uniform spread but I am not sure about the correctness of this approach.

Answer 1

Indeed, the distribution function of the eigenphases of a random matrix in $\operatorname{SO}(n)$ has a peak at 0 and at $\pm\pi$. It only becomes uniform for large $n$. The joint distribution function for $n=2m$ even is given by (Girko, 1985) $$p(\theta_1, \cdots, \theta_m) = C \prod_{1 \leq k < j \leq m} (\cos\theta_k - \cos\theta_j)^2~.$$ As a simple example, for $n=4$ this gives upon integration over $\theta_2$ at fixed $\theta_1$ the density profile $$p(\theta)=\frac{\cos 2 \theta+2}{4 \pi }.$$ If you sample from $O(n)$, you can either discard the matrices with determinant $-1$, or multiply the first column by $-1$.

---

By way of illustration, this is the density profile for $n=100$ (computed by averaging over 500 random matrices), when the peaks at $0$ and $\pm\pi$ have become very small.

The reason why the Haar measure for SO$(n)$ does not give a rotationally invariant density of eigenvalues on the unit circle is that the eigenvalues must come in complex conjugate pairs $e^{\pm i\theta}$. This constraint is not rotationally invariant.

Answer 2

There's an asymmetry because members of O(n) include rotation with random angle and flips with random axis, for small $n$, the fraction of flips is large, which create peaks at $-1$ and $1$ peaks. Restricting to $SO(n)$ removes the peaks at $-1$, and $1$ remaining density follows cosine formula from Girko

negateFirstRow[A_] := {-A[[1, All]]}~Join~A[[2 ;;, All]];
sampleO[n_] := RandomVariate[CircularRealMatrixDistribution[n]];
sampleSO[n_] := 
  With[{mat = sampleO@n}, If[Det[mat] > 0, mat, negateFirstRow@mat]];

n = 4;
numSamples = 100000;
label = StringForm["Eigenvalues of O(``)", n];
angles = Arg /@ Flatten[Table[Eigenvalues@sampleO@n, {numSamples}]];
Histogram[angles, Automatic, PDF, PlotLabel -> label, 
 AxesLabel -> {"arg", "density"}]

label = StringForm["Eigenvalues of SO(``)", n];
angles = Arg /@ Flatten[Table[Eigenvalues@sampleSO@n, {numSamples}]];
observedPlot = 
  Histogram[angles, {Pi/15}, PDF, PlotLabel -> label, 
   AxesLabel -> {"arg", "density"}];
predictedPlot = 
  Plot[(Cos[2 t] + 2)/(4 Pi), {t, -Pi, Pi}, 
   PlotRange -> {0, 3/(4 Pi)}, 
   PlotLegends -> {TraditionalForm[(Cos[2 t] + 2)/(4 Pi)]}];
Show[observedPlot, predictedPlot]