Что можно¶

Современные алгоритмы позволяют достаточно эффективно решать следующие задачи:

Проверить является ли данный многочлен SOS
Для многочлена $p(x_1, \ldots, x_n, c_1, \ldots, c_m)$ линейного относительно $c_1, \ldots, c_m$ и линейной функции $\ell(c_1, \ldots, c_m)$ найти максимум $\ell(c_1, \ldots, c_m)$ среди всех наборов чисел $c_1, \ldots, c_m$ таких, что $p(x_1, \ldots, x_n, c_1, \ldots, c_m)$ является SOS относительно $c_1, \ldots, c_m$ .
Все это ещё можно делать не на всем пространстве, а при наличии условий на $x_1, \ldots, x_n$ (про это завтра).

Одной из библиотек, которая умеет делать такие вещи является SumOfSquares.jl на языке Julia.

using SumOfSquares
using DynamicPolynomials
using Mosek
using MosekTools

┌ Info: Precompiling SumOfSquares [4b9e565b-77fc-50a5-a571-1244f986bda1]
└ @ Base loading.jl:1260
┌ Info: Precompiling MosekTools [1ec41992-ff65-5c91-ac43-2df89e9693a4]
└ @ Base loading.jl:1260

Проверка SOS¶

Пример 1. Неравенство о средних для двух чисел $x^2 + y^2 - 2 xy \geqslant 0$

@polyvar x y
opt = optimizer_with_attributes(Mosek.Optimizer, "QUIET" => true);
model = SOSModel(opt)
con_ref = @constraint(model, x^2 - 2 * x * y + y^2 >= 0)
optimize!(model)
@show termination_status(model)
round(sos_decomposition(con_ref)[1], digits=3)

termination_status(model) = MathOptInterface.OPTIMAL

Пример 2. Неравенство о средних для трех чисел $x^6 + y^6 + z^6 - 3 x^2 y^2 z^2 \geqslant 0$

@polyvar x y z
opt = optimizer_with_attributes(Mosek.Optimizer, "QUIET" => true);
model = SOSModel(opt)
con_ref = @constraint(model, x^6 + y^6 + z^6 - 3 * x^2 * y^2 * z^2 >= 0)
optimize!(model)
@show termination_status(model)
for p = sos_decomposition(con_ref)
    print(round(p, digits=3), "\n")
end

termination_status(model) = MathOptInterface.OPTIMAL
0.056*x^3 + 0.031*x^2*y + 0.032*x^2*z - 0.892*x*y^2 + 0.836*x*z^2 - 0.007*y^3 - 0.022*y^2*z - 0.024*y*z^2 - 0.01*z^3
-0.936*x^2*y - 0.042*x^2*z - 0.029*x*y^2 + 0.03*x*z^2 + 0.167*y^3 + 0.051*y^2*z + 0.768*y*z^2 - 0.008*z^3
0.014*x^3 + 0.057*x^2*y - 0.923*x^2*z - 0.033*x*y^2 + 0.019*x*z^2 - 0.029*y^3 + 0.791*y^2*z - 0.028*y*z^2 + 0.132*z^3
0.706*x^3 - 0.112*x^2*y + 0.237*x^2*z - 0.316*x*y^2 - 0.391*x*z^2 + 0.338*y^3 + 0.375*y^2*z - 0.226*y*z^2 - 0.612*z^3
-0.705*x^3 - 0.1*x^2*y + 0.227*x^2*z + 0.321*x*y^2 + 0.384*x*z^2 + 0.314*y^3 + 0.397*y^2*z - 0.214*y*z^2 - 0.624*z^3
0.019*x^3 + 0.312*x^2*y + 0.196*x^2*z - 0.004*x*y^2 - 0.015*x*z^2 - 0.871*y^3 + 0.272*y^2*z + 0.558*y*z^2 - 0.468*z^3
0.0
0.0
0.0
0.0

Заметим, что разложение отличается от того, которое было в анонсе курса, их, вообще говоря, бесконечно много разных.

Пример 3. Многочлен Моцкина $x^4 y^2 + x^2 y^4 - 3 x^2 y^2 + 1$

opt = optimizer_with_attributes(Mosek.Optimizer, "QUIET" => true);
model = SOSModel(opt)
moz = x^4*y^2 + x^2*y^4 - 3*x^2*y^2 + 1
@constraint(model, moz >= 0)
optimize!(model)
@show termination_status(model)

termination_status(model) = MathOptInterface.INFEASIBLE

INFEASIBLE::TerminationStatusCode = 2

Это означает, что разложения в сумму квадратов нет.

Пример 4. А теперь умножим на $(x^2 + y^2)$

opt = optimizer_with_attributes(Mosek.Optimizer, "QUIET" => true);
model = SOSModel(opt)
moz = x^4*y^2 + x^2*y^4 - 3*x^2*y^2 + 1
con_ref = @constraint(model, moz * (x^2 + y^2) >= 0)
optimize!(model)
@show termination_status(model)
for p = sos_decomposition(con_ref)
    print(round(p, digits=3), "\n")
end

Утверждается, что многочлен Моцкина умноженный на $x^2 + y^2$ является суммой квадратов многочленов выше. Заметим, что это приближенное равенство (и нулевые многочлены на самом деле не нулевые, а очень близкие к нулю). Однако, мы можем округлить числа и проверить, вдруг получится точное разложение:

((-x^3*y - x*y^3 + 2*x*y)^2 + (-x^2*y + y)^2 + (x*y^2 - x)^2)  - moz * (x^2 + y^2)

И действительно получилось! То, как округлять, чтобы получать верные равенства - отдельная история.

SOS оптимизация¶

Простешим примером использования SOS оптимизации является поиск нижних оценко на значение многочлена. Например, если мы хотим найти нижнюю оценку на $p(x_1, \ldots, x_n)$ , мы можем поискать максимум $c$ такой, что $p(x_1, \ldots, x_n) - c$ является SOS

opt = optimizer_with_attributes(Mosek.Optimizer, "QUIET" => true);
model = SOSModel(opt)
@variable(model, c)
moz = x^4*y^2 + x^2*y^4 - 3*x^2*y^2 + 1
@constraint(model, moz - c >= 0)
@objective(model, Max, c)
optimize!(model)
@show termination_status(model)
@show round(value(c), digits = 3)

termination_status(model) = MathOptInterface.INFEASIBLE
round(value(c), digits = 3) = -0.652

-0.652

Так мы выяснили, что многочлен Моцкина всегда $\geqslant -0.655$ . Но ведь мы же знаем, что он неотрицателен! Чтобы получить оценку получше, можно прибегнуть к теореме Артина (а точнее Пойа-Резника)

opt = optimizer_with_attributes(Mosek.Optimizer, "QUIET" => true);
model = SOSModel(opt)
@variable(model, c)
moz = x^4 * y^2 + x^2 * y^4 - 3 * x^2 * y^2 + 1
@constraint(model, moz * (x^2 + y^2) - c * (x^2 + y^2) >= 0)
@objective(model, Max, c)
optimize!(model)
@show termination_status(model)
@show round(value(c), digits = 3)

termination_status(model) = MathOptInterface.OPTIMAL
round(value(c), digits = 3) = -0.0

-0.0