I was wondering if there was already a package out there for generating random matrices from SU(N) and SO(N). This mostly boils down to generating random points on n-spheres and is not difficult to implement for particular cases, but I was wondering if someone had already done this as it would be nice to have a more efficient and properly parallelized implementation than I am likely to write.
RandomMatrices.jl lets you sample from Haar measure.
This seems to have quite a lot of deprecation warnings at the moment. That aside, my understanding of this is that this can generate me random U(N) matrices… the only way I can think off the top of my head of producing SU(N) matrices from these (essentially dividing out the determinant) seems rather inefficient to me (not to mention the fact that I’m then essentially generating an extra S^{1} variable only to divide it out again in a reall complicated way!)… Obviously this isn’t really my area of expertise, so perhaps I’m missing something stupidly obvious.
I suppose Q^2 probably has the right distribution, though it requires a dense matrix multiplication, as far as I know.
@alanedelman might have a better suggestion.
I’ll fix the deprecation warnings soon (I thought I already had…).
Hmm, having looked at the code HaarMeasure.jl the default version just calls QR on a Ginibre (random Gaussian) matrix. If there is a qrfact call that uses Given’s instead of Householder, then it should give you the equivalent for SO(N).
There is also an O(N^2) implementation of Stewart’s algorithm. Perhaps that can be modified?
Maybe you can pullback your sampling distribution to its Lie Algebra which is quite simple (? pun intended).
Actually, I just realized something that makes me incredibly happy (though, bear in mind, I don’t know how efficient this is, probably not very). I can just use expm of the generators. Why does this make me very happy? Because look how beautifully this can be written in Julia. For example, for SU(2), after defining the Pauli matrices and a random unit vector n on \mathbb{R}^{N^{2}-1} and a random scalar in [0, 2\pi)
σ⃗ = [σ₁, σ₂, σ₃]
U = expm(im*φ*n'*σ⃗)
How cool is that? Unfortunately this probably involves dense matrix multiplication, so despite looking beautiful it’s probably incredibly inefficient, so alternatively I can do
\cos(\phi) + i(n_{a}T_{a})\sin(\phi)
where T_{a} are the generators. This should involve only scalar multiplication and dense matrix addition. These are both essentially what @rveltz suggested.
4 Likes