jump-dev
Sparse LDL for QuadtoSOC
QuadtoSOC computes a cholesky factorization, which fails when the Q matrix is positive semidefinite but not positive definite:
https://github.com/jump-dev/MathOptInterface.jl/blob/8c8636f47cbcd0bda8494ac9dcb852df75e2bcc5/src/Bridges/Constraint/bridges/quad_to_soc.jl#L71-L85
Per @dpo we could use LDLFactorizations with a pivot tolerance of zero (or HSL/MUMPS, if we want external binary dependencies).
Can you expand on why? If you are just using the package, does it affect your license?
IANAL, but under a conservative interpretation, LGPL is equivalent to GPL for Julia code (the header and linking exceptions for LGPL are very unclear in the context of Julia). GPL is viral and would force MOI and packages that depend on it to be distributed under GPL.
IANAL either and I'm not so sure how much I could change that license. The package essentially started as a line by line translation of Tim Davis's C code (part of UMFPACK). It was later improved and expanded.
Complicated stuff. I did not realize that Julia made LGPL ambiguous, but it kind of makes sense.
IANAL either and I'm not so sure how much I could change that license. The package essentially started as a line by line translation of Tim Davis's C code (part of UMFPACK). It was later improved and expanded.
That makes sense, you don't have choice of licenses in this case.
I don't think we want to add a binary dependency just to compute this factorization, and the license issue is a problem. So not sure what is actionable here?
https://github.com/oxfordcontrol/QDLDL.jl is another possible solution. It's a pure julia library under a permissive license. We'd need to check if it can handle zero pivots.
This license issue is surprising: the LDLt code - I believe this one? - is part of the Julia standard library in LinearAlgebra.ldlt. And yet a translation of the code into Julia becomes incompatibly licensed with Julia?
There also seems to be a positive semi-definite cholesky(A, Val(True)).
https://github.com/JuliaLang/julia/blob/master/THIRDPARTY.md notes that SuiteSparse (and by extension any Julia translations) are GPL and LGPL. If you compile julia with USE_GPL_LIBS=0 then these libraries should disappear from stdlib.
We haven't gotten any requests for this, but it would be reasonable for someone to ask for a version of JuMP that works (possibly with some features disabled) when Julia is compiled with USE_GPL_LIBS=0.
This is actually less obvious than it seems.
As one problem, we're already calling a SuiteSparse method:
julia> Q
2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4 stored entries:
2.0 -2.0
-2.0 2.0
help?> LinearAlgebra.cholesky(LinearAlgebra.Symmetric(Q))
cholesky(A::SparseMatrixCSC; shift = 0.0, check = true, perm = nothing) -> CHOLMOD.Factor
Compute the Cholesky factorization of a sparse positive definite matrix A. A must be a SparseMatrixCSC or a Symmetric/Hermitian view of a
SparseMatrixCSC. Note that even if A doesn't have the type tag, it must still be symmetric or Hermitian. If perm is not given, a
fill-reducing permutation is used. F = cholesky(A) is most frequently used to solve systems of equations with F\b, but also the methods diag,
det, and logdet are defined for F. You can also extract individual factors from F, using F.L. However, since pivoting is on by default, the
factorization is internally represented as A == P'*L*L'*P with a permutation matrix P; using just L without accounting for P will give
incorrect answers. To include the effects of permutation, it's typically preferable to extract "combined" factors like PtL = F.PtL (the
equivalent of P'*L) and LtP = F.UP (the equivalent of L'*P).
When check = true, an error is thrown if the decomposition fails. When check = false, responsibility for checking the decomposition's
validity (via issuccess) lies with the user.
Setting the optional shift keyword argument computes the factorization of A+shift*I instead of A. If the perm argument is provided, it should
be a permutation of 1:size(A,1) giving the ordering to use (instead of CHOLMOD's default AMD ordering).
Examples
≡≡≡≡≡≡≡≡≡≡
In the following example, the fill-reducing permutation used is [3, 2, 1]. If perm is set to 1:3 to enforce no permutation, the number of
nonzero elements in the factor is 6.
julia> A = [2 1 1; 1 2 0; 1 0 2]
3×3 Matrix{Int64}:
2 1 1
1 2 0
1 0 2
julia> C = cholesky(sparse(A))
SuiteSparse.CHOLMOD.Factor{Float64}
type: LLt
method: simplicial
maxnnz: 5
nnz: 5
success: true
julia> C.p
3-element Vector{Int64}:
3
2
1
julia> L = sparse(C.L);
julia> Matrix(L)
3×3 Matrix{Float64}:
1.41421 0.0 0.0
0.0 1.41421 0.0
0.707107 0.707107 1.0
julia> L * L' ≈ A[C.p, C.p]
true
julia> P = sparse(1:3, C.p, ones(3))
3×3 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
⋅ ⋅ 1.0
⋅ 1.0 ⋅
1.0 ⋅ ⋅
julia> P' * L * L' * P ≈ A
true
julia> C = cholesky(sparse(A), perm=1:3)
SuiteSparse.CHOLMOD.Factor{Float64}
type: LLt
method: simplicial
maxnnz: 6
nnz: 6
success: true
julia> L = sparse(C.L);
julia> Matrix(L)
3×3 Matrix{Float64}:
1.41421 0.0 0.0
0.707107 1.22474 0.0
0.707107 -0.408248 1.1547
julia> L * L' ≈ A
true
│ Note
│
│ This method uses the CHOLMOD library from SuiteSparse, which only supports doubles or complex doubles. Input matrices not of those
│ element types will be converted to SparseMatrixCSC{Float64} or SparseMatrixCSC{ComplexF64} as appropriate.
│
│ Many other functions from CHOLMOD are wrapped but not exported from the Base.SparseArrays.CHOLMOD module.
And even though the docs say sparse positive definite matrix A, it seems to work for semi- as well:
julia> using ECOS
julia> model = Model(ECOS.Optimizer);
julia> @variable(model, x >= 1)
x
julia> @variable(model, y <= 2)
y
julia> @objective(model, Min, (x - y)^2)
x² - 2 y*x + y²
julia> optimize!(model)
ECOS 2.0.8 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS
It pcost dcost gap pres dres k/t mu step sigma IR | BT
0 +4.616e-01 +7.177e+00 +1e+01 4e-01 8e-01 1e+00 3e+00 --- --- 1 1 - | - -
1 -5.039e+00 +2.375e+00 +8e-01 2e-01 6e-01 5e+00 2e-01 0.9180 5e-04 1 1 1 | 0 0
2 -4.954e-01 -6.934e-02 +3e-01 4e-02 5e-02 3e-02 7e-02 0.8408 1e-01 1 2 2 | 0 0
3 -2.618e-02 -9.109e-03 +1e-02 2e-03 2e-03 3e-03 3e-03 0.9615 5e-03 1 1 1 | 0 0
4 -2.968e-04 -9.368e-05 +1e-04 2e-05 2e-05 4e-05 3e-05 0.9882 1e-04 1 1 1 | 0 0
5 -3.335e-06 -6.835e-07 +2e-06 2e-07 2e-07 5e-07 4e-07 0.9871 1e-04 2 1 1 | 0 0
6 -5.621e-08 -4.092e-10 +3e-08 5e-09 5e-09 1e-08 9e-09 0.9791 1e-04 1 1 1 | 0 0
7 -1.060e-09 +1.403e-10 +7e-10 1e-10 1e-10 3e-10 2e-10 0.9787 1e-04 1 1 1 | 0 0
OPTIMAL (within feastol=1.1e-10, reltol=7.0e-01, abstol=7.4e-10).
Runtime: 0.000410 seconds.
julia> print(backend(model).optimizer.model.model_cache)
Minimize ScalarAffineFunction{Float64}:
0.0 + 1.0 v[3]
Subject to:
VectorAffineFunction{Float64}-in-Nonnegatives
┌ ┐
│-1.0 + 1.0 x│
└ ┘ ∈ Nonnegatives(1)
┌ ┐
│0.0 - 1.0 v[2]│
└ ┘ ∈ Nonnegatives(1)
VectorAffineFunction{Float64}-in-SecondOrderCone
┌ ┐
│-2.1213203435596424 + 2.82842712474619 x - 2.82842712474619 v[2] + 0.7071067811865475 v[3]│
│3.5355339059327373 - 2.82842712474619 x + 2.82842712474619 v[2] - 0.7071067811865475 v[3] │
│0.0 + 1.4142135623730951 x - 1.414213562373095 v[2] │
│0.0 + 2.1073424255447017e-8 v[2] │
└ ┘ ∈ SecondOrderCone(4)
Do we actually have an example where Cholesky fails (due too numerics) and it should pass? And a LDL would be better?