Добро пожаловать! Как подопечный в программе наставничества ChiPy, я напишу несколько сообщений в блоге о своем проекте, который должен был узнать, как реализовать пару алгоритмов машинного обучения для выполнения на видеокарте. В этом сообщении блога я представлю несколько фундаментальных концепций машинного обучения на примере логистической регрессии, а также код с простой реализацией на Python и OpenCL, взаимодействующий с PyOpenCL. Этот пост предназначен для широкой аудитории; Если вы новичок в машинном обучении, возможно, стоит прочитать начало и самостоятельно повозиться с кодом. Если вы здесь ради OpenCL, просто прыгайте вниз. И если у вас есть отзывы, хорошие или плохие, дайте мне знать! Если вы хотите увидеть код, который я использую, вы можете заглянуть в мой репозиторий на github. Поехали!



Модель логистической регрессии

Есть много случаев, когда кажется естественным попытаться смоделировать результат события как вероятность того, что событие произойдет, поскольку мы часто записываем результат какого-либо события в виде двоичных данных. Кроме того, у нас могут быть некоторые данные, которые, по нашему предположению, должны коррелировать с наступлением указанного события. В медицине такие сценарии возникают часто.

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

Для логистической регрессии (которая является типом регрессионной модели в классе обобщенных линейных моделей или GLM) нам нужно несколько вещей. Первый - это распределение вероятностей, которое отражает структуру выходных данных, которые мы пытаемся оценить, и способ связать наши входные данные со средним значением этого распределения. Я собираюсь просто дать вам эти уравнения здесь, но я настоятельно рекомендую прочитать приложение в конце этого сообщения в блоге для более полного вывода модели. Для системы, которая генерирует двоичные результаты с некоторой вероятностью от 0 до 1, что «положительное» событие произойдет, естественным выбором распределения вероятностей является распределение Бернулли. Это функция плотности вероятности (pdf) выглядит так:

Где пи - вероятность того, что у = 1. Однако нам нужно найти способ сопоставить произвольно большое (конечное) значение с диапазоном (0,0, 1,0), чтобы оно могло отражать значение числа пи или нашу веру в то, что y = 1. Функция, которую мы здесь будем использовать, - это логистическая функция (или логит, или сигмоид). Его форма выглядит так:

Подумайте, что происходит, когда бета очень большая или очень маленькая. Если бета большая, дробь приближается к 1, а если бета мала, дробь приближается к 0, что делает это отличным способом сопоставить некоторые входные данные (сохраненные как значение бета) с диапазоном (0,0, 1,0), который мы можем интерпретировать как вероятность того, что y = 1. Давайте формализуем это и представим, какой должна быть бета-версия:

где X - наши данные, а тета - линейные коэффициенты каждой характеристики данных. С этой структурой мы можем начать думать о том, как мы на самом деле подойдем к вычислению значений тета, учитывая набор входных и выходных данных, поскольку это позволит нам «изучить» карту между нашими входными функциями (X) и выходные данные (y).

Подгонка и оптимизация модели

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

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

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

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

Как видите, в строке 2 алгоритма «Learn Something» первое, что нам нужно сделать, это попытаться понять, каким должен быть ответ. В случае логистической регрессии у нас уже есть это из первой части сообщения в блоге. Сейчас нам нужно определить, как мы 1) собираемся измерять ошибку и 2) предлагать корректировки модели, чтобы уменьшить ошибку в следующий раз. Обычно мы это делаем, настраивая функцию, которую хотим оптимизировать. В литературе по машинному обучению это часто называют функцией потерь или функцией стоимости. В идеале эта функция является выпуклой (или вогнутой), что означает, что она гарантированно имеет минимальное (или максимальное) значение - эта структура особенно полезна, поскольку позволяет нам использовать производную этой функции для поиска максимума или минимума. Напомним, что в исчислении минимум выпуклой функции для некоторого параметра можно найти, установив частную производную функции по рассматриваемому параметру равной нулю, а затем решив для этого параметра. Итак, какую функцию мы должны использовать для измерения нашей производительности. Хорошее место для начала - функция правдоподобия:

В основном то, что мы хотим сделать здесь, - это максимизировать эту функцию по набору обучающих выборок (пары (x, y) от 1 до N), что означает, что мы корректируем значения тета, которые увеличили бы общее измерение правдоподобия с учетом нашего данные. Однако первое, что нам нужно сделать, это записать правую часть уравнения через X и тета; К счастью, у нас есть результат, полученный ранее, так что давайте просто подключим его.

Хорошо, это выглядит лучше, но продукт по-прежнему затрудняет вычисление и решение проблемы. Однако мы можем упростить это, логарифмируя обе части, превращая произведение в сумму.

И для удобства превратим это в задачу минимизации, отрицая каждую часть уравнения:

Большой! А теперь давайте займемся алгеброй и вычислим производную по тэте для этого плохого парня. Шучу, я думаю, что уже достаточно замучил вас сегодня (и я не хочу печатать все это латексом: P), поэтому приведу уравнение для производной:

Следует отметить, что я еще не упомянул, что я использую матричную нотацию для X и theta, поэтому для большей ясности, если бы у вас было j функций (столбцов) в X, мы бы представили theta как вектор j +1 числа (мы добавляем единицу как термин «перехват» или «смещение»). Следовательно, приведенное выше уравнение для каждого из наших j + 1 тета-параметров будет выглядеть так:

Итак, все, что нам нужно сделать сейчас, это установить это значение равным нулю и решить для теты, верно? К сожалению, не все так просто. Учитывая нелинейность, присутствующую в уравнении, мы не можем аналитически решить для theta. Но мы можем использовать свойства производной, чтобы «подтолкнуть» значения тета в правильном направлении на протяжении многих итераций. Формальное название этого процесса - градиентный спуск (так как мы спускаемся по градиенту к оптимуму). Алгоритм градиентного спуска в этом случае будет выглядеть так:

В приведенном выше псевдокоде альфа относится к «скорости обучения» или тому, насколько велика тета на каждом шаге. На практике это может быть сложно установить, и существует множество схем, в которых используются адаптивные методы, чтобы поддерживать настройку теты достаточно большой, чтобы мы снизили стоимость, но не настолько большой, чтобы мы могли выйти за пределы глобального оптимума. В целом алгоритм довольно прост и поразительно похож на глупый алгоритм «Выучить что-то», не так ли? В следующем разделе мы посмотрим, как это реализовано в Python.

Установка данных на ЦП

Прежде чем приступить к какой-либо подгонке, нам действительно нужны некоторые данные. Поскольку я хочу убедиться, что код работает, я буду генерировать синтетические данные, которые определенно можно смоделировать с помощью логистической регрессии. Код для подгонки модели на ЦП прямо размещен ниже, однако не похоже, что у Medium есть какие-либо хорошие решения, чтобы обойти ужасное форматирование кода. Ссылку на блокнот jupyter можно найти здесь.

import os
import pyopencl as pcl
import numpy as np
import scipy.stats as ss
import pandas as pd
import math
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
x0_1 = ss.norm(loc=10.0, scale=2.0)
x0_0 = ss.norm(loc=7.0, scale=2.0)
x1_1 = ss.norm(loc=5.0, scale=3.0)
x1_0 = ss.norm(loc=-5.0, scale=3.0)
nsamps=1024
X_1 = pd.DataFrame(index=range(nsamps), 
                   columns=['x0','x1', 'y'])
X_0 = pd.DataFrame(index=range(nsamps), 
                   columns=['x0','x1', 'y'])
X_1.loc[:, 'x0'] = x0_1.rvs(size=(nsamps,)).astype(np.float32)
X_1.loc[:, 'x1'] = x1_1.rvs(size=(nsamps,)).astype(np.float32)
X_1.loc[:, 'y'] = np.ones(shape=(nsamps,)).astype(np.float32)
X_0.loc[:, 'x0'] = x0_0.rvs(size=(nsamps,)).astype(np.float32)
X_0.loc[:, 'x1'] = x1_0.rvs(size=(nsamps,)).astype(np.float32)
X_0.loc[:, 'y'] = np.zeros(shape=(nsamps,)).astype(np.float32)
X_all = pd.concat((X_1, X_0), ignore_index=True)
X_all = X_all.reindex(np.random.permutation(X_all.index))
X = X_all.loc[:, ['x0', 'x1']]
y = X_all.loc[:,'y']
def sig(X, theta):
    lin = X.dot(theta)
    sig = 1.0 / (1.0 + np.exp(-lin))
    return sig
 
    
def lr_cost(X, theta, y):
    est = sig(X, theta)
    log_est = np.log(est)
    cost = y*log_est + (1-y) * (1-log_est)
    cost *= -1
    return cost
def cost_ss(est, actual):
    cost = ((est - actual)**2.0)
    cost /= 2.0
    return cost
def grad(X, theta, y):
    diff = y - sig(X, theta)
    g = np.zeros(X.shape)
    for i in range(X.shape[0]):
        for j  in range(X.shape[1]):
            g[i,j] = diff[i] * X[i,j]
    return g
def fit_params(X, y, theta):
    tol = 1e-5
    learning_rate = 1e-2
    costs = [np.inf]
    for i in range(0, 1000):
        my_cost = lr_cost(X, theta, y)
        my_cost = my_cost.sum()
        costs.append(my_cost)
        if abs(my_cost - costs[-2]) < tol:
            break
        else:
            my_gradient = grad(X, theta, y)
            my_gradient = my_gradient.mean(axis=0)
            theta = theta + learning_rate * my_gradient
            
    return theta, costs
X_train, X_test, y_train, y_test = train_test_split(X, 
                                                    y, 
                                                    train_size=0.6)
new_t = np.random.normal(size=(X_train.shape[1], )).astype(np.float32)
fitted_theta, outcost = fit_params(X_train.values, 
                                   y_train.values, 
                                   new_t)
print('cost on test set: ', lr_cost(X_test.values, fitted_theta, y_test.values).sum())
print('average sum of square error on test set: ',cost_ss(sig(X_test.values, fitted_theta), y_test.values).sum()/X_test.shape[0])
plt.figure(figsize=(10,6))
plt.plot(outcost[1:])
plt.xlabel('Iteration')
plt.ylabel('Cost')
plt.title('Cost During Fitting - CPU')
plt.show()

Установка данных на GPU

Итак, ребята, мы наконец-то здесь, и это суть моего проекта Программа наставничества ChiPy. Адаптация моделей машинного обучения для выполнения на графическом процессоре привела к значительному увеличению производительности для ряда проблем в области искусственного интеллекта, компьютерного зрения и обработки естественного языка, и интерес к этой области растет. Не верите мне? Посмотрите на курс акций NVIDIA за последний год. Я пишу код на Python около 6 лет, но я никогда не интегрировал код с других языков с Python, и я подумал, что это даст мне возможность узнать об этом процессе, а также о программировании на GPU.

У меня есть хорошо прокомментированный код, которым я могу поделиться для вычислений на GPU с OpenCL и PyOpenCL для логистической регрессии, но, учитывая форматирование Medium, я думаю, что было бы лучше оставить ссылку на этот блокнот здесь. Однако я использую это место, чтобы немного рассказать о том, как работает OpenCL.

OpenCL - это сокращение от «Open Computing Language» - он очень похож на C99 и специально разработан для выполнения в высокой степени параллелизма. Это может значительно ускорить большое количество вычислений, в основном любой алгоритм, большая часть работы которого не зависит от ранее синхронизированных результатов. Например, большие матричные умножения, обработка отдельных элементов в большом списке, попарные сравнения длинных списков взаимодействующих элементов и классические задачи «N тел» в физическом моделировании, таком как молекулярная динамика. Распределяя рабочую нагрузку между большим количеством рабочих и позволяя всем им работать параллельно, мы можем сэкономить много времени. Так как же на самом деле выглядит структура программы OpenCL?

Для любой программы OpenCL в основном есть два фрагмента кода: один, который работает на «хосте», и ядро, которое запускается на «устройстве». Хост-программа написана на любом языке, который может взаимодействовать с библиотекой OpenCL, что означает, что это обычно программа C, C ++ или даже Python. В этом случае ядро ​​представляет собой специализированный код OpenCL, который выполняется на устройстве (которое в моем случае является выделенной видеокартой на моем ноутбуке, но это может быть другой процессор). Код хоста управляет общим потоком и логикой программы, которая будет считывать и записывать данные, отправлять код и данные в GPU, обеспечивая выполнение графическим процессором заданных инструкций и копирование результатов из памяти GPU обратно в память хоста. OpenCL также предоставляет чрезвычайно гибкий интерфейс для синхронизации данных и инструкций на множестве различных устройств. Хорошая диаграмма для понимания структуры OpenCL показана ниже:

Иерархия рабочего пространства OpenCL следующая: контекст (хост-программа может запускать несколько контекстов), устройство (контекст может содержать несколько устройств), очередь команд (устройство может иметь несколько очередей команд) и ядро ​​(несколько ядер могут быть последовательно добавлен в очередь команд).

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

__kernel void square(__global float * in, __global float * out)
{
   int gid = get_global_id(0);
   out[gid] = in[gid] * in[gid];
}

Обратите внимание, что мы никогда не настраиваем цикл for, который делает что-то вроде получения индекса массива. Когда мы настраиваем ядро, указав глобальный размер рабочего пространства (который в данном случае может быть длиной массива), OpenCL запустит большое количество этих ядер на графическом процессоре, которые будут работать независимо в своем собственном подразделе. массива.

Я надеюсь, что в этом сообщении в блоге было что предложить всем, независимо от их уровня знаний о машинном обучении, Python и OpenCL. Код, который я предоставляю в записной книжке jupyter, работает, но ни в коем случае не оптимизирован, но я буду продолжать улучшать его в ходе программы наставничества по python. В своих следующих сообщениях в блоге я планирую показать код OpenCL, который оптимизирован и может обрабатывать большие наборы данных с высокой производительностью. Будьте на связи!

Приложение: Вывод модели логистической регрессии

Я рада, что вы решили приехать сюда и посмотреть, почему все мое махание руками оправдано. История здесь начинается с семейства вероятностных распределений, называемого экспоненциальным семейством. Все распределения в экспоненциальном семействе имеют несколько общих черт, а именно то, что все они имеют функцию плотности вероятности (pdf) в форме:

где beta - параметр локализации (~ mean), phi - параметр масштаба (~ дисперсия), а y - наши данные. Многие из наиболее часто используемых распределений вероятностей следуют этой форме, например, распределение Гаусса, распределение Пуассона и биномиальное распределение. У экспоненциального семейства распределений есть несколько важных свойств, касающихся среднего и дисперсии. Я не собираюсь выводить их сейчас, но нам нужно будет использовать их позже, когда мы построим модель логистической регрессии:

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

Итак, какое распределение мы должны использовать в сценарии, в котором мы хотим смоделировать двоичный результат происходящего (или не происходящего) события? Хорошим выбором для этого будет распределение Бернулли, у которого есть файл в формате pdf, который выглядит так:

где пи - наша вера в то, что событие y, которое может принимать значение 1 или 0, произойдет или будет равно 1. Например, в случаях, когда мы думаем, что вероятность того, что наше интересующее событие произойдет, такая же, как и у события не происходит (например, подбрасывание честной монеты) пи будет равно 0,5. Хорошо, это здорово, но не похоже, что это часть семейства экспоненциальных распределений! По крайней мере, пока. Давай займемся алгеброй.

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

Следовательно, решение этого уравнения относительно пи дает нам:

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

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

а также

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

Подумайте, что происходит, когда мы меняем значение беты: когда оно действительно велико, дробь приближается к значению 1, а когда она очень мала, дробь приближается к значению нуля. Логистическая (или логит, или сигмоидальная) функция дает нам именно то, что мы хотим, способ сопоставить любое сколь угодно большое или малое число с диапазоном (0,1). Что касается более технических замечаний в литературе, вы можете часто встретить обобщенные линейные модели, описанные как имеющие функцию «связи», которая присуща каждому распределению, часто записываемая как таковая:

Идея здесь заключается в том, что функция g предоставляет способ описать карту между некоторым значением eta и ожидаемым значением X. Однако в нашем случае мы действительно хотим обратное к g - это позволит нам ввести некоторое значение для eta и получите ожидаемое значение распределения X (которое в нашем случае является распределением Бернулли). Но подождите, он у нас уже есть для распределения Бернулли! У нас есть функция, которая позволяет нам выразить пи как функцию от бета (предпоследнее уравнение). Последнее, что нам нужно сделать, это как использовать это с нашими данными. Самый простой способ - построить линейную модель, которая связывает характеристики наших данных и веса этих характеристик с ожидаемым значением распределения. В нашем случае это будет выглядеть так:

Взгляните на уравнение 5. Мы в основном уже решили его, когда выяснили, какой была производная кумулянтной функции для распределения Бернулли. Замена, которая дает нам:

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

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

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