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
  16
  17 def is_symmetric_psd(A):
  18     """
  19     Determine whether or not the matrix ``A`` is symmetric
  20     positive-semidefinite.
  21
  22     INPUT:
  23
  24     - ``A`` - The matrix in question
  25
  26     OUTPUT:
  27
  28     Either ``True`` if ``A`` is symmetric positive-semidefinite, or
  29     ``False`` otherwise.
  30
  31     EXAMPLES:
  32
  33     Every completely positive matrix is symmetric
  34     positive-semidefinite::
  35
  36         sage: v = vector(map(abs, random_vector(ZZ, 10)))
  37         sage: A = v.column() * v.row()
  38         sage: is_symmetric_psd(A)
  39         True
  40
  41     The following matrix is symmetric but not positive semidefinite::
  42
  43         sage: A = matrix(ZZ, [[1, 2], [2, 1]])
  44         sage: is_symmetric_psd(A)
  45         False
  46
  47     This matrix isn't even symmetric::
  48
  49         sage: A = matrix(ZZ, [[1, 2], [3, 4]])
  50         sage: is_symmetric_psd(A)
  51         False
  52
  53     """
  54
  55     if A.base_ring() == SR:
  56         msg = 'The matrix ``A`` cannot be symbolic.'
  57         raise ValueError.new(msg)
  58
  59     # First make sure that ``A`` is symmetric.
  60     if not A.is_symmetric():
  61         return False
  62
  63     # If ``A`` is symmetric, we only need to check that it is positive
  64     # semidefinite. For that we can consult its minimum eigenvalue,
  65     # which should be zero or greater. Since ``A`` is symmetric, its
  66     # eigenvalues are guaranteed to be real.
  67     return min(A.eigenvalues()) >= 0
  68
  69
  70 def unit_eigenvectors(A):
  71     """
  72     Return the unit eigenvectors of a symmetric positive-definite matrix.
  73
  74     INPUT:
  75
  76     - ``A`` - The matrix whose eigenvectors we want to compute.
  77
  78     OUTPUT:
  79
  80     A list of (eigenvalue, eigenvector) pairs where each eigenvector is
  81     associated with its paired eigenvalue of ``A`` and has norm `1`.
  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     EXAMPLES:
 149
 150     Create a symmetric positive-semidefinite matrix over the symbolic
 151     ring and factor it::
 152
 153         sage: A = matrix(SR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
 154         sage: X = factor_psd(A)
 155         sage: A2 = (X*X.transpose()).simplify_full()
 156         sage: A == A2
 157         True
 158
 159     Attempt to factor the same matrix over ``RR`` which won't work
 160     because ``RR`` isn't exact::
 161
 162         sage: A = matrix(RR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
 163         sage: factor_psd(A)
 164         Traceback (most recent call last):
 165         ...
 166         ValueError: The base ring of ``A`` must be either exact or symbolic.
 167
 168     Attempt to factor the same matrix over ``ZZ`` which won't work
 169     because ``ZZ`` isn't a field::
 170
 171         sage: A = matrix(ZZ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
 172         sage: factor_psd(A)
 173         Traceback (most recent call last):
 174         ...
 175         ValueError: The base ring of ``A`` must be a field.
 176
 177     """
 178
 179     if not A.base_ring().is_exact() and not A.base_ring() is SR:
 180         msg = 'The base ring of ``A`` must be either exact or symbolic.'
 181         raise ValueError(msg)
 182
 183     if not A.base_ring().is_field():
 184         raise ValueError('The base ring of ``A`` must be a field.')
 185
 186     if not A.base_ring() is SR:
 187         # Change the base field of ``A`` so that we are sure we can take
 188         # roots. The symbolic ring has no algebraic_closure method.
 189         A = A.change_ring(A.base_ring().algebraic_closure())
 190
 191
 192     # Get the eigenvectors, and filter out the ones that correspond to
 193     # the eigenvalue zero.
 194     all_evs = unit_eigenvectors(A)
 195     evs = [ (val,vec) for (val,vec) in all_evs if not val == 0 ]
 196
 197     d = [ sqrt(val) for (val,vec) in evs ]
 198     root_D = diagonal_matrix(d).change_ring(A.base_ring())
 199
 200     Q = matrix(A.base_ring(), [ vec for (val,vec) in evs ]).transpose()
 201
 202     return Q*root_D*Q.transpose()
 203
 204
 205 def random_psd(V, accept_zero=True, rank=None):
 206     """
 207     Generate a random symmetric positive-semidefinite matrix over the
 208     vector space ``V``. That is, the returned matrix will be a linear
 209     transformation on ``V``, with the same base ring as ``V``.
 210
 211     We take a very loose interpretation of "random," here. Otherwise we
 212     would never (for example) choose a matrix on the boundary of the
 213     cone (with a zero eigenvalue).
 214
 215     INPUT:
 216
 217     - ``V`` - The vector space on which the returned matrix will act.
 218
 219     - ``accept_zero`` - Do you want to accept the zero matrix (which
 220                         is symmetric PSD? Defaults to ``True``.
 221
 222     - ``rank`` - Require the returned matrix to have the given rank
 223                  (optional).
 224
 225     OUTPUT:
 226
 227     A random symmetric positive semidefinite matrix, i.e. a linear
 228     transformation from ``V`` to itself.
 229
 230     ALGORITHM:
 231
 232     The matrix is constructed from some number of spectral projectors,
 233     which in turn are created at "random" from the underlying vector
 234     space ``V``.
 235
 236     If no particular ``rank`` is desired, we choose the number of
 237     projectors at random. Otherwise, we keep adding new projectors until
 238     the desired rank is achieved.
 239
 240     Finally, before returning, we check if the matrix is zero. If
 241     ``accept_zero`` is ``False``, we restart the process from the
 242     beginning.
 243
 244     EXAMPLES:
 245
 246     Well, it doesn't crash at least::
 247
 248         sage: V = VectorSpace(QQ, 2)
 249         sage: A = random_psd(V)
 250         sage: A.matrix_space()
 251         Full MatrixSpace of 2 by 2 dense matrices over Rational Field
 252         sage: is_symmetric_psd(A)
 253         True
 254
 255     A matrix with the desired rank is returned::
 256
 257         sage: V = VectorSpace(QQ, 5)
 258         sage: A = random_psd(V,False,1)
 259         sage: A.rank()
 260         1
 261         sage: A = random_psd(V,False,2)
 262         sage: A.rank()
 263         2
 264         sage: A = random_psd(V,False,3)
 265         sage: A.rank()
 266         3
 267         sage: A = random_psd(V,False,4)
 268         sage: A.rank()
 269         4
 270         sage: A = random_psd(V,False,5)
 271         sage: A.rank()
 272         5
 273
 274     If the user asks for a rank that's too high, we fail::
 275
 276         sage: V = VectorSpace(QQ, 2)
 277         sage: A = random_psd(V,False,3)
 278         Traceback (most recent call last):
 279         ...
 280         ValueError: The ``rank`` must be between 0 and the dimension of ``V``.
 281
 282     """
 283
 284     # We construct the matrix from its spectral projectors. Since
 285     # there can be at most ``n`` of them, where ``n`` is the dimension
 286     # of our vector space, we want to choose a random integer between
 287     # ``0`` and ``n`` and then construct that many random elements of
 288     # ``V``.
 289     n = V.dimension()
 290
 291     rank_A = 0
 292     if rank is None:
 293         # Choose one randomly
 294         rank_A = ZZ.random_element(n+1)
 295     elif (rank < 0) or (rank > n):
 296         # The rank of ``A`` can be at most ``n``.
 297         msg  = 'The ``rank`` must be between 0 and the dimension of ``V``.'
 298         raise ValueError(msg)
 299     else:
 300         # Use the one the user gave us.
 301         rank_A = rank
 302
 303     # Begin with the zero matrix, and add projectors to it if we have any.
 304     A = V.zero().column()*V.zero().row()
 305
 306     # Careful, begin at idx=1 so that we only generate a projector
 307     # when rank_A is greater than zero.
 308     while A.rank() < rank_A:
 309         v = V.random_element()
 310         A += v.column()*v.row()
 311
 312     if accept_zero or not A.is_zero():
 313         # We either don't care what ``A`` is, or it's non-zero, so
 314         # just return it.
 315         return A
 316     else:
 317         # Uh oh, we need to generate a new one.
 318         return random_psd(V, accept_zero, rank)