Sparse LDL for QuadtoSOC

[mlubin]

QuadtoSOC computes a cholesky factorization, which fails when the Q matrix is positive semidefinite but not positive definite:

MathOptInterface.jl/src/Bridges/Constraint/bridges/quad_to_soc.jl

Lines 71 to 85 in 8c8636f

	F = try
	LinearAlgebra.cholesky(LinearAlgebra.Symmetric(Q))
	catch
	throw(
	MOI.UnsupportedConstraint{typeof(func),typeof(set)}(
	"Unable to transform a quadratic constraint into a " *
	"second-order cone constraint because the quadratic " *
	"constraint is not strongly convex.\n\nConvex constraints " *
	"that are not strongly convex (i.e., the matrix is positive " *
	"semidefinite but not positive definite) are not supported " *
	"yet.\n\nNote that a quadratic equality constraint is " *
	"non-convex.",
	),
	)
	end

Per @dpo we could use LDLFactorizations with a pivot tolerance of zero (or HSL/MUMPS, if we want external binary dependencies).

[mlubin]

LDLFactorizations is LGPL, so we can't use it in MOI 😢

[abelsiqueira]

Can you expand on why? If you are just using the package, does it affect your license?

[mlubin]

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.

[dpo]

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.

[abelsiqueira]

Complicated stuff. I did not realize that Julia made LGPL ambiguous, but it kind of makes sense.

[mlubin]

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.

[odow]

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?

[mlubin]

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.

[jd-foster]

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)).

[mlubin]

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.

[mlubin]

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.

[odow]

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?