Skip to content

Commit 1393f70

Browse files
committed
src/sage/matrix/special.py: new random_unitary_matrix() function
Add a new function to generate random unitary matrices based on the paper "The Generation of All Rational Orthogonal Matrices" by Liebeck and Osborne (1991).
1 parent 08c0989 commit 1393f70

File tree

1 file changed

+166
-0
lines changed

1 file changed

+166
-0
lines changed

src/sage/matrix/special.py

Lines changed: 166 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@
4040
:meth:`~sage.matrix.special.random_rref_matrix`
4141
:meth:`~sage.matrix.special.random_subspaces_matrix`
4242
:meth:`~sage.matrix.special.random_unimodular_matrix`
43+
:meth:`~sage.matrix.special.random_unitary_matrix`
4344
:meth:`~sage.matrix.special.toeplitz`
4445
:meth:`~sage.matrix.special.vandermonde`
4546
:meth:`~sage.matrix.special.vector_on_axis_rotation_matrix`
@@ -241,6 +242,9 @@ def random_matrix(ring, nrows, ncols=None, algorithm='randomize', implementation
241242
242243
- ``'unimodular'`` -- creates a matrix of determinant 1
243244
245+
- ``'unitary'`` -- creates a (square) unitary matrix over a
246+
subfield of the complex numbers
247+
244248
- ``'diagonalizable'`` -- creates a diagonalizable matrix. if the
245249
base ring is ``QQ`` creates a diagonalizable matrix whose eigenvectors,
246250
if computed by hand, will have only integer entries. See the
@@ -665,6 +669,8 @@ def random_matrix(ring, nrows, ncols=None, algorithm='randomize', implementation
665669
return random_subspaces_matrix(parent, *args, **kwds)
666670
elif algorithm == 'unimodular':
667671
return random_unimodular_matrix(parent, *args, **kwds)
672+
elif algorithm == 'unitary':
673+
return random_unitary_matrix(parent, *args, **kwds)
668674
else:
669675
raise ValueError('random matrix algorithm "%s" is not recognized' % algorithm)
670676

@@ -3056,6 +3062,166 @@ def random_unimodular_matrix(parent, upper_bound=None, max_tries=100):
30563062
return random_matrix(ring, size,algorithm='echelonizable',rank=size, upper_bound=upper_bound, max_tries=max_tries)
30573063

30583064

3065+
@matrix_method
3066+
def random_unitary_matrix(parent):
3067+
r"""
3068+
Generate a random (square) unitary matrix of a given size
3069+
with entries in a subfield of the complex numbers.
3070+
3071+
INPUT:
3072+
3073+
- ``parent`` -- :class:`sage.matrix.matrix_space.MatrixSpace`; a
3074+
square matrix space over a subfield of the complex numbers
3075+
having characteristic zero.
3076+
3077+
OUTPUT:
3078+
3079+
A random unitary matrix in ``parent``. A :exc:`ValueError` is
3080+
raised if ``parent`` is not an appropriate matrix space (is not
3081+
square, or is not over a usable field).
3082+
3083+
ALGORITHM:
3084+
3085+
The method described by Liebeck and Osborne [LieOs1991]_ is used
3086+
almost verbatim.
3087+
3088+
.. WARNING:
3089+
3090+
Inexact rings may violate your expectations; in particular,
3091+
the rings ``RR`` and ``CC`` are accepted by this method but
3092+
the resulting matrix will usually fail the ``is_unitary()``
3093+
check due to numerical issues.
3094+
3095+
EXAMPLES:
3096+
3097+
This function is not in the global namespace and must be imported
3098+
before being used::
3099+
3100+
sage: from sage.matrix.constructor import random_unitary_matrix
3101+
3102+
Matrices with rational entries::
3103+
3104+
sage: n = ZZ.random_element(10)
3105+
sage: MS = MatrixSpace(QQ, n)
3106+
sage: U = random_unitary_matrix(MS)
3107+
sage: U.is_unitary() and U in MS
3108+
True
3109+
3110+
Matrices over a quadraric field::
3111+
3112+
sage: n = ZZ.random_element(10)
3113+
sage: K = QuadraticField(-1,'i')
3114+
sage: MS = MatrixSpace(K, n)
3115+
sage: U = random_unitary_matrix(MS)
3116+
sage: U.is_unitary() and U in MS
3117+
True
3118+
3119+
Matrices with entries in the algebraic real field (slow)::
3120+
3121+
sage: # long time
3122+
sage: n = ZZ.random_element(4)
3123+
sage: MS = MatrixSpace(AA, n)
3124+
sage: U = random_unitary_matrix(MS)
3125+
sage: U.is_unitary() and U in MS
3126+
True
3127+
3128+
Matrices with entries in the algebraic field (slower yet)::
3129+
3130+
sage: # long time
3131+
sage: n = ZZ.random_element(2)
3132+
sage: MS = MatrixSpace(QQbar, n)
3133+
sage: U = random_unitary_matrix(MS)
3134+
sage: U.is_unitary() and U in MS
3135+
True
3136+
3137+
Double-precision real/complex matrices can be constructed as
3138+
well::
3139+
3140+
sage: n = ZZ.random_element(10)
3141+
sage: MS = MatrixSpace(RDF, n)
3142+
sage: U = random_unitary_matrix(MS)
3143+
sage: U.is_unitary() and U in MS
3144+
True
3145+
sage: MS = MatrixSpace(CDF, n)
3146+
sage: U = random_unitary_matrix(MS)
3147+
sage: U.is_unitary() and U in MS
3148+
True
3149+
3150+
TESTS:
3151+
3152+
We raise a :exc:`ValueError` if the supplied matrix space is not
3153+
square::
3154+
3155+
sage: MS = MatrixSpace(QQ, 3, 4)
3156+
sage: random_unitary_matrix(MS)
3157+
Traceback (most recent call last):
3158+
...
3159+
ValueError: parent must be square
3160+
3161+
Likewise if its base ring is not a field::
3162+
3163+
sage: MS = MatrixSpace(ZZ, 3)
3164+
sage: random_unitary_matrix(MS)
3165+
Traceback (most recent call last):
3166+
...
3167+
ValueError: base ring of parent must be a field
3168+
3169+
Likewise if that field is not of characteristic zero::
3170+
3171+
sage: MS = MatrixSpace(GF(5), 3)
3172+
sage: random_unitary_matrix(MS)
3173+
Traceback (most recent call last):
3174+
...
3175+
ValueError: base ring of parent must have characteristic zero
3176+
3177+
Likewise if that field is not a subfield of the complex numbers::
3178+
3179+
sage: R = Qp(7)
3180+
sage: R.characteristic()
3181+
0
3182+
sage: MS = MatrixSpace(R, 3)
3183+
sage: random_unitary_matrix(MS)
3184+
Traceback (most recent call last):
3185+
...
3186+
ValueError: base ring of parent must be a subfield of the complex
3187+
numbers
3188+
"""
3189+
n = parent.nrows()
3190+
if n != parent.ncols():
3191+
raise ValueError("parent must be square")
3192+
3193+
F = parent.base_ring()
3194+
if not F.is_field():
3195+
raise ValueError("base ring of parent must be a field")
3196+
elif not F.characteristic().is_zero():
3197+
# This is probably covered by checking for subfields of
3198+
# RLF/CLF, but it's mentioned explicitly in the paper and my
3199+
# faith in the subfield check is not 100%, so we verify the
3200+
# characteristic separately.
3201+
raise ValueError("base ring of parent must have characteristic zero")
3202+
3203+
from sage.rings.real_lazy import RLF, CLF
3204+
if not (RLF.has_coerce_map_from(F) or
3205+
F.has_coerce_map_from(RLF) or
3206+
CLF.has_coerce_map_from(F) or
3207+
F.has_coerce_map_from(CLF)):
3208+
# The implementation of SR.random_element() currently just
3209+
# returns a random integer coerced into SR, so there is no
3210+
# benefit to allowing SR here when QQ is available.
3211+
raise ValueError("base ring of parent must be a subfield "
3212+
"of the complex numbers")
3213+
3214+
I = identity_matrix(F,n)
3215+
A = random_matrix(F,n)
3216+
S = A - A.conjugate_transpose()
3217+
U = (S-I).inverse()*(S+I)
3218+
D = identity_matrix(F,n)
3219+
for i in range(n):
3220+
if ZZ.random_element(2).is_zero():
3221+
D[i,i] *= F(-1)
3222+
return D*U
3223+
3224+
30593225
@matrix_method
30603226
def random_diagonalizable_matrix(parent, eigenvalues=None, dimensions=None):
30613227
"""

0 commit comments

Comments
 (0)