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