X-Git-Url: https://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fdoubly_nonnegative.py;h=4f7f950bb2c897509c2cb4b41d5200bccfd67c49;hb=a5700d65325514a505d24fb96b75c2b0f2f6e94f;hp=f705e5e98b083d6cf85fac39148d171cdf452d99;hpb=2d95c4e34d085c9647c73b73b2957f937cfee26b;p=sage.d.git diff --git a/mjo/cone/doubly_nonnegative.py b/mjo/cone/doubly_nonnegative.py index f705e5e..4f7f950 100644 --- a/mjo/cone/doubly_nonnegative.py +++ b/mjo/cone/doubly_nonnegative.py @@ -14,13 +14,13 @@ 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 "mjo.cone" -# resolves. +# 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 def is_doubly_nonnegative(A): @@ -29,7 +29,7 @@ def is_doubly_nonnegative(A): INPUT: - - ``A`` - The matrix in question + - ``A`` - The matrix in question OUTPUT: @@ -67,9 +67,240 @@ def is_doubly_nonnegative(A): +def has_admissible_extreme_rank(A): + """ + The extreme matrices of the doubly-nonnegative cone have some + restrictions on their ranks. This function checks to see whether or + not ``A`` could be extreme based on its rank. + + INPUT: + + - ``A`` - The matrix in question + + OUTPUT: + + ``False`` if the rank of ``A`` precludes it from being an extreme + matrix of the doubly-nonnegative cone, ``True`` otherwise. + + REFERENCE: + + Hamilton-Jester, Crista Lee; Li, Chi-Kwong. Extreme Vectors of + Doubly Nonnegative Matrices. Rocky Mountain Journal of Mathematics + 26 (1996), no. 4, 1371--1383. doi:10.1216/rmjm/1181071993. + http://projecteuclid.org/euclid.rmjm/1181071993. + + EXAMPLES: + + The zero matrix has rank zero, which is admissible:: + + sage: A = zero_matrix(QQ, 5, 5) + sage: has_admissible_extreme_rank(A) + True + + """ + if not A.is_symmetric(): + # This function is more or less internal, so blow up if passed + # something unexpected. + raise ValueError('The matrix ``A`` must be symmetric.') + + 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) + + +def E(matrix_space, i,j): + """ + Return the ``i``,``j``th element of the standard basis in + ``matrix_space``. + + INPUT: + + - ``matrix_space`` - The underlying matrix space of whose basis + the returned matrix is an element + + - ``i`` - The row index of the single nonzero entry + + - ``j`` - The column index of the single nonzero entry + + OUTPUT: + + A basis element of ``matrix_space``. It has a single \"1\" in the + ``i``,``j`` row,column and zeros elsewhere. + + EXAMPLES:: + + sage: M = MatrixSpace(ZZ, 2, 2) + sage: E(M,0,0) + [1 0] + [0 0] + sage: E(M,0,1) + [0 1] + [0 0] + sage: E(M,1,0) + [0 0] + [1 0] + sage: E(M,1,1) + [0 0] + [0 1] + sage: E(M,2,1) + Traceback (most recent call last): + ... + IndexError: Index `i` is out of bounds. + sage: E(M,1,2) + Traceback (most recent call last): + ... + IndexError: Index `j` is out of bounds. + + """ + # We need to check these ourselves, see below. + if i >= matrix_space.nrows(): + raise IndexError('Index `i` is out of bounds.') + if j >= matrix_space.ncols(): + raise IndexError('Index `j` is out of bounds.') + + # The basis here is returned as a one-dimensional list, so we need + # to compute the offset into it based on ``i`` and ``j``. Since we + # compute the index ourselves, we need to do bounds-checking + # manually. Otherwise for e.g. a 2x2 matrix space, the index (0,2) + # 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] + + + def is_extreme_doubly_nonnegative(A): """ Returns ``True`` if the given matrix is an extreme matrix of the doubly-nonnegative cone, and ``False`` otherwise. + + REFERENCES: + + 1. Hamilton-Jester, Crista Lee; Li, Chi-Kwong. Extreme Vectors of + Doubly Nonnegative Matrices. Rocky Mountain Journal of Mathematics + 26 (1996), no. 4, 1371--1383. doi:10.1216/rmjm/1181071993. + http://projecteuclid.org/euclid.rmjm/1181071993. + + 2. Berman, Abraham and Shaked-Monderer, Naomi. Completely Positive + Matrices. World Scientific, 2003. + + EXAMPLES: + + The zero matrix is an extreme matrix:: + + sage: A = zero_matrix(QQ, 5, 5) + sage: is_extreme_doubly_nonnegative(A) + True + + Any extreme vector of the completely positive cone is an extreme + vector of the doubly-nonnegative cone:: + + sage: v = vector([1,2,3,4,5,6]) + sage: A = v.column() * v.row() + sage: A = A.change_ring(QQ) + sage: is_extreme_doubly_nonnegative(A) + True + + We should be able to generate the extreme completely positive + vectors randomly:: + + sage: v = vector(map(abs, random_vector(ZZ, 4))) + sage: A = v.column() * v.row() + sage: A = A.change_ring(QQ) + sage: is_extreme_doubly_nonnegative(A) + True + sage: v = vector(map(abs, random_vector(ZZ, 10))) + sage: A = v.column() * v.row() + sage: A = A.change_ring(QQ) + sage: is_extreme_doubly_nonnegative(A) + True + + The following matrix is completely positive but has rank 3, so by a + remark in reference #1 it is not extreme:: + + sage: A = matrix(QQ, [[1,2,1],[2,6,3],[1,3,5]]) + sage: is_extreme_doubly_nonnegative(A) + False + + The following matrix is completely positive (diagonal) with rank 2, + so it is also not extreme:: + + sage: A = matrix(QQ, [[1,0,0],[2,0,0],[0,0,0]]) + sage: is_extreme_doubly_nonnegative(A) + False + """ - raise NotImplementedError() + + if not A.base_ring().is_exact() and not A.base_ring() is SR: + 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 method. + A = A.change_ring(A.base_ring().algebraic_closure()) + + # Step 1 (see reference #1) + k = A.rank() + + if k == 0: + # Short circuit, we know the zero matrix is extreme. + return True + + if not is_symmetric_psd(A): + return False + + # Step 1.5, appeal to Theorem 3.1 in reference #1 to short + # circuit. + if not has_admissible_extreme_rank(A): + return False + + # Step 2 + X = factor_psd(A) + + # Step 3 + # + # Begin with an empty spanning set, and add a new matrix to it + # 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): + if A[i,j] == 0: + M = A.matrix_space() + S = X.transpose() * (E(M,i,j) + E(M,j,i)) * X + spanning_set.append(S) + + # The spanning set that we have at this point is of matrices. We + # only care about the dimension of the spanned space, and Sage + # 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()) + vectors = map(phi,spanning_set) + + V = span(vectors, A.base_ring()) + d = V.dimension() + + # Needed to safely divide by two here (we don't want integer + # division). We ensured that the base ring of ``A`` is a field + # earlier. + two = A.base_ring()(2) + return d == (k*(k + 1)/two - 1)