Code Monkey home page Code Monkey logo

limitedldlfactorizations.jl's Introduction

Limited-Memory LDLᵀ Factorization

DOI CI Build Status codecov Documentation/stable Documentation/dev

A Port of LLDL to Julia. See https://github.com/optimizers/lldl.

Please cite this repository if you use LimitedLDLFactorizations.jl in your work: see CITATION.cff.

LimitedLDLFactorizations is a limited-memory LDLᵀ factorization for symmetric matrices. Given a symmetric matrix A, we search for a unit lower triangular L, a diagonal D and a diagonal ∆ such that LDLᵀ is an incomplete factorization of A+∆. The elements of the diagonal matrix ∆ have the form ±α, where α ≥ 0 is a shift.

It is possible to only supply the lower triangle of A and/or a prescribed permutation that attempts to diminish fill-in. AMD.jl and Metis.jl are recommended packages for computing fill-reducing orderings of sparse matrices.

Installing

julia> ]
pkg> add LimitedLDLFactorizations

Brief Description

The only functions exported are lldl, \, ldiv! and nnz. Supply a dense array or sparse matrix to lldl. Dense arrays will be converted to sparse. The strict lower triangle and diagonal of sparse matrices will be extracted.

Optionally, supply

  • a memory parameter to allow more fill in the L factor;
  • a drop tolerance to discard small elements in the L factor;
  • an initial shift to speed up the process in case the factorization does not exist without shift.

Using a memory parameter larger than or equal to the size of A will yield an exact factorization provided one exists with the permutation supplied. In particular, the full factorization exists for any symmetric permutation of a symmetric quasi-definite matrix.

lldl returns a factorization in the form of a LimitedLDLFactorization object. The \ and ldiv! methods are implemented for objects of type LimitedLDLFactorization

More Examples

See examples/example.jl and tests/runtest.jl.

Complete Description

[1] C.-J. Lin and J. J. Moré. Incomplete Cholesky factorizations with limited memory. SIAM Journal on Scientific Computing, 21(1):24--45, 1999. DOI 10.1137/S1064827597327334.
[2] http://www.mcs.anl.gov/~more/icfs
[3] D. Orban. Limited-Memory LDLᵀ Factorization of Symmetric Quasi-Definite Matrices with Application to Constrained Optimization. Numerical Algorithms 70(1):9--41, 2015. DOI 10.1007/s11075-014-9933-x.
[4] https://github.com/optimizers/lldl

Bug reports and discussions

If you think you found a bug, feel free to open an issue. Focused suggestions and requests can also be opened as issues. Before opening a pull request, start an issue or a discussion on the topic, please.

If you want to ask a question not suited for a bug report, feel free to start a discussion here. This forum is for general discussion about this repository and the JuliaSmoothOptimizers, so questions about any of our packages are welcome.

limitedldlfactorizations.jl's People

Contributors

abelsiqueira avatar amontoison avatar dpo avatar geoffroyleconte avatar github-actions[bot] avatar jsobot avatar juliatagbot avatar monssaftoukal avatar t-bltg avatar tgalizzi avatar tmigot avatar

Stargazers

avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar

Watchers

avatar avatar avatar avatar

Forkers

abelsiqueira tgalizzi stjordanis amontoison geoffroyleconte t-bltg proofofconceptforjulismoothoptimizers web5530top jefffessler

limitedldlfactorizations.jl's Issues

Add Support for Matrices with `Int32` and `UInt32` as Indices Vectors

It seems the code fails when the input sparse matrix has indices with types which are not Int64:

using SparseArrays;
using LimitedLDLFactorizations;

mA = sprand(10, 10, 0.1);
vI, vJ, vK = findnz(mA);
mB = sparse(Int32.(vI), Int32.(vJ), vK, mA.m, mA.n);

lldl(mA); #<! Works
lldl(mB); #<! Fails

I get the following error:

julia> lldl(mB)
ERROR: MethodError: no method matching abspermute!(::Vector{Float64}, ::SubArray{Int32, 1, Vector{Int32}, Tuple{UnitRange{Int64}}, true}, ::Int64)

I think the code should support any integer as indices and any float as value. Wasn't that the design goal?

Performance wise, if the input data is 32 Bit having all other computations in 32 but might assist with performance (Memory and calculation throughput).

For the UInt32 case:

mB = sparse(UInt32.(vI), UInt32.(vJ), vK, mA.m, mA.n);
lldl(mB); #<! Fails

The error:

julia> lldl(mB); #<! Fails
ERROR: MethodError: no method matching amd(::SparseMatrixCSC{Float64, UInt32})

It also fails for UInt64.

Remark: I didn't make the matrix symmetric as it is has no significance for the error.

Issue with the operator \ in Julia 1.2.0

Based on the tests provided by the package, I tried the following script using Julia version 1.2.0:

using LinearAlgebra, SparseArrays, Test

using AMD, Metis, LimitedLDLFactorizations

  A = [
    1.7 0 0 0 0 0 0 0 0.13 0
    0 1.0 0 0 0.02 0 0 0 0 0.01
    0 0 1.5 0 0 0 0 0 0 0
    0 0 0 1.1 0 0 0 0 0 0
    0 0.02 0 0 2.6 0 0.16 0.09 0.52 0.53
    0 0 0 0 0 1.2 0 0 0 0
    0 0 0 0 0.16 0 1.3 0 0 0.56
    0 0 0 0 0.09 0 0 1.6 0.11 0
    0.13 0 0 0 0.52 0 0 0.11 1.4 0
    0 0.01 0 0 0.53 0 0.56 0 0 3.1
  ]
  A = sparse(A)

    LLDL = lldl(A, memory = 10)

    sol = ones(A.n)
    b = A * sol
    x = LLDL \ b

Error message:

$ julia test.jl 
ERROR: LoadError: MethodError: no method matching adjoint(::Tuple{SparseMatrixCSC{Float64,Int64},Array{Float64,1},Float64,Array{Int64,1}})
Closest candidates are:
  adjoint(!Matched::Missing) at missing.jl:79
  adjoint(!Matched::LightGraphs.DefaultDistance) at /xxx/.julia/packages/LightGraphs/IgJif/src/distance.jl:22
  adjoint(!Matched::Number) at number.jl:193
  ...
Stacktrace:
 [1] \(::Tuple{SparseMatrixCSC{Float64,Int64},Array{Float64,1},Float64,Array{Int64,1}}, ::Array{Float64,1}) at ./operators.jl:563
 [2] top-level scope at /xxxx/test.jl:23
 [3] include at ./boot.jl:328 [inlined]
 [4] include_relative(::Module, ::String) at ./loading.jl:1094
 [5] include(::Module, ::String) at ./Base.jl:31
 [6] exec_options(::Base.JLOptions) at ./client.jl:295
 [7] _start() at ./client.jl:464
in expression starting at /xxx/test.jl:23

I cannot understand the issue with the script. Help would be welcome. Thank you.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Package name

@dpo, what do you wanna do about the package names of LLDL and LDL?

Update README

The README needs to be updated to say that implicit ordering is now supported and should give an example.

Support for Modified Incomplete Choleksy Decomposition

In the book Yousef Saad - Iterative Methods for Sparse Linear Systems (2nd Edition):

image

This method, originally for the ILU decomposition, can improve the performance of the conditioner greatly.
Basically it enforces the conditioner to comply A e = L L' e where e is a vector of ones.
Namely the sum of values of a row matches the input matrix.
The way it is implemented is that any component which is dropped in the incomplete decomposition is compensated by the diagonal component.

Since component of the row is computed anyway it should come almost freely (Unless I'm wrong and the values are not computed).

Error with the permutation

K2_5bug.csv

using DelimitedFiles, SparseArrays, LinearAlgebra
path_bug = "K2_5bug.csv"
K2_5 = sparse(readdlm(path_bug, ',', Float64))
using LimitedLDLFactorizations
K2_5_lowtri = tril(K2_5)
F = lldl(K2_5_lowtri, memory=10)

This closes julia. The matrix that I am trying to factorize is SQD and nearly singular. However this works fine:

F = lldl(K2_5_lowtri, 1: K2_5_lowtri.n, memory=10)

Removing the @inbounds and @simd macros should help to locate the issue but I could not manage to solve it.

Testing fails on 1.7.0

Running ] test LimitedLDLFactorizations on 1.7.0 fails:

[...]
     Testing Running tests...
SQD no shift: Error During Test at ~/.julia/packages/LimitedLDLFactorizations/r21rs/test/runtests.jl:5
  Got exception outside of a @test
  ArgumentError: Illegal buffers for SparseMatrixCSC construction 10 [1, 2, 4, 4, 4, 8, 8, 11, 13, 14, 14] [9, 5, 10, 7, 8, 9, 10, 8, 9, 10, 9, 10, 10, 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, 9, 5, 10, 7, 8, 9, 10, 10, 9] [0.07647058823529412, 0.02, 0.01, 0.06154793045083859, 0.0346207108785967, 0.20003077396522545, 0.2038005847053393, -0.00429353465240259, -0.02480708910277052, 0.4087826636611924, 0.05752526570865536, -0.010068305077340346, -0.07185227820756272, 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, 0.0, 0.0, 0.0, 0.0, 0.0, 0.08120589331564292, 0.012147079105600677, 0.005594312599767557, 0.08143008846051507, 0.0431345532201147, 0.25763282748712707, 0.1801245921491028, 0.26251723244756425, 0.07079149130272516]
  Stacktrace:
    [1] SparseMatrixCSC
      @ $SRC_DIR/usr/share/julia/stdlib/v1.7/SparseArrays/src/sparsematrix.jl:29 [inlined]
    [2] lldl(T::SparseMatrixCSC{Float64, Int64}, adiag::SparseVector{Float64, Int64}, P::UnitRange{Int64}; memory::Int64, α::Float64, droptol::Float64)
      @ LimitedLDLFactorizations ~/.julia/packages/LimitedLDLFactorizations/r21rs/src/LimitedLDLFactorizations.jl:229
    [3] #lldl#3
      @ ~/.julia/packages/LimitedLDLFactorizations/r21rs/src/LimitedLDLFactorizations.jl:65 [inlined]
    [4] macro expansion
      @ ~/.julia/packages/LimitedLDLFactorizations/r21rs/test/runtests.jl:27 [inlined]
    [5] macro expansion
      @ $SRC_DIR/usr/share/julia/stdlib/v1.7/Test/src/Test.jl:1283 [inlined]
    [6] top-level scope
      @ ~/.julia/packages/LimitedLDLFactorizations/r21rs/test/runtests.jl:7
    [7] include(fname::String)
      @ Base.MainInclude ./client.jl:451
    [8] top-level scope
      @ none:6
    [9] eval
      @ ./boot.jl:373 [inlined]
   [10] exec_options(opts::Base.JLOptions)
      @ Base ./client.jl:268
   [11] _start()
      @ Base ./client.jl:495
Test Summary: | Pass  Error  Total
SQD no shift  |    2      1      3
ERROR: LoadError: Some tests did not pass: 2 passed, 0 failed, 1 errored, 0 broken.
in expression starting at ~/.julia/packages/LimitedLDLFactorizations/r21rs/test/runtests.jl:5
ERROR: Package LimitedLDLFactorizations errored during testing

Support multi-vector b 

LLDL = lldl(A, memory = 10)
x = LLDL \ b

currently only vector b is supported. I am wondering if multi-vector (b as a matrix) could be supported (so that x will be 2d as well). A naive multi-threading seems to have too much overhead.

Support for Unitful sparse arrays? [feature request]

The current code for lldl fails for Unitful sparse (or dense) arrays. MWE:

using SparseArrays: spdiagm
using LimitedLDLFactorizations: lldl
using Unitful: s
A = spdiagm(Float32.(1:3)*s^2) # unitful sparse array
lldl(A) # fails

ERROR: MethodError: no method matching lldl(::SparseArrays.SparseMatrixCSC{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}, Int64}, ::Type{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}})
Closest candidates are:
  lldl(::SparseArrays.SparseMatrixCSC{Tv, Ti}; kwargs...) where {Tv<:Number, Ti<:Integer} at ~/.julia/packages/LimitedLDLFactorizations/usbOW/src/LimitedLDLFactorizations.jl:493
  lldl(::SparseArrays.SparseMatrixCSC{Tv, Ti}, ::Type{Tf}; P, memory, droptol, α, α_increase_factor, check_tril) where {Tv<:Number, Ti<:Integer, Tf<:Real} at ~/.julia/packages/LimitedLDLFactorizations/usbOW/src/LimitedLDLFactorizations.jl:471

Some of the lldl code uses Number which is general enough to include unitful values, such as here:

LimitedLDLFactorizations.jl/src/LimitedLDLFactorizations.jl

Line 493 in 51eec8d

lldl(A::SparseMatrixCSC{Tv, Ti}; kwargs...) where {Tv <: Number, Ti <: Integer} =

but then just a bit deeper into the code the types get restricted to Real, like here:

LimitedLDLFactorizations.jl/src/LimitedLDLFactorizations.jl

Line 480 in 51eec8d

) where {Tv <: Number, Ti <: Integer, Tf <: Real}

I realize that Number is general enough to include complex number types and perhaps you don't want to support those? (Though that could be useful too, right?) My Unitful values are reinterpretable as real numbers, but currently not when stored in a sparse matrix: JuliaSparse/SparseArrays.jl#289

I figure it's a long shot, but I thought I'd ask to see if there is any possibility of supporting more general Number types, especially Unitful values and maybe eventually complex values too? Since you have a pure Julia version, it seems like it would be possible!

Work with lower triangle instead of strict lower triangle and diagonal

Is it possible to change the algorithm so that the user can work with a lower triangular matrix without having to call tril(A, -1)? I find this not very convenient since most factorization algorithms require the lower or upper triangle (with the diagonal). We could still allocate adiag = diag(A) if needed.
@dpo I can try to make these changes if you agree.

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.