mjo/cone/symmetric_psd.py

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