Matrices in Julia with Orthogonality constraints

Question

How do you construct random unitary matrices in Julia

An Orthogonal matrix is defined as $QQ^T=I_n=Q^TQ$ with $Q \in \mathbb{R}^{n \times n}$. This can be easily found in Julia from Manifolds.jl by the following command:

Q=rand((OrthogonalMatrices(n))).

Likewise, a unitary matrix is defined as $QQ^{} = I_n = Q^{} Q$ with * being the complex conjugate transposed. In Manifolds.jl, they have a function called UnitaryMatrices(n). However, when I write the following command: Q=rand((UnitaryMatrices(n))).

I get the following error message: "and!(::ManifoldsBase.AbstractManifold, ::Matrix{ComplexF64}) is ambiguous."

Answer 1

Although I didn't find how to avoid the error from the Manifolds package, I did find a way to create the matrix you wished with a different package: RandomMatrices.

The code looks like this:

julia> using LinearAlgebra, RandomMatrices

julia> const n = 4
4

julia> M = rand(RandomMatrices.Haar(2), n)
4×4 Matrix{ComplexF64}:
 0.0491609+0.547619im    0.337504+0.0588388im   0.107132-0.635825im     -0.252285+0.317691im
   0.37847+0.212615im    0.252793-0.512421im    0.558867+0.412853im    0.00980966+0.0468962im
 -0.257792+0.510858im   -0.327271-0.00628619im   0.10582-0.00169264im    0.742162+0.0584425im
  0.424879+0.0360062im  -0.324348+0.586346im    0.146708+0.262194im    -0.0767287+0.522515im

julia> M * M'
4×4 Matrix{ComplexF64}:
          1.0+0.0im           1.54646e-16-4.51345e-17im  -1.81615e-16-1.63497e-16im  -8.33993e-17+6.35032e-17im
  1.54646e-16+4.51345e-17im           1.0+0.0im          -6.82757e-18+5.95489e-17im  -2.02409e-16-7.38223e-17im
 -1.81615e-16+1.63497e-16im  -6.82757e-18-5.95489e-17im           1.0+0.0im           4.49412e-18+2.613e-16im
 -8.33993e-17-6.35032e-17im  -2.02409e-16+7.38223e-17im   4.49412e-18-2.613e-16im             1.0+0.0im

julia> M * M' ≈ I
true

The "magic" parameter 2 to Haar is necessary (1 gives an orthogonal matrix). This package may have some installation issues with other packages (at least in the environment I've tested it).

Hope this helps and maybe someone else will figure out how to use Manifolds for this.

Answer 2

Generate a random complex matrix, skew it, then apply the Cayley Transform.

Here is a Julia one-liner

n=5;  C = rand(ComplexF64,n,n);  Q = (I-C+C') / (I+C-C');  Q'Q ≈ I
true