Pivoted sparse Cholesky for QuadtoSOCBridge.

Project.toml

@@ -4,6 +4,7 @@ version = "1.49.0"
 
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8"
 CodecBzip2 = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
 CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -22,10 +23,12 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b"
 
 [extensions]
+MathOptInterfaceCliqueTreesExt = "CliqueTrees"
 MathOptInterfaceLDLFactorizationsExt = "LDLFactorizations"
 
 [compat]
 BenchmarkTools = "1"
+CliqueTrees = "1.17"
 CodecBzip2 = "0.6, 0.7, 0.8"
 CodecZlib = "0.6, 0.7"
 ForwardDiff = "1"
@@ -45,9 +48,10 @@ Test = "1"
 julia = "1.10"
 
 [extras]
-LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b"
+CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8"
 JSONSchema = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
+LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b"
 ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc"
 
 [targets]
-test = ["LDLFactorizations", "JSONSchema", "ParallelTestRunner"]
+test = ["CliqueTrees", "LDLFactorizations", "JSONSchema", "ParallelTestRunner"]

ext/MathOptInterfaceCliqueTreesExt.jl

@@ -0,0 +1,35 @@
+# 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 MathOptInterfaceCliqueTreesExt
+
+import CliqueTrees.Multifrontal: ChordalCholesky
+import LinearAlgebra: RowMaximum, cholesky!
+import MathOptInterface as MOI
+import SparseArrays: findnz, sparse
+
+function MOI.Bridges.Constraint.compute_sparse_sqrt_fallback(
+    Q::AbstractMatrix,
+    ::F,
+    ::S,
+) where {F<:MOI.ScalarQuadraticFunction,S<:MOI.AbstractSet}
+    G = cholesky!(ChordalCholesky(Q), RowMaximum())
+    U = sparse(G.U) * G.P
+
+    # Verify the factorization reconstructs Q. This catches indefinite
+    # matrices where the diagonal is all zeros (e.g., [0 -1; -1 0]).
+    if !isapprox(Q, U' * U; atol = 1e-10)
+        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
+
+    return findnz(U)
+end
+
+end  # module

test/Bridges/Constraint/test_QuadtoSOCBridge_CliqueTrees.jl

@@ -0,0 +1,108 @@
+# 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 TestConstraintQuadToSOCCliqueTrees
+
+import SparseArrays
+using Test
+
+import CliqueTrees
+import MathOptInterface as MOI
+
+function runtests()
+    for name in names(@__MODULE__; all = true)
+        if startswith("$(name)", "test_")
+            @testset "$(name)" begin
+                getfield(@__MODULE__, name)()
+            end
+        end
+    end
+    return
+end
+
+function test_compute_sparse_sqrt_edge_cases()
+    for A in Any[
+        # 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 eigenvalue. Pivoted Cholesky works.
+        [1.0 1.0; 1.0 1.0],
+        # Cholesky succeeds, even though 0 eigenvalue
+        [2.0 2.0; 2.0 2.0],
+        # Cholesky fails because of 0 column/row. Pivoted Cholesky works.
+        [2.0 0.0; 0.0 0.0],
+        # Early zero pivot - this case breaks LDLFactorizations but works
+        # with CliqueTrees' pivoted Cholesky.
+        [1.0 1.0 0.0; 1.0 1.0 0.0; 0.0 0.0 1.0],
+    ]
+        B = SparseArrays.sparse(A)
+        f = zero(MOI.ScalarQuadraticFunction{eltype(A)})
+        s = MOI.GreaterThan(zero(eltype(A)))
+        I, J, V = MOI.Bridges.Constraint.compute_sparse_sqrt(B, f, s)
+        U = zeros(eltype(A), 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
+    # Test failures
+    for A in Any[
+        [-1.0 0.0; 0.0 1.0],
+        # Indefinite matrix
+        [0.0 -1.0; -1.0 0.0],
+        # BigFloat not supported
+        BigFloat[-1.0 0.0; 0.0 1.0],
+        BigFloat[1.0 0.0; 0.0 2.0],
+        BigFloat[1.0 1.0; 1.0 1.0],
+    ]
+        B = SparseArrays.sparse(A)
+        f = zero(MOI.ScalarQuadraticFunction{eltype(A)})
+        s = MOI.GreaterThan(zero(eltype(A)))
+        @test_throws(
+            MOI.UnsupportedConstraint{typeof(f),typeof(s)},
+            MOI.Bridges.Constraint.compute_sparse_sqrt(B, f, s),
+        )
+    end
+    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_early_zero_pivot()
+    # This matrix has an early zero pivot that causes LDLFactorizations to
+    # halt early, but CliqueTrees' pivoted Cholesky handles it correctly.
+    inner = MOI.Utilities.Model{Float64}()
+    model = MOI.Bridges.Constraint.QuadtoSOC{Float64}(inner)
+    x = MOI.add_variables(model, 3)
+    # (x[1] + x[2])^2 + x[3]^2 = x[1]^2 + 2*x[1]*x[2] + x[2]^2 + x[3]^2
+    # Q = [1 1 0; 1 1 0; 0 0 1]
+    f = 0.5 * x[1] * x[1] + 1.0 * x[1] * x[2] + 0.5 * x[2] * x[2] + 0.5 * x[3] * x[3]
+    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)
+    # Verify the constraint was created successfully
+    @test MOI.output_dimension(g) == 5  # [1, rhs, Ux...]
+    return
+end
+
+end  # module
+
+TestConstraintQuadToSOCCliqueTrees.runtests()