Аппроксимация эллипсоидом
artpgchВведение
Представьте, что у вас есть набор точек, разбросанных в пространстве, и вы хотите найти идеальную эллиптическую границу, которая охватит их все, имея при этом наименьший возможный объем. Эта классическая задача оптимизации — поиск эллипсоида минимального объёма — находит применение в машинном обучении, робототехнике и анализе данных.
В этом примере мы не только решим эту задачу, но и детально исследуем внутренние механизмы процесса оптимизации. Используя функции пакета JuMP в Julia, мы покажем, как современные инструменты автоматически преобразуют сложные математические формулировки в формы, понятные для решателей. Мы изучим, какие преобразования используются, и исследуем, как разные формулировки одной и той же задачи влияют на скорость и эффективность решения.
using JuMP, Clarabel, LinearAlgebra, Plots, Random, Test
Постановка задачи
Предположим, что нам дано множество , состоящее из точек в -мерном пространстве:
Наша цель — найти оптимальный вектор и оптимальную симметричную матрицу , такие, чтобы эллипс:
,
содержал множество S и имел наименее возможный объём.
Оптимальные значения и задаются следующей задачей оптимизации:
где .
и это оптимальные значения и , полученные в результате решения задачи оптимизации.
Данные
Сгенерируем несколько точек:
function создать_набор_точек(
m; # количество 2-мерных точек
a = 10, # масштаб по оси x
b = 2, # масштаб по оси y
ρ = π / 6, # вращение точек вокруг начала координат
random_seed = 1,
)
rng = Random.MersenneTwister(random_seed)
P = randn(rng, Float64, m, 2)
Ф = [a*cos(ρ) a*sin(ρ); -b*sin(ρ) b*cos(ρ)]
S = P * Ф
return S
end
Для простоты этого примера возьмём :
S = создать_набор_точек(100);
Визуализируем точки:
r = 1.1 * maximum(abs.(S))
plot = Plots.scatter(S[:, 1], S[:, 2]; xlim = (-r, r), ylim = (-r, r), label = nothing, c = :green, shape = :x, size = (600, 600))
Модель JuMP
Теперь давайте построим модель JuMP. В результате решения задачи будет вычислены значения и .
модель = Model(Clarabel.Optimizer)
set_silent(модель)
m, n = size(S)
@variable(модель, z[1:n])
@variable(модель, Z[1:n, 1:n], PSD)
@variable(модель, s)
@variable(модель, t)
@constraint(модель, [s z'; z Z] >= 0, PSDCone())
@constraint(
модель,
[i in 1:m],
S[i, :]' * Z * S[i, :] - 2 * S[i, :]' * z + s <= 1,
)
@constraint(модель, [t; vec(Z)] in MOI.RootDetConeSquare(n))
@objective(модель, Max, t)
optimize!(модель)
assert_is_solved_and_feasible(модель)
solution_summary(модель)
Результаты вычислений
В результате оптимизации модели, получаем значения и :
D = value.(Z)
c = D \ value.(z)
Мы можем убедиться, что каждая точка лежит внутри эллипсоида, проверив, что величина наибольшего нормализованного радиуса не превышает единицу:
наибольший_радиус = maximum(map(x -> (x - c)' * D * (x - c), eachrow(S)))
Совмещая полученное решение с графиком, мы видим эллипсоид, минимизирующий объём аппроксимации:
P = sqrt(D)
q = -P * c
data = [tuple(P \ [cos(θ) - q[1], sin(θ) - q[2]]...) for θ in 0:0.05:(2pi+0.05)]
Plots.plot!(plot, data; c = :crimson, label = nothing)
Альтернативные постановки задачи
Модель использует специализированный конус MOI.RootDetConeSquare, который не поддерживается решателем Clarabel. Для обеспечения совместимости JuMP автоматически преобразует задачу в эквивалентную формулировку, используя доступные решателю конструкции. Анализ применённых преобразований доступен через команду print_active_bridges, отображающую цепочку математических транcформаций.
print_active_bridges(модель)
Анализ вывода показывает, что решатель Clarabel не поддерживает целевые функции типа MOI.VariableIndex. Для обеспечения совместимости JuMP автоматически применил преобразование FunctionConversionBridge, трансформировав целевую функцию в поддерживаемый формат MOI.ScalarAffineFunction{Float64}.
Для оптимизации взаимодействия с решателем можно либо использовать автоматическую переформулировку JuMP, либо явно задать целевую функцию в изначально поддерживаемом формате.
@objective(модель, Max, 1.0 * t + 0.0);
Активные мосты (active bridges) — это концепция в контексте JuMP и MOI (MathOptInterface), которая используется для описания механизмов преобразования или адаптации между различными типами данных или структур в процессе оптимизации.
С помощью команды print_active_bridges, снова выведем активные мосты:
print_active_bridges(модель)
Для оптимизации модели под решатель Clarabel можно выполнить следующие эквивалентные преобразования:
-
Переформулировать ограничения положительной определённости, заменив прямую декларацию переменных как положительно определённых матриц на явные ограничения, которые решатель понимает лучше.
-
Преобразовать ограничение с блочной матрицей, используя её симметричную форму.
-
Объединить скалярные неравенства в одно векторное ограничение на неотрицательность.
-
Изменить тип ограничения на определитель, используя его альтернативное представление, с которым решатель работает эффективнее.
Эти изменения позволяют сократить количество автоматических преобразований в модели и повысить эффективность её решения.
модель = Model(Clarabel.Optimizer)
set_silent(модель)
@variable(модель, z[1:n])
@variable(модель, s)
@variable(модель, t)
@variable(модель, Z[1:n, 1:n], Symmetric)
@constraint(модель, Z >= 0, PSDCone())
@constraint(модель, LinearAlgebra.Symmetric([s z'; z Z]) >= 0, PSDCone())
f = [1 - S[i, :]' * Z * S[i, :] + 2 * S[i, :]' * z - s for i in 1:m]
@constraint(модель, f in MOI.Nonnegatives(m))
@constraint(модель, 1 * [t; triangle_vec(Z)] .+ 0 in MOI.RootDetConeTriangle(n))
@objective(модель, Max, 1 * t + 0)
optimize!(модель)
assert_is_solved_and_feasible(модель)
solve_time_1 = solve_time(модель)
Данная формулировка выводит наиболее краткий результат:
print_active_bridges(модель)
Анализ текущей модели показывает, что ограничения типа MOI.PositiveSemidefiniteConeTriangle остаются неподдерживаемыми решателем Clarabel. Для достижения полной совместимости требуется их дополнительное преобразование в эквивалентные конструкции, изначально распознаваемые решателем. Это позволит минимизировать цепочку автоматических преобразований и оптимизировать процесс решения.
модель = Model(Clarabel.Optimizer)
set_silent(модель)
@variable(модель, z[1:n])
@variable(модель, s)
@variable(модель, t)
@variable(модель, Z[1:n, 1:n], Symmetric)
f = triangle_vec(Z)
scale_f = [1.0, sqrt(2), 1.0]
@constraint(модель, scale_f .* f in MOI.Scaled(MOI.PositiveSemidefiniteConeTriangle(n)))
g = triangle_vec(LinearAlgebra.Symmetric([s z'; z Z]))
scale_g = [1.0, sqrt(2), 1.0, sqrt(2), sqrt(2), 1.0]
@constraint(модель, scale_g .* g in MOI.Scaled(MOI.PositiveSemidefiniteConeTriangle(1 + n)))
f = [1 - S[i, :]' * Z * S[i, :] + 2 * S[i, :]' * z - s for i in 1:m]
@constraint(модель, f in MOI.Nonnegatives(m))
@constraint(модель, 1 * [t; triangle_vec(Z)] .+ 0 in MOI.RootDetConeTriangle(n))
@objective(модель, Max, 1 * t + 0)
optimize!(модель)
assert_is_solved_and_feasible(модель)
solve_time_2 = solve_time(модель)
Данная формулировка выводит ещё более краткий результат:
print_active_bridges(модель)
В текущей модели осталось одно неподдерживаемое ограничение. JuMP преобразовало его в комбинацию условий положительной определённости и специального геометрического конуса.
Поскольку решатель Clarabel не работает с геометрическим конусом напрямую, система автоматически заменила его набором степенных конусов.
Для исследования альтернативных подходов можно отключить это автоматическое преобразование с помощью функции remove_bridge, что позволит протестировать другие способы переформулировки ограничения.
remove_bridge(модель, MOI.Bridges.Constraint.GeoMeanToPowerBridge)
optimize!(модель)
assert_is_solved_and_feasible(модель)
Теперь время, затраченное на решение, составляет:
solve_time_3 = solve_time(модель)
При этом, ранее оно составляло:
solve_time_2
Почему время решения разное?
print_active_bridges(модель)
После отключения предыдущего преобразования JuMP применил альтернативный подход: ограничение было переформулировано через конус относительной энтропии, который затем преобразован в поддерживаемые решателем экспоненциальные конусы. Это демонстрирует гибкость системы в выборе эквивалентных математических представлений.
Выводы
Процесс нахождения оптимального эллипсоида продемонстрировал ключевой принцип: одна и та же математическая задача может быть представлена в нескольких эквивалентных формах, и выбор формы может существенно влиять на производительность.
Гибкость JuMP и его способность автоматически подбирать «мосты» (bridges) для преобразования задач освобождает исследователя от рутинной работы. Однако, как мы убедились, понимание этого механизма позволяет вручную выбирать наиболее эффективные пути, отказываясь от автоматических преобразований в пользу тех, что лучше подходят для конкретного решателя.
Таким образом, этот пример служит не только иллюстрацией метода аппроксимации эллипсоидом, но и ценным уроком о том, что глубокое понимание инструментария позволяет управлять процессом оптимизации, а не просто наблюдать за его результатом. Это открывает путь к решению более сложных и масштабных задач.