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,,nj=1,\cdots,n
   δ=dj\delta=d_{j}
   p=yjp=y_{j}
   dj=dj+αp2d_{j}=d_{j}+\alpha p^{2}
   β=αpdj\beta=\displaystyle\frac{\alpha p}{d_{j}}
   α=αδdj\alpha=\displaystyle\frac{\alpha\delta}{d_{j}}
   for i=j+1,,ni=j+1,\cdots,n
      yi=yi-plijy_{i}=y_{i}-pl_{ij}
      lij=lij+βyil_{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,,ni=1,\cdots,n
   δ=uii\delta=u_{ii}
   p=yip=y_{i}
   uii=uii+αp2u_{ii}=u_{ii}+\alpha p^{2}
   β=αp\beta=\alpha p
   p=pδp=\displaystyle\frac{p}{\delta}
   δ=uiiδ\delta=\displaystyle\frac{u_{ii}}{\delta}
   α=αδ\alpha=\displaystyle\frac{\alpha}{\delta}
   for j=i+1,,nj=i+1,\cdots,n
      yj=yj-puijy_{j}=y_{j}-pu_{ij}
      uij=δuij+βyyu_{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

  1. The double subscript notation is replaced by the indexing rule defined in Equation 71.

  2. The outer loop counter i ranges from zero to n-1.

  3. 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 yi 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 yi 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 yi is encountered.

For a fuller discussion of the derivation and implementation of LU factor update, see Section 6.2.