# How to generate a random unitary matrix perfectly in Julia

**URL:** <https://discourse.julialang.org/t/how-to-generate-a-random-unitary-matrix-perfectly-in-julia/34102>\
**Category:** Numerics\
**Created:** [February 3, 2020, 3:45am UTC](https://discourse.julialang.org/t/how-to-generate-a-random-unitary-matrix-perfectly-in-julia/34102 "2020-02-03T03:45:31Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![JunjieWu](https://avatars.discourse-cdn.com/v4/letter/j/aca169/32.png) [@JunjieWu](https://discourse.julialang.org/u/JunjieWu)\
**Post date:** [February 3, 2020, 3:45am UTC](https://discourse.julialang.org/t/how-to-generate-a-random-unitary-matrix-perfectly-in-julia/34102/1 "2020-02-03T03:45:31Z")

</div>

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?

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

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

```

---

<div class="post-metadata">

**Author:** ![johnczito](https://avatars.discourse-cdn.com/v4/letter/j/53a042/32.png) [@johnczito](https://discourse.julialang.org/u/johnczito)\
**Post date:** [February 3, 2020, 5:13am UTC](https://discourse.julialang.org/t/how-to-generate-a-random-unitary-matrix-perfectly-in-julia/34102/2 "2020-02-03T05:13:52Z")

</div>

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

It would pass a test at least:

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

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

```

---

<div class="post-metadata">

**Author:** ![JunjieWu](https://avatars.discourse-cdn.com/v4/letter/j/aca169/32.png) [@JunjieWu](https://discourse.julialang.org/u/JunjieWu)\
**Post date:** [February 3, 2020, 5:15am UTC](https://discourse.julialang.org/t/how-to-generate-a-random-unitary-matrix-perfectly-in-julia/34102/3 "2020-02-03T05:15:23Z")

</div>

OK. 😄 Thank you very much.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [February 3, 2020, 7:41am UTC](https://discourse.julialang.org/t/how-to-generate-a-random-unitary-matrix-perfectly-in-julia/34102/4 "2020-02-03T07:41:59Z")

</div>

> [@johnczito](#):
>
> I’m not sure “perfectly” is possible

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.

---

<div class="post-metadata">

**Author:** ![johnczito](https://avatars.discourse-cdn.com/v4/letter/j/53a042/32.png) [@johnczito](https://discourse.julialang.org/u/johnczito)\
**Post date:** [February 3, 2020, 9:30am UTC](https://discourse.julialang.org/t/how-to-generate-a-random-unitary-matrix-perfectly-in-julia/34102/6 "2020-02-03T09:30:17Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [February 3, 2020, 10:05am UTC](https://discourse.julialang.org/t/how-to-generate-a-random-unitary-matrix-perfectly-in-julia/34102/7 "2020-02-03T10:05:58Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [February 3, 2020, 10:22am UTC](https://discourse.julialang.org/t/how-to-generate-a-random-unitary-matrix-perfectly-in-julia/34102/8 "2020-02-03T10:22:40Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [February 3, 2020, 11:06am UTC](https://discourse.julialang.org/t/how-to-generate-a-random-unitary-matrix-perfectly-in-julia/34102/9 "2020-02-03T11:06:38Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![thofma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thofma/32/1691_2.png) [@thofma](https://discourse.julialang.org/u/thofma)\
**Post date:** [February 3, 2020, 11:57am UTC](https://discourse.julialang.org/t/how-to-generate-a-random-unitary-matrix-perfectly-in-julia/34102/10 "2020-02-03T11:57:26Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [February 3, 2020, 12:47pm UTC](https://discourse.julialang.org/t/how-to-generate-a-random-unitary-matrix-perfectly-in-julia/34102/11 "2020-02-03T12:47:07Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [February 3, 2020, 4:18pm UTC](https://discourse.julialang.org/t/how-to-generate-a-random-unitary-matrix-perfectly-in-julia/34102/12 "2020-02-03T16:18:00Z")

</div>

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