From: Michael Orlitzky Date: Thu, 30 Oct 2014 18:04:33 +0000 (-0400) Subject: Refactor symmetric_pds/doubly_nonnegative. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;ds=sidebyside;h=2d95c4e34d085c9647c73b73b2957f937cfee26b;p=sage.d.git Refactor symmetric_pds/doubly_nonnegative. Ensure that we're working over an amenable base ring in factor_psd. New module: completely_positive (cone). --- diff --git a/mjo/cone/completely_positive.py b/mjo/cone/completely_positive.py new file mode 100644 index 0000000..3353c24 --- /dev/null +++ b/mjo/cone/completely_positive.py @@ -0,0 +1,68 @@ +""" +The completely positive cone `$\mathcal{K}$` over `\mathbb{R}^{n}$` is +the set of all matrices `$A$`of the form `$\sum uu^{T}$` for `$u \in +\mathbb{R}^{n}_{+}$`. Equivalently, `$A = XX{T}$` where all entries of +`$X$` are nonnegative. +""" + +# Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we +# have to explicitly mangle our sitedir here so that "mjo.cone" +# resolves. +from os.path import abspath +from site import addsitedir +addsitedir(abspath('../../')) +from mjo.cone.symmetric_psd import factor_psd, is_symmetric_psd + + +def is_completely_positive(A): + """ + Determine whether or not the matrix ``A`` is completely positive. + + INPUT: + + - ``A`` - The matrix in question + + OUTPUT: + + Either ``True`` if ``A`` is completely positive, or ``False`` + otherwise. + + EXAMPLES: + + Generate an extreme completely positive matrix and check that we + identify it correctly:: + + sage: v = vector(map(abs, random_vector(ZZ, 10))) + sage: A = v.column() * v.row() + sage: A = A.change_ring(QQ) + sage: is_completely_positive(A) + True + + The following matrix isn't positive semidefinite, so it can't be + completely positive:: + + sage: A = matrix(QQ, [[1, 2], [2, 1]]) + sage: is_completely_positive(A) + False + + This matrix isn't even symmetric, so it can't be completely + positive:: + + sage: A = matrix(QQ, [[1, 2], [3, 4]]) + sage: is_completely_positive(A) + False + + """ + + if A.base_ring() == SR: + msg = 'The matrix ``A`` cannot be the symbolic.' + raise ValueError.new(msg) + + if not is_symmetric_psd(A): + return False + + # Would crash if we hadn't ensured that ``A`` was symmetric + # positive-semidefinite. + X = factor_psd(A) + return X.rank() == 1 + diff --git a/mjo/cone/doubly_nonnegative.py b/mjo/cone/doubly_nonnegative.py index b071e41..f705e5e 100644 --- a/mjo/cone/doubly_nonnegative.py +++ b/mjo/cone/doubly_nonnegative.py @@ -19,7 +19,7 @@ from sage.all import * from os.path import abspath from site import addsitedir addsitedir(abspath('../../')) -from mjo.cone.symmetric_psd import factor_psd +from mjo.cone.symmetric_psd import factor_psd, is_symmetric_psd @@ -54,22 +54,16 @@ def is_doubly_nonnegative(A): """ if A.base_ring() == SR: - msg = 'The base ring of ``A`` cannot be the Symbolic Ring' + msg = 'The matrix ``A`` cannot be the symbolic.' raise ValueError.new(msg) - # First make sure that ``A`` is symmetric. - if not A.is_symmetric(): - return False - # Check that all of the entries of ``A`` are nonnegative. if not all([ a >= 0 for a in A.list() ]): return False - # If ``A`` is symmetric and non-negative, we only need to check - # that it is positive semidefinite. For that we can consult its - # minimum eigenvalue, which should be zero or greater. Since ``A`` - # is symmetric, its eigenvalues are guaranteed to be real. - return min(A.eigenvalues()) >= 0 + # It's nonnegative, so all we need to do is check that it's + # symmetric positive-semidefinite. + return is_symmetric_psd(A) diff --git a/mjo/cone/symmetric_psd.py b/mjo/cone/symmetric_psd.py index fae0ae2..d85ac24 100644 --- a/mjo/cone/symmetric_psd.py +++ b/mjo/cone/symmetric_psd.py @@ -15,6 +15,59 @@ addsitedir(abspath('../../')) from mjo.symbolic import matrix_simplify_full +def is_symmetric_psd(A): + """ + Determine whether or not the matrix ``A`` is symmetric + positive-semidefinite. + + INPUT: + + - ``A`` - The matrix in question + + OUTPUT: + + Either ``True`` if ``A`` is symmetric positive-semidefinite, or + ``False`` otherwise. + + EXAMPLES: + + Every completely positive matrix is symmetric + positive-semidefinite:: + + sage: v = vector(map(abs, random_vector(ZZ, 10))) + sage: A = v.column() * v.row() + sage: is_symmetric_psd(A) + True + + The following matrix is symmetric but not positive semidefinite:: + + sage: A = matrix(ZZ, [[1, 2], [2, 1]]) + sage: is_symmetric_psd(A) + False + + This matrix isn't even symmetric:: + + sage: A = matrix(ZZ, [[1, 2], [3, 4]]) + sage: is_symmetric_psd(A) + False + + """ + + if A.base_ring() == SR: + msg = 'The matrix ``A`` cannot be symbolic.' + raise ValueError.new(msg) + + # First make sure that ``A`` is symmetric. + if not A.is_symmetric(): + return False + + # If ``A`` is symmetric, we only need to check that it is positive + # semidefinite. For that we can consult its minimum eigenvalue, + # which should be zero or greater. Since ``A`` is symmetric, its + # eigenvalues are guaranteed to be real. + return min(A.eigenvalues()) >= 0 + + def unit_eigenvectors(A): """ Return the unit eigenvectors of a symmetric positive-definite matrix. @@ -61,11 +114,15 @@ def factor_psd(A): INPUT: - - ``A`` - The matrix to factor. + - ``A`` - The matrix to factor. The base ring of ``A`` must be either + exact or the symbolic ring (to compute eigenvalues), and it + must be a field so that we can take its algebraic closure + (necessary to e.g. take square roots). OUTPUT: - A matrix ``X`` such that `A = XX^{T}`. + A matrix ``X`` such that `A = XX^{T}`. The base field of ``X`` will + be the algebraic closure of the base field of ``A``. ALGORITHM: @@ -89,7 +146,10 @@ def factor_psd(A): `$D$` will have dimension `$k \times k$`. In the end everything works out the same. - EXAMPLES:: + EXAMPLES: + + Create a symmetric positive-semidefinite matrix over the symbolic + ring and factor it:: sage: A = matrix(SR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]]) sage: X = factor_psd(A) @@ -97,8 +157,38 @@ def factor_psd(A): sage: A == A2 True + Attempt to factor the same matrix over ``RR`` which won't work + because ``RR`` isn't exact:: + + sage: A = matrix(RR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]]) + sage: factor_psd(A) + Traceback (most recent call last): + ... + ValueError: The base ring of ``A`` must be either exact or symbolic. + + Attempt to factor the same matrix over ``ZZ`` which won't work + because ``ZZ`` isn't a field:: + + sage: A = matrix(ZZ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]]) + sage: factor_psd(A) + Traceback (most recent call last): + ... + ValueError: The base ring of ``A`` must be a field. + """ + if not A.base_ring().is_exact() and not A.base_ring() is SR: + raise ValueError('The base ring of ``A`` must be either exact or symbolic.') + + if not A.base_ring().is_field(): + raise ValueError('The base ring of ``A`` must be a field.') + + if not A.base_ring() is SR: + # Change the base field of ``A`` so that we are sure we can take + # roots. The symbolic ring has no algebraic closure. + A = A.change_ring(A.base_ring().algebraic_closure()) + + # Get the eigenvectors, and filter out the ones that correspond to # the eigenvalue zero. all_evs = unit_eigenvectors(A)