Ticket #10837: trac_10837-norms-condition-CDF.patch

File trac_10837-norms-condition-CDF.patch, 16.9 KB (added by Rob Beezer, 12 years ago)

  • sage/matrix/matrix2.pyx 

    # HG changeset patch
# User Rob Beezer <beezer@ups.edu>
# Date 1298510872 28800
# Node ID 5c31405d7399ff39fbf5cc3cdb966c0793704908
# Parent  3c67675a1df921b9a9fa0ef007fc54c1a1d1e8f6
10837: matrix and vector norms, condition number over RDF/CDF

diff -r 3c67675a1df9 -r 5c31405d7399 sage/matrix/matrix2.pyx
    a b  
    68546854            sage: A.norm()
    68556855            15.0
    68566856       
    6857         ::
    6858        
    6859             sage: A = matrix(CDF, 2, 3, [3*I,4,1-I,1,2,0])
     6857        Norms of numerical matrices over high-precision reals are computed by this routine.
     6858        Faster routines for double precision entries from `RDF` or `CDF` are provided by
     6859        the :class:`~sage.matrix.matrix_double_dense.Matrix_double_dense` class.  ::
     6860       
     6861            sage: A = matrix(RR, 2, 3, [3*I,4,1-I,1,2,0])
    68606862            sage: A.norm('frob')
    68616863            5.65685424949
    68626864            sage: A.norm(2)

  • sage/matrix/matrix_double_dense.pyx 

    diff -r 3c67675a1df9 -r 5c31405d7399 sage/matrix/matrix_double_dense.pyx
    a b  
    3636
    3737import sage.rings.real_double
    3838import sage.rings.complex_double
     39
    3940from matrix cimport Matrix
    4041from sage.structure.element cimport ModuleElement,Vector
    4142from constructor import matrix
    … …  
    537538    #    compute_LU(self)   
    538539    #
    539540    ########################################################################
    540            
     541
     542    def condition(self, p='frob'):
     543        r"""
     544        Returns the condition number of a square nonsingular matrix.
     545
     546        Roughly speaking, this is a measure of how sensitive
     547        the matrix is to round-off errors in numerical computations.
     548        The minimum possible value is 1.0, and larger numbers indicate
     549        greater sensitivity.
     550
     551        INPUT:
     552
     553        - ``p`` - default: 'frob' - controls which norm is used
     554        to compute the condition number, allowable values are
     555        'frob' (for the Frobenius norm), integers -2, -1, 1, 2,
     556        positive and negative infinity. See output discussion
     557        for specifics.
     558
     559        OUTPUT:
     560
     561        The condition number of a matrix is the product of a norm
     562        of the matrix times the norm of the inverse of the matrix.
     563        This requires that the matrix be square and invertible
     564        (nonsingular, full rank).
     565
     566        Returned value is a double precision floating point value
     567        in ``RDF``, or ``Infinity``.  Row and column sums described below are
     568        sums of the absolute values of the entries, where the
     569        absolute value of the complex number `a+bi` is `\sqrt{a^2+b^2}`.
     570        Singular values are the "diagonal" entries of the "S" matrix in
     571        the singular value decomposition.
     572
     573        - ``p = 'frob'``: the default norm employed in computing
     574          the condition number, the Frobenius norm, which for a
     575          matrix `A=(a_{ij})` computes
     576
     577          .. math::
     578
     579                \left(\sum_{i,j}\left\lvert{a_{i,j}}\right\rvert^2\right)^{1/2}
     580
     581        - ``p = Infinity`` or ``p = oo``: the maximum row sum.
     582        - ``p = -Infinity`` or ``p = -oo``: the minimum colum sum.
     583        - ``p = 1``: the maximum column sum.
     584        - ``p = -1``: the minimum column sum.
     585        - ``p = 2``: the 2-norm, equal to the maximum singular value.
     586        - ``p = -2``: the minimum singular value.
     587
     588        ALGORITHM:
     589
     590        Computation is performed by the ``cond()`` function of
     591        the SciPy/NumPy library.
     592
     593        EXAMPLES:
     594
     595        First over the reals.  ::
     596
     597            sage: A = matrix(RDF, 4, [(1/4)*x^3 for x in range(16)]); A
     598            [   0.0   0.25    2.0   6.75]
     599            [  16.0  31.25   54.0  85.75]
     600            [ 128.0 182.25  250.0 332.75]
     601            [ 432.0 549.25  686.0 843.75]
     602            sage: A.condition()
     603            9923.88955...
     604            sage: A.condition(p='frob')
     605            9923.88955...
     606            sage: A.condition(p=Infinity)
     607            22738.5
     608            sage: A.condition(p=-Infinity)
     609            17.5
     610            sage: A.condition(p=1)
     611            12139.21...
     612            sage: A.condition(p=-1)
     613            550.0
     614            sage: A.condition(p=2)
     615            9897.8088...
     616            sage: A.condition(p=-2)
     617            0.000101032462...
     618
     619        And over the complex numbers.  ::
     620
     621            sage: B = matrix(CDF, 3, [x + x^2*I for x in range(9)]); B
     622            [           0  1.0 + 1.0*I  2.0 + 4.0*I]
     623            [ 3.0 + 9.0*I 4.0 + 16.0*I 5.0 + 25.0*I]
     624            [6.0 + 36.0*I 7.0 + 49.0*I 8.0 + 64.0*I]
     625            sage: B.condition()
     626            203.851798...
     627            sage: B.condition(p='frob')
     628            203.851798...
     629            sage: B.condition(p=Infinity)
     630            369.55630...
     631            sage: B.condition(p=-Infinity)
     632            5.46112969...
     633            sage: B.condition(p=1)
     634            289.251481...
     635            sage: B.condition(p=-1)
     636            20.4566639...
     637            sage: B.condition(p=2)
     638            202.653543...
     639            sage: B.condition(p=-2)
     640            0.00493453005...
     641
     642        Hilbert matrices are famously ill-conditioned, while
     643        an identity matrix can hit the minimum with the right norm.  ::
     644
     645            sage: A = matrix(RDF, 10, [1/(i+j+1) for i in range(10) for j in range(10)])
     646            sage: A.condition()
     647            1.63346888329e+13
     648            sage: id = identity_matrix(CDF, 10)
     649            sage: id.condition(p=1)
     650            1.0
     651
     652        Return values are in `RDF`.  ::
     653
     654            sage: A = matrix(CDF, 2, range(1,5))
     655            sage: A.condition() in RDF
     656            True
     657
     658        Rectangular and singular matrices raise errors.  ::
     659
     660            sage: A = matrix(RDF, 2, 3, range(6))
     661            sage: A.condition()
     662            Traceback (most recent call last):
     663            ...
     664            TypeError: matrix must be square, not 2 x 3
     665
     666            sage: A = matrix(QQ, 5, range(25))
     667            sage: A.is_singular()
     668            True
     669            sage: B = A.change_ring(CDF)
     670            sage: B.condition()
     671            Traceback (most recent call last):
     672            ...
     673            LinAlgError: Singular matrix
     674
     675        Improper values of ``p`` are caught.  ::
     676
     677            sage: A = matrix(CDF, 2, range(1,5))
     678            sage: A.condition(p='bogus')
     679            Traceback (most recent call last):
     680            ...
     681            ValueError: condition number 'p' must be +/- infinity, 'frob' or an integer, not bogus
     682            sage: A.condition(p=632)
     683            Traceback (most recent call last):
     684            ...
     685            ValueError: condition number integer values of 'p' must be -2, -1, 1 or 2, not 632
     686        """
     687        if not self.is_square():
     688            raise TypeError("matrix must be square, not %s x %s" % (self.nrows(), self.ncols()))
     689        global numpy
     690        if numpy is None:
     691            import numpy
     692        import sage.rings.infinity
     693        import sage.rings.integer
     694        import sage.rings.real_double
     695        if p == sage.rings.infinity.Infinity:
     696            p = numpy.inf
     697        elif p == -sage.rings.infinity.Infinity:
     698            p = -numpy.inf
     699        elif p == 'frob':
     700            p = 'fro'
     701        else:
     702            try:
     703                p = sage.rings.integer.Integer(p)
     704            except:
     705                raise ValueError("condition number 'p' must be +/- infinity, 'frob' or an integer, not %s" % p)
     706            if not p in [-2,-1,1,2]:
     707                raise ValueError("condition number integer values of 'p' must be -2, -1, 1 or 2, not %s" % p)
     708        # may raise a LinAlgError if matrix is singular
     709        c = numpy.linalg.cond(self._matrix_numpy, p=p)
     710        if c == numpy.inf:
     711            return sage.rings.infinity.Infinity
     712        else:
     713            return sage.rings.real_double.RDF(c)
     714
     715    def norm(self, p='frob'):
     716        r"""
     717        Returns the norm of the matrix.
     718
     719        INPUT:
     720
     721        - ``p`` - default: 'frob' - controls which norm is computed,
     722          allowable values are 'frob' (for the Frobenius norm),
     723          integers -2, -1, 1, 2, positive and negative infinity.
     724          See output discussion for specifics.
     725
     726        OUTPUT:
     727
     728        Returned value is a double precision floating point value
     729        in ``RDF``.  Row and column sums described below are
     730        sums of the absolute values of the entries, where the
     731        absolute value of the complex number `a+bi` is `\sqrt{a^2+b^2}`.
     732        Singular values are the "diagonal" entries of the "S" matrix in
     733        the singular value decomposition.
     734
     735        - ``p = 'frob'``: the default, the Frobenius norm, which for
     736          a matrix `A=(a_{ij})` computes
     737
     738          .. math::
     739
     740                \left(\sum_{i,j}\left\lvert{a_{i,j}}\right\rvert^2\right)^{1/2}
     741
     742        - ``p = Infinity`` or ``p = oo``: the maximum row sum.
     743        - ``p = -Infinity`` or ``p = -oo``: the minimum colum sum.
     744        - ``p = 1``: the maximum column sum.
     745        - ``p = -1``: the minimum column sum.
     746        - ``p = 2``: the 2-norm, equal to the maximum singular value.
     747        - ``p = -2``: the minimum singular value.
     748
     749        ALGORITHM:
     750
     751        Computation is performed by the ``norm()`` function of
     752        the SciPy/NumPy library.
     753
     754        EXAMPLES:
     755
     756        First over the reals.  ::
     757
     758            sage: A = matrix(RDF, 3, range(-3, 6)); A
     759            [-3.0 -2.0 -1.0]
     760            [ 0.0  1.0  2.0]
     761            [ 3.0  4.0  5.0]
     762            sage: A.norm()
     763            8.30662386...
     764            sage: A.norm(p='frob')
     765            8.30662386...
     766            sage: A.norm(p=Infinity)
     767            12.0
     768            sage: A.norm(p=-Infinity)
     769            3.0
     770            sage: A.norm(p=1)
     771            8.0
     772            sage: A.norm(p=-1)
     773            6.0
     774            sage: A.norm(p=2)
     775            7.99575670...
     776            sage: A.norm(p=-2)
     777            3.84592537...e-16
     778
     779        And over the complex numbers.  ::
     780
     781            sage: B = matrix(CDF, 2, [[1+I, 2+3*I],[3+4*I,3*I]]); B
     782            [1.0 + 1.0*I 2.0 + 3.0*I]
     783            [3.0 + 4.0*I       3.0*I]
     784            sage: B.norm()
     785            7.0
     786            sage: B.norm(p='frob')
     787            7.0
     788            sage: B.norm(p=Infinity)
     789            8.0
     790            sage: B.norm(p=-Infinity)
     791            5.01976483...
     792            sage: B.norm(p=1)
     793            6.60555127...
     794            sage: B.norm(p=-1)
     795            6.41421356...
     796            sage: B.norm(p=2)
     797            6.66189877...
     798            sage: B.norm(p=-2)
     799            2.14921023...
     800
     801        Since it is invariant under unitary multiplication, the
     802        Frobenius norm is equal to the square root of the sum of
     803        squares of the singular values.  ::
     804
     805            sage: A = matrix(RDF, 5, range(1,26))
     806            sage: f = A.norm(p='frob')
     807            sage: U, S, V = A.SVD()
     808            sage: s = sqrt(sum([S[i,i]^2 for i in range(5)]))
     809            sage: abs(f-s) < 1.0e-12
     810            True
     811
     812        Return values are in `RDF`.
     813
     814            sage: A = matrix(CDF, 2, range(4))
     815            sage: A.norm() in RDF
     816            True
     817
     818        Improper values of ``p`` are caught.  ::
     819
     820            sage: A.norm(p='bogus')
     821            Traceback (most recent call last):
     822            ...
     823            ValueError: matrix norm 'p' must be +/- infinity, 'frob' or an integer, not bogus
     824            sage: A.norm(p=632)
     825            Traceback (most recent call last):
     826            ...
     827            ValueError: matrix norm integer values of 'p' must be -2, -1, 1 or 2, not 632
     828        """
     829        global numpy
     830        if numpy is None:
     831            import numpy
     832        import sage.rings.infinity
     833        import sage.rings.integer
     834        import sage.rings.real_double
     835        if p == sage.rings.infinity.Infinity:
     836            p = numpy.inf
     837        elif p == -sage.rings.infinity.Infinity:
     838            p = -numpy.inf
     839        elif p == 'frob':
     840            p = 'fro'
     841        else:
     842            try:
     843                p = sage.rings.integer.Integer(p)
     844            except:
     845                raise ValueError("matrix norm 'p' must be +/- infinity, 'frob' or an integer, not %s" % p)
     846            if not p in [-2,-1,1,2]:
     847                raise ValueError("matrix norm integer values of 'p' must be -2, -1, 1 or 2, not %s" % p)
     848        return sage.rings.real_double.RDF(numpy.linalg.norm(self._matrix_numpy, ord=p))
     849
    541850    def LU(self):
    542851        """
    543852        Computes the LU decomposition of a matrix.

  • sage/modules/vector_double_dense.pyx 

    diff -r 3c67675a1df9 -r 5c31405d7399 sage/modules/vector_double_dense.pyx
    a b  
    3939
    4040from sage.structure.element cimport Element, ModuleElement, RingElement, Vector
    4141
     42from sage.rings.real_double import RDF
     43
    4244from sage.rings.complex_double import CDF
    4345from sage.rings.complex_double cimport ComplexDoubleElement, new_ComplexDoubleElement
    4446
    … …  
    571573        return self.change_ring(CDF)
    572574
    573575
     576    def norm(self, p=2):
     577        r"""
     578        Returns the norm (or related computations) of the vector.
     579
     580        INPUT:
     581
     582        - ``p`` - default: 2 - controls which norm is computed,
     583          allowable values are any real number and positive and
     584          negative infinity.  See output discussion for specifics.
     585
     586        OUTPUT:
     587
     588        Returned values is a double precision floating point value
     589        in ``RDF`` (or an integer when ``p=0``.  The default value
     590        of ``p = 2`` is the "usual" Euclidean norm.  For other values:
     591
     592        - ``p = Infinity`` or ``p = oo``: the maximum of the
     593          absolute values of the entries, where the absolute value
     594          of the complex number `a+bi` is `\sqrt{a^2+b^2}`.
     595        - ``p = -Infinity`` or ``p = -oo``: the minimum of the
     596          absolute values of the entries.
     597        - ``p = 0`` : the number of nonzero entries in the vector.
     598        - ``p`` is any other real number: for a vector `\vec{x}`
     599          this method computes
     600
     601          .. math::
     602
     603                \left(\sum_i x_i^p\right)^{1/p}
     604
     605          For ``p < 0`` this function is not a norm, but the above
     606          computation may be useful for other purposes.
     607
     608        ALGORITHM:
     609
     610        Computation is performed by the ``norm()`` function of
     611        the SciPy/NumPy library.
     612
     613        EXAMPLES:
     614
     615        First over the reals.  ::
     616
     617            sage: v = vector(RDF, range(9))
     618            sage: v.norm()
     619            14.28285685...
     620            sage: v.norm(p=2)
     621            14.28285685...
     622            sage: v.norm(p=6)
     623            8.744039097...
     624            sage: v.norm(p=Infinity)
     625            8.0
     626            sage: v.norm(p=-oo)
     627            0.0
     628            sage: v.norm(p=0)
     629            8
     630            sage: v.norm(p=0.3)
     631            4099.153615...
     632
     633        And over the complex numbers.  ::
     634
     635            sage: w = vector(CDF, [3-4*I, 0, 5+12*I])
     636            sage: w.norm()
     637            13.9283882...
     638            sage: w.norm(p=2)
     639            13.9283882...
     640            sage: w.norm(p=0)
     641            2
     642            sage: w.norm(p=4.2)
     643            13.0555695...
     644            sage: w.norm(p=oo)
     645            13.0
     646
     647        Negative values of ``p`` are allowed and will
     648        provide the same computation as for positive values.
     649        A zero entry in the vector will raise a warning and return
     650        zero. ::
     651
     652            sage: v = vector(CDF, range(1,10))
     653            sage: v.norm(p=-3.2)
     654            0.953760808...
     655            sage: w = vector(CDF, [-1,0,1])
     656            sage: w.norm(p=-1.6)
     657            0.0
     658
     659        Return values are in ``RDF``, or an integer when ``p = 0``.  ::
     660
     661            sage: v = vector(RDF, [1,2,4,8])
     662            sage: v.norm() in RDF
     663            True
     664            sage: v.norm(p=0) in ZZ
     665            True
     666
     667        Improper values of ``p`` are caught.  ::
     668
     669            sage: w = vector(CDF, [-1,0,1])
     670            sage: w.norm(p='junk')
     671            Traceback (most recent call last):
     672            ...
     673            ValueError: vector norm 'p' must be +/- infinity or a real number, not junk
     674        """
     675        global numpy
     676        if numpy is None:
     677            import numpy
     678        import sage.rings.infinity
     679        import sage.rings.integer
     680        if p == sage.rings.infinity.Infinity:
     681            p = numpy.inf
     682        elif p == -sage.rings.infinity.Infinity:
     683            p = -numpy.inf
     684        else:
     685            try:
     686                p = RDF(p)
     687            except:
     688                raise ValueError("vector norm 'p' must be +/- infinity or a real number, not %s" % p)
     689        n = numpy.linalg.norm(self._vector_numpy, ord=p)
     690        # p = 0 returns integer *count* of non-zero entries
     691        if n.dtype == numpy.int64:
     692            return sage.rings.integer.Integer(n)
     693        else:
     694            return RDF(n)
     695
     696
    574697    #############################
    575698    # statistics
    576699    #############################