""" 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}$` """ from sage.all import * 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. SETUP:: sage: from mjo.cone.symmetric_psd import is_symmetric_psd 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. OUTPUT: A list of (eigenvalue, eigenvector) pairs where each eigenvector is associated with its paired eigenvalue of ``A`` and has norm `1`. 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: bool(unit_evs[0][1].norm() == 1) True sage: bool(unit_evs[1][1].norm() == 1) True sage: bool(unit_evs[2][1].norm() == 1) 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, []) def factor_psd(A): """ Factor a symmetric positive-semidefinite matrix ``A`` into `XX^{T}`. INPUT: - ``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}`. The base field of ``X`` will be the algebraic closure of the base field of ``A``. ALGORITHM: Since ``A`` is symmetric and positive-semidefinite, we can diagonalize it by some matrix `$Q$` whose columns are orthogonal eigenvectors of ``A``. Then, `$A = QDQ^{T}$` From this representation we can take the square root of `$D$` (since all eigenvalues of ``A`` are nonnegative). If we then let `$X = Q*sqrt(D)*Q^{T}$`, we have, `$XX^{T} = Q*sqrt(D)*Q^{T}Q*sqrt(D)*Q^{T} = Q*D*Q^{T} = A$` as desired. In principle, this is the algorithm used, although we ignore the eigenvectors corresponding to the eigenvalue zero. Thus if `$rank(A) = k$`, the matrix `$Q$` will have dimention `$n \times k$`, and `$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 ring and factor it:: sage: A = matrix(SR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]]) sage: X = factor_psd(A) 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) evs = [ (val,vec) for (val,vec) in all_evs if not val == 0 ] d = [ sqrt(val) 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() 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. SETUP:: sage: from mjo.cone.symmetric_psd import is_symmetric_psd, random_psd 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)