Создание программы OpenCL

Работу по переносу вычислений мы начнем с создания программы OpenCL. Выбор такого подхода весьма очевиден. Мы уже организовали процесс средствами MQL5. Следовательно, весь процесс операций нам ясен и понятен. Операции вычислений будут осуществляться в программе OpenCL. В основной программе нам предстоит организовать процесс передачи данных и вызов программы OpenCL. Последний процесс проще организовать, когда нам уже будет известно, какие данные и в какой kernel необходимо передать.

Код программы запишем в отдельный файл opencl_program.cl. Имя и расширения файла может быть любым, так как в последующем мы будем подгружать ее в коде основной программы в качестве ресурса. Я же использую расширение *.cl как общепринятое расширение для обозначения программ OpenCL. В целом использование стандартных расширений облегчает чтение проектов со сложной структурой файлов.

Для прямого прохода по аналогии с реализацией MQL5 создадим kernel PerceptronFeedForward. Для проведения операций нам потребуются в качестве исходных данных вектор результатов предыдущего слоя (inputs) и матрица весовых коэффициентов (weights). Результат операций запишем в вектор результатов (outputs). Кроме массивов данных нам потребуются количество нейронов в предыдущем слое для контроля выхода за пределы массива.

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

Матрица весовых коэффициентов в нашем буфере представляется в виде вектора, в котором последовательно идут N весовых коэффициентов первого нейрона, за ними — N весовых коэффициентов второго нейрона и так далее. N на 1 элемент больше количества нейронов в предыдущем слое, так как мы используем элемент bias-смещения. У нас уже есть порядковый номер нейрона и количество нейронов в предыдущем слое, поэтому мы можем определить смещение в векторе весовых коэффициентов до первого весового коэффициента нашего нейрона.

Далее создадим локальную переменную для сбора накопительной суммы произведений и инициируем ее значением коэффициента bias-смещения. После этого организовываем цикл и считаем сумму произведений исходного сигнала на весовые коэффициенты для конкретного нейрона нашего потока. Здесь следует обратить внимание, что в OpenCL существует поддержка векторных операций. Это особый вид операций, который позволяет микропроцессору за один цикл производить одну операцию сразу над несколькими значениями. В предложенной реализации я использовал векторы типа double4. Это вектор из четырех элементов типа double. При этом для перевода данных из массива дискретных значений в вектор была создана функция ToVect4, которую мы рассмотрим немного позже. Для получения суммы произведений я использовал встроенную функцию dot. Она относится к векторным операциям и позволяет получить дискретное значение произведения двух векторов. Это позволило нам использовать в цикле шаг 4 и тем самым в 4 раза сократить количество итераций.

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

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

Полный код кернела представлен ниже.

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

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

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

То же можно сказать и о реализации функции активации Swish.

Но есть некоторые сложности при реализации функции Softmax. Это связано со сложностью передачи данных между потоками для подсчета суммы всех значений вектора экспоненты в нейронном слое. Для решения вопроса было решено разделить процесс на несколько этапов. Вдобавок мы воспользуемся возможностью Work-group использовать общие массивы в локальной памяти.

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

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

Обратите внимание, что OpenCL не предоставляет возможности создания динамических массивов. Поэтому размер локального массива должен быть определен до компиляции программы. Значит, нам надо искать некий компромисс. Излишний размер ведет к неэффективному использованию памяти. А слишком малый размер массива ограничивает количество активных параллельных потоков. Конечно, решать такую задачу гораздо легче, когда известны параметры используемого устройства и архитектура модели. Следовательно, нам нужен механизм, позволяющий легко и быстро изменить данный параметр перед непосредственным компилированием программы. И для этого нет ничего лучше, как воспользоваться макроподстановкой. Точно также, как мы сделали для типа данных. В коде мы указываем константу LOCAL_SIZE, значение которой мы присваиваем в нашем файле констант defines.mqh.

Далее мы организовываем цикл, в котором каждый поток суммирует свою часть экспоненциальных значений. И записывает полученное значение в соответствующий элемент локального массива.

После завершения операций цикла мы устанавливаем barrier, который позволяет синхронизировать все потоки локальной группы. Данная операция приостанавливает выполнение каждого потока в ожидание завершения итераций цикла всех потоков локальной группы.

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

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

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

За прямым проходом идет обратный проход. Начинается он с определения градиентов ошибки на выходе нейронного слоя. Эту функцию выполняет кернел CalcOutputGradient. В параметрах кернел получает указатели на три вектора: первые два вектора целевых значений и результатов прямого прохода являются исходными данными для функции, а третий служит для записи результатов вычислений. Также в параметрах указывается используемая функция потерь. Код кернела полностью повторяет алгоритм ранее рассмотренного аналогичного метода, написанного средствами MQL5. В нем мы также видим разветвление алгоритма в зависимости от используемой функции потерь.

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

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

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

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

Значение производных S-образных функций, таких как сигмоида и гиперболический тангенс, легко посчитать через результат функции активации.

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

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

Как можно заметить, производная функции выражается через значение функции активации, значения сигмоиды и параметра функции активации β. Подставив формулу сигмоиды, получим нижеследующее выражение.

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

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

Наиболее сложным выглядит кернел вычисления производной функции Softmax. Как и в случае S-образных функций, для вычисления производной функции Softmax достаточно значений самой функции активации. Только, как вы помните, для вычисления значения одного элемента вектора при прямом проходе мы использовали значения всех элементов вектора до функции активации (вычисление общей суммы экспонент). Следовательно, значение каждого элемента после активации зависит от всех элементов вектора до активации. А значит, каждый элемент до активации должен получить свою долю ошибки с каждого элемента на выходе нейронного слоя. В общем виде производная функции Softmax вычисляется по формуле ниже.

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

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

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

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

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

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

Распределение градиента ошибки через нейронный слой будет осуществляться в кернеле CalcHiddenGradient. На вход кернелу в параметрах подаются указатели на 3 массива:

• weights — матрица весовых коэффициентов;
• gradients — массив градиентов, скорректированных на производную функции активации;
• gradient_inputs — массив для записи градиентов предшествующего слоя.

Кроме того, в параметрах указывается количество нейронов в верхнем слое (размер массива gradients). Алгоритм построения кернела сильно напоминает метод прямого прохода — мы так же используем функции dot и ToVect4. Отличие в используемых массивах: если при прямом проходе мы брали входной сигнал и умножали на весовые коэффициенты, то сейчас мы градиент ошибки умножаем на весовые коэффициенты. Есть еще один момент в использовании функции ToVect4 для матрицы весовых коэффициентов. Когда мы рассматривали эту функцию для прямого прохода, то говорили о схожей функции параметров первого элемента для копирования start и сдвига shift. Тогда мы использовали шаг в 1 элемент. Сейчас, перебирая элементы по массиву градиентов, мы будем подбирать соответствующие весовые коэффициенты. Но если при прямом проходе у нас нейроны и весовые коэффициенты шли по порядку, то при обратном проходе мы берем весовые коэффициенты поперек матрицы весов. В векторном выражении матрицы весов мы будем использовать шаг между двумя элементами для копирования на 1 элемент больше количества нейронов в предыдущем слое (элемент смещения bias). В то же время сдвиг будет равен порядковому номеру обрабатываемого нейрона нижнего слоя.

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

В остальном, мы также используем векторные операции, которые позволяют нам использовать цикл с шагом в 4 элемента.

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

1. Накопление градиентов ошибки на некотором интервале.
2. Усреднение накопленного градиента и корректировка матрицы весов.

Собирать градиенты ошибок для каждого весового коэффициента будем в кернеле CalcDeltaWeights. В параметрах кернелу передаются указатели на массив данных на выходе предыдущего слоя, массив градиентов и массив для накопления дельт необходимых корректировок весовых коэффициентов.

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

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

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

Напомню, что мы используем дополнительный элемент bias-смещения. Для этого элемента используется постоянное значение входящего элемента равное единице. Мы его не учли при создании пространства задач, но мы просто обязаны также накапливать градиент ошибки по нему. С математической точки зрения производная умножения на 1 равна единице. А значит, градиент ошибки для данного элемента равен градиенту ошибки перед функцией активации соответствующего нейрона. Чтобы не повторять итерацию записи градиента ошибки веса bias-смещения, данную итерацию мы будем осуществлять только для потока с индексом 0 во втором измерении.

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

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

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

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

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

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

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

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

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

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

В остальном алгоритм кернела повторяет ранее рассмотренные методы обновления весовых коэффициентов.

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

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

Рассмотрим реализацию данного метода в кернеле AdaDeltaUpdate. В параметрах кернелу передаются указатели на четыре массива данных:

• delta_weights — массив накопленных градиентов ошибки;
• weights — матрица весовых коэффициентов;
• momentumW — матрица экспоненциальных средних квадратов весовых коэффициентов;
• momentumG — матрица экспоненциальных средних квадратов градиентов ошибки.

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

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

Последним мы изучали метод адаптивной оценки моментов Adam. Напомню математические формулы этого метода.

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

Рассмотрим реализацию метода в кернеле AdamUpdate. В параметрах кернела, как всегда, мы будем передавать указатели на массивы данных:

• delta_weights — накопленные дельты градиентов;
• weights — матрица весовых коэффициентов;
• momentumM — матрица накопленных градиентов;
• momentumV — матрица накопленных квадратов градиентов.

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

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

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

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

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

И в заключение кернела сохраним полученные значения в соответствующие массивы. Не забываем обнулить массив накопленных дельт.

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