This is not a functional implementation (as it stands, it's rather MATLAB-ish), but I'll leave this snippet around and hope somebody could make something purely functional out of this outer product form of Cholesky decomposition:

m = Array[Min, {4, 4}]; (* example matrix *)

Do[
  m[[k, k]] = b = Sqrt[a = m[[k, k]]];
  m[[k, k + 1 ;;]] = (v = m[[k, k + 1 ;;]])/b;
  m[[k + 1 ;;, k + 1 ;;]] -= Outer[Times, v, v]/a;
  , {k, Length[m]}];
UpperTriangularize[m]

{{1, 1, 1, 1}, {0, 1, 1, 1}, {0, 0, 1, 1}, {0, 0, 0, 1}}

This is not a functional implementation (as it stands, it's rather MATLAB-ish), but I'll leave this snippet around and hope somebody could make something purely functional out of this outer product form of Cholesky decomposition:

m = Array[Min, {4, 4}]; (* example matrix *)

Do[
  m[[k, k]] = b = Sqrt[a = m[[k, k]]];
  m[[k, k + 1 ;;]] = (v = m[[k, k + 1 ;;]])/b;
  m[[k + 1 ;;, k + 1 ;;]] -= Outer[Times, v, v]/a;
  , {k, Length[m]}];
UpperTriangularize[m]

{{1, 1, 1, 1}, {0, 1, 1, 1}, {0, 0, 1, 1}, {0, 0, 0, 1}}

(added 6/12/2012)

As expected, the method given above can be made purely functional, but to me the functional version looks a lot less elegant:

k = 0; n = Length[m];
MapAt[Sqrt,
 Nest[Function[m, ++k;
   Block[{a = m[[k, k]], v = m[[k, k + 1 ;;]], s},
    s = Sqrt[m[[k, k]]];
    ArrayFlatten[{{ReplacePart[m[[1 ;; k, 1 ;; k]], {k, k} -> s], 
       DiagonalMatrix[Append[ConstantArray[1, k - 1], 1/s]].m[[1 ;; k, k + 1 ;;]]},
                 {0, m[[k + 1 ;;, k + 1 ;;]] - Outer[Times, v, v]/a}}]]],
  m, n - 1],
 {n, n}]