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         raise ValueError('The base ring of ``A`` must be either exact or symbolic.')
 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.
 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()