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