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