X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fsymmetric_psd.py;h=89eba53c9687823a0882304a49d9ddefc1d009db;hb=7af2b9d146a6bf2fb8acc3c342983de577b417ce;hp=d85ac247e12ca2a4936851182081415bba688061;hpb=2d95c4e34d085c9647c73b73b2957f937cfee26b;p=sage.d.git diff --git a/mjo/cone/symmetric_psd.py b/mjo/cone/symmetric_psd.py index d85ac24..89eba53 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}$` @@ -6,85 +6,30 @@ 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 - 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. 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) @@ -93,17 +38,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 ) @@ -146,6 +83,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 +94,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 @@ -178,14 +119,15 @@ def factor_psd(A): """ 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.') + msg = 'The base ring of ``A`` must be either exact or symbolic.' + raise ValueError(msg) 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. + # roots. The symbolic ring has no algebraic_closure method. A = A.change_ring(A.base_ring().algebraic_closure()) @@ -194,9 +136,136 @@ 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() + + +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 + 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 random_symmetric_psd + + EXAMPLES: + + Well, it doesn't crash at least:: + + sage: set_random_seed() + sage: V = VectorSpace(QQ, 2) + sage: A = random_symmetric_psd(V) + sage: A.matrix_space() + Full MatrixSpace of 2 by 2 dense matrices over Rational Field + sage: A.is_positive_semidefinite() + True + + A matrix with the desired rank is returned:: + + sage: set_random_seed() + sage: V = VectorSpace(QQ, 5) + sage: A = random_symmetric_psd(V,False,1) + sage: A.rank() + 1 + sage: A = random_symmetric_psd(V,False,2) + sage: A.rank() + 2 + sage: A = random_symmetric_psd(V,False,3) + sage: A.rank() + 3 + sage: A = random_symmetric_psd(V,False,4) + sage: A.rank() + 4 + 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_symmetric_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 + + 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