From 9bf405b589c4cfde0f456de4e4bfb3f63fa20986 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Tue, 17 Mar 2026 17:02:22 +1300 Subject: [PATCH 1/8] Revert "Revert LDLFactorizations package extension (#2972)" This reverts commit def69e15f8189f4a13685a24ab0a13d10bf9d92f. --- Project.toml | 6 ++- ext/MathOptInterfaceLDLFactorizationsExt.jl | 50 +++++++++++++++++++ .../Constraint/bridges/QuadtoSOCBridge.jl | 35 ++++++++++--- .../Constraint/test_QuadtoSOCBridge.jl | 38 +++++++++----- 4 files changed, 109 insertions(+), 20 deletions(-) create mode 100644 ext/MathOptInterfaceLDLFactorizationsExt.jl diff --git a/Project.toml b/Project.toml index 25fe2d1f43..55037da0c7 100644 --- a/Project.toml +++ b/Project.toml @@ -19,9 +19,11 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [weakdeps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b" [extensions] MathOptInterfaceBenchmarkToolsExt = "BenchmarkTools" +MathOptInterfaceLDLFactorizationsExt = "LDLFactorizations" [compat] BenchmarkTools = "1" @@ -30,6 +32,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" @@ -44,8 +47,9 @@ julia = "1.10" [extras] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b" JSONSchema = "7d188eb4-7ad8-530c-ae41-71a32a6d4692" ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc" [targets] -test = ["BenchmarkTools", "JSONSchema", "ParallelTestRunner"] +test = ["BenchmarkTools", "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 1270f3f31f..810d6ca03f 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, @@ -399,6 +399,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 8c24eb4d68d6110931a4be98b2318de7dcaaee2f Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Fri, 6 Mar 2026 01:06:10 -0500 Subject: [PATCH 2/8] Pivoted sparse Cholesky for QuadtoSOCBridge. --- Project.toml | 11 +- ext/MathOptInterfaceCliqueTreesExt.jl | 35 ++++++ .../test_QuadtoSOCBridge_CliqueTrees.jl | 108 ++++++++++++++++++ 3 files changed, 153 insertions(+), 1 deletion(-) 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 55037da0c7..ccac8a3fbe 100644 --- a/Project.toml +++ b/Project.toml @@ -3,6 +3,11 @@ uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" version = "1.50.1" [deps] +<<<<<<< HEAD +======= +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" +>>>>>>> 26563e3e (Pivoted sparse Cholesky for QuadtoSOCBridge.) CodecBzip2 = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd" CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" @@ -23,10 +28,12 @@ LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b" [extensions] MathOptInterfaceBenchmarkToolsExt = "BenchmarkTools" +MathOptInterfaceCliqueTreesExt = "CliqueTrees" MathOptInterfaceLDLFactorizationsExt = "LDLFactorizations" [compat] BenchmarkTools = "1" +CliqueTrees = "1.17" CodecBzip2 = "0.6, 0.7, 0.8" CodecZlib = "0.6, 0.7" ForwardDiff = "1" @@ -47,9 +54,11 @@ julia = "1.10" [extras] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b" JSONSchema = "7d188eb4-7ad8-530c-ae41-71a32a6d4692" +LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b" ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc" [targets] -test = ["BenchmarkTools", "LDLFactorizations", "JSONSchema", "ParallelTestRunner"] +test = ["BenchmarkTools", "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 d9f1245a0446cb3a9801dbb1217dd6802c2dff27 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Thu, 12 Mar 2026 15:28:24 +1300 Subject: [PATCH 3/8] Refactor how extensions for sparse square root are handled --- Project.toml | 8 +- 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 | 109 +++++++++++++++--- .../test_QuadtoSOCBridge_CliqueTrees.jl | 108 ----------------- 8 files changed, 228 insertions(+), 207 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 ccac8a3fbe..6e65fc0ed6 100644 --- a/Project.toml +++ b/Project.toml @@ -3,11 +3,6 @@ uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" version = "1.50.1" [deps] -<<<<<<< HEAD -======= -BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" -CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" ->>>>>>> 26563e3e (Pivoted sparse Cholesky for QuadtoSOCBridge.) CodecBzip2 = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd" CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" @@ -24,6 +19,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [weakdeps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b" [extensions] @@ -43,7 +39,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 810d6ca03f..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,10 +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 1.0; 1.0 1.0], ] B = SparseArrays.sparse(A) f = zero(MOI.ScalarQuadraticFunction{eltype(A)}) @@ -399,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 da7e747c368c692845768a43cb229b9279ce9c17 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Wed, 18 Mar 2026 13:19:20 +1300 Subject: [PATCH 4/8] 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 876a20fa73b4a1081ac070ce7af167d5b8ac50a0 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Wed, 18 Mar 2026 14:17:19 +1300 Subject: [PATCH 5/8] 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) From 8fb6daa2c0238c0c3e3618037244f9ff4c5022de Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Fri, 27 Mar 2026 14:47:42 +1300 Subject: [PATCH 6/8] Remove LDLFactorizations --- Project.toml | 7 +- ext/MathOptInterfaceCliqueTreesExt.jl | 6 +- ext/MathOptInterfaceLDLFactorizationsExt.jl | 30 ---- .../Constraint/bridges/QuadtoSOCBridge.jl | 129 +++++++++++++----- src/Utilities/Utilities.jl | 1 - src/Utilities/sparse_square_root.jl | 72 ---------- .../Constraint/test_QuadtoSOCBridge.jl | 73 +++++----- 7 files changed, 130 insertions(+), 188 deletions(-) delete mode 100644 ext/MathOptInterfaceLDLFactorizationsExt.jl delete mode 100644 src/Utilities/sparse_square_root.jl diff --git a/Project.toml b/Project.toml index 6e65fc0ed6..f6343aa3fc 100644 --- a/Project.toml +++ b/Project.toml @@ -20,12 +20,10 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [weakdeps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" -LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b" [extensions] MathOptInterfaceBenchmarkToolsExt = "BenchmarkTools" MathOptInterfaceCliqueTreesExt = "CliqueTrees" -MathOptInterfaceLDLFactorizationsExt = "LDLFactorizations" [compat] BenchmarkTools = "1" @@ -35,7 +33,6 @@ CodecZlib = "0.6, 0.7" ForwardDiff = "1" JSON = "0.21, 1" JSONSchema = "1" -LDLFactorizations = "0.10" LinearAlgebra = "1" MutableArithmetics = "1" NaNMath = "0.3, 1" @@ -51,10 +48,8 @@ julia = "1.10" [extras] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" -LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b" JSONSchema = "7d188eb4-7ad8-530c-ae41-71a32a6d4692" -LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b" ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc" [targets] -test = ["BenchmarkTools", "CliqueTrees", "LDLFactorizations", "JSONSchema", "ParallelTestRunner"] +test = ["BenchmarkTools", "CliqueTrees", "JSONSchema", "ParallelTestRunner"] diff --git a/ext/MathOptInterfaceCliqueTreesExt.jl b/ext/MathOptInterfaceCliqueTreesExt.jl index e1b34e58dd..f0e76b9d56 100644 --- a/ext/MathOptInterfaceCliqueTreesExt.jl +++ b/ext/MathOptInterfaceCliqueTreesExt.jl @@ -11,10 +11,10 @@ import LinearAlgebra import MathOptInterface as MOI import SparseArrays -MOI.Utilities.is_defined(::MOI.Utilities.CliqueTreesExt) = true +MOI.Bridges.Constraint.is_defined(::MOI.Bridges.Constraint.CliqueTrees) = true -function MOI.Utilities.compute_sparse_sqrt( - ::MOI.Utilities.CliqueTreesExt, +function MOI.Bridges.Constraint.compute_sparse_sqrt( + ::MOI.Bridges.Constraint.CliqueTrees, Q::AbstractMatrix, ) G = LinearAlgebra.cholesky!( diff --git a/ext/MathOptInterfaceLDLFactorizationsExt.jl b/ext/MathOptInterfaceLDLFactorizationsExt.jl deleted file mode 100644 index c1bd13ea87..0000000000 --- a/ext/MathOptInterfaceLDLFactorizationsExt.jl +++ /dev/null @@ -1,30 +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 MathOptInterfaceLDLFactorizationsExt - -import LDLFactorizations -import LinearAlgebra -import MathOptInterface as MOI -import SparseArrays - -MOI.Utilities.is_defined(::MOI.Utilities.LDLFactorizationsExt) = true - -function MOI.Utilities.compute_sparse_sqrt( - ::MOI.Utilities.LDLFactorizationsExt, - Q::AbstractMatrix, -) - n = LinearAlgebra.checksquare(Q) - factor = LDLFactorizations.ldl(Q) - if !LDLFactorizations.factorized(factor) || minimum(factor.D) < 0 - return nothing - end - L = sqrt.(factor.D) * LinearAlgebra.UnitLowerTriangular(factor.L) - J, I, V = SparseArrays.findnz(SparseArrays.sparse(L)) - 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 26db3bb468..eb2f40ef53 100644 --- a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl +++ b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl @@ -60,6 +60,97 @@ end const QuadtoSOC{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{QuadtoSOCBridge{T},OT} +abstract type AbstractExt end + +is_defined(::AbstractExt) = false + +struct CliqueTrees <: 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 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(), CliqueTrees()) + if is_defined(ext) + try + if (ret = compute_sparse_sqrt(ext, Q)) !== nothing + return ret + end + catch + end + end + end + return nothing +end + +function _get_sqrt_error_message(is_clique_trees_defined::Bool) + 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. + """ + if is_clique_trees_defined + return msg + end + clique_trees = """ + + ## CliqueTrees.jl + + 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. + """ + return msg * clique_trees +end + function bridge_constraint( ::Type{QuadtoSOCBridge{T}}, model, @@ -84,43 +175,9 @@ function bridge_constraint( MOI.ScalarAffineTerm(scale * term.coefficient, term.variable), ) for term in func.affine_terms ] - sqrt_ret = MOI.Utilities.compute_sparse_sqrt(LinearAlgebra.Symmetric(Q)) + sqrt_ret = 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. - """ + msg = _get_sqrt_error_message(is_defined(CliqueTrees())) 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]) diff --git a/src/Utilities/Utilities.jl b/src/Utilities/Utilities.jl index 505c91a79b..34d28429f8 100644 --- a/src/Utilities/Utilities.jl +++ b/src/Utilities/Utilities.jl @@ -70,6 +70,5 @@ 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 deleted file mode 100644 index 1cd7182c02..0000000000 --- a/src/Utilities/sparse_square_root.jl +++ /dev/null @@ -1,72 +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. - -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 642af819ff..2a523c067e 100644 --- a/test/Bridges/Constraint/test_QuadtoSOCBridge.jl +++ b/test/Bridges/Constraint/test_QuadtoSOCBridge.jl @@ -9,7 +9,6 @@ module TestConstraintQuadToSOC using Test import CliqueTrees -import LDLFactorizations import LinearAlgebra import MathOptInterface as MOI import SparseArrays @@ -342,12 +341,13 @@ function test_semidefinite_cholesky_fail() end function test_linear_algebra_compute_sparse_sqrt_edge_cases() - ext = MOI.Utilities.LinearAlgebraExt() + ext = MOI.Bridges.Constraint.LinearAlgebraExt() 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], ] - I, J, V = MOI.Utilities.compute_sparse_sqrt(ext, SparseArrays.sparse(A)) + Q = SparseArrays.sparse(A) + I, J, V = MOI.Bridges.Constraint.compute_sparse_sqrt(ext, Q) U = zeros(eltype(A), size(A)) for (i, j, v) in zip(I, J, V) U[i, j] += v @@ -361,42 +361,13 @@ function test_linear_algebra_compute_sparse_sqrt_edge_cases() 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, A) === nothing - end - return -end - -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 + @test MOI.Bridges.Constraint.compute_sparse_sqrt(ext, A) === nothing end return end function test_clique_trees_compute_sparse_sqrt_edge_cases() - ext = MOI.Utilities.CliqueTreesExt() + ext = MOI.Bridges.Constraint.CliqueTrees() 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], @@ -407,7 +378,8 @@ function test_clique_trees_compute_sparse_sqrt_edge_cases() 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)) + Q = SparseArrays.sparse(A) + I, J, V = MOI.Bridges.Constraint.compute_sparse_sqrt(ext, Q) U = zeros(eltype(A), size(A)) for (i, j, v) in zip(I, J, V) U[i, j] += v @@ -420,8 +392,8 @@ function test_clique_trees_compute_sparse_sqrt_edge_cases() [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 + Q = SparseArrays.sparse(A) + @test MOI.Bridges.Constraint.compute_sparse_sqrt(ext, Q) === nothing end return end @@ -437,7 +409,8 @@ function test_compute_sparse_sqrt_edge_cases() 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)) + Q = SparseArrays.sparse(A) + I, J, V = MOI.Bridges.Constraint.compute_sparse_sqrt(Q) U = zeros(eltype(A), size(A)) for (i, j, v) in zip(I, J, V) U[i, j] += v @@ -450,8 +423,8 @@ function test_compute_sparse_sqrt_edge_cases() [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 + Q = SparseArrays.sparse(A) + @test MOI.Bridges.Constraint.compute_sparse_sqrt(Q) === nothing end return end @@ -489,6 +462,26 @@ function test_clique_trees_early_zero_pivot() return end +function test_clique_trees_error_message() + for flag in (true, false) + msg = MOI.Bridges.Constraint._get_sqrt_error_message(flag) + @test occursin("CliqueTrees", msg) == !flag + end + return +end + +struct _DummyCliqueTrees <: MOI.Bridges.Constraint.AbstractExt end + +function test_is_defined_default_fallback() + @test !MOI.Bridges.Constraint.is_defined(_DummyCliqueTrees()) + Q = ones(2, 2) + @test_throws( + MethodError, + MOI.Bridges.Constraint.compute_sparse_sqrt(_DummyCliqueTrees(), Q), + ) + return +end + end # module TestConstraintQuadToSOC.runtests() From 99ae2ff113971c0734de2e7b7552db6344938a67 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Fri, 27 Mar 2026 16:12:59 +1300 Subject: [PATCH 7/8] Update --- ext/MathOptInterfaceCliqueTreesExt.jl | 6 ++--- .../Constraint/bridges/QuadtoSOCBridge.jl | 24 +++++++++---------- .../Constraint/test_QuadtoSOCBridge.jl | 20 ++++++++-------- 3 files changed, 25 insertions(+), 25 deletions(-) diff --git a/ext/MathOptInterfaceCliqueTreesExt.jl b/ext/MathOptInterfaceCliqueTreesExt.jl index f0e76b9d56..75cad5b44f 100644 --- a/ext/MathOptInterfaceCliqueTreesExt.jl +++ b/ext/MathOptInterfaceCliqueTreesExt.jl @@ -11,10 +11,10 @@ import LinearAlgebra import MathOptInterface as MOI import SparseArrays -MOI.Bridges.Constraint.is_defined(::MOI.Bridges.Constraint.CliqueTrees) = true +MOI.Bridges.Constraint.is_defined(::MOI.Bridges.Constraint._CliqueTrees) = true -function MOI.Bridges.Constraint.compute_sparse_sqrt( - ::MOI.Bridges.Constraint.CliqueTrees, +function MOI.Bridges.Constraint._compute_sparse_sqrt( + ::MOI.Bridges.Constraint._CliqueTrees, Q::AbstractMatrix, ) G = LinearAlgebra.cholesky!( diff --git a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl index eb2f40ef53..895af0cc7a 100644 --- a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl +++ b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl @@ -60,17 +60,17 @@ end const QuadtoSOC{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{QuadtoSOCBridge{T},OT} -abstract type AbstractExt end +abstract type _AbstractExt end -is_defined(::AbstractExt) = false +is_defined(::_AbstractExt) = false -struct CliqueTrees <: AbstractExt end +struct _CliqueTrees <: _AbstractExt end -struct LinearAlgebraExt <: AbstractExt end +struct _LinearAlgebra <: _AbstractExt end -is_defined(::LinearAlgebraExt) = true +is_defined(::_LinearAlgebra) = true -function compute_sparse_sqrt(::LinearAlgebraExt, Q::AbstractMatrix) +function _compute_sparse_sqrt(::_LinearAlgebra, Q::AbstractMatrix) factor = LinearAlgebra.cholesky(Q; check = false) if !LinearAlgebra.issuccess(factor) return nothing @@ -85,7 +85,7 @@ function compute_sparse_sqrt(::LinearAlgebraExt, Q::AbstractMatrix) end """ - compute_sparse_sqrt(Q::AbstractMatrix) + _compute_sparse_sqrt(Q::AbstractMatrix) Attempts to compute a sparse square root such that `Q = A' * A`. @@ -101,20 +101,20 @@ If unsuccessful, this function returns `nothing`. 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`. +These extension packages must be loaded before calling `_compute_sparse_sqrt`. The extensions currently supported are: * The pivoted Cholesky in `CliqueTrees.jl` """ -function compute_sparse_sqrt(Q::AbstractMatrix) +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(), CliqueTrees()) + for ext in (_LinearAlgebra(), _CliqueTrees()) if is_defined(ext) try - if (ret = compute_sparse_sqrt(ext, Q)) !== nothing + if (ret = _compute_sparse_sqrt(ext, Q)) !== nothing return ret end catch @@ -175,7 +175,7 @@ function bridge_constraint( MOI.ScalarAffineTerm(scale * term.coefficient, term.variable), ) for term in func.affine_terms ] - sqrt_ret = compute_sparse_sqrt(LinearAlgebra.Symmetric(Q)) + sqrt_ret = _compute_sparse_sqrt(LinearAlgebra.Symmetric(Q)) if sqrt_ret === nothing msg = _get_sqrt_error_message(is_defined(CliqueTrees())) return throw(MOI.UnsupportedConstraint{typeof(func),typeof(set)}(msg)) diff --git a/test/Bridges/Constraint/test_QuadtoSOCBridge.jl b/test/Bridges/Constraint/test_QuadtoSOCBridge.jl index 2a523c067e..d51737d273 100644 --- a/test/Bridges/Constraint/test_QuadtoSOCBridge.jl +++ b/test/Bridges/Constraint/test_QuadtoSOCBridge.jl @@ -341,13 +341,13 @@ function test_semidefinite_cholesky_fail() end function test_linear_algebra_compute_sparse_sqrt_edge_cases() - ext = MOI.Bridges.Constraint.LinearAlgebraExt() + ext = MOI.Bridges.Constraint._LinearAlgebra() 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], ] Q = SparseArrays.sparse(A) - I, J, V = MOI.Bridges.Constraint.compute_sparse_sqrt(ext, Q) + I, J, V = MOI.Bridges.Constraint._compute_sparse_sqrt(ext, Q) U = zeros(eltype(A), size(A)) for (i, j, v) in zip(I, J, V) U[i, j] += v @@ -361,13 +361,13 @@ function test_linear_algebra_compute_sparse_sqrt_edge_cases() 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.Bridges.Constraint.compute_sparse_sqrt(ext, A) === nothing + @test MOI.Bridges.Constraint._compute_sparse_sqrt(ext, A) === nothing end return end function test_clique_trees_compute_sparse_sqrt_edge_cases() - ext = MOI.Bridges.Constraint.CliqueTrees() + ext = MOI.Bridges.Constraint._CliqueTrees() 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], @@ -379,7 +379,7 @@ function test_clique_trees_compute_sparse_sqrt_edge_cases() BigFloat[1.0 1.0; 1.0 1.0], ] Q = SparseArrays.sparse(A) - I, J, V = MOI.Bridges.Constraint.compute_sparse_sqrt(ext, Q) + I, J, V = MOI.Bridges.Constraint._compute_sparse_sqrt(ext, Q) U = zeros(eltype(A), size(A)) for (i, j, v) in zip(I, J, V) U[i, j] += v @@ -393,7 +393,7 @@ function test_clique_trees_compute_sparse_sqrt_edge_cases() BigFloat[-1.0 0.0; 0.0 1.0], ] Q = SparseArrays.sparse(A) - @test MOI.Bridges.Constraint.compute_sparse_sqrt(ext, Q) === nothing + @test MOI.Bridges.Constraint._compute_sparse_sqrt(ext, Q) === nothing end return end @@ -410,7 +410,7 @@ function test_compute_sparse_sqrt_edge_cases() BigFloat[1.0 1.0; 1.0 1.0], ] Q = SparseArrays.sparse(A) - I, J, V = MOI.Bridges.Constraint.compute_sparse_sqrt(Q) + I, J, V = MOI.Bridges.Constraint._compute_sparse_sqrt(Q) U = zeros(eltype(A), size(A)) for (i, j, v) in zip(I, J, V) U[i, j] += v @@ -424,7 +424,7 @@ function test_compute_sparse_sqrt_edge_cases() BigFloat[-1.0 0.0; 0.0 1.0], ] Q = SparseArrays.sparse(A) - @test MOI.Bridges.Constraint.compute_sparse_sqrt(Q) === nothing + @test MOI.Bridges.Constraint._compute_sparse_sqrt(Q) === nothing end return end @@ -470,14 +470,14 @@ function test_clique_trees_error_message() return end -struct _DummyCliqueTrees <: MOI.Bridges.Constraint.AbstractExt end +struct _DummyCliqueTrees <: MOI.Bridges.Constraint._AbstractExt end function test_is_defined_default_fallback() @test !MOI.Bridges.Constraint.is_defined(_DummyCliqueTrees()) Q = ones(2, 2) @test_throws( MethodError, - MOI.Bridges.Constraint.compute_sparse_sqrt(_DummyCliqueTrees(), Q), + MOI.Bridges.Constraint._compute_sparse_sqrt(_DummyCliqueTrees(), Q), ) return end From 8049e9f831eebffdcf9329ee296625d29722f078 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Fri, 27 Mar 2026 16:27:17 +1300 Subject: [PATCH 8/8] Update --- src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl index 895af0cc7a..58fbd1f538 100644 --- a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl +++ b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl @@ -177,7 +177,7 @@ function bridge_constraint( ] sqrt_ret = _compute_sparse_sqrt(LinearAlgebra.Symmetric(Q)) if sqrt_ret === nothing - msg = _get_sqrt_error_message(is_defined(CliqueTrees())) + msg = _get_sqrt_error_message(is_defined(_CliqueTrees())) 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])