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