Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix various bugs with Real-valued Hermitian matrices #3805

Merged
merged 7 commits into from
Aug 12, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 57 additions & 5 deletions src/operators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -333,11 +333,9 @@ end
function Base.:+(a::GenericAffExpr, q::GenericQuadExpr)
return GenericQuadExpr(a + q.aff, copy(q.terms))
end
function Base.:-(a::GenericAffExpr, q::GenericQuadExpr)
result = -q
# This makes an unnecessary copy of aff, but it's important for a to appear
# first.
result.aff = a + result.aff
function Base.:-(a::GenericAffExpr{S}, q::GenericQuadExpr{T}) where {S,T}
result = -_copy_convert_coef(_MA.promote_operation(-, S, T), q)
add_to_expression!(result.aff, a)
return result
end

Expand Down Expand Up @@ -571,3 +569,57 @@ function Base.complex(
)
return r + im * i
end

# These methods exist in LinearAlgebra for subtypes of Real. Without them, we
# return a `Matrix` which looses the Hermitian information.
function Base.:+(
Copy link
Member

@blegat blegat Aug 12, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add a comment that the fault returns a Matrix and we overload to return a Hermitian. This is already done in LinearAlgebra for subtypes of Real

odow marked this conversation as resolved.
Show resolved Hide resolved
A::LinearAlgebra.Symmetric{V,Matrix{V}},
B::LinearAlgebra.Hermitian,
) where {
V<:Union{
GenericVariableRef{<:Real},
GenericAffExpr{<:Real},
GenericQuadExpr{<:Real},
},
}
return LinearAlgebra.Hermitian(A) + B
end

function Base.:+(
A::LinearAlgebra.Hermitian,
B::LinearAlgebra.Symmetric{V,Matrix{V}},
) where {
V<:Union{
GenericVariableRef{<:Real},
GenericAffExpr{<:Real},
GenericQuadExpr{<:Real},
},
}
return A + LinearAlgebra.Hermitian(B)
end

function Base.:-(
A::LinearAlgebra.Symmetric{V,Matrix{V}},
B::LinearAlgebra.Hermitian,
) where {
V<:Union{
GenericVariableRef{<:Real},
GenericAffExpr{<:Real},
GenericQuadExpr{<:Real},
},
}
return LinearAlgebra.Hermitian(A) - B
end

function Base.:-(
A::LinearAlgebra.Hermitian,
B::LinearAlgebra.Symmetric{V,Matrix{V}},
) where {
V<:Union{
GenericVariableRef{<:Real},
GenericAffExpr{<:Real},
GenericQuadExpr{<:Real},
},
}
return A - LinearAlgebra.Hermitian(B)
end
16 changes: 16 additions & 0 deletions src/sd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -703,6 +703,22 @@ function build_constraint(
return VectorConstraint(x, MOI.Zeros(length(x)), shape)
end

# If we have a real-valued Hermitian matrix, then it is actually Symmetric, and
# not Complex-valued Hermitian.
function build_constraint(
error_fn::Function,
H::LinearAlgebra.Hermitian{V},
set::Zeros,
) where {
V<:Union{
GenericVariableRef{<:Real},
GenericAffExpr{<:Real},
GenericQuadExpr{<:Real},
},
}
return build_constraint(error_fn, LinearAlgebra.Symmetric(H), set)
end

reshape_set(s::MOI.Zeros, ::HermitianMatrixShape) = Zeros()

function build_constraint(error_fn::Function, ::AbstractMatrix, ::Nonnegatives)
Expand Down
19 changes: 19 additions & 0 deletions test/test_constraint.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2088,4 +2088,23 @@ function test_abstract_vector_orthants()
return
end

function test_real_hermitian_in_zeros()
model = Model()
@variable(model, x[1:2, 1:2], Symmetric)
c = @constraint(model, LinearAlgebra.Hermitian(x) in Zeros())
obj = constraint_object(c)
@test obj.func == [x[1, 1], x[1, 2], x[2, 2]]
@test obj.shape == SymmetricMatrixShape(2; needs_adjoint_dual = true)
H = LinearAlgebra.Hermitian([1 2; 2 3])
c = @constraint(model, x == H)
obj = constraint_object(c)
@test obj.func == [x[1, 1] - 1, x[1, 2] - 2, x[2, 2] - 3]
@test obj.shape == SymmetricMatrixShape(2; needs_adjoint_dual = true)
c = @constraint(model, LinearAlgebra.Hermitian(x .^ 2) in Zeros())
obj = constraint_object(c)
@test obj.func == [x[1, 1]^2, x[1, 2]^2, x[2, 2]^2]
@test obj.shape == SymmetricMatrixShape(2; needs_adjoint_dual = true)
return
end

end # module
38 changes: 36 additions & 2 deletions test/test_operator.jl
Original file line number Diff line number Diff line change
Expand Up @@ -469,7 +469,7 @@ function test_extension_basic_operators_affexpr(
@test_expression_with_string aff - aff "0 x"
# 4-4 AffExpr--QuadExpr
@test_expression_with_string aff2 + q "2.5 y*z + 1.2 y + 7.1 x + 3.7"
@test_expression_with_string aff2 - q "-2.5 y*z + 1.2 y - 7.1 x - 1.3"
@test_expression_with_string aff2 - q "-2.5 y*z - 7.1 x + 1.2 y - 1.3"
@test_expression_with_string aff2 * q "(1.2 y + 1.2) * (2.5 y*z + 7.1 x + 2.5)"
@test_expression_with_string aff2 / q "(1.2 y + 1.2) / (2.5 y*z + 7.1 x + 2.5)"
@test transpose(aff) === aff
Expand Down Expand Up @@ -499,7 +499,7 @@ function test_extension_basic_operators_quadexpr(
@test_expression_with_string q * 2 "5 y*z + 14.2 x + 5"
@test_expression_with_string q / 2 "1.25 y*z + 3.55 x + 1.25"
@test q == q
@test_expression_with_string aff2 - q "-2.5 y*z + 1.2 y - 7.1 x - 1.3"
@test_expression_with_string aff2 - q "-2.5 y*z - 7.1 x + 1.2 y - 1.3"
# 4-2 QuadExpr--Variable
@test_expression_with_string q + w "2.5 y*z + 7.1 x + w + 2.5"
@test_expression_with_string q - w "2.5 y*z + 7.1 x - w + 2.5"
Expand Down Expand Up @@ -689,4 +689,38 @@ function test_base_complex()
return
end

function test_aff_minus_quad()
model = Model()
@variable(model, x)
a, b = 1.0 * x, (2 + 3im) * x^2
@test a - b == -(b - a)
@test b - a == -(a - b)
a, b = (1.0 + 2im) * x, 3 * x^2 + 4 * x
@test a - b == -(b - a)
@test b - a == -(a - b)
return
end

function test_hermitian_and_symmetric()
model = Model()
@variable(model, A[1:2, 1:2], Symmetric)
@variable(model, B[1:2, 1:2], Hermitian)
for (x, y) in (
(A, B),
(B, A),
(1.0 * A, B),
(B, 1.0 * A),
(1.0 * A, 1.0 * B),
(1.0 * B, 1.0 * A),
(1.0 * LinearAlgebra.Symmetric(A .* A), 1.0 * B),
(1.0 * B, 1.0 * LinearAlgebra.Symmetric(A .* A)),
)
@test x + y isa LinearAlgebra.Hermitian
@test x + y == x .+ y
@test x - y isa LinearAlgebra.Hermitian
@test x - y == x .- y
end
return
end

end
Loading