Ticket #9754: trac_9754-random-subspace-unimodular-matrix-v2.patch

sage/matrix/constructor.py

@@ -808 +808 @@
 
        -  ``echelon_form`` - creates a matrix in echelon form
 
-       -  ``echelonizable`` - creates a matrix that has a predictable echelon form
+       -  ``echelonizable`` - creates a matrix that has a predictable
+          echelon form
+
+       - ``subspaces`` - creates a matrix whose four subspaces, when
+         explored, have reasonably sized, integral valued, entries.
+
+       - ``unimodular`` - creates a matrix of determinant 1.
 
     -  ``*args, **kwds`` - arguments and keywords to describe additional properties.
        See more detailed documentation below.
@@ -1004 +1010 @@
     matrix also has integer entries.  Other exact rings may be also
     specified, but there is no notion of controlling the size.  Square matrices
     of full rank generated by this function always have determinant one, and
-    can be constructed with the upcoming ``unimodular`` keyword. ::
+    can be constructed with the ``unimodular`` keyword. ::
 
         sage: A=random_matrix(QQ, 4, 8, algorithm='echelonizable', rank=3, upper_bound=60); A # random
         sage: A.base_ring()
@@ -1026 +1032 @@
 
         from sage.matrix.constructor import random_echelonizable_matrix
 
+    Random matrices with predictable subspaces.  The ``algorithm='subspaces'``
+    keyword, along with an optional nullity (``nullity``) will return 
+    a matrix whose natural basis vectors for its four fundamental subspaces, if computed as
+    described in the documentation of the :func:`~sage.matrix.constructor.random_subspaces_matrix`
+    contain only integer entries.  If ``nullity``, is not set, the 
+    nullity of the matrix will be generated randomly. ::
+
+        sage: B=random_matrix(QQ, 5, 6, algorithm='subspaces', nullity=2); B #random
+        sage: B_expanded=B.augment(identity_matrix(5)).rref()
+        sage: all([x in ZZ for x in B_expanded.list()])
+        True
+        sage: C=B_expanded.submatrix(0,0,B.nrows()-B.nullity(),B.ncols())
+        sage: L=B_expanded.submatrix(B.nrows()-B.nullity(),B.ncols())
+        sage: B.right_kernel()==C.right_kernel()
+        True
+        sage: B.row_space()==C.row_space()
+        True
+        sage: B.column_space()==L.right_kernel()
+        True
+        sage: B.left_kernel()==L.row_space()
+        True
+
+    For more, see the documentation of the :func:`~sage.matrix.constructor.random_subspaces_matrix`
+    function.  In the notebook or at the Sage command-line, first execute the following to make
+    this further documentation available::
+
+        from sage.matrix.constructor import random_subspaces_matrix
+
+    Random unimodular matrices.  The ``algorithm='unimodular'``
+    keyword, along with an optional entry size control (``upper_bound``)
+    will return a matrix of determinant 1. When the base ring is ``ZZ``
+    or ``QQ`` the result has integer entries, whose magnitudes
+    can be limited by the value of ``upper_bound``. ::
+
+        sage: C=random_matrix(QQ, 5, algorithm='unimodular', upper_bound=70); C # random
+        sage: det(C)
+        1
+        sage: C.base_ring()
+        Rational Field
+        sage: (C.nrows(), C.ncols())
+        (5, 5)
+        sage: all([abs(x)<70 for x in C.list()])
+        True
+
+    For more, see the documentation of the :func:`~sage.matrix.constructor.random_unimodular_matrix`
+    function.  In the notebook or at the Sage command-line, first execute the following to make
+    this further documentation available::
+
+        from sage.matrix.constructor import random_unimodular_matrix
+
     TESTS:
 
     We return an error for a bogus value of ``algorithm``::
@@ -1059 +1115 @@
         return random_rref_matrix(parent, *args, **kwds)
     elif algorithm == 'echelonizable':
         return random_echelonizable_matrix(parent, *args, **kwds)
+    elif algorithm == 'subspaces':
+        return random_subspaces_matrix(parent, *args, **kwds)
+    elif algorithm == 'unimodular':
+        return random_unimodular_matrix(parent, *args, **kwds)
     else:
         raise ValueError('random matrix algorithm "%s" is not recognized' % algorithm)
 
@@ -1709 +1769 @@
         sage: B.rank()
         4
 
-    Square matrices over ZZ or QQ with full rank are unimodular. ::
+    Square matrices over ZZ or QQ with full rank are always unimodular. ::
 
         sage: E=random_matrix(QQ, 7, 7, algorithm='echelonizable', rank=7); E # random
         [  1  -1   5  12 -24 -41  47]
@@ -1767 +1827 @@
                 row_index=0
                 while row_index<rows:
                     # To each row in a pivot column add a scalar multiple of the pivot row.
+                    # for full rank, square matrices, using only this row operation preserves the determinant of 1.
                     if pivots!=row_index:
                     # To ensure a leading one is not removed by the addition of the pivot row by its
                     # additive inverse.
@@ -1811 +1872 @@
         if rows>1:
             matrix.add_multiple_of_row(0,randint(1,rows-1),ring.random_element())
     return matrix
+
+def random_subspaces_matrix(parent, nullity=None):
+    """
+    Create a matrix of the designated size and nullity whose right and
+    left null spaces, column space, and row space have desirable
+    properties that simplify the subspaces. 
+
+    INPUT:
+
+    - ``parent`` - A matrix space specifying the base ring, dimensions, and
+      representation (dense/sparse) for the result.  The base ring must be exact.
+
+    - ``nullity`` - The desired nullity of the return matrix. (Default: None)
+
+    OUTPUT:
+
+    A matrix whose natrual basis vectors for its four subspaces, when 
+    computed, have reasonably sized, integral valued, entries.
+
+    .. note::
+
+        It is easiest to use this function via a call to the
+        :func:`~sage.matrix.constructor.random_matrix`
+        function with the ``algorithm='subspaces'`` keyword.  We provide
+        one example accessing this function directly, while the remainder will
+        use this more general function.
+
+    EXAMPLES:
+
+    A 6x8 matrix with designated nullity of 3.  The four subspaces are
+    determined using one simple routine in which we augment the
+    original matrix with the equal row dimension identity matrix.  The
+    resulting matrix is then put in reduced row-echelon form and the
+    subspaces can then be determined by analyzing subdivisions of this
+    matrix. See the four subspaces routine in [BEEZER]_ for more. ::
+
+        sage: from sage.matrix.constructor import random_subspaces_matrix
+        sage: matrix_space = sage.matrix.matrix_space.MatrixSpace(QQ, 6, 8)
+        sage: B=random_subspaces_matrix(matrix_space, nullity=3); B # random
+        [ 113  339  -46  218 -243  641 -269 -306]
+        [ -33  -99   13  -63   69 -185   77   90]
+        [  35  105  -14   67  -74  197  -82  -95]
+        [ -18  -54    7  -34   37 -100   41   49]
+        [ -26  -78   10  -49   53 -144   59   71]
+        [  -8  -24    3  -15   16  -44   18   22]
+        sage: B.nullity()
+        3
+        sage: B_expanded=B.augment(identity_matrix(6)).rref()
+        sage: B_expanded # random
+        [  1   3   0   0   1   1   3  -2   0   0  -3   0  -9  16]
+        [  0   0   1   0   3  -2  -1  -3   0   0   2   0  11 -27]
+        [  0   0   0   1  -1   2  -3  -1   0   0   2   0   7 -14]
+        [  0   0   0   0   0   0   0   0   1   0  -5   0  -3   2]
+        [  0   0   0   0   0   0   0   0   0   1   1   0   1  -3]
+        [  0   0   0   0   0   0   0   0   0   0   0   1  -1   1]
+        sage: B_expanded.subdivide(B.nrows()-B.nullity(),B.ncols());B_expanded # random
+        [  1   3   0   0   1   1   3  -2|  0   0  -3   0  -9  16]
+        [  0   0   1   0   3  -2  -1  -3|  0   0   2   0  11 -27]
+        [  0   0   0   1  -1   2  -3  -1|  0   0   2   0   7 -14]
+        [-------------------------------+-----------------------]
+        [  0   0   0   0   0   0   0   0|  1   0  -5   0  -3   2]
+        [  0   0   0   0   0   0   0   0|  0   1   1   0   1  -3]
+        [  0   0   0   0   0   0   0   0|  0   0   0   1  -1   1]
+        sage: C=B_expanded.subdivision(0,0)
+        sage: C # random
+        [ 1  3  0  0  1  1  3 -2]
+        [ 0  0  1  0  3 -2 -1 -3]
+        [ 0  0  0  1 -1  2 -3 -1]
+        sage: L=B_expanded.subdivision(1,1)
+        sage: L # random
+        [ 1  0 -5  0 -3  2]
+        [ 0  1  1  0  1 -3]
+        [ 0  0  0  1 -1  1]
+        sage: B.right_kernel()==C.right_kernel()
+        True
+        sage: B.row_space()==C.row_space()
+        True
+        sage: B.column_space()==L.right_kernel()
+        True
+        sage: B.left_kernel()==L.row_space()
+        True
+
+    A matrix to show that the null space of the L matrix is the column space of the starting matrix. ::
+
+        sage: A=random_matrix(QQ, 5, 7, algorithm='subspaces', nullity=None); A # random
+        [-31  12  -9 -27  21   2 -15]
+        [105 -24   6 103 -30 -34  79]
+        [ 29  -9   5  26 -14  -5  17]
+        [233 -55  16 228 -71 -73 173]
+        [-42  10  -3 -41  13  13 -31]
+        sage: A.nullity() # random
+        1
+        sage: A_expanded=A.augment(identity_matrix(5)).rref()
+        sage: A_expanded # random
+        [  1   0   0   0   0   1   0   0   3   7  25 151]
+        [  0   1   0   0   1   2   1   0   5  21  84 493]
+        [  0   0   1   0  -1   2   0   0   2  13  53 308]
+        [  0   0   0   1   0  -1   1   0  -2  -3  -9 -57]
+        [  0   0   0   0   0   0   0   1  -3   1   1  -2]
+        sage: C=A_expanded.submatrix(0,0,A.nrows()-A.nullity(),A.ncols())
+        sage: L=A_expanded.submatrix(A.nrows()-A.nullity(),A.ncols())
+        sage: A.right_kernel()==C.right_kernel()
+        True
+        sage: A.row_space()==C.row_space()
+        True
+        sage: A.column_space()==L.right_kernel()
+        True
+        sage: A.left_kernel()==L.row_space()
+        True
+
+    TESTS:
+
+    The designated rank of the L matrix cannot be greater than the number of desired rows. ::
+
+        sage: random_matrix(QQ, 19, 20, algorithm='subspaces', nullity=21)
+        Traceback (most recent call last):
+        ...
+        ValueError: nullity cannot exceed the number of rows or columns.
+
+    REFERENCES:
+
+        .. [BEEZER] `A First Course in Linear Algebra <http://linear.ups.edu/>`_.
+           Robert A. Beezer, accessed 15 July 2010.
+
+    AUTHOR:
+
+    Billy Wonderly (2010-07)
+    """
+
+    import sage.gsl.probability_distribution as pd
+
+    ring = parent.base_ring()
+    rows = parent.nrows()
+    columns = parent.nrows()
+
+    if nullity>rows:
+        raise ValueError("nullity cannot exceed the number of rows or columns.")
+    # If nullity is not designated, generate using probability distribution skewing to smaller numbers, always at least 1.
+    if nullity==None:
+        nullity_generator=pd.RealDistribution("beta",[1.4,5.5])
+        nullity=int(nullity_generator.get_random_element()*(rows-1)+1)
+    else:
+        nullity=nullity
+    rank=rows-nullity
+    B=random_matrix(ring, rows, columns, algorithm='echelon_form', num_pivots=rank)
+    # Create a nonsingular matrix whose columns will be used to stack a matrix over the L matrix, forming a nonsingular matrix.
+    K_nonzero_columns=random_matrix(ring, rank, rank, algorithm='echelonizable', rank=rank) 
+    K=matrix(QQ,rank,rows)
+    L=random_matrix(ring, nullity, rows, algorithm='echelon_form', num_pivots=nullity)
+    for column in range(len(L.nonpivots())):
+        for entry in range(rank):
+            K[entry,L.nonpivots()[column]]=K_nonzero_columns[entry,column]
+    J=K.stack(L)
+    J_inverse=J.inverse()
+    N=B.augment(J)
+    M=J_inverse*N # Entries remain whole values as J is unimodular.
+    return_matrix=M.submatrix(0,0,M.nrows(),M.ncols()-M.nrows())
+    return return_matrix
+
+def random_unimodular_matrix(parent, upper_bound=None):
+    """
+    Generate a random unimodular (determinant 1) matrix of a desired size over a desired ring.
+
+    INPUT:
+
+    - ``parent`` - A matrix space specifying the base ring, dimensions
+      and representation (dense/sparse) for the result.  The base ring
+      must be exact.
+
+    - ``upper_bound`` - For large matrices over QQ or ZZ,
+      ``upper_bound`` is the largest value matrix entries can achieve.
+
+    OUTPUT:
+
+    An invertible matrix with the desired properties and determinant 1.
+
+    .. note::
+
+        It is easiest to use this function via a call to the
+        :func:`~sage.matrix.constructor.random_matrix`
+        function with the ``algorithm='unimodular'`` keyword.  We provide
+        one example accessing this function directly, while the remainder will
+        use this more general function.
+
+    EXAMPLES:
+
+    A matrix size 5 over QQ. ::
+
+        sage: from sage.matrix.constructor import random_unimodular_matrix
+        sage: matrix_space = sage.matrix.matrix_space.MatrixSpace(QQ, 5)
+        sage: A=random_unimodular_matrix(matrix_space); A # random
+        [  -8   31   85  148 -419]
+        [   2   -9  -25  -45  127]
+        [  -1   10   30   65 -176]
+        [  -3   12   33   58 -164]
+        [   5  -21  -59 -109  304]
+        sage: det(A)
+        1
+
+    A matrix size 6 with entries no larger than 50. ::
+
+        sage: B=random_matrix(ZZ, 7, algorithm='unimodular', upper_bound=50);B # random
+        [ -1   0   3   1  -2   2   9]
+        [  1   2  -5   0  14  19 -49]
+        [ -3  -2  12   5  -6  -4  24]
+        [  1   2  -9  -3   3   4  -7]
+        [ -2  -1   7   2  -8  -5  31]
+        [  2   2  -6  -3   8  16 -32]
+        [  1   2  -9  -2   5   6 -12]
+        sage: det(B)
+        1
+
+    A matrix over the number Field in `y` with defining polynomial `y^2-2y-2`. ::
+
+        sage: y = var('y')
+        sage: K=NumberField(y^2-2*y-2,'y')
+        sage: C=random_matrix(K, 3, algorithm='unimodular');C # random
+        [   2*y - 33 681*y - 787   31*y - 37]
+        [      y + 6 -155*y + 83    -7*y + 4]
+        [         -y   24*y + 51       y + 3]
+        sage: det(C)
+        1
+
+    TESTS:
+
+    Unimodular matrices are square. ::
+
+        sage: random_matrix(QQ, 5, 6, algorithm='unimodular')
+        Traceback (most recent call last):
+        ...
+        TypeError: a unimodular matrix must be square.
+
+    Only matrices over ZZ and QQ can have size control. ::
+
+        sage: F.<a>=GF(5^7)
+        sage: random_matrix(F, 5, algorithm='unimodular', upper_bound=20)
+        Traceback (most recent call last):
+        ...
+        TypeError: only matrices over ZZ or QQ can have size control.
+
+    AUTHOR:
+
+    Billy Wonderly (2010-07)
+    """
+
+    ring=parent.base_ring()
+    size=parent.nrows()
+    if parent.nrows()!=parent.ncols():
+        raise TypeError("a unimodular matrix must be square.")
+    if upper_bound!=None and (ring!=ZZ and ring!=QQ):
+        raise TypeError("only matrices over ZZ or QQ can have size control.")
+    if upper_bound==None:
+        # random_echelonizable_matrix() always returns a determinant one matrix if given full rank.
+        return random_matrix(ring, size, algorithm='echelonizable', rank=size)
+    elif upper_bound!=None and (ring==ZZ or ring==QQ):
+        return random_matrix(ring, size,algorithm='echelonizable',rank=size, upper_bound=upper_bound)
+