Positive seems to fail for larger matrices

[oschulz]

Not sure if it's to be expected: Positive seems to fail for larger matrices.

using PositiveFactorizations, LinearAlgebra
isposdef(Matrix(cholesky(Positive, Hermitian(rand(10,10)))))

seems to (almost?) always returns true, but

isposdef(Matrix(cholesky(Positive, Hermitian(rand(20,20)))))

typically returns false.

We (me and @lmh91) have observed this in actual use cases with empirical 300x300 covariance matrices. Sometimes, Positive succeeds in producing a posdef matrix, but sometimes it fails.

[oschulz]

@timholy if it would be helpful for diagnosis (in cases is an actual a bug, not just a known limitation of the algorithm used by Positive), I have one of those troublesome 300x300 covariance matrices as an HDF5 file.

[timholy]

Why are you converting into a Matrix? The factorization itself is (should be, anyway!) rigorously posdef, you're probably seeing roundoff error.

[oschulz]

Ah, right, of course! Thanks, Tim - we just always had that conversion in the code (the's some intermediate function that expects an Array in the workflow), and we just never ran into the rounding error before, even with a 300x300 matrix - just luck, I guess. We should of course go directly to a PDMat.

[timholy]

You can use L = F.L and do everything that you would have done with the original matrix in terms of L*L'. Demo using the famously difficult Hilbert matrix:

julia> H = [1/(i+j-1) for i = 1:20, j=1:20]
20×20 Array{Float64,2}:
 1.0        0.5        0.333333   …  0.0555556  0.0526316  0.05     
 0.5        0.333333   0.25          0.0526316  0.05       0.047619 
 0.333333   0.25       0.2           0.05       0.047619   0.0454545
 0.25       0.2        0.166667      0.047619   0.0454545  0.0434783
 0.2        0.166667   0.142857      0.0454545  0.0434783  0.0416667
 0.166667   0.142857   0.125      …  0.0434783  0.0416667  0.04     
 0.142857   0.125      0.111111      0.0416667  0.04       0.0384615
 0.125      0.111111   0.1           0.04       0.0384615  0.037037 
 0.111111   0.1        0.0909091     0.0384615  0.037037   0.0357143
 0.1        0.0909091  0.0833333     0.037037   0.0357143  0.0344828
 0.0909091  0.0833333  0.0769231  …  0.0357143  0.0344828  0.0333333
 0.0833333  0.0769231  0.0714286     0.0344828  0.0333333  0.0322581
 0.0769231  0.0714286  0.0666667     0.0333333  0.0322581  0.03125  
 0.0714286  0.0666667  0.0625        0.0322581  0.03125    0.030303 
 0.0666667  0.0625     0.0588235     0.03125    0.030303   0.0294118
 0.0625     0.0588235  0.0555556  …  0.030303   0.0294118  0.0285714
 0.0588235  0.0555556  0.0526316     0.0294118  0.0285714  0.0277778
 0.0555556  0.0526316  0.05          0.0285714  0.0277778  0.027027 
 0.0526316  0.05       0.047619      0.0277778  0.027027   0.0263158
 0.05       0.047619   0.0454545     0.027027   0.0263158  0.025641 

julia> cholesky(H)
ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
 [1] checkpositivedefinite at /home/tim/src/julia-1/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/factorization.jl:18 [inlined]
 [2] #cholesky!#124(::Bool, ::typeof(cholesky!), ::Hermitian{Float64,Array{Float64,2}}, ::Val{false}) at /home/tim/src/julia-1/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/cholesky.jl:226
 [3] #cholesky! at ./none:0 [inlined]
 [4] #cholesky!#125(::Bool, ::typeof(cholesky!), ::Array{Float64,2}, ::Val{false}) at /home/tim/src/julia-1/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/cholesky.jl:258
 [5] #cholesky#129 at ./none:0 [inlined]
 [6] cholesky at /home/tim/src/julia-1/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/cholesky.jl:348 [inlined] (repeats 2 times)
 [7] top-level scope at REPL[22]:1

julia> F = cholesky(Positive, H);

julia> L = F.L;

julia> H - L*L'
20×20 Array{Float64,2}:
 0.0  0.0  0.0  0.0           0.0          …   0.0           0.0        
 0.0  0.0  0.0  0.0           0.0              0.0           0.0        
 0.0  0.0  0.0  0.0           0.0              0.0           0.0        
 0.0  0.0  0.0  0.0           0.0              0.0           0.0        
 0.0  0.0  0.0  0.0          -1.38778e-17     -6.93889e-18  -6.93889e-18
 0.0  0.0  0.0  0.0           1.38778e-17  …   6.93889e-18   0.0        
 0.0  0.0  0.0  0.0           0.0              6.93889e-18   0.0        
 0.0  0.0  0.0  0.0           0.0             -6.93889e-18   0.0        
 0.0  0.0  0.0  0.0           0.0              6.93889e-18   0.0        
 0.0  0.0  0.0  1.38778e-17   0.0              0.0           6.93889e-18
 0.0  0.0  0.0  0.0          -1.38778e-17  …   0.0           0.0        
 0.0  0.0  0.0  0.0           1.38778e-17      6.93889e-18  -6.93889e-18
 0.0  0.0  0.0  0.0          -6.93889e-18      4.45848e-12   7.01655e-12
 0.0  0.0  0.0  0.0          -6.93889e-18      0.0           0.0        
 0.0  0.0  0.0  0.0           0.0              1.54326e-12   2.79415e-12
 0.0  0.0  0.0  0.0          -6.93889e-18  …   1.04083e-17   0.0        
 0.0  0.0  0.0  0.0           0.0             -1.04139e-13  -1.98858e-13
 0.0  0.0  0.0  0.0           6.93889e-18     -7.02091e-13  -1.29227e-12
 0.0  0.0  0.0  0.0          -6.93889e-18     -1.0          -2.84616e-12
 0.0  0.0  0.0  0.0          -6.93889e-18     -2.84616e-12  -5.45897e-12

Of course, near-equality is only possible when the original is actually positive-definite. When it is not, there will be (often large) differences. Here you can see it was forced to cheat on the [end-1,end-1] element due to roundoff error. But that looks like the only one.

[oschulz]

Nice to know, thanks Tim! In my specific case, I only have to wrap it in a PDMat in the end, to use it for an MvTDist. But we do some calculations before, I'll have a look an see if we can use that earlier on. It's just weighting and adding covariance matrices though, so comparatively primitive.