Решение уравнения Пуассона на массивнопараллельном графическом процессоре с использованием нейросетевой оптимизации Текст научной статьи по специальности «Компьютерные и информационные науки»

ИЗВЕСТИЯ

ПЕНЗЕНСКОГО ГОСУДАРСТВЕННОГО ПЕДАГОГИЧЕСКОГО УНИВЕРСИТЕТА имени В. Г. БЕЛИНСКОГО ФИЗИКО-МАТЕМАТИЧЕСКИЕ И ТЕХНИЧЕСКИЕ НАУКИ № 18 (22)2010

IZVESTIA

PENZENSKOGO GOSUDARSTVENNOGO PEDAGOGICHESKOGO UNIVERSITETA imeni V. G. BELINSKOGO PHYSICAL, MATHEMATICAL AND TECHNICAL SCIENCES № 18 (22) 2010

УДК 004.272:004.032.26

РЕШЕНИЕ УРАВНЕНИЯ ПУАССОНА НА МАССИВНО-ПАРАЛЛЕЛЬНОМ ГРАФИЧЕСКОМ ПРОЦЕССОРЕ С ИСПОЛЬЗОВАНИЕМ НЕЙРОСЕТЕВОЙ ОПТИМИЗАЦИИ

© Н. О. МАТВЕЕВА, Е. В. АРТЮХИНА Пензенский государственный педагогический университет им. В. Г. Белинского, кафедра вычислительных систем и моделирования e-mail: tigermno@rambler.ru

Матвеева Н.О., Е. В. Артюхина - Решение уравнения Пуассона на массивно-параллельном графическом процессоре с использованием нейросетевой оптимизации // Известия ПГПУ им. В. Г. Белинского. 2010. № 18 (22). С. 189-195. - Рассматривается возможность эффективной реализация алгоритма обучения радиально-базисной нейронной сети (RBFNN), решающей уравнение Пуассона, на процессоре с массивно-параллельной архитектурой. Предлагаются методы распараллеливания основных вычислений алгоритма обучения RBFNN с использованием технологии CUDA. Проводится сравнение времени обучения RBFNNна центральном и массивно-параллельном графическом процессорах.

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

Matveyeva N. O., Artyukhina E. V. - Solution of the differential Poisson equation on massive parallel Graphics Processing Unit using the neuronet optimization // Izv. Penz. gos. pedagog. univ. im. V. G. Belinskogo. 2010. № 18 (22). P. 189-195. - The article describes the possibility of the effective realization of the training algorithm of the radial basis neural network (RBFNN) on massive parallel Graphics Processing Unit (GPU) that solves the differential Poisson equation. The parallel realization of the training RBFNN algorithm are proposed, that uses CUDA. Execution time using CPU versus GPU is compared.

Keywords: massive parallel processing, Graphics Processing Unit, parallelism, CUDA, radial basis neural network, differential equation.

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

Вычислительные мощности современных суперкомпьютеров на данный момент доступны далеко не всем, но созданные относительно недавно и активно развивающиеся многопроцессорные системы, первоначально появившиеся как видеокарты новой архитектуры, предоставляют реальную возможность ускорения вычислений путем их распараллеливания. Графическим процессорам первоначально был свойственен параллелизм, поэтому неудивительно, что видеокарты стали многопроцессорными. Еще до появления графических процессоров (GPU) нового поколения были попытки использовать вычислительные ресурсы видео карты для решения задач общего назначения, а после выхода удобной технологии программирования новых GPU nVidia CUDA (Compute Unified Device Architecture) [1], они стали использоваться для решения самых разных задач, требующих больших вычислительных затрат.

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

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

Альтернативным подходом является использование различных вариантов метода взвешенных невязок, когда в качестве базисных функций применяются радиально-базисные функции (RBF - radial basis function) [2-6]. Использование RBF-функций рассматривается как бессеточный метод (meshless, meshfree). Применение радиально-базисных функций позволяет обойтись без трудоемкого построения сетки и позволяет получить приближенное дифференцируемое решение в произвольных точках области. Недостатком метода является сложность определения параметров радиально-базисных функций.

Бессеточные методы эффективно реализуются на радиально-базисных нейронных сетях (RBFNN) [6-9]. RBFNN отличаются простотой, так как содержат только один скрытый слой, поэтому исключается неформализуе-мый подбор структуры сети, необходимый для большинства типов нейронных сетей. RBFNN могут быть реализованы на параллельных системах. Главное достоинство RBFNN состоит в использовании принципов обучения для формирования оптимальных параметров радиально-базисных функций.

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

Структура радиально-базисной нейронной сети

Радиально-базисная нейронная сеть (рис. 1) представляет собой сеть с двумя слоями [6-8].

> ц

Рис. 1. Радиально-базисная нейронная сеть

Первый слой осуществляет преобразование входного вектора х с использованием радиально-базисных функций (RBF). Практически используются различные радиально-базисные функции. В дальнейшем будем использовать наиболее часто употребляемую функцию - Гауссиан, имеющий вид для k -го нейрона

ф* (х ) = ехР (-г2/а2) (1)

где х - входной вектор, гк = ||х - с* || - радиус к -го нейрона, ак - ширина к -го нейрона, с* - центр к -го нейрона.

Выход сети описывается выражением

т

и=Х ф* (х) > (2)

к =1

где wk - вес, связывающий выходной нейрон с к -ым нейроном первого слоя; т - число нейронов первого слоя.

Алгоритм обучения радиально-базисной нейронной сети

Рассмотрим градиентный алгоритм обучения радиально-базисной нейронной сети на примере решения двумерного уравнения Пуассона

+д7=/(х,у"), (х,у)е°, (3)

и = Р (х, У), (X, у)е50, (4)

где 30. - граница области; / и р - известные функции (х, у).

_ Г2

Выбирая в качестве радиально-базисной функции Гауссиан (1), определяемый как фк (х,у) = е ак, рассмотрим RBF-сеть как аппроксиматор функции решения:

Кх )=£ wk Фк (гк), (5)

к=1

где т - число радиально-базисных функций (скрытых нейронов), гк = ^(х _ схк)2 + (у _ сук ) , (схк сук ) - координаты центра нейрона к, ак - ширина нейрона к.

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

ных точках

( , 1 ^ д ли (х, у,) дли (х,., у,)

(V, с, а

1 N

I(м>,с,а) = - £

2 /=1

д2и (х, у,) д2и (х, у,) )

---------- +------/(х, у,)

+■2 £[и (х, у})-р1 ], (6)

21+1

дх2 ду2

где X - штрафной множитель, р1 - значение граничных условий первого рода в точке 1 границы, N и К - количество внутренних и граничных контрольных точек.

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

Очень важно соблюдать при обучении соотношение между оптимальным количеством нейронов т и количеством контрольных точек [7, С. 389]

т ж 3N + К , (7)

где ж - знак пропорциональности, N - число внутренних контрольных точек, К - число контрольных точек на границе.

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

Градиентный алгоритм спуска, одновременно оптимизирующий веса, центры и ширину может быть построен в виде последовательности двух шагов [9].

Шаг 1. Зафиксировав ^ки), находим центры и ширину, минимизирующие функционал ошибки (6), используя метод градиентного спуска

д1 (с(и_1) а'"_1) ^(и)) с(и) = с(и-1) _ в(и-1) У к , к , к !

Ск ск в дсГ» ,

2

д1 (0(п) а(и-1) “(и))

(п) (И-1) (п-1) V - ’ - ’ - —

а- = а; '-а' '— -------т—т------,

- - да[пе)

где — и 2^(Ди, - /-) “-е а

да- -=1

—4 - 3а2г2 + -4

а!

+ и 2 ХД. -р, а-20-.

^о-игад.,-С)““ Є •(',-о-)• + і.2^(«;-л)е- •(х^.

д0Х- ,= 1 а- 1=1 а-

лт ( Г'2 2 л 2 К - — (х о ^

-О- и2]Т(Д[(-/,)“-е а іг,-о,.У^ + Ь2]Г(.1 -р,)• е а2

,-и1 (Л,

у- -и1 и- ,=1 М-

Направление изменения параметров сети выбираются противоположным градиенту функционала I.

Шаг 2. Зафиксировав центрыи ширину, находим веса, минимизирующие функционал ошибки (6).

Для обучения весов сети использовался алгоритм сопряженных градиентов для минимизации квадратичного функционала, предложенный в [6]. Так как RBF-сеть линейна относительно весов, то метод сопряженных градиентов позволяет найти глобальный минимум функционала ошибки. Данный алгоритм эквивалентен при зафиксированных центрах и ширине решению системы линейных алгебраических уравнений Aw = s, где А =1 (М ТМ + ) - симметричная положительно определенная матрица, s =1 (М^ + Я^Тр ), М - матри-

-І ..2 „2

ца N х т с элементами т,к = 4е ак Г‘к а , N - матрица К х т с элементами пк = е ак , f - вектор значений функ-

аь

ции правой части во внутренних точках, р - вектор граничных условий в граничных контрольных точках. Метод сопряженных градиентов для обучения сети представлен на рисунке 2.

го = s - А« 0

Ро = го

Для 1 = 1,2,...

а 1 = О;, Г )/(Ар,, р;)

« М = « ; + а ;Р ;

1 = Г -а; АР;

Р; = 1 Г1+1)/(Г1, Г)

Если||г,+^|/||^| <6, то КОНЕЦ

Р1+1 = Г+1 +Р ;Р 1 увеличить1

Рис. 2. Метод сопряженных градиентов для обучения сети.

Алгоритмы распараллеливания вычислений на графическом процессоре в технологии CUDA

Как было описано выше, для минимизации весов сети использовался алгоритм сопряженных градиентов. Распараллеливание данного алгоритма заключается в основном в параллельной реализации матрично-векторных операций как описано в [10]. Но кроме коррекции весов большая доля времени приходится на вычисление матриц М и N, расстояний гк, коррекцию векторов ширины и центров. При этом элементы данных матриц и векторов находятся не с помощью матрично-векторных операций, а в циклах, что делает невозможным использование стандартной библиотеки CUDA CUBLAS. Например, функция вычисления матрицы М и г1 в МА^АВ представлена на рис. 3.

function rim global w a nv m r2v f R1 M R1=zeros(1,nv);

M=zeros(nv,m); for i=1:nv

for k=1:m

r2=r2_v(i,k);

M(i,k)=4*exp(-r2/(a(k)A2))*((r2-

a(k)A2)/(a(k)A4));

R1(i)=R1(i)+M(i,k)*w(k);

end

Ri(i)=Ri(i)-f(i);

end

end

Рис. 3. Функция вычисления матрицы М и вектора R1 в MATLAB

global void r1m(float* M, float* R1,float*

r2 v,float* a, float* w, float* f, int n)

{

const unsigned int tx = threadldx.x;

const unsigned int by = blockldx.y;

int i = tx; int j = by*n + i;

extern shared float sdata[];

float r = r2 v[j]; float ak2 = powf(a[i],2);

float sd = 4*expf(-r/ak2)*(r-ak2)/powf(ak2,2);

M[j]= sd;

sdata[i] = sd*w[i];

for(int stride = BLOCK SIZE/2; stride > 0; stride >>=

1)

{

syncthreads(); if(tx < stride)

sdata[tx] += sdata[tx + stride];

}

syncthreads();

if(tx == 0)

R1[by] = sdata[0] - f[by];

}

Рис. 4. Функция вычисления матрицы М и вектора R1 в CUDA

Для вычисления одного элемента матрицы М требуется прочитать из памяти гл и ак и выполнить все опе-—2

-—- 2 2 2 — — а

рации по формуле т- = 4е а- -- - , т.е. необходимо небольшое количество чтений из памяти и большое коли' ак

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

На графическом процессоре каждый элемент матрицы М вычисляется независимо от других в отдельном

т

потоке. Вычисление вектора г, = ^ т л“к - /, реализуется в этом том же ядре, что и матрицы М, для уменьшения

-=1

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

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

Анализ результатов экспериментов

Алгоритм обучения радиально-базисной нейронной сети реализовывался с помощью технологии CUDA [1, 11] на графическом процессоре nVidia GeForce 8800 GTX, на центральном процессоре Intel Pentium 4 3 GHz - в системе MATLAB 7.7.0. Экспериментальное исследование проводилось на примере модельной задачи (3) - (4), для f (х, y) = sin(nx)• sin (ny), p (x, y) = 0. Анализ результатов экспериментов показал, что часть алгоритма, не включающая обучение весов сети, распараллеливается очень эффективно, значительно уменьшая время вычислений. Метод сопряженных градиентов для обучения весов становится более эффективным при увеличении числа нейронов, что можно видеть на рисунке 5. Для решаемой задачи увеличивать количество нейронов не имело смысла, поэтому ускорения представленные в таблице 1 были получены в основном за счет распараллеливания шага 1, вычисления матриц M, N и др. При решении более сложного уравнения эффективность распараллеливания на графическом процессоре будет значительней, но и для решаемой задачи было получено существенное ускорение.

18 |---------Т-----------Т-----------Т----------1-----------

16 -.........;...........<..........:• i..........-

14 -.........:•..........•:........../..........;..........-

12 [■........г............f /........;..........|..........-

(И т ■

■XI

а.

о

:£ о

В -

4 -

2 -

о ■ о

Рис. 5. График зависимости величины ускорения от числа нейронов для метода сопряженных градиентов

Таблица 1.

Результаты экспериментов

Параметры сети GPU, мс CPU, мс T Ускорение S = -CPU- TGPU

Число контр. точек = 224 , нейронов = 64 1,67 64,59 38,58

Число контр. точек = 524, нейронов = 64 2,11 176,3 83,42

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

Заключение

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

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

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

Благодарности. Работа выполнена по тематическому плану научно-исследовательских работ Пензенского государственного педагогического университета, проводимых по заданию Федерального агентства по образованию.

СПИСОК ЛИТЕРАТУРЫ

1. NVIDIA CUDA Homepage [http://www.nvidia.ru/object/cuda.html].

2. Buhmann M. D. Radial Basis Functions: Theory and Implementations. Cambridge: Cambridge University Press, 2004. 259 pp.

3. Liu G. R., Gu Y. T. An Introduction to Meshfree Methods and Their Programming. NY: Springer, 2005. 479 pp.

4. Meshfree Methods for Partial Differential Equations / Ed. M. Griebel, M.A. Schweitzer. NY: Springer, 2008. 412 pp.

5. Толстых А. И., Широбоков Д. А. Бессеточный метод на основе радиальных базисных функций // Журнал вычислительной математики и математической физики. 2005. Т. 45. № 8. С. 1498-1505.

6. Горбаченко В. И., Артюхина Е. В. Исследование градиентных алгоритмов обучения весов радиально-базисных нейронных сетей для решения краевых задач математической физики // Известия ПГПУ им. Белинского. 2009. №13(17). С. 134-143.

7. Хайкин С. Нейронные сети: полный курс. М.: Вильямс, 2006. 1104 с.

8. Осовский С. Нейронные сети для обработки информации. М.: Финансы и статистика, 2002. 344 с.

9. Jianyu L., Siwei L., Yingjiana Q., Yapinga H.. Numerical solution of elliptic partial differential equation using radial basis function neural networks // Neural Networks. 2003. Vol. 16. № (5/6). P. 729-734.

10. Горбаченко В.И., Матвеева Н.О. Реализация итерационных алгоритмов решения систем линейных алгебраических уравнений на графических процессорах в технологии CUDA // Вопросы радиоэлектроники. Серия: электронная вычислительная техника. Вып. 5. М., 2008. С. 65-75.

11. NVIDIA CUDA Compute Unified Device Architecture. Programming Guide. Version 2.0. 06.07.2008.