-"""
+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}$`
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('../../'))
-
-
-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):
"""
Return the unit eigenvectors of a symmetric positive-definite matrix.
INPUT:
- - ``A`` - The matrix whose eigenvectors we want to compute.
+ - ``A`` -- The matrix whose unit 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`.
+ associated with its paired eigenvalue of ``A`` and has norm `1`. If
+ the base ring of ``A`` is not algebraically closed, then returned
+ eigenvectors may (necessarily) be over its algebraic closure and not
+ the base ring of ``A`` itself.
+
+ 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: unit_evs = list(unit_eigenvectors(A))
sage: bool(unit_evs[0][1].norm() == 1)
True
sage: bool(unit_evs[1][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, [])
-
+ return ( (val,vec.normalized())
+ for (val,vecs,multiplicity) in A.eigenvectors_right()
+ for vec in vecs )
`$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
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 ]
+ d = ( val.sqrt() 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()
+ 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):
+def random_symmetric_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
``accept_zero`` is ``False``, we restart the process from the
beginning.
+ SETUP::
+
+ sage: from mjo.cone.symmetric_psd import random_symmetric_psd
+
EXAMPLES:
Well, it doesn't crash at least::
+ sage: set_random_seed()
sage: V = VectorSpace(QQ, 2)
- sage: A = random_psd(V)
+ sage: A = random_symmetric_psd(V)
sage: A.matrix_space()
Full MatrixSpace of 2 by 2 dense matrices over Rational Field
- sage: is_symmetric_psd(A)
+ sage: A.is_positive_semidefinite()
True
A matrix with the desired rank is returned::
+ sage: set_random_seed()
sage: V = VectorSpace(QQ, 5)
- sage: A = random_psd(V,False,1)
+ sage: A = random_symmetric_psd(V,False,1)
sage: A.rank()
1
- sage: A = random_psd(V,False,2)
+ sage: A = random_symmetric_psd(V,False,2)
sage: A.rank()
2
- sage: A = random_psd(V,False,3)
+ sage: A = random_symmetric_psd(V,False,3)
sage: A.rank()
3
- sage: A = random_psd(V,False,4)
+ sage: A = random_symmetric_psd(V,False,4)
sage: A.rank()
4
- sage: A = random_psd(V,False,5)
+ sage: A = random_symmetric_psd(V,False,5)
sage: A.rank()
5
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_psd(V,False,3)
+ sage: A = random_symmetric_psd(V,False,3)
Traceback (most recent call last):
...
ValueError: The ``rank`` must be between 0 and the dimension of ``V``.
# 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_element().column()*V.zero_element().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)
+ 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
projects / sage.d.git / blobdiff