JuliaStats / PDMats.jl

Uniform Interface for positive definite matrices of various structures
Other
104 stars 43 forks source link

Account for pivoting in sparse Cholesky decompositions #175

Closed devmotion closed 11 months ago

devmotion commented 1 year ago

As explained in the docstring, one should be careful when working with Cholesky decompositions of sparse arrays since the property :L does not take into account the pivoting - in contrast to :PtL which also includes the necessary permutations.

This PR changes chol_lower and chol_upper for such sparse Cholesky decompositions such that pivoting is taken into account. Unfortunately, since SparseArrays does not support multiplication of such FactorComponents, code that multiplies the output of chol_lower/chol_upper requires a different, explicitly constructed sparse equivalent of these factors. Alternatively, one could define internal functions for multiplying the output of chol_lower/chol_upper, let them fall back to * by default, and add special implementations for these sparse factors. I've not done this for now because I was not completely sure whether the sparse special case justifies such changes to the (internal) API.

Additionally, the PR adds a seemingly missing definition of chol_upper(::Matrix) and fixes a few test errors in test/pdmtypes.jl and test/specialarrays.jl caused by upstream bugs.

Fixes #120.

Edit: I suggest merging #180 first. That PR fixes the test issues on the master branch properly whereas this PR here only contains a workaround.

devmotion commented 1 year ago

Test failures on Julia 1.6 seem to be caused by an upstream issue in ArrayLayouts: https://github.com/JuliaLinearAlgebra/ArrayLayouts.jl/pull/161

dkarrasch commented 1 year ago

Would it make sense to make an extra package for the PDSparseMat stuff? According to juliahub.com, PDMats.jl has 1k+ dependents, which, starting from v1.10, will all need to load SparseArrays.jl (no longer in the sysimage) by default. That is a heavy load for the ecosystem, which would be good to avoid, if possible.

devmotion commented 1 year ago

Maybe, I'm not sure. It would require a breaking change though so it's also unclear whether/how quickly the package ecosystem would adopt the breaking release.

In any case, I think this consideration is unrelated to the PR. The PR just fixes bugs in the current release of PDMats.

codecov-commenter commented 11 months ago

Codecov Report

Patch coverage is 100.00% of modified lines.

Files Changed Coverage
src/chol.jl 100.00%
src/pdsparsemat.jl 100.00%

:loudspeaker: Thoughts on this report? Let us know!.

devmotion commented 11 months ago

This PR is ready for review as well @andreasnoack 🙂