Unify PDMat and PDSparseMat + move SparseArrays support to an extension
#188
Conversation
Codecov ReportAttention:
Additional details and impacted files@@ Coverage Diff @@
## master #188 +/- ##
==========================================
+ Coverage 91.61% 92.18% +0.56%
==========================================
Files 9 11 +2
Lines 680 640 -40
==========================================
- Hits 623 590 -33
+ Misses 57 50 -7 ☔ View full report in Codecov by Sentry. |
There was a problem hiding this comment.
Thank you so much for taking up the challenge and moving this forward! I think this is a great service to the ecosystem. Perhaps we can make this work even for older Julia versions by using SuiteSparse or SparseArrays.SuiteSparse in the extension? On v1.9+, this should be a no-op, but on older versions (without Pkg extensions) this will just get loaded by including the extension code files. I'm not sure though if that would require "version-dependent" Project.toml files or other hassle.
AFAICT on older Julia versions we would have to keep both SparseArrays and SuiteSparse as dependencies. However, if only Julia >= 1.9 the extension only depends on SparseArrays, then SuiteSparse would still be a regular dependency and since it depends on SparseArrays the extension would be loaded unconditionally. Alternatively, the extension could depend on both SparseArrays and SuiteSparse (or only SuiteSparse) - but then users would have to load both SparseArrays and SuiteSparse (or only SuiteSparse). Both alternatives seem suboptimal. However, I just thought of another option that I will try later: Making SparseArrays and SuiteSparse both a weak dependency but letting the extension only depend on SparseArrays. |
|
Seems to work 🙂 On Julia 1.0 and 1.6, PDMats depends on SparseArrays and SuiteSparse (https://github.com/JuliaStats/PDMats.jl/actions/runs/6404900221/job/17386347115?pr=188#step:5:36, https://github.com/JuliaStats/PDMats.jl/actions/runs/6404900221/job/17386347711?pr=188#step:5:27) and on Julia 1.9 and nightly neither SparseArrays nor SuiteSparse are dependencies (https://github.com/JuliaStats/PDMats.jl/actions/runs/6404900221/job/17386348265?pr=188#step:5:24, https://github.com/JuliaStats/PDMats.jl/actions/runs/6404900221/job/17386348761?pr=188#step:5:24) and the tests pass successfully with installing and loading only SparseArrays. |
|
I'll update the PR once the other PRs (non-breaking bugfixes and possibly a new feature) are merged. |
There was a problem hiding this comment.
Nice work!
One other thing we might want to look at, quite a few of the methods starting with ### (un)whitening rely on the Choleksy factorization. If there isn't time to review/generalize those algorithms, perhaps it might make sense to define a PDMatCholesky{T,S<:AbstractMatrix{T}} = PDMat{T,S,<:Cholesky} and then restrict the argument of those methods to be a PDMatCholesky? That way we're at least not promising things we can't deliver.
| end | ||
|
|
||
| function PDMat(mat::AbstractMatrix,chol::Cholesky{T,S}) where {T,S} | ||
| function PDMat(mat::AbstractMatrix{T}, chol::Cholesky) where {T<:Real} |
There was a problem hiding this comment.
| function PDMat(mat::AbstractMatrix{T}, chol::Cholesky) where {T<:Real} | |
| function PDMat(mat::AbstractMatrix{T}, fmat::Factorization) where {T<:Real} |
There was a problem hiding this comment.
Maybe we should wait with generalizing the constructor until we actually support and test other factorizations.
Co-authored-by: Tim Holy <[email protected]>
|
Bump |
|
The PR is the result of some discussions I had with @andreasnoack about |
|
Right, I wasn't bumping you 🙂 |
| mat::S # input matrix | ||
| chol::Cholesky{T,S} # Cholesky factorization of mat | ||
| fact::F # factorization of mat |
There was a problem hiding this comment.
One needs an inner constructor that checks isposdef(fact) before calling new. Otherwise I could call
E = eigen(A)
PDMat{eltype(A),typeof(A),typeof(E)}(A, E)even if A is not posdef.
There was a problem hiding this comment.
Currently none of the other constructors enforces positive semi-definite matrices: #22
So ideally this should be addressed more generally, to ensure that such checks are performed by all AbstractPDMat types. But I wonder if it should be left for a separate PR?
There was a problem hiding this comment.
Sure. But it's not a hard change, we should just do it once other dust has settled. Given that the PDMat wrapper promises positive-definiteness, it is kinda important to enforce this, especially if we widen the type definition so that it can no longer be guaranteed from the input types alone.
| @test M isa PDMat{Int,<:SymTridiagonal,<:Cholesky} | ||
| @test M == S | ||
| end | ||
| end |
There was a problem hiding this comment.
Probably want to add tests for what happens if one tries to construct with factorizations besides cholesky.
This PR unifies
PDMatandPDSparseMatand makesPDSparseMatan (unexported) alias ofPDMat. This makes it also possible to move support forSparseArrays to an extension, without any workarounds such as non-functional types in the main package (see #174).Initially,
PDMatandPDSparseMatonly supportedMatrix{Float64}andSparseMatrixCSC{Float64}matrices. This restriction was lifted in #37 but the separate types with support forCholeskyandCHOLMOD.Factors were kept.Due to the refactoring of SuiteSparse and SparseArrays (it was integrated in SparseArrays in Julia 1.9), on Julia >= 1.9 only the SparseArrays dependency is needed and SuiteSparse only loads functionality from SparseArrays (https://github.com/JuliaSparse/SuiteSparse.jl/blob/e8285dd13a6d5b5cf52d8124793fc4d622d07554/src/SuiteSparse.jl). It seemed inconvenient to demand that users on Julia >= 1.9 load both SparseArrays and SuiteSparse for the sparse extension, and hence I decided to let the extension only depend on SparseArrays. IIUC, however, this means that we can't support Julia < 1.9 anymore since there SparseArrays does not contain the requiredCHOLMODfunctionality. @dkarrasch, is this correct?Since this means we drop support for the current LTS version, my suggestion would be to backport bugfixes to older Julia versions (at least >= 1.6) until another Julia version >= 1.9 becomes the new LTS version.Edit: The PR keeps support for Julia 1.0-, I found a workaround: #188 (comment)
Fixes #174.