7 Symmetric Matrices 7.7 Backward Substitution for Symmetric Systems 8 Sparse Matrices

7.8 Symmetric Factor Update

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.

7.8.1 Symmetric LDU Factor Update

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.

Algorithm 20: Symmetric LDU Factor Update

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.

7.8.2 Symmetric LU Factor Update

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.

Algorithm 21: Symmetric LU Factor Update

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 87.
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.