diff --git a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl index 1031949636..2274cb854c 100644 --- a/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl +++ b/src/Bridges/Constraint/bridges/QuadtoSOCBridge.jl @@ -73,8 +73,20 @@ function bridge_constraint( if !less_than LinearAlgebra.rmul!(Q, -1) end - F = try - LinearAlgebra.cholesky(LinearAlgebra.Symmetric(Q)) + # Construct the VectorAffineFunction. We're aiming for: + # | 1 | + # | -a^T x + ub | ∈ RotatedSecondOrderCone() + # | Ux + 0 | + # Start with the -a^T x terms... + vector_terms = MOI.VectorAffineTerm{T}[ + MOI.VectorAffineTerm( + 2, + MOI.ScalarAffineTerm(scale * term.coefficient, term.variable), + ) for term in func.affine_terms + ] + # Compute a sparse U such that U' × U == Q. + I, J, V = try + _compute_sparse_U(Q) catch throw( MOI.UnsupportedConstraint{typeof(func),typeof(set)}( @@ -88,35 +100,19 @@ function bridge_constraint( ), ) end - # Construct the VectorAffineFunction. We're aiming for: - # | 1 | - # | -a^T x + ub | ∈ RotatedSecondOrderCone() - # | Ux + 0 | - # Start with the -a^T x terms... - vector_terms = MOI.VectorAffineTerm{T}[ - MOI.VectorAffineTerm( - 2, - 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]]] + 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)) @@ -129,6 +125,31 @@ function bridge_constraint( ) end +function _compute_sparse_U(Q::SparseArrays.SparseMatrixCSC) + return _compute_sparse_U(LinearAlgebra.Symmetric(Q)) +end + +function _compute_sparse_U(Q::LinearAlgebra.Symmetric) + F = LinearAlgebra.cholesky(Q; check = false) + if LinearAlgebra.issuccess(F) + L, p = SparseArrays.sparse(F.L), F.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 + # Cholesky failed. Use instead an eigen value decomposition. This is not + # triangular, but nothing says that it has to be. + E = LinearAlgebra.eigen(LinearAlgebra.Symmetric(Matrix(Q))) + if !all(isreal, E.values) || minimum(E.values) < 0 + error("Matrix is not PSD") + end + D = LinearAlgebra.Diagonal(sqrt.(E.values)) * E.vectors + return SparseArrays.findnz(SparseArrays.sparse(D)) +end + function _matrix_from_quadratic_terms( terms::Vector{MOI.ScalarQuadraticTerm{T}}, ) where {T} diff --git a/test/Bridges/Constraint/QuadtoSOCBridge.jl b/test/Bridges/Constraint/QuadtoSOCBridge.jl index 20db1a09df..0ee2feff31 100644 --- a/test/Bridges/Constraint/QuadtoSOCBridge.jl +++ b/test/Bridges/Constraint/QuadtoSOCBridge.jl @@ -339,6 +339,51 @@ 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, 0.0, sum_y])) + return +end + +function test_compute_sparse_U_edge_cases() + 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_U(B) + 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( + ErrorException("Matrix is not PSD"), + MOI.Bridges.Constraint._compute_sparse_U(B), + ) + return +end + end # module TestConstraintQuadToSOC.runtests()