From 871568d9018f6054017377d1f4f6c4726857bf35 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 11:12:51 +0200 Subject: [PATCH 01/10] Sampling basis --- .../src/tutorials/Getting started/sampling.jl | 23 ++++++++++++ src/Certificate/ideal.jl | 11 +++--- src/constraints.jl | 36 +++++++++++-------- 3 files changed, 52 insertions(+), 18 deletions(-) create mode 100644 docs/src/tutorials/Getting started/sampling.jl diff --git a/docs/src/tutorials/Getting started/sampling.jl b/docs/src/tutorials/Getting started/sampling.jl new file mode 100644 index 000000000..66e22ad55 --- /dev/null +++ b/docs/src/tutorials/Getting started/sampling.jl @@ -0,0 +1,23 @@ +# # Sampling basis + +#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Getting started/sampling.ipynb) +#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Getting started/sampling.ipynb) +# **Contributed by**: Benoît Legat + +using Test #src +using DynamicPolynomials +using SumOfSquares +import Hypatia + +# In this tutorial, we show how to use a different polynomial basis +# for enforcing the equality between the polynomial and its Sum-of-Squares decomposition. + +@polyvar x +p = x^4 - 4x^3 - 2x^2 + 12x + 9 + +model = Model(Hypatia.Optimizer) +set_silent(model) +@constraint(model, p in SOSCone(), zero_basis = BoxSampling([-1], [1])) +optimize!(model) +solution_summary(model) +backend(model) diff --git a/src/Certificate/ideal.jl b/src/Certificate/ideal.jl index 51f98c4fc..9f57c6a4d 100644 --- a/src/Certificate/ideal.jl +++ b/src/Certificate/ideal.jl @@ -60,15 +60,18 @@ function _reduce_with_domain(basis::MB.SubBasis{B}, domain) where {B} end function zero_basis( - ::AbstractIdealCertificate, + cert::AbstractIdealCertificate, basis, domain, gram_bases, weights, ) - return _reduce_with_domain( - _combine_with_gram(basis, gram_bases, weights), - domain, + return MB.explicit_basis_covering( + cert.zero_basis, + _reduce_with_domain( + _combine_with_gram(basis, gram_bases, weights), + domain, + ), ) end diff --git a/src/constraints.jl b/src/constraints.jl index 4b6c3aebd..c11290018 100644 --- a/src/constraints.jl +++ b/src/constraints.jl @@ -367,7 +367,7 @@ end _promote_coef_type(::Type{V}, ::Type) where {V<:JuMP.AbstractVariableRef} = V _promote_coef_type(::Type{F}, ::Type{T}) where {F,T} = promote_type(F, T) -function _default_basis( +function _default_gram_basis( coeffs, basis::MB.SubBasis{B}, gram_basis::MB.MonomialIndexedBasis{G}, @@ -375,7 +375,7 @@ function _default_basis( if B === G return coeffs, basis else - return _default_basis( + return _default_gram_basis( SA.SparseCoefficients(basis.monomials, coeffs), MB.implicit_basis(basis), gram_basis, @@ -383,20 +383,20 @@ function _default_basis( end end -function _default_basis( +function _default_gram_basis( p::SA.AbstractCoefficients, basis::MB.FullBasis{B}, gram_basis::MB.MonomialIndexedBasis{G}, ) where {B,G} if B === G - return _default_basis( + return _default_gram_basis( collect(SA.values(p)), MB.SubBasis{B}(collect(SA.keys(p))), gram_basis, ) else new_basis = MB.FullBasis{G,MP.monomial_type(typeof(basis))}() - return _default_basis( + return _default_gram_basis( SA.coeffs(p, basis, new_basis), new_basis, gram_basis, @@ -422,14 +422,14 @@ function _default_gram_basis(_, basis::MB.MonomialIndexedBasis) return basis end -function _default_basis(a::SA.AlgebraElement, basis) +function _default_gram_basis(a::SA.AlgebraElement, basis) b = SA.basis(a) gram = _default_gram_basis(b, basis) - return _default_basis(SA.coeffs(a), b, gram)..., gram + return _default_gram_basis(SA.coeffs(a), b, gram)..., gram end -function _default_basis(p::MP.AbstractPolynomialLike, basis) - return _default_basis( +function _default_gram_basis(p::MP.AbstractPolynomialLike, basis) + return _default_gram_basis( MB.algebra_element( MB.sparse_coefficients(p), MB.FullBasis{MB.Monomial,MP.monomial_type(p)}(), @@ -438,6 +438,11 @@ function _default_basis(p::MP.AbstractPolynomialLike, basis) ) end +_default_zero_basis(basis, ::Nothing) = MB.implicit_basis(basis) +function _default_zero_basis(_, nodes::MB.AbstractNodes) + return MB.ImplicitLagrangeBasis(nodes) +end + function JuMP.build_constraint( _error::Function, p, @@ -446,11 +451,14 @@ function JuMP.build_constraint( zero_basis = nothing, kws..., ) - __coefs, basis, gram_basis = _default_basis(p, basis) - if isnothing(zero_basis) - zero_basis = MB.implicit_basis(basis) - end - set = JuMP.moi_set(cone, basis, gram_basis, zero_basis; kws...) + __coefs, basis, gram_basis = _default_gram_basis(p, basis) + set = JuMP.moi_set( + cone, + basis, + gram_basis, + _default_zero_basis(basis, zero_basis); + kws..., + ) _coefs = PolyJuMP.non_constant(__coefs) # If a polynomial with real coefficients is used with the Hermitian SOS # cone, we want to promote the coefficients to complex From 242646cd0cb91768ad6457e76b2a4992c3760ffe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 14:21:53 +0200 Subject: [PATCH 02/10] Fixes --- src/Certificate/ideal.jl | 16 ++++++++++++---- src/constraints.jl | 4 ++-- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/src/Certificate/ideal.jl b/src/Certificate/ideal.jl index 9f57c6a4d..64ffd667c 100644 --- a/src/Certificate/ideal.jl +++ b/src/Certificate/ideal.jl @@ -273,15 +273,23 @@ function _quotient_basis_type( } end +function _zero_basis_type(::Type{<:SimpleIdealCertificate{C,G,Z}}) where {C,G,Z} + return Z +end + +function _zero_basis_type(::Type{Remainder{C}}) where {C} + return _zero_basis_type(C) +end + function MA.promote_operation( ::typeof(zero_basis), - ::Type{<:Union{SimpleIdealCertificate,Remainder}}, - ::Type{B}, + C::Type{<:Union{SimpleIdealCertificate,Remainder}}, + ::Type, ::Type{SemialgebraicSets.FullSpace}, ::Type, ::Type, -) where {B} - return B +) + return MB.explicit_basis_type(_zero_basis_type(C)) end function MA.promote_operation( diff --git a/src/constraints.jl b/src/constraints.jl index c11290018..c80ff0f17 100644 --- a/src/constraints.jl +++ b/src/constraints.jl @@ -439,8 +439,8 @@ function _default_gram_basis(p::MP.AbstractPolynomialLike, basis) end _default_zero_basis(basis, ::Nothing) = MB.implicit_basis(basis) -function _default_zero_basis(_, nodes::MB.AbstractNodes) - return MB.ImplicitLagrangeBasis(nodes) +function _default_zero_basis(basis, nodes::MB.AbstractNodes) + return MB.ImplicitLagrangeBasis(MP.variables(basis), nodes) end function JuMP.build_constraint( From d31d85fc85e41ca363964065fe3e7083c8b7830f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 16:40:57 +0200 Subject: [PATCH 03/10] Fixes --- src/Certificate/ideal.jl | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/Certificate/ideal.jl b/src/Certificate/ideal.jl index 64ffd667c..fd3999b10 100644 --- a/src/Certificate/ideal.jl +++ b/src/Certificate/ideal.jl @@ -40,9 +40,11 @@ function _combine_with_gram( return MB.SubBasis{B}(keys(SA.coeffs(p))) end -_reduce_with_domain(basis::MB.SubBasis, ::FullSpace) = basis +function _reduce_with_domain(basis::MB.SubBasis, zero_basis, ::FullSpace) + return MB.explicit_basis_covering(zero_basis, basis) +end -function _reduce_with_domain(basis::MB.SubBasis{B}, domain) where {B} +function _reduce_with_domain(basis::MB.SubBasis{B}, zero_basis::MB.FullBasis{B}, domain) where {B} if B !== MB.Monomial error("Only Monomial basis support with an equalities in domain") end @@ -59,6 +61,10 @@ function _reduce_with_domain(basis::MB.SubBasis{B}, domain) where {B} ) end +_zero_basis(c::SimpleIdealCertificate) = c.zero_basis + +_zero_basis(c::Remainder) = _zero_basis(c.gram_certificate) + function zero_basis( cert::AbstractIdealCertificate, basis, @@ -66,12 +72,10 @@ function zero_basis( gram_bases, weights, ) - return MB.explicit_basis_covering( - cert.zero_basis, - _reduce_with_domain( - _combine_with_gram(basis, gram_bases, weights), - domain, - ), + return _reduce_with_domain( + _combine_with_gram(basis, gram_bases, weights), + _zero_basis(cert), + domain, ) end From 9ff438c4cdb83b83f4b856ffc2306c5a5687f9f5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 16:42:28 +0200 Subject: [PATCH 04/10] Fix format --- src/Certificate/ideal.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/Certificate/ideal.jl b/src/Certificate/ideal.jl index fd3999b10..c8afac692 100644 --- a/src/Certificate/ideal.jl +++ b/src/Certificate/ideal.jl @@ -44,7 +44,11 @@ function _reduce_with_domain(basis::MB.SubBasis, zero_basis, ::FullSpace) return MB.explicit_basis_covering(zero_basis, basis) end -function _reduce_with_domain(basis::MB.SubBasis{B}, zero_basis::MB.FullBasis{B}, domain) where {B} +function _reduce_with_domain( + basis::MB.SubBasis{B}, + ::MB.FullBasis{B}, + domain, +) where {B} if B !== MB.Monomial error("Only Monomial basis support with an equalities in domain") end From c526218c5867f870cb3493167cbefd4b757c3f79 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 17:05:37 +0200 Subject: [PATCH 05/10] Fix --- src/Certificate/ideal.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Certificate/ideal.jl b/src/Certificate/ideal.jl index c8afac692..dad1e7e32 100644 --- a/src/Certificate/ideal.jl +++ b/src/Certificate/ideal.jl @@ -65,10 +65,6 @@ function _reduce_with_domain( ) end -_zero_basis(c::SimpleIdealCertificate) = c.zero_basis - -_zero_basis(c::Remainder) = _zero_basis(c.gram_certificate) - function zero_basis( cert::AbstractIdealCertificate, basis, @@ -281,10 +277,14 @@ function _quotient_basis_type( } end +_zero_basis(c::SimpleIdealCertificate) = c.zero_basis + function _zero_basis_type(::Type{<:SimpleIdealCertificate{C,G,Z}}) where {C,G,Z} return Z end +_zero_basis(c::Remainder) = _zero_basis(c.gram_certificate) + function _zero_basis_type(::Type{Remainder{C}}) where {C} return _zero_basis_type(C) end From cbcabf0ac2caf051919425d5c4eb3fc53b9e96d5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 17:55:00 +0200 Subject: [PATCH 06/10] Fix --- src/Certificate/ideal.jl | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/Certificate/ideal.jl b/src/Certificate/ideal.jl index dad1e7e32..8c2d07995 100644 --- a/src/Certificate/ideal.jl +++ b/src/Certificate/ideal.jl @@ -44,14 +44,18 @@ function _reduce_with_domain(basis::MB.SubBasis, zero_basis, ::FullSpace) return MB.explicit_basis_covering(zero_basis, basis) end -function _reduce_with_domain( - basis::MB.SubBasis{B}, - ::MB.FullBasis{B}, +function _reduce_with_domain(basis, zero_basis, domain) + return __reduce_with_domain(basis, zero_basis, domain) +end + +function __reduce_with_domain(_, _, _) + error("Only Monomial basis support with an equalities in domain") +end +function __reduce_with_domain( + basis::MB.SubBasis{MB.Monomial}, + ::MB.FullBasis{MB.Monomial}, domain, -) where {B} - if B !== MB.Monomial - error("Only Monomial basis support with an equalities in domain") - end +) I = ideal(domain) # set of standard monomials that are hit standard = Set{eltype(basis.monomials)}() From 9ad01fcb209d2a9dd3b7b96880196dcc1d12d802 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 17:55:21 +0200 Subject: [PATCH 07/10] Fix format --- src/Certificate/ideal.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Certificate/ideal.jl b/src/Certificate/ideal.jl index 8c2d07995..4dcb8c32f 100644 --- a/src/Certificate/ideal.jl +++ b/src/Certificate/ideal.jl @@ -49,7 +49,7 @@ function _reduce_with_domain(basis, zero_basis, domain) end function __reduce_with_domain(_, _, _) - error("Only Monomial basis support with an equalities in domain") + return error("Only Monomial basis support with an equalities in domain") end function __reduce_with_domain( basis::MB.SubBasis{MB.Monomial}, From 70027f127da93dad25663243d10f4d82366da103 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 17:58:41 +0200 Subject: [PATCH 08/10] Fix --- src/constraints.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/constraints.jl b/src/constraints.jl index c80ff0f17..866a29775 100644 --- a/src/constraints.jl +++ b/src/constraints.jl @@ -367,7 +367,7 @@ end _promote_coef_type(::Type{V}, ::Type) where {V<:JuMP.AbstractVariableRef} = V _promote_coef_type(::Type{F}, ::Type{T}) where {F,T} = promote_type(F, T) -function _default_gram_basis( +function _default_basis( coeffs, basis::MB.SubBasis{B}, gram_basis::MB.MonomialIndexedBasis{G}, @@ -375,7 +375,7 @@ function _default_gram_basis( if B === G return coeffs, basis else - return _default_gram_basis( + return _default_basis( SA.SparseCoefficients(basis.monomials, coeffs), MB.implicit_basis(basis), gram_basis, @@ -383,20 +383,20 @@ function _default_gram_basis( end end -function _default_gram_basis( +function _default_basis( p::SA.AbstractCoefficients, basis::MB.FullBasis{B}, gram_basis::MB.MonomialIndexedBasis{G}, ) where {B,G} if B === G - return _default_gram_basis( + return _default_basis( collect(SA.values(p)), MB.SubBasis{B}(collect(SA.keys(p))), gram_basis, ) else new_basis = MB.FullBasis{G,MP.monomial_type(typeof(basis))}() - return _default_gram_basis( + return _default_basis( SA.coeffs(p, basis, new_basis), new_basis, gram_basis, @@ -422,14 +422,14 @@ function _default_gram_basis(_, basis::MB.MonomialIndexedBasis) return basis end -function _default_gram_basis(a::SA.AlgebraElement, basis) +function _default_basis(a::SA.AlgebraElement, basis) b = SA.basis(a) gram = _default_gram_basis(b, basis) - return _default_gram_basis(SA.coeffs(a), b, gram)..., gram + return _default_basis(SA.coeffs(a), b, gram)..., gram end -function _default_gram_basis(p::MP.AbstractPolynomialLike, basis) - return _default_gram_basis( +function _default_basis(p::MP.AbstractPolynomialLike, basis) + return _default_basis( MB.algebra_element( MB.sparse_coefficients(p), MB.FullBasis{MB.Monomial,MP.monomial_type(p)}(), From 2f2c32ea593e7213ea05841ba96073a6f8bc3b2e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 18:13:29 +0200 Subject: [PATCH 09/10] Update src/constraints.jl --- src/constraints.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/constraints.jl b/src/constraints.jl index 866a29775..9f6b096f8 100644 --- a/src/constraints.jl +++ b/src/constraints.jl @@ -451,7 +451,7 @@ function JuMP.build_constraint( zero_basis = nothing, kws..., ) - __coefs, basis, gram_basis = _default_gram_basis(p, basis) + __coefs, basis, gram_basis = _default_basis(p, basis) set = JuMP.moi_set( cone, basis, From 39dc9ebe5f3b8a7e924a724cdbd3392e86fd04b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 22:26:48 +0200 Subject: [PATCH 10/10] Delete docs/src/tutorials/Getting started/sampling.jl --- .../src/tutorials/Getting started/sampling.jl | 23 ------------------- 1 file changed, 23 deletions(-) delete mode 100644 docs/src/tutorials/Getting started/sampling.jl diff --git a/docs/src/tutorials/Getting started/sampling.jl b/docs/src/tutorials/Getting started/sampling.jl deleted file mode 100644 index 66e22ad55..000000000 --- a/docs/src/tutorials/Getting started/sampling.jl +++ /dev/null @@ -1,23 +0,0 @@ -# # Sampling basis - -#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Getting started/sampling.ipynb) -#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Getting started/sampling.ipynb) -# **Contributed by**: Benoît Legat - -using Test #src -using DynamicPolynomials -using SumOfSquares -import Hypatia - -# In this tutorial, we show how to use a different polynomial basis -# for enforcing the equality between the polynomial and its Sum-of-Squares decomposition. - -@polyvar x -p = x^4 - 4x^3 - 2x^2 + 12x + 9 - -model = Model(Hypatia.Optimizer) -set_silent(model) -@constraint(model, p in SOSCone(), zero_basis = BoxSampling([-1], [1])) -optimize!(model) -solution_summary(model) -backend(model)