From cf9724b9fc63ed903893fdc390d2e1d5575ab6da Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Fri, 6 Mar 2026 01:06:10 -0500 Subject: [PATCH 1/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 db4a366813b821a4e738a9f2038f1c1976c84f12 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Thu, 12 Mar 2026 15:28:24 +1300 Subject: [PATCH 2/5] Refactor how extensions for sparse square root are handled --- Project.toml | 2 +- ext/MathOptInterfaceCliqueTreesExt.jl | 32 +++-- 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, 225 insertions(+), 201 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..68383bd483 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] diff --git a/ext/MathOptInterfaceCliqueTreesExt.jl b/ext/MathOptInterfaceCliqueTreesExt.jl index 2b13acebf2..093c71ad57 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 +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. 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)) + 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..178e11177e 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..768211e788 --- /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 `LinearAlgebra.cholesky`. 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..7b32b460b7 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 = 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 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 60e83ed22d4595f945f006c8381a32da84dbf12d Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Fri, 13 Mar 2026 15:02:18 +1300 Subject: [PATCH 3/5] Apply suggestions from code review Co-authored-by: Oscar Dowson --- Project.toml | 2 +- src/Utilities/sparse_square_root.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 68383bd483..645b8fc854 100644 --- a/Project.toml +++ b/Project.toml @@ -28,7 +28,7 @@ MathOptInterfaceLDLFactorizationsExt = "LDLFactorizations" [compat] BenchmarkTools = "1" -CliqueTrees = "1.17" +CliqueTrees = "1.13" CodecBzip2 = "0.6, 0.7, 0.8" CodecZlib = "0.6, 0.7" ForwardDiff = "1" diff --git a/src/Utilities/sparse_square_root.jl b/src/Utilities/sparse_square_root.jl index 768211e788..1cd7182c02 100644 --- a/src/Utilities/sparse_square_root.jl +++ b/src/Utilities/sparse_square_root.jl @@ -44,7 +44,7 @@ If unsuccessful, this function returns `nothing`. ## Extensions -By default, this function attempts to use a `LinearAlgebra.cholesky`. If that +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`. @@ -52,7 +52,7 @@ 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` + * 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 From 3661f0ef03dec450f5493c912612b1c92762613a Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Fri, 13 Mar 2026 15:34:55 +1300 Subject: [PATCH 4/5] Update CliqueTrees and OrderedCollections versions --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 645b8fc854..1ea74639d2 100644 --- a/Project.toml +++ b/Project.toml @@ -28,7 +28,7 @@ MathOptInterfaceLDLFactorizationsExt = "LDLFactorizations" [compat] BenchmarkTools = "1" -CliqueTrees = "1.13" +CliqueTrees = "1.17" CodecBzip2 = "0.6, 0.7, 0.8" CodecZlib = "0.6, 0.7" ForwardDiff = "1" @@ -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" From efe77a92cabb9f364e76da8f7ef83c822eb253e4 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Fri, 13 Mar 2026 17:07:23 +1300 Subject: [PATCH 5/5] Update --- ext/MathOptInterfaceCliqueTreesExt.jl | 4 ++-- src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl | 4 ++-- test/Bridges/Constraint/test_QuadtoSOCBridge.jl | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/ext/MathOptInterfaceCliqueTreesExt.jl b/ext/MathOptInterfaceCliqueTreesExt.jl index 093c71ad57..e1b34e58dd 100644 --- a/ext/MathOptInterfaceCliqueTreesExt.jl +++ b/ext/MathOptInterfaceCliqueTreesExt.jl @@ -22,8 +22,8 @@ function MOI.Utilities.compute_sparse_sqrt( LinearAlgebra.RowMaximum(), ) U = SparseArrays.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]). + # Verify the factorization reconstructs Q. We don't have something like + # LinearAlgebra.issuccess(G) if !isapprox(Q, U' * U; atol = 1e-10) return nothing end diff --git a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl index 178e11177e..cbc328fe17 100644 --- a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl +++ b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl @@ -100,10 +100,10 @@ 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()) + 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()) + msg = _error_msg(MOI.Utilities.CliqueTreesExt()) return throw(MOI.AddConstraintNotAllowed{F,S}(msg)) end msg = """ diff --git a/test/Bridges/Constraint/test_QuadtoSOCBridge.jl b/test/Bridges/Constraint/test_QuadtoSOCBridge.jl index 7b32b460b7..21ede98582 100644 --- a/test/Bridges/Constraint/test_QuadtoSOCBridge.jl +++ b/test/Bridges/Constraint/test_QuadtoSOCBridge.jl @@ -480,7 +480,7 @@ function test_clique_trees_early_zero_pivot() 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] + 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}()))