Strange (and wrong) factorization

[bvdmitri]

Consider this matrix:

F = [42491.1429254459 1.0544416413649244e6 64.9016820609457 1712.2779951809016; 1.0544416413649244e6 2.616823794441869e7 1610.468694700484 42488.422800411565; 64.9016820609457 1610.468694700484 0.10421453600353446 2.6155294717625517; 1712.2779951809016 42488.422800411565 2.6155294717625517 69.0045838263577]

It is not particularly good, but also not awfully bad, eg

det(F) # 0.020309321582317352 which is close to 0, but OK

Here is what I get with the package against the Julia base

julia> cholesky(Positive, F).L
4×4 LowerTriangular{Float64, Matrix{Float64}}:
  206.134       ⋅          ⋅           ⋅ 
 5115.33      40.9241      ⋅           ⋅ 
    0.314852  -0.0025192  0.071248     ⋅ 
    8.30663   -0.0664632  2.98325e-8  1.0

julia> cholesky(F).L
4×4 LowerTriangular{Float64, Matrix{Float64}}:
  206.134       ⋅          ⋅           ⋅ 
 5115.33      40.9241      ⋅           ⋅ 
    0.314852  -0.0025192  0.071248     ⋅ 
    8.30663   -0.0664632  2.98325e-8  0.000236804

Basically the factorizations are almost identical with the exception of the last element. It also gives the wrong inverse

julia> inv(cholesky(Positive, F)) * F
4×4 Matrix{Float64}:
 1.0          -5.05104e-11  -1.11159e-15   0.0805992
 5.60538e-14   1.0          -1.24247e-16  -0.00162406
 1.14508e-13  -7.95464e-11   1.0           4.18713e-7
 2.27374e-13  -7.27596e-12  -1.33227e-15   5.60763e-8

[bvdmitri]

What is the rationale behind this line? https://github.com/timholy/PositiveFactorizations.jl/blob/master/src/cholesky.jl#L96 It basically makes the factorization completely wrong in marginal cases, where the diagonal entry is slightly below the threshold.

[dpascall]

@bvdmitri I have a problem where doing a modified Cholesky is definitely the way to go. Is this (currently) unresolved issue drastic enough mean that this package shouldn't be used, and I should try coding it up for myself rather than trusting the output here?