Ticket #11608: trac_11608-eigenvalues-symmetric-multiplicity-rebase.patch

File trac_11608-eigenvalues-symmetric-multiplicity-rebase.patch, 11.4 KB (added by rbeezer, 9 years ago)

Rebased on 5.0-beta10

  • sage/matrix/matrix_double_dense.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1336024761 25200
    # Node ID 498d7d78a4d02f5bf82ec0bfa9233305bb3b8dae
    # Parent  a6beb86b32aaf7e43db585455b87827fb29afd74
    11608: RDF/CDF eigenvalues, symmetric matrices, multiplicities
    
    diff --git a/sage/matrix/matrix_double_dense.pyx b/sage/matrix/matrix_double_dense.pyx
    a b  
    14481448
    14491449    right_eigenspaces = eigenspaces_right
    14501450
    1451     def eigenvalues(self):
     1451    def eigenvalues(self, algorithm='default', tol=None):
    14521452        r"""
    1453         Returns a list of the eigenvalues (with multiplicity)
    1454         of this matrix.  The returned eigenvalues are elements of CDF.
     1453        Returns a list of eigenvalues.
     1454
     1455
     1456        INPUT:
     1457
     1458        - ``self`` - a square matrix
     1459
     1460        - ``algorithm`` - default: ``'default'``
     1461
     1462          - ``'default'`` - applicable to any matrix
     1463            with double-precision floating point entries.
     1464            Uses the :meth:`~scipy.linalg.eigvals` method from SciPy.
     1465
     1466          - ``'symmetric'`` - checks that any complex matrix
     1467            (i.e. with entries from :class:`~sage.rings.complex_double.CDF`)
     1468            can be coerced into a real matrix (i.e. with entries from
     1469            :class:`~sage.rings.real_double.RDF`), then applies the
     1470            algorithm for Hermitian matrices.  This algorithm can
     1471            be significantly faster than the ``'default'`` algorithm.
     1472
     1473          - ``'hermitian'`` - uses the :meth:`~scipy.linalg.eigh` method
     1474            from SciPy, which applies only to Hermitian matrices.  Since
     1475            Hermitian is defined as a matrix equaling its conjugate-transpose,
     1476            for a matrix with real entries this property is equivalent to
     1477            being symmetric.  This algorithm can be significantly faster
     1478            than the ``'default'`` algorithm.
     1479
     1480        - ``'tol'`` - default: ``None`` - if set to a value other than
     1481          ``None`` this is interpreted as a small real number used to aid in
     1482          grouping eigenvalues that are numerically similar.  See the output
     1483          description for more information.
     1484
     1485        .. warning::
     1486
     1487        When using the ``'symmetric'`` or ``'hermitian'`` algorithms,
     1488        no check is made on the input matrix, and only the entries below,
     1489        and on, the main diagonal are employed in the computation.
     1490
     1491        OUTPUT:
     1492
     1493        Default output for a square matrix of size $n$ is a list of $n$
     1494        eigenvalues from the complex double field,
     1495        :class:`~sage.rings.complex_double.CDF`.  If the ``'symmetric'``
     1496        or ``'hermitian'`` algorithms are chosen, the returned eigenvalues
     1497        are from the real double field,
     1498        :class:`~sage.rings.real_double.RDF`.
     1499
     1500        If a tolerance is specified, an attempt is made to group eigenvalues
     1501        that are numerically similar.  The return is then a list of pairs,
     1502        where each pair is an eigenvalue followed by its multiplicity.
     1503        The eigenvalue reported is the mean of the eigenvalues computed,
     1504        and these eigenvalues are contained in an interval (or disk) whose
     1505        radius is less than ``5*tol`` for $n < 10,000$ in the worst case.
     1506
     1507        More precisely, for an $n\times n$ matrix, the diameter of the
     1508        interval containing similar eigenvalues could be as large as sum
     1509        of the reciprocals of the first $n$ integers times ``tol``.
     1510
     1511        .. warning::
     1512
     1513        Use caution when using the  ``tol`` parameter to group eigenvalues.
     1514        See the examples below to see how this can go wrong.
    14551515
    14561516        EXAMPLES::
    14571517
    14581518            sage: m = matrix(RDF, 2, 2, [1,2,3,4])
    1459             sage: m.eigenvalues()
    1460             [-0.372281323269, 5.37228132327]
    1461             sage: parent(m.eigenvalues()[0])
     1519            sage: ev = m.eigenvalues(); ev
     1520            [-0.372281323..., 5.37228132...]
     1521            sage: ev[0].parent()
    14621522            Complex Double Field
    14631523
    14641524            sage: m = matrix(RDF, 2, 2, [0,1,-1,0])
    1465             sage: m.eigenvalues()
     1525            sage: m.eigenvalues(algorithm='default')
    14661526            [1.0*I, -1.0*I]
    14671527
    14681528            sage: m = matrix(CDF, 2, 2, [I,1,-I,0])
    14691529            sage: m.eigenvalues()
    1470             [-0.624810533844 + 1.30024259022*I, 0.624810533844 - 0.30024259022*I]
     1530            [-0.624810533... + 1.30024259...*I, 0.624810533... - 0.30024259...*I]
     1531
     1532        The adjacency matrix of a graph will be symmetric, and the
     1533        eigenvalues will be real.  ::
     1534
     1535            sage: A = graphs.PetersenGraph().adjacency_matrix()
     1536            sage: A = A.change_ring(RDF)
     1537            sage: ev = A.eigenvalues(algorithm='symmetric'); ev
     1538            [-2.0, -2.0, -2.0, -2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 3.0]
     1539            sage: ev[0].parent()
     1540            Real Double Field
     1541
     1542        The matrix ``A`` is "random", but the construction of ``B``
     1543        provides a positive-definite Hermitian matrix.  Note that
     1544        the eigenvalues of a Hermitian matrix are real, and the
     1545        eigenvalues of a positive-definite matrix will be positive.  ::
     1546
     1547            sage: A = matrix([[ 4*I + 5,  8*I + 1,  7*I + 5, 3*I + 5],
     1548            ...               [ 7*I - 2, -4*I + 7, -2*I + 4, 8*I + 8],
     1549            ...               [-2*I + 1,  6*I + 6,  5*I + 5,  -I - 4],
     1550            ...               [ 5*I + 1,  6*I + 2,    I - 4, -I + 3]])
     1551            sage: B = (A*A.conjugate_transpose()).change_ring(CDF)
     1552            sage: ev = B.eigenvalues(algorithm='hermitian'); ev
     1553            [2.68144025..., 49.5167998..., 274.086188..., 390.71557...]
     1554            sage: ev[0].parent()
     1555            Real Double Field
     1556
     1557        A tolerance can be given to aid in grouping eigenvalues that
     1558        are similar numerically.  However, if the parameter is too small
     1559        it might split too finely.  Too large, and it can go wrong very
     1560        badly.  Use with care.  ::
     1561
     1562            sage: G = graphs.PetersenGraph()
     1563            sage: G.spectrum()
     1564            [3, 1, 1, 1, 1, 1, -2, -2, -2, -2]
     1565
     1566            sage: A = G.adjacency_matrix().change_ring(RDF)
     1567            sage: A.eigenvalues(algorithm='symmetric', tol=1.0e-5)
     1568            [(-2.0, 4), (1.0, 5), (3.0, 1)]
     1569
     1570            sage: A.eigenvalues(algorithm='symmetric', tol=1.0e-20)
     1571            [(-2.0, 1), (-2.0, 2), (-2.0, 1), (1.0, 1), (1.0, 1), (1.0, 1), (1.0, 1), (1.0, 1), (3.0, 1)]
     1572
     1573            sage: A.eigenvalues(algorithm='symmetric', tol=2.5)
     1574            [(-2.0, 4), (1.33333333333, 6)]
     1575
     1576        An (extreme) example of properly grouping similar eigenvalues.  ::
     1577
     1578            sage: G = graphs.HigmanSimsGraph()
     1579            sage: A = G.adjacency_matrix().change_ring(RDF)
     1580            sage: A.eigenvalues(algorithm='symmetric', tol=1.0e-5)
     1581            [(-8.0, 22), (2.0, 77), (22.0, 1)]
     1582
     1583        TESTS:
     1584
     1585        Testing bad input.  ::
     1586
     1587            sage: A = matrix(CDF, 2, range(4))
     1588            sage: A.eigenvalues(algorithm='junk')
     1589            Traceback (most recent call last):
     1590            ...
     1591            ValueError: algorithm must be 'default', 'symmetric', or 'hermitian', not junk
     1592
     1593            sage: A = matrix(CDF, 2, 3, range(6))
     1594            sage: A.eigenvalues()
     1595            Traceback (most recent call last):
     1596            ...
     1597            ValueError: matrix must be square, not 2 x 3
     1598
     1599            sage: A = matrix(CDF, 2, [1, 2, 3, 4*I])
     1600            sage: A.eigenvalues(algorithm='symmetric')
     1601            Traceback (most recent call last):
     1602            ...
     1603            TypeError: cannot apply symmetric algorithm to matrix with complex entries
     1604
     1605            sage: A = matrix(CDF, 2, 2, range(4))
     1606            sage: A.eigenvalues(tol='junk')
     1607            Traceback (most recent call last):
     1608            ...
     1609            TypeError: tolerance parameter must be a real number, not junk
     1610
     1611            sage: A = matrix(CDF, 2, 2, range(4))
     1612            sage: A.eigenvalues(tol=-0.01)
     1613            Traceback (most recent call last):
     1614            ...
     1615            ValueError: tolerance parameter must be positive, not -0.01
     1616
     1617        A very small matrix.  ::
    14711618
    14721619            sage: matrix(CDF,0,0).eigenvalues()
    14731620            []
    14741621        """
     1622        import sage.rings.real_double
     1623        import sage.rings.complex_double
     1624        import numpy
     1625        if not algorithm in ['default', 'symmetric', 'hermitian']:
     1626            msg = "algorithm must be 'default', 'symmetric', or 'hermitian', not {0}"
     1627            raise ValueError(msg.format(algorithm))
    14751628        if not self.is_square():
    1476             raise ArithmeticError("self must be a square matrix")
     1629            msg = 'matrix must be square, not {0} x {1}'
     1630            raise ValueError(msg.format(self.nrows(), self.ncols()))
     1631        if algorithm == 'symmetric' and self.base_ring() == sage.rings.complex_double.CDF:
     1632            try:
     1633                self = self.change_ring(sage.rings.real_double.RDF)  # check side effect
     1634            except TypeError:
     1635                raise TypeError('cannot apply symmetric algorithm to matrix with complex entries')
     1636        if algorithm == 'symmetric':
     1637            algorithm = 'hermitian'
     1638        multiplicity = not tol is None
     1639        if multiplicity:
     1640            try:
     1641                tol = float(tol)
     1642            except (ValueError, TypeError):
     1643                msg = 'tolerance parameter must be a real number, not {0}'
     1644                raise TypeError(msg.format(tol))
     1645            if tol < 0:
     1646                msg = 'tolerance parameter must be positive, not {0}'
     1647                raise ValueError(msg.format(tol))
     1648
    14771649        if self._nrows == 0:
    14781650            return []
    14791651        global scipy
    14801652        if scipy is None:
    14811653            import scipy
    14821654        import scipy.linalg
    1483         return [sage.rings.complex_double.CDF(x) for x in scipy.linalg.eigvals(self._matrix_numpy)]
     1655        if self._nrows == 0:
     1656            return []
     1657        global scipy
     1658        if scipy is None:
     1659            import scipy
     1660        import scipy.linalg
     1661        global numpy
     1662        if numpy is None:
     1663            import numpy
     1664        # generic eigenvalues, or real eigenvalues for Hermitian
     1665        if algorithm == 'default':
     1666            return_class = sage.rings.complex_double.CDF
     1667            evalues = scipy.linalg.eigvals(self._matrix_numpy)
     1668        elif algorithm=='hermitian':
     1669            return_class = sage.rings.real_double.RDF
     1670            evalues = scipy.linalg.eigh(self._matrix_numpy, eigvals_only=True)
     1671        if not multiplicity:
     1672            return [return_class(e) for e in evalues]
     1673        else:
     1674            # pairs in ev_group are
     1675            #   slot 0: the sum of "equal" eigenvalues, "s"
     1676            #   slot 1: number of eigenvalues in this sum, "m"
     1677            #   slot 2: average of these eigenvalues, "avg"
     1678            # we test if "new" eigenvalues are close to the group average
     1679            ev_group = []
     1680            for e in evalues:
     1681                location = None
     1682                best_fit = tol
     1683                for i in range(len(ev_group)):
     1684                    s, m, avg = ev_group[i]
     1685                    d = numpy.abs(avg - e)
     1686                    if d < best_fit:
     1687                        best_fit = d
     1688                        location = i
     1689                if location is None:
     1690                    ev_group.append([e, 1, e])
     1691                else:
     1692                    ev_group[location][0] += e
     1693                    ev_group[location][1] += 1
     1694                    ev_group[location][2] = ev_group[location][0]/ev_group[location][1]
     1695            return [(return_class(avg), m) for _, m, avg in ev_group]
    14841696
    14851697    def left_eigenvectors(self):
    14861698        r"""