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