From 6bb77d76f9a639a94b836fc24f067a37a6186f42 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Tue, 4 Nov 2014 09:57:53 -0500
Subject: [PATCH] Add the random_psd() function.

---
 mjo/cone/symmetric_psd.py | 117 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 117 insertions(+)

diff --git a/mjo/cone/symmetric_psd.py b/mjo/cone/symmetric_psd.py
index c9f2a4d..bb46fc9 100644
--- a/mjo/cone/symmetric_psd.py
+++ b/mjo/cone/symmetric_psd.py
@@ -201,3 +201,120 @@ def factor_psd(A):
     Q = matrix(A.base_ring(), [ vec for (val,vec) in evs ]).transpose()
 
     return Q*root_D*Q.transpose()
+
+
+def random_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
+    transformation on ``V``, with the same base ring as ``V``.
+
+    We take a very loose interpretation of "random," here. Otherwise we
+    would never (for example) choose a matrix on the boundary of the
+    cone (with a zero eigenvalue).
+
+    INPUT:
+
+    - ``V`` - The vector space on which the returned matrix will act.
+
+    - ``accept_zero`` - Do you want to accept the zero matrix (which
+                        is symmetric PSD? Defaults to ``True``.
+
+    - ``rank`` - Require the returned matrix to have the given rank
+                 (optional).
+
+    OUTPUT:
+
+    A random symmetric positive semidefinite matrix, i.e. a linear
+    transformation from ``V`` to itself.
+
+    ALGORITHM:
+
+    The matrix is constructed from some number of spectral projectors,
+    which in turn are created at "random" from the underlying vector
+    space ``V``.
+
+    If no particular ``rank`` is desired, we choose the number of
+    projectors at random. Otherwise, we keep adding new projectors until
+    the desired rank is achieved.
+
+    Finally, before returning, we check if the matrix is zero. If
+    ``accept_zero`` is ``False``, we restart the process from the
+    beginning.
+
+    EXAMPLES:
+
+    Well, it doesn't crash at least::
+
+        sage: V = VectorSpace(QQ, 2)
+        sage: A = random_psd(V)
+        sage: A.matrix_space()
+        Full MatrixSpace of 2 by 2 dense matrices over Rational Field
+        sage: is_symmetric_psd(A)
+        True
+
+    A matrix with the desired rank is returned::
+
+        sage: V = VectorSpace(QQ, 5)
+        sage: A = random_psd(V,False,1)
+        sage: A.rank()
+        1
+        sage: A = random_psd(V,False,2)
+        sage: A.rank()
+        2
+        sage: A = random_psd(V,False,3)
+        sage: A.rank()
+        3
+        sage: A = random_psd(V,False,4)
+        sage: A.rank()
+        4
+        sage: A = random_psd(V,False,5)
+        sage: A.rank()
+        5
+
+    If the user asks for a rank that's too high, we fail::
+
+        sage: V = VectorSpace(QQ, 2)
+        sage: A = random_psd(V,False,3)
+        Traceback (most recent call last):
+        ...
+        ValueError: The ``rank`` must be between 0 and the dimension of ``V``.
+
+    """
+
+    # We construct the matrix from its spectral projectors. Since
+    # there can be at most ``n`` of them, where ``n`` is the dimension
+    # of our vector space, we want to choose a random integer between
+    # ``0`` and ``n`` and then construct that many random elements of
+    # ``V``.
+    n = V.dimension()
+
+    rank_A = 0
+    if rank is None:
+        # Choose one randomly
+        rank_A = ZZ.random_element(n+1)
+    elif (rank < 0) or (rank > n):
+        # The rank of ``A`` can be at most ``n``.
+        msg  = 'The ``rank`` must be between 0 and the dimension of ``V``.'
+        raise ValueError(msg)
+    else:
+        # 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_element().column()*V.zero_element().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)
-- 
2.43.2