# HG changeset patch
# User Billy Wonderly <bwonderly@ups.edu>
# Date 1282951283 25200
# Node ID df08ce5b6045f88976e9b9f1d9d35dcfe8529e20
# Parent  7f94d055da42997b9edb78b851c7aee9f341d282
Trac 9754: random subspaces and unimodular  matrices

diff -r 7f94d055da42 -r df08ce5b6045 sage/matrix/constructor.py
--- a/sage/matrix/constructor.py	Fri Aug 27 17:51:38 2010 -0700
+++ b/sage/matrix/constructor.py	Fri Aug 27 16:21:23 2010 -0700
@@ -808,7 +808,13 @@
 
        -  ``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.
 
        - ``diagonalizable`` - creates a diagonalizable matrix whose
          eigenvectors, if computed by hand, will have only integer
@@ -1008,7 +1014,7 @@
     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()
@@ -1061,6 +1067,58 @@
 
         from sage.matrix.constructor import random_diagonalizable_matrix
 
+    Random matrices with predictable subspaces.  The ``algorithm='subspaces'``
+    keyword, along with an optional rank (``rank``) 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 ``rank``, is not set, the 
+    rank of the matrix will be generated randomly. ::
+
+        sage: B=random_matrix(QQ, 5, 6, algorithm='subspaces', rank=3); B #random
+        sage: B_expanded=B.augment(identity_matrix(5)).rref()
+        sage: (B.nrows(), B.ncols())
+        (5, 6)
+        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``::
@@ -1096,6 +1154,10 @@
         return random_echelonizable_matrix(parent, *args, **kwds)
     elif algorithm == 'diagonalizable':
         return random_diagonalizable_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)
 
@@ -1590,10 +1652,10 @@
 
     Matrices must have the number of pivot columns be greater than zero. ::
 
-        sage: random_matrix(QQ, 5, 4, algorithm='echelon_form', num_pivots=0)
+        sage: random_matrix(QQ, 5, 4, algorithm='echelon_form', num_pivots=-1)
         Traceback (most recent call last):
         ...
-        ValueError: the number of pivots must be greater than zero.
+        ValueError: the number of pivots must be zero or greater.
 
     AUTHOR:
 
@@ -1607,8 +1669,8 @@
         num_pivots=ZZ(num_pivots)
     except TypeError:
         raise TypeError("the number of pivots must be an integer.")
-    if num_pivots<=0:
-        raise ValueError("the number of pivots must be greater than zero.")
+    if num_pivots<0:
+        raise ValueError("the number of pivots must be zero or greater.")
     ring = parent.base_ring()
     if not ring.is_exact():
         raise TypeError("the base ring must be exact.")
@@ -1620,7 +1682,7 @@
         one=ring.one()
         # Create a matrix of the desired size to be modified and then returned.
         return_matrix=copy(parent.zero_matrix())
-        pivots=[0] #Force first column to be a pivot.
+        pivots=[0] #Force first column to be a pivot. No harm if no pivots at all.
         # Probability distribution for the placement of leading one's.
         pivot_generator=pd.RealDistribution("beta",[1.6,4.3]) 
         while len(pivots)<num_pivots: 
@@ -1746,7 +1808,7 @@
         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]
@@ -1764,10 +1826,10 @@
     Matrices must have a rank greater than zero.
     (For zero rank, just create a zero matrix.) ::
 
-        sage: random_matrix(QQ, 3, 4, algorithm='echelonizable', rank=0)
+        sage: random_matrix(QQ, 3, 4, algorithm='echelonizable', rank=-1)
         Traceback (most recent call last):
         ...
-        ValueError: matrices must have rank greater than zero.
+        ValueError: matrices must have rank zero or greater.
 
     The base ring must be exact. ::
 
@@ -1785,8 +1847,8 @@
 
     ring = parent.base_ring()
     rows = parent.nrows()
-    if rank<=0:
-        raise ValueError("matrices must have rank greater than zero.")
+    if rank<0:
+        raise ValueError("matrices must have rank zero or greater.")
     matrix = random_rref_matrix(parent, rank)
     # Entries of matrices over the ZZ or QQ can get large, entry size is regulated by finding the largest
     # entry of the resultant matrix after addition of scalar multiple of a row.
@@ -1804,6 +1866,7 @@
                 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.
@@ -1849,6 +1912,277 @@
             matrix.add_multiple_of_row(0,randint(1,rows-1),ring.random_element())
     return matrix
 
+def random_subspaces_matrix(parent, rank=None):
+    """
+    Create a matrix of the designated size and rank 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.
+
+    - ``rank`` - The desired rank 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 rank 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, rank=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.rank()
+        3
+        sage: B.nullity()
+        3
+        sage: (B.nrows(), B.ncols())
+        (6, 8)
+        sage: all([x in ZZ for x in B.list()])
+        True
+        sage: B_expanded=B.augment(identity_matrix(6)).rref()
+        sage: all([x in ZZ for x in B_expanded.list()])
+        True
+        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', rank=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.nrows(), A.ncols())
+        (5, 7)
+        sage: all([x in ZZ for x in A.list()])
+        True
+        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: all([x in ZZ for x in A_expanded.list()])
+        True
+        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', rank=21)
+        Traceback (most recent call last):
+        ...
+        ValueError: rank 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.ncols()
+
+    if rank>rows or rank>columns:
+        raise ValueError("rank cannot exceed the number of rows or columns.")
+    # If rank is not designated, generate using probability distribution skewing to smaller numbers, always at least 1.
+    if rank==None:
+        left_nullity_generator=pd.RealDistribution("beta",[1.4,5.5])
+        nullity=int(left_nullity_generator.get_random_element()*(rows-1)+1)
+        rank=rows-nullity
+    else:
+        rank=rank
+    nullity=rows-rank
+    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)
+    # by multiplying the B matrix by J.inverse() we hide the B matrix
+    # of the solution using row operations required to change the solution
+    # K matrix to the identity matrix.
+    return J.inverse()*B
+
+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)
+
+
 def random_diagonalizable_matrix(parent,eigenvalues=None,dimensions=None):
     """
     Create a random matrix that diagonalizes nicely.