Add LDLFactorizations.jl as a package extension

Project.toml

@@ -18,13 +18,20 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
+[weakdeps]
+LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b"
+
+[extensions]
+MathOptInterfaceLDLFactorizationsExt = "LDLFactorizations"
+
 [compat]
 BenchmarkTools = "1"
 CodecBzip2 = "0.6, 0.7, 0.8"
 CodecZlib = "0.6, 0.7"
 ForwardDiff = "1"
 JSON = "0.21, 1"
 JSONSchema = "1"
+LDLFactorizations = "0.10"
 LinearAlgebra = "1"
 MutableArithmetics = "1"
 NaNMath = "0.3, 1"
@@ -37,7 +44,8 @@ Test = "1"
 julia = "1.10"
 
 [extras]
+LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b"
 JSONSchema = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
 
 [targets]
-test = ["JSONSchema"]
+test = ["LDLFactorizations", "JSONSchema"]

ext/MathOptInterfaceLDLFactorizationsExt.jl

@@ -0,0 +1,46 @@
+# Copyright (c) 2017: Miles Lubin and contributors
+# Copyright (c) 2017: Google Inc.
+#
+# Use of this source code is governed by an MIT-style license that can be found
+# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
+
+module MathOptInterfaceLDLFactorizationsExt
+
+import LDLFactorizations
+import LinearAlgebra
+import MathOptInterface as MOI
+import SparseArrays
+
+# The type signature of this function is not important, so long as it is more
+# specific than the (untyped) generic fallback with the error pointing to
+# LDLFactorizations.jl
+function MOI.Bridges.Constraint.compute_sparse_sqrt_fallback(
+    Q::AbstractMatrix,
+    ::F,
+    ::S,
+) where {F<:MOI.ScalarQuadraticFunction,S<:MOI.AbstractSet}
+    n = LinearAlgebra.checksquare(Q)
+    factor = LDLFactorizations.ldl(Q)
+    if minimum(factor.D) < 0
+        msg = """
+        Unable to transform a quadratic constraint into a SecondOrderCone
+        constraint because the quadratic constraint is not convex.
+        """
+        throw(MOI.UnsupportedConstraint{F,S}(msg))
+    end
+    # We have Q = P' * L * D * L' * P. We want to find Q = U' * U, so
+    # U = sqrt(D) * L' * P. First, compute L'. Note I and J are reversed:
+    J, I, V = SparseArrays.findnz(factor.L)
+    # Except L doesn't include the identity along the diagonal. Add it back.
+    append!(J, 1:n)
+    append!(I, 1:n)
+    append!(V, ones(n))
+    # Now scale by sqrt(D)
+    for (k, i) in enumerate(I)
+        V[k] *= sqrt(factor.D[i, i])
+    end
+    # Finally, permute the columns of L'. The rows stay in the same order.
+    return I, factor.P[J], V
+end
+
+end  # module

src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl

@@ -60,6 +60,39 @@ end
 const QuadtoSOC{T,OT<:MOI.ModelLike} =
     SingleBridgeOptimizer{QuadtoSOCBridge{T},OT}
 
+function compute_sparse_sqrt_fallback(Q, ::F, ::S) where {F,S}
+    msg = """
+    Unable to transform a quadratic constraint into a SecondOrderCone
+    constraint because the quadratic constraint is not strongly convex and
+    our Cholesky decomposition failed.
+
+    If the constraint is convex but not strongly convex, you can work-around
+    this issue by manually installing and loading `LDLFactorizations.jl`:
+    ```julia
+    import Pkg; Pkg.add("LDLFactorizations")
+    using LDLFactorizations
+    ```
+
+    LDLFactorizations.jl is not included by default because it is licensed
+    under the LGPL.
+    """
+    return throw(MOI.AddConstraintNotAllowed{F,S}(msg))
+end
+
+function compute_sparse_sqrt(Q, func, set)
+    factor = LinearAlgebra.cholesky(Q; check = false)
+    if !LinearAlgebra.issuccess(factor)
+        return compute_sparse_sqrt_fallback(Q, func, set)
+    end
+    L, p = SparseArrays.sparse(factor.L), factor.p
+    # We have Q = P' * L * L' * P. We want to find Q = U' * U, so U = L' * P
+    # First, compute L'. Note I and J are reversed
+    J, I, V = SparseArrays.findnz(L)
+    # Then, we want to permute the columns of L'. The rows stay in the same
+    # order.
+    return I, p[J], V
+end
+
 function bridge_constraint(
     ::Type{QuadtoSOCBridge{T}},
     model,
@@ -73,21 +106,6 @@ function bridge_constraint(
     if !less_than
         LinearAlgebra.rmul!(Q, -1)
     end
-    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 (that is, the matrix is positive " *
-                "semidefinite but not positive definite) are not supported " *
-                "yet.\n\nNote that a quadratic equality constraint is " *
-                "non-convex.",
-            ),
-        )
-    end
     # Construct the VectorAffineFunction. We're aiming for:
     #  |          1  |
     #  | -a^T x + ub | ∈ RotatedSecondOrderCone()
@@ -99,24 +117,20 @@ function bridge_constraint(
             MOI.ScalarAffineTerm(scale * term.coefficient, term.variable),
         ) for term in func.affine_terms
     ]
-    # Note that `L = U'` so we add the terms of L'x
-    L, p = SparseArrays.sparse(F.L), F.p
-    I, J, V = SparseArrays.findnz(L)
-    for i in 1:length(V)
-        # Cholesky is a pivoted decomposition, so L × L' == Q[p, p].
-        # To get the variable, we need the row of L, I[i], then to map that
-        # through the permutation vector, so `p[I[i]]`, and then through the
-        # index_to_variable_map.
-        xi = index_to_variable_map[p[I[i]]]
+    I, J, V = compute_sparse_sqrt(LinearAlgebra.Symmetric(Q), func, set)
+    for (i, j, v) in zip(I, J, V)
         push!(
             vector_terms,
-            MOI.VectorAffineTerm(J[i] + 2, MOI.ScalarAffineTerm(V[i], xi)),
+            MOI.VectorAffineTerm(
+                i + 2,
+                MOI.ScalarAffineTerm(v, index_to_variable_map[j]),
+            ),
         )
     end
     # This is the [1, ub, 0] vector...
     set_constant = MOI.constant(set)
     MOI.throw_if_scalar_and_constant_not_zero(func, typeof(set))
-    vector_constant = vcat(one(T), -scale * set_constant, zeros(T, size(L, 1)))
+    vector_constant = vcat(one(T), -scale * set_constant, zeros(T, size(Q, 1)))
     f = MOI.VectorAffineFunction(vector_terms, vector_constant)
     dimension = MOI.output_dimension(f)
     soc = MOI.add_constraint(model, f, MOI.RotatedSecondOrderCone(dimension))

test/Bridges/Constraint/QuadtoSOCBridge.jl

@@ -10,6 +10,7 @@ import LinearAlgebra
 import SparseArrays
 using Test
 
+import LDLFactorizations
 import MathOptInterface as MOI
 
 function runtests()
@@ -339,6 +340,67 @@ function test_constraint_primal_no_quad_terms()
     return
 end
 
+function test_semidefinite_cholesky_fail()
+    inner = MOI.Utilities.Model{Float64}()
+    model = MOI.Bridges.Constraint.QuadtoSOC{Float64}(inner)
+    x = MOI.add_variables(model, 2)
+    f = 0.5 * x[1] * x[1] + 1.0 * x[1] * x[2] + 0.5 * x[2] * x[2]
+    c = MOI.add_constraint(model, f, MOI.LessThan(1.0))
+    F, S = MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone
+    ci = only(MOI.get(inner, MOI.ListOfConstraintIndices{F,S}()))
+    g = MOI.get(inner, MOI.ConstraintFunction(), ci)
+    y = MOI.get(inner, MOI.ListOfVariableIndices())
+    sum_y = 1.0 * y[1] + 1.0 * y[2]
+    @test isapprox(g, MOI.Utilities.vectorize([1.0, 1.0, sum_y, 0.0]))
+    return
+end
+
+function test_compute_sparse_sqrt_edge_cases()
+    f = zero(MOI.ScalarQuadraticFunction{Float64})
+    s = MOI.GreaterThan(0.0)
+    for A in [
+        # Trivial Cholesky
+        [1.0 0.0; 0.0 2.0],
+        # Cholesky works, with pivoting
+        [1.0 0.0 1.0; 0.0 1.0 1.0; 1.0 1.0 3.0],
+        # Cholesky fails due to 0 eigen value
+        [1.0 1.0; 1.0 1.0],
+        # Cholesky succeeds, even though 0 eigen value
+        [2.0 2.0; 2.0 2.0],
+        # Cholesky fails because of 0 column/row
+        [2.0 0.0; 0.0 0.0],
+    ]
+        B = SparseArrays.sparse(A)
+        I, J, V = MOI.Bridges.Constraint.compute_sparse_sqrt(B, f, s)
+        U = zeros(size(A))
+        for (i, j, v) in zip(I, J, V)
+            U[i, j] += v
+        end
+        @test isapprox(A, U' * U; atol = 1e-10)
+    end
+    A = [-1.0 0.0; 0.0 1.0]
+    B = SparseArrays.sparse(A)
+    @test_throws(
+        MOI.UnsupportedConstraint{typeof(f),typeof(s)},
+        MOI.Bridges.Constraint.compute_sparse_sqrt(B, f, s),
+    )
+    return
+end
+
+function test_compute_sparse_sqrt_fallback()
+    # Test the default fallback that is hit when LDLFactorizations isn't loaded.
+    # We could put the test somewhere else so it runs before this file is
+    # loaded, but that's pretty flakey for a long-term solution. Instead, we're
+    # going to abuse the lack of a strong type signature to hit it:
+    f = zero(MOI.ScalarAffineFunction{Float64})
+    A = SparseArrays.sparse([-1.0 0.0; 0.0 1.0])
+    @test_throws(
+        MOI.AddConstraintNotAllowed{typeof(f),MOI.GreaterThan{Float64}},
+        MOI.Bridges.Constraint.compute_sparse_sqrt(A, f, MOI.GreaterThan(0.0)),
+    )
+    return
+end
+
 end  # module
 
 TestConstraintQuadToSOC.runtests()