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

File trac_10837-norms-condition-CDF-v2.patch, 16.7 KB (added by Rob Beezer, 11 years ago)

Rebase

  • sage/matrix/matrix2.pyx 

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

diff --git a/sage/matrix/matrix2.pyx b/sage/matrix/matrix2.pyx
    a b  
    88158815            sage: A.norm()
    88168816            15.0
    88178817       
    8818         ::
    8819        
    8820             sage: A = matrix(CDF, 2, 3, [3*I,4,1-I,1,2,0])
     8818        Norms of numerical matrices over high-precision reals are computed by this routine.
     8819        Faster routines for double precision entries from `RDF` or `CDF` are provided by
     8820        the :class:`~sage.matrix.matrix_double_dense.Matrix_double_dense` class.  ::
     8821       
     8822            sage: A = matrix(RR, 2, 3, [3*I,4,1-I,1,2,0])
    88218823            sage: A.norm('frob')
    88228824            5.65685424949
    88238825            sage: A.norm(2)

  • sage/matrix/matrix_double_dense.pyx 

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

  • sage/modules/vector_double_dense.pyx 

    diff --git a/sage/modules/vector_double_dense.pyx b/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    #############################