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'

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

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.

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.

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.”