= Solution
The following <R (programming language)> code implements the given <recursion> for the <Cholesky decomposition>. The input is assumed to be a real <symmetric matrix> that is a <positive-definite matrix>.
``
mychol <- function(A) {
n <- nrow(A)
if (n == 1L) return(matrix(sqrt(A[1L, 1L]), 1L, 1L))
j <- 2L:n
L <- matrix(0, n, n)
L[1L, 1L] <- sqrt(A[1L, 1L])
b <- A[j, 1L, drop = FALSE] / L[1L, 1L]
L[j, 1L] <- b
S <- A[j, j, drop = FALSE] - tcrossprod(b)
L[j, j] <- mychol(S)
L
}
``
Here `tcrossprod(b)` is the <outer product> $bb^T$. The `drop = FALSE` arguments preserve <matrix> dimensions, including the final one-by-one subproblem, and the <base case> also returns a <matrix>.
To verify correctness, partition $A=\begin{pmatrix}a&r^T\\r&B\end{pmatrix}$. The code sets $\ell=\sqrt a$, $b=r/\ell$ and recursively factors $S=B-bb^T$. This <Schur complement> is a <positive-definite matrix>: for $z\ne0$,
$$
z^TSz=
\begin{pmatrix}-r^Tz/a\\z\end{pmatrix}^{\!T}
A
\begin{pmatrix}-r^Tz/a\\z\end{pmatrix}>0.
$$
If $CC^T=S$ is the recursive answer, then the returned <lower triangular matrix> satisfies
$$
L=\begin{pmatrix}\ell&0\\b&C\end{pmatrix},
\qquad
LL^T=\begin{pmatrix}a&\ell b^T\\\ell b&bb^T+CC^T\end{pmatrix}
=\begin{pmatrix}a&r^T\\r&B\end{pmatrix}=A.
$$
The one-dimensional identity starts this <mathematical induction> proof. Each recursive call reduces the <matrix> size by one and therefore terminates.
Back to article page
