In [1]:
using LinearAlgebra
function cholesky(A::Matrix)
    n = size(A,1)
    L = copy(A)
    d = ones(n)
    D = Diagonal(d)
    Lt = LowerTriangular(L) # This maps just the lower-tri portion
    # get the process started
    for i=1:n-1
        d[i] = L[i,i]
        L[i:n,i] ./= d[i]
        L[i+1:n,i+1:n] -= d[i]*L[i+1:n,i]*L[i+1:n,i]'
    end
    d[n] = L[n,n]
    L[n,n] = 1.0
    return (L=Lt, D=D, M=L)
end
A = [3 -1.0 -1.0; -1 3.0 -1; -1 -1 3.0]
F = cholesky(A)
display(F.L)

@show F.L*F.D*(F.L)' - A

3×3 LowerTriangular{Float64, Matrix{Float64}}:
  1.0         ⋅    ⋅ 
 -0.333333   1.0   ⋅ 
 -0.333333  -0.5  1.0
F.L * F.D * (F.L)' - A = [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0]
Out[1]:
3×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0
In [2]:
A = randn(10,10)
A = A'*A
F = cholesky(A)
@show norm(F.L*F.D*(F.L)' - A)

norm(F.L * F.D * (F.L)' - A) = 6.175100983377661e-15
Out[2]:
6.175100983377661e-15
In [3]:
A = [1 -1.0 -1.0; -1 0 -1; -1 -1 0]
#eigvals(A)
F = cholesky(A)
display(F.L)
display(F.D)
display(F.L*F.D*F.L' - A)

3×3 LowerTriangular{Float64, Matrix{Float64}}:
  1.0   ⋅    ⋅ 
 -1.0  1.0   ⋅ 
 -1.0  2.0  1.0
3×3 Diagonal{Float64, Vector{Float64}}:
 1.0    ⋅    ⋅ 
  ⋅   -1.0   ⋅ 
  ⋅     ⋅   3.0
3×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0
In [1]:
using LinearAlgebra
function modified_cholesky(A::Matrix,delta::Real,beta::Real)
    n = size(A,1)
    L = copy(A)
    d = ones(n)
    D = Diagonal(d)
    Lt = LowerTriangular(L) # This maps just the upper-tri portion
    # get the process started
    for i=1:n-1
        theta = maximum(abs.(L[i+1:n,i]))
        d[i] = max(abs(L[i,i]), (theta/beta)^2, delta)
        #d[i] = 100.0
        L[i:n,i] ./= d[i]
        L[i,i] = 1.0
        L[i+1:n,i+1:n] -= d[i]*L[i+1:n,i]*L[i+1:n,i]'
    end
    d[n] = max(abs.(L[n,n]), delta)
    L[n,n] = 1.0
    return (L=Lt, D=D, M=L)
end
A = [3 -1.0 -1.0; -1 3.0 -1; -1 -1 3.0]
F = modified_cholesky(A, 0.1, 5.0)
display(F.L)
display(F.L*F.D*F.L')
@show F.L*F.D*(F.L)' - A

3×3 LowerTriangular{Float64, Matrix{Float64}}:
  1.0         ⋅    ⋅ 
 -0.333333   1.0   ⋅ 
 -0.333333  -0.5  1.0
3×3 Matrix{Float64}:
  3.0  -1.0  -1.0
 -1.0   3.0  -1.0
 -1.0  -1.0   3.0
F.L * F.D * (F.L)' - A = [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0]
Out[1]:
3×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0
In [2]:
A = randn(10,10)
A = A'*A
F = modified_cholesky(A, 0.1, 5.0)
@show norm(F.L*F.D*(F.L)' - A)

norm(F.L * F.D * (F.L)' - A) = 0.09996486748890554
Out[2]:
0.09996486748890554
In [3]:
A = randn(10,10)
A = A'*A-0.1*I
F = modified_cholesky(A, 0.1, 5.0)
@show norm(F.L*F.D*(F.L)' - A)

norm(F.L * F.D * (F.L)' - A) = 9.304803091120704
Out[3]:
9.304803091120704
In [4]:
display(F.L*F.D*(F.L)' - A)

10×10 Matrix{Float64}:
 0.0  0.0   0.0           0.0          …   0.0           0.0
 0.0  0.0  -4.44089e-16   0.0             -2.22045e-16   0.0
 0.0  0.0   0.0           0.0             -2.77556e-17   0.0
 0.0  0.0   0.0           0.0             -2.22045e-16   0.0
 0.0  0.0   0.0           0.0              8.88178e-16   2.22045e-16
 0.0  0.0   0.0           0.0          …   0.0           0.0
 0.0  0.0   4.44089e-16   0.0             -4.44089e-16  -5.55112e-17
 0.0  0.0   0.0           2.22045e-16     -1.66533e-16  -4.44089e-16
 0.0  0.0   5.55112e-17  -2.22045e-16      4.44089e-16   2.22045e-16
 0.0  0.0  -1.11022e-16   1.11022e-16      0.0           9.3048
In [5]:
Diagonal(F.L*F.D*(F.L)' - A)
Out[5]:
10×10 Diagonal{Float64, Vector{Float64}}:
 0.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅            ⋅            ⋅ 
  ⋅   0.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅            ⋅            ⋅ 
  ⋅    ⋅   0.0   ⋅    ⋅    ⋅    ⋅    ⋅            ⋅            ⋅ 
  ⋅    ⋅    ⋅   0.0   ⋅    ⋅    ⋅    ⋅            ⋅            ⋅ 
  ⋅    ⋅    ⋅    ⋅   0.0   ⋅    ⋅    ⋅            ⋅            ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅   0.0   ⋅    ⋅            ⋅            ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   0.0   ⋅            ⋅            ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   8.88178e-16   ⋅            ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅           4.44089e-16   ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅            ⋅           9.3048
In [ ]: