Логистическая регрессия

На этой странице приводится пошаговое руководство по использованию алгоритма логистической регрессии в Julia с применением Flux. После этого создается простая модель логистической регрессии без использования Flux, и ее функциональные компоненты сравниваются с реализацией Flux.

Сначала импортируем необходимые пакеты Julia.

julia> using Flux, Statistics, MLDatasets, DataFrames, OneHotArrays

Набор данных

Сначала импортируем набор данных из MLDatasets.jl. Мы будем использовать набор Iris, который содержит данные о трех разных видах ириса (Iris). Набор состоит из 150 точек данных (x), у каждой из которых четыре признака. Каждый x сопоставлен с y, названием вида ириса (Iris). Следующий код скачивает набор данных Iris при первом запуске.

julia> Iris()
dataset Iris:
  metadata   =>    Dict{String, Any} with 4 entries
  features   =>    150×4 DataFrame
  targets    =>    150×1 DataFrame
  dataframe  =>    150×5 DataFrame

julia> x, y = Iris(as_df=false)[:];

Давайте посмотрим на набор данных:

julia> y
1×150 Matrix{InlineStrings.String15}:
 "Iris-setosa"  "Iris-setosa"  …  "Iris-virginica"  "Iris-virginica"

julia> x |> summary
"4×150 Matrix{Float64}"

Значения y здесь соответствуют типу ириса; всего имеется 150 точек данных. Значения x описывают длину чашелистика, ширину чашелистика, длину лепестка и ширину лепестка (все в см (cm)) 150 растений ириса (поэтому матрица имеет размер 4×150). У разных типов ириса разные длина и ширина чашелистиков и лепестков, и эта связь подчиняется определенной закономерности. С учетом этого можно обучить простой классификатор, который выводит вид ириса на основе длины и ширины чашелистиков и лепестков в качестве входов.

Нашим следующим шагом будет преобразование этих данных в форму, пригодную для подачи в модель машинного обучения. Значения x выстроены в виде матрицы и в идеале должны быть преобразованы в тип Float32 (см. раздел Советы по производительности), но метки должны быть в прямой унитарной кодировке. В этой отличной ветке обсуждения можно найти информацию о различных приемах прямой унитарной кодировки данных с использованием внешних пакетов Julia и без них.

julia> x = Float32.(x);

julia> y = vec(y);

julia> custom_y_onehot = unique(y) .== permutedims(y)
3×150 BitMatrix:
 1  1  1  1  1  1  1  1  1  1  1  1  1  …  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     1  1  1  1  1  1  1  1  1  1  1  1

Ту же операцию можно выполнить с помощью функции onehotbatch из пакета OneHotArrays. Мы используем оба этих выхода параллельно, чтобы показать, насколько интуитивно понятна работа с FluxML.

julia> const classes = ["Iris-setosa", "Iris-versicolor", "Iris-virginica"];

julia> flux_y_onehot = onehotbatch(y, classes)
3×150 OneHotMatrix(::Vector{UInt32}) with eltype Bool:
 1  1  1  1  1  1  1  1  1  1  1  1  1  …  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅     ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅     1  1  1  1  1  1  1  1  1  1  1  1

Данные готовы. На следующем шаге необходимо создать для них классификатор.

Построение модели

Модель логистической регрессии математически выражается так:

где W — матрица весов, b — вектор смещений, а σ — какая-либо функция активации. В нашем случае мы воспользуемся функцией активации softmax, так как будем выполнять задачу многоклассовой классификации.

julia> m(W, b, x) = W*x .+ b
m (generic function with 1 method)

Обратите внимание, что у данной модели нет функции активации, но мы вернемся к этому позже.

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

julia> W = rand(Float32, 3, 4);

julia> b = [0.0f0, 0.0f0, 0.0f0];

Теперь наша модель может принять весь набор данных и спрогнозировать класс каждого x за один проход. Но нам нужно убедиться в том, что модель выводит вероятности принадлежности входов к соответствующим классам. Так как у модели три входа, определяется вероятность принадлежности каждого из них к определенному классу.

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

Функция softmax пропорционально уменьшает выходы до значений вероятности так, что сумма всех конечных выходов равна 1. Давайте реализуем это в Julia.

julia> custom_softmax(x) = exp.(x) ./ sum(exp.(x), dims=1)
custom_softmax (generic function with 1 method)

Реализация выглядит достаточно просто. Обратите внимание, что мы указываем dims=1 в функции sum для вычисления суммы вероятностей в каждом столбце. Напомним, что выходом модели является матрица 3 × 150 (спрогнозированные y), в которой каждый столбец представляет собой выход для соответствующего входа.

Давайте объединим эту функцию softmax с моделью, чтобы построить полную модель custom_model.

julia> custom_model(W, b, x) = m(W, b, x) |> custom_softmax
custom_model (generic function with 1 method)

Проверим, работает ли модель.

julia> custom_model(W, b, x) |> size
(3, 150)

Работает! Проверим, работает ли функция softmax.

julia> all(0 .<= custom_model(W, b, x) .<= 1)
true

julia> sum(custom_model(W, b, x), dims=1)
1×150 Matrix{Float32}:
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  …  1.0  1.0  1.0  1.0  1.0  1.0  1.0

Каждое выходное значение находится в диапазоне от 0 до 1, и сумма значений в каждом столбце равна 1.

Давайте преобразуем нашу модель custom_model в модель Flux. Flux предоставляет пользователям очень стройный API, пользоваться которым практически так же удобно, как писать собственный код.

Обратите внимание, что переменные вида flux_* в этом руководстве являются универсальными, то есть их можно использовать и с другими сходными наборами данных, а переменные вида custom_* применимы только в рамках задачи данного руководства.

julia> flux_model = Chain(Dense(4 => 3), softmax)
Chain(
  Dense(4 => 3),                        # 15 параметров
  NNlib.softmax,
)

Dense(4 => 3) — это слой с четырьмя входами (четыре признака у каждой точки данных) и тремя выходами (три класса или метки). Этот слой соответствует определенной выше математической модели. На внутреннем уровне Flux вычисляет выход с использованием того же выражения, но на этот раз нам не нужно самим инициализировать параметры — Flux сделает это за нас.

Используемая здесь функция softmax из библиотеки NNLib.jl повторно экспортируется Flux. Наконец, Flux предоставляет пользователям структуру Chain, позволяющую легко накладывать уровни друг на друга.

Доступ к весам и смещениям модели можно получить следующим образом.

julia> flux_model[1].weight, flux_model[1].bias
(Float32[0.78588694 -0.45968163 -0.77409476 0.2358028; -0.9049773 -0.58643705 0.466441 -0.79523873; 0.82426906 0.4143493 0.7630932 0.020588955], Float32[0.0, 0.0, 0.0])

Теперь можно передать сразу все данные, у каждой точки которых четыре признака (входа).

Потери и точность

На следующем шаге следует определить некоторые количественные значения для модели, которые будут максимизироваться или минимизироваться в процессе обучения. Этими значениями будут функция потерь и показатель точности.

Начнем с определения функции потерь: logitcrossentropy.

julia> custom_logitcrossentropy(ŷ, y) = mean(.-sum(y .* logsoftmax(ŷ; dims = 1); dims = 1));

Теперь можно включить custom_logitcrossentropy в функцию, которая принимает параметры модели, x и y, и возвращает значение потерь.

julia> function custom_loss(W, b, x, y)
           ŷ = custom_model(W, b, x)
           custom_logitcrossentropy(ŷ, y)
       end;

julia> custom_loss(W, b, x, custom_y_onehot)
1.1714406827505623

Функция потерь работает!

Flux предоставляет множество полезных функций потерь в лаконичной форме. По сути, определенная выше функция custom_logitcrossentropy была взята напрямую из Flux. Имеющиеся во Flux функции включают в себя проверки исправности, обеспечивают высокую производительность и хорошо взаимодействуют со всей экосистемой FluxML.

julia> function flux_loss(flux_model, x, y)
           ŷ = flux_model(x)
           Flux.logitcrossentropy(ŷ, y)
       end;

julia> flux_loss(flux_model, x, flux_y_onehot)
1.2156688659673647

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

Эту задачу можно разделить на две задачи:

  1. определение индекса максимального элемента в каждом столбце выходной матрицы;

  2. преобразование этого индекса в имя класса.

Индекс максимального элемента должен определяться по столбцам (напомним, что каждый столбец — это выход одной точки данных x). С этой целью можно воспользоваться функцией Julia argmax.

julia> argmax(custom_y_onehot, dims=1)  # вычисляем декартов индекс максимального элемента в каждом столбце
1×150 Matrix{CartesianIndex{2}}:
 CartesianIndex(1, 1)  CartesianIndex(1, 2)  …  CartesianIndex(3, 150)

julia> max_idx = [x[1] for x in argmax(custom_y_onehot; dims=1)]
1×150 Matrix{Int64}:
 1  1  1  1  1  1  1  1  1  1  1  1  1  …  3  3  3  3  3  3  3  3  3  3  3  3

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

julia> function custom_onecold(custom_y_onehot)
           max_idx = [x[1] for x in argmax(custom_y_onehot; dims=1)]
           vec(classes[max_idx])
       end;

julia> custom_onecold(custom_y_onehot)
150-element Vector{String}:
 "Iris-setosa"
 "Iris-setosa"
 "Iris-setosa"
 "Iris-setosa"
 "Iris-setosa"
 "Iris-setosa"
 "Iris-setosa"
 "Iris-setosa"
 "Iris-setosa"
 "Iris-setosa"
 ⋮
 "Iris-virginica"
 "Iris-virginica"
 "Iris-virginica"
 "Iris-virginica"
 "Iris-virginica"
 "Iris-virginica"
 "Iris-virginica"
 "Iris-virginica"
 "Iris-virginica"

Работает!

Flux предоставляет пользователям функцию onecold, так что писать ее самостоятельно не нужно. Сравним нашу функцию custom_onecold с Flux.onecold.

julia> istrue = Flux.onecold(flux_y_onehot, classes) .== custom_onecold(custom_y_onehot);

julia> all(istrue)
true

Обе функции работают одинаково!

Теперь перейдем к показателю точности (accuracy) и запустим с ним необученную модель custom_model.

julia> custom_accuracy(W, b, x, y) = mean(custom_onecold(custom_model(W, b, x)) .== y);

julia> custom_accuracy(W, b, x, y)
0.3333333333333333

Для определения функции точности можно было бы также воспользоваться встроенными возможностями Flux.

julia> flux_accuracy(x, y) = mean(Flux.onecold(flux_model(x), classes) .== y);

julia> flux_accuracy(x, y)
0.24

Обучение модели

Давайте обучим нашу модель с помощью классического алгоритма градиентного спуска. Согласно алгоритму градиентного спуска веса и смещения должны итеративно обновляться по следующим математическим уравнениям:

Здесь W — матрица весов, b — вектор смещений,  — скорость обучения,  — производная функции потерь относительно веса, а  — производная функции потерь относительно смещения.

Производные вычисляются с помощью инструмента автоматического дифференцирования. Во Flux для этой цели применяется Zygote.jl. Так как Zygote.jl — это автономный пакет Julia, его можно использовать и без Flux. Дополнительные сведения см. в документации по Zygote.jl.

В первую очередь следует получить градиент функции потерь относительно весов и смещений. Flux повторно экспортирует функцию gradient из Zygote, поэтому импортировать Zygote явным образом для ее использования не требуется. gradient принимает функцию и ее аргументы и возвращает кортеж, содержащий ∂f/∂x для каждого аргумента x. Давайте передадим custom_loss и аргументы, которые требуются для custom_loss, в gradient. Для выполнения градиентного спуска потребуются производные функции потерь (custom_loss) относительно весов (∂f/∂w) и смещения (∂f/∂b), но мы можем игнорировать частные производные функции потерь (custom_loss) относительно x (∂f/∂x) и y в прямой унитарной кодировке (∂f/∂y).

julia> dLdW, dLdb, _, _ = gradient(custom_loss, W, b, x, custom_y_onehot);

Теперь можно обновить параметры по алгоритму градиентного спуска:

julia> W .= W .- 0.1 .* dLdW;

julia> b .= b .- 0.1 .* dLdb;

Параметры обновлены. Теперь можно проверить значение пользовательской функции потерь:

julia> custom_loss(W, b, x, custom_y_onehot)
1.164742997664842

Потери уменьшились! Давайте включим логику обучения в функцию.

julia> function train_custom_model()
           dLdW, dLdb, _, _ = gradient(custom_loss, W, b, x, custom_y_onehot)
           W .= W .- 0.1 .* dLdW
           b .= b .- 0.1 .* dLdb
       end;

Мы можем ввести функцию обучения в цикл и обучить модель в течение большего числа эпох. Цикл можно адаптировать под потребности пользователя, а условия указывать в виде обычного кода Julia. Здесь мы обучим модель в течение максимум 500 эпох, но во избежание переобучения модели цикл будет прерван, как только значение точности достигнет или превысит 0.98.

julia> for i = 1:500
            train_custom_model();
            custom_accuracy(W, b, x, y) >= 0.98 && break
       end

julia> @show custom_accuracy(W, b, x, y);
custom_accuracy(W, b, x, y) = 0.98

Все работает! Наша модель достигла точности 0.98! Давайте посмотрим на потери.

julia> custom_loss(W, b, x, custom_y_onehot)
0.6520349798243569

Как и ожидалось, потери также уменьшились! Повторим те же действия с моделью flux_model.

Мы можем написать похожий цикл для модели flux_model и обучить ее аналогичным образом.

julia> flux_loss(flux_model, x, flux_y_onehot)
1.215731131385928

julia> function train_flux_model()
           dLdm, _, _ = gradient(flux_loss, flux_model, x, flux_y_onehot)
           @. flux_model[1].weight = flux_model[1].weight - 0.1 * dLdm[:layers][1][:weight]
           @. flux_model[1].bias = flux_model[1].bias - 0.1 * dLdm[:layers][1][:bias]
       end;

julia> for i = 1:500
            train_flux_model();
            flux_accuracy(x, y) >= 0.98 && break
       end

Посмотрим на точность и значение потерь:

julia> @show flux_accuracy(x, y);
flux_accuracy(x, y) = 0.98

julia> flux_loss(flux_model, x, flux_y_onehot)
0.6952386604624324

Итоговые потери и точность очень близки.


Подведем итог: в этом руководстве мы увидели, как выполнить алгоритм логистической регрессии в Julia с использованием пакета Flux и без него. Мы начали с импорта классического набора данных Iris и перевели метки в прямую унитарную кодировку. Затем мы самостоятельно определили модель, функцию потерь и точность.

Наконец, мы обучили модель, вручную написав алгоритм градиентного спуска и оптимизировав потери. Что любопытно, мы реализовали большинство функций самостоятельно, а затем параллельно сравнили их с функциями из Flux!

Впервые опубликовано 1 апреля 2023 г. Автор: Саранш Чопра (Saransh Chopra).