Ticket #10876: trac_10876-elementary-matrices-v4.patch

File trac_10876-elementary-matrices-v4.patch, 16.0 KB (added by rbeezer, 10 years ago)
  • sage/matrix/all.py

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1299281755 28800
    # Node ID 5ca241a9df64be62b44138e3ec13b61cf476ae79
    # Parent  6a679959b54b7975af1841df4dd77d28b0faef28
    10876: constructor for elementary matrices
    
    diff -r 6a679959b54b -r 5ca241a9df64 sage/matrix/all.py
    a b  
    11from matrix_space import MatrixSpace, is_MatrixSpace
    2 from constructor import matrix, Matrix, column_matrix, random_matrix, diagonal_matrix, identity_matrix, block_matrix, block_diagonal_matrix, jordan_block, zero_matrix, ones_matrix
     2from constructor import matrix, Matrix, column_matrix, random_matrix, diagonal_matrix, identity_matrix, block_matrix, block_diagonal_matrix, jordan_block, zero_matrix, ones_matrix, elementary_matrix
    33from matrix import is_Matrix
    44from berlekamp_massey import berlekamp_massey
    55
  • sage/matrix/constructor.py

    diff -r 6a679959b54b -r 5ca241a9df64 sage/matrix/constructor.py
    a b  
    16411641    one = ring(1)
    16421642    return matrix_space.MatrixSpace(ring, nrows, ncols, sparse).matrix([one]*nents)
    16431643
     1644def elementary_matrix(arg0, arg1=None, **kwds):
     1645    r"""
     1646    Creates a square matrix that corresponds to a row operation or a column operation.
     1647
     1648    FORMATS:
     1649
     1650    In each case, ``R`` is the base ring, and is optional. ``n`` is the size
     1651    of the square matrix created.  Any call may include the ``sparse`` keyword
     1652    to determine the representation used.  The default is ``False`` which
     1653    leads to a dense representation.  We describe the matrices by their
     1654    associated row operation, see the output description for more.
     1655
     1656    - ``elementary_matrix(R, n, row1=i, row2=j)``
     1657
     1658      The matrix which swaps rows ``i`` and ``j``.
     1659
     1660    - ``elementary_matrix(R, n, row1=i, scale=s)``
     1661
     1662      The matrix which multiplies row ``i`` by ``s``.
     1663
     1664    - ``elementary_matrix(R, n, row1=i, row2=j, scale=s)``
     1665
     1666      The matrix which multiplies row ``j`` by ``s``
     1667      and adds it to row ``i``.
     1668
     1669    Elementary matrices representing column operations are created
     1670    in an entirely analogous way, replacing ``row1`` by ``col1`` and
     1671    replacing ``row2`` by ``col2``.
     1672
     1673    Specifying the ring for entries of the matrix is optional.  If it
     1674    is not given, and a scale parameter is provided, then a ring containing
     1675    the value of ``scale`` will be used.  Otherwise, the ring defaults
     1676    to the integers.
     1677
     1678    OUTPUT:
     1679
     1680    An elementary matrix is a square matrix that is very close to being
     1681    an identity matrix.  If ``E`` is an elementary matrix and ``A`` is any
     1682    matrix with the same number of rows, then ``E*A`` is the result of
     1683    applying a row operation to ``A``.  This is how the three types
     1684    created by this function are described.  Similarly, an elementary matrix
     1685    can be associated with a column operation, so if ``E`` has the same number
     1686    of columns as ``A`` then ``A*E`` is the result of performing a column
     1687    operation on ``A``.
     1688
     1689    An elementary matrix representing a row operation is created if ``row1``
     1690    is specified, while an elementary matrix representing a column operation
     1691    is created if ``col1`` is specified.
     1692
     1693    EXAMPLES:
     1694
     1695    Over the integers, creating row operations. Recall that row
     1696    and column numbering begins at zero.  ::
     1697
     1698        sage: A = matrix(ZZ, 4, 10, range(40)); A
     1699        [ 0  1  2  3  4  5  6  7  8  9]
     1700        [10 11 12 13 14 15 16 17 18 19]
     1701        [20 21 22 23 24 25 26 27 28 29]
     1702        [30 31 32 33 34 35 36 37 38 39]
     1703
     1704        sage: E = elementary_matrix(4, row1=1, row2=3); E
     1705        [1 0 0 0]
     1706        [0 0 0 1]
     1707        [0 0 1 0]
     1708        [0 1 0 0]
     1709        sage: E*A
     1710        [ 0  1  2  3  4  5  6  7  8  9]
     1711        [30 31 32 33 34 35 36 37 38 39]
     1712        [20 21 22 23 24 25 26 27 28 29]
     1713        [10 11 12 13 14 15 16 17 18 19]
     1714
     1715        sage: E = elementary_matrix(4, row1=2, scale=10); E
     1716        [ 1  0  0  0]
     1717        [ 0  1  0  0]
     1718        [ 0  0 10  0]
     1719        [ 0  0  0  1]
     1720        sage: E*A
     1721        [  0   1   2   3   4   5   6   7   8   9]
     1722        [ 10  11  12  13  14  15  16  17  18  19]
     1723        [200 210 220 230 240 250 260 270 280 290]
     1724        [ 30  31  32  33  34  35  36  37  38  39]
     1725
     1726        sage: E = elementary_matrix(4, row1=2, row2=1, scale=10); E
     1727        [ 1  0  0  0]
     1728        [ 0  1  0  0]
     1729        [ 0 10  1  0]
     1730        [ 0  0  0  1]
     1731        sage: E*A
     1732        [  0   1   2   3   4   5   6   7   8   9]
     1733        [ 10  11  12  13  14  15  16  17  18  19]
     1734        [120 131 142 153 164 175 186 197 208 219]
     1735        [ 30  31  32  33  34  35  36  37  38  39]
     1736
     1737    Over the rationals, now as column operations. Recall that row
     1738    and column numbering begins at zero.  Checks now have the
     1739    elementary matrix on the right.  ::
     1740
     1741        sage: A = matrix(QQ, 5, 4, range(20)); A
     1742        [ 0  1  2  3]
     1743        [ 4  5  6  7]
     1744        [ 8  9 10 11]
     1745        [12 13 14 15]
     1746        [16 17 18 19]
     1747
     1748        sage: E = elementary_matrix(QQ, 4, col1=1, col2=3); E
     1749        [1 0 0 0]
     1750        [0 0 0 1]
     1751        [0 0 1 0]
     1752        [0 1 0 0]
     1753        sage: A*E
     1754        [ 0  3  2  1]
     1755        [ 4  7  6  5]
     1756        [ 8 11 10  9]
     1757        [12 15 14 13]
     1758        [16 19 18 17]
     1759
     1760        sage: E = elementary_matrix(QQ, 4, col1=2, scale=1/2); E
     1761        [  1   0   0   0]
     1762        [  0   1   0   0]
     1763        [  0   0 1/2   0]
     1764        [  0   0   0   1]
     1765        sage: A*E
     1766        [ 0  1  1  3]
     1767        [ 4  5  3  7]
     1768        [ 8  9  5 11]
     1769        [12 13  7 15]
     1770        [16 17  9 19]
     1771
     1772        sage: E = elementary_matrix(QQ, 4, col1=2, col2=1, scale=10); E
     1773        [ 1  0  0  0]
     1774        [ 0  1 10  0]
     1775        [ 0  0  1  0]
     1776        [ 0  0  0  1]
     1777        sage: A*E
     1778        [  0   1  12   3]
     1779        [  4   5  56   7]
     1780        [  8   9 100  11]
     1781        [ 12  13 144  15]
     1782        [ 16  17 188  19]
     1783
     1784    An elementary matrix is always nonsingular.  Then repeated row
     1785    operations can be represented by products of elementary matrices,
     1786    and this product is again nonsingular.  If row operations are to
     1787    preserve fundamental properties of a matrix (like rank), we do not
     1788    allow scaling a row by zero.  Similarly, the corresponding elementary
     1789    matrix is not constructed.  Also, we do not allow adding a multiple
     1790    of a row to itself, since this could also lead to a new zero row.  ::
     1791
     1792        sage: A = matrix(QQ, 4, 10, range(40)); A
     1793        [ 0  1  2  3  4  5  6  7  8  9]
     1794        [10 11 12 13 14 15 16 17 18 19]
     1795        [20 21 22 23 24 25 26 27 28 29]
     1796        [30 31 32 33 34 35 36 37 38 39]
     1797
     1798        sage: E1 = elementary_matrix(QQ, 4, row1=0, row2=1)
     1799        sage: E2 = elementary_matrix(QQ, 4, row1=3, row2=0, scale=100)
     1800        sage: E = E2*E1
     1801        sage: E.is_singular()
     1802        False
     1803        sage: E*A
     1804        [  10   11   12   13   14   15   16   17   18   19]
     1805        [   0    1    2    3    4    5    6    7    8    9]
     1806        [  20   21   22   23   24   25   26   27   28   29]
     1807        [1030 1131 1232 1333 1434 1535 1636 1737 1838 1939]
     1808
     1809        sage: E3 = elementary_matrix(QQ, 4, row1=3, scale=0)
     1810        Traceback (most recent call last):
     1811        ...
     1812        ValueError: scale parameter of row of elementary matrix must be non-zero
     1813
     1814        sage: E4 = elementary_matrix(QQ, 4, row1=3, row2=3, scale=12)
     1815        Traceback (most recent call last):
     1816        ...
     1817        ValueError: cannot add a multiple of a row to itself
     1818
     1819    If the ring is not specified, and a scale parameter is given, the
     1820    base ring for the matrix is chosen to contain the scale parameter.
     1821    Otherwise, if no ring is given, the default is the integers. ::
     1822
     1823        sage: E = elementary_matrix(4, row1=1, row2=3)
     1824        sage: E.parent()
     1825        Full MatrixSpace of 4 by 4 dense matrices over Integer Ring
     1826
     1827        sage: E = elementary_matrix(4, row1=1, scale=4/3)
     1828        sage: E.parent()
     1829        Full MatrixSpace of 4 by 4 dense matrices over Rational Field
     1830
     1831        sage: E = elementary_matrix(4, row1=1, scale=I)
     1832        sage: E.parent()
     1833        Full MatrixSpace of 4 by 4 dense matrices over Symbolic Ring
     1834
     1835        sage: E = elementary_matrix(4, row1=1, scale=CDF(I))
     1836        sage: E.parent()
     1837        Full MatrixSpace of 4 by 4 dense matrices over Complex Double Field
     1838
     1839        sage: E = elementary_matrix(4, row1=1, scale=QQbar(I))
     1840        sage: E.parent()
     1841        Full MatrixSpace of 4 by 4 dense matrices over Algebraic Field
     1842
     1843    Returned matrices have a dense implementation by default,
     1844    but a sparse implementation may be requested.  ::
     1845
     1846        sage: E = elementary_matrix(4, row1=0, row2=1)
     1847        sage: E.is_dense()
     1848        True
     1849
     1850        sage: E = elementary_matrix(4, row1=0, row2=1, sparse=True)
     1851        sage: E.is_sparse()
     1852        True
     1853
     1854    And the ridiculously small cases.  The zero-row matrix cannot be built
     1855    since then there are no rows to manipulate. ::
     1856
     1857        sage: elementary_matrix(QQ, 1, row1=0, row2=0)
     1858        [1]
     1859        sage: elementary_matrix(QQ, 0, row1=0, row2=0)
     1860        Traceback (most recent call last):
     1861        ...
     1862        ValueError: size of elementary matrix must be 1 or greater, not 0
     1863
     1864    TESTS::
     1865
     1866        sage: E = elementary_matrix('junk', 5, row1=3, row2=1, scale=12)
     1867        Traceback (most recent call last):
     1868        ...
     1869        TypeError: optional first parameter must be a ring, not junk
     1870
     1871        sage: E = elementary_matrix(5, row1=3, scale='junk')
     1872        Traceback (most recent call last):
     1873        ...
     1874        TypeError: scale must be an element of some ring, not junk
     1875
     1876        sage: E = elementary_matrix(ZZ, 5, row1=3, col2=3, scale=12)
     1877        Traceback (most recent call last):
     1878        ...
     1879        ValueError: received an unexpected keyword: col2=3
     1880
     1881        sage: E = elementary_matrix(QQ, row1=3, scale=12)
     1882        Traceback (most recent call last):
     1883        ...
     1884        ValueError: size of elementary matrix must be given
     1885
     1886        sage: E = elementary_matrix(ZZ, 4/3, row1=3, row2=1, scale=12)
     1887        Traceback (most recent call last):
     1888        ...
     1889        TypeError: size of elementary matrix must be an integer, not 4/3
     1890
     1891        sage: E = elementary_matrix(ZZ, -3, row1=3, row2=1, scale=12)
     1892        Traceback (most recent call last):
     1893        ...
     1894        ValueError: size of elementary matrix must be 1 or greater, not -3
     1895
     1896        sage: E = elementary_matrix(ZZ, 5, row2=1, scale=12)
     1897        Traceback (most recent call last):
     1898        ...
     1899        ValueError: row1 or col1 must be specified
     1900
     1901        sage: E = elementary_matrix(ZZ, 5, row1=3, col1=3, scale=12)
     1902        Traceback (most recent call last):
     1903        ...
     1904        ValueError: cannot specify both row1 and col1
     1905
     1906        sage: E = elementary_matrix(ZZ, 5, row1=4/3, row2=1, scale=12)
     1907        Traceback (most recent call last):
     1908        ...
     1909        TypeError: row of elementary matrix must be an integer, not 4/3
     1910
     1911        sage: E = elementary_matrix(ZZ, 5, col1=5, col2=1, scale=12)
     1912        Traceback (most recent call last):
     1913        ...
     1914        ValueError: column of elementary matrix must be positive and smaller than 5, not 5
     1915
     1916        sage: E = elementary_matrix(ZZ, 5, col1=3, col2=4/3, scale=12)
     1917        Traceback (most recent call last):
     1918        ...
     1919        TypeError: column of elementary matrix must be an integer, not 4/3
     1920
     1921        sage: E = elementary_matrix(ZZ, 5, row1=3, row2=-1, scale=12)
     1922        Traceback (most recent call last):
     1923        ...
     1924        ValueError: row of elementary matrix must be positive and smaller than 5, not -1
     1925
     1926        sage: E = elementary_matrix(ZZ, 5, row1=3, row2=1, scale=4/3)
     1927        Traceback (most recent call last):
     1928        ...
     1929        TypeError: scale parameter of elementary matrix must an element of Integer Ring, not 4/3
     1930
     1931        sage: E = elementary_matrix(ZZ, 5, row1=3)
     1932        Traceback (most recent call last):
     1933        ...
     1934        ValueError: insufficient parameters provided to construct elementary matrix
     1935
     1936        sage: E = elementary_matrix(ZZ, 5, row1=3, row2=3, scale=12)
     1937        Traceback (most recent call last):
     1938        ...
     1939        ValueError: cannot add a multiple of a row to itself
     1940
     1941        sage: E = elementary_matrix(ZZ, 5, col1=3, scale=0)
     1942        Traceback (most recent call last):
     1943        ...
     1944        ValueError: scale parameter of column of elementary matrix must be non-zero
     1945
     1946    AUTHOR:
     1947
     1948    - Rob Beezer (2011-03-04)
     1949    """
     1950    import sage.structure.element
     1951    # determine ring and matrix size
     1952    if not arg1 is None and not rings.is_Ring(arg0):
     1953        raise TypeError('optional first parameter must be a ring, not {0}'.format(arg0))
     1954    scale = kwds.pop('scale', None)
     1955    if rings.is_Ring(arg0):
     1956        R = arg0
     1957        arg0 = arg1
     1958    elif scale is not None:
     1959        if not sage.structure.element.is_RingElement(scale):
     1960            raise TypeError('scale must be an element of some ring, not {0}'.format(scale))
     1961        R = scale.parent()
     1962    else:
     1963        R = rings.ZZ
     1964    if arg0 is None:
     1965        raise ValueError('size of elementary matrix must be given')
     1966    try:
     1967        n = rings.Integer(arg0)
     1968    except:
     1969        raise TypeError('size of elementary matrix must be an integer, not {0}'.format(arg0))
     1970    if n <= 0:
     1971        raise ValueError('size of elementary matrix must be 1 or greater, not {0}'.format(n))
     1972    # row operations or column operations?
     1973    # column operation matrix will be transpose of a row operation matrix
     1974    row1 = kwds.pop('row1', None)
     1975    col1 = kwds.pop('col1', None)
     1976    if row1 is None and col1 is None:
     1977        raise ValueError('row1 or col1 must be specified')
     1978    if not row1 is None and not col1 is None:
     1979        raise ValueError('cannot specify both row1 and col1')
     1980    rowop = not row1 is None
     1981    if rowop:
     1982        opstring = "row"
     1983        row2 = kwds.pop('row2', None)
     1984    else:
     1985        opstring = "column"
     1986        row1 = col1
     1987        row2 = kwds.pop('col2', None)
     1988    sparse = kwds.pop('sparse', False)
     1989    if kwds:
     1990        extra = kwds.popitem()
     1991        raise ValueError('received an unexpected keyword: {0}={1}'.format(extra[0], extra[1]))
     1992
     1993    # analyze parameters to determine matrix type
     1994    try:
     1995        row1 = rings.Integer(row1)
     1996    except:
     1997        raise TypeError('{0} of elementary matrix must be an integer, not {1}'.format(opstring, row1))
     1998    if row1 < 0 or row1 >= n :
     1999        raise ValueError('{0} of elementary matrix must be positive and smaller than {1}, not {2}'.format(opstring, n, row1))
     2000    if not row2 is None:
     2001        try:
     2002            row2 = rings.Integer(row2)
     2003        except:
     2004            raise TypeError('{0} of elementary matrix must be an integer, not {1}'.format(opstring, row2))
     2005        if row2 < 0 or row2 >= n :
     2006            raise ValueError('{0} of elementary matrix must be positive and smaller than {1}, not {2}'.format(opstring, n, row2))
     2007    if not scale is None:
     2008        try:
     2009            scale = R(scale)
     2010        except:
     2011            raise TypeError('scale parameter of elementary matrix must an element of {0}, not {1}'.format(R, scale))
     2012
     2013    # determine type of matrix and adjust an identity matrix
     2014    # return row operation matrix or the transpose as a column operation matrix
     2015    elem = identity_matrix(R, n, sparse=sparse)
     2016    if row2 is None and scale is None:
     2017        raise ValueError('insufficient parameters provided to construct elementary matrix')
     2018    elif not row2 is None and not scale is None:
     2019        if row1 == row2:
     2020            raise ValueError('cannot add a multiple of a {0} to itself'.format(opstring))
     2021        elem[row1, row2] = scale
     2022    elif not row2 is None and scale is None:
     2023        elem[row1, row1] = 0
     2024        elem[row2, row2] = 0
     2025        elem[row1, row2] = 1
     2026        elem[row2, row1] = 1
     2027    elif row2 is None and not scale is None:
     2028        if scale == 0:
     2029            raise ValueError('scale parameter of {0} of elementary matrix must be non-zero'.format(opstring))
     2030        elem[row1, row1] = scale
     2031    if rowop:
     2032        return elem
     2033    else:
     2034        return elem.transpose()
     2035
    16442036def _determine_block_matrix_grid(sub_matrices):
    16452037    r"""
    16462038    For internal use. This tries to determine the dimensions