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

File trac_10876-elementary-matrices-v3.patch, 15.7 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 55747d5dc011fd72be72940cc10d5ef5924d3d9a
    # Parent  6a679959b54b7975af1841df4dd77d28b0faef28
    10876: constructor for elementary matrices
    
    diff -r 6a679959b54b -r 55747d5dc011 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 55747d5dc011 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 fails
     1855    since 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: row of elementary matrix must be positive and smaller than 0, not 0
     1863
     1864    TESTS::
     1865
     1866        sage: E = elementary_matrix(5, row1=3, scale='junk')
     1867        Traceback (most recent call last):
     1868        ...
     1869        TypeError: scale must be an element of some ring, not junk
     1870
     1871        sage: E = elementary_matrix(ZZ, 5, row1=3, col2=3, scale=12)
     1872        Traceback (most recent call last):
     1873        ...
     1874        ValueError: received an unexpected keyword: col2=3
     1875
     1876        sage: E = elementary_matrix(QQ, row1=3, scale=12)
     1877        Traceback (most recent call last):
     1878        ...
     1879        ValueError: size of elementary matrix must be given
     1880
     1881        sage: E = elementary_matrix(ZZ, 4/3, row1=3, row2=1, scale=12)
     1882        Traceback (most recent call last):
     1883        ...
     1884        TypeError: size of elementary matrix must be an integer, not 4/3
     1885
     1886        sage: E = elementary_matrix(ZZ, -3, row1=3, row2=1, scale=12)
     1887        Traceback (most recent call last):
     1888        ...
     1889        ValueError: size of elementary matrix must be positive, not -3
     1890
     1891        sage: E = elementary_matrix(ZZ, 5, row2=1, scale=12)
     1892        Traceback (most recent call last):
     1893        ...
     1894        ValueError: row1 or col1 must be specified
     1895
     1896        sage: E = elementary_matrix(ZZ, 5, row1=3, col1=3, scale=12)
     1897        Traceback (most recent call last):
     1898        ...
     1899        ValueError: cannot specify both row1 and col1
     1900
     1901        sage: E = elementary_matrix(ZZ, 5, row1=4/3, row2=1, scale=12)
     1902        Traceback (most recent call last):
     1903        ...
     1904        TypeError: row of elementary matrix must be an integer, not 4/3
     1905
     1906        sage: E = elementary_matrix(ZZ, 5, col1=5, col2=1, scale=12)
     1907        Traceback (most recent call last):
     1908        ...
     1909        ValueError: column of elementary matrix must be positive and smaller than 5, not 5
     1910
     1911        sage: E = elementary_matrix(ZZ, 5, col1=3, col2=4/3, scale=12)
     1912        Traceback (most recent call last):
     1913        ...
     1914        TypeError: column of elementary matrix must be an integer, not 4/3
     1915
     1916        sage: E = elementary_matrix(ZZ, 5, row1=3, row2=-1, scale=12)
     1917        Traceback (most recent call last):
     1918        ...
     1919        ValueError: row of elementary matrix must be positive and smaller than 5, not -1
     1920
     1921        sage: E = elementary_matrix(ZZ, 5, row1=3, row2=1, scale=4/3)
     1922        Traceback (most recent call last):
     1923        ...
     1924        TypeError: scale parameter of elementary matrix must an element of Integer Ring, not 4/3
     1925
     1926        sage: E = elementary_matrix(ZZ, 5, row1=3)
     1927        Traceback (most recent call last):
     1928        ...
     1929        ValueError: insufficient parameters provided to construct elementary matrix
     1930
     1931        sage: E = elementary_matrix(ZZ, 5, row1=3, row2=3, scale=12)
     1932        Traceback (most recent call last):
     1933        ...
     1934        ValueError: cannot add a multiple of a row to itself
     1935
     1936        sage: E = elementary_matrix(ZZ, 5, col1=3, scale=0)
     1937        Traceback (most recent call last):
     1938        ...
     1939        ValueError: scale parameter of column of elementary matrix must be non-zero
     1940
     1941    AUTHOR:
     1942
     1943    - Rob Beezer (2011-03-04)
     1944    """
     1945    import sage.structure.element
     1946    # determine ring and matrix size
     1947    scale = kwds.pop('scale', None)
     1948    if rings.is_Ring(arg0):
     1949        R = arg0
     1950        arg0 = arg1
     1951    elif scale is not None:
     1952        if not sage.structure.element.is_RingElement(scale):
     1953            raise TypeError('scale must be an element of some ring, not {0}'.format(scale))
     1954        R = scale.parent()
     1955    else:
     1956        R = rings.ZZ
     1957    if arg0 is None:
     1958        raise ValueError('size of elementary matrix must be given')
     1959    try:
     1960        n = rings.Integer(arg0)
     1961    except:
     1962        raise TypeError('size of elementary matrix must be an integer, not {0}'.format(arg0))
     1963    if n < 0:
     1964        raise ValueError('size of elementary matrix must be positive, not {0}'.format(n))
     1965    # row operations or column operations?
     1966    # column operation matrix will be transpose of a row operation matrix
     1967    row1 = kwds.pop('row1', None)
     1968    col1 = kwds.pop('col1', None)
     1969    if row1 is None and col1 is None:
     1970        raise ValueError('row1 or col1 must be specified')
     1971    if not row1 is None and not col1 is None:
     1972        raise ValueError('cannot specify both row1 and col1')
     1973    rowop = not row1 is None
     1974    if rowop:
     1975        opstring = "row"
     1976        row2 = kwds.pop('row2', None)
     1977    else:
     1978        opstring = "column"
     1979        row1 = col1
     1980        row2 = kwds.pop('col2', None)
     1981    sparse = kwds.pop('sparse', False)
     1982    if kwds:
     1983        extra = kwds.popitem()
     1984        raise ValueError('received an unexpected keyword: {0}={1}'.format(extra[0], extra[1]))
     1985
     1986    # analyze parameters to determine matrix type
     1987    try:
     1988        row1 = rings.Integer(row1)
     1989    except:
     1990        raise TypeError('{0} of elementary matrix must be an integer, not {1}'.format(opstring, row1))
     1991    if row1 < 0 or row1 >= n :
     1992        raise ValueError('{0} of elementary matrix must be positive and smaller than {1}, not {2}'.format(opstring, n, row1))
     1993    if not row2 is None:
     1994        try:
     1995            row2 = rings.Integer(row2)
     1996        except:
     1997            raise TypeError('{0} of elementary matrix must be an integer, not {1}'.format(opstring, row2))
     1998        if row2 < 0 or row2 >= n :
     1999            raise ValueError('{0} of elementary matrix must be positive and smaller than {1}, not {2}'.format(opstring, n, row2))
     2000    if not scale is None:
     2001        try:
     2002            scale = R(scale)
     2003        except:
     2004            raise TypeError('scale parameter of elementary matrix must an element of {0}, not {1}'.format(R, scale))
     2005
     2006    # determine type of matrix and adjust an identity matrix
     2007    # return row operation matrix or the transpose as a column operation matrix
     2008    elem = identity_matrix(R, n, sparse=sparse)
     2009    if row2 is None and scale is None:
     2010        raise ValueError('insufficient parameters provided to construct elementary matrix')
     2011    elif not row2 is None and not scale is None:
     2012        if row1 == row2:
     2013            raise ValueError('cannot add a multiple of a {0} to itself'.format(opstring))
     2014        elem[row1, row2] = scale
     2015    elif not row2 is None and scale is None:
     2016        elem[row1, row1] = 0
     2017        elem[row2, row2] = 0
     2018        elem[row1, row2] = 1
     2019        elem[row2, row1] = 1
     2020    elif row2 is None and not scale is None:
     2021        if scale == 0:
     2022            raise ValueError('scale parameter of {0} of elementary matrix must be non-zero'.format(opstring))
     2023        elem[row1, row1] = scale
     2024    if rowop:
     2025        return elem
     2026    else:
     2027        return elem.transpose()
     2028
    16442029def _determine_block_matrix_grid(sub_matrices):
    16452030    r"""
    16462031    For internal use. This tries to determine the dimensions