Классификаторы

Цель данного руководства — продемонстрировать формулирование задач классификации с помощью JuMP.

Задачи классификации связаны с построением функций, называемых классификаторами, которые позволяют эффективно классифицировать данные по двум или более разным наборам. Распространенный вариант их применения — классификация новых точек данных после обучения классификатора на известных данных.

Теоретическая информация и модели для этого руководства взяты из раздела 9.4 работы Ferris2007.

Требуемые пакеты

В этом руководстве используются следующие пакеты.

using JuMP
import DelimitedFiles
import Ipopt
import LinearAlgebra
import Plots
import Random
import Test

Данные и визуализация

Для начала сгенерируем некоторое количество точек для тестирования. Аргумент — это количество двухмерных точек:

function generate_test_points(m; random_seed = 1)
    rng = Random.MersenneTwister(random_seed)
    return 2.0 .* rand(rng, Float64, m, 2)
end

generate_test_points (generic function with 1 method)

Для примера возьмем :

P = generate_test_points(100);

Точки представлены построчно в матрице P. Давайте визуализируем точки с помощью пакета Plots:

plot = Plots.scatter(
    P[:, 1],
    P[:, 2];
    xlim = (0, 2.02),
    ylim = (0, 2.02),
    color = :white,
    size = (600, 600),
    legend = false,
)

Нам нужно разделить точки на два отдельных множества, находящихся по разные стороны разделительной линии. Затем мы маркируем каждую точку в зависимости от того, по какую сторону от линии она находится. На основе меток точек мы покажем, как создать классификатор с использованием модели JuMP. После этого мы сможем проверить, насколько хорошо классификатор воспроизводит исходные метки и границы между ними.

Давайте проведем линию, разделяющую точки на два множества, определив градиент и константу:

w_0, g_0 = [5, 3], 8
line(v::AbstractArray; w = w_0, g = g_0) = w' * v - g
line(x::Real; w = w_0, g = g_0) = -(w[1] * x - g) / w[2];

Функция множественной диспетчеризации Julia позволяет определить векторную форму и форму с одной переменной функции line с одним и тем же именем.

Добавим линию на график:

Plots.plot!(plot, line; linewidth = 5)

Теперь промаркируем точки с учетом того, по какую сторону от линии они находятся. Для предстоящей формулировки JuMP с числовой точки зрения удобно использовать метки +1 и --1.

labels = ifelse.(line.(eachrow(P)) .>= 0, 1, -1)
Plots.scatter!(
    plot,
    P[:, 1],
    P[:, 2];
    shape = ifelse.(labels .== 1, :cross, :xcross),
    markercolor = ifelse.(labels .== 1, :blue, :crimson),
    markersize = 8,
)

Наша цель — показать, что мы можем восстановить линию на основе лишь точек и меток.

Формулировка: линейная машина опорных векторов

Классификатор, известный как линейная машина опорных векторов (SVM), ищет аффинную функцию , которая удовлетворяет условию для всех точек с меткой -1 и условию для всех точек с меткой +1.

Реализующая его квадратичная программа в линейных ограничениях выглядит так:

где — это диагональная матрица меток.

Для параметра положительного штрафа необходимо значение по умолчанию:

C_0 = 100.0;

Формулировка на языке JuMP

Модель JuMP выглядит так:

function solve_SVM_classifier(P::Matrix, labels::Vector; C::Float64 = C_0)
    m, n = size(P)
    model = Model(Ipopt.Optimizer)
    set_silent(model)
    @variable(model, w[1:n])
    @variable(model, g)
    @variable(model, y[1:m] >= 0)
    @objective(model, Min, 1 / 2 * w' * w + C * sum(y))
    D = LinearAlgebra.Diagonal(labels)
    @constraint(model, D * (P * w .- g) .+ y .>= 1)
    optimize!(model)
    Test.@test is_solved_and_feasible(model)
    slack = extrema(value.(y))
    println("Minimum slack: ", slack[1], "\nMaximum slack: ", slack[2])
    classifier(x) = line(x; w = value.(w), g = value(g))
    return model, classifier
end

solve_SVM_classifier (generic function with 1 method)

Результаты

Давайте восстановим значения, определяющие классификатор, решив модель:

_, classifier = solve_SVM_classifier(P, labels)

(A JuMP Model
├ solver: Ipopt
├ objective_sense: MIN_SENSE
│ └ objective_function_type: QuadExpr
├ num_variables: 103
├ num_constraints: 200
│ ├ AffExpr in MOI.GreaterThan{Float64}: 100
│ └ VariableRef in MOI.GreaterThan{Float64}: 100
└ Names registered in the model
  └ :g, :w, :y, Main.classifier)

Получив решение, мы можем спросить: было ли значение штрафной константы «достаточно большим» для этого набора данных? Об этом можно судить отчасти по диапазону ослабляющих переменных. Если он слишком велик, необходимо увеличить штрафную константу.

Давайте построим график решения и посмотрим, как мы справились:

Plots.plot!(plot, classifier; linewidth = 5, linestyle = :dashdotdot)

Как видите, нам удалось восстановить разделительную линию, используя только информацию о точках и их метках.

Неразделимые классы точек

А что если множества точек нельзя четко разделить линией (или гиперплоскостью в более высоких размерностях)? Сработает ли такой метод? Давайте повторим процесс, но на этот раз смоделируем неразделимые классы точек, смешивая несколько соседних точек по разные стороны линии.

nearby_indices = abs.(line.(eachrow(P))) .< 1.1
labels_new = ifelse.(nearby_indices, -labels, labels)
model, classifier = solve_SVM_classifier(P, labels_new)
plot = Plots.scatter(
    P[:, 1],
    P[:, 2];
    xlim = (0, 2.02),
    ylim = (0, 2.02),
    color = :white,
    size = (600, 600),
    legend = false,
)
Plots.scatter!(
    plot,
    P[:, 1],
    P[:, 2];
    shape = ifelse.(labels_new .== 1, :cross, :xcross),
    markercolor = ifelse.(labels_new .== 1, :blue, :crimson),
    markersize = 8,
)
Plots.plot!(plot, classifier; linewidth = 5, linestyle = :dashdotdot)

Наша формулировка JuMP по-прежнему создает классификатор, но он неправильно классифицирует некоторые неразделимые точки.

Мы можем выяснить, какие точки влияют на форму линии, сравнив двойственные значения аффинных ограничений с константой штрафа :

affine_cons = all_constraints(model, AffExpr, MOI.GreaterThan{Float64})
active_cons = findall(isapprox.(dual.(affine_cons), C_0; atol = 0.001))
findall(nearby_indices) ⊆ active_cons

true

Последний оператор говорит о том, что неразделимые точки активно участвуют в определении классификатора. Остальные точки представляют интерес и выделены на графике:

P_active = P[setdiff(active_cons, findall(nearby_indices)), :]
Plots.scatter!(
    plot,
    P_active[:, 1],
    P_active[:, 2];
    shape = :hexagon,
    markersize = 8,
    markeropacity = 0.5,
)

Усложнение задачи: двойственность и ядерный метод

Теперь рассмотрим альтернативную формулировку линейной машины SVM, решив двойственную задачу.

Двойственная программа

Программа, двойственная по отношению к линейной программе SVM, также является квадратичной программой в линейных ограничениях:

Модель JuMP выглядит так:

function solve_dual_SVM_classifier(P::Matrix, labels::Vector; C::Float64 = C_0)
    m, n = size(P)
    model = Model(Ipopt.Optimizer)
    set_silent(model)
    @variable(model, 0 <= u[1:m] <= C)
    D = LinearAlgebra.Diagonal(labels)
    @objective(model, Min, 1 / 2 * u' * D * P * P' * D * u - sum(u))
    @constraint(model, con, sum(D * u) == 0)
    optimize!(model)
    Test.@test is_solved_and_feasible(model)
    w = P' * D * value.(u)
    g = dual(con)
    classifier(x) = line(x; w = w, g = g)
    return classifier
end

solve_dual_SVM_classifier (generic function with 1 method)

Мы восстанавливаем градиентный вектор линии , задавая , а константу линии — как двойственное значение одиночного аффинного ограничения.

В двойственной задаче меньше переменных и ограничений, поэтому зачастую решить двойственную форму может быть проще.

Мы можем проверить, удалось ли двойственной форме восстановить классификатор:

classifier = solve_dual_SVM_classifier(P, labels)
Plots.plot!(plot, classifier; linewidth = 5, linestyle = :dash)

Ядерный метод

Линейные методы SVM не ограничиваются нахождением разделяющих гиперплоскостей в исходном пространстве набора данных. Можно сначала преобразовать обучающие данные с помощью нелинейного отображения, применить метод, а затем отобразить гиперплоскость обратно в исходное пространство.

Фактические данные, описывающие набор точек, хранятся в матрице , но в случае двойственной программы на самом деле важна матрица Грама , выражающая попарное сравнение (внутреннее произведение) каждого вектора точек. Отсюда следует, что любое отображение множества точек достаточно определить на уровне попарных отображений между точками. Такие отображения называются ядерными функциями:

где правая часть применяет некоторое преобразование , после чего находится внутреннее произведение в этом пространстве.

На практике можно избежать явного указания , а вместо этого определить ядерную функцию непосредственно между парами векторов. Такой вариант использования ядерной функции без знания отображения называется ядерным методом (или иногда ядерным трюком).

Классификатор с использованием гауссова ядра

Продемонстрируем применение гауссова ядра или радиальной базисной функции:

для некоторого положительного параметра .

k_gauss(s::Vector, t::Vector; μ = 0.5) = exp(-μ * LinearAlgebra.norm(s - t)^2)

k_gauss (generic function with 1 method)

Если дана матрица точек, выраженных построчно, и ядро, следующая функция возвращает преобразованную матрицу , которая заменяет :

function pairwise_transform(kernel::Function, P::Matrix{T}) where {T}
    m, n = size(P)
    K = zeros(T, m, m)
    for j in 1:m, i in 1:j
        K[i, j] = K[j, i] = kernel(P[i, :], P[j, :])
    end
    return LinearAlgebra.Symmetric(K)
end

pairwise_transform (generic function with 1 method)

Теперь все готово для определения задачи оптимизации. Нам необходимо предоставить ядерную функцию, которая будет использоваться в задаче. Обратите внимание, что все дополнительные именованные аргументы (например, значения параметров) передаются в ядро.

function solve_kernel_SVM_classifier(
    kernel::Function,
    P::Matrix,
    labels::Vector;
    C::Float64 = C_0,
    kwargs...,
)
    m, n = size(P)
    K = pairwise_transform(kernel, P)
    model = Model(Ipopt.Optimizer)
    set_silent(model)
    @variable(model, 0 <= u[1:m] <= C)
    D = LinearAlgebra.Diagonal(labels)
    con = @constraint(model, sum(D * u) == 0)
    @objective(model, Min, 1 / 2 * u' * D * K * D * u - sum(u))
    optimize!(model)
    Test.@test is_solved_and_feasible(model)
    u_sol, g_sol = value.(u), dual(con)
    function classifier(v::Vector)
        return sum(
            D[i, i] * u_sol[i] * kernel(P[i, :], v; kwargs...) for i in 1:m
        ) - g_sol
    end
    return classifier
end

solve_kernel_SVM_classifier (generic function with 1 method)

На этот раз мы не восстанавливаем градиентный вектор линии напрямую. Вместо этого мы вычисляем классификатор с помощью функции:

где — это вектор строки из .

Набор данных шахматной доски

Продемонстрируем этот нелинейный метод на наборе данных шахматной доски.

filename = joinpath(@__DIR__, "data", "checker", "checker.txt")
checkerboard = DelimitedFiles.readdlm(filename, ' ', Int)
labels = ifelse.(iszero.(checkerboard[:, 1]), -1, 1)
B = checkerboard[:, 2:3] ./ 100.0  # изменяем масштаб квадратов на [0,2] x [0,2].
plot = Plots.scatter(
    B[:, 1],
    B[:, 2];
    color = ifelse.(labels .== 1, :white, :black),
    markersize = ifelse.(labels .== 1, 4, 2),
    size = (600, 600),
    legend = false,
)

Способен ли этот метод сгенерировать явно нелинейную поверхность? Давайте решим квадратичную задачу на основе гауссова ядра со следующими параметрами:

classifier = solve_kernel_SVM_classifier(k_gauss, B, labels; C = 1e5, μ = 10.0)
grid = [[x, y] for x in 0:0.01:2, y in 0:0.01:2]
grid_pos = [Tuple(g) for g in grid if classifier(g) >= 0]
Plots.scatter!(plot, grid_pos; markersize = 0.2)

Оказывается, ядерный метод может хорошо работать в качестве нелинейного классификатора.

Результат довольно сильно зависит от выбора параметров: большие значения допускают более сложную границу, тогда как меньшие значения приводят к более гладкой границе для классификатора. Определение более эффективной ядерной функции и выбор параметров осуществляются в процессе перекрестной проверки по отношению к набору данных: для проверки оптимального выбора параметров относительно статистической величины погрешности используются различные наборы для тестирования, обучения и настройки.