Cholesky factorization

using LinearAlgebra

A randomly chosen matrix is extremely unlikely to be symmetric. However, there is a simple way to symmetrize one.

A = rand(1.0:9.0,4,4)
B = A + A'
4×4 Array{Float64,2}:
 12.0  10.0  11.0  12.0
 10.0   4.0   9.0  10.0
 11.0   9.0  12.0   5.0
 12.0  10.0   5.0   4.0

Similarly, a random symmetric matrix is unlikely to be positive definite. The Cholesky algorithm always detects a non-PD matrix by quitting with an error.

# cholesky(B)   throws a PosDefException

It’s not hard to manufacture an SPD matrix to try out the Cholesky factorization.

B = A'*A
cf = cholesky(B)
Cholesky{Float64,Array{Float64,2}}
U factor:
4×4 UpperTriangular{Float64,Array{Float64,2}}:
 13.1909  5.68574  11.3715     5.68574
   ⋅      6.37749  -0.416335   0.419038
   ⋅       ⋅        1.58629   -4.08545
   ⋅       ⋅         ⋅         2.96748

What’s returned is a “factorization object.” (This allows it to be used efficiently in various contexts.) Another step is required to extract the factor as a matrix.

R = Matrix(cf.U)
4×4 Array{Float64,2}:
 13.1909  5.68574  11.3715     5.68574
  0.0     6.37749  -0.416335   0.419038
  0.0     0.0       1.58629   -4.08545
  0.0     0.0       0.0        2.96748
opnorm(R'*R - B)
0.0