Solution

ID: past-exam-of-the-simulation-and-statistical-programming-course-of-the-university-of-oxford/2023/3/a/solution

The following R 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 . 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 . The code sets , and recursively factors . This Schur complement is a positive-definite matrix: for ,
If is the recursive answer, then the returned lower triangular matrix satisfies
The one-dimensional identity starts this mathematical induction proof. Each recursive call reduces the matrix size by one and therefore terminates.

New to topics? Read the docs here!