Решение задачи поиска наименьшего расстояния до политопа с использованием графических ускорителей Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 16, № 5, 2011

Решение задачи поиска наименьшего расстояния до политопа с использованием графических ускорителей

Е. А. Нурминский, П.Л. Поздняк Институт автоматики и процессов управления ДВО РАН, Владивосток, Россия

e-mail: nurmi@dvo . ru, pozpl@dvo . ru

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

Ключевые слова: проекция, графические ускорители.

Введение

Решение задачи проекции начала координат на n-мерный выпуклый многогранник X, заданный своими вершинами, широко используется при построении алгоритмов негладкой оптимизации [1], распознавании образов и автоматической классификации [2], в компьютерной томографии и радиационной терапии [3]. Для практических приложений характерны большие размерности исходного пространства и большое количество X

требования по вычислительной эффективности и точности. Эти обстоятельства мотивировали поиск новых алгоритмов решения задачи проекции и способов ее реализации с использованием нестандартных вычислительных архитектур. В работе исследована возможность ускорения вычислений с помощью графических ускорителей GPU (Graphics Processing Unit). Применение GPU для решения задач большой размерности хорошо зарекомендовало себя в таких областях как молекулярная динамика [4], магнитно-резонансная томография [5] и т. д. В нашей задаче использование GPU при решении задач высокой размерности позволяет увеличить скорость вычислений до 1011 раз.

1. Метод подходящих аффинных подпространств

X

n-мерном евклидовом пространстве E со скалярным произведением xy и норм ой ||x|| = у/хх. Другими словами, речь идет о нахождении точки х* g X такой, что

x

* II2

mm x

x€X

2

(1)

В данной задаче предполагается, что множество X является политопом, т. е. выпуклой оболочкой конечного набора точек X пространства E:

X = {X\x2,...,Xm} = [хг,г е M},

где M = {1, 2,..., m} — индексное множество вершин политопа X, Для множества Z С X с индексным множеством ind (Z) С M аффинную (äff (Z)) и выпуклую (со (Z)) оболочки определим традиционным образом:

aff(Z) = {х = Л^хг, Л = 1}, (2)

xsz,ism d (Z) ¿sind (Z)

eo(X) = {x = ЛiXi,

x&z,i&md (z)

Y^ л, = 1,Л > 0,г е ind (Z)}. ieind (z)

X

набора точек SI = {xi,i е I С M} таких, что

min ||z||2 = min \\z\\2.

Подходящим нуль-мерным симплексом является, например, любая точка xi, г е M.

Метод подходящих аффинных подпространств состоит в построении последовательности множеств Ik С M таких, что aff(xi, г е Ik) представляет собой подходящее аффинное подпространство и расстояние от нуля

rk = min ||х||2

aff{xi,ieifc}

до этих подпространств строго монотонно убывает. Множества Ik и соответетвенно Xk = co{xi,i е Ik} будем называть базисом,

В методе подходящих аффинных подпространств на вход каждой итерации подается базис, построенный на предыдущих шагах, В результате выполнения итерации алгоритма строится новый базис со строго меньшим расстоянием до начала координат, чем предыдущий. Сама итерация содержит внутренний цикл, в котором, возможно, происходит удаление части вершин из базиса. Для размерности пространства n трудоемкость итерации внутреннего цикла составляет не более O(n2) при добавлении вектора в базис и не более 0(кп2) при удалении вершин из базиса, где к — количество векторов в получившемся базисе.

Описание алгоритма выглядит следующим образом. Положим, что в начале k-й итерации имеются базис Ik, соответствующий ему подходящий подеимплеке Sk и аффинное подпространство Hk = äff {Sk}, Дальнейшие действия описывают выполнение k

Внешний цикл алгоритма:

1) решаем задачу нахождения элемента минимальной нормы для аффинного подпространства Hk

min ||z||2 = ||zk||2; (4)

zSHk

2) если хггк > ||гк||2 для всех г 6 М \ 1к, то гкгг > ||гк ||2 для всех г 6 М и гк — решение задачи (1), и работа алгоритма на этом прекращается;

3) при невыполнении п. 2 выбираем хгн из условия

||хН гк ||2 = шт ||хггк ||2,

хг хк <\\хк\\2

где гн 6 М\1. Образуем новое индексное множество 1а = 1к и {гн}, полагаем и>а = гк и,

занулив счетчик внутренних итераций в = 0, выполняем следующий внутренний цикл:

a) образуем аффинное подпространство Щ = {хг,г 6 1а};

b) решаем задачу проекции на аффинное подпространство Щ

и 11 2 и Н 11 2

шт ||г|| = ||г || ;

г&На

c) если г 5 6 X, то полагаем Нк+1 = Н£ ,1к+1 = 1а и выходим из внутреннего цикла. Иначе полагаем

и^ = + (1 —

где находим максимальное ß такое, что uß G X, т. е.

ßmax = max ß.

Легко видеть, что найденное таким образом ßmax принадлежит некоторой грани симплекса Sk, т. е.

ußmax = ws+1 = Y^ eiX, Y^ 9i = 1 9i > 0' i G Is, H С Is.

ieis ieis

Положим, что я — индекс вектора в индексном пространстве I s, соответствующий выбранному ßmax, Удалим этот вектор из базиса, В результате получим некоторый подходящий подеимплеке Sk = Sk\x?, После нахождения такого подсимплекса полагаем Is+1 = 1Дя, Н1к+1 = Щ\x? и переходим к следующей итерации внутреннего цикла; (I) завершим выполнение итерации внутреннего цикла и вернемся к его началу; 4) завершим выполнение итерации внешнего цикла и вернемся к его началу. Данный алгоритм находит решение задачи проекции за конечное число итераций и для него показана глобальная "лучше чем линейная" скорость сходимости [6]. Как видно из приведенного описания, большую долю в алгоритме составляют матрично-векторные операции, что дает возможность использования графических ускорителей для ускорения его работы на задачах больших размерностей.

1.1. Вычислительные возможности GPU

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

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

NVIDIA (http://www.nvidia.com) содержат от 8 до 240 ядер на одном чипе. Это позволяет эффективно решать задачи с параллелизмом по данным с большим количеством арифметических операций по отношению к количеству обращений к памяти.

С целью реализации высокоуровневого доступа к аппаратным возможностям GPU для разработчиков неграфических приложений компанией NVIDIA была создана программно-аппаратная архитектура CUDA (Compute Unified Device Architechure). В состав CIJDA входят API (Application Programming Interface) для работы с GPU, специальный SDK (Software Development Kit) и компилятор языка С. Использование этой технологии позволяет писать масштабируемые параллельные программы, используя простое расширение языка С, при этом не подстраивая код под особенности конкретной модели GPU.

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

Для поддержки CUDA на аппаратном уровне NVIDIA было разработано семейство продуктов GeForce, Quadro и Tesla, построенных на вычислительной архитектуре Fermi. Основой этой архитектуры является массив процессоров, разделенный на группы многопоточных мультипроцессоров SM (Streaming Multiprocessors), каждый из которых в свою очередь содержит по восемь последовательных SP (Scalar Processor) процессорных ядер.

В отличие от ранних графических ускорителей архитектура Fermi предоставляет набор инструкций, воспринимаемых SP-ядром, близкий к тому набору, который программист ожидает от ядра обычного процессора. В частности, это полная поддержка целочисленной арифметики и арифметики с плавающей точкой. Используемая в наших экспериментах GeForce 280GTX содержит 240 последовательных процессорных ядер SP, сгруппированных в 30 мультипроцессоров SM.

Каждое многопоточное последовательное процессорное ядро SP может поддерживать до 128 потоков, причем создание потока, планирование и распределение ресурсов происходят полностью на устройстве. Каждый поток программы на CUDA отображается непосредственно на физический поток, находящийся на GPU. Таким образом, используемый GPU поддерживает до 30 720 параллельных потоков.

Мультипроцессор SM работает согласно модели SIMD (Single Instruction, Multiple Data), т.е. все потоки мультипроцессора выполняют одну и ту же инструкцию. Обмен данными между SP-ядрами мультипроцессора может происходить при помощи высокоскоростной области памяти (shared memory), объем которой для GF280 составляет 16 КБ на каждый мультипроцессор. Следует отметить, что для анонсируемой GF300 эта цифра увеличится до 64 КБ. Для хранения больших объемов данных на устройстве имеется встроенная (global) память объемом 1 ГБ и скоростью доступа 140 ГБ/с. Для сравнения — скорость доступа к оперативной памяти для DDR3 составляет порядка 20 ГБ/с (для модулей РСЗ-19200), Через общую память предполагается также обмен данными между различными мультипроцессорами.

Для сравнительных экспериментов по приросту производительности CPU/GPU были использованы машина на базе процессора Intel Xeon Е7330 с тактовой частотой 2.40 ГГц и многопоточные версии библиотек BLAS-ATLAS и LAPACK-ATLAS,

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

Основной концепцией реализации метода подходящих аффинных подпространств для GPU является перенос матрично-векторных операций на графический ускоритель, В алгоритме можно выделить несколько шагов, использующих такие операции. Во-первых, это проверка оптимальности построенного на г-м шаге внешнего цикла алгоритма приближения zH, которая производится с помощью системы неравенств ХгzH > ||zH||2,г = 1,..., N Введя матрицу X, столбцами которой являются векторы Хг, данное условие можно переписать в виде zHX > ||zH||2е, где e = (1,..., 1), Реализация такой проверки сводится к умножению матрицы на вектор и применению к полученному результату векторной версии алгоритма редукции [7] для выбора наибольшего значения элемента массива в полученном векторе.

Самой ресурсоемкой в вычислительном плане частью алгоритма является проекция пуля па базисное подпространство Ив, получаемое во внешнем цикле, где

Hb = {х = Y ^хг, ^^г = 1}

гев гев

и В С {1, 2,..., N} описывает набор базисных векторов. Количество элементов множества В обозначим через щ. Вычисление барицентрических координат нового приближения в имеющемся базисе Хв = {Хг, г € В} проводится для случая пь < п + 1 по формуле

^ = (Х|Хв )-1е/ет (ХТХв )-1е; (5)

для случая пь = п + 1 барицентрические координаты являются решением системы

Хв^ = 0, ет^ = 0. (6)

Последнее вызвано тем, что при ||В|| = п +1 матрица Ав — Х^Хв является вырожденной, Основной по трудоемкости операцией здесь является обращение матрицы Грамма Х^Хв, Исходя из особенностей алгоритма для реализации этой процедуры был выбран алгоритм блочного обращения [8, 9], в котором для блочной матрицы

А = Мц А12

А И21 А22

с невырожденной подматрицей А11 и невырожденной обратной матрицей Д = А22 — А21А-11А12 обратная матрица А-1 вычисляется как

A-1 + A-11A12D"1A2IA-11 -A-11A12D-1

-1 11 11 12 21 11 - 11 12 A = _n-1 4-1 n-1 . vJ

Таким образом, вычисление обратной матрицы Грамма происходит следующим образом: па к-м шаге алгоритма предполагается, что уже есть обратная левая верхняя подматрица матрицы Грамма Г— размерности (к — 1) х (к — 1), Применяя формулу (7) к подматрице Грамма Гк, получим обратную ей матрицу Г-1, затем переходим к шагу к + 1.

Использование блочного алгоритма позволяет снизить затраты на пересчет обратной матрицы Грамма, поскольку при добавлении нового вектора V в базис к уже существующей матрице Грамма Гк добавляются справа и снизу соответственно векторы

т = XT ■ v и тт соответственно. Этот алгоритм также использует операции BLAS третьего уровня [10], эффективные для исполнения на GPU для матриц большой размерности, В ходе вычислительных экспериментов было установлено, что для данной подзадачи предпочтительно применять смешанную модель вычислений CPU-GPU, Причиной этого могут быть относительно низкая производительность GPU на задачах малой размерности и большая начальная латентность вычислений (в среднем 8 мс для GF280), Поэтому пока подпространство Н имеет размерность менее некоторого порогового значения вычисления ведутся полностью на CPU, Так как основной операцией для подзадачи нахождения проекции нуля на Н является перемножение матриц, пороговое значение размерности подпространства оценивается для конкретной сборки из соображений прироста скорости операции матричного умножения и времени, затрачиваемого на копирование данных из оперативной памяти на память устройства. Для рассматриваемой сборки при размерности вектора 105 это составляет 128 векторов,

2. Численные эксперименты

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

2.1. Тестирование на центрированных симплексах 2.1.1. Тест ПОЛНЫЙ- СИМПЛЕКС- С

Первый тип тестовых множеств представлял собой (n — 1)-мерные симплексы с центром тяжести в начале координат, переведенные в n-мерное пространство сдвигом по n-й координате. Ниже приводится код на языке матрично-векторных вычислений octave (http://www.gnu.org/software/octave/), реализующий генерацию таких множеств:

rand("seed", 12345); dim=1000, nvec=dim; eps=0.01;

S=diag(100*rand(l,(dim-1) ));

S=[S zeros(dim -1, 1)];

c=sum(S,2)/dim;

S -= с * ones(1,dim);

S = [S ; eps*ones(l, dim)];

При решении задачи для построенных таким образом симплексов характерно последовательное включение всех вершин в базис. Поэтому в алгоритме при решении задач (5) и (6) эффективно применяется блочное обращение матрицы Грамма XTX, использующее полученную на предыдущем шаге обратную матрицу при добавлении новых векторов в базис. Для данного типа тестовых множеств были поставлены две серии численных экспериментов. Первая серия проводилась для множеств размерности n х n, где n

n

Видна степенная зависимость времени исполнения как CPU-, так и GPU-реализаций с общими показателями степени nCPU = 1.96 (кривая I), nGPU = 1.25 (кривая 2). Улучшение сходимости в GPU-реализации можно объяснить заменой скалярных вычислений

1е+07

| 1е+06 w

| 100000 х ч о с

cq С*

(D он

pq

10000 1000 100

10

GPU-

CPU —

1

1е+0.6

100 1000

Размерность вектора

10000

0 £

к

к <D К 100000

ч

о

с

¡3 CQ

t* 10000

<D

Он

1000

1000 10000 100000 Размерность вектора

Рис. 1. Центрированные симплексы с изменением количества векторов и размерности пространства (тест ПОЛНЫЙ-СИМПЛЕКС-С)

Рис. 2. Центрированные симплексы с фиксированным количеством векторов и изменением размерности пространства (тест ПОДСИМПЛЕКС-С)

векторными. Уменьшение показателя степени для GPU на 0.71 а не на 1 происходит в силу того, что физически на мультипроцессоре в один такт обсчитывается не весь вектор, а только его часть. Также видно, что для задач размерности меньше 1000 х 1000 более эффективной является CPU-реализация. Это можно объяснить небольшой утилизацией вычислительных ресурсов GPU и его латентностью вычислений.

2.1.2. Тест ПОДСИМПЛЕКС-С

Вторая серия численных экспериментов проведена для множеств, построенных по тому же алгоритму, но содержащих фиксированное количество векторов 103 (рис. 2). Изменялась только размерность векторов, которая варьировалась от 103 до 105. На этой серии тестов алгоритм показал практически линейную зависимость времени нахождения проекции от размерности векторов. Наибольший прирост производительности, составляющий примерно 9 раз, получен для векторов максимальной размерности порядка 105, что позволяет предположить дальнейшее увеличение преимуществ GPU с ростом размерности. Общие показатели степени для данной серии алгоритмов составляют 1.4 для GPU и 2 для CPU.

2.2. Тестирование на случайном множестве

Это множество тестов строилось таким образом, что исходный набор данных состоял из п векторов £ = (£1, £2, • • •, £ш) построенных следующим образом:

{вса1е 1 • ^ — , если г = 1, 2, ..., п — 1, вса1е2 • я + в^г/£, если г = п,

где я — независимые равномерно распределенные на [0,1] случайные величины, вса1е1, вса1е2 — константы масштабирования, в^г/£ — константа сдвига. Для каждого вектора были проведены также масштабирование последней координаты с коэффициентом 10-3 и две серии тестов — с изменением размерности и количества векторов (ПОЛНЫИ-СИМПЛЕКС-С) и с изменением только размерности (ПОДСИМПЛЕКС-С) (рис. 3).

Размерность вектора Размерность вектора

Рис. 3. Случайно сгенерированные симплексы с изменением размерности пространства и количества векторов — ПОЛНЫЙ-СИМПЛЕКС-С (а) и изменением размерности пространства и фиксированным количеством векторов — ПОДСИМПЛЕКС-С (б)

Код для octave, описывающий генерацию множеств, показан ниже:

randO'seed", 12345); scale = [100,0.1]; shift = 500;

S = rand(dim - 1, nvec);

S = [scale(l) * (S - sum(S,2)/nvec * ones(l, nvec)); scale(2)*rand(l,nvec)+shift ];

Работа алгоритма на сгенерированных таким образом множествах характеризуется присутствием в подавляющем большинстве тестов фазы выведения векторов из базиса. При этом для матрицы Грамма, получившейся при изъятии вектора из базиса, становится невозможным применение ускорения обращения с использованием блочного алгоритма. Несмотря на это, алгоритм демонстрирует характер сходимости, схожий с предыдущими измерениями. Показатель сходимости для тестов ПОЛНЫЙ-СИМПЛЕКС-С равен 1.85 для GPU и 3.4 для CPU. Для тестов ПОДСИМПЛЕКС-С этот параметр составляет 1.87 для GPU и 3.14 для CPU.

Заключение

Результаты проведенных численных экспериментов подтверждают возможность практического применения GPU-реализации метода подходящих аффинных подпространств для решения задач проекции большой размерности (105). Основное ограничение на размерность задачи накладывает количество оперативной памяти на графическом ускорителе, что проявляется в том, что для тестовых задач, решение которых требовало использования всей памяти графического ускорителя, точка увеличения показателя степени в характеристике сходимости алгоритма так и не была достигнута. Применение реализации метода подходящих аффинных подпространств для графических ускорителей неэффективно для задач проекции малых размеров, что объясняется низкой утилизацией вычислительных ресурсов GPU. Как показали эксперименты, преимущества в производительности над CPU-версией начинается для задач с размерностью более 1500 х 1500.

Список литературы

[1] Нурминский Е.А. Численные методы выпуклой оптимизации. М.: Наука, 1991. 160 с.

[2] Boyd S., Vandenberghe S. Convex Optimization. Cambridge Univ. Press, 2004.

[31 CensorY., Jiang M., Louis A.K. Mathematical Methods in Biomedical Imaging and Intensity-Modulated Radiation Therapy (IMRT). Pisa, Italy: Birkhauser-Verlag, 2008. 522 p.

[4] Boyarchenkov A.S., Potashnikov S.I. Using graphic accelerators and CUDA technology for solving problems of molecular dynamic // Comput. Methods and Programming. 2009. Vol. 10, No. 1. P. 9-23.

[5] Stone S., Haldar J. Accelerating advanced MRI reconstructions on GPUs // 5th Conf. on Computing Frontiers. ACM. N.Y. USA, 2008. P. 261-272.

[6] Нурминский E.A. О сходимости метода подходящих аффинных подпространств для решения задачи о наименьшем расстоянии до симплекса // Журн. вычисл. математики и матем. физики. 2005. Т. 45, вып. 11. С. 1996-2004.

[7] Lin Н.Х., Sips H.J. Parallel vector reduction algorithms and architectures // J. of Parallel and Distributed Comput. 1988. Vol. 5, No. 2. P. 103^130.

[81 Castillo M., Chan E., Francisco D. Making Programming Synonymous With Programming for Linear Algebra Libraries. Technical Report, 2008.

[9] Нейдеккер X., Магнус Я.P. Матричное дифференциальное исчисление с приложениями к статистике и эконометрике. М.: Физматлит, 2002. 495 с.

[101 Barrachina S., Castillo М. Evaluation and tuning of the level 3 CUBLAS for graphics processors // IPDPS. 2008. P. 1-8.

Поступила в редакцию 10 октября 2010 г., с доработки — 31 мая 2011 г.