6.2 LU Factor Update

Similar algorithms exist for updating LU decompositions of A. If U is upper triangular and L is unit lower triangular, an element uij from U is related to the elements uiju^{\prime}_{ij} and did^{\prime}_{i} of the LDU decomposition as follows.

uij=uijdiu_{ij}=u^{\prime}_{ij}d^{\prime}_{i} (50)

The recurrence relations of the inner loop of Algorithm 8 must change to reflect this relationship. The following statement updates the z vector for a unit upper triangular U.

zj=zj-quijz_{j}=z_{j}-q{u}^{\prime}_{ij} (51)

If U is upper triangular, the statement becomes

zj=zj-puijδz_{j}=z_{j}-p\frac{u_{ij}}{\delta} (52)

where δ\delta is the value of uii before it was changed during stage i of the procedure. Along the same lines, the factor update statement

uij=uij+β2zju^{\prime}_{ij}=u^{\prime}_{ij}+{\beta}_{2}{z}_{j} (53)

becomes

uijuii=uijδ+β2zj\frac{u_{ij}}{u_{ii}}=\frac{{u}_{ij}}{\delta}+{\beta}_{2}{z}_{j} (54)

Solving for the updated value of uij yields

uij=uii(uijδ+β2zj)u_{ij}=u_{ii}(\frac{{u}_{ij}}{\delta}+{\beta}_{2}{z}_{j}) (55)

Taking these observations into consideration and pulling operations on constants out of the inner loop, Algorithm 9 updates U based on a rank one change to A.

Algorithm 9: LU Factor Update
for i=1,,ni=1,\cdots,n
    δ=uii\delta=u_{ii}
    p=yip=y_{i}
    q=ziq=z_{i}
    uii=uii+αpqu_{ii}=u_{ii}+\alpha pq
    β1=αuii\beta_{1}=\displaystyle\frac{\alpha}{u_{ii}}
    β2=αq\beta_{2}=\alpha q
    q=qδq=\displaystyle\frac{q}{\delta}
    δ=uiiδ\delta=\displaystyle\frac{u_{ii}}{\delta}
    α=αδ\alpha=\displaystyle\frac{\alpha}{\delta}
    for j=i+1,,nj=i+1,\cdots,n
        yj=yj-pljiy_{j}=y_{j}-pl_{ji}
        zj=zj-quijz_{j}=z_{j}-qu_{ij}
        lji=lji+β1yjl_{ji}=l_{ji}+\beta_{1}y_{j}
        uij=δuij+β2zju_{ij}=\delta u_{ij}+\beta_{2}z_{j}

If U is unit upper triangular and L is lower triangular, a similar algorithm is derived from the observation that lij of L and lijl^{\prime}_{ij}, did^{\prime}_{i} of the LDU decomposition are related as follows.

lij=lijdjl_{ij}=l^{\prime}_{ij}d^{\prime}_{j} (56)

The resulting algorithm deviates from Algorithm 8 in a manner that parallels Algorithm 9.