How to generate a random unitary matrix perfectly in Julia

Specific Domains Numerics
1

Dear friends,

I am trying to implement a function that generates a random unitary matrix. However, I found that I met a precision issue. Is any idea how to solve this problem?

using LinearAlgebra

function RandomUnitaryMatrix(N::Int)
    x = (rand(N,N) + rand(N,N)*im) / sqrt(2)
    f = qr(x)
    diagR = sign.(real(diag(f.R)))
    diagR[diagR.==0] .= 1
    diagRm = diagm(diagR)
    u = f.Q * diagRm
    
    return u
end

u = RandomUnitaryMatrix(2)
u * u'

The output is

2×2 Array{Complex{Float64},2}:
         1.0+0.0im          1.66533e-16-1.11022e-16im
 1.66533e-16+1.11022e-16im          1.0+0.0im
2

Why isn’t what you have “good enough”?

It would pass a test at least:

using LinearAlgebra

X = [         1.0+0.0im          1.66533e-16-1.11022e-16im;
        1.66533e-16+1.11022e-16im          1.0+0.0im  ]

julia> isapprox(X, I)
true

Furthermore, I’m not sure “perfectly” is possible. If you use the RandomMatrices.jl package to simulate a real, unitary (orthogonal) matrix, you can get the same behavior:

using Random, RandomMatrices
Random.seed!(8675309)

d = Haar(1)
R = rand(d, 2)

julia> R * R'
2×2 Array{Float64,2}:
 1.0          5.55112e-17
 5.55112e-17  1.0
3 Likes
3

OK. :smile: Thank you very much.

4

There do exist Matrices that are exactly unitary. Bumping the entries in the example up or down using nextfloat or prevfloat may find one. But probably you don’t actually need this and it’s probably exponentially expensive with the dimension.

3 Likes
6

Sure. I don’t mean that you couldn’t ever get u * u' exactly equal to I . But writing a sampler where it’s guaranteed for every draw is a tall order, and as you say probably not necessary.

7

Isn’t the matrix found that way not exactly unitary, but only that it appears exactly unitary within Float64 calculations?

8

Correct. I don’t think there any rational points on the circle besides poles so the only true unitary matrices are diagonal

1 Like
9

Yeah, in that sense, “perfectly unitary” is indeed too much to ask for.

1 Like
10

There are actually infinitely many rational points on the circle, so there are non-diagonal “true” unitary matrices:

julia> a = 3//5
3//5

julia> b = 4//5
4//5

julia> u = [a -b; im*b im*a]
2×2 Array{Complex{Rational{Int64}},2}:
 3//5+0//1*im  -4//5+0//1*im
 0//1+4//5*im   0//1+3//5*im

julia> u * u'
2×2 Array{Complex{Rational{Int64}},2}:
 1//1+0//1*im  0//1+0//1*im
 0//1+0//1*im  1//1+0//1*im

julia> u * u' == I
true
4 Likes
11

But not in Float64 representation :wink: “The hypotenuse of a primitive Pythagorean triple is an odd integer larger than 1, as its square is the sum of an odd and an even number.”

2 Likes
12

If you want precision, maybe use a Rational{Int64} to calc the Random unitary matrix, and converting to float on demand

1 Like