X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fsymmetric_psd.py;h=ef5d477ce67aa4f40a3466c48904888ddfaa840b;hb=14b952fe32ba20b7bac685da07633b2ce5b5a5af;hp=c9f2a4df4ac87ffca9e44a956516af794a934ee7;hpb=68df5e61c6a34482ec8cb3bd24a289593101d5e0;p=sage.d.git diff --git a/mjo/cone/symmetric_psd.py b/mjo/cone/symmetric_psd.py index c9f2a4d..ef5d477 100644 --- a/mjo/cone/symmetric_psd.py +++ b/mjo/cone/symmetric_psd.py @@ -6,15 +6,6 @@ all symmetric positive-semidefinite matrices (as a subset of from sage.all import * -# Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we -# have to explicitly mangle our sitedir here so that "mjo.symbolic" -# resolves. -from os.path import abspath -from site import addsitedir -addsitedir(abspath('../../')) -from mjo.symbolic import matrix_simplify_full - - def is_symmetric_psd(A): """ Determine whether or not the matrix ``A`` is symmetric @@ -29,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 @@ -81,6 +76,10 @@ def unit_eigenvectors(A): A list of (eigenvalue, eigenvector) pairs where each eigenvector is associated with its paired eigenvalue of ``A`` and has norm `1`. + SETUP:: + + sage: from mjo.cone.symmetric_psd import unit_eigenvectors + EXAMPLES:: sage: A = matrix(QQ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]]) @@ -146,6 +145,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 @@ -153,7 +156,7 @@ def factor_psd(A): sage: A = matrix(SR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]]) sage: X = factor_psd(A) - sage: A2 = matrix_simplify_full(X*X.transpose()) + sage: A2 = (X*X.transpose()).simplify_full() sage: A == A2 True @@ -201,3 +204,123 @@ 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. + + SETUP:: + + sage: from mjo.cone.symmetric_psd import is_symmetric_psd, random_psd + + 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().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)