Editing Cholesky decomposition (section)

===Updating the decomposition===
A task that often arises in practice is that one needs to update a Cholesky decomposition. In more details, one has already computed the Cholesky decomposition <math display=inline>\mathbf{A} = \mathbf{L}\mathbf{L}^*</math> of some matrix <math display=inline>\mathbf{A}</math>, then one changes the matrix <math display=inline>\mathbf{A}</math> in some way into another matrix, say <math display=inline> \tilde{\mathbf{A}} </math>, and one wants to compute the Cholesky decomposition of the updated matrix: <math display=inline> \tilde{\mathbf{A}} = \tilde{\mathbf{L}} \tilde{\mathbf{L}}^* </math>. The question is now whether one can use the Cholesky decomposition of <math display=inline>\mathbf{A}</math> that was computed before to compute the Cholesky decomposition of <math display=inline> \tilde{\mathbf{A}} </math>.

==== Rank-one update ====
The specific case, where the updated matrix <math display=inline> \tilde{\mathbf{A}} </math> is related to the matrix <math display=inline>\mathbf{A}</math> by <math display=inline> \tilde{\mathbf{A}} = \mathbf{A} + \mathbf{x} \mathbf{x}^* </math>, is known as a ''rank-one update''.

Here is a function<ref>Based on: {{cite book|last=Stewart|first=G. W.|title=Basic decompositions|year=1998|publisher=Soc. for Industrial and Applied Mathematics|location=Philadelphia|isbn=0-89871-414-1}}</ref> written in [[Matlab]] syntax that realizes a rank-one update:
<syntaxhighlight lang="matlab">
function [L] = cholupdate(L, x)
    n = length(x);
    for k = 1:n
        r = sqrt(L(k, k)^2 + x(k)^2);
        c = r / L(k, k);
        s = x(k) / L(k, k);
        L(k, k) = r;
        if k < n
            L((k+1):n, k) = (L((k+1):n, k) + s * x((k+1):n)) / c;
            x((k+1):n) = c * x((k+1):n) - s * L((k+1):n, k);
        end
    end
end
</syntaxhighlight>
A ''rank-n update'' is one where for a matrix <math display=inline>\mathbf{M}</math> one updates the decomposition such that <math display=inline> \tilde{\mathbf{A}} = \mathbf{A} + \mathbf{M} \mathbf{M}^* </math>.  This can be achieved by successively performing rank-one updates for each of the columns of <math display=inline>\mathbf{M}</math>.

==== Rank-one downdate ====
A ''rank-one downdate'' is similar to a rank-one update, except that the addition is replaced by subtraction: <math display=inline> \tilde{\mathbf{A}} = \mathbf{A} - \mathbf{x} \mathbf{x}^* </math>. This only works if the new matrix <math display=inline> \tilde{\mathbf{A}} </math> is still positive definite.

The code for the rank-one update shown above can easily be adapted to do a rank-one downdate: one merely needs to replace the two additions in the assignment to <code>r</code> and <code>L((k+1):n, k)</code> by subtractions.

==== Adding and removing rows and columns ====
If a symmetric and positive definite matrix <math display=inline> \mathbf A </math> is represented in block form as

<math display=block>
\mathbf{A} = 
\begin{pmatrix}
 \mathbf A_{11} & \mathbf A_{13} \\
 \mathbf A_{13}^{\mathrm{T}} & \mathbf A_{33} \\
\end{pmatrix}
</math>

and its upper Cholesky factor
<math display=block>
\mathbf{L} = 
\begin{pmatrix}
 \mathbf L_{11} &  \mathbf L_{13} \\
  0      & \mathbf L_{33} \\
\end{pmatrix},
</math>

then for a new matrix <math display=inline> \tilde{\mathbf{A}} </math>, which is the same as <math display=inline> \mathbf A </math> but with the insertion of new rows and columns,
<math display=block>\begin{align}
\tilde{\mathbf{A}} &= 
\begin{pmatrix}
 \mathbf A_{11} & \mathbf A_{12} & \mathbf A_{13} \\
 \mathbf A_{12}^{\mathrm{T}} & \mathbf A_{22} & \mathbf A_{23} \\
 \mathbf A_{13}^{\mathrm{T}} & \mathbf A_{23}^{\mathrm{T}} &  \mathbf A_{33} \\
\end{pmatrix}
\end{align}
</math>

Now there is an interest in finding the Cholesky factorization of <math display=inline> \tilde{\mathbf{A}} </math>, which can be called <math display=inline> \tilde{\mathbf S} </math>, without directly computing the entire decomposition. 
<math display=block>\begin{align}
\tilde{\mathbf{S}} &= 
\begin{pmatrix}
 \mathbf S_{11} & \mathbf S_{12} & \mathbf S_{13} \\
 0 & \mathbf S_{22} & \mathbf S_{23} \\
 0 & 0 &  \mathbf S_{33} \\
\end{pmatrix}.
\end{align}
</math>

Writing <math display=inline> \mathbf A \setminus \mathbf{b}</math> for the solution of <math display=inline> \mathbf A \mathbf x = \mathbf b</math>, which can be found easily for triangular matrices, and <math display=inline> \text{chol} (\mathbf M)</math> for the Cholesky decomposition of <math display=inline> \mathbf M </math>, the following relations can be found:
<math display=block>\begin{align}
\mathbf S_{11} &= \mathbf L_{11}, \\
\mathbf S_{12} &= \mathbf L_{11}^{\mathrm{T}} \setminus \mathbf A_{12}, \\
\mathbf S_{13} &= \mathbf L_{13}, \\
\mathbf S_{22} &= \mathrm{chol} \left(\mathbf A_{22}  - \mathbf S_{12}^{\mathrm{T}} \mathbf S_{12}\right), \\
\mathbf S_{23} &= \mathbf S_{22}^{\mathrm{T}} \setminus \left(\mathbf A_{23} - \mathbf S_{12}^{\mathrm{T}} \mathbf S_{13}\right), \\
\mathbf S_{33} &= \mathrm{chol} \left(\mathbf L_{33}^{\mathrm{T}} \mathbf L_{33} - \mathbf S_{23}^{\mathrm{T}} \mathbf S_{23}\right).
\end{align}
</math>

These formulas may be used to determine the Cholesky factor after the insertion of rows or columns in any position, if the row and column dimensions are appropriately set (including to zero). The inverse problem, 

<math display=block>\begin{align}
\tilde{\mathbf{A}} &= 
\begin{pmatrix}
 \mathbf A_{11} & \mathbf A_{12} & \mathbf A_{13} \\
 \mathbf A_{12}^{\mathrm{T}} & \mathbf A_{22} & \mathbf A_{23} \\
 \mathbf A_{13}^{\mathrm{T}} & \mathbf A_{23}^{\mathrm{T}} &  \mathbf A_{33} \\
\end{pmatrix}
\end{align}
</math>
with known Cholesky decomposition 
<math display=block>\begin{align}
\tilde{\mathbf{S}} &= 
\begin{pmatrix}
 \mathbf S_{11} & \mathbf S_{12} & \mathbf S_{13} \\
 0 & \mathbf S_{22} & \mathbf S_{23} \\
 0 & 0 &  \mathbf S_{33} \\
\end{pmatrix}
\end{align}
</math>

and the desire to determine the Cholesky factor
<math display=block>\begin{align}
\mathbf{L} &= 
\begin{pmatrix}
 \mathbf L_{11} &  \mathbf L_{13} \\
  0      & \mathbf L_{33} \\
\end{pmatrix}
\end{align}
</math>

of the matrix <math display=inline> \mathbf A </math> with rows and columns removed,
<math display=block>\begin{align}
\mathbf{A} &= 
\begin{pmatrix}
 \mathbf A_{11} & \mathbf A_{13} \\
 \mathbf A_{13}^{\mathrm{T}} & \mathbf A_{33} \\
\end{pmatrix},
\end{align}
</math>

yields the following rules:
<math display=block>\begin{align}
\mathbf L_{11} &= \mathbf S_{11}, \\
\mathbf L_{13} &= \mathbf S_{13}, \\
\mathbf L_{33} &= \mathrm{chol} \left(\mathbf S_{33}^{\mathrm{T}} \mathbf S_{33} + \mathbf S_{23}^{\mathrm{T}} \mathbf S_{23}\right).
\end{align}
</math>

Notice that the equations above that involve finding the Cholesky decomposition of a new matrix are all of the form <math display=inline> \tilde{\mathbf{A}} = \mathbf{A} \pm \mathbf{x} \mathbf{x}^* </math>, which allows them to be efficiently calculated using the update and downdate procedures detailed in the previous section.<ref>Osborne, M. (2010), Appendix B.</ref>