mjo/ldlt.py

   1 from sage.all import *
   2
   3 def is_positive_semidefinite_naive(A):
   4     r"""
   5     A naive positive-semidefinite check that tests the eigenvalues for
   6     nonnegativity.  We follow the sage convention that positive
   7     (semi)definite matrices must be symmetric or Hermitian.
   8
   9     SETUP::
  10
  11         sage: from mjo.ldlt import is_positive_semidefinite_naive
  12
  13     TESTS:
  14
  15     The trivial matrix is vaciously positive-semidefinite::
  16
  17         sage: A = matrix(QQ, 0)
  18         sage: A
  19         []
  20         sage: is_positive_semidefinite_naive(A)
  21         True
  22
  23     """
  24     if A.nrows() == 0:
  25         return True # vacuously
  26     return A.is_hermitian() and all( v >= 0 for v in A.eigenvalues() )
  27
  28
  29 def ldlt_naive(A):
  30     r"""
  31     Perform a pivoted `LDL^{T}` factorization of the Hermitian
  32     positive-semidefinite matrix `A`.
  33
  34     This is a naive, recursive implementation that is inefficient due
  35     to Python's lack of tail-call optimization. The pivot strategy is
  36     to choose the largest diagonal entry of the matrix at each step,
  37     and to permute it into the top-left position. Ultimately this
  38     results in a factorization `A = PLDL^{T}P^{T}`, where `P` is a
  39     permutation matrix, `L` is unit-lower-triangular, and `D` is
  40     diagonal decreasing from top-left to bottom-right.
  41
  42     ALGORITHM:
  43
  44     The algorithm is based on the discussion in Golub and Van Loan, but with
  45     some "typos" fixed.
  46
  47     OUTPUT:
  48
  49     A triple `(P,L,D)` such that `A = PLDL^{T}P^{T}` and where,
  50
  51       * `P` is a permutaiton matrix
  52       * `L` is unit lower-triangular
  53       * `D` is a diagonal matrix whose entries are decreasing from top-left
  54         to bottom-right
  55
  56     SETUP::
  57
  58         sage: from mjo.ldlt import ldlt_naive, is_positive_semidefinite_naive
  59
  60     EXAMPLES:
  61
  62     All three factors should be the identity when the original matrix is::
  63
  64         sage: I = matrix.identity(QQ,4)
  65         sage: P,L,D = ldlt_naive(I)
  66         sage: P == I and L == I and D == I
  67         True
  68
  69     TESTS:
  70
  71     Ensure that a "random" positive-semidefinite matrix is factored correctly::
  72
  73         sage: set_random_seed()
  74         sage: n = ZZ.random_element(5)
  75         sage: A = matrix.random(QQ, n)
  76         sage: A = A*A.transpose()
  77         sage: is_positive_semidefinite_naive(A)
  78         True
  79         sage: P,L,D = ldlt_naive(A)
  80         sage: A == P*L*D*L.transpose()*P.transpose()
  81         True
  82
  83     """
  84     n = A.nrows()
  85
  86     # Use the fraction field of the given matrix so that division will work
  87     # when (for example) our matrix consists of integer entries.
  88     ring = A.base_ring().fraction_field()
  89
  90     if n == 0 or n == 1:
  91         # We can get n == 0 if someone feeds us a trivial matrix.
  92         P = matrix.identity(ring, n)
  93         L = matrix.identity(ring, n)
  94         D = A
  95         return (P,L,D)
  96
  97     A1 = A.change_ring(ring)
  98     diags = A1.diagonal()
  99     s = diags.index(max(diags))
 100     P1 = copy(A1.matrix_space().identity_matrix())
 101     P1.swap_rows(0,s)
 102     A1 = P1.T * A1 * P1
 103     alpha1 = A1[0,0]
 104
 105     # Golub and Van Loan mention in passing what to do here. This is
 106     # only sensible if the matrix is positive-semidefinite, because we
 107     # are assuming that we can set everything else to zero as soon as
 108     # we hit the first on-diagonal zero.
 109     if alpha1 == 0:
 110         P = A1.matrix_space().identity_matrix()
 111         L = P
 112         D = A1.matrix_space().zero()
 113         return (P,L,D)
 114
 115     v1 = A1[1:n,0]
 116     A2 = A1[1:,1:]
 117
 118     P2, L2, D2 = ldlt_naive(A2 - (v1*v1.transpose())/alpha1)
 119
 120     P1 = P1*block_matrix(2,2, [[ZZ(1), ZZ(0)],
 121                                [0*v1,  P2]])
 122     L1 = block_matrix(2,2, [[ZZ(1),                    ZZ(0)],
 123                             [P2.transpose()*v1/alpha1, L2]])
 124     D1 = block_matrix(2,2, [[alpha1, ZZ(0)],
 125                             [0*v1,   D2]])
 126
 127     return (P1,L1,D1)
 128
 129
 130
 131 def ldlt_fast(A):
 132     r"""
 133     Perform a fast, pivoted `LDL^{T}` factorization of the Hermitian
 134     positive-semidefinite matrix `A`.
 135
 136     This function is much faster than ``ldlt_naive`` because the
 137     tail-recursion has been unrolled into a loop.
 138     """
 139     ring = A.base_ring().fraction_field()
 140     A = A.change_ring(ring)
 141
 142     # Keep track of the permutations in a vector rather than in a
 143     # matrix, for efficiency.
 144     n = A.nrows()
 145     p = list(range(n))
 146
 147     for k in range(n):
 148         # We need to loop once for every diagonal entry in the
 149         # matrix. So, as many times as it has rows/columns. At each
 150         # step, we obtain the permutation needed to put things in the
 151         # right place, then the "next" entry (alpha) of D, and finally
 152         # another column of L.
 153         diags = A.diagonal()[k:n]
 154         alpha = max(diags)
 155
 156         # We're working *within* the matrix ``A``, so every index is
 157         # offset by k. For example: after the second step, we should
 158         # only be looking at the lower 3-by-3 block of a 5-by-5 matrix.
 159         s = k + diags.index(alpha)
 160
 161         # Move the largest diagonal element up into the top-left corner
 162         # of the block we're working on (the one starting from index k,k).
 163         # Presumably this is faster than hitting the thing with a
 164         # permutation matrix.
 165         #
 166         # Since "L" is stored in the lower-left "half" of "A", it's a
 167         # good thing that we need to permuts "L," too. This is due to
 168         # how P2.T appears in the recursive algorithm applied to the
 169         # "current" column of L There, P2.T is computed recusively, as
 170         # 1 x P3.T, and P3.T = 1 x P4.T, etc, from the bottom up. All
 171         # are eventually applied to "v" in order.  Here we're working
 172         # from the top down, and rather than keep track of what
 173         # permutations we need to perform, we just perform them as we
 174         # go along. No recursion needed.
 175         A.swap_columns(k,s)
 176         A.swap_rows(k,s)
 177
 178         # Update the permutation "matrix" with the swap we just did.
 179         p_k = p[k]
 180         p[k] = p[s]
 181         p[s] = p_k
 182
 183         # Now the largest diagonal is in the top-left corner of the
 184         # block below and to the right of index k,k. When alpha is
 185         # zero, we can just leave the rest of the D/L entries
 186         # zero... which is exactly how they start out.
 187         if alpha != 0:
 188             # Update the "next" block of A that we'll work on during
 189             # the following iteration. I think it's faster to get the
 190             # entries of a row than a column here?
 191             for i in range(n-k-1):
 192                 for j in range(i+1):
 193                     A[k+1+j,k+1+i] = A[k+1+j,k+1+i] - A[k,k+1+j]*A[k,k+1+i]/alpha
 194                     A[k+1+i,k+1+j] = A[k+1+j,k+1+i] # keep it symmetric!
 195
 196             for i in range(n-k-1):
 197                 # Store the "new" (kth) column of L, being sure to set
 198                 # the lower-left "half" from the upper-right "half"
 199                 A[k+i+1,k] = A[k,k+1+i]/alpha
 200
 201     MS = A.matrix_space()
 202     P = MS.matrix(lambda i,j: p[j] == i)
 203     D = MS.diagonal_matrix(A.diagonal())
 204
 205     for i in range(n):
 206         A[i,i] = 1
 207         for j in range(i+1,n):
 208             A[i,j] = 0
 209
 210     return P,A,D
 211
 212
 213 def block_ldlt_naive(A, check_hermitian=False):
 214     r"""
 215     Perform a block-`LDL^{T}` factorization of the Hermitian
 216     matrix `A`.
 217
 218     This is a naive, recursive implementation akin to
 219     ``ldlt_naive()``, where the pivots (and resulting diagonals) are
 220     either `1 \times 1` or `2 \times 2` blocks. The pivots are chosen
 221     using the Bunch-Kaufmann scheme that is both fast and numerically
 222     stable.
 223
 224     OUTPUT:
 225
 226     A triple `(P,L,D)` such that `A = PLDL^{T}P^{T}` and where,
 227
 228       * `P` is a permutation matrix
 229       * `L` is unit lower-triangular
 230       * `D` is a block-diagonal matrix whose entries are decreasing
 231         from top-left to bottom-right and whose blocks are of size
 232         one or two.
 233     """
 234     n = A.nrows()
 235
 236     # Use the fraction field of the given matrix so that division will work
 237     # when (for example) our matrix consists of integer entries.
 238     ring = A.base_ring().fraction_field()
 239
 240     if n == 0 or n == 1:
 241         # We can get n == 0 if someone feeds us a trivial matrix.
 242         P = matrix.identity(ring, n)
 243         L = matrix.identity(ring, n)
 244         D = A
 245         return (P,L,D)
 246
 247     alpha = (1 + ZZ(17).sqrt()) * ~ZZ(8)
 248     A1 = A.change_ring(ring)
 249
 250     # Bunch-Kaufmann step 1, Higham step "zero." We use Higham's
 251     # "omega" notation instead of Bunch-Kaufman's "lamda" because
 252     # lambda means other things in the same context.
 253     column_1_subdiag = A1[1:,0].list()
 254     omega_1 = max([ a_i1.abs() for a_i1 in column_1_subdiag ])
 255
 256     if omega_1 == 0:
 257         # "There's nothing to do at this step of the algorithm,"
 258         # which means that our matrix looks like,
 259         #
 260         #   [ 1 0 ]
 261         #   [ 0 B ]
 262         #
 263         # We could still do a pivot_one_by_one() here, but it would
 264         # pointlessly subract a bunch of zeros and multiply by one.
 265         B = A1[1:,1:]
 266         P2, L2, D2 = block_ldlt_naive(B)
 267         P1 = block_matrix(2,2, [[ZZ(1), ZZ(0)],
 268                                 [ZZ(0), P2]])
 269         L1 = block_matrix(2,2, [[ZZ(1), ZZ(0)],
 270                                 [ZZ(0), L2]])
 271         D1 = block_matrix(2,2, [[ZZ(1), ZZ(0)],
 272                                 [ZZ(0), D2]])
 273         return (P1,L1,D1)
 274
 275     def pivot_one_by_one(M, c=None):
 276         # Perform a one-by-one pivot on "M," swapping row/columns "c".
 277         # If "c" is None, no swap is performed.
 278         if c is not None:
 279             P1 = copy(M.matrix_space().identity_matrix())
 280             P1.swap_rows(0,c)
 281             M = P1.T * M * P1
 282
 283         # The top-left entry is now our 1x1 pivot.
 284         C = M[1:n,0]
 285         B = M[1:,1:]
 286
 287         P2, L2, D2 = block_ldlt_naive(B - (C*C.transpose())/M[0,0])
 288
 289         if c is None:
 290             P1 = block_matrix(2,2, [[ZZ(1), ZZ(0)],
 291                                     [0*C,   P2]])
 292         else:
 293             P1 = P1*block_matrix(2,2, [[ZZ(1), ZZ(0)],
 294                                        [0*C,   P2]])
 295
 296         L1 = block_matrix(2,2, [[ZZ(1),                    ZZ(0)],
 297                                 [P2.transpose()*C/M[0,0], L2]])
 298         D1 = block_matrix(2,2, [[M[0,0], ZZ(0)],
 299                                 [0*C,     D2]])
 300
 301         return (P1,L1,D1)
 302
 303
 304     if A1[0,0].abs() > alpha*omega_1:
 305         return pivot_one_by_one(A1)
 306
 307     r = 1 + column_1_subdiag.index(omega_1)
 308
 309     # If the matrix is Hermitian, we need only look at the above-
 310     # diagonal entries to find the off-diagonal of maximal magnitude.
 311     omega_r = max( a_rj.abs() for a_rj in A1[:r,r].list() )
 312
 313     if A1[0,0].abs()*omega_r >= alpha*(omega_1^2):
 314         return pivot_one_by_one(A1)
 315
 316     if A1[r,r].abs() > alpha*omega_r:
 317         # Higham step (3)
 318         # Another 1x1 pivot, but this time swapping indices 0,r.
 319         return pivot_one_by_one(A1,r)
 320
 321     # Higham step (4)
 322     # If we made it here, we have to do a 2x2 pivot.
 323     return None