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

diff --git a/mjo/cone/symmetric_psd.py b/mjo/cone/symmetric_psd.py
index e4be629..1b3dd8a 100644
--- a/mjo/cone/symmetric_psd.py
+++ b/mjo/cone/symmetric_psd.py
@@ -20,6 +20,10 @@ def is_symmetric_psd(A):
     Either ``True`` if ``A`` is symmetric positive-semidefinite, or
     ``False`` otherwise.
 
+    SETUP::
+
+        sage: from mjo.cone.symmetric_psd import is_symmetric_psd
+
     EXAMPLES:
 
     Every completely positive matrix is symmetric
@@ -65,17 +69,24 @@ def unit_eigenvectors(A):
 
     INPUT:
 
-    - ``A`` - The matrix whose eigenvectors we want to compute.
+    - ``A`` -- The matrix whose unit eigenvectors we want to compute.
 
     OUTPUT:
 
     A list of (eigenvalue, eigenvector) pairs where each eigenvector is
-    associated with its paired eigenvalue of ``A`` and has norm `1`.
+    associated with its paired eigenvalue of ``A`` and has norm `1`. If
+    the base ring of ``A`` is not algebraically closed, then returned
+    eigenvectors may (necessarily) be over its algebraic closure and not
+    the base ring of ``A`` itself.
+
+    SETUP::
+
+        sage: from mjo.cone.symmetric_psd import unit_eigenvectors
 
     EXAMPLES::
 
         sage: A = matrix(QQ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
-        sage: unit_evs = unit_eigenvectors(A)
+        sage: unit_evs = list(unit_eigenvectors(A))
         sage: bool(unit_evs[0][1].norm() == 1)
         True
         sage: bool(unit_evs[1][1].norm() == 1)
@@ -84,17 +95,9 @@ def unit_eigenvectors(A):
         True
 
     """
-    # This will give us a list of lists whose elements are the
-    # eigenvectors we want.
-    ev_lists = [ (val,vecs) for (val,vecs,multiplicity)
-                            in A.eigenvectors_right() ]
-
-    # Pair each eigenvector with its eigenvalue and normalize it.
-    evs = [ [(l, vec/vec.norm()) for vec in vecs] for (l,vecs) in ev_lists ]
-
-    # Flatten the list, abusing the fact that "+" is overloaded on lists.
-    return sum(evs, [])
-
+    return ( (val,vec.normalized())
+             for (val,vecs,multiplicity) in A.eigenvectors_right()
+             for vec in vecs )
 
 
 
@@ -137,6 +140,10 @@ def factor_psd(A):
     `$D$` will have dimension `$k \times k$`. In the end everything
     works out the same.
 
+    SETUP::
+
+        sage: from mjo.cone.symmetric_psd import factor_psd
+
     EXAMPLES:
 
     Create a symmetric positive-semidefinite matrix over the symbolic
@@ -186,10 +193,10 @@ def factor_psd(A):
     all_evs = unit_eigenvectors(A)
     evs = [ (val,vec) for (val,vec) in all_evs if not val == 0 ]
 
-    d = [ sqrt(val) for (val,vec) in evs ]
+    d = ( val.sqrt() for (val,vec) in evs )
     root_D = diagonal_matrix(d).change_ring(A.base_ring())
 
-    Q = matrix(A.base_ring(), [ vec for (val,vec) in evs ]).transpose()
+    Q = matrix(A.base_ring(), ( vec for (val,vec) in evs )).transpose()
 
     return Q*root_D*Q.transpose()
 
@@ -233,6 +240,10 @@ def random_psd(V, accept_zero=True, rank=None):
     ``accept_zero`` is ``False``, we restart the process from the
     beginning.
 
+    SETUP::
+
+        sage: from mjo.cone.symmetric_psd import is_symmetric_psd, random_psd
+
     EXAMPLES:
 
     Well, it doesn't crash at least::