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]
3×3 Matrix{Float64}:
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
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) = 3.4582155199852684e-15
3.4582155199852684e-15
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
cholesky(A; throw_indefinite_error=true)
ArgumentError: A is not positive definite
Stacktrace:
[1] cholesky(A::Matrix{Float64}; throw_indefinite_error::Bool)
@ Main ./In[25]:12
[2] top-level scope
@ In[32]:1
function is_posdef(A)
try
F = cholesky(A; throw_indefinite_error=true)
catch ArgumentError
return false
end
return true
end
is_posdef(A)
false
is_posdef(rand(10,10) |> A->A'*A)
true
"""
bisection_to_positive_definite(A)
Return a value tau such that A + tau*I
"""
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)
2.0
eigvals(A + 2*I)
3-element Vector{Float64}:
0.2679491924311227
2.9999999999999996
3.7320508075688776
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]
3×3 Matrix{Float64}:
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
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.012848604931466667
0.012848604931466667
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) = 4.12429806439097
4.12429806439097
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
Diagonal(F.L*F.D*(F.L)' - A)
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