X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fdoubly_nonnegative.py;h=e8cacda5cd492100096eff2cf557c8992a34d9b3;hb=1bbade9f41ffbfe366b15d0db657f666bc1f025d;hp=ebb5b685fd3c5e1843fb4379e7920316dc159a15;hpb=7a884aa96ce490425b82d977d3242954511efbdf;p=sage.d.git diff --git a/mjo/cone/doubly_nonnegative.py b/mjo/cone/doubly_nonnegative.py index ebb5b68..e8cacda 100644 --- a/mjo/cone/doubly_nonnegative.py +++ b/mjo/cone/doubly_nonnegative.py @@ -1,4 +1,4 @@ -""" +r""" The doubly-nonnegative cone in `S^{n}` is the set of all such matrices that both, @@ -13,14 +13,10 @@ It is represented typically by either `\mathcal{D}^{n}` or 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 our module names -# resolve. -from os.path import abspath -from site import addsitedir -addsitedir(abspath('../../')) -from mjo.cone.symmetric_psd import factor_psd, is_symmetric_psd -from mjo.matrix_vector import isomorphism +from mjo.cone.symmetric_psd import (factor_psd, + is_symmetric_psd, + random_symmetric_psd) +from mjo.basis_repr import basis_repr def is_doubly_nonnegative(A): @@ -36,6 +32,10 @@ def is_doubly_nonnegative(A): Either ``True`` if ``A`` is doubly-nonnegative, or ``False`` otherwise. + SETUP:: + + sage: from mjo.cone.doubly_nonnegative import is_doubly_nonnegative + EXAMPLES: Every completely positive matrix is doubly-nonnegative:: @@ -58,7 +58,7 @@ def is_doubly_nonnegative(A): raise ValueError.new(msg) # Check that all of the entries of ``A`` are nonnegative. - if not all([ a >= 0 for a in A.list() ]): + if not all( a >= 0 for a in A.list() ): return False # It's nonnegative, so all we need to do is check that it's @@ -68,7 +68,7 @@ def is_doubly_nonnegative(A): def is_admissible_extreme_rank(r, n): - """ + r""" 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. @@ -85,6 +85,10 @@ def is_admissible_extreme_rank(r, n): the doubly-nonnegative cone in `$\mathbb{R}^{n}$`, or ``False`` otherwise. + SETUP:: + + sage: from mjo.cone.doubly_nonnegative import is_admissible_extreme_rank + EXAMPLES: For dimension 5, only ranks zero, one, and three are admissible:: @@ -150,6 +154,10 @@ def has_admissible_extreme_rank(A): 26 (1996), no. 4, 1371--1383. doi:10.1216/rmjm/1181071993. http://projecteuclid.org/euclid.rmjm/1181071993. + SETUP:: + + sage: from mjo.cone.doubly_nonnegative import has_admissible_extreme_rank + EXAMPLES: The zero matrix has rank zero, which is admissible:: @@ -196,7 +204,7 @@ def has_admissible_extreme_rank(A): return is_admissible_extreme_rank(r,n) -def E(matrix_space, i,j): +def stdE(matrix_space, i,j): """ Return the ``i``,``j``th element of the standard basis in ``matrix_space``. @@ -215,26 +223,30 @@ def E(matrix_space, i,j): A basis element of ``matrix_space``. It has a single \"1\" in the ``i``,``j`` row,column and zeros elsewhere. + SETUP:: + + sage: from mjo.cone.doubly_nonnegative import stdE + EXAMPLES:: sage: M = MatrixSpace(ZZ, 2, 2) - sage: E(M,0,0) + sage: stdE(M,0,0) [1 0] [0 0] - sage: E(M,0,1) + sage: stdE(M,0,1) [0 1] [0 0] - sage: E(M,1,0) + sage: stdE(M,1,0) [0 0] [1 0] - sage: E(M,1,1) + sage: stdE(M,1,1) [0 0] [0 1] - sage: E(M,2,1) + sage: stdE(M,2,1) Traceback (most recent call last): ... IndexError: Index `i` is out of bounds. - sage: E(M,1,2) + sage: stdE(M,1,2) Traceback (most recent call last): ... IndexError: Index `j` is out of bounds. @@ -253,7 +265,7 @@ def E(matrix_space, i,j): # would be computed as offset 3 into a four-element list and we # would succeed incorrectly. idx = matrix_space.ncols()*i + j - return matrix_space.basis()[idx] + return list(matrix_space.basis())[idx] @@ -272,6 +284,10 @@ def is_extreme_doubly_nonnegative(A): 2. Berman, Abraham and Shaked-Monderer, Naomi. Completely Positive Matrices. World Scientific, 2003. + SETUP:: + + sage: from mjo.cone.doubly_nonnegative import is_extreme_doubly_nonnegative + EXAMPLES: The zero matrix is an extreme matrix:: @@ -355,11 +371,11 @@ def is_extreme_doubly_nonnegative(A): # whenever we come across an index pair `$(i,j)$` with # `$A_{ij} = 0$`. spanning_set = [] - for j in range(0, A.ncols()): - for i in range(0,j): + for j in range(A.ncols()): + for i in range(j): if A[i,j] == 0: M = A.matrix_space() - S = X.transpose() * (E(M,i,j) + E(M,j,i)) * X + S = X.transpose() * (stdE(M,i,j) + stdE(M,j,i)) * X spanning_set.append(S) # The spanning set that we have at this point is of matrices. We @@ -367,7 +383,7 @@ def is_extreme_doubly_nonnegative(A): # can't compute the dimension of a set of matrices anyway, so we # convert them all to vectors and just ask for the dimension of the # resulting vector space. - (phi, phi_inverse) = isomorphism(A.matrix_space()) + (phi, phi_inverse) = basis_repr(A.matrix_space()) vectors = map(phi,spanning_set) V = span(vectors, A.base_ring()) @@ -378,3 +394,151 @@ 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. + + SETUP:: + + sage: from mjo.cone.doubly_nonnegative import (is_doubly_nonnegative, + ....: random_doubly_nonnegative) + + 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_symmetric_psd(V, accept_zero, rank) + + while not all( x >= 0 for x in A.list() ): + A = random_symmetric_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. + + SETUP:: + + sage: from mjo.cone.doubly_nonnegative import (is_extreme_doubly_nonnegative, + ....: random_extreme_doubly_nonnegative) + + 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 rank is not None and 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