From ba36ab4b7ae4103c1528c7150ac0dfd22c5791df Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 28 Oct 2023 23:31:51 +0200 Subject: [PATCH] Projection in minimizer extraction tuto --- .../extracting_minimizers.jl | 131 +++++++++++++++++- 1 file changed, 126 insertions(+), 5 deletions(-) diff --git a/docs/src/tutorials/Polynomial Optimization/extracting_minimizers.jl b/docs/src/tutorials/Polynomial Optimization/extracting_minimizers.jl index 8d795bb38..571b04e60 100644 --- a/docs/src/tutorials/Polynomial Optimization/extracting_minimizers.jl +++ b/docs/src/tutorials/Polynomial Optimization/extracting_minimizers.jl @@ -7,6 +7,7 @@ # [L09] Laurent, Monique. # *Sums of squares, moment matrices and optimization over polynomials.* # Emerging applications of algebraic geometry (2009): 157-270. +# World Scientific, **2009**. # ## Introduction @@ -16,20 +17,95 @@ # $x_2 \le 2x_1^4 - 8x_1^3 + 8x_1^2 + 2$, # $x_2 \le 4x_1^4 - 32x_1^3 + 88x_1^2 - 96x_1 + 36$ and the box constraints # $0 \le x_1 \le 3$ and $0 \le x_2 \le 4$, -# World Scientific, **2009**. using Test #src using DynamicPolynomials -@polyvar x[1:2] +DynamicPolynomials.@polyvar x[1:2] p = -sum(x) using SumOfSquares K = @set x[1] >= 0 && x[1] <= 3 && x[2] >= 0 && x[2] <= 4 && x[2] <= 2x[1]^4 - 8x[1]^3 + 8x[1]^2 + 2 && x[2] <= 4x[1]^4 - 32x[1]^3 + 88x[1]^2 - 96x[1] + 36 +f1 = 2x[1]^4 - 8x[1]^3 + 8x[1]^2 + 2 +f2 = 4x[1]^4 - 32x[1]^3 + 88x[1]^2 - 96x[1] + 36 + +import MultivariatePolynomials as MP + +function nonlinear_vars!(vars, p) + for mono in MP.monomials(p) + if MP.degree(mono) > 1 + for var in MP.effective_variables(mono) + push!(vars, var) + end + end + end +end + +function nonlinear_vars(K::BasicSemialgebraicSet{T,P}) where {T,P} + vars = Set{MP.variable_union_type(P)}() + for p in inequalities(K) + nonlinear_vars!(vars, p) + end + for p in equalities(K) + nonlinear_vars!(vars, p) + end + return sort(collect(vars), rev=true) +end + +function polyhedral_proj!(vars, v, K) + nl_vars = nonlinear_vars(K) + lin_vars = setdiff(vars, nl_vars) +end + +function proj!(v) + v[1] = min(max(v[1], 0), 3) + v[2] = min(v[2], f1(v[1])) + v[2] = min(v[2], f2(v[1])) + v[2] = min(max(v[2], 0), 4) + return v +end + +proj!([3, 4]) + +function laurent(ν) + μ = measure(ν) + v = [moment_value(μ, var) for var in x] + @show v + proj!(v) + return v, p(x => v) +end + +import MultivariatePolynomials as MP +function gaussian(ν, leading::Bool = true) + monos = MP.monomials(MP.variables(ν.basis.monomials), 0:1) + I = [MultivariateMoments._index(ν.basis, mono) for mono in monos] + Q = ν.Q[I, I] + F = eigen(Q) + display(F) + best_v = nothing + best_obj = Inf + J = leading ? size(F.vectors, 2) : axes(F.vectors, 2) + for i in J + v = F.vectors[:, i] + @show v + v /= v[1] + v = v[2:end] + @show v + proj!(v) + @show v + obj = p(x => v) + @show obj + if obj < best_obj + best_v = v + best_obj = obj + end + end + return best_v, best_obj +end # We will now see how to find the optimal solution using Sum of Squares Programming. # We first need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v1.12/installation/#Supported-solvers) for a list of the available choices. -import CSDP -solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true) +import Clarabel +solver = Clarabel.Optimizer # A Sum-of-Squares certificate that $p \ge \alpha$ over the domain `S`, ensures that $\alpha$ is a lower bound to the polynomial optimization problem. # The following function searches for the largest lower bound and finds zero using the `d`th level of the hierarchy`. @@ -49,12 +125,31 @@ end model4 = solve(4) @test objective_value(model4) ≈ -7 rtol=1e-4 #src @test termination_status(model4) == MOI.OPTIMAL #src +ν4 = moment_matrix(model4[:c]) +laurent(ν4) + +import MultivariateMoments as MM +function _vec(Q) + n = size(Q, 2) + return [Q[i, j] for j in 1:n for i in 1:j] +end +function truncate(ν, d) + vars = MP.variables(ν.basis) + monos = MP.monomials(vars, 0:d) + I = [MM._index(ν.basis, mono) for mono in monos] + return MM.MomentMatrix( + MM.SymMatrix(_vec(ν.Q[I, I]), length(I)), + MonomialBasis(monos), + ) +end # The second level improves the lower bound model5 = solve(5) @test objective_value(model5) ≈ -20/3 rtol=1e-4 #src @test termination_status(model5) == MOI.OPTIMAL #src +ν5 = moment_matrix(model5[:c]) +laurent(ν5) # The third level finds the optimal objective value as lower bound... @@ -65,12 +160,38 @@ model7 = solve(7) # ...and proves it by exhibiting the minimizer. ν7 = moment_matrix(model7[:c]) -η = atomic_measure(ν7, 1e-3) # Returns nothing as the dual is not atomic +laurent(ν7) +η = atomic_measure(ν7, 1e-3) @test length(η.atoms) == 1 #src @test η.atoms[1].center ≈ [2.3295, 3.1785] rtol=1e-4 #src # We can indeed verify that the objective value at `x_opt` is equal to the lower bound. +opt = [2.3295, 3.1785] x_opt = η.atoms[1].center @test x_opt ≈ [2.3295, 3.1785] rtol=1e-4 #src p(x_opt) +#x_opt = + +#atomic_measure(ν4, UserRank()) + +compute_support!(ν5, FixedRank(1)) +ν5.support + +ν5_2 = truncate(ν5, 1) +compute_support!(ν5_2, FixedRank(2)) +ν5_2.support + +[p(x => x_opt) for p in ν5.support.I.p] +[p(x => [2.6666, 1.23456]) for p in ν5.support.I.p] +compute_support!(ν5, FixedRank(2)) +ν5.support +pp = -6.945530612461534 - 2.2636172932906113*x[2] + x[2]^2 + +lag, system = PolyJuMP.lagrangian_kkt(MOI.MIN_SENSE, 1.0 * p, K) + +using Macaulay +νmax3 = moment_matrix(system.I.p, Clarabel.Optimizer, 3) +laurent(νmax3) +νmax4 = moment_matrix(system.I.p, Clarabel.Optimizer, 4) +laurent(νmax4)