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}$$ | (94) |
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}}$$ | (95) |
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 directlyb effected by the
change in A. Papers by
Tinney
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.
| while $i\neq$ finished |
| $j=$ scan( $\mathbf{U},$row $i,i+1,n,u$ ) |
| $k=$ scan( $\mathbf{L},$column $i,i+1,n,l$ ) |
| $i=$ min( $j,k$ ) |
| $^\dagger$ Operations and data types are defined in Section 8.2.1. |
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.
| while $i\neq$ finished |
| $i=$ scan( $\mathbf{U},$row $i,i+1,n,u$ ) |
| $^\dagger$ Operations and data types are defined in Section 8.2.1. |
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.
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.
| while $i\neq$ finished |
| get( $\mathbf{U},i,i,u$ ) |
| $\delta=u$ |
| $p=y_{i}$ |
| $q=z_{i}$ |
| $u=u+\alpha pq$ |
| put( $\mathbf{U},i,i,u$ ) |
| $\beta_{1}=\displaystyle\frac{\alpha p}{u}$ |
| $\beta_{2}=\alpha q$ |
| $q=\displaystyle\frac{q}{\delta}$ |
| $\delta=\displaystyle\frac{u}{\delta}$ |
| $\alpha=\displaystyle\frac{\alpha}{\delta}$ |
| $k=$ finished |
| while [ $j=$ scan( $\mathbf{U},$row $i,i+1,n,u$ ) ] $\neq$ finished |
| if $j$ is first element on row $i$ |
| $k=j$ |
| $z_{j}=z_{j}-qu$ |
| $u=\delta u+\beta_{2}z_{j}$ |
| put( $\mathbf{U},i,j,u$ ) |
| get( $\mathbf{L},j,i,l$ ) |
| $y_{j}=y_{j}-pl$ |
| $l=l+\beta_{1}y_{j}$ |
| put( $\mathbf{L},j,i,l$ ) |
| $i=k$ |
| $^\dagger$ Operations and data types are defined in Section 8.2.1. |
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.
| while $i\neq$ finished |
| get( $\mathbf{U},i,i,u$ ) |
| $\delta=u$ |
| $p=y_{i}$ |
| $u=u+\alpha p^{2}$ |
| put( $\mathbf{U},i,i,u$ ) |
| $\beta=\displaystyle\frac{\alpha p}{u}$ |
| $q=\displaystyle\frac{q}{\delta}$ |
| $\delta=\displaystyle\frac{u}{\delta}$ |
| $\alpha=\displaystyle\frac{\alpha}{\delta}$ |
| $k=$ finished |
| while [ $j=$ scan( $\mathbf{U},$row $i,i+1,n,u$ ) ] $\neq$ finished |
| if $j$ is first element on row $i$ |
| $k=j$ |
| $y_{j}=y_{j}-pu$ |
| $u=\delta u+\beta z_{j}$ |
| put( $\mathbf{U},i,j,u$ ) |
| $i=k$ |
| $^\dagger$ Operations and data types are defined in Section 8.2.1. |