Ticket #13018: trac_13018-is-positive-definite-v2.patch

File trac_13018-is-positive-definite-v2.patch, 9.0 KB (added by rbeezer, 10 years ago)
  • sage/matrix/matrix2.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1338050285 25200
    # Node ID 479229c6ef975b07734c626612766cc6a46e242b
    # Parent  548a894e33b20ed79fdd778dfbf386b2c7748367
    Positive definite check for exact matrices
    
    diff --git a/sage/matrix/matrix2.pyx b/sage/matrix/matrix2.pyx
    a b  
    1034610346        from sage.modules.free_module_element import vector
    1034710347        L, d = self._indefinite_factorization(algorithm, check=check)
    1034810348        return L, vector(L.base_ring(), d)
     10349
     10350    def is_positive_definite(self):
     10351        r"""
     10352        Determines if a real or symmetric matrix is positive definite.
     10353
     10354        A square matrix `A` is postive definite if it is
     10355        symmetric with real entries or Hermitan with complex entries,
     10356        and for every non-zero vector `\vec{x}`
     10357
     10358        .. math::
     10359
     10360            \vec{x}^\ast A\vec{x} > 0
     10361
     10362        Here `\vec{x}^\ast` is the conjugate-transpose, which can be
     10363        simplified to just the transpose in the real case.
     10364
     10365        ALGORITHM:
     10366
     10367        A matrix is positive definite if and only if the
     10368        diagonal entries from the indefinite factorization
     10369        are all positive (see :meth:`indefinite_factorization`).
     10370        So this algorithm is of order ``n^3/3`` and may be applied
     10371        to matrices with elements of any ring that has a fraction
     10372        field contained within the reals or complexes.
     10373       
     10374        INPUT:
     10375
     10376        Any square matrix.
     10377
     10378        OUTPUT:
     10379
     10380        This routine will return ``True`` if the matrix is square,
     10381        symmetric or Hermitian, and meets the condition above
     10382        for the quadratic form.
     10383
     10384        The base ring for the elements of the matrix needs to
     10385        have a fraction field implemented and the computations
     10386        that result from the indefinite factorization must be
     10387        convertable to real numbers that are comparable to zero.
     10388
     10389        EXAMPLES:
     10390
     10391        A real symmetric matrix that is positive definite,
     10392        as evidenced by the positive entries for the diagonal
     10393        matrix of the indefinite factorization and the postive
     10394        determinants of the leading principal submatrices. ::
     10395
     10396            sage: A = matrix(QQ, [[ 4, -2,  4,  2],
     10397            ...                   [-2, 10, -2, -7],
     10398            ...                   [ 4, -2,  8,  4],
     10399            ...                   [ 2, -7,  4,  7]])
     10400            sage: A.is_positive_definite()
     10401            True
     10402            sage: _, d = A.indefinite_factorization(algorithm='symmetric')
     10403            sage: d
     10404            (4, 9, 4, 1)
     10405            sage: [A[:i,:i].determinant() for i in range(1,A.nrows()+1)]
     10406            [4, 36, 144, 144]
     10407
     10408        A real symmetric matrix which is not positive definite, along
     10409        with a vector that makes the quadratic form negative. ::
     10410
     10411            sage: A = matrix(QQ, [[ 3,  -6,   9,   6,  -9],
     10412            ...                   [-6,  11, -16, -11,  17],
     10413            ...                   [ 9, -16,  28,  16, -40],
     10414            ...                   [ 6, -11,  16,   9, -19],
     10415            ...                   [-9,  17, -40, -19,  68]])
     10416            sage: A.is_positive_definite()
     10417            False
     10418            sage: _, d = A.indefinite_factorization(algorithm='symmetric')
     10419            sage: d
     10420            (3, -1, 5, -2, -1)
     10421            sage: [A[:i,:i].determinant() for i in range(1,A.nrows()+1)]
     10422            [3, -3, -15, 30, -30]
     10423            sage: u = vector(QQ, [2, 2, 0, 1, 0])
     10424            sage: u.row()*A*u
     10425            (-3)
     10426           
     10427        A real symmetric matrix with a singular leading
     10428        principal submatrix, that is therefore not positive definite.
     10429        The vector ``u`` makes the quadratic form zero.  ::
     10430
     10431            sage: A = matrix(QQ, [[21, 15, 12, -2],
     10432            ...                   [15, 12,  9,  6],
     10433            ...                   [12,  9,  7,  3],
     10434            ...                   [-2,  6,  3,  8]])
     10435            sage: A.is_positive_definite()
     10436            False
     10437            sage: [A[:i,:i].determinant() for i in range(1,A.nrows()+1)]
     10438            [21, 27, 0, -75]
     10439            sage: u = vector(QQ, [1,1,-3,0])
     10440            sage: u.row()*A*u
     10441            (0)
     10442           
     10443        An Hermitian matrix that is positive definite.  ::
     10444
     10445            sage: C.<I> = NumberField(x^2 + 1)
     10446            sage: A = matrix(C, [[        23,  17*I + 3,  24*I + 25,     21*I],
     10447            ...                  [ -17*I + 3,        38, -69*I + 89, 7*I + 15],
     10448            ...                  [-24*I + 25, 69*I + 89,        976, 24*I + 6],
     10449            ...                  [     -21*I, -7*I + 15,  -24*I + 6,       28]])
     10450            sage: A.is_positive_definite()
     10451            True
     10452            sage: _, d = A.indefinite_factorization(algorithm='hermitian')
     10453            sage: d
     10454            (23, 576/23, 89885/144, 142130/17977)
     10455            sage: [A[:i,:i].determinant() for i in range(1,A.nrows()+1)]
     10456            [23, 576, 359540, 2842600]
     10457           
     10458        An Hermitian matrix that is not positive definite.
     10459        The vector ``u`` makes the quadratic form negative.  ::
     10460       
     10461            sage: C.<I> = QuadraticField(-1)
     10462            sage: B = matrix(C, [[      2, 4 - 2*I, 2 + 2*I],
     10463            ...                  [4 + 2*I,       8,    10*I],
     10464            ...                  [2 - 2*I,   -10*I,      -3]])
     10465            sage: B.is_positive_definite()
     10466            False
     10467            sage: _, d = B.indefinite_factorization(algorithm='hermitian')
     10468            sage: d
     10469            (2, -2, 3)
     10470            sage: [B[:i,:i].determinant() for i in range(1,B.nrows()+1)]
     10471            [2, -4, -12]
     10472            sage: u = vector(C, [-5 + 10*I, 4 - 3*I, 0])
     10473            sage: u.row().conjugate()*B*u
     10474            (-50)
     10475           
     10476        A positive definite matrix over an algebraically closed field.  ::
     10477
     10478            sage: A = matrix(QQbar, [[        2,   4 + 2*I,   6 - 4*I],
     10479            ...                      [ -2*I + 4,        11, 10 - 12*I],
     10480            ...                      [  4*I + 6, 10 + 12*I,        37]])
     10481            sage: A.is_positive_definite()
     10482            True
     10483            sage: [A[:i,:i].determinant() for i in range(1,A.nrows()+1)]
     10484            [2, 2, 6]
     10485           
     10486        TESTS:
     10487
     10488        If the base ring lacks a ``conjugate`` method, it
     10489        will be assumed to not be Hermitian and thus symmetric.
     10490        If the base ring does not make sense as a subfield of
     10491        the reals, then this routine will fail since comparison
     10492        to zero is meaningless.  ::
     10493
     10494            sage: F.<a> = FiniteField(5^3)
     10495            sage: a.conjugate()
     10496            Traceback (most recent call last):
     10497            ...
     10498            AttributeError: 'sage.rings.finite_rings.element_givaro.FiniteField_givaroElement'
     10499            object has no attribute 'conjugate'
     10500            sage: A = matrix(F,
     10501            ...       [[      a^2 + 2*a, 4*a^2 + 3*a + 4,       3*a^2 + a, 2*a^2 + 2*a + 1],
     10502            ...        [4*a^2 + 3*a + 4,       4*a^2 + 2,             3*a, 2*a^2 + 4*a + 2],
     10503            ...        [      3*a^2 + a,             3*a,       3*a^2 + 2, 3*a^2 + 2*a + 3],
     10504            ...        [2*a^2 + 2*a + 1, 2*a^2 + 4*a + 2, 3*a^2 + 2*a + 3, 3*a^2 + 2*a + 4]])
     10505            sage: A.is_positive_definite()
     10506            Traceback (most recent call last):
     10507            ...
     10508            TypeError: cannot convert computations from
     10509            Finite Field in a of size 5^3 into real numbers
     10510       
     10511        AUTHOR:
     10512
     10513        - Rob Beezer (2012-05-24)
     10514        """
     10515        from sage.rings.all import RR
     10516        # check to see if the Hermitian routine is usable
     10517        # otherwise we will assume the symmetric case
     10518        imaginary = True
     10519        a = self.base_ring().an_element()
     10520        try:
     10521            a.conjugate()
     10522        except AttributeError:
     10523            imaginary = False
     10524        if imaginary:
     10525            if not self.is_hermitian():
     10526                return False
     10527        else:
     10528            if not self.is_symmetric():
     10529                return False
     10530        try:
     10531            if imaginary:
     10532                _, d = self._indefinite_factorization('hermitian', check=False)
     10533            else:
     10534                _, d = self._indefinite_factorization('symmetric', check=False)
     10535        except ValueError as e:
     10536            # a zero singular leading principal submatrix is one
     10537            # indicator that the matrix is not positive definite
     10538            if e.message.find('leading principal submatrix is singular') != -1:
     10539                return False
     10540            else:
     10541                raise ValueError(e)
     10542        # Now have diagonal entries (hopefully real) and so can
     10543        # test with a generator (which will short-circuit)
     10544        # positive definite iff all entries of d are positive
     10545        try:
     10546            posdef = all( RR(x) > 0 for x in d )
     10547        except TypeError:
     10548            universe = Sequence(d).universe()
     10549            msg = "cannot convert computations from {0} into real numbers"
     10550            raise TypeError(msg.format(universe))
     10551        return posdef
    1034910552       
    1035010553    def hadamard_bound(self):
    1035110554        r"""