From cf9724b9fc63ed903893fdc390d2e1d5575ab6da Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Fri, 6 Mar 2026 01:06:10 -0500 Subject: [PATCH] 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()