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

diff --git a/mjo/cone/symmetric_psd.py b/mjo/cone/symmetric_psd.py
index 1b3dd8a..fd6f950 100644
--- a/mjo/cone/symmetric_psd.py
+++ b/mjo/cone/symmetric_psd.py
@@ -1,4 +1,4 @@
-"""
+r"""
 The positive semidefinite cone `$S^{n}_{+}$` is the cone consisting of
 all symmetric positive-semidefinite matrices (as a subset of
 `$\mathbb{R}^{n \times n}$`
@@ -29,6 +29,7 @@ def is_symmetric_psd(A):
     Every completely positive matrix is symmetric
     positive-semidefinite::
 
+        sage: set_random_seed()
         sage: v = vector(map(abs, random_vector(ZZ, 10)))
         sage: A = v.column() * v.row()
         sage: is_symmetric_psd(A)
@@ -46,6 +47,13 @@ def is_symmetric_psd(A):
         sage: is_symmetric_psd(A)
         False
 
+    The trivial matrix in a trivial space is trivially symmetric and
+    positive-semidefinite::
+
+        sage: A = matrix(QQ, 0,0)
+        sage: is_symmetric_psd(A)
+        True
+
     """
 
     if A.base_ring() == SR:
@@ -60,7 +68,11 @@ def is_symmetric_psd(A):
     # 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
+    if A.is_zero():
+        # A is trivial... so trivially positive-semudefinite.
+        return True
+    else:
+        return min(A.eigenvalues()) >= 0
 
 
 def unit_eigenvectors(A):
@@ -201,7 +213,7 @@ def factor_psd(A):
     return Q*root_D*Q.transpose()
 
 
-def random_psd(V, accept_zero=True, rank=None):
+def random_symmetric_psd(V, accept_zero=True, rank=None):
     """
     Generate a random symmetric positive-semidefinite matrix over the
     vector space ``V``. That is, the returned matrix will be a linear
@@ -242,14 +254,16 @@ def random_psd(V, accept_zero=True, rank=None):
 
     SETUP::
 
-        sage: from mjo.cone.symmetric_psd import is_symmetric_psd, random_psd
+        sage: from mjo.cone.symmetric_psd import (is_symmetric_psd,
+        ....:                                     random_symmetric_psd)
 
     EXAMPLES:
 
     Well, it doesn't crash at least::
 
+        sage: set_random_seed()
         sage: V = VectorSpace(QQ, 2)
-        sage: A = random_psd(V)
+        sage: A = random_symmetric_psd(V)
         sage: A.matrix_space()
         Full MatrixSpace of 2 by 2 dense matrices over Rational Field
         sage: is_symmetric_psd(A)
@@ -257,27 +271,29 @@ def random_psd(V, accept_zero=True, rank=None):
 
     A matrix with the desired rank is returned::
 
+        sage: set_random_seed()
         sage: V = VectorSpace(QQ, 5)
-        sage: A = random_psd(V,False,1)
+        sage: A = random_symmetric_psd(V,False,1)
         sage: A.rank()
         1
-        sage: A = random_psd(V,False,2)
+        sage: A = random_symmetric_psd(V,False,2)
         sage: A.rank()
         2
-        sage: A = random_psd(V,False,3)
+        sage: A = random_symmetric_psd(V,False,3)
         sage: A.rank()
         3
-        sage: A = random_psd(V,False,4)
+        sage: A = random_symmetric_psd(V,False,4)
         sage: A.rank()
         4
-        sage: A = random_psd(V,False,5)
+        sage: A = random_symmetric_psd(V,False,5)
         sage: A.rank()
         5
 
     If the user asks for a rank that's too high, we fail::
 
+        sage: set_random_seed()
         sage: V = VectorSpace(QQ, 2)
-        sage: A = random_psd(V,False,3)
+        sage: A = random_symmetric_psd(V,False,3)
         Traceback (most recent call last):
         ...
         ValueError: The ``rank`` must be between 0 and the dimension of ``V``.
@@ -303,19 +319,23 @@ def random_psd(V, accept_zero=True, rank=None):
         # Use the one the user gave us.
         rank_A = rank
 
-    # Begin with the zero matrix, and add projectors to it if we have any.
-    A = V.zero().column()*V.zero().row()
-
-    # Careful, begin at idx=1 so that we only generate a projector
-    # when rank_A is greater than zero.
-    while A.rank() < rank_A:
-        v = V.random_element()
-        A += v.column()*v.row()
-
-    if accept_zero or not A.is_zero():
-        # We either don't care what ``A`` is, or it's non-zero, so
-        # just return it.
-        return A
-    else:
-        # Uh oh, we need to generate a new one.
-        return random_psd(V, accept_zero, rank)
+    if n == 0 and not accept_zero:
+        # We're gonna loop forever trying to satisfy this...
+        raise ValueError('You must have accept_zero=True when V is trivial')
+
+    # Loop until we find a suitable "A" that will then be returned.
+    while True:
+        # Begin with the zero matrix, and add projectors to it if we
+        # have any.
+        A = matrix.zero(V.base_ring(), n, n)
+
+        # Careful, begin at idx=1 so that we only generate a projector
+        # when rank_A is greater than zero.
+        while A.rank() < rank_A:
+            v = V.random_element()
+            A += v.column()*v.row()
+
+        if accept_zero or not A.is_zero():
+            # We either don't care what ``A`` is, or it's non-zero, so
+            # just return it.
+            return A