-"""
+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('../../'))
-from mjo.symbolic import matrix_simplify_full
-
-
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 )
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:
`$D$` will have dimension `$k \times k$`. In the end everything
works out the same.
- EXAMPLES::
+ 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 = 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)
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_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
+ 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 random_symmetric_psd
+
+ EXAMPLES:
+
+ 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()
+ Full MatrixSpace of 2 by 2 dense matrices over Rational Field
+ 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_symmetric_psd(V,False,1)
+ sage: A.rank()
+ 1
+ sage: A = random_symmetric_psd(V,False,2)
+ sage: A.rank()
+ 2
+ sage: A = random_symmetric_psd(V,False,3)
+ sage: A.rank()
+ 3
+ sage: A = random_symmetric_psd(V,False,4)
+ sage: A.rank()
+ 4
+ 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_symmetric_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
+
+ 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