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

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

Self conatined

  • sage/matrix/matrix_double_dense.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1264132569 28800
    # Node ID 4e6e9d0a0909019ad3823de75c2f9e872fe3b128
    # Parent  4498ef979fc42b03037970223b0873b07d63b96b
    [mq]: eigen
    
    diff -r 4498ef979fc4 -r 4e6e9d0a0909 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, [[-5, 3, 2, 8],[10, 2, 4, -2],[-1, -10, -10, -17],[-2, 7, 6, 13]])
     633            sage: spectrum = m.eigenspaces_left()
     634            sage: spectrum[0]
     635            (2.0, Vector space of degree 4 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.5 0.5 0.5 0.5])
    637638
    638             sage: e, v = es[0]
    639             sage: v = v.basis()[0]
    640             sage: a = v * m
    641             sage: b = e * v
    642             sage: diff = a.change_ring(CDF) - b
    643             sage: abs(abs(diff)) < 1e-10
     639            sage: e, V = spectrum[2]
     640            sage: v = V.basis()[0]
     641            sage: diff = (v*m).change_ring(CDF) - e*v
     642            sage: abs(abs(diff)) < 1e-14
    644643            True
    645             sage: diff # random -- very small numbers
    646             (-2.6645352591e-15, -7.1054273576e-15, -3.5527136788e-15)
    647644
    648         TESTS:
     645        TESTS::
     646       
    649647            sage: m.eigenspaces_left(algebraic_multiplicity=True)
    650648            Traceback (most recent call last):
    651649            ...
    652             ValueError: algebraic_multiplicity can only be False for double precision matrices           
     650            ValueError: algebraic_multiplicity must be set to False for double precision matrices
    653651        """
    654         # Raise an error if algebraic_multiplicity is True since that
    655         # would normally change the format of the return value.
     652        # For numerical values we leave decisions about
     653        # multiplicity to the calling routine
    656654        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()
     655            raise ValueError, "algebraic_multiplicity must be set to False for double precision matrices"
     656        spectrum = self.left_eigenvectors()
    661657        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)))
     658        for evalue,evectors,_ in spectrum:
     659            evector = evectors[0]
     660            espace = evector.parent().span_of_basis([evector],check=False)
     661            pairs.append((evalue, espace))
    665662        return pairs
    666663
    667664    left_eigenspaces = eigenspaces_left
     665   
     666    def eigenspaces_right(self, var='a', algebraic_multiplicity=False):
     667        r"""
     668        Computes the right eigenspaces of a matrix of double precision
     669        real or complex numbers (i.e. RDF or CDF).
     670
     671        INPUT:
     672
     673        - ``var`` - ignored for numerical matrices
     674        - ``algebraic_multiplicity`` - must be set to ``False``
     675          for numerical matrices, and will raise an error otherwise.
     676
     677        OUTPUT:
     678
     679        Return a list of pairs ``(e, V)`` where ``e`` is a (complex)
     680        eigenvalue and ``V`` is the associated right eigenspace as a
     681        vector space.
     682
     683        No attempt is made to determine if an eigenvalue has multiplicity
     684        greater than one, so all the eigenspaces returned have dimension one.
     685
     686        EXAMPLES::
     687
     688            sage: m = matrix(RDF, [[-9, -14, 19, -74],[-1, 2, 4, -11],[-4, -12, 6, -32],[0, -2, -1, 1]])
     689            sage: m
     690            [ -9.0 -14.0  19.0 -74.0]
     691            [ -1.0   2.0   4.0 -11.0]
     692            [ -4.0 -12.0   6.0 -32.0]
     693            [  0.0  -2.0  -1.0   1.0]
     694            sage: spectrum = m.eigenspaces_right()
     695            sage: spectrum[0]
     696            (2.0, Vector space of degree 4 and dimension 1 over Real Double Field
     697            User basis matrix:
     698            [ 0.258198889747 -0.516397779494  0.774596669241  0.258198889747])
     699
     700            sage: e, V = spectrum[2]
     701            sage: v = V.basis()[0]
     702            sage: diff = (m*v).change_ring(CDF) - e*v
     703            sage: abs(abs(diff)) < 1e-14
     704            True
     705
     706        TESTS::
     707       
     708            sage: m.eigenspaces_right(algebraic_multiplicity=True)
     709            Traceback (most recent call last):
     710            ...
     711            ValueError: algebraic_multiplicity must be set to False for double precision matrices
     712        """
     713        # For numerical values we leave decisions about
     714        # multiplicity to the calling routine
     715        if algebraic_multiplicity:
     716            raise ValueError, "algebraic_multiplicity must be set to False for double precision matrices"
     717        spectrum = self.right_eigenvectors()
     718        pairs = []
     719        for evalue,evectors,_ in spectrum:
     720            evector = evectors[0]
     721            espace = evector.parent().span_of_basis([evector],check=False)
     722            pairs.append((evalue, espace))
     723        return pairs
     724
     725    right_eigenspaces = eigenspaces_right
    668726
    669727    def eigenvalues(self):
    670         """
     728        r"""
    671729        Returns a list of the eigenvalues (with multiplicity)
    672730        of this matrix.  The returned eigenvalues are elements of CDF.
    673731
     
    685743            sage: m = matrix(CDF, 2, 2, [I,1,-I,0])
    686744            sage: m.eigenvalues()
    687745            [-0.624810533844 + 1.30024259022*I, 0.624810533844 - 0.30024259022*I]
     746
    688747            sage: matrix(CDF,0,0).eigenvalues()
    689748            []
    690749        """
     
    699758        return [sage.rings.complex_double.CDF(x) for x in scipy.linalg.eigvals(self._matrix_numpy)]
    700759
    701760    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.
     761        r"""
     762        Compute the left eigenvectors of a matrix of double precision
     763        real or complex numbers (i.e. RDF or CDF).
    706764
    707765        OUTPUT:
    708              eigenvalues -- as a list
    709              corresponding eigenvectors -- as an RDF matrix whose rows
    710                            are the eigenvectors.
     766        Returns a list of triples, each of the form ``(e,[v],1)``,
     767        where ``e`` is the eigenvalue, and ``v`` is an associated
     768        left eigenvector.  If the matrix is of size `n`, then there are
     769        `n` triples.  Values are computed with the SciPy library.
    711770
    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
     771        The format of this output is designed to match the format
     772        for exact results.  However, since matrices here have numerical
     773        entries, the resulting eigenvalues will also be numerical.  No
     774        attempt is made to determine if two eigenvalues are equal, or if
     775        eigenvalues might actually be zero.  So the algebraic multiplicity
     776        of each eigenvalue is reported as 1.  Decisions about equal
     777        eigenvalues or zero eigenvalues should be addressed in the
     778        calling routine.
     779
     780        EXAMPLES::
     781
     782            sage: m = matrix(RDF, [[-5, 3, 2, 8],[10, 2, 4, -2],[-1, -10, -10, -17],[-2, 7, 6, 13]])
     783            sage: m
     784            [ -5.0   3.0   2.0   8.0]
     785            [ 10.0   2.0   4.0  -2.0]
     786            [ -1.0 -10.0 -10.0 -17.0]
     787            [ -2.0   7.0   6.0  13.0]
     788            sage: spectrum = m.left_eigenvectors()
     789            sage: spectrum[0]
     790            (2.0, [(0.5, 0.5, 0.5, 0.5)], 1)
     791            sage: spectrum[1]
     792            (1.0, [(-0.615457454897, -0.492365963917, -0.492365963917, -0.369274472938)], 1)
     793            sage: spectrum[2]
     794            (-2.0, [(-0.800640769025, -0.32025630761, -0.480384461415, -0.160128153805)], 1)
     795            sage: spectrum[3]
     796            (-1.0, [(0.316227766017, 0.316227766017, 0.632455532034, 0.632455532034)], 1)
    723797        """
    724798        if not self.is_square():
    725799            raise ArithmeticError, "self must be a square matrix"
     
    730804            import scipy
    731805        import scipy.linalg
    732806        v,eig = scipy.linalg.eig(self._matrix_numpy, right=False, left=True)
     807        # scipy puts eigenvectors in columns, we will extract from rows
    733808        eig = matrix(eig.T)
    734         return ([sage.rings.complex_double.CDF(x) for x in v], eig)
     809        return [(sage.rings.complex_double.CDF(v[i]), [eig[i]], 1) for i in range(len(v))]
    735810
    736811    eigenvectors_left = left_eigenvectors
    737812   
    738813    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.
     814        r"""
     815        Compute the right eigenvectors of a matrix of double precision
     816        real or complex numbers (i.e. RDF or CDF).
    743817
    744818        OUTPUT:
    745              eigenvalues -- as a list
    746              corresponding eigenvectors -- as an RDF matrix whose columns
    747                            are the eigenvectors.
     819        Returns a list of triples, each of the form ``(e,[v],1)``,
     820        where ``e`` is the eigenvalue, and ``v`` is an associated
     821        right eigenvector.  If the matrix is of size `n`, then there
     822        are `n` triples.  Values are computed with the SciPy library.
    748823
    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         """
     824        The format of this output is designed to match the format
     825        for exact results.  However, since matrices here have numerical
     826        entries, the resulting eigenvalues will also be numerical.  No
     827        attempt is made to determine if two eigenvalues are equal, or if
     828        eigenvalues might actually be zero.  So the algebraic multiplicity
     829        of each eigenvalue is reported as 1.  Decisions about equal
     830        eigenvalues or zero eigenvalues should be addressed in the
     831        calling routine.
     832
     833        EXAMPLES::
     834       
     835            sage: m = matrix(RDF, [[-9, -14, 19, -74],[-1, 2, 4, -11],[-4, -12, 6, -32],[0, -2, -1, 1]])
     836            sage: m
     837            [ -9.0 -14.0  19.0 -74.0]
     838            [ -1.0   2.0   4.0 -11.0]
     839            [ -4.0 -12.0   6.0 -32.0]
     840            [  0.0  -2.0  -1.0   1.0]
     841            sage: spectrum = m.right_eigenvectors()
     842            sage: spectrum[0]
     843            (2.0, [(0.258198889747, -0.516397779494, 0.774596669241, 0.258198889747)], 1)
     844            sage: spectrum[1]
     845            (1.0, [(-0.547722557505, 0.36514837167, -0.73029674334, -0.182574185835)], 1)
     846            sage: spectrum[2]
     847            (-2.0, [(0.693375245282, -0.138675049056, 0.693375245282, 0.138675049056)], 1)
     848            sage: spectrum[3]
     849            (-1.0, [(0.426401432711, -0.213200716356, 0.852802865422, 0.213200716356)], 1)
     850            """
    761851        if not self.is_square():
    762852            raise ArithmeticError, "self must be a square matrix"
    763853        if self._nrows == 0:
     
    767857            import scipy
    768858        import scipy.linalg
    769859        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))
     860        # scipy puts eigenvectors in columns, we will extract from rows
     861        eig = matrix(eig.T)
     862        return [(sage.rings.complex_double.CDF(v[i]), [eig[i]], 1) for i in range(len(v))]
    771863
    772864    eigenvectors_right = right_eigenvectors
    773865
  • sage/modular/modform/numerical.py

    diff -r 4498ef979fc4 -r 4e6e9d0a0909 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.