Refactor symmetric_pds/doubly_nonnegative.

[sage.d.git] / mjo / cone / symmetric_psd.py

diff --git a/mjo/cone/symmetric_psd.py b/mjo/cone/symmetric_psd.py
index fae0ae2307bb4e8b0583cf5b2200951f66a7e256..d85ac247e12ca2a4936851182081415bba688061 100644 (file)
--- 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)