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 unit_eigenvectors(A):
  19     """
  20     Return the unit eigenvectors of a symmetric positive-definite matrix.
  21
  22     INPUT:
  23
  24       - ``A`` - The matrix whose eigenvectors we want to compute.
  25
  26     OUTPUT:
  27
  28     A list of (eigenvalue, eigenvector) pairs where each eigenvector is
  29     associated with its paired eigenvalue of ``A`` and has norm `1`.
  30
  31     EXAMPLES::
  32
  33         sage: A = matrix(QQ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
  34         sage: unit_evs = unit_eigenvectors(A)
  35         sage: bool(unit_evs[0][1].norm() == 1)
  36         True
  37         sage: bool(unit_evs[1][1].norm() == 1)
  38         True
  39         sage: bool(unit_evs[2][1].norm() == 1)
  40         True
  41
  42     """
  43     # This will give us a list of lists whose elements are the
  44     # eigenvectors we want.
  45     ev_lists = [ (val,vecs) for (val,vecs,multiplicity)
  46                             in A.eigenvectors_right() ]
  47
  48     # Pair each eigenvector with its eigenvalue and normalize it.
  49     evs = [ [(l, vec/vec.norm()) for vec in vecs] for (l,vecs) in ev_lists ]
  50
  51     # Flatten the list, abusing the fact that "+" is overloaded on lists.
  52     return sum(evs, [])
  53
  54
  55
  56
  57 def factor_psd(A):
  58     """
  59     Factor a symmetric positive-semidefinite matrix ``A`` into
  60     `XX^{T}`.
  61
  62     INPUT:
  63
  64       - ``A`` - The matrix to factor.
  65
  66     OUTPUT:
  67
  68     A matrix ``X`` such that `A = XX^{T}`.
  69
  70     ALGORITHM:
  71
  72     Since ``A`` is symmetric and positive-semidefinite, we can
  73     diagonalize it by some matrix `$Q$` whose columns are orthogonal
  74     eigenvectors of ``A``. Then,
  75
  76       `$A = QDQ^{T}$`
  77
  78     From this representation we can take the square root of `$D$`
  79     (since all eigenvalues of ``A`` are nonnegative). If we then let
  80     `$X = Q*sqrt(D)*Q^{T}$`, we have,
  81
  82       `$XX^{T} = Q*sqrt(D)*Q^{T}Q*sqrt(D)*Q^{T} = Q*D*Q^{T} = A$`
  83
  84     as desired.
  85
  86     In principle, this is the algorithm used, although we ignore the
  87     eigenvectors corresponding to the eigenvalue zero. Thus if `$rank(A)
  88     = k$`, the matrix `$Q$` will have dimention `$n \times k$`, and
  89     `$D$` will have dimension `$k \times k$`. In the end everything
  90     works out the same.
  91
  92     EXAMPLES::
  93
  94         sage: A = matrix(SR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
  95         sage: X = factor_psd(A)
  96         sage: A2 = matrix_simplify_full(X*X.transpose())
  97         sage: A == A2
  98         True
  99
 100     """
 101
 102     # Get the eigenvectors, and filter out the ones that correspond to
 103     # the eigenvalue zero.
 104     all_evs = unit_eigenvectors(A)
 105     evs = [ (val,vec) for (val,vec) in all_evs if not val == 0 ]
 106
 107     d = [ sqrt(val) for (val,vec) in evs ]
 108     root_D = diagonal_matrix(d).change_ring(A.base_ring())
 109
 110     Q = matrix(A.base_ring(), [ vec for (val,vec) in evs ]).transpose()
 111
 112     return Q*root_D*Q.transpose()