From cd7bc4d83a4e06df55980a6177c2fabe383f7b51 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Sun, 6 Jul 2025 15:24:26 -0400 Subject: [PATCH 1/4] src/sage/matrix/special.py: add some "# long time" tags I'm getting a few warnings about tests in this file that take longer than 5s to run (without --long). Here we add "# long time" tags. --- src/sage/matrix/special.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/sage/matrix/special.py b/src/sage/matrix/special.py index f1de08e87d5..fe98679a855 100644 --- a/src/sage/matrix/special.py +++ b/src/sage/matrix/special.py @@ -303,7 +303,7 @@ def random_matrix(ring, nrows, ncols=None, algorithm='randomize', implementation sage: expected(100) 1/25250 sage: add_samples(ZZ, 5, 5) - sage: while not all(abs(dic[a]/total_count - expected(a)) < 0.001 for a in dic): + sage: while not all(abs(dic[a]/total_count - expected(a)) < 0.001 for a in dic): # long time ....: add_samples(ZZ, 5, 5) The ``distribution`` keyword set to ``uniform`` will limit values @@ -313,7 +313,7 @@ def random_matrix(ring, nrows, ncols=None, algorithm='randomize', implementation sage: total_count = 0 sage: dic = defaultdict(Integer) sage: add_samples(ZZ, 5, 5, distribution='uniform') - sage: while not all(abs(dic[a]/total_count - expected(a)) < 0.001 for a in dic): + sage: while not all(abs(dic[a]/total_count - expected(a)) < 0.001 for a in dic): # long time ....: add_samples(ZZ, 5, 5, distribution='uniform') The ``x`` and ``y`` keywords can be used to distribute entries uniformly. @@ -324,14 +324,14 @@ def random_matrix(ring, nrows, ncols=None, algorithm='randomize', implementation sage: total_count = 0 sage: dic = defaultdict(Integer) sage: add_samples(ZZ, 4, 8, x=70, y=100) - sage: while not all(abs(dic[a]/total_count - expected(a)) < 0.001 for a in dic): + sage: while not all(abs(dic[a]/total_count - expected(a)) < 0.001 for a in dic): # long time ....: add_samples(ZZ, 4, 8, x=70, y=100) sage: expected = lambda n : 1/10 if n in range(-5, 5) else 0 sage: total_count = 0 sage: dic = defaultdict(Integer) sage: add_samples(ZZ, 3, 7, x=-5, y=5) - sage: while not all(abs(dic[a]/total_count - expected(a)) < 0.001 for a in dic): + sage: while not all(abs(dic[a]/total_count - expected(a)) < 0.001 for a in dic): # long time ....: add_samples(ZZ, 3, 7, x=-5, y=5) If only ``x`` is given, then it is used as the upper bound of a range @@ -341,7 +341,7 @@ def random_matrix(ring, nrows, ncols=None, algorithm='randomize', implementation sage: total_count = 0 sage: dic = defaultdict(Integer) sage: add_samples(ZZ, 5, 5, x=25) - sage: while not all(abs(dic[a]/total_count - expected(a)) < 0.001 for a in dic): + sage: while not all(abs(dic[a]/total_count - expected(a)) < 0.001 for a in dic): # long time ....: add_samples(ZZ, 5, 5, x=25) To control the number of nonzero entries, use the ``density`` keyword From 35c6c2857b9e2d9a0e3caaebd00960b6d2b334ab Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Sun, 6 Jul 2025 20:39:21 -0400 Subject: [PATCH 2/4] src/doc/en/reference/references: add the LieOs1991 reference New reference (Liebeck and Osborne, 1991) to be used in a new method for generating random unitary matrices. --- src/doc/en/reference/references/index.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/doc/en/reference/references/index.rst b/src/doc/en/reference/references/index.rst index af0c2e7afad..ae8868e6a85 100644 --- a/src/doc/en/reference/references/index.rst +++ b/src/doc/en/reference/references/index.rst @@ -4485,6 +4485,10 @@ REFERENCES: Mathematics. Springer-Verlag, New York, 1997. ISBN 0-387-98254-X +.. [LieOs1991] Hans Liebeck and Anthony Osborne. + The Generation of All Rational Orthogonal Matrices. + The American Mathematical Monthly, 98(2):131-133, 1991. + .. [Lim] \C. H. Lim, *CRYPTON: A New 128-bit Block Cipher*; available at https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=89a685057c2ce4f632b105dcaec264b9162a1800 From 1ceab98239bde9c1fbd6d0f16a9ab04e3703d766 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Sun, 6 Jul 2025 20:46:23 -0400 Subject: [PATCH 3/4] 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). --- src/sage/matrix/special.py | 166 +++++++++++++++++++++++++++++++++++++ 1 file changed, 166 insertions(+) diff --git a/src/sage/matrix/special.py b/src/sage/matrix/special.py index fe98679a855..d6229b18928 100644 --- a/src/sage/matrix/special.py +++ b/src/sage/matrix/special.py @@ -40,6 +40,7 @@ :meth:`~sage.matrix.special.random_rref_matrix` :meth:`~sage.matrix.special.random_subspaces_matrix` :meth:`~sage.matrix.special.random_unimodular_matrix` + :meth:`~sage.matrix.special.random_unitary_matrix` :meth:`~sage.matrix.special.toeplitz` :meth:`~sage.matrix.special.vandermonde` :meth:`~sage.matrix.special.vector_on_axis_rotation_matrix` @@ -241,6 +242,9 @@ def random_matrix(ring, nrows, ncols=None, algorithm='randomize', implementation - ``'unimodular'`` -- creates a matrix of determinant 1 + - ``'unitary'`` -- creates a (square) unitary matrix over a + subfield of the complex numbers + - ``'diagonalizable'`` -- creates a diagonalizable matrix. if the base ring is ``QQ`` creates a diagonalizable matrix whose eigenvectors, 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 return random_subspaces_matrix(parent, *args, **kwds) elif algorithm == 'unimodular': return random_unimodular_matrix(parent, *args, **kwds) + elif algorithm == 'unitary': + return random_unitary_matrix(parent, *args, **kwds) else: raise ValueError('random matrix algorithm "%s" is not recognized' % algorithm) @@ -3056,6 +3062,166 @@ def random_unimodular_matrix(parent, upper_bound=None, max_tries=100): return random_matrix(ring, size,algorithm='echelonizable',rank=size, upper_bound=upper_bound, max_tries=max_tries) +@matrix_method +def random_unitary_matrix(parent): + r""" + Generate a random (square) unitary matrix of a given size + with entries in a subfield of the complex numbers. + + INPUT: + + - ``parent`` -- :class:`sage.matrix.matrix_space.MatrixSpace`; a + square matrix space over a subfield of the complex numbers + having characteristic zero. + + OUTPUT: + + A random unitary matrix in ``parent``. A :exc:`ValueError` is + raised if ``parent`` is not an appropriate matrix space (is not + square, or is not over a usable field). + + ALGORITHM: + + The method described by Liebeck and Osborne [LieOs1991]_ is used + almost verbatim. + + .. WARNING: + + Inexact rings may violate your expectations; in particular, + the rings ``RR`` and ``CC`` are accepted by this method but + the resulting matrix will usually fail the ``is_unitary()`` + check due to numerical issues. + + EXAMPLES: + + This function is not in the global namespace and must be imported + before being used:: + + sage: from sage.matrix.constructor import random_unitary_matrix + + Matrices with rational entries:: + + sage: n = ZZ.random_element(10) + sage: MS = MatrixSpace(QQ, n) + sage: U = random_unitary_matrix(MS) + sage: U.is_unitary() and U in MS + True + + Matrices over a quadraric field:: + + sage: n = ZZ.random_element(10) + sage: K = QuadraticField(-1,'i') + sage: MS = MatrixSpace(K, n) + sage: U = random_unitary_matrix(MS) + sage: U.is_unitary() and U in MS + True + + Matrices with entries in the algebraic real field (slow):: + + sage: # long time + sage: n = ZZ.random_element(4) + sage: MS = MatrixSpace(AA, n) + sage: U = random_unitary_matrix(MS) + sage: U.is_unitary() and U in MS + True + + Matrices with entries in the algebraic field (slower yet):: + + sage: # long time + sage: n = ZZ.random_element(2) + sage: MS = MatrixSpace(QQbar, n) + sage: U = random_unitary_matrix(MS) + sage: U.is_unitary() and U in MS + True + + Double-precision real/complex matrices can be constructed as + well:: + + sage: n = ZZ.random_element(10) + sage: MS = MatrixSpace(RDF, n) + sage: U = random_unitary_matrix(MS) + sage: U.is_unitary() and U in MS + True + sage: MS = MatrixSpace(CDF, n) + sage: U = random_unitary_matrix(MS) + sage: U.is_unitary() and U in MS + True + + TESTS: + + We raise a :exc:`ValueError` if the supplied matrix space is not + square:: + + sage: MS = MatrixSpace(QQ, 3, 4) + sage: random_unitary_matrix(MS) + Traceback (most recent call last): + ... + ValueError: parent must be square + + Likewise if its base ring is not a field:: + + sage: MS = MatrixSpace(ZZ, 3) + sage: random_unitary_matrix(MS) + Traceback (most recent call last): + ... + ValueError: base ring of parent must be a field + + Likewise if that field is not of characteristic zero:: + + sage: MS = MatrixSpace(GF(5), 3) + sage: random_unitary_matrix(MS) + Traceback (most recent call last): + ... + ValueError: base ring of parent must have characteristic zero + + Likewise if that field is not a subfield of the complex numbers:: + + sage: R = Qp(7) + sage: R.characteristic() + 0 + sage: MS = MatrixSpace(R, 3) + sage: random_unitary_matrix(MS) + Traceback (most recent call last): + ... + ValueError: base ring of parent must be a subfield of the complex + numbers + """ + n = parent.nrows() + if n != parent.ncols(): + raise ValueError("parent must be square") + + F = parent.base_ring() + if not F.is_field(): + raise ValueError("base ring of parent must be a field") + elif not F.characteristic().is_zero(): + # This is probably covered by checking for subfields of + # RLF/CLF, but it's mentioned explicitly in the paper and my + # faith in the subfield check is not 100%, so we verify the + # characteristic separately. + raise ValueError("base ring of parent must have characteristic zero") + + from sage.rings.real_lazy import RLF, CLF + if not (RLF.has_coerce_map_from(F) or + F.has_coerce_map_from(RLF) or + CLF.has_coerce_map_from(F) or + F.has_coerce_map_from(CLF)): + # The implementation of SR.random_element() currently just + # returns a random integer coerced into SR, so there is no + # benefit to allowing SR here when QQ is available. + raise ValueError("base ring of parent must be a subfield " + "of the complex numbers") + + I = identity_matrix(F,n) + A = random_matrix(F,n) + S = A - A.conjugate_transpose() + U = (S-I).inverse()*(S+I) + D = identity_matrix(F,n) + for i in range(n): + if ZZ.random_element(2).is_zero(): + D[i,i] *= F(-1) + return D*U + + @matrix_method def random_diagonalizable_matrix(parent, eigenvalues=None, dimensions=None): """ From 49ac0cdc3cba96f83ea3f04c06b8a337ad750685 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Mon, 7 Jul 2025 12:35:25 -0400 Subject: [PATCH 4/4] src/sage/matrix/special.py: optimize random_unitary_matrix() Use the lower-level methods Matrix.set_row_to_multiple_of_row() and random.random() to implement random_unitary_matrix(). This should bring with it a small speed increase, and actually makes the code simpler. Thanks to John Cremona for the suggestion. --- src/sage/matrix/special.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/sage/matrix/special.py b/src/sage/matrix/special.py index d6229b18928..76705b3d970 100644 --- a/src/sage/matrix/special.py +++ b/src/sage/matrix/special.py @@ -3215,11 +3215,16 @@ def random_unitary_matrix(parent): A = random_matrix(F,n) S = A - A.conjugate_transpose() U = (S-I).inverse()*(S+I) - D = identity_matrix(F,n) + + # Scale the rows of U by plus/minus one with equal probability. + # This generates the equivalence class of U according to the + # Liebeck/Osborne paper. + from random import random for i in range(n): - if ZZ.random_element(2).is_zero(): - D[i,i] *= F(-1) - return D*U + if random() < 0.5: + U.set_row_to_multiple_of_row(i, i, -1) + + return U @matrix_method