X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fsymmetric_psd.py;h=e4be629447987bca42a6d3e795f1b9926bc378c8;hb=a4109ab945f5d3ed94207c936e60b5b187ae450b;hp=fae0ae2307bb4e8b0583cf5b2200951f66a7e256;hpb=bc32004c617931e80f5c87775086e3f0e05f502e;p=sage.d.git diff --git a/mjo/cone/symmetric_psd.py b/mjo/cone/symmetric_psd.py index fae0ae2..e4be629 100644 --- a/mjo/cone/symmetric_psd.py +++ b/mjo/cone/symmetric_psd.py @@ -6,13 +6,57 @@ 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): @@ -21,7 +65,7 @@ def unit_eigenvectors(A): INPUT: - - ``A`` - The matrix whose eigenvectors we want to compute. + - ``A`` - The matrix whose eigenvectors we want to compute. OUTPUT: @@ -61,11 +105,15 @@ def factor_psd(A): INPUT: - - ``A`` - The matrix to factor. + - ``A`` - The matrix to factor. The base ring of ``A`` must be either + exact or the symbolic ring (to compute eigenvalues), and it + must be a field so that we can take its algebraic closure + (necessary to e.g. take square roots). OUTPUT: - A matrix ``X`` such that `A = XX^{T}`. + A matrix ``X`` such that `A = XX^{T}`. The base field of ``X`` will + be the algebraic closure of the base field of ``A``. ALGORITHM: @@ -89,16 +137,50 @@ def factor_psd(A): `$D$` will have dimension `$k \times k$`. In the end everything works out the same. - EXAMPLES:: + EXAMPLES: + + Create a symmetric positive-semidefinite matrix over the symbolic + ring and factor it:: 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 + Attempt to factor the same matrix over ``RR`` which won't work + because ``RR`` isn't exact:: + + sage: A = matrix(RR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]]) + sage: factor_psd(A) + Traceback (most recent call last): + ... + ValueError: The base ring of ``A`` must be either exact or symbolic. + + Attempt to factor the same matrix over ``ZZ`` which won't work + because ``ZZ`` isn't a field:: + + sage: A = matrix(ZZ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]]) + sage: factor_psd(A) + Traceback (most recent call last): + ... + ValueError: The base ring of ``A`` must be a field. + """ + 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()) + + # Get the eigenvectors, and filter out the ones that correspond to # the eigenvalue zero. all_evs = unit_eigenvectors(A) @@ -110,3 +192,119 @@ def factor_psd(A): 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): + """ + 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. + + EXAMPLES: + + Well, it doesn't crash at least:: + + sage: V = VectorSpace(QQ, 2) + sage: A = random_psd(V) + sage: A.matrix_space() + Full MatrixSpace of 2 by 2 dense matrices over Rational Field + sage: is_symmetric_psd(A) + True + + A matrix with the desired rank is returned:: + + sage: V = VectorSpace(QQ, 5) + sage: A = random_psd(V,False,1) + sage: A.rank() + 1 + sage: A = random_psd(V,False,2) + sage: A.rank() + 2 + sage: A = random_psd(V,False,3) + sage: A.rank() + 3 + sage: A = random_psd(V,False,4) + sage: A.rank() + 4 + sage: A = random_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) + 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 + + # 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)