X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fsymmetric_psd.py;h=d711773034666dfb8fedcfb1dabdd9bc34aa344e;hb=928b7d49fda98ff105c92293b5797bb7a2b9873a;hp=b3f05f62c707e5b6c7e02d8c2806e23f543f7d03;hpb=af8292334d767865c22bd2f1002ab1f214e9f1ae;p=sage.d.git diff --git a/mjo/cone/symmetric_psd.py b/mjo/cone/symmetric_psd.py index b3f05f6..d711773 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,84 +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('../../')) - - -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) @@ -92,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 ) @@ -145,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 @@ -194,15 +136,15 @@ 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_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 @@ -241,40 +183,44 @@ 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 random_symmetric_psd + EXAMPLES: Well, it doesn't crash at least:: 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) + sage: A.is_positive_semidefinite() True A matrix with the desired rank is returned:: 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: 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``. @@ -300,19 +246,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