From 65d9193c142074701dca76f67d33bf7f60f60657 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Mon, 2 Feb 2026 08:46:57 +1300 Subject: [PATCH 1/3] Add LDLFactorizations.jl as a package extension --- Project.toml | 10 ++- ext/MathOptInterfaceLDLFactorizationsExt.jl | 46 +++++++++++++ .../Constraint/bridges/QuadtoSOCBridge.jl | 66 +++++++++++-------- test/Bridges/Constraint/QuadtoSOCBridge.jl | 48 ++++++++++++++ 4 files changed, 143 insertions(+), 27 deletions(-) create mode 100644 ext/MathOptInterfaceLDLFactorizationsExt.jl diff --git a/Project.toml b/Project.toml index 20f2b66f3e..c1da77ff95 100644 --- a/Project.toml +++ b/Project.toml @@ -18,6 +18,12 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +[weakdeps] +LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b" + +[extensions] +MathOptInterfaceLDLFactorizationsExt = "LDLFactorizations" + [compat] BenchmarkTools = "1" CodecBzip2 = "0.6, 0.7, 0.8" @@ -25,6 +31,7 @@ CodecZlib = "0.6, 0.7" ForwardDiff = "1" JSON = "0.21, 1" JSONSchema = "1" +LDLFactorizations = "0.10" LinearAlgebra = "1" MutableArithmetics = "1" NaNMath = "0.3, 1" @@ -37,7 +44,8 @@ Test = "1" julia = "1.10" [extras] +LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b" JSONSchema = "7d188eb4-7ad8-530c-ae41-71a32a6d4692" [targets] -test = ["JSONSchema"] +test = ["LDLFactorizations", "JSONSchema"] diff --git a/ext/MathOptInterfaceLDLFactorizationsExt.jl b/ext/MathOptInterfaceLDLFactorizationsExt.jl new file mode 100644 index 0000000000..7eeb8ad9d6 --- /dev/null +++ b/ext/MathOptInterfaceLDLFactorizationsExt.jl @@ -0,0 +1,46 @@ +# 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_root_fallback( + Q::AbstractMatrix, + ::F, + ::S, +) where {F<:MOI.ScalarQuadraticFunction,S<:MOI.AbstractSet} + n = LinearAlgebra.checksquare(Q) + factor = LDLFactorizations.ldl(Q) + if minimum(factor.D) < 0 + 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 1031949636..9f3e167e95 100644 --- a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl +++ b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl @@ -60,6 +60,39 @@ end const QuadtoSOC{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{QuadtoSOCBridge{T},OT} +function compute_sparse_sqrt_root_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. + + In some cases 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_root(Q, func, set) + factor = LinearAlgebra.cholesky(Q; check = false) + if !LinearAlgebra.issuccess(factor) + return compute_sparse_sqrt_root_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 +end + function bridge_constraint( ::Type{QuadtoSOCBridge{T}}, model, @@ -73,21 +106,6 @@ function bridge_constraint( if !less_than LinearAlgebra.rmul!(Q, -1) end - F = try - LinearAlgebra.cholesky(LinearAlgebra.Symmetric(Q)) - catch - throw( - MOI.UnsupportedConstraint{typeof(func),typeof(set)}( - "Unable to transform a quadratic constraint into a " * - "second-order cone constraint because the quadratic " * - "constraint is not strongly convex.\n\nConvex constraints " * - "that are not strongly convex (that is, the matrix is positive " * - "semidefinite but not positive definite) are not supported " * - "yet.\n\nNote that a quadratic equality constraint is " * - "non-convex.", - ), - ) - end # Construct the VectorAffineFunction. We're aiming for: # | 1 | # | -a^T x + ub | ∈ RotatedSecondOrderCone() @@ -99,24 +117,20 @@ function bridge_constraint( MOI.ScalarAffineTerm(scale * term.coefficient, term.variable), ) for term in func.affine_terms ] - # Note that `L = U'` so we add the terms of L'x - L, p = SparseArrays.sparse(F.L), F.p - I, J, V = SparseArrays.findnz(L) - for i in 1:length(V) - # Cholesky is a pivoted decomposition, so L × L' == Q[p, p]. - # To get the variable, we need the row of L, I[i], then to map that - # through the permutation vector, so `p[I[i]]`, and then through the - # index_to_variable_map. - xi = index_to_variable_map[p[I[i]]] + I, J, V = compute_sparse_sqrt_root(LinearAlgebra.Symmetric(Q), func, set) + for (i, j, v) in zip(I, J, V) push!( vector_terms, - MOI.VectorAffineTerm(J[i] + 2, MOI.ScalarAffineTerm(V[i], xi)), + MOI.VectorAffineTerm( + i + 2, + MOI.ScalarAffineTerm(v, index_to_variable_map[j]), + ), ) end # This is the [1, ub, 0] vector... set_constant = MOI.constant(set) MOI.throw_if_scalar_and_constant_not_zero(func, typeof(set)) - vector_constant = vcat(one(T), -scale * set_constant, zeros(T, size(L, 1))) + vector_constant = vcat(one(T), -scale * set_constant, zeros(T, size(Q, 1))) f = MOI.VectorAffineFunction(vector_terms, vector_constant) dimension = MOI.output_dimension(f) soc = MOI.add_constraint(model, f, MOI.RotatedSecondOrderCone(dimension)) diff --git a/test/Bridges/Constraint/QuadtoSOCBridge.jl b/test/Bridges/Constraint/QuadtoSOCBridge.jl index 20db1a09df..323450f4e3 100644 --- a/test/Bridges/Constraint/QuadtoSOCBridge.jl +++ b/test/Bridges/Constraint/QuadtoSOCBridge.jl @@ -10,6 +10,7 @@ import LinearAlgebra import SparseArrays using Test +import LDLFactorizations import MathOptInterface as MOI function runtests() @@ -339,6 +340,53 @@ function test_constraint_primal_no_quad_terms() 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_compute_sparse_U_edge_cases() + f = zero(MOI.ScalarQuadraticFunction{Float64}) + s = MOI.GreaterThan(0.0) + for A in [ + # 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 eigen value + [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 + [2.0 0.0; 0.0 0.0], + ] + B = SparseArrays.sparse(A) + I, J, V = MOI.Bridges.Constraint.compute_sparse_sqrt_root(B, f, s) + U = zeros(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 + A = [-1.0 0.0; 0.0 1.0] + B = SparseArrays.sparse(A) + @test_throws( + MOI.UnsupportedConstraint{typeof(f),typeof(s)}, + MOI.Bridges.Constraint.compute_sparse_sqrt_root(B, f, s), + ) + return +end + end # module TestConstraintQuadToSOC.runtests() From d298200f24d3304678be25e5b4f95ed01da5edab Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Mon, 2 Feb 2026 10:11:02 +1300 Subject: [PATCH 2/3] Update --- src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl index 9f3e167e95..9e8b6bf68b 100644 --- a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl +++ b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl @@ -66,8 +66,8 @@ function compute_sparse_sqrt_root_fallback(Q, ::F, ::S) where {F,S} constraint because the quadratic constraint is not strongly convex and our Cholesky decomposition failed. - In some cases you can work-around this issue by manually installing and - loading `LDLFactorizations.jl`: + 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 From de7c66d7b15c0023faa65e98db7e377bebda614f Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Mon, 2 Feb 2026 13:39:48 +1300 Subject: [PATCH 3/3] Update --- ext/MathOptInterfaceLDLFactorizationsExt.jl | 2 +- .../Constraint/bridges/QuadtoSOCBridge.jl | 8 ++++---- test/Bridges/Constraint/QuadtoSOCBridge.jl | 20 ++++++++++++++++--- 3 files changed, 22 insertions(+), 8 deletions(-) diff --git a/ext/MathOptInterfaceLDLFactorizationsExt.jl b/ext/MathOptInterfaceLDLFactorizationsExt.jl index 7eeb8ad9d6..9211069776 100644 --- a/ext/MathOptInterfaceLDLFactorizationsExt.jl +++ b/ext/MathOptInterfaceLDLFactorizationsExt.jl @@ -14,7 +14,7 @@ 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_root_fallback( +function MOI.Bridges.Constraint.compute_sparse_sqrt_fallback( Q::AbstractMatrix, ::F, ::S, diff --git a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl index 9e8b6bf68b..944308ce9f 100644 --- a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl +++ b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl @@ -60,7 +60,7 @@ end const QuadtoSOC{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{QuadtoSOCBridge{T},OT} -function compute_sparse_sqrt_root_fallback(Q, ::F, ::S) where {F,S} +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 @@ -79,10 +79,10 @@ function compute_sparse_sqrt_root_fallback(Q, ::F, ::S) where {F,S} return throw(MOI.AddConstraintNotAllowed{F,S}(msg)) end -function compute_sparse_sqrt_root(Q, func, set) +function compute_sparse_sqrt(Q, func, set) factor = LinearAlgebra.cholesky(Q; check = false) if !LinearAlgebra.issuccess(factor) - return compute_sparse_sqrt_root_fallback(Q, func, set) + 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 @@ -117,7 +117,7 @@ function bridge_constraint( MOI.ScalarAffineTerm(scale * term.coefficient, term.variable), ) for term in func.affine_terms ] - I, J, V = compute_sparse_sqrt_root(LinearAlgebra.Symmetric(Q), func, set) + I, J, V = compute_sparse_sqrt(LinearAlgebra.Symmetric(Q), func, set) for (i, j, v) in zip(I, J, V) push!( vector_terms, diff --git a/test/Bridges/Constraint/QuadtoSOCBridge.jl b/test/Bridges/Constraint/QuadtoSOCBridge.jl index 323450f4e3..00b64a3ff1 100644 --- a/test/Bridges/Constraint/QuadtoSOCBridge.jl +++ b/test/Bridges/Constraint/QuadtoSOCBridge.jl @@ -355,7 +355,7 @@ function test_semidefinite_cholesky_fail() return end -function test_compute_sparse_U_edge_cases() +function test_compute_sparse_sqrt_edge_cases() f = zero(MOI.ScalarQuadraticFunction{Float64}) s = MOI.GreaterThan(0.0) for A in [ @@ -371,7 +371,7 @@ function test_compute_sparse_U_edge_cases() [2.0 0.0; 0.0 0.0], ] B = SparseArrays.sparse(A) - I, J, V = MOI.Bridges.Constraint.compute_sparse_sqrt_root(B, f, s) + I, J, V = MOI.Bridges.Constraint.compute_sparse_sqrt(B, f, s) U = zeros(size(A)) for (i, j, v) in zip(I, J, V) U[i, j] += v @@ -382,7 +382,21 @@ function test_compute_sparse_U_edge_cases() B = SparseArrays.sparse(A) @test_throws( MOI.UnsupportedConstraint{typeof(f),typeof(s)}, - MOI.Bridges.Constraint.compute_sparse_sqrt_root(B, f, s), + MOI.Bridges.Constraint.compute_sparse_sqrt(B, f, s), + ) + 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