Integrate CliqueTrees.jl - #586
Conversation
|
Your PR requires formatting changes to meet the project's style guidelines. Click here to view the suggested changes.diff --git a/lib/MadNLPCliqueTrees/src/MadNLPCliqueTrees.jl b/lib/MadNLPCliqueTrees/src/MadNLPCliqueTrees.jl
index 045d4e9a..1360db60 100644
--- a/lib/MadNLPCliqueTrees/src/MadNLPCliqueTrees.jl
+++ b/lib/MadNLPCliqueTrees/src/MadNLPCliqueTrees.jl
@@ -55,35 +55,35 @@ mutable struct CliqueTreesSolver{T, F <: Factorization{T}} <: AbstractLinearSolv
logger::MadNLPLogger
end
-function _build_factorization(tril::SparseMatrixCSC{T, Int32}, opt::CliqueTreesOptions, ::Val{LDL}) where T
+function _build_factorization(tril::SparseMatrixCSC{T, Int32}, opt::CliqueTreesOptions, ::Val{LDL}) where {T}
S = Symmetric(tril, :L)
F = ChordalLDLt{:L}(S; alg = opt.cliquetrees_ordering)::FChordalLDLt{:L, T, Int32}
return F
end
-function _build_factorization(tril::SparseMatrixCSC{T, Int32}, opt::CliqueTreesOptions, ::Val{CHOLESKY}) where T
+function _build_factorization(tril::SparseMatrixCSC{T, Int32}, opt::CliqueTreesOptions, ::Val{CHOLESKY}) where {T}
S = Symmetric(tril, :L)
F = ChordalCholesky{:L}(S; alg = opt.cliquetrees_ordering)::FChordalCholesky{:L, T, Int32}
return F
end
function CliqueTreesSolver(
- tril::SparseMatrixCSC{T, Int32};
- opt = CliqueTreesOptions(),
- logger = MadNLPLogger(),
- pos = size(tril, 1),
-) where T
+ tril::SparseMatrixCSC{T, Int32};
+ opt = CliqueTreesOptions(),
+ logger = MadNLPLogger(),
+ pos = size(tril, 1),
+ ) where {T}
signs = fill(-1, size(tril, 1)); signs[1:pos] .= 1
F = _build_factorization(tril, opt, Val(opt.cliquetrees_algorithm))
return CliqueTreesSolver{T, typeof(F)}(tril, F, signs, opt, logger)
end
-function factorize!(M::CliqueTreesSolver{T, <:ChordalLDLt}) where T
- ldlt!(copy!(M.F, M.tril), M.opt.cliquetrees_strategy; signs=M.signs, reg=M.opt.cliquetrees_regularization, check=false)
+function factorize!(M::CliqueTreesSolver{T, <:ChordalLDLt}) where {T}
+ ldlt!(copy!(M.F, M.tril), M.opt.cliquetrees_strategy; signs = M.signs, reg = M.opt.cliquetrees_regularization, check = false)
return M
end
-function solve_linear_system!(M::CliqueTreesSolver{T, <:ChordalLDLt}, rhs::Vector{T}) where T
+function solve_linear_system!(M::CliqueTreesSolver{T, <:ChordalLDLt}, rhs::Vector{T}) where {T}
if issuccess(M.F)
ldiv!(M.F, rhs)
end
@@ -91,7 +91,7 @@ function solve_linear_system!(M::CliqueTreesSolver{T, <:ChordalLDLt}, rhs::Vecto
return rhs
end
-function inertia(M::CliqueTreesSolver{T, <:ChordalLDLt}) where T
+function inertia(M::CliqueTreesSolver{T, <:ChordalLDLt}) where {T}
d = M.F.d
pos = 0; zer = 0; neg = 0
@inbounds for di in d
@@ -106,15 +106,15 @@ function inertia(M::CliqueTreesSolver{T, <:ChordalLDLt}) where T
return pos, zer, neg
end
-introduce(::CliqueTreesSolver{T, <:ChordalLDLt}) where T =
+introduce(::CliqueTreesSolver{T, <:ChordalLDLt}) where {T} =
"CliqueTrees/LDLᵀ v$(pkgversion(CliqueTrees))"
-function factorize!(M::CliqueTreesSolver{T, <:ChordalCholesky}) where T
- cholesky!(copy!(M.F, M.tril), M.opt.cliquetrees_strategy; check=false)
+function factorize!(M::CliqueTreesSolver{T, <:ChordalCholesky}) where {T}
+ cholesky!(copy!(M.F, M.tril), M.opt.cliquetrees_strategy; check = false)
return M
end
-function solve_linear_system!(M::CliqueTreesSolver{T, <:ChordalCholesky}, rhs::Vector{T}) where T
+function solve_linear_system!(M::CliqueTreesSolver{T, <:ChordalCholesky}, rhs::Vector{T}) where {T}
if issuccess(M.F)
ldiv!(M.F, rhs)
end
@@ -122,7 +122,7 @@ function solve_linear_system!(M::CliqueTreesSolver{T, <:ChordalCholesky}, rhs::V
return rhs
end
-function inertia(M::CliqueTreesSolver{T, <:ChordalCholesky}) where T
+function inertia(M::CliqueTreesSolver{T, <:ChordalCholesky}) where {T}
n = size(M.tril, 1)
if issuccess(M.F)
@@ -132,14 +132,14 @@ function inertia(M::CliqueTreesSolver{T, <:ChordalCholesky}) where T
end
end
-introduce(::CliqueTreesSolver{T, <:ChordalCholesky}) where T =
+introduce(::CliqueTreesSolver{T, <:ChordalCholesky}) where {T} =
"CliqueTrees/Cholesky v$(pkgversion(CliqueTrees))"
is_inertia(::CliqueTreesSolver) = true
improve!(::CliqueTreesSolver) = false
input_type(::Type{<:CliqueTreesSolver}) = :csc
default_options(::Type{<:CliqueTreesSolver}) = CliqueTreesOptions()
-is_supported(::Type{<:CliqueTreesSolver}, ::Type{T}) where T <: AbstractFloat = true
+is_supported(::Type{<:CliqueTreesSolver}, ::Type{T}) where {T <: AbstractFloat} = true
export CliqueTreesSolver, CliqueTreesOptions
diff --git a/lib/MadNLPCliqueTrees/test/runtests.jl b/lib/MadNLPCliqueTrees/test/runtests.jl
index a8ae9fd4..f9545cad 100644
--- a/lib/MadNLPCliqueTrees/test/runtests.jl
+++ b/lib/MadNLPCliqueTrees/test/runtests.jl
@@ -10,12 +10,12 @@ using SparseArrays
# Cholesky algorithm on positive definite matrix
@testset "Cholesky algorithm" begin
for T in (Float32, Float64)
- row = Int32[1,2,2]; col = Int32[1,1,2]; val = T[1., .1, 2.]
+ row = Int32[1, 2, 2]; col = Int32[1, 1, 2]; val = T[1.0, 0.1, 2.0]
b = T[1.0, 3.0]
sol = T[0.8542713567839195, 1.4572864321608041]
csc = sparse(row, col, val, 2, 2)
- opt = CliqueTreesOptions(cliquetrees_algorithm=MadNLP.CHOLESKY)
- M = CliqueTreesSolver(csc; opt=opt)
+ opt = CliqueTreesOptions(cliquetrees_algorithm = MadNLP.CHOLESKY)
+ M = CliqueTreesSolver(csc; opt = opt)
MadNLP.factorize!(M)
@test MadNLP.is_inertia(M)
@test MadNLP.inertia(M) == (2, 0, 0)
diff --git a/src/LinearSolvers/cholmod.jl b/src/LinearSolvers/cholmod.jl
index 868bf762..4de47861 100644
--- a/src/LinearSolvers/cholmod.jl
+++ b/src/LinearSolvers/cholmod.jl
@@ -17,7 +17,7 @@ end
function CHOLMODSolver(
csc::SparseMatrixCSC{T};
- opt=CHOLMODOptions(), logger=MadNLPLogger(), pos=nothing,
+ opt = CHOLMODOptions(), logger = MadNLPLogger(), pos = nothing,
) where T
p = Vector{Float64}(undef,csc.n)
d = Vector{Float64}(undef,csc.n)
diff --git a/src/LinearSolvers/lapack.jl b/src/LinearSolvers/lapack.jl
index 88a2e0b1..f2fb29d1 100644
--- a/src/LinearSolvers/lapack.jl
+++ b/src/LinearSolvers/lapack.jl
@@ -22,7 +22,7 @@ mutable struct LapackCPUSolver{T, MT} <: AbstractLinearSolver{T}
A::MT;
opt=LapackOptions(),
logger=MadNLPLogger(),
- pos=nothing,
+ pos = nothing,
) where {MT <: AbstractMatrix}
T = eltype(A)
m,n = size(A)
diff --git a/src/LinearSolvers/ldl.jl b/src/LinearSolvers/ldl.jl
index bf287b8b..87097fbe 100644
--- a/src/LinearSolvers/ldl.jl
+++ b/src/LinearSolvers/ldl.jl
@@ -13,7 +13,7 @@ end
function LDLSolver(
tril::SparseMatrixCSC{T};
- opt=LDLFactorizationsOptions(), logger=MadNLPLogger(), pos=nothing,
+ opt = LDLFactorizationsOptions(), logger = MadNLPLogger(), pos = nothing,
) where T
# TODO: convert tril to triu, not full
full, tril_to_full_view = get_tril_to_full(T,tril)
diff --git a/src/LinearSolvers/mumps.jl b/src/LinearSolvers/mumps.jl
index 8d78a931..86bd01e1 100644
--- a/src/LinearSolvers/mumps.jl
+++ b/src/LinearSolvers/mumps.jl
@@ -164,7 +164,7 @@ end
# ---------------------------------------------------------------------------------------
function MumpsSolver(csc::SparseMatrixCSC{T,Int32};
- opt=MumpsOptions(), logger=MadNLPLogger(), pos=nothing,
+ opt = MumpsOptions(), logger = MadNLPLogger(), pos = nothing,
) where T
I,J = findIJ(csc)
diff --git a/src/LinearSolvers/umfpack.jl b/src/LinearSolvers/umfpack.jl
index 244eb24b..1a2a7f3c 100644
--- a/src/LinearSolvers/umfpack.jl
+++ b/src/LinearSolvers/umfpack.jl
@@ -21,7 +21,7 @@ end
function UmfpackSolver(
csc::SparseMatrixCSC{T};
- opt=UmfpackOptions(), logger=MadNLPLogger(), pos=nothing,
+ opt = UmfpackOptions(), logger = MadNLPLogger(), pos = nothing,
) where T
p = Vector{Float64}(undef,csc.n)
d = Vector{Float64}(undef,csc.n) |
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #586 +/- ##
==========================================
+ Coverage 84.16% 89.85% +5.68%
==========================================
Files 52 37 -15
Lines 4680 3677 -1003
==========================================
- Hits 3939 3304 -635
+ Misses 741 373 -368 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
|
Thanks for working on that @samuelsonric!
julia> A = sparse(Float128[1 1; 1 0])
2×2 SparseMatrixCSC{Float128, Int64} with 3 stored entries:
1.00000000000000000000000000000000000e+00 1.00000000000000000000000000000000000e+00
1.00000000000000000000000000000000000e+00 ⋅
julia> F = ldlt!(ChordalLDLt(A))
2×2 FChordalLDLt{:L, Float128, Int64} with 3 stored entries:
1.00000000000000000000000000000000000e+00 ⋅
1.00000000000000000000000000000000000e+00 1.00000000000000000000000000000000000e+00
1.00000000000000000000000000000000000e+00 ⋅
⋅ -1.00000000000000000000000000000000000e+00
julia> A = sparse(Float64[0 1; 1 1])
2×2 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
⋅ 1.0
1.0 1.0
julia> F = ldlt!(ChordalLDLt(A))
ERROR: SingularException(1)
Stacktrace:
[1] ldlt!(F::FChordalLDLt{:L, Float64, Int64}; check::Bool, reg::Nothing)
@ CliqueTrees.Multifrontal ~/.julia/packages/CliqueTrees/wgHkM/src/Multifrontal.jl/src/ldlt.jl:66
[2] ldlt!(F::FChordalLDLt{:L, Float64, Int64})
@ CliqueTrees.Multifrontal ~/.julia/packages/CliqueTrees/wgHkM/src/Multifrontal.jl/src/ldlt.jl:46
[3] top-level scope
@ REPL[15]:1
julia> A = sparse(Float64[0 1; 1 0])
2×2 SparseMatrixCSC{Float64, Int64} with 2 stored entries:
⋅ 1.0
1.0 ⋅
julia> F = ldlt!(ChordalLDLt(A))
ERROR: SingularException(1)
Stacktrace:
[1] ldlt!(F::FChordalLDLt{:L, Float128, Int64}; check::Bool, reg::Nothing)
@ CliqueTrees.Multifrontal ~/.julia/packages/CliqueTrees/wgHkM/src/Multifrontal.jl/src/ldlt.jl:66
[2] ldlt!(F::FChordalLDLt{:L, Float128, Int64})
@ CliqueTrees.Multifrontal ~/.julia/packages/CliqueTrees/wgHkM/src/Multifrontal.jl/src/ldlt.jl:46
[3] top-level scope
@ REPL[11]:1It may be slower but we need its robustness for our KKT systems. |
|
I think @samuelsonric Something that will be very useful for us is to regularize on the fly the pivots in the sparse Cholesky if the inertia is wrong (matrix not SPD). |
|
This is great! Looks like the solver performs quite robustly with |
|
Hello @amontoison. Element TypesThe algorithm can handle matrices with any element type, since we have defined fallbacks for each BLAS/LAPACK kernel we use.
These fallbacks are sometimes a little primitive, but I can give them attention if needed. 1x1 PivotingPivoting in sparse matrix factorization is tricky, since numerical pivoting can destroy sparsity. The simplest approach is to perform 1x1 pivoting within dense blocks (supernodes). This is done if you pass julia> A = sparse(Float64[0 1; 1 1]);
2×2 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
⋅ 1.0
1.0 1.0
julia> ldlt!(ChordalLDLt(A), RowMaximum())
2×2 FChordalLDLt{:L, Float64, Int64} with 3 stored entries:
1.0 ⋅
1.0 1.0
1.0 ⋅
⋅ -1.0Beware! This changes Dynamic RegularizationA very simple type of dynamic regularization is already implemented. It is the same strategy used in QDLDL.jl and LDLFactorizations.jl: small pivots are replaced with julia> ldlt!(ChordalLDLt(A); signs=[1, 1], reg=DynamicRegularization())
2×2 FChordalLDLt{:L, Float64, Int64} with 3 stored entries:
1.0 ⋅
1.0e6 1.0
1.0e-6 ⋅
⋅ 1.0e-6I am going to try implementing one of the strategies in the blog post as well. |
|
@amontoison Try this: Gill-Murray-Wright 1981. julia> using CliqueTrees.Multifrontal, LinearAlgebra
julia> A = [
1 1 2
1 1 + 1e-20 3
2 3 1
];
julia> F = ChordalLDLt(A);
julia> ldlt!(F; signs=[1, 1, 1], reg=GMW81(F));
julia> F.P' * F.L * F.D * F.L' * F.P
3×3 Matrix{Float64}:
3.77124 1.0 2.0
1.0 6.01561 3.0
2.0 3.0 3.24264 |
|
You can use it with and without pivoting; this is triggered by passing |
|
Wow! Thanks a lot @samuelsonric !!! |
|
I only implemented LDL, but Cholesky is very similar, and I could implement at as well. |
|
Ok, it makes sense. 👍 |
|
@samuelsonric |
|
@samuelsonric
|
I think that these algorithms are used in LU and symmetric indefinite factorizations, not Cholesky / LDLt. If I implement a symmetric indefinite factorization, then I suppose that I will also have to implement matching. By symmetric indefinite factorization, I mean LBLt with 2x2 pivoting. |
|
@amontoison I added an option Additionally, if the information is known, it would be good to pass a keyword argument |
|
Is it not possible to define a function |
|
The failing tests seem to be unrelated to my contribution. |
|
@amontoison My initial implementation of GMW81 was very buggy; that should no longer be the case. I implemented a second modified Cholesky algorithm as well -- SE99 -- with better theoretical properties. These can be used both with julia> using CliqueTrees.Multifrontal, LinearAlgebra
julia> A = [
1 1 2
1 1 3
2 3 1
];
julia> isposdef(A)
false
julia> F = ChordalCholesky(A);
julia> cholesky!(F; reg=SE99(F));
julia> B = Matrix(F)
3×3 Matrix{Float64}:
3.0 1.0 2.0
1.0 3.0 3.0
2.0 3.0 3.375
julia> isposdef(B)
trueWhere were you thinking about using modified Cholesky? |
|
Thank you @samuelsonric ! |
|
@amontoison I updated the workflow. Can you explain what you mean by this?
CliqueTrees needs to be imported somwhere. |
|
@samuelsonric does CliqueTrees have an API for uniform batch factorization/solve? It's not required for this PR/MadNLP, but it would be super cool to try madsuite-org/MadIPM.jl#91 on CPU with a proper batched solver. |
|
@klamike No but I can try writing it. Did you mean to say GPU? |
|
I meant CPU, on NVIDIA GPUs we already have CUDSS. But a backend-agnostic CliqueTrees on GPU would be very cool too 😄 |
|
@klamike Cholesky or LDLt or both? |
|
LDLt :) |
|
It's also relevant for #598 |
This PR adds integrates CliqueTrees.jl with MadNLP.jl.
CliqueTrees.jl contains a supernodal LDLt factorization routine. Unlike LDLFactorizations.jl and QDLDL.jl, which process the matrix column-by-coumn, the CliqueTrees.jl uses BLAS level 3 kernels to handle dense chunks. This can significantly improve factorization speed. Here is an example.
Here are the results.
CliqueTreesSolverLDLSolver