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