Cholesky updates of A’A

Updating the Cholesky factorization of A'A=LL' when one or more columns are added to or removed from matrix A can be done very efficiently obviating the re-factorization from scratch.

Removing columns

Assume we know the Cholesky factorisation of A'A=LL' and we remove a column from matrix A which can be written as follows

A = \begin{bmatrix} A_{1:\kappa} & c & A_{\kappa+1:n} \end{bmatrix}

Then,

\begin{aligned} A'A &= \begin{bmatrix} A_{1:\kappa} & c & A_{\kappa+1:n} \end{bmatrix}'\begin{bmatrix} A_{1:\kappa} & c & A_{\kappa+1:n} \end{bmatrix}\\ &=\begin{bmatrix} A_{1:\kappa}'A_{1:\kappa} & A_{1:\kappa}'c & A_{1:\kappa}'A_{\kappa+1:n}\\ c'A_{1:\kappa} & c'c & c'A_{\kappa+1:n}\\ A_{\kappa+1:n}'A_{1:\kappa} & A_{1:\kappa}c & A_{\kappa+1:n}'A_{\kappa+1:n} \end{bmatrix}\\ &=LL'\\ &=\begin{bmatrix} L_{11}\\ \ell_{12} & \ell_{22}\\ L_{31} & \ell_{32} & L_{33} \end{bmatrix} \begin{bmatrix} L_{11}\\ \ell_{12} & \ell_{22}\\ L_{31} & \ell_{32} & L_{33} \end{bmatrix}' \end{aligned}

Once we delete the column c at \kappa+1 we have

\begin{bmatrix} A_{1:\kappa}'A_{1:\kappa} & A_{1:\kappa}'A_{\kappa+1:n}\\ A_{\kappa+1:n}'A_{1:\kappa} & A_{\kappa+1:n}'A_{\kappa+1:n} \end{bmatrix}= \begin{bmatrix} \bar{L}_{11}\\ \bar{L}_{31} & \bar{L}_{33} \end{bmatrix} \begin{bmatrix} \bar{L}_{11}\\ \bar{L}_{31} & \bar{L}_{33} \end{bmatrix}'

from which we have that

\begin{aligned} \bar{L}_{11} &= L_{11}\\ \bar{L}_{31} &= L_{31} \end{aligned}

and

\bar{L}_{33}\bar{L}_{33}' = L_{33}L_{33}' + \ell_{32}\ell_{32}'.

The last equation is a rank-1 update (see Golub and Van Loan) and it comes at the
cost of 2(n-\kappa)^2+4(n-\kappa).

Appending columns using permutation matrix

Having computed a Cholesky factorisation with permutation for the matrix A'A, that is A'A = PLL'P', we need to compute the Cholesky factorisation of \bar{A}'\bar{A}, where \bar{A}=[A\ c].

We have

\begin{aligned} \bar{A}'\bar{A} &= \begin{bmatrix}A & c\end{bmatrix}'\begin{bmatrix}A & c\end{bmatrix}\\ &= \begin{bmatrix}A'A & A'c\\c'A & c'c\end{bmatrix}\\ &= \begin{bmatrix}P \\ & 1\end{bmatrix} \begin{bmatrix}L_1\\\ell_1' & \ell_2\end{bmatrix} \begin{bmatrix}L_1\\\ell_1' & \ell_2\end{bmatrix}' \begin{bmatrix}P \\ & 1\end{bmatrix}' \end{aligned}

From which we have that

A'A = PL_1L_1'P' \Rightarrow L_1 = L,

and $\ell_1$ is computed from the following nice linear system

PL\ell_1 = A'c \Leftrightarrow L\ell_1 = P'A'c,

and, provided that c'c - \ell_1'\ell_1>0, \bar{A}'\bar{A} is positive definite and

\ell_2^2 = c'c - \ell_1'\ell_1.

Remark. to solve A'Ax=b we do PLL'P'x= b and we set x=Py so we then need to solve

PLL'y = b \Leftrightarrow LL'y = P'b.

The permuted Cholesky factor of \bar{A}'\bar{A}=\bar{P}\bar{L}\bar{L}'\bar{P}' is

\bar{L}=\begin{bmatrix}L_1\\\ell_1' & \ell_2\end{bmatrix},

and the corresponding permutation matrix is

\bar{P}=\begin{bmatrix}P \\ & 1\end{bmatrix}.

The computational cost of this update is n^2+3n.

Inserting columns

In this example we shall insert a column in A, so we shall define the matrix

\tilde{A} = \begin{bmatrix} A_{1:\kappa} & c & A_{\kappa+1:n} \end{bmatrix},

where \kappa\in\mathbb{N}_{[1,n]}.

There is then a permutation matrix \tilde{P} so that \tilde{A}\tilde{P}' = [A\ c].

It is then easy to update the factorisation
\begin{aligned} &(\tilde{A}\tilde{P})'(\tilde{A}\tilde{P})=\tilde{L}\tilde{L}', \\ \Leftrightarrow &\ \tilde{A}'\tilde{A} = \tilde{P}\tilde{L}\tilde{L}'\tilde{P}' \end{aligned}

which is the updated Cholesky factorisation of \tilde{A}'\tilde{A} with permutation matrix \tilde{P}.

Inserting many columns

Now we are going to insert many columns, recursively, at various positions in \tilde{A}. Essentially, we are going to update a permuted Cholesky factorisation by inserting a column in any position in A, so the updated matrix becomes

\bar{A} = \begin{bmatrix} A_{1:\kappa} & c & A_{\kappa+1:n} \end{bmatrix}

There is a permutation matrix \tilde{P} so that

\bar{A}\tilde{P} = \begin{bmatrix} A & c \end{bmatrix} \triangleq \tilde{A}

Say we have A'A=PLL'P'. We can then compute a permuted Cholesky factorisation for \tilde{A}'\tilde{A},

\begin{aligned} &\tilde{A}'\tilde{A}= \bar{P}\tilde{L}\tilde{L}'\bar{P}'\\ \Leftrightarrow &\ \tilde{P}'\bar{A}'\bar{A}\tilde{P} = \bar{P}\tilde{L}\tilde{L}'\bar{P}'\\ \Leftrightarrow &\ \bar{A}'\bar{A} = \tilde{P}\bar{P}\tilde{L}\tilde{L}'\bar{P}'\tilde{P}', \end{aligned}
which is the permuted Cholesky factorisation of \bar{A}'\bar{A} with permutation matrix \tilde{P}\bar{P}.

Let now A\in\Re^{m\times n} be a given matrix and \alpha be a collection of column indexes of A, and A_\alpha\in\Re^{m\times |\alpha|} and we know the Cholesky factorisation of A_{\alpha}'A_{\alpha}=L_{\alpha}L_{\alpha}'.

Let \bar{\alpha} = \alpha\cup \{\alpha^*\}, i.e., the \alpha^*-th column of A is added to A_{\alpha} to form the new matrix

A_{\bar\alpha} = \begin{bmatrix} A_{\alpha} & A_{\alpha^*} \end{bmatrix}.

Having augmented A_{\alpha} we can now update the factorisation of A_{\alpha}'A_{\alpha} and compute an L_{\bar\alpha} such that

A_{\bar\alpha}'A_{\bar\alpha} = L_{\bar\alpha}L_{\bar\alpha}'.

 

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

mathbabe

Exploring and venting about quantitative issues

Look at the corners!

The math blog of Dmitry Ostrovsky

The Unapologetic Mathematician

Mathematics for the interested outsider

Almost Sure

A random mathematical blog

Mathematix

Mathematix is the mathematician of the village of Asterix

%d bloggers like this: