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.

[araujoms]

I'm not a maintainer of this package, but I have recently been studying the same problem in order to sample random unitaries for my own package, and I believe I understand the issue.

First of all, in order to correct the sign of the matrix, you need to transform it from a bunch of Householder reflections to a regular matrix, but this transformation alone is O(d^3), which kind of defeats the purpose of Stewart's algorithm in offering a O(d^2) method. In his paper he proposes storing the unitary as the set of householder reflections together with the diagonal matrix to correct the signs. Which of course works, but some type other than QRPackedQ would be needed. Nevertheless, Stewart's algorithm is still much faster than doing the full QR decomposition, so it's still worth using it even without going for exotic types.

Secondly, I don't think what you are doing is correct. You can't correct the sign via a scalar. And furthermore Julia's qr makes the diagonal of R real even for complex matrices, so there's no need to use complex phases, +-1 suffices.

[dlfivefifty]

Note this package doesn’t really have an active maintainer so if you wish to contribute and eventually take over as maintainer it would be really appreciated

[araujoms]

I appreciate the vote of confidence, but I'm afraid this package needs more love than I can give. I can fix the implementation of Stewart's algorithm, but I have neither the time nor the knowledge needed for something deeper.

Maybe it's a good idea to advertise the package in JuliaCon? There must be plenty of people familiar with random matrix theory there, and I don't think this is a case where we need to be afraid of Jia Tan.