-"""
+r"""
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}$`
Every completely positive matrix is symmetric
positive-semidefinite::
+ sage: set_random_seed()
sage: v = vector(map(abs, random_vector(ZZ, 10)))
sage: A = v.column() * v.row()
sage: is_symmetric_psd(A)
sage: is_symmetric_psd(A)
False
+ The trivial matrix in a trivial space is trivially symmetric and
+ positive-semidefinite::
+
+ sage: A = matrix(QQ, 0,0)
+ sage: is_symmetric_psd(A)
+ True
+
"""
if A.base_ring() == SR:
# 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
+ if A.is_zero():
+ # A is trivial... so trivially positive-semudefinite.
+ return True
+ else:
+ return min(A.eigenvalues()) >= 0
def unit_eigenvectors(A):
Well, it doesn't crash at least::
+ sage: set_random_seed()
sage: V = VectorSpace(QQ, 2)
sage: A = random_symmetric_psd(V)
sage: A.matrix_space()
A matrix with the desired rank is returned::
+ sage: set_random_seed()
sage: V = VectorSpace(QQ, 5)
sage: A = random_symmetric_psd(V,False,1)
sage: A.rank()
If the user asks for a rank that's too high, we fail::
+ sage: set_random_seed()
sage: V = VectorSpace(QQ, 2)
sage: A = random_symmetric_psd(V,False,3)
Traceback (most recent call last):
# 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_symmetric_psd(V, accept_zero, rank)
+ if n == 0 and not accept_zero:
+ # We're gonna loop forever trying to satisfy this...
+ raise ValueError('You must have accept_zero=True when V is trivial')
+
+ # Loop until we find a suitable "A" that will then be returned.
+ while True:
+ # Begin with the zero matrix, and add projectors to it if we
+ # have any.
+ A = matrix.zero(V.base_ring(), n, n)
+
+ # 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