If the LU decomposition of the n × n symmetric matrix A exists and the factorization of a related matrix
$\mathbf{A^{\prime}=A+\Delta A}$ |
is desired, factor update is often the procedure of choice.
Section 6 examines factor update techniques for dense, asymmetric matrices. The current section examines techniques that exploit computational efficiencies introduced by symmetry. Symmetry reduces the work required to update the factorization of A by half, just as it reduces the work required to decompose A in the first place.
More specifically, the current section examines procedures for updating the factors of A following a symmetric rank one modification
$\mathbf{A^{\prime}=A}+\alpha\mathbf{y{y}^{T}}$ |
where $\alpha$ is a scalar and y is an n vector.
Algorithm C1 of Gill et al. [7] updates the factors L and D of a LDU decomposition of A. The algorithm assumes that the upper triangular factors U are implicit. Algorithm 20 mimics Algorithm C1. The scalar $\alpha$ and the y vector are destroyed by this procedure. The factors of A are overwritten by their new values.
for $j=1,\cdots,n$ |
$\delta=d_{j}$ |
$p=y_{j}$ |
$d_{j}=d_{j}+\alpha p^{2}$ |
$\beta=\displaystyle\frac{\alpha p}{d_{j}}$ |
$\alpha=\displaystyle\frac{\alpha\delta}{d_{j}}$ |
for $i=j+1,\cdots,n$ |
$y_{i}=y_{i}-pl_{ij}$ |
$l_{ij}=l_{ij}+\beta y_{i}$ |
See Section 6.1 for a discussion of updating asymmetric LDU factorizations following rank one changes to a matrix.
If the LU decomposition of the symmetric matrix A exists and U is stored explicitly, the recurrence relations of the inner loop of Algorithm 20 must change. Following a train of thought similar to the derivation of Algorithm 9 (see Section 6.2 for details) results in Algorithm 21 which updates U based on a rank one change to A.
for $i=1,\cdots,n$ |
$\delta=u_{ii}$ |
$p=y_{i}$ |
$u_{ii}=u_{ii}+\alpha p^{2}$ |
$\beta=\alpha p$ |
$p=\displaystyle\frac{p}{\delta}$ |
$\delta=\displaystyle\frac{u_{ii}}{\delta}$ |
$\alpha=\displaystyle\frac{\alpha}{\delta}$ |
for $j=i+1,\cdots,n$ |
$y_{j}=y_{j}-pu_{ij}$ |
$u_{ij}=\delta u_{ij}+\beta y_{y}$ |
If U is maintained in a zero-based linear array, Algorithm 21 changes in the normal manner, that is
The double subscript notation is replaced by the indexing rule defined in Equation 71.
The outer loop counter i ranges from zero to n-1.
The inner loop counter j ranges from i+1 to n-1.
The outer loops of the symmetric factor update algorithms do not have to begin at one unless y is full. If U has any leading zeros, the initial value of i (or j in Algorithm 20) should be the index of y_{i} the first nonzero element of y. If there is no a priori information about the structure of y but there is a high probability of leading zeros, testing y_{i} for zero at the beginning of the loop might save a lot of work. However, you must remember to suspend the test as soon as the first nonzero value of y_{i} is encountered.
For a fuller discussion of the derivation and implementation of LU factor update, see Section 6.2.