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.
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=\text{ }$scan($\mathbf{U},$row $i,i+1,n,u$) |
$k=\text{ }$scan($\mathbf{L},$column $i,i+1,n,l$) |
$i=\text{ }$min($j,k$) |
In words, Algorithm 32 starts at u_{ii} (the diagonal element in U corresponding to the nonzero element in y) and looks to the right until it finds the first nonzero element u_{ij}. It then starts at element l_{ii} and looks down column i of L until it finds a nonzero element l_{ik}. The value of i is then reset to the smaller of j or k and the process is repeated for the next u_{ii} and l_{ii}. 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 PAP^{T} has a symmetric sparsity pattern Algorithm 32 simplifies to Algorithm 33.
while $i\neq$finished |
$i=\text{ }$scan($\mathbf{U},$row $i,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.
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 y_{i}.
The z vector is full.
The product yz^{T} has the same sparsity structure as A.
The y and z^{T} vectors have been permuted by P and P^{T} 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$ |
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$ |