From e0780f914caad47eba957f1ba598c4136f205e70 Mon Sep 17 00:00:00 2001
From: Oscar Dowson <o.dowson@gmail.com>
Date: Tue, 17 Mar 2026 17:02:22 +1300
Subject: [PATCH 1/5] Revert "Revert LDLFactorizations package extension
 (#2972)"

This reverts commit def69e15f8189f4a13685a24ab0a13d10bf9d92f.
---
 Project.toml                                  | 10 +++-
 ext/MathOptInterfaceLDLFactorizationsExt.jl   | 50 +++++++++++++++++++
 .../Constraint/bridges/QuadtoSOCBridge.jl     | 35 ++++++++++---
 .../Constraint/test_QuadtoSOCBridge.jl        | 38 +++++++++-----
 4 files changed, 113 insertions(+), 20 deletions(-)
 create mode 100644 ext/MathOptInterfaceLDLFactorizationsExt.jl

diff --git a/Project.toml b/Project.toml
index 6cf49918e7..3693b96bc5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -18,6 +18,12 @@ 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"
@@ -25,6 +31,7 @@ CodecZlib = "0.6, 0.7"
 ForwardDiff = "1"
 JSON = "0.21, 1"
 JSONSchema = "1"
+LDLFactorizations = "0.10"
 LinearAlgebra = "1"
 MutableArithmetics = "1"
 NaNMath = "0.3, 1"
@@ -38,8 +45,9 @@ Test = "1"
 julia = "1.10"
 
 [extras]
+LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b"
 JSONSchema = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
 ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc"
 
 [targets]
-test = ["JSONSchema", "ParallelTestRunner"]
+test = ["LDLFactorizations", "JSONSchema", "ParallelTestRunner"]
diff --git a/ext/MathOptInterfaceLDLFactorizationsExt.jl b/ext/MathOptInterfaceLDLFactorizationsExt.jl
new file mode 100644
index 0000000000..c3b5c7a31c
--- /dev/null
+++ b/ext/MathOptInterfaceLDLFactorizationsExt.jl
@@ -0,0 +1,50 @@
+# 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)
+    # Ideally we should use `LDLFactorizations.factorized(factor)` here, but it
+    # has some false negatives. Instead we check that the factorization appeared
+    # to work. This is a heuristic. There might be other cases where check is
+    # insufficient.
+    if minimum(factor.D) < 0 || any(issubnormal, factor.D)
+        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
diff --git a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl
index a4f4f0b6d6..437547c006 100644
--- a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl
+++ b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl
@@ -60,6 +60,25 @@ 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)
     # There's a big try-catch here because Cholesky can fail even if
     # `check = false`. As one example, it currently (v1.12) fails with
@@ -69,7 +88,10 @@ function compute_sparse_sqrt(Q, func, set)
     # The try-catch isn't a performance concern because the alternative is not
     # being able to reformulate the problem.
     try
-        factor = LinearAlgebra.cholesky(Q)
+        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
@@ -77,12 +99,11 @@ function compute_sparse_sqrt(Q, func, set)
         # Then, we want to permute the columns of L'. The rows stay in the same
         # order.
         return I, p[J], V
-    catch
-        msg = """
-        Unable to transform a quadratic constraint into a SecondOrderCone
-        constraint because the quadratic constraint is not strongly convex and
-        our Cholesky decomposition failed.
-        """
+    catch err
+        if err isa MOI.AddConstraintNotAllowed
+            rethrow(err)
+        end
+        msg = "There was an error computing a matrix square root"
         throw(MOI.UnsupportedConstraint{typeof(func),typeof(set)}(msg))
     end
 end
diff --git a/test/Bridges/Constraint/test_QuadtoSOCBridge.jl b/test/Bridges/Constraint/test_QuadtoSOCBridge.jl
index a32cde5480..9ca703a8e4 100644
--- a/test/Bridges/Constraint/test_QuadtoSOCBridge.jl
+++ b/test/Bridges/Constraint/test_QuadtoSOCBridge.jl
@@ -10,6 +10,7 @@ import LinearAlgebra
 import SparseArrays
 using Test
 
+import LDLFactorizations
 import MathOptInterface as MOI
 
 function runtests()
@@ -344,16 +345,13 @@ function test_semidefinite_cholesky_fail()
     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]
-    @test_throws(
-        MOI.UnsupportedConstraint,
-        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]))
+    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
 
@@ -363,8 +361,12 @@ function test_compute_sparse_sqrt_edge_cases()
         [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. LDL works
+        [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. LDL works
+        [2.0 0.0; 0.0 0.0],
     ]
         B = SparseArrays.sparse(A)
         f = zero(MOI.ScalarQuadraticFunction{eltype(A)})
@@ -379,8 +381,6 @@ function test_compute_sparse_sqrt_edge_cases()
     # Test failures
     for A in Any[
         [-1.0 0.0; 0.0 1.0],
-        [1.0 1.0; 1.0 1.0],
-        [2.0 0.0; 0.0 0.0],
         # Found from test_quadratic_nonconvex_constraint_basic
         [0.0 -1.0; -1.0 0.0],
         # Different element type. We could potentially make this work in future,
@@ -400,6 +400,20 @@ function test_compute_sparse_sqrt_edge_cases()
     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()

From bdb3c0f4297fa7d9d0269f96bc874bc006f1aea4 Mon Sep 17 00:00:00 2001
From: Richard Samuelson <samuelson.rick@gmail.com>
Date: Fri, 6 Mar 2026 01:06:10 -0500
Subject: [PATCH 2/5] Pivoted sparse Cholesky for QuadtoSOCBridge.

---
 Project.toml                                  |   8 +-
 ext/MathOptInterfaceCliqueTreesExt.jl         |  35 ++++++
 .../test_QuadtoSOCBridge_CliqueTrees.jl       | 108 ++++++++++++++++++
 3 files changed, 149 insertions(+), 2 deletions(-)
 create mode 100644 ext/MathOptInterfaceCliqueTreesExt.jl
 create mode 100644 test/Bridges/Constraint/test_QuadtoSOCBridge_CliqueTrees.jl

diff --git a/Project.toml b/Project.toml
index 3693b96bc5..e7ebd54581 100644
--- a/Project.toml
+++ b/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"]
diff --git a/ext/MathOptInterfaceCliqueTreesExt.jl b/ext/MathOptInterfaceCliqueTreesExt.jl
new file mode 100644
index 0000000000..2b13acebf2
--- /dev/null
+++ b/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
diff --git a/test/Bridges/Constraint/test_QuadtoSOCBridge_CliqueTrees.jl b/test/Bridges/Constraint/test_QuadtoSOCBridge_CliqueTrees.jl
new file mode 100644
index 0000000000..56496f6e92
--- /dev/null
+++ b/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()

From 8bcb841c479bf63467a7cd29b0f5d50f351bfbb2 Mon Sep 17 00:00:00 2001
From: Oscar Dowson <o.dowson@gmail.com>
Date: Thu, 12 Mar 2026 15:28:24 +1300
Subject: [PATCH 3/5] Refactor how extensions for sparse square root are
 handled

---
 Project.toml                                  |   4 +-
 ext/MathOptInterfaceCliqueTreesExt.jl         |  36 +++---
 ext/MathOptInterfaceLDLFactorizationsExt.jl   |  38 ++----
 .../Constraint/bridges/QuadtoSOCBridge.jl     |  63 +++++-----
 src/Utilities/Utilities.jl                    |   1 +
 src/Utilities/sparse_square_root.jl           |  72 ++++++++++++
 .../Constraint/test_QuadtoSOCBridge.jl        | 110 +++++++++++++++---
 .../test_QuadtoSOCBridge_CliqueTrees.jl       | 108 -----------------
 8 files changed, 228 insertions(+), 204 deletions(-)
 create mode 100644 src/Utilities/sparse_square_root.jl
 delete mode 100644 test/Bridges/Constraint/test_QuadtoSOCBridge_CliqueTrees.jl

diff --git a/Project.toml b/Project.toml
index e7ebd54581..1ea74639d2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,7 +4,6 @@ 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"
@@ -20,6 +19,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [weakdeps]
+CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8"
 LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b"
 
 [extensions]
@@ -38,7 +38,7 @@ LDLFactorizations = "0.10"
 LinearAlgebra = "1"
 MutableArithmetics = "1"
 NaNMath = "0.3, 1"
-OrderedCollections = "1"
+OrderedCollections = "1.1"
 ParallelTestRunner = "2.4.1"
 PrecompileTools = "1"
 Printf = "1"
diff --git a/ext/MathOptInterfaceCliqueTreesExt.jl b/ext/MathOptInterfaceCliqueTreesExt.jl
index 2b13acebf2..e1b34e58dd 100644
--- a/ext/MathOptInterfaceCliqueTreesExt.jl
+++ b/ext/MathOptInterfaceCliqueTreesExt.jl
@@ -6,30 +6,28 @@
 
 module MathOptInterfaceCliqueTreesExt
 
-import CliqueTrees.Multifrontal: ChordalCholesky
-import LinearAlgebra: RowMaximum, cholesky!
+import CliqueTrees
+import LinearAlgebra
 import MathOptInterface as MOI
-import SparseArrays: findnz, sparse
+import SparseArrays
 
-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
+MOI.Utilities.is_defined(::MOI.Utilities.CliqueTreesExt) = true
 
-    # Verify the factorization reconstructs Q. This catches indefinite
-    # matrices where the diagonal is all zeros (e.g., [0 -1; -1 0]).
+function MOI.Utilities.compute_sparse_sqrt(
+    ::MOI.Utilities.CliqueTreesExt,
+    Q::AbstractMatrix,
+)
+    G = LinearAlgebra.cholesky!(
+        CliqueTrees.Multifrontal.ChordalCholesky(Q),
+        LinearAlgebra.RowMaximum(),
+    )
+    U = SparseArrays.sparse(G.U) * G.P
+    # Verify the factorization reconstructs Q. We don't have something like
+    # LinearAlgebra.issuccess(G)
     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))
+        return nothing
     end
-
-    return findnz(U)
+    return SparseArrays.findnz(U)
 end
 
 end  # module
diff --git a/ext/MathOptInterfaceLDLFactorizationsExt.jl b/ext/MathOptInterfaceLDLFactorizationsExt.jl
index c3b5c7a31c..c1bd13ea87 100644
--- a/ext/MathOptInterfaceLDLFactorizationsExt.jl
+++ b/ext/MathOptInterfaceLDLFactorizationsExt.jl
@@ -11,39 +11,19 @@ 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(
+MOI.Utilities.is_defined(::MOI.Utilities.LDLFactorizationsExt) = true
+
+function MOI.Utilities.compute_sparse_sqrt(
+    ::MOI.Utilities.LDLFactorizationsExt,
     Q::AbstractMatrix,
-    ::F,
-    ::S,
-) where {F<:MOI.ScalarQuadraticFunction,S<:MOI.AbstractSet}
+)
     n = LinearAlgebra.checksquare(Q)
     factor = LDLFactorizations.ldl(Q)
-    # Ideally we should use `LDLFactorizations.factorized(factor)` here, but it
-    # has some false negatives. Instead we check that the factorization appeared
-    # to work. This is a heuristic. There might be other cases where check is
-    # insufficient.
-    if minimum(factor.D) < 0 || any(issubnormal, factor.D)
-        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])
+    if !LDLFactorizations.factorized(factor) || minimum(factor.D) < 0
+        return nothing
     end
-    # Finally, permute the columns of L'. The rows stay in the same order.
+    L = sqrt.(factor.D) * LinearAlgebra.UnitLowerTriangular(factor.L)
+    J, I, V = SparseArrays.findnz(SparseArrays.sparse(L))
     return I, factor.P[J], V
 end
 
diff --git a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl
index 437547c006..cbc328fe17 100644
--- a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl
+++ b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl
@@ -60,8 +60,8 @@ end
 const QuadtoSOC{T,OT<:MOI.ModelLike} =
     SingleBridgeOptimizer{QuadtoSOCBridge{T},OT}
 
-function compute_sparse_sqrt_fallback(Q, ::F, ::S) where {F,S}
-    msg = """
+function _error_msg(::MOI.Utilities.LDLFactorizationsExt)
+    return """
     Unable to transform a quadratic constraint into a SecondOrderCone
     constraint because the quadratic constraint is not strongly convex and
     our Cholesky decomposition failed.
@@ -76,36 +76,41 @@ function compute_sparse_sqrt_fallback(Q, ::F, ::S) where {F,S}
     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)
-    # There's a big try-catch here because Cholesky can fail even if
-    # `check = false`. As one example, it currently (v1.12) fails with
-    # `BigFloat`. Similarly, we want to guard against errors in
-    # `compute_sparse_sqrt_fallback`.
-    #
-    # The try-catch isn't a performance concern because the alternative is not
-    # being able to reformulate the problem.
-    try
-        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
-    catch err
-        if err isa MOI.AddConstraintNotAllowed
-            rethrow(err)
-        end
-        msg = "There was an error computing a matrix square root"
-        throw(MOI.UnsupportedConstraint{typeof(func),typeof(set)}(msg))
+function _error_msg(::MOI.Utilities.CliqueTreesExt)
+    return """
+    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 `CliqueTrees.jl`:
+    ```julia
+    import Pkg; Pkg.add("CliqueTrees")
+    using CliqueTrees
+    ```
+
+    CliqueTrees.jl is not included by default because it contains a number of
+    heavy dependencies.
+    """
+end
+
+function compute_sparse_sqrt(Q, ::F, ::S) where {F,S}
+    if (ret = MOI.Utilities.compute_sparse_sqrt(Q)) !== nothing
+        return ret
+    elseif !MOI.Utilities.is_defined(MOI.Utilities.LDLFactorizationsExt())
+        msg = _error_msg(MOI.Utilities.LDLFactorizationsExt())
+        return throw(MOI.AddConstraintNotAllowed{F,S}(msg))
+    elseif !MOI.Utilities.is_defined(MOI.Utilities.CliqueTreesExt())
+        msg = _error_msg(MOI.Utilities.CliqueTreesExt())
+        return throw(MOI.AddConstraintNotAllowed{F,S}(msg))
     end
+    msg = """
+    Unable to transform a quadratic constraint into a SecondOrderCone
+    constraint because the quadratic constraint is not convex.
+    """
+    return throw(MOI.UnsupportedConstraint{F,S}(msg))
 end
 
 function bridge_constraint(
diff --git a/src/Utilities/Utilities.jl b/src/Utilities/Utilities.jl
index 34d28429f8..505c91a79b 100644
--- a/src/Utilities/Utilities.jl
+++ b/src/Utilities/Utilities.jl
@@ -70,5 +70,6 @@ include("lazy_iterators.jl")
 include("set_dot.jl")
 
 include("distance_to_set.jl")
+include("sparse_square_root.jl")
 
 end # module
diff --git a/src/Utilities/sparse_square_root.jl b/src/Utilities/sparse_square_root.jl
new file mode 100644
index 0000000000..1cd7182c02
--- /dev/null
+++ b/src/Utilities/sparse_square_root.jl
@@ -0,0 +1,72 @@
+# 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.
+
+abstract type AbstractExt end
+
+is_defined(::AbstractExt) = false
+
+struct LDLFactorizationsExt <: AbstractExt end
+
+struct CliqueTreesExt <: AbstractExt end
+
+struct LinearAlgebraExt <: AbstractExt end
+
+is_defined(::LinearAlgebraExt) = true
+
+function compute_sparse_sqrt(::LinearAlgebraExt, Q::AbstractMatrix)
+    factor = LinearAlgebra.cholesky(Q; check = false)
+    if !LinearAlgebra.issuccess(factor)
+        return nothing
+    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
+
+"""
+    compute_sparse_sqrt(Q::AbstractMatrix)
+
+Attempts to compute a sparse square root such that `Q = A' * A`.
+
+## Return value
+
+If successful, this function returns an `(I, J, V)` triplet of the sparse `A`
+matrix.
+
+If unsuccessful, this function returns `nothing`.
+
+## Extensions
+
+By default, this function attempts to use a Cholesky decomposition. If that
+fails, it may optionally use various extension packages.
+
+These extension packages must be loaded before calling `compute_sparse_sqrt`.
+
+The extensions currently supported are:
+
+ * The LDL routine in `LDLFactorizations.jl`
+ * The pivoted Cholesky in `CliqueTrees.jl`
+"""
+function compute_sparse_sqrt(Q::AbstractMatrix)
+    # There's a big try-catch here because Cholesky can fail even if
+    # `check = false`. The try-catch isn't a performance concern because the
+    # alternative is not being able to reformulate the problem.
+    for ext in (LinearAlgebraExt(), LDLFactorizationsExt(), CliqueTreesExt())
+        if is_defined(ext)
+            try
+                if (ret = compute_sparse_sqrt(ext, Q)) !== nothing
+                    return ret
+                end
+            catch
+            end
+        end
+    end
+    return nothing
+end
diff --git a/test/Bridges/Constraint/test_QuadtoSOCBridge.jl b/test/Bridges/Constraint/test_QuadtoSOCBridge.jl
index 9ca703a8e4..21ede98582 100644
--- a/test/Bridges/Constraint/test_QuadtoSOCBridge.jl
+++ b/test/Bridges/Constraint/test_QuadtoSOCBridge.jl
@@ -6,12 +6,13 @@
 
 module TestConstraintQuadToSOC
 
-import LinearAlgebra
-import SparseArrays
 using Test
 
+import CliqueTrees
 import LDLFactorizations
+import LinearAlgebra
 import MathOptInterface as MOI
+import SparseArrays
 
 function runtests()
     for name in names(@__MODULE__; all = true)
@@ -383,11 +384,7 @@ function test_compute_sparse_sqrt_edge_cases()
         [-1.0 0.0; 0.0 1.0],
         # Found from test_quadratic_nonconvex_constraint_basic
         [0.0 -1.0; -1.0 0.0],
-        # Different element type. We could potentially make this work in future,
-        # but it first requires https://github.com/JuliaSmoothOptimizers/LDLFactorizations.jl/pull/142
         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)})
@@ -400,17 +397,96 @@ function test_compute_sparse_sqrt_edge_cases()
     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)),
-    )
+function test_ldlfactorizations_compute_sparse_sqrt_edge_cases()
+    ext = MOI.Utilities.LDLFactorizationsExt()
+    for A in AbstractMatrix[
+        [1.0 0.0; 0.0 2.0],
+        [1.0 0.0 1.0; 0.0 1.0 1.0; 1.0 1.0 3.0],
+        # [1.0 1.0; 1.0 1.0],
+        # [2.0 2.0; 2.0 2.0],
+        # [2.0 0.0; 0.0 0.0],
+    ]
+        I, J, V = MOI.Utilities.compute_sparse_sqrt(ext, SparseArrays.sparse(A))
+        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],
+        [0.0 -1.0; -1.0 0.0],
+        BigFloat[-1.0 0.0; 0.0 1.0],
+        [1.0 1.0 0.0; 1.0 1.0 0.0; 0.0 0.0 1.0],
+    ]
+        @test MOI.Utilities.compute_sparse_sqrt(ext, SparseArrays.sparse(A)) ===
+              nothing
+    end
+    return
+end
+
+function test_clique_trees_compute_sparse_sqrt_edge_cases()
+    ext = MOI.Utilities.CliqueTreesExt()
+    for A in AbstractMatrix[
+        [1.0 0.0; 0.0 2.0],
+        [1.0 0.0 1.0; 0.0 1.0 1.0; 1.0 1.0 3.0],
+        [1.0 1.0; 1.0 1.0],
+        [2.0 2.0; 2.0 2.0],
+        [2.0 0.0; 0.0 0.0],
+        [1.0 1.0 0.0; 1.0 1.0 0.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],
+    ]
+        I, J, V = MOI.Utilities.compute_sparse_sqrt(ext, SparseArrays.sparse(A))
+        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],
+        [0.0 -1.0; -1.0 0.0],
+        BigFloat[-1.0 0.0; 0.0 1.0],
+    ]
+        @test MOI.Utilities.compute_sparse_sqrt(ext, SparseArrays.sparse(A)) ===
+              nothing
+    end
+    return
+end
+
+function test_clique_trees_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_clique_trees_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 = sum(0.5 * x[i] * x[i] for i in 1:3) + 1.0 * x[1] * 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)
+    # Verify the constraint was created successfully
+    @test MOI.output_dimension(g) == 5  # [1, rhs, Ux...]
     return
 end
 
diff --git a/test/Bridges/Constraint/test_QuadtoSOCBridge_CliqueTrees.jl b/test/Bridges/Constraint/test_QuadtoSOCBridge_CliqueTrees.jl
deleted file mode 100644
index 56496f6e92..0000000000
--- a/test/Bridges/Constraint/test_QuadtoSOCBridge_CliqueTrees.jl
+++ /dev/null
@@ -1,108 +0,0 @@
-# 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()

From 109977b684e88a6a43f46d604a3c7fd49985e761 Mon Sep 17 00:00:00 2001
From: Oscar Dowson <o.dowson@gmail.com>
Date: Wed, 18 Mar 2026 13:19:20 +1300
Subject: [PATCH 4/5] Update

---
 .../Constraint/bridges/QuadtoSOCBridge.jl     | 95 ++++++++-----------
 .../Constraint/test_QuadtoSOCBridge.jl        | 59 +++---------
 2 files changed, 54 insertions(+), 100 deletions(-)

diff --git a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl
index cbc328fe17..26db3bb468 100644
--- a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl
+++ b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl
@@ -60,59 +60,6 @@ end
 const QuadtoSOC{T,OT<:MOI.ModelLike} =
     SingleBridgeOptimizer{QuadtoSOCBridge{T},OT}
 
-function _error_msg(::MOI.Utilities.LDLFactorizationsExt)
-    return """
-    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.
-    """
-end
-
-function _error_msg(::MOI.Utilities.CliqueTreesExt)
-    return """
-    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 `CliqueTrees.jl`:
-    ```julia
-    import Pkg; Pkg.add("CliqueTrees")
-    using CliqueTrees
-    ```
-
-    CliqueTrees.jl is not included by default because it contains a number of
-    heavy dependencies.
-    """
-end
-
-function compute_sparse_sqrt(Q, ::F, ::S) where {F,S}
-    if (ret = MOI.Utilities.compute_sparse_sqrt(Q)) !== nothing
-        return ret
-    elseif !MOI.Utilities.is_defined(MOI.Utilities.LDLFactorizationsExt())
-        msg = _error_msg(MOI.Utilities.LDLFactorizationsExt())
-        return throw(MOI.AddConstraintNotAllowed{F,S}(msg))
-    elseif !MOI.Utilities.is_defined(MOI.Utilities.CliqueTreesExt())
-        msg = _error_msg(MOI.Utilities.CliqueTreesExt())
-        return throw(MOI.AddConstraintNotAllowed{F,S}(msg))
-    end
-    msg = """
-    Unable to transform a quadratic constraint into a SecondOrderCone
-    constraint because the quadratic constraint is not convex.
-    """
-    return throw(MOI.UnsupportedConstraint{F,S}(msg))
-end
-
 function bridge_constraint(
     ::Type{QuadtoSOCBridge{T}},
     model,
@@ -137,8 +84,46 @@ function bridge_constraint(
             MOI.ScalarAffineTerm(scale * term.coefficient, term.variable),
         ) for term in func.affine_terms
     ]
-    I, J, V = compute_sparse_sqrt(LinearAlgebra.Symmetric(Q), func, set)
-    for (i, j, v) in zip(I, J, V)
+    sqrt_ret = MOI.Utilities.compute_sparse_sqrt(LinearAlgebra.Symmetric(Q))
+    if sqrt_ret === nothing
+        msg = """
+        ## SecondOrderCone reformulation
+
+        We tried to reformulate the quadratic constraint into a SecondOrderCone,
+        but this failed because the quadratic constraint is not strongly convex
+        and our matrix factorization failed.
+
+        ## Package extensions
+
+        If the constraint is convex but not strongly convex, you can work-around
+        this issue by manually installing and loading one of the following
+        packages.
+
+        ### LDLFactorizations.jl
+
+        Currently active: $(MOI.Utilities.is_defined(MOI.Utilities.LDLFactorizationsExt()))
+
+        ```julia
+        import Pkg; Pkg.add("LDLFactorizations")
+        using LDLFactorizations
+        ```
+        LDLFactorizations.jl is not included by default because it is licensed
+        under the LGPL.
+
+        ### CliqueTrees.jl
+
+        Currently active: $(MOI.Utilities.is_defined(MOI.Utilities.CliqueTreesExt()))
+
+        ```julia
+        import Pkg; Pkg.add("CliqueTrees")
+        using CliqueTrees
+        ```
+        CliqueTrees.jl is not included by default because it contains a number of
+        heavy dependencies.
+        """
+        return throw(MOI.UnsupportedConstraint{typeof(func),typeof(set)}(msg))
+    end
+    for (i, j, v) in zip(sqrt_ret[1], sqrt_ret[2], sqrt_ret[3])
         push!(
             vector_terms,
             MOI.VectorAffineTerm(
diff --git a/test/Bridges/Constraint/test_QuadtoSOCBridge.jl b/test/Bridges/Constraint/test_QuadtoSOCBridge.jl
index 21ede98582..3e8da39f62 100644
--- a/test/Bridges/Constraint/test_QuadtoSOCBridge.jl
+++ b/test/Bridges/Constraint/test_QuadtoSOCBridge.jl
@@ -28,32 +28,17 @@ end
 include("../utilities.jl")
 
 function test_error_for_nonconvex_quadratic_constraints()
-    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
-    bridged_mock = MOI.Bridges.Constraint.QuadtoSOC{Float64}(mock)
-    x = MOI.add_variable(bridged_mock)
+    inner = MOI.Utilities.Model{Float64}()
+    model = MOI.Bridges.Constraint.QuadtoSOC{Float64}(inner)
+    x = MOI.add_variable(model)
+    F = MOI.ScalarQuadraticFunction{Float64}
     @test_throws(
-        MOI.UnsupportedConstraint,
-        MOI.add_constraint(
-            bridged_mock,
-            MOI.ScalarQuadraticFunction(
-                [MOI.ScalarQuadraticTerm(1.0, x, x)],
-                MOI.ScalarAffineTerm{Float64}[],
-                0.0,
-            ),
-            MOI.GreaterThan(0.0),
-        )
+        MOI.UnsupportedConstraint{F,MOI.GreaterThan{Float64}},
+        MOI.add_constraint(model, 1.0 * x * x, MOI.GreaterThan(0.0))
     )
     @test_throws(
-        MOI.UnsupportedConstraint,
-        MOI.add_constraint(
-            bridged_mock,
-            MOI.ScalarQuadraticFunction(
-                [MOI.ScalarQuadraticTerm(-1.0, x, x)],
-                MOI.ScalarAffineTerm{Float64}[],
-                0.0,
-            ),
-            MOI.LessThan(0.0),
-        )
+        MOI.UnsupportedConstraint{F,MOI.LessThan{Float64}},
+        MOI.add_constraint(model, -1.0 * x * x, MOI.LessThan(0.0))
     )
     return
 end
@@ -356,23 +341,13 @@ function test_semidefinite_cholesky_fail()
     return
 end
 
-function test_compute_sparse_sqrt_edge_cases()
-    for A in Any[
-        # Trivial Cholesky
+function test_linear_algebra_compute_sparse_sqrt_edge_cases()
+    ext = MOI.Utilities.LinearAlgebraExt()
+    for A in AbstractMatrix[
         [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. LDL works
-        [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. LDL works
-        [2.0 0.0; 0.0 0.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)
+        I, J, V = MOI.Utilities.compute_sparse_sqrt(ext, SparseArrays.sparse(A))
         U = zeros(eltype(A), size(A))
         for (i, j, v) in zip(I, J, V)
             U[i, j] += v
@@ -382,17 +357,11 @@ function test_compute_sparse_sqrt_edge_cases()
     # Test failures
     for A in Any[
         [-1.0 0.0; 0.0 1.0],
-        # Found from test_quadratic_nonconvex_constraint_basic
         [0.0 -1.0; -1.0 0.0],
         BigFloat[-1.0 0.0; 0.0 1.0],
+        [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)))
-        @test_throws(
-            MOI.UnsupportedConstraint{typeof(f),typeof(s)},
-            MOI.Bridges.Constraint.compute_sparse_sqrt(B, f, s),
-        )
+        @test MOI.Utilities.compute_sparse_sqrt(ext, A) === nothing
     end
     return
 end

From 8cd38f3dc354e3e851bf7d9666ae6e7e2d412320 Mon Sep 17 00:00:00 2001
From: Oscar Dowson <o.dowson@gmail.com>
Date: Wed, 18 Mar 2026 14:17:19 +1300
Subject: [PATCH 5/5] Update

---
 .../Constraint/test_QuadtoSOCBridge.jl        | 30 +++++++++++++++++++
 1 file changed, 30 insertions(+)

diff --git a/test/Bridges/Constraint/test_QuadtoSOCBridge.jl b/test/Bridges/Constraint/test_QuadtoSOCBridge.jl
index 3e8da39f62..642af819ff 100644
--- a/test/Bridges/Constraint/test_QuadtoSOCBridge.jl
+++ b/test/Bridges/Constraint/test_QuadtoSOCBridge.jl
@@ -426,6 +426,36 @@ function test_clique_trees_compute_sparse_sqrt_edge_cases()
     return
 end
 
+function test_compute_sparse_sqrt_edge_cases()
+    for A in AbstractMatrix[
+        [1.0 0.0; 0.0 2.0],
+        [1.0 0.0 1.0; 0.0 1.0 1.0; 1.0 1.0 3.0],
+        [1.0 1.0; 1.0 1.0],
+        [2.0 2.0; 2.0 2.0],
+        [2.0 0.0; 0.0 0.0],
+        [1.0 1.0 0.0; 1.0 1.0 0.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],
+    ]
+        I, J, V = MOI.Utilities.compute_sparse_sqrt(SparseArrays.sparse(A))
+        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],
+        [0.0 -1.0; -1.0 0.0],
+        BigFloat[-1.0 0.0; 0.0 1.0],
+    ]
+        @test MOI.Utilities.compute_sparse_sqrt(SparseArrays.sparse(A)) ===
+              nothing
+    end
+    return
+end
+
 function test_clique_trees_semidefinite_cholesky_fail()
     inner = MOI.Utilities.Model{Float64}()
     model = MOI.Bridges.Constraint.QuadtoSOC{Float64}(inner)