8.8 Sparse LU Factor Update

If the LU decomposition of the sparse product PAQ exists and the factorization of a related matrix

𝐏𝐀𝐐=𝐏(𝐀+𝚫𝐀)𝐐\mathbf{PA^{\prime}Q=P(A+\Delta A)Q} (76)

is desired, factor update is at times the procedure of choice. The current section examines procedures for updating the sparse factorization of PAQ following a rank one modification of A, that is

𝐀=𝐀+α𝐲𝐳𝐓\mathbf{A^{\prime}=A}+\alpha\mathbf{y{z}^{T}} (77)

where α\alpha is a scalar, y and z are n vectors, and A{}^{\prime} has the same sparsity pattern as A.

The condition on the structure of A{}^{\prime} is not imposed by the factor update process, but is instead a comment on the utility of factor update in a sparse environment. If the modification to A introduces new elements into the matrix, the pivot sequence determined during symbolic factorization may no longer apply. The sparsity degradation introduced by an inappropriate pivot sequence may outweigh the benefits gained from the updating the existing factorization.

The performance of factor update algorithms is often enhanced by restricting pivot operations to the portions of L and U that are directly effected by the change in A. Papers by Tinney et al. [16] and Chan and Brandwajn [2] describe a systematic methodology for determining this subset of LU. The rows of U and columns in L that are changed during factor update are referred to as the factorization path. The fundamental operation is to determine the factorization path associated with a vector y with just one nonzero element. Such a vector is called a singleton. Its factorization path is called a singleton path. If more than one of the elements in y are nonzero, the composite factorization path is simply the union of the singleton paths.

8.8.1 Factorization Path of a Singleton Update

If the LU decomposition of PAQ is stored as an asymmetric sparse matrix and a singleton vector y has one nonzero element at location i, its factorization path through LU is determined by Algorithm 32. Each value ascribed to i by Algorithm 32 is a vertex on the factorization path of y.

Algorithm 32: Factorization Path
while ii\neqfinished
   j= j=\text{ }scan(𝐔,\mathbf{U},row i,i+1,n,ui,i+1,n,u)
   k= k=\text{ }scan(𝐋,\mathbf{L},column i,i+1,n,li,i+1,n,l)
   i= i=\text{ }min(j,kj,k)

In words, Algorithm 32 starts at uii (the diagonal element in U corresponding to the nonzero element in y) and looks to the right until it finds the first nonzero element uij. It then starts at element lii and looks down column i of L until it finds a nonzero element lik. The value of i is then reset to the smaller of j or k and the process is repeated for the next uii and lii. The procedure ends when there are no elements to the right of the diagonal in U or below the diagonal in L for some vertex on the factorization path. Obviously, it is assumes that a column scan is independent of a row scan but works in a similar manner.

If the LU decomposition of PAPT has a symmetric sparsity pattern Algorithm 32 simplifies to Algorithm 33.

Algorithm 33: Symmetric Factorization Path
while ii\neqfinished
   i= i=\text{ }scan(𝐔,\mathbf{U},row i,i+1,n,ui,i+1,n,u)

The symmetry makes the search through column i of L unnecessary. In Algorithm 33, the use of U to determine the factorization path was arbitrary. Either L or U can anchor the process.

8.8.2 Revising LU after a Singleton Update

Algorithm 34 updates a structurally symmetric LU factorization after a rank one modification to A. It assumes:

  • L is unit lower triangular and maintained in a sparse data structure that communicates using the buffer l.

  • U is upper triangular and maintained in a sparse data structure that communicates using the buffer u.

  • The sparse data structure does not permit a column scan.

  • The y vector is full and all entries are zero except yi.

  • The z vector is full.

  • The product yzT has the same sparsity structure as A.

  • The y and zT vectors have been permuted by P and PT so they are in the same frame of reference as L and U.

Algorithm 34: Structurally Symmetric Sparse LU Factor Update
while ii\neqfinished
   get(𝐔,i,i,u\mathbf{U},i,i,u)
   δ=u\delta=u
   p=yip=y_{i}
   q=ziq=z_{i}
   u=u+αpqu=u+\alpha pq
   put(𝐔,i,i,u\mathbf{U},i,i,u)
   β1=αpu\beta_{1}=\displaystyle\frac{\alpha p}{u}
   β2=αq\beta_{2}=\alpha q
   q=qδq=\displaystyle\frac{q}{\delta}
   δ=uδ\delta=\displaystyle\frac{u}{\delta}
   α=αδ\alpha=\displaystyle\frac{\alpha}{\delta}
   k=k=finished
   while[j=j=scan(𝐔,\mathbf{U},row i,i+1,n,ui,i+1,n,u)] \neq finished
      if jj is first element on row ii
         k=jk=j
      zj=zj-quz_{j}=z_{j}-qu
      u=δu+β2zju=\delta u+\beta_{2}z_{j}
      put(𝐔,i,j,u\mathbf{U},i,j,u)
      get(𝐋,j,i,l\mathbf{L},j,i,l)
      yj=yj-ply_{j}=y_{j}-pl
      l=l+β1yjl=l+\beta_{1}y_{j}
      put(𝐋,j,i,l\mathbf{L},j,i,l)
   i=ki=k

The requirement that L and U occupy different data structures in Algorithm 34 is pedantic. In practice, L will occupy the subdiagonal portion of a sparse matrix and U will occupy the diagonal and superdiagonal portions of the same matrix.

If A is symmetric, y = z, and A{}^{\prime} has the same sparsity pattern as A, Algorithm 34 simplifies to Algorithm 35.

Algorithm 35: Symmetric Sparse LU Factor Update
while ii\neqfinished
   get(𝐔,i,i,u\mathbf{U},i,i,u)
   δ=u\delta=u
   p=yip=y_{i}
   u=u+αp2u=u+\alpha p^{2}
   put(𝐔,i,i,u\mathbf{U},i,i,u)
   β=αpu\beta=\displaystyle\frac{\alpha p}{u}
   q=qδq=\displaystyle\frac{q}{\delta}
   δ=uδ\delta=\displaystyle\frac{u}{\delta}
   α=αδ\alpha=\displaystyle\frac{\alpha}{\delta}
   k=k=finished
   while[j=j=scan(𝐔,\mathbf{U},row i,i+1,n,ui,i+1,n,u)] \neq finished
      if jj is first element on row ii
         k=jk=j
      yj=yj-puy_{j}=y_{j}-pu
      u=δu+βzju=\delta u+\beta z_{j}
      put(𝐔,i,j,u\mathbf{U},i,j,u)
   i=ki=k