X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fdoubly_nonnegative.py;h=56a655e52c5cd60cbf1073c351e65833e5f54500;hb=c2d15c0884c8b92483f58826747887bd2bcdcdeb;hp=4f7f950bb2c897509c2cb4b41d5200bccfd67c49;hpb=a5700d65325514a505d24fb96b75c2b0f2f6e94f;p=sage.d.git diff --git a/mjo/cone/doubly_nonnegative.py b/mjo/cone/doubly_nonnegative.py index 4f7f950..56a655e 100644 --- a/mjo/cone/doubly_nonnegative.py +++ b/mjo/cone/doubly_nonnegative.py @@ -19,7 +19,7 @@ from sage.all import * from os.path import abspath from site import addsitedir addsitedir(abspath('../../')) -from mjo.cone.symmetric_psd import factor_psd, is_symmetric_psd +from mjo.cone.symmetric_psd import factor_psd, is_symmetric_psd, random_psd from mjo.matrix_vector import isomorphism @@ -67,6 +67,67 @@ def is_doubly_nonnegative(A): +def is_admissible_extreme_rank(r, n): + """ + The extreme matrices of the doubly-nonnegative cone have some + restrictions on their ranks. This function checks to see whether the + rank ``r`` would be an admissible rank for an ``n``-by-``n`` matrix. + + INPUT: + + - ``r`` - The rank of the matrix. + + - ``n`` - The dimension of the vector space on which the matrix acts. + + OUTPUT: + + Either ``True`` if a rank ``r`` matrix could be an extreme vector of + the doubly-nonnegative cone in `$\mathbb{R}^{n}$`, or ``False`` + otherwise. + + EXAMPLES: + + For dimension 5, only ranks zero, one, and three are admissible:: + + sage: is_admissible_extreme_rank(0,5) + True + sage: is_admissible_extreme_rank(1,5) + True + sage: is_admissible_extreme_rank(2,5) + False + sage: is_admissible_extreme_rank(3,5) + True + sage: is_admissible_extreme_rank(4,5) + False + sage: is_admissible_extreme_rank(5,5) + False + + When given an impossible rank, we just return false:: + + sage: is_admissible_extreme_rank(100,5) + False + + """ + if r == 0: + # Zero is in the doubly-nonnegative cone. + return True + + if r > n: + # Impossible, just return False + return False + + # See Theorem 3.1 in the cited reference. + if r == 2: + return False + + if n.mod(2) == 0: + # n is even + return r <= max(1, n-3) + else: + # n is odd + return r <= max(1, n-2) + + def has_admissible_extreme_rank(A): """ The extreme matrices of the doubly-nonnegative cone have some @@ -93,9 +154,35 @@ def has_admissible_extreme_rank(A): The zero matrix has rank zero, which is admissible:: - sage: A = zero_matrix(QQ, 5, 5) - sage: has_admissible_extreme_rank(A) - True + sage: A = zero_matrix(QQ, 5, 5) + sage: has_admissible_extreme_rank(A) + True + + Likewise, rank one is admissible for dimension 5:: + + sage: v = vector(QQ, [1,2,3,4,5]) + sage: A = v.column()*v.row() + sage: has_admissible_extreme_rank(A) + True + + But rank 2 is never admissible:: + + sage: v1 = vector(QQ, [1,0,0,0,0]) + sage: v2 = vector(QQ, [0,1,0,0,0]) + sage: A = v1.column()*v1.row() + v2.column()*v2.row() + sage: has_admissible_extreme_rank(A) + False + + In dimension 5, three is the only other admissible rank:: + + sage: v1 = vector(QQ, [1,0,0,0,0]) + sage: v2 = vector(QQ, [0,1,0,0,0]) + sage: v3 = vector(QQ, [0,0,1,0,0]) + sage: A = v1.column()*v1.row() + sage: A += v2.column()*v2.row() + sage: A += v3.column()*v3.row() + sage: has_admissible_extreme_rank(A) + True """ if not A.is_symmetric(): @@ -106,20 +193,7 @@ def has_admissible_extreme_rank(A): r = rank(A) n = ZZ(A.nrows()) # Columns would work, too, since ``A`` is symmetric. - if r == 0: - # Zero is in the doubly-nonnegative cone. - return True - - # See Theorem 3.1 in the cited reference. - if r == 2: - return False - - if n.mod(2) == 0: - # n is even - return r <= max(1, n-3) - else: - # n is odd - return r <= max(1, n-2) + return is_admissible_extreme_rank(r,n) def E(matrix_space, i,j): @@ -304,3 +378,141 @@ def is_extreme_doubly_nonnegative(A): # earlier. two = A.base_ring()(2) return d == (k*(k + 1)/two - 1) + + +def random_doubly_nonnegative(V, accept_zero=True, rank=None): + """ + Generate a random doubly nonnegative 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. + + INPUT: + + - ``V`` - The vector space on which the returned matrix will act. + + - ``accept_zero`` - Do you want to accept the zero matrix (which + is doubly nonnegative)? Default to ``True``. + + - ``rank`` - Require the returned matrix to have the given rank + (optional). + + OUTPUT: + + A random doubly nonnegative matrix, i.e. a linear transformation + from ``V`` to itself. + + EXAMPLES: + + Well, it doesn't crash at least:: + + sage: V = VectorSpace(QQ, 2) + sage: A = random_doubly_nonnegative(V) + sage: A.matrix_space() + Full MatrixSpace of 2 by 2 dense matrices over Rational Field + sage: is_doubly_nonnegative(A) + True + + A matrix with the desired rank is returned:: + + sage: V = VectorSpace(QQ, 5) + sage: A = random_doubly_nonnegative(V,False,1) + sage: A.rank() + 1 + sage: A = random_doubly_nonnegative(V,False,2) + sage: A.rank() + 2 + sage: A = random_doubly_nonnegative(V,False,3) + sage: A.rank() + 3 + sage: A = random_doubly_nonnegative(V,False,4) + sage: A.rank() + 4 + sage: A = random_doubly_nonnegative(V,False,5) + sage: A.rank() + 5 + + """ + + # Generate random symmetric positive-semidefinite matrices until + # one of them is nonnegative, then return that. + A = random_psd(V, accept_zero, rank) + + while not all([ x >= 0 for x in A.list() ]): + A = random_psd(V, accept_zero, rank) + + return A + + + +def random_extreme_doubly_nonnegative(V, accept_zero=True, rank=None): + """ + Generate a random extreme doubly nonnegative 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. + + INPUT: + + - ``V`` - The vector space on which the returned matrix will act. + + - ``accept_zero`` - Do you want to accept the zero matrix + (which is extreme)? Defaults to ``True``. + + - ``rank`` - Require the returned matrix to have the given rank + (optional). WARNING: certain ranks are not possible + in any given dimension! If an impossible rank is + requested, a ValueError will be raised. + + OUTPUT: + + A random extreme doubly nonnegative matrix, i.e. a linear + transformation from ``V`` to itself. + + EXAMPLES: + + Well, it doesn't crash at least:: + + sage: V = VectorSpace(QQ, 2) + sage: A = random_extreme_doubly_nonnegative(V) + sage: A.matrix_space() + Full MatrixSpace of 2 by 2 dense matrices over Rational Field + sage: is_extreme_doubly_nonnegative(A) + True + + Rank 2 is never allowed, so we expect an error:: + + sage: V = VectorSpace(QQ, 5) + sage: A = random_extreme_doubly_nonnegative(V, False, 2) + Traceback (most recent call last): + ... + ValueError: Rank 2 not possible in dimension 5. + + Rank 4 is not allowed in dimension 5:: + + sage: V = VectorSpace(QQ, 5) + sage: A = random_extreme_doubly_nonnegative(V, False, 4) + Traceback (most recent call last): + ... + ValueError: Rank 4 not possible in dimension 5. + + """ + + if not is_admissible_extreme_rank(rank, V.dimension()): + msg = 'Rank %d not possible in dimension %d.' + raise ValueError(msg % (rank, V.dimension())) + + # Generate random doubly-nonnegative matrices until + # one of them is extreme, then return that. + A = random_doubly_nonnegative(V, accept_zero, rank) + + while not is_extreme_doubly_nonnegative(A): + A = random_doubly_nonnegative(V, accept_zero, rank) + + return A