X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fsymmetric_psd.py;h=c6b9f69dd179a9f9950cbcd3b96a0c6bfd4e8436;hb=1e8a8f3916b47dbd059928d0e7d684364fb817dd;hp=fae0ae2307bb4e8b0583cf5b2200951f66a7e256;hpb=bc32004c617931e80f5c87775086e3f0e05f502e;p=sage.d.git

diff --git a/mjo/cone/symmetric_psd.py b/mjo/cone/symmetric_psd.py
index fae0ae2..c6b9f69 100644
--- a/mjo/cone/symmetric_psd.py
+++ b/mjo/cone/symmetric_psd.py
@@ -15,13 +15,66 @@ 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.
 
     INPUT:
 
-      - ``A`` - The matrix whose eigenvectors we want to compute.
+    - ``A`` - The matrix whose eigenvectors we want to compute.
 
     OUTPUT:
 
@@ -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)