In [1]:
using LinearAlgebra
function cholesky(A::Matrix; throw_indefinite_error::Bool = false)
    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]
        if throw_indefinite_error && d[i] < 0 
            throw(ArgumentError("A is not positive definite"))
        end 
        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
    
    if throw_indefinite_error && d[n] < 0 
        throw(ArgumentError("A is not positive definite"))
    end 
    
    return (L=Lt, D=D, M=L)
end
A = [3 -1.0 -1.0; -1 3.0 -1; -1 -1 3.0]
display(A)
F = cholesky(A)
display(F.L)

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

3×3 Matrix{Float64}:
  3.0  -1.0  -1.0
 -1.0   3.0  -1.0
 -1.0  -1.0   3.0
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) = 4.252041856088627e-15
Out[2]:
4.252041856088627e-15
In [3]:
A = [1 -1.0 -1.0; -1 0 -1; -1 -1 0]
@show eigvals(A)
F = cholesky(A)
display(F.L)
display(F.D)
display(F.L*F.D*F.L' - A)

eigvals(A) = [-1.7320508075688772, 1.0000000000000002, 1.7320508075688774]

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 [4]:
cholesky(A; throw_indefinite_error=true)

ArgumentError: A is not positive definite

Stacktrace:
 [1] cholesky(A::Matrix{Float64}; throw_indefinite_error::Bool)
   @ Main ./In[1]:12
 [2] top-level scope
   @ In[4]:1

Bisection to perturb¶

In [5]:
function is_posdef(A)
    try
        F = cholesky(A; throw_indefinite_error=true)
    catch ArgumentError
        return false
    end
    return true
end 
is_posdef(A)
Out[5]:
false
In [6]:
is_posdef(rand(10,10) |> A->A'*A)
Out[6]:
true
In [7]:
""" 
    bisection_to_positive_definite(A)

Return a value tau such that A + tau*I is sym. pos. def. 
"""
function bisection_to_positive_definite(A; taumax=1e9)
    # check the input 
    tau = 1.0 
    if is_posdef(A)
        return 0.0 
    else
        while is_posdef(A + tau*I) == false 
            tau = 2*tau 
            if tau > taumax 
                throw(ArgumentError("the matrix has a strong negative direction"))
            end 
        end 
        return tau 
    end
end 
bisection_to_positive_definite(A)
Out[7]:
2.0
In [8]:
eigvals(A + 2*I)
Out[8]:
3-element Vector{Float64}:
 0.2679491924311227
 2.9999999999999996
 3.7320508075688776

Modified Cholesky¶

In [9]:
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)
        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]
@show eigvals(A)
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

eigvals(A) = [0.9999999999999991, 4.0, 4.0]

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[9]:
3×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0
In [29]:
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) = 6.282611876161638e-15
Out[29]:
6.282611876161638e-15
In [34]:
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) = 5.81109460383906
Out[34]:
5.81109460383906
In [35]:
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           0.0           0.0              0.0          -2.22045e-16
 0.0   0.0          -8.88178e-16  -1.11022e-16     -2.22045e-16   0.0
 0.0   0.0           0.0           0.0             -2.22045e-16  -5.55112e-17
 0.0   0.0           0.0          -4.44089e-16      4.44089e-16   1.249e-16
 0.0   1.11022e-16   5.55112e-17   8.88178e-16  …  -3.33067e-16   0.0
 0.0   0.0           0.0          -4.44089e-16      2.22045e-16  -1.11022e-16
 0.0   0.0           0.0          -4.44089e-16      4.44089e-16   0.0
 0.0  -2.22045e-16  -5.55112e-17  -2.22045e-16      8.88178e-16  -6.66134e-16
 0.0   0.0           0.0           1.66533e-16     -6.66134e-16   5.81109
In [36]:
Diagonal(F.L*F.D*(F.L)' - A)
Out[36]:
10×10 Diagonal{Float64, Vector{Float64}}:
 0.0   ⋅     ⋅            ⋅    ⋅   …   ⋅            ⋅            ⋅ 
  ⋅   0.0    ⋅            ⋅    ⋅       ⋅            ⋅            ⋅ 
  ⋅    ⋅   -8.88178e-16   ⋅    ⋅       ⋅            ⋅            ⋅ 
  ⋅    ⋅     ⋅           0.0   ⋅       ⋅            ⋅            ⋅ 
  ⋅    ⋅     ⋅            ⋅   0.0      ⋅            ⋅            ⋅ 
  ⋅    ⋅     ⋅            ⋅    ⋅   …   ⋅            ⋅            ⋅ 
  ⋅    ⋅     ⋅            ⋅    ⋅       ⋅            ⋅            ⋅ 
  ⋅    ⋅     ⋅            ⋅    ⋅      8.88178e-16   ⋅            ⋅ 
  ⋅    ⋅     ⋅            ⋅    ⋅       ⋅           8.88178e-16   ⋅ 
  ⋅    ⋅     ⋅            ⋅    ⋅       ⋅            ⋅           5.81109
In [ ]: