Unitary matrices produced by `Stewart` do not seem to be distributed with Haar measure

[pohwanwuzan]

Hi there,

First of all thanks for this nice package. I had a problem when I tried to use the function Stewart to generate unitary matrices in my simulation though. The function rests in the file HaarMeasure.jl, so I assume that it's intended to produce unitary matrices distributed with Haar measure. According to [1], the normalized eigenvalue density is $\rho(\theta) = \frac{1}{2 \pi}$, where $\theta$ is the argument of the eigenvalues of the generated matrices.

I tried to reproduce Figure 2 in [1] with the following script:

using LinearAlgebra: eigvals
using Random: Xoshiro
using RandomMatrices: Stewart

rng = Xoshiro(1013)
T = ComplexF64

m = 50        # order of the matrix
N = 10000     # number of samples

theta = Vector{Float64}(undef, m * N)

get_theta(z) = imag(log(z))

# simulation
for i = 1:N

    V = Stewart(T, m)

    ev = eigvals(V)

    @views copy!(theta[((i - 1) * m + 1):(i * m)], get_theta.(ev))

end

# write to csv
open("./theta.csv", "w") do io

    println(io, "theta")

    for i = 1:length(theta)
        println(io, theta[i])
    end

end

The resulting density does not seem to be constant; instead it looks like the wrong distribution in Figure 1 in [1]:

Looking into the source code, it seems that the function calls larfg from LAPACK to generate Householder reflectors, but the reflectors themselves are not rotated as suggested in [1] (specifically equation (7.22)). To test out I just made a simple change to the source code:

function Stewart(::Type{ComplexF64}, n)
    τ = Array{ComplexF64}(undef, n)
    H = complex.(randn(n, n), randn(n,n))

    pτ = pointer(τ)
    pβ = pointer(H)
    pH = pointer(H, 2)

    r = one(ComplexF64)  # store the rotation

    for i = 0:n-2
        r *= - abs(H[i + 1, i + 1]) / H[i + 1, i + 1]   # update the rotation
        larfg!(n - i, pβ + (n + 1)*16i, pH + (n + 1)*16i, 1, pτ + 16i)
    end

    return (LinearAlgebra.QRPackedQ(H, τ), r)

end

Then I modified the script and carried out the small experiment again:

using LinearAlgebra: eigvals
using Random: Xoshiro
using RandomMatrices: Stewart

rng = Xoshiro(1013)
T = ComplexF64

m = 50        # order of the matrix
N = 10000     # number of samples

theta = Vector{Float64}(undef, m * N)
U = Matrix{T}(undef, m, m)

get_theta(z) = imag(log(z))

# simulation
for i = 1:N

    V, r = Stewart(T, m)
    copy!(U, V)

    @. U = r * U

    ev = eigvals(U)

    @views copy!(theta[((i - 1) * m + 1):(i * m)], get_theta.(ev))

end

# write to csv
open("./theta.csv", "w") do io

    println(io, "theta")

    for i = 1:length(theta)
        println(io, theta[i])
    end

end

The density now appears constant:

Perhaps I'm missing something here, and it would be great to hear your opinion. Cheers!

References

[1] Mezzadri, F., 2006. How to generate random matrices from the classical compact groups. arXiv preprint math-ph/0609050.