mjo/cone/symmetric_psd.py

   1 r"""
   2 The positive semidefinite cone `$S^{n}_{+}$` is the cone consisting of
   3 all symmetric positive-semidefinite matrices (as a subset of
   4 `$\mathbb{R}^{n \times n}$`
   5 """
   6
   7 from sage.all import *
   8
   9 def unit_eigenvectors(A):
  10     """
  11     Return the unit eigenvectors of a symmetric positive-definite matrix.
  12
  13     INPUT:
  14
  15     - ``A`` -- The matrix whose unit eigenvectors we want to compute.
  16
  17     OUTPUT:
  18
  19     A list of (eigenvalue, eigenvector) pairs where each eigenvector is
  20     associated with its paired eigenvalue of ``A`` and has norm `1`. If
  21     the base ring of ``A`` is not algebraically closed, then returned
  22     eigenvectors may (necessarily) be over its algebraic closure and not
  23     the base ring of ``A`` itself.
  24
  25     SETUP::
  26
  27         sage: from mjo.cone.symmetric_psd import unit_eigenvectors
  28
  29     EXAMPLES::
  30
  31         sage: A = matrix(QQ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
  32         sage: unit_evs = list(unit_eigenvectors(A))
  33         sage: bool(unit_evs[0][1].norm() == 1)
  34         True
  35         sage: bool(unit_evs[1][1].norm() == 1)
  36         True
  37         sage: bool(unit_evs[2][1].norm() == 1)
  38         True
  39
  40     """
  41     return ( (val,vec.normalized())
  42              for (val,vecs,multiplicity) in A.eigenvectors_right()
  43              for vec in vecs )
  44
  45
  46
  47 def factor_psd(A):
  48     """
  49     Factor a symmetric positive-semidefinite matrix ``A`` into
  50     `XX^{T}`.
  51
  52     INPUT:
  53
  54     - ``A`` - The matrix to factor. The base ring of ``A`` must be either
  55               exact or the symbolic ring (to compute eigenvalues), and it
  56               must be a field so that we can take its algebraic closure
  57               (necessary to e.g. take square roots).
  58
  59     OUTPUT:
  60
  61     A matrix ``X`` such that `A = XX^{T}`. The base field of ``X`` will
  62     be the algebraic closure of the base field of ``A``.
  63
  64     ALGORITHM:
  65
  66     Since ``A`` is symmetric and positive-semidefinite, we can
  67     diagonalize it by some matrix `$Q$` whose columns are orthogonal
  68     eigenvectors of ``A``. Then,
  69
  70       `$A = QDQ^{T}$`
  71
  72     From this representation we can take the square root of `$D$`
  73     (since all eigenvalues of ``A`` are nonnegative). If we then let
  74     `$X = Q*sqrt(D)*Q^{T}$`, we have,
  75
  76       `$XX^{T} = Q*sqrt(D)*Q^{T}Q*sqrt(D)*Q^{T} = Q*D*Q^{T} = A$`
  77
  78     as desired.
  79
  80     In principle, this is the algorithm used, although we ignore the
  81     eigenvectors corresponding to the eigenvalue zero. Thus if `$rank(A)
  82     = k$`, the matrix `$Q$` will have dimention `$n \times k$`, and
  83     `$D$` will have dimension `$k \times k$`. In the end everything
  84     works out the same.
  85
  86     SETUP::
  87
  88         sage: from mjo.cone.symmetric_psd import factor_psd
  89
  90     EXAMPLES:
  91
  92     Create a symmetric positive-semidefinite matrix over the symbolic
  93     ring and factor it::
  94
  95         sage: A = matrix(SR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
  96         sage: X = factor_psd(A)
  97         sage: A2 = (X*X.transpose()).simplify_full()
  98         sage: A == A2
  99         True
 100
 101     Attempt to factor the same matrix over ``RR`` which won't work
 102     because ``RR`` isn't exact::
 103
 104         sage: A = matrix(RR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
 105         sage: factor_psd(A)
 106         Traceback (most recent call last):
 107         ...
 108         ValueError: The base ring of ``A`` must be either exact or symbolic.
 109
 110     Attempt to factor the same matrix over ``ZZ`` which won't work
 111     because ``ZZ`` isn't a field::
 112
 113         sage: A = matrix(ZZ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
 114         sage: factor_psd(A)
 115         Traceback (most recent call last):
 116         ...
 117         ValueError: The base ring of ``A`` must be a field.
 118
 119     """
 120
 121     if not A.base_ring().is_exact() and not A.base_ring() is SR:
 122         msg = 'The base ring of ``A`` must be either exact or symbolic.'
 123         raise ValueError(msg)
 124
 125     if not A.base_ring().is_field():
 126         raise ValueError('The base ring of ``A`` must be a field.')
 127
 128     if not A.base_ring() is SR:
 129         # Change the base field of ``A`` so that we are sure we can take
 130         # roots. The symbolic ring has no algebraic_closure method.
 131         A = A.change_ring(A.base_ring().algebraic_closure())
 132
 133
 134     # Get the eigenvectors, and filter out the ones that correspond to
 135     # the eigenvalue zero.
 136     all_evs = unit_eigenvectors(A)
 137     evs = [ (val,vec) for (val,vec) in all_evs if not val == 0 ]
 138
 139     d = ( val.sqrt() for (val,vec) in evs )
 140     root_D = diagonal_matrix(d).change_ring(A.base_ring())
 141
 142     Q = matrix(A.base_ring(), ( vec for (val,vec) in evs )).transpose()
 143
 144     return Q*root_D*Q.transpose()
 145
 146
 147 def random_symmetric_psd(V, accept_zero=True, rank=None):
 148     """
 149     Generate a random symmetric positive-semidefinite matrix over the
 150     vector space ``V``. That is, the returned matrix will be a linear
 151     transformation on ``V``, with the same base ring as ``V``.
 152
 153     We take a very loose interpretation of "random," here. Otherwise we
 154     would never (for example) choose a matrix on the boundary of the
 155     cone (with a zero eigenvalue).
 156
 157     INPUT:
 158
 159     - ``V`` - The vector space on which the returned matrix will act.
 160
 161     - ``accept_zero`` - Do you want to accept the zero matrix (which
 162                         is symmetric PSD? Defaults to ``True``.
 163
 164     - ``rank`` - Require the returned matrix to have the given rank
 165                  (optional).
 166
 167     OUTPUT:
 168
 169     A random symmetric positive semidefinite matrix, i.e. a linear
 170     transformation from ``V`` to itself.
 171
 172     ALGORITHM:
 173
 174     The matrix is constructed from some number of spectral projectors,
 175     which in turn are created at "random" from the underlying vector
 176     space ``V``.
 177
 178     If no particular ``rank`` is desired, we choose the number of
 179     projectors at random. Otherwise, we keep adding new projectors until
 180     the desired rank is achieved.
 181
 182     Finally, before returning, we check if the matrix is zero. If
 183     ``accept_zero`` is ``False``, we restart the process from the
 184     beginning.
 185
 186     SETUP::
 187
 188         sage: from mjo.cone.symmetric_psd import random_symmetric_psd
 189
 190     EXAMPLES:
 191
 192     Well, it doesn't crash at least::
 193
 194         sage: V = VectorSpace(QQ, 2)
 195         sage: A = random_symmetric_psd(V)
 196         sage: A.matrix_space()
 197         Full MatrixSpace of 2 by 2 dense matrices over Rational Field
 198         sage: A.is_positive_semidefinite()
 199         True
 200
 201     A matrix with the desired rank is returned::
 202
 203         sage: V = VectorSpace(QQ, 5)
 204         sage: A = random_symmetric_psd(V,False,1)
 205         sage: A.rank()
 206         1
 207         sage: A = random_symmetric_psd(V,False,2)
 208         sage: A.rank()
 209         2
 210         sage: A = random_symmetric_psd(V,False,3)
 211         sage: A.rank()
 212         3
 213         sage: A = random_symmetric_psd(V,False,4)
 214         sage: A.rank()
 215         4
 216         sage: A = random_symmetric_psd(V,False,5)
 217         sage: A.rank()
 218         5
 219
 220     If the user asks for a rank that's too high, we fail::
 221
 222         sage: V = VectorSpace(QQ, 2)
 223         sage: A = random_symmetric_psd(V,False,3)
 224         Traceback (most recent call last):
 225         ...
 226         ValueError: The ``rank`` must be between 0 and the dimension of ``V``.
 227
 228     """
 229
 230     # We construct the matrix from its spectral projectors. Since
 231     # there can be at most ``n`` of them, where ``n`` is the dimension
 232     # of our vector space, we want to choose a random integer between
 233     # ``0`` and ``n`` and then construct that many random elements of
 234     # ``V``.
 235     n = V.dimension()
 236
 237     rank_A = 0
 238     if rank is None:
 239         # Choose one randomly
 240         rank_A = ZZ.random_element(n+1)
 241     elif (rank < 0) or (rank > n):
 242         # The rank of ``A`` can be at most ``n``.
 243         msg  = 'The ``rank`` must be between 0 and the dimension of ``V``.'
 244         raise ValueError(msg)
 245     else:
 246         # Use the one the user gave us.
 247         rank_A = rank
 248
 249     if n == 0 and not accept_zero:
 250         # We're gonna loop forever trying to satisfy this...
 251         raise ValueError('You must have accept_zero=True when V is trivial')
 252
 253     # Loop until we find a suitable "A" that will then be returned.
 254     while True:
 255         # Begin with the zero matrix, and add projectors to it if we
 256         # have any.
 257         A = matrix.zero(V.base_ring(), n, n)
 258
 259         # Careful, begin at idx=1 so that we only generate a projector
 260         # when rank_A is greater than zero.
 261         while A.rank() < rank_A:
 262             v = V.random_element()
 263             A += v.column()*v.row()
 264
 265         if accept_zero or not A.is_zero():
 266             # We either don't care what ``A`` is, or it's non-zero, so
 267             # just return it.
 268             return A