For orthogonal matrices, I use this procedure:
function random_orthogonal(n, θ)
X = randn(n, n)
Y = X - X'
Y *= sqrt(2) / frobenius_norm(Y)
return expm(θ*Y)
end
For unitary matrices, this is what I tried:
function random_unitary(n, θ)
X = randn(n, n) + rand(n, n)*im
Y = X - X'
Y *= sqrt(2) / frobenius_norm(Y)
return expm(θ*Y)
end
(' is the conjugate transpose function.)
If we let $\theta$ be a random value from $[0, 2 \pi]$, then random_orthogonal produces the expected distribution of eigenvalues. However, random_unitary does not. What am I doing wrong?
(I realize this isn't the most computationally efficient way of doing this; I'm just trying to understand it from a theoretical perspective.)