Алгоритм увеличения точности нейронных сетей и его приложения
 Скачать 3.63 Mb.
|
|
что каждый столбец 𝐵 будет считан из GPU RAM как минимум |𝑖| раз, и каждая стро- ка 𝐴 считана |𝑗| раз. Это значительно уменьшает соотношение между арифметическими операциями и доступами к памяти. Проблема решается разбиением 𝐶 на подматрицы, каж- дая из которых хранится в общей памяти блоков, скорость которой значительно превыша- ет скорость GPU RAM. Существует большое количество оптимизаций [127], применимость которых зависит как от размера матриц, так и от конкретной версии платформы CUDA. Вместо самостоятельного отслеживания изменений и модификации кода, практически все сопутствующие улучшения производительности можно получить используя стороннюю биб- лиотеку CUBLAS [128] (Basic Linear Algebra Subroutine), а конкретно, функцию cublasSgemm(cublasHandle handle, cublasOperation_t transa, cublasOperation_t transb, int m,int n,int k, const float* alpha, const float* A,int lda, const float* B,int ldb, const float* beta, float* C, int ldc), которая выполняет операцию C = aplha * op (A) × op (B) + beta * C. (4.4.1) 98 имя назначение handle определяет поток, выполняющий умножение transa флаг транспонирования A transb флаг транспонирования B alpha, beta константы в выражении 4.4.1 m, n, k число столбцов C , строк C и строк op (B) A, B, C адреса памяти для матриц lda, ldb, lcd чисто строк в A , B и C Таблица 4. Аргументы cublasSgemm. Так, нередко после обновления библиотеки, производительность увеличивалась на 10 ∼ 15% Минусом подхода конечно же является потеря гибкости в силу отсутствия исходного кода. Назначения аргументов представлены в Таблице 4. Подразумевается хранение массивов по столбцам. Эффективность cublasSgemm сильно зависит от размера матриц. Продемонстри- руем её для NVIDIA RTX 2060, CUDA 10.2. Устройство имеет частоту 1275Mhz и теорети- ческую производительность в 4.9 · 10 12 операций с плавающей точкой в секунду (FLOPS – floating point operations per second). Рассмотрим матричное умножение, которое требуется для прямого распространения |𝑎| паттернов по сети между слоями со 128 нейронами: 𝐶 128×|𝑎| = 𝐴 128×128 × 𝐵 128×|𝑎| Результаты приведены на Рисунке 4.4.1. Как можно видеть, перемножение небольших мат- риц крайне неэффективно 4 , однако, при увеличении |𝑎| эффективность вычислений заметно возрастает. Например, если разбить матрицу 𝐵 из 4096 векторов на блоки по 64 и перемно- жить их с 𝐴 отдельно, то итоговое время будет в ∼ 10 раз больше, чем при их одновременной обработке. Разумеется, смысла в таком разбиении нет, однако в следующем пункте мы ука- жем способ «склеивания» матриц для разных производных, приводящий к схожему увели- чению |𝑎| . Производительность в рассмотренном случае даже для большого числа паттернов не дотягивает до теоретической из-за относительно малого количества строк, определяемого числом нейронов в скрытых слоях. Так, при увеличении их до 512, уже для 6144 входных векторов можно получить 4.5 · 10 9 FLOPS, то есть 92% от теоретически возможной. 4.4.3. Особенности реализации. Основная идея реализации алгоритма довольно про- ста – вычисления должны быть устроены так, чтобы задержки между вызовами cublasSgemm и поэлементными ядрами были минимальны, а все массивы, на которых они работают, име- ли максимальный размер. Так как количество строк в матрицах определяется архитекту- рой сети, можно лишь попытаться увеличить число столбцов. Очевидно, что все паттерны лучше всего обрабатывать одновременно. Число столбцов также можно увеличить с помо- щью склейки матриц различных производных. В формуле (4.1.5) матричное произведение 𝑊 𝜃𝜅 × [𝜕 𝑠 𝜎 (𝑧 𝜅𝑎 )] вычисляется для всех 𝑠 независимо. Пусть его нужно выполнить 𝑛 раз. 4 по мере развития CUDA и CUBLAS эта ситуация исправляется 99 64 128 256 512 1448 4096 11 585 32 768 0.2 1 2 3 4 5 a GFlops Рис. 4.4.1. Производительность матричного умножения на RTX 2060 в зави- симости от количества векторов, обрабатываемых одновременно (масштаб по оси абсцисс логарифмический). Прямая сверху – теоретическая производитель- ность на данной частоте. Если вместо этого сначала вычислить 𝜕 𝑠 𝜎 (𝑧 𝜅𝑎 ) для всех 𝑠 , а результаты записать в одну матрицу с новым индексом столбца ˜ 𝑎 ∈ [1, . . . , |𝑎|, . . . , 𝑛 · |𝑎|] , то вместо 𝑛 умножений мат- рицы |𝜃| × |𝜅| на матрицу |𝜅| × |𝑎| можно всего один раз умножить её на |𝜅| × 𝑛 · |𝑎| Такой «трюк» вынуждает хранить матрицы только по столбцам, однако, в случае большого количества вычисляемых производных, он позволяет значительно повысить эффективность вычислений даже для малого количества паттернов. Формула обратного распространения (4.2.2) допускает аналогичную модификацию. Что касается задержки между вызовами функций, возможны два варианта: либо па- мять для всех матриц активностей и их производных выделена на GPU, либо существуют промежуточные операции загрузки/сохранения, которые выполняются во время других вы- числений. Первый вариант подходит для GPU с достаточным объемом памяти. Второй ва- риант может быть применён при подходящей скорости обмена данными между GPU RAM и CPU RAM. Как правило, у персональных компьютеров она велика, тогда как в облачных сервисах её недостаточно. Опишем метод уменьшения количества промежуточной памяти, необходимой для вы- полнения (4.2.2). Вместо предварительного вычисления произведений [︀𝑊 𝜃𝜅 ]︀ 𝑇 ×𝜕𝐸/𝜕 (︀𝜕 𝑠 𝑧 𝜃𝑎 )︀ и их хранения в промежуточной памяти, их можно записывать непосредственно в память, выделенную для 𝜕𝐸/𝜕 (𝜕 𝑝 𝑧 𝜅𝑎 ) при 𝑝 = 𝑠 , а затем вычислять 𝜕𝐸/𝜕 (𝜕 𝑝 𝑧 𝜅𝑎 ) в правильном порядке. Если при таком хранении начать вычисление для некоторого 𝑝 , то матричное про- изведение, записанное по этому адресу, должно быть поэлементно умножено на 𝜎 ′ (𝑧 𝜅𝑎 ) , следовательно, его нельзя будет использовать где-либо ещё. Чтобы обойти это ограниче- ние, нужно начать с 𝑝 = (0, 0, . . . , 0) , что изменит матрицу [︀𝑊 𝜃𝜅 ]︀ 𝑇 × 𝜕𝐸/𝜕 (︀𝜕 𝑠 𝑧 𝜃𝑎 )︀ с 𝑠 = (0, 0, . . . , 0) , которая, однако, не используется ни для каких других 𝑝 , потому что 𝑠−𝑝 долж- но быть положительным. Следом можно вычислить любую первую производную, например 𝑝 = (1, 0, 0, . . . , 0) , что изменит матрицу [︀𝑊 𝜃𝜅 ]︀ 𝑇 × 𝜕𝐸/𝜕 (︀𝜕 𝑠 𝑧 𝜃𝑎 )︀ для 𝑠 = (1, 0, 0, . . . , 0) , 100 однако, единственный другой случай, где (1, 0, 0 . . . , 0) − 𝑝 не имеет отрицательных компо- нент есть 𝑝 = (0, 0, . . . , 0) , но он уже был вычислен. Таким образом могут быть вычислены все первые, затем все вторые производные, и так далее. 4.4.4. cpde-2. С учётом вышеизложенных оптимизаций был написан программный ком- плекс «cpde-2». Опишем начальные данные и дополнительный код, с помощью которых мож- но воспользоваться «cpde-2» для обучения нейронных сетей. Для целей этого параграфа ну- мерация для всех индексов будет начинаться с нуля, то есть если ранее 𝑎 ∈ [1, . . . , |𝑎|] то сейчас тот же индекс 𝑎 ∈ [0, . . . , |𝑎| − 1] . Итак, для дифференциального уравнения или си- стемы уравнений существует список встречающихся производных. Мы начнем с введения на множестве мульти-индексов алфавитного порядка, с помощью которого любой такой список запишется единственным образом. Два примера приведены в Таблице 5. 0 0 1 0 0 1 2 0 1 1 0 2 3 0 2 1 1 2 0 3 0 0 0 1 0 0 0 1 0 0 0 1 2 0 0 1 1 0 1 0 1 0 2 0 0 1 1 0 0 2 Таблица 5. Мульти-индексы всех производных до 3 порядка для 2 переменных и до 2 порядка для 3 переменных, записанные в алфавитном порядке. Согласно Лемме 2, если в уравнении встречается некоторая производная 𝜕 𝑠 , то в список необходимо включить все другие 𝜕 𝑝 и 𝜕 𝑞 для которых 𝜕 𝑠 = 𝜕 𝑝 𝜕 𝑞 , например 𝜕 3 /𝜕𝑥 1 𝜕𝑥 2 2 потребует вычисления 𝜕/𝜕𝑥 1 , 𝜕/𝜕𝑦 2 , 𝜕 2 /𝜕𝑥 2 2 , 𝜕 2 /𝜕𝑥 1 𝜕𝑥 2 , но не 𝜕 2 /𝜕𝑥 2 1 . Таким образом, в списках производных могут отсутствовать некоторые строки, присутствующие в полном списке производных. Алфавитный порядок аргументов был также выбран для аргументов функций вида (4.1.8). Пусть имеется дифференциальное уравнение, в котором встречаются k производных включая нулевую 𝑉 (𝑥 𝛼𝑎 , 𝜕 𝑠[0] 𝑣, 𝜕 𝑠[1] 𝑣, 𝜕 𝑠[2] 𝑣, . . . , 𝜕 𝑠[k−1] 𝑣) = 0, (4.4.2) где мульти-индексами 𝑠 [1] , . . . , 𝑠 [k − 1] перечислены все используемые производные в ал- фавитном порядке (нулевой мульти-индекс 𝑣 ≡ 𝜕 𝑠[0] 𝑣 введён для удобства). Для решения (4.4.2) с помощью «cpde-2» необходимо задать архитектуру сети, список производных, мас- сивы 𝜕 𝑠[𝑖] 𝑥 𝛼𝑎 , а также написать ядро для вычисления невязки. Рассмотрим все эти пункты по порядку. 4.4.4.1. Параметры задачи. Архитектура сети и список производных указываются в тек- стовом файле SconfS.txt, который расположен в одной директории с исполняемым файлом 101 8 2 64 64 64 64 64 64 1 1 0 0 0 0 0 0 1 28 3 0 0 0 1 1 0 1 0 1 2 2 0 2 1 1 2 0 2 3 3 0 3 2 1 3 1 2 3 0 3 4 4 0 4 3 1 4 2 2 4 1 3 4 0 4 5 5 0 5 4 1 5 3 2 5 2 3 5 1 4 5 0 5 Таблица 6. Файл SconfS.txt, использованный для решения уравнения Пуас- сона методом обучения с исключением. «cpde-2». Пример SconfS.txt для краевой задачи уравнения Пуассона, разобранной в Па- раграфе 2.3, представлен в Таблице 6. Самое первое число в файле есть количество слоёв сети, далее следуют количества нейронов в каждом слое. Так как рассматривается скалярная функция от двух переменных, то входной слой имеет размерность 2, а выходной размерность 1. В следующей строке перечислены типы функций активации каждого слоя. Ноль «0» со- ответствует логистической сигмоиде 𝜎 (𝑥) = 1/ (1 + exp (−𝑥)) , а единица «1» – линейной функции активации 𝜎 (𝑥) = 𝑥 . Далее идёт количество строк и столбцов в таблице производ- ных. Число столбцов на единицу больше, чем количество независимых переменных, так как в каждой строке перед мульти-индексом дополнительно указан порядок данной производной. Тип разделителей между числами роли не играет. Несмотря на то, что в уравнении Пуас- сона встречаются производные не выше второй, из-за дополнительного дифференцирования невязки возникают производные вплоть до пятого порядка. 4.4.4.2. Начальные данные. Матрицы начальных данных 𝜕 𝑠[𝑖] 𝑥 𝛼𝑎 , где 𝑖 ∈ [0, . . . , k − 1] , записываются в файл in.dat в двоичном виде. Его структура представлена в Таблице 7. Первое число есть общее количество матриц k , далее следует текстовый разделитель «#» и два целых числа, задающие размер первой матрицы. Далее идет сама матрица 𝜕 𝑠[0] 𝑥 𝛼𝑎 = 102 Переменная Хранимые данные Тип переменной k количество матриц Integer32 # разделитель Character8 h 0 число строк 𝑥 𝛼𝑎 Integer32 w 0 число столбцов 𝑥 𝛼𝑎 Integer32 𝜕 𝑠[0] 𝑥 𝛼𝑎 матрица 𝑥 𝛼𝑎 Float32, h 0 * w 0 # разделитель Character8 h 1 число строк 𝜕 𝑠[1] 𝑥 𝛼𝑎 Integer32 w 1 число столбцов 𝜕 𝑠[1] 𝑥 𝛼𝑎 Integer32 𝜕 𝑠[1] 𝑥 𝛼𝑎 матрица 𝜕 𝑠[1] 𝑥 𝛼𝑎 Float32, h 1 * w 1 𝜕 𝑠[𝑘−1] 𝑥 𝛼𝑎 матрица 𝜕 𝑠[k−1] 𝑥 𝛼𝑎 Float32, h k−1 * w k−1 Таблица 7. Структура файла in.dat 𝑥 𝛼𝑎 записанная линейно по столбцам в формате float32. Затем идет следующий текстовый разделитель «#» и так далее. 4.4.4.3. Вычисление невязки. Обсудим устройство ядра, вычисляющего невязку. Оно но- сит название e_universal, находится в файле initd.txt, и имеет следующую структуру __global__ void e_universal ( f l o a t ** ADDR, f l o a t ** ADDR_dE, f l o a t ** ADDRX, f l o a t **ADDRe, i n t h , i n t xh , i n t w){ long unsigned i n t t i d = threadIdx . x + blockIdx . x * blockDim . x ; // объявление остальных переменных while ( t i d < w){ //вычисления для фиксированного a = t i d t i d += blockDim . x * dridDim . x ; } } Мы будем считать невязку локальной величиной, а значит треды могут работают с разными 𝑎 параллельно. Конкретный номер столбца 𝑎 определяется переменной tid (thread id), цикл while требуется при несовпадающем числе тредов и столбцов (то есть почти всегда). Доступ к производным входа осуществляется через массив указателей float** ADDRX 𝜕 𝑠[𝑖] 𝑥 𝛼 [𝑎] → ADDRX[i][a * hX + 𝛼] = ADDRX[i][tid * hX + 𝛼]. 103 Здесь через hX обозначено число строк матрицы 𝑥 𝛼𝑎 (размерность входного пространства). В рассматриваемом примере hX=2. Аналогично предоставляется доступ к производным вы- хода 𝜕 𝑠[𝑖] 𝑣 𝜔𝑎 → ADDR[i][a * hV + 𝜔] = ADDR[i][tid * hV + 𝜔], здесь hV есть число выходов сети (количество неизвестных функций). В рассматриваемом примере hV=1. Индексы 𝛼 и 𝜔 перебираются внутри кода для каждого треда. При xH>1 или hV>1 доступ к памяти, вообще говоря, будет осуществляется неоптимальным образом, однако, в силу небольшого числа строк, это не влияет хоть сколько-нибудь заметным обра- зом на производительность алгоритма в целом. Фактически, в массивах ADDR и ADDRX находятся все необходимые данные для вычисления невязки 𝑉 в каждой точке сетки. Об автоматической генерации кода для её вычисления написано в Пункте 4.4.5. После вычисления невязок нужно сложить их в глобальную ошибку 𝐸 = ∑︀ 𝑎 𝑉 2 и приступить к обратному проходу, который вычислит градиент 𝐸 по весам сети и позволит её минимизировать. Как было установлено в Главе 2, результаты можно значительно улучшить если начать минимизировать ошибку с дополнительными производными 𝑉 , например до второго порядка 𝐸 = ∑︁ 𝑎 [︀𝑉 2 + 𝑉 2 𝑥 + 𝑉 2 𝑦 + 𝑉 2 𝑥𝑥 + 𝑉 2 𝑥𝑦 + 𝑉 2 𝑦𝑦 ]︀ , и в процессе постепенно отбрасывать старшие производные. В результате такого дополнения, список производных заметно вырастет. Согласно процедуре из Параграфа 4.2, обратный проход начинается с вычисления градиента 𝐸 по всем матрицам полученным на выходе. Мы предполагаем, что ядро сперва вычисляет невязку и её производные, для которых выделены локальные переменные с именем « e » и индексом согласно алфавитному порядку, в данном случае e0 = 𝑉 [𝑎], e1 = 𝜕𝑉 [𝑎] 𝜕𝑥 , e2 = 𝜕𝑉 [𝑎] 𝜕𝑦 , e3 = 𝜕 2 𝑉 [𝑎] 𝜕𝑥 2 , e4 = 𝜕 2 𝑉 [𝑎] 𝜕𝑥𝜕𝑦 , e5 = 𝜕 2 𝑉 [𝑎] 𝜕𝑦 2 Так как для каждого треда tid = a , то индекс 𝑎 у переменных e можно не писать. Величина для минимизации есть сумма их квадратов 5 в каждой точке 𝐸 = 1 2 ∑︁ 𝑎 [︀e0 2 + e1 2 + e2 2 + e3 2 + e4 2 + e5 2 ]︀ . (4.4.3) От 𝐸 необходимо взять градиент по всем элементам матриц ADDR (напомним, что мы рас- сматриваем случай hV=1). В силу линейности (4.4.3), при дифференцировании по каждому 𝑎 = tid сумма исчезает и остаются только локальные слагаемые 𝜕𝐸 𝜕 ADDR[0][tid] = 𝜕𝐸 𝜕𝑣 = 𝑒0 𝜕𝑒0 𝜕v + 𝑒1 𝜕𝑒1 𝜕v + 𝑒2 𝜕𝑒2 𝜕v + 𝑒3 𝜕𝑒3 𝜕v + 𝑒4 𝜕𝑒4 𝜕v + 𝑒5 𝜕𝑒5 𝜕v , 5 нормировка будет рассмотрена в следующем пункте 104 𝜕𝐸 𝜕 ADDR[1][tid] = 𝜕E 𝜕 (𝜕v /𝜕x 1 ) = e0 𝜕e0 𝜕 (𝜕v /𝜕x 1 ) + e1 𝜕e1 𝜕 (𝜕v /𝜕x 1 ) + e2 𝜕e2 𝜕 (𝜕v /𝜕x 1 ) + +e3 𝜕e 3 𝜕 (𝜕v /𝜕x 1 ) + e4 𝜕e4 𝜕 (𝜕v /𝜕x 1 ) + e5 𝜕e5 𝜕 (𝜕v /𝜕x 1 ) В общем случае 𝜕𝐸 𝜕 ADDR[i][tid] = 𝜕𝐸 𝜕 (︀𝜕 𝑠[𝑖] 𝑣 )︀ = e0 𝜕e0 𝜕 (︀𝜕 s[i ] v )︀ + e1 𝜕e1 𝜕 (︀𝜕 s[i ] v )︀ + e2 𝜕e2 𝜕 (︀𝜕 s[i ] v )︀ + e3 𝜕e3 𝜕 (︀𝜕 s[i ] v )︀ + e4 𝜕e4 𝜕 (︀𝜕 s[i ] v )︀ + e5 𝜕e5 𝜕 (︀𝜕 s[i ] v )︀ . Удобно воспользоваться промежуточным векторами De_DZ (вместо «_» будет стоять поряд- ковый номер вектора) с размерностью k . Они также являются локальными переменными для каждого треда, и индекс 𝑎 можно не писать (︃ 𝜕e0 𝜕𝑣 , 𝜕e0 𝜕 (︀𝜕 𝑠[1] 𝑣 )︀ , 𝜕e0 𝜕 (︀𝜕 𝑠[2] 𝑣 )︀ , . . . )︃ → float De0DZ[k], (4.4.4) (︃ 𝜕e1 𝜕𝑣 , 𝜕e1 𝜕 (︀𝜕 𝑠[1] 𝑣 )︀ , 𝜕e1 𝜕 (︀𝜕 𝑠[2] 𝑣 )︀ , . . . )︃ → float De1DZ[k], (︃ 𝜕e2 𝜕𝑣 , 𝜕e2 𝜕 (︀𝜕 𝑠[1] 𝑣 )︀ , 𝜕e2 𝜕 (︀𝜕 𝑠[2] 𝑣 )︀ , . . . )︃ → float De2DZ[k]. Теперь выражения для градиента 𝜕𝐸/𝜕 ADDR выражаются в простом виде. Они должны быть записаны в массив ADDR_dE: 𝜕𝐸 𝜕𝑢 = e0 · De0DZ[0] + e1 · De1DZ[0] + e2 · De2DZ[0] + . . . → ADDR _ dE[0][tid], 𝜕𝐸 𝜕 (𝜕𝑣/𝜕𝑥 1 ) ≡ 𝜕𝐸 𝜕 (︀𝜕 𝑠[1] 𝑣 )︀ = e0 · De0DZ[1] + e1 * De1DZ[1]+ + e2 * De2DZ[1] + . . . → ADDR _ dE[1][tid], 𝜕𝐸 𝜕 (︀𝜕 𝑠[𝑖] 𝑣 )︀ = e0 * De0DZ[i] + e1 * De1DZ[i]+ +e2 * De2DZ[i] + . . . → ADDR _ dE[i][tid]. (4.4.5) После этого комплекс «cpde-2» может осуществить обратный проход и вычислить гради- ент 𝐸 по всем весам сети. При высоких порядках обучения с исключением, суммы могут состоять из довольно большого числа слагаемых. Для уменьшения количества кода стоит 105 0 0 1 0 0 1 2 0 1 1 Таблица 8. Алфавитный список производных для выражений (4.4.6), (4.4.7). заранее определить какие из компонент векторов типа (4.4.4) будут нулями и заранее убрать соответствующие члены из сумм вида (4.4.5). 4.4.4.4. Нормировка. Очень часто при решении возникает ситуация, когда функция и её производные имеют разный порядок величины, то же самое может случиться и с невязкой. В таком случае нельзя будет минимизировать выражение (4.4.3). Универсальную нормиров- ку введём с помощью следующих соображений – итоговый вклад в градиент 𝐸 от каждой производной невязки полностью описывается вектором De_DZ длины |