Ticket #4756: trac_4756_eigenmatrices_double.patch

File trac_4756_eigenmatrices_double.patch, 7.6 KB (added by Rob Beezer, 13 years ago)
  • sage/matrix/matrix_double_dense.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1263866033 28800
    # Node ID 5488e7f6e3109d5adedea1527e9be0a8f0d82f89
    # Parent  33c04bb5aea0255c4fac02782a9974deaf33f1ce
    [mq]: eigen
    
    diff -r 33c04bb5aea0 -r 5488e7f6e310 sage/matrix/matrix_double_dense.pyx
    a b  
    611611
    612612    def eigenspaces_left(self, var='a', algebraic_multiplicity=False):
    613613        r"""
    614         Return a list of pairs (e, V) where e runs through all complex
    615         eigenvalues of this matrix, and V is the corresponding
    616         left eigenspace (always a 1-dimensional complex vector space).
     614        Return a list of pairs ``(e, V)`` where ``e`` is a (complex)
     615        eigenvalue and ``V`` is the associated left eigenspace.
    617616
    618         INPUT
    619             Both the var and the algebraic_multiplicity arguments are
    620             ignored.  They do not make much sense in a numerical
    621             situation.
     617        INPUT:
     618        - ``var`` - ignored for numerical matrices
     619        - ``algebraic_multiplicity`` - must be set to ``False``
     620          for numerical matrices, and will raise an error otherwise
    622621
    623622        EXAMPLES:
    624623            sage: m = matrix(RDF, 3, range(9)); m
     
    651650            sage: m.eigenspaces_left(algebraic_multiplicity=True)
    652651            Traceback (most recent call last):
    653652            ...
    654             ValueError: algebraic_multiplicity can only be False for double precision matrices           
     653            ValueError: algebraic_multiplicity must be set to False for double precision matrices
    655654        """
    656655        # Raise an error if algebraic_multiplicity is True since that
    657656        # would normally change the format of the return value.
    658657        if algebraic_multiplicity:
    659             raise ValueError, "algebraic_multiplicity can only be False for double precision matrices"
     658            raise ValueError, "algebraic_multiplicity must be set to False for double precision matrices"
    660659           
    661         e, v = self.left_eigenvectors()
    662         v = v.rows()
     660        spectrum = self.left_eigenvectors()
    663661        pairs = []
    664         for l from 0<=l<len(e):
    665             c = v[l]
    666             pairs.append((e[l], c.parent().span_of_basis([c], check=False)))
     662        for k in range(len(spectrum)):
     663            evalue = spectrum[k][0]
     664            evector = spectrum[k][1][0]
     665            pairs.append((evalue, evector.parent().span_of_basis([evector],check=False)))
    667666        return pairs
    668667
    669668    left_eigenspaces = eigenspaces_left
     
    701700        return [sage.rings.complex_double.CDF(x) for x in scipy.linalg.eigvals(self._matrix_numpy)]
    702701
    703702    def left_eigenvectors(self):
    704         """
    705         Computes the eigenvalues and *left* eigenvectors of this
    706         matrix m acting *from the right*.  I.e., vectors v such that
    707         v*m = lambda*v.
     703        r"""
     704        Compute the left eigenvectors of a matrix over CDF.
    708705
    709706        OUTPUT:
    710              eigenvalues -- as a list
    711              corresponding eigenvectors -- as an RDF matrix whose rows
    712                            are the eigenvectors.
     707        Returns a list of triples, each of the form (e,[v],1),
     708        where ``e`` is the eigenvalue, and ``v`` is an associated
     709        eigenvector.  If the matrix is of size `n`, then there are
     710        that many triples.
     711
     712        ..note
     713        The format of this output is designed to match the format
     714        for exact results.  However, since matrices here have numerical
     715        entries, the resulting eigenvalues will also be numerical.  No
     716        attempt is made to determine if two eigenvalues are equal, or if
     717        eigenvalues might actually be zero.  So the algebraic multiplicity
     718        of each eigenvalue is reported as 1.
    713719
    714720        EXAMPLES:
    715721            sage: m = Matrix(RDF, 3, range(9))
    716722            sage: es = m.left_eigenvectors()
    717723            sage: es # random-ish platform-dependent output (low order digits)
    718             ([13.3484692283, -1.34846922835, -9.10854412047e-16],
    719             [-0.440242867236 -0.567868371314 -0.695493875393]
    720             [-0.897878732262 -0.278434036822  0.341010658618]
    721             [ 0.408248290464 -0.816496580928  0.408248290464])
    722             sage: e, v = es; e = e[0]; v = v[0]
    723             sage: abs(abs(e*v - v*m)) < 1e-10
    724             True
     724            [(13.3484692283, [(0.440242867236, 0.567868371314, 0.695493875393)], 1), (-1.34846922835, [(0.897878732262, 0.278434036822, -0.341010658618)], 1), (-1.1327908706e-15, [(0.408248290464, -0.816496580928, 0.408248290464)], 1)]
    725725        """
    726726        if not self.is_square():
    727727            raise ArithmeticError, "self must be a square matrix"
     
    733733        import scipy.linalg
    734734        v,eig = scipy.linalg.eig(self._matrix_numpy, right=False, left=True)
    735735        eig = matrix(eig.T)
    736         return ([sage.rings.complex_double.CDF(x) for x in v], eig)
     736        return [(sage.rings.complex_double.CDF(v[i]), [eig[i]], 1) for i in range(len(v))]
    737737
    738738    eigenvectors_left = left_eigenvectors
    739739   
    740740    def right_eigenvectors(self):
    741         """
    742         Computes the eigenvalues and *right* eigenvectors of this
    743         matrix m acting *from the left*.  I.e., vectors v such that
    744         m * v = lambda*v, where v is viewed as a column vector.
     741        r"""
     742        Compute the right eigenvectors of a matrix over CDF.
    745743
    746744        OUTPUT:
    747              eigenvalues -- as a list
    748              corresponding eigenvectors -- as an RDF matrix whose columns
    749                            are the eigenvectors.
     745        Returns a list of triples, each of the form (e,[v],1),
     746        where ``e`` is the eigenvalue, and ``v`` is an associated
     747        eigenvector.  If the matrix is of size `n`, then there are
     748        that many triples.
     749
     750        ..note
     751        The format of this output is designed to match the format
     752        for exact results.  However, since matrices here have numerical
     753        entries, the resulting eigenvalues will also be numerical.  No
     754        attempt is made to determine if two eigenvalues are equal, or if
     755        eigenvalues might actually be zero.  So the algebraic multiplicity
     756        of each eigenvalue is reported as 1.
    750757
    751758        EXAMPLES:
    752759            sage: m = Matrix(RDF, 3, range(9))
    753             sage: evals,evecs = m.right_eigenvectors()
    754             sage: evals # random low-order bits
    755             [13.3484692283, -1.34846922835, -1.1327908706e-15]
    756             sage: evecs # random low-order bits
    757             [ 0.164763817282  0.799699663112  0.408248290464]
    758             [ 0.505774475901  0.104205787719 -0.816496580928]
    759             [ 0.846785134519 -0.591288087674  0.408248290464]
    760             sage: max([max(m*evec - evec*eval) for eval,evec in zip(evals, evecs.columns())]) < 1e-14
    761             True
    762         """
     760            sage: m.right_eigenvectors()
     761            [(13.3484692283, [(0.164763817282, 0.799699663112, 0.408248290464)], 1), (-1.34846922835, [(0.505774475901, 0.104205787719, -0.816496580928)], 1), (-1.1327908706e-15, [(0.846785134519, -0.591288087674, 0.408248290464)], 1)]
     762            """
    763763        if not self.is_square():
    764764            raise ArithmeticError, "self must be a square matrix"
    765765        if self._nrows == 0:
     
    769769            import scipy
    770770        import scipy.linalg
    771771        v,eig = scipy.linalg.eig(self._matrix_numpy, right=True, left=False)
    772         return ([sage.rings.complex_double.CDF(x) for x in v], matrix(eig))
     772        eig = matrix(eig)
     773        return [(sage.rings.complex_double.CDF(v[i]), [eig[i]], 1) for i in range(len(v))]
    773774
    774775    eigenvectors_right = right_eigenvectors
    775776