The simplest method is based on the QR decomposition of a random complex matrix with independent elements following a Normal distribution.
In Mathematica one can generate a random complex number with Normal distribution for both the real and imaginary parts with the following function
NormalRandomComplex[] := (#[[1]] + #[[2]]*I) &@
RandomReal[MultinormalDistribution[{0, 0}, {{1, 0}, {0, 1}}]]
where {0, 0} specifies the origin and {{1, 0}, {0, 1}} the covariance matrix (standard deviation 1)
The function that generates the random n x n unitary matrix is
RandomUnitaryMatrix[n_] := Module[{z, q, r},
z = Table[NormalRandomComplex[], {n}, {n}];
{q, r} = QRDecomposition[z];
q.DiagonalMatrix[Sign[Tr[r, List]]]
] z = Table[NormalRandomComplex[], {n}, {n}];
{q, r} = QRDecomposition[z];
q.DiagonalMatrix[Sign[Tr[r, List]]]
In this code, the QR decomposition of the complex random matrix with normal distribution outputs a unitary matrix "q" and a triangular matrix "r". The sought random unitary matrix is a correction of "q", which must be multiplied by a diagonal matrix filled with 1 and -1 according to the sign of the corresponding diagonal element of "r".
No comments:
Post a Comment