# How to generate a random unitary matrix perfectly in Julia

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
``````

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

OK. Thank you very much.

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

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.

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

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

1 Like

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

1 Like

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

But not in Float64 representation “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

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

1 Like