Ticket #4756: trac_4756-double-eigen.2.patch

File trac_4756-double-eigen.2.patch, 13.6 KB (added by Rob Beezer, 13 years ago)

Self-contained patch, apply only this.

  • sage/matrix/matrix_double_dense.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1264132569 28800
    # Node ID 6c4190c4ba571ea4615296cce5a44a3d87c2cd0c
    # Parent  4498ef979fc42b03037970223b0873b07d63b96b
    [mq]: eigen
    
    diff -r 4498ef979fc4 -r 6c4190c4ba57 sage/matrix/matrix_double_dense.pyx
    a b  
    609609
    610610    def eigenspaces_left(self, var='a', algebraic_multiplicity=False):
    611611        r"""
    612         Return a list of pairs (e, V) where e runs through all complex
    613         eigenvalues of this matrix, and V is the corresponding
    614         left eigenspace (always a 1-dimensional complex vector space).
     612        Computes the left eigenspaces of a matrix of double precision
     613        real or complex numbers (i.e. RDF or CDF).
    615614
    616         INPUT
    617             Both the var and the algebraic_multiplicity arguments are
    618             ignored.  They do not make much sense in a numerical
    619             situation.
     615        INPUT:
    620616
    621         EXAMPLES:
    622             sage: m = matrix(RDF, 3, range(9)); m
    623             [0.0 1.0 2.0]
    624             [3.0 4.0 5.0]
    625             [6.0 7.0 8.0]
    626             sage: es = m.eigenspaces_left()
    627             sage: es # random
    628             [(13.3484692283, Vector space of degree 3 and dimension 1 over Real Double Field
     617        - ``var`` - ignored for numerical matrices
     618        - ``algebraic_multiplicity`` - must be set to ``False``
     619          for numerical matrices, and will raise an error otherwise.
     620
     621        OUTPUT:
     622
     623        Return a list of pairs ``(e, V)`` where ``e`` is a (complex)
     624        eigenvalue and ``V`` is the associated left eigenspace as a
     625        vector space.
     626
     627        No attempt is made to determine if an eigenvalue has multiplicity
     628        greater than one, so all the eigenspaces returned have dimension one.
     629
     630        EXAMPLES::
     631       
     632            sage: m = matrix(RDF, 3, range(9))
     633            sage: spectrum = m.eigenspaces_left()
     634            sage: spectrum[0]
     635            (13.3484692283, Vector space of degree 3 and dimension 1 over Real Double Field
    629636            User basis matrix:
    630             [-0.440242867236 -0.567868371314 -0.695493875393]),
    631             (-1.34846922835, Vector space of degree 3 and dimension 1 over Real Double Field
    632             User basis matrix:
    633             [-0.897878732262 -0.278434036822  0.341010658618]),
    634             (-9.10854412047e-16, Vector space of degree 3 and dimension 1 over Real Double Field
    635             User basis matrix:
    636             [ 0.408248290464 -0.816496580928  0.408248290464])]
     637            [0.440242867236 0.567868371314 0.695493875393])
    637638
    638             sage: e, v = es[0]
     639            sage: e, v = spectrum[0]
    639640            sage: v = v.basis()[0]
    640641            sage: a = v * m
    641642            sage: b = e * v
     
    645646            sage: diff # random -- very small numbers
    646647            (-2.6645352591e-15, -7.1054273576e-15, -3.5527136788e-15)
    647648
    648         TESTS:
     649        TESTS::
     650       
    649651            sage: m.eigenspaces_left(algebraic_multiplicity=True)
    650652            Traceback (most recent call last):
    651653            ...
    652             ValueError: algebraic_multiplicity can only be False for double precision matrices           
     654            ValueError: algebraic_multiplicity must be set to False for double precision matrices
    653655        """
    654         # Raise an error if algebraic_multiplicity is True since that
    655         # would normally change the format of the return value.
     656        # For numerical values we leave decisions about
     657        # multiplicity to the calling routine
    656658        if algebraic_multiplicity:
    657             raise ValueError, "algebraic_multiplicity can only be False for double precision matrices"
    658            
    659         e, v = self.left_eigenvectors()
    660         v = v.rows()
     659            raise ValueError, "algebraic_multiplicity must be set to False for double precision matrices"
     660        spectrum = self.left_eigenvectors()
    661661        pairs = []
    662         for l from 0<=l<len(e):
    663             c = v[l]
    664             pairs.append((e[l], c.parent().span_of_basis([c], check=False)))
     662        for evalue,evectors,_ in spectrum:
     663            evector = evectors[0]
     664            espace = evector.parent().span_of_basis([evector],check=False)
     665            pairs.append((evalue, espace))
    665666        return pairs
    666667
    667668    left_eigenspaces = eigenspaces_left
     669   
     670    def eigenspaces_right(self, var='a', algebraic_multiplicity=False):
     671        r"""
     672        Computes the right eigenspaces of a matrix of double precision
     673        real or complex numbers (i.e. RDF or CDF).
     674
     675        INPUT:
     676
     677        - ``var`` - ignored for numerical matrices
     678        - ``algebraic_multiplicity`` - must be set to ``False``
     679          for numerical matrices, and will raise an error otherwise.
     680
     681        OUTPUT:
     682
     683        Return a list of pairs ``(e, V)`` where ``e`` is a (complex)
     684        eigenvalue and ``V`` is the associated right eigenspace as a
     685        vector space.
     686
     687        No attempt is made to determine if an eigenvalue has multiplicity
     688        greater than one, so all the eigenspaces returned have dimension one.
     689
     690        EXAMPLES::
     691
     692            sage: m = matrix(RDF, 3, range(9))
     693            sage: spectrum = m.eigenspaces_right()
     694            sage: spectrum[0]
     695            (13.3484692283, Vector space of degree 3 and dimension 1 over Real Double Field
     696            User basis matrix:
     697            [0.164763817282 0.505774475901 0.846785134519])
     698
     699            sage: e, v = spectrum[0]
     700            sage: v = v.basis()[0]
     701            sage: a = m * v
     702            sage: b = e * v
     703            sage: diff = a.change_ring(CDF) - b
     704            sage: abs(abs(diff)) < 1e-10
     705            True
     706            sage: diff # random -- very small numbers
     707            (2.22044604925e-15, 9.7699626167e-15, 7.1054273576e-15)
     708
     709        TESTS::
     710       
     711            sage: m.eigenspaces_right(algebraic_multiplicity=True)
     712            Traceback (most recent call last):
     713            ...
     714            ValueError: algebraic_multiplicity must be set to False for double precision matrices
     715        """
     716        # For numerical values we leave decisions about
     717        # multiplicity to the calling routine
     718        if algebraic_multiplicity:
     719            raise ValueError, "algebraic_multiplicity must be set to False for double precision matrices"
     720        spectrum = self.right_eigenvectors()
     721        pairs = []
     722        for evalue,evectors,_ in spectrum:
     723            evector = evectors[0]
     724            espace = evector.parent().span_of_basis([evector],check=False)
     725            pairs.append((evalue, espace))
     726        return pairs
     727
     728    right_eigenspaces = eigenspaces_right
    668729
    669730    def eigenvalues(self):
    670         """
     731        r"""
    671732        Returns a list of the eigenvalues (with multiplicity)
    672733        of this matrix.  The returned eigenvalues are elements of CDF.
    673734
     
    685746            sage: m = matrix(CDF, 2, 2, [I,1,-I,0])
    686747            sage: m.eigenvalues()
    687748            [-0.624810533844 + 1.30024259022*I, 0.624810533844 - 0.30024259022*I]
     749
    688750            sage: matrix(CDF,0,0).eigenvalues()
    689751            []
    690752        """
     
    699761        return [sage.rings.complex_double.CDF(x) for x in scipy.linalg.eigvals(self._matrix_numpy)]
    700762
    701763    def left_eigenvectors(self):
    702         """
    703         Computes the eigenvalues and *left* eigenvectors of this
    704         matrix m acting *from the right*.  I.e., vectors v such that
    705         v*m = lambda*v.
     764        r"""
     765        Compute the left eigenvectors of a matrix of double precision
     766        real or complex numbers (i.e. RDF or CDF).
    706767
    707768        OUTPUT:
    708              eigenvalues -- as a list
    709              corresponding eigenvectors -- as an RDF matrix whose rows
    710                            are the eigenvectors.
     769        Returns a list of triples, each of the form ``(e,[v],1)``,
     770        where ``e`` is the eigenvalue, and ``v`` is an associated
     771        left eigenvector.  If the matrix is of size `n`, then there are
     772        `n` triples.  Values are computed with the SciPy library.
    711773
    712         EXAMPLES:
    713             sage: m = Matrix(RDF, 3, range(9))
    714             sage: es = m.left_eigenvectors()
    715             sage: es # random-ish platform-dependent output (low order digits)
    716             ([13.3484692283, -1.34846922835, -9.10854412047e-16],
    717             [-0.440242867236 -0.567868371314 -0.695493875393]
    718             [-0.897878732262 -0.278434036822  0.341010658618]
    719             [ 0.408248290464 -0.816496580928  0.408248290464])
    720             sage: e, v = es; e = e[0]; v = v[0]
    721             sage: abs(abs(e*v - v*m)) < 1e-10
    722             True
     774        The format of this output is designed to match the format
     775        for exact results.  However, since matrices here have numerical
     776        entries, the resulting eigenvalues will also be numerical.  No
     777        attempt is made to determine if two eigenvalues are equal, or if
     778        eigenvalues might actually be zero.  So the algebraic multiplicity
     779        of each eigenvalue is reported as 1.  Decisions about equal
     780        eigenvalues or zero eigenvalues should be addressed in the
     781        calling routine.
     782
     783        EXAMPLES::
     784
     785            sage: m = matrix(RDF, 3, range(9))
     786            sage: m.left_eigenvectors()
     787            [(13.3484692283, [(0.440242867236, 0.567868371314, 0.695493875393)], 1), (-1.34846922835, [(0.897878732262, 0.278434036822, -0.341010658618)], 1), (-6.2265089409e-16, [(0.408248290464, -0.816496580928, 0.408248290464)], 1)]
    723788        """
    724789        if not self.is_square():
    725790            raise ArithmeticError, "self must be a square matrix"
     
    730795            import scipy
    731796        import scipy.linalg
    732797        v,eig = scipy.linalg.eig(self._matrix_numpy, right=False, left=True)
     798        # scipy puts eigenvectors in columns, we will extract from rows
    733799        eig = matrix(eig.T)
    734         return ([sage.rings.complex_double.CDF(x) for x in v], eig)
     800        return [(sage.rings.complex_double.CDF(v[i]), [eig[i]], 1) for i in range(len(v))]
    735801
    736802    eigenvectors_left = left_eigenvectors
    737803   
    738804    def right_eigenvectors(self):
    739         """
    740         Computes the eigenvalues and *right* eigenvectors of this
    741         matrix m acting *from the left*.  I.e., vectors v such that
    742         m * v = lambda*v, where v is viewed as a column vector.
     805        r"""
     806        Compute the right eigenvectors of a matrix of double precision
     807        real or complex numbers (i.e. RDF or CDF).
    743808
    744809        OUTPUT:
    745              eigenvalues -- as a list
    746              corresponding eigenvectors -- as an RDF matrix whose columns
    747                            are the eigenvectors.
     810        Returns a list of triples, each of the form ``(e,[v],1)``,
     811        where ``e`` is the eigenvalue, and ``v`` is an associated
     812        right eigenvector.  If the matrix is of size `n`, then there
     813        are `n` triples.  Values are computed with the SciPy library.
    748814
    749         EXAMPLES:
    750             sage: m = Matrix(RDF, 3, range(9))
    751             sage: evals,evecs = m.right_eigenvectors()
    752             sage: evals # random low-order bits
    753             [13.3484692283, -1.34846922835, -1.1327908706e-15]
    754             sage: evecs # random low-order bits
    755             [ 0.164763817282  0.799699663112  0.408248290464]
    756             [ 0.505774475901  0.104205787719 -0.816496580928]
    757             [ 0.846785134519 -0.591288087674  0.408248290464]
    758             sage: max([max(m*evec - evec*eval) for eval,evec in zip(evals, evecs.columns())]) < 1e-14
    759             True
    760         """
     815        The format of this output is designed to match the format
     816        for exact results.  However, since matrices here have numerical
     817        entries, the resulting eigenvalues will also be numerical.  No
     818        attempt is made to determine if two eigenvalues are equal, or if
     819        eigenvalues might actually be zero.  So the algebraic multiplicity
     820        of each eigenvalue is reported as 1.  Decisions about equal
     821        eigenvalues or zero eigenvalues should be addressed in the
     822        calling routine.
     823
     824        EXAMPLES::
     825       
     826            sage: m = matrix(RDF, 3, range(9))
     827            sage: m.right_eigenvectors()
     828            [(13.3484692283, [(0.164763817282, 0.505774475901, 0.846785134519)], 1), (-1.34846922835, [(0.799699663112, 0.104205787719, -0.591288087674)], 1), (-6.2265089409e-16, [(0.408248290464, -0.816496580928, 0.408248290464)], 1)]
     829            """
    761830        if not self.is_square():
    762831            raise ArithmeticError, "self must be a square matrix"
    763832        if self._nrows == 0:
     
    767836            import scipy
    768837        import scipy.linalg
    769838        v,eig = scipy.linalg.eig(self._matrix_numpy, right=True, left=False)
    770         return ([sage.rings.complex_double.CDF(x) for x in v], matrix(eig))
     839        # scipy puts eigenvectors in columns, we will extract from rows
     840        eig = matrix(eig.T)
     841        return [(sage.rings.complex_double.CDF(v[i]), [eig[i]], 1) for i in range(len(v))]
    771842
    772843    eigenvectors_right = right_eigenvectors
    773844
  • sage/modular/modform/numerical.py

    diff -r 4498ef979fc4 -r 6c4190c4ba57 sage/modular/modform/numerical.py
    a b  
    206206            t += randint(-50,50)*M.T(p).matrix()
    207207
    208208        self._hecke_matrix = t
    209         evals, B = t.change_ring(CDF).right_eigenvectors()
     209        # evals, B = t.change_ring(CDF).right_eigenvectors()
     210        from sage.matrix.constructor import matrix
     211        spectrum = t.change_ring(CDF).right_eigenvectors()
     212        evals = [spectrum[i][0] for i in range(len(spectrum))]
     213        B = matrix(CDF, [spectrum[i][1][0] for i in range(len(spectrum))]).transpose()
    210214       
    211215        # Find the eigenvalues that occur with multiplicity 1 up
    212216        # to the given eps.