implement random_unitary

[orialb]

Needed for https://github.com/ITensor/ITensors.jl/issues/377 (and useful also for other use-cases). I implemented it in such a way that random_unitary(n,n) will sample according to the Haar measure (following https://arxiv.org/abs/math-ph/0609050, page 11).

[emstoudenmire]

Looks good! Thanks for adding

[mtfishman]

For the sake of generic code, as well as how it will be used in ITensors.jl for randomMPS, should we combine random_orthog and random_unitary into one function random_unitary(::Type{<:Number}, ::Int, ::Int) where you can specify the type? This would allow it to be used for other types like single-precision as well. I think also it should default to Float64, since we basically always default to Float64.

We could always make an alias random_orthog(::Type{ElT}, n::Int, m::Int) where {ElT <: Real} = random_unitary(ElT, n, m).

[orialb]

I think having a function called random_unitary which defaults to returning only matrices in O(n) (say in the square case) can be quite surprising. I would expect it to draw uniformly from U(n). So in my opinion if we combine these two functions, random_unitary should return complex matrices by default when an eltype is not specified.

Having the alias would be also good in my opinion, as having random_orthog might be clearer in some cases where one is explicitly restricting to orthogonal matrices.

[mtfishman]

Yeah that is a good point, from a random matrix theory perspective it could be surprising. So let's have:

random_unitary(::Type{ElT}, n, m) where {ElT <: Number}
random_unitary(n, m) = random_unitary(ComplexF64, n, m)
random_orthog(::Type{ElT}, n, m) where {ElT <: Real} = random_unitary(ElT, n, m)
random_orthog(n, m) = random_orthog(Float64, n, m)