CC BY
27
УДК 004.272:519.67:004.942 В.Ю. Воронов1
МЕТОД АВТОМАТИЧЕСКОГО ВЫБОРА И НАСТРОЙКИ ПАРАМЕТРОВ РЕШАТЕЛЯ РАЗРЕЖЕННОЙ СЛАУ
В работе представлен метод автоматического выбора метода решения разреженной СЛАУ и значений его параметров на основе анализа матрицы коэффициентов СЛАУ. Выбор осуществляется с целью максимизации метрики производительности. Представлена формальная модель производительности решателя СЛАУ, на основе которой сформулирован алгоритм построения метода, его обучения на прецедентах и применения. Выполнена апробация метода для итерационных решателей СЛАУ из параллельной библиотеки PETSc на платформе МВС-100К МСЦ РАН.
Ключевые слова: метод опорных векторов, автоматическая настройка производительности приложений, решение разреженных систем линейных алгебраических уравнений, PETSc.
Введение. Вычислительный эксперимент с использованием высокопроизводительных платформ является важным элементом прикладной науки. Требование эффективного использования ресурсов вычислительных платформ делает актуальной проблематику автоматической настройки производительности (automatic performance tuning) параллельных приложений. Эта тематика подразумевает создание алгоритмов настройки приложений на специфику решаемой задачи и особенности целевой высокопроизводительной архитектуры. Решение разреженных систем линейных алгебраических уравнений (разреженных СЛАУ) является одной из важных задач, лежащих в основе широкого круга высокопроизводительных научных приложений. Известно, что на производительность численного решения задачи существенно влияет выбор методов решения и специфика конкретной вычислительной платформы [1, 2]. Эти свойства характерны при использовании итерационных методов решения разреженных СЛАУ. В ряде случаев определение наилучшего метода решения конкретного класса задач априори недоступно и требует экспериментальной настройки. Таким образом, проблема определения наилучшего метода решения заданного класса разреженных СЛАУ для целевой архитектуры является одной из важных задач настройки производительности.
Существующие методы настройки производительности методов решения разреженных СЛАУ можно разделить на две группы — на методы, настраивающие алгоритмы на основе выявления архитектурных особенностей платформы и на основе анализа матрицы коэффициентов СЛАУ. Наиболее известные инструментальные системы, реализующие методы первой группы, — пакеты ATLAS [3] и OSKI [4]. В них реализована автоматическая настройка функций пакета базовых операций линейной алгебры BLAS, настройке подвергаются операции вида "матрица-вектор" и "матрица-матрица" с целью более эффективного использования структуры кэш-памяти и других особенностей платформы. В пакете ATLAS настройка осуществляется двумя основными способами — варьирования параметров функций с выбором оптимальных настроек на этапе выполнения приложения и создания множества нескольких реализаций функции с дальнейшим выбором определенной реализации функции для данной платформы. Пакет OSKI может учитывать априорную информацию пользователя о структуре заполнения матрицы, о характере и потоке вычислений в приложении. Образец подхода настройки производительности решателя на основе анализа матрицы коэффициентов СЛАУ — работа [5], авторы которой решают задачу прогнозирования числа обусловленности матрицы. Известно, что алгоритмы его вычисления требуют порядка 0(п3) операций, вычисление оценки на число обусловленности снижает сложность до 0(п2) операций (такова, например, группа алгоритмов ***CON в пакете LAPACK), что тем не менее является высоким показателем, в особенности если требуется определить, является матрица плохо обусловленной или нет. В работе описан алгоритм построения регрессии на обучающей выборке СЛАУ, алгоритм прогнозирования, а также процедура отбора признаков в заданном пространстве признаков СЛАУ. Данный подход находит применение в системе IPRS, разрабатываемой авторами в Университете Кентуки для автоматического выбора предобусловливателя СЛАУ.
Предлагаемый в данной работе алгоритм автоматического выбора метода решения разреженной СЛАУ и настройки его параметров основан на применении методов машинного обучения.
1Факультет ВМиК МГУ, асп., emaikbasravQangel.eme.msu.ru.
* Работа выполнена при поддержке грантов РФФИ № 08-07-00445-а, 08-07-12081.
1. Модель производительности решателя разреженной СЛАУ. Интуитивно представляя задачей выбор наилучшего метода решений и значений его набора параметров в зависимости от характеристик входных данных, представим формальную модель предметной области и дадим ее интерпретацию для процесса решения разреженной СЛАУ с действительными коэффициентами А^х = Ь. Для анализа производительности в предлагаемой модели использование метода решения разреженной СЛАУ представляется функцией сложности ^):
Я '■
где ¿^ = (/1, /2,..., ¡п)т — вектор признаков, соответствующих решаемой задаче А¿; Sj = («1, «2, • • • ... , «то)т — вектор параметров используемого метода решения; Нк = (/11, /г,2, • • •, — вектор характеристик вычислительной платформы; Р= (р) — значение метрики производительности выбранного метода решения заданной СЛАУ на данной платформе. Рассмотрим более подробно представленную модель.
• Переход от матрицы СЛАУ к признаковому пространству необходим для понижения размерности данных и позволяет оперировать со СЛАУ различной размерности. Под признаком понимается функция от матрицы коэффициентов СЛАУ / : Жп —> Ж.
• Вектор параметров метода — набор значений, кодирующий применяемые для решения СЛАУ тип итерационного метода, тип предобусловливателя, значения параметров итерационного метода и предобусловливателя.
• Вектор характеристик вычислительной платформы — набор параметров, описывающих вычислительную платформу. В общем случае каждый параметр является значением некоторой функции, описывающей производительность заданной платформы.
• Метрика производительности — действительная величина, характеризующая производительность решения СЛАУ с помощью выбранного метода. В простейшем случае это время решения разреженной СЛАУ, на практике широко применяется метрика произведения времени решения СЛАУ на пиковый объем занимаемой оперативной памяти.
Представленная модель производительности позволяет формулировать задачу определения наилучшего метода решения заданной разреженной СЛАУ из некоторого множества методов В, на котором достигается минимум функции сложности:
Б* =a,тgmmQ(FuSj,Hk). (1)
Аналитическое задание функции позволяет применять для решения (1) оптимизационные методы. В общем случае вид функции С} априори неизвестен, а задается набором значений д^, полученных при решении некоторого набора СЛАУ. Для решения (1) предлагается применить метод обучения на прецедентах [6].
2. Алгоритм выбора решателя СЛАУ и настройки его параметров. Решение (1) будет реализовано в несколько этапов. На первом этапе алгоритма выполняется построение функции, аппроксимирующей С} на некотором обучающем наборе СЛАУ Т. С использованием построенной зависимости на втором этапе определяется производительность методов решения для рассматриваемой СЛАУ, из которых выбирается наилучший.
Для того чтобы увеличить точность определения метода решения СЛАУ, необходимо учитывать специфику применения итерационных методов решения разреженных СЛАУ. Прежде всего, необходимо сузить множество поиска методов за счет прогнозирования сходимости данного метода при решении заданной разреженной СЛАУ. Для этого фиксируются условия сходимости метода (требуемый допуск нормы невязки, максимальное количество итераций метода) и задача классификации методов по сходимости формулируется в терминах бинарной классификации, для решения которой применяется метод опорных векторов [7, 8]. Задача бинарной классификации векторов {х| х^ € Жп, г = 1,...,/} для классов у^ = { — 1,1} может быть решена методом С-классификации опорных векторов (С-ЭУС)
путем решения задачи квадратичного программирования
тт + <7
У^ШТф(щ) + Ъ) ^ 1 " Сг,
Сг^О,
либо решения двойственной задачи
1 Т г1 т тт-а ц/а — е а,
а 2
УГ« = 0, (2)
О < а,- < С,
где = фТ(щ)ф(ху) — ядро, (¿^ = х^-), С ^ 0, е — единичный вектор. Результат
решения (2) — функция
/(х) = sgní + Ъ\, (3)
4=1 '
определяющая принадлежность элемента х к классу у. Для построения классификатора может использоваться также метод ^-классификации опорных векторов {и-ЭУС).
Для построения аппроксимации функции <3 на обучающем наборе задач применяются методы регрессии опорных векторов: для заданного набора значений функции
((XI, дх),..., (хп,дп)), х^!1, д* € Ж,
— метод е-регрессии опорных векторов (е-ЗУИ) в форме задачи квадратичного программирования
I I
mm -wvw ■ ■
i= 1 i= 1
T,
1
w'b'l'k* 2
wтф(щ) + b- +
Zi - wтф(щ) - b < e + i*,
% — 1
и в форме двойственной задачи
1
min - (а — a*)TQ(a — а*) + е V^ Zi(a — а*), а,а* 2
%=1
I
= (4)
г= 1
«j ^0, а* ^ С, i — 1 ►>>>►/►
где аналогично с классификацией Ä"(xj,xj) = <f>(x.J)<f)(x.j). В результате решения (4) искомая аппроксимирующая функция может быть записана в виде
I
f(x) = ^(^«г + Q!*)K(Xj,x) + Ь. i= 1
По аналогии с методом I/-SVC введением в модели (2) параметра v определяется г/-регрессия опорных векторов (V-SVR).
В применении методов опорных векторов важную роль играет использование функции-ядра. Применение ядра вместо скалярного произведения позволяет решать при классификации проблему линейной разделимости и получать нелинейную аппроксимирующую функцию в результате регрессии. В качестве функции-ядра наиболее часто применяются следующие функции:
• линейная функция = х^х^ как частный случай линейного метода опорных векторов;
• радиальная базисная функция = И2, 7 > 0;
• полиномиальная функция = + г)а, 7 > 0;
• сигмоидная функция = Ш^х^х^- + г).
В качестве параметров функции-ядра тут выступают г, й, 7. Для численного решения задач (2), (4) используются алгоритмы, основанные на декомпозиции, в частности ЭМО-метод, применимость, сложность и критерии сходимости которого для методов опорных векторов исследованы в работах [9, 10].
Для повышения точности построения аппроксимирующей функции предложено использовать следующий метод. Пусть набор параметров метода решения в г представлен в виде объединения непересекающихся множеств:
^ = У ■
Тогда аппроксимирующая функция может быть найдена в виде композиции функций, определенных на подмножествах Sit:
я= Е (5)
Решение (4) тем самым сведено к построению множества аппроксимирующих функций д^ меньшей размерности, что означает уменьшение вычислительной сложности построения аппроксимации и ведет к повышению точности. Разбиение на подмножества Sjt проводится так, что каждое подмножество соответствует комбинации параметров одного итерационного метода и предобусловливателя.
Таким образом, общая последовательность этапов обучения предлагаемого алгоритма следующая.
• Определить обучающую выборку СЛАУ Т = {¿<¿1, множество методов решения СЛАУ В и сформировать набор значений функции сложности д = {Р^^ ///,.. -Ру/г) путем решения СЛАУ на целевой платформе.
• Каждому элементу из набора значений функции сложности сопоставить класс
{1, если для СЛАУ г метод решения ] сходится, — 1, если для СЛАУ i метод решения ] расходится.
• Создать и обучить классификатор С : (Рг, Я*;) —> (су/г). Применение классификатора позволяет определить множество сходящихся методов Б, для которого заведомо имеет смысл искать решение, тем самым сузить пространство поиска методов решения СЛАУ.
• Построить аппроксимирующую функцию (Fi, Sj, Нк) —> {Рг:1ь) на множестве сходящихся методов решения Б в виде суммы аппроксимирующих функций на подпространствах .
Последовательность этапов алгоритма определения наилучшего метода решения СЛАУ, не входящей в тестовую выборку, следующая:
• формирование из входной СЛАУ набора признаков {¿<¿1;
• определение для набора методов решения Sj подмножества Sj, для которого выполняются условия сходимости путем применения построенного классификатора;
• определение производительности Р.на множестве методов решения Б* с использованием построенной аппроксимации функции качества, выбор наиболее производительного решателя.
При использовании методов машинного обучения высокая точность метода на обучающем наборе в общем случае не гарантирует высокой точности решения на произвольных СЛАУ. Для увеличения обобщающей способности построения аппроксимации применяются следующие методы:
• предобработка данных (в самом распространенном случае ограничиваются приведением признаков выборки к масштабу [0,1] или [—1,1]);
• использование в качестве критерия точности для построения классификатора и аппроксимирующей функции метода скользящего контроля;
• поиск параметров методов (2) и (4) (предпочитая те параметры метода, которые максимизируют точность скользящего контроля на обучающей выборке СЛАУ). На практике используются две распространенные статегии отыскания параметров 7 и С (либо 7 и ь>, в зависимости от типа метода опорных векторов) — перебор параметров на некоторой сетке значений = (2~По+г, 2~По+:г) либо иерархический поиск параметров на измельчающейся сетке значений параметров (поиск состоит из нескольких итераций, на каждой итерации параметры ищутся в более густой сетке в окрестности оптимальной точки, найденной на предыдущей итерации).
3. Апробация предложенного алгоритма. Для оценки применимости разработанного метода выполнена апробация на основе методов решения из параллельной библиотеки научных вычислений РЕТЭс [11]. Множество методов решения из метода решения соответствуют группе итерационных методов, предобусловливателей и значений параметров. В качестве итерационных решателей рассматривались методы минимальных обобщенных невязок (GMR.ES) и бисопряженных градиентов со стабилизацией (BiCGStab). В качестве предобусловливателей использовались методы неполной декомпозиции (1Ы1), обычный и блочный методы Якоби, аддитивный метод Шварца. В качестве параметров методов рассматривались:
Перечень основных характеристик СЛАУ
СЛАУ n nnz p(A)
GHS_indef/exdata_l 6001 2269500 6,742507e+06
FEMLAB/sme3Db 29067 2081063 l,128736e+05
Boeing/bcsstk39 46772 2060662 2,783841e+04
Bova/rmalO 46835 2329092 8,735202e+04
GHS_psdef/vanbody 47072 2329056 3,642314e+07
Boeing/ct20stif 52329 2600295 l,983118e+06
Nasa/nasasrb 54870 2677324 l,403127e+05
Rudnyi/water_tank 60740 2035281 l,003749e+05
PAESEC/H20 67024 2216736 3,344466e+03
GHS_psdef/oilpan 73752 2148558 5,907112e+04
FEMLAB/poisson3Db 85623 2374949 2,783841e+04
Oberwolfach/filter3D 106437 2707179 l,570727e+04
Schenk_ISEI/barrier2-l 113076 2129496 3,952821e+02
Schenk JSEI/barrier2-2 113076 2129496 4,991775e+02
Schenk JSEI/barrier2-3 113076 2129496 5,497282e+02
Schenk_ISEI/barrier2-4 113076 2129496 3,884241e+02
Schenk JSEI/barrier2-10 115625 2158759 9,190972e+02
Schenk_ISEI/barrier2-ll 115625 2158759 9,256846e+02
Schenk JSEI/barrier2-12 115625 2158759 9,206099e+02
Schenk_ISEI/barrier2-9 115625 2158759 9,275252e+02
vanHeukelum/cagel2 130228 2032536 9,891092e+00
Schenk_ISEI/para-4 153226 2930882 3,981658e+03
Schenk JSEI/рагаЛО 155924 2094873 7,843124e+04
Schenk_ISEI/para-5 155924 2094873 9,282488e+04
Schenk_ISEI/para-6 155924 2094873 l,560355e+04
Schenk_ISEI/para-7 155924 2094873 5,128946e+03
Schenk_ISEI/para-8 155924 2094873 2,873521e+03
Schenk_ISEI/para-9 155924 2094873 3,082653e+04
Schenk JBMNA/c-big 345241 2340859 2,819281y+05
GHS_indef/darcy003 389874 2097566 4,359254e+04
GHS_indef/mario002 389874 2097566 4,359254e+04
GHS_indef/helm2d03 392257 2741935 6,550870e+02
Sandia/ASIC_680k 682862 2638997 l,515031e+07
Примечание, п — размерность СЛАУ, ппг — число ненулевых элементов матрицы СЛАУ, р(А) — число обусловленности матрицы СЛАУ.
• параметр перезапуска для итерационного метода ОМИБЭ;
• тип алгоритма переупорядочения, число уровней заполнения, коэффициент заполнения и порог отсечения ведущего элемента в предобусловливателях, основанных на неполной факторизации;
• число блоков, тип метода построения блочного предобусловливателя Якоби, тип блочного решателя и его настройки;
• число блоков, тип метода, параметр перекрытия блоков для предобусловливателя аддитивного Шварца, а также тип блочного решателя и его настройки, аналогичные настройкам в неполной факторизации.
Отметим, что реализация переборной стратегии для такого поискового пространства не может дать качественного результата, поскольку, во-первых, некоторые параметры вещественные и, во-вторых, в общем случае не наблюдается непрерывной зависимости производительности от изменений значения параметра или малого изменения вещественного параметра.
Для построения признакового пространства СЛАУ использован пакет АпаМос! [12]. Следует отметить, что переход к признаковому пространству является трудоемкой вычислительной операцией, в данной работе выбраны те признаки, общая вычислительная сложность определения которых заведомо не превышает сложности решения разреженной СЛАУ.
В качестве метрики производительности Р^к использовалось произведение времени работы решателя и объема использованной оперативной памяти.
Для реализации методов опорных векторов применялся пакет операций машинного обучения ЫЬЭУМ [13]. Выборка СЛАУ для вычислительного эксперимента была взята из репозитория разреженных матриц университета Флориды [14], в котором собраны СЛАУ, возникающие из ряда научных и промышленных проблем; выборка составляет 35 СЛАУ, размерность матриц выборки в диапазоне 103-105, разреженность порядка 0,004-5%. Общие характеристики матриц выборки приведены в таблице. Вся выборка была случайным образом разделена на обучающую, в которую вошло 23 СЛАУ, и тестовую, куда вошли оставшиеся 12 СЛАУ.
Рис. 1. График зависимости точности скользящего контроля классификатора на обучающем (сплошная линия) и тестовом (пунктирная линия) наборах от
параметров классификатора
Для всех СЛАУ правая часть была получена умножением матрицы на известное решение х = = (1,...,1)т, начальное приближение решения ж о = (0,..., 0)т, критерий сходимости — достижение нормы невязки решения ||г||2 < Ю-9, либо итерационный процесс прерывается по достижении
1000 итераций. Вычислительный эксперимент проводился на платформе МВС-100К МЦС РАН, использована библиотека РЕТЭс версии 2.3.3-р15.
Обратимся к результатам апробации метода. На рис. 1 представлены результаты обучения и настройки параметров классификаторов сходимости. Представленный график иллюстрирует, что классификаторы с высокой точностью скользящего контроля для обучающей выборки дают также высокий результат точности на тестовой выборке. Высокую точность полученного классификатора иллюстрирует построенная ИОС-кривая [15], представленная на рис. 2. Величина площади под кривой данного классификатора составляет 0,9681. Таким образом, предложенный метод с высокой степенью точности выделяет из множества методов решения сходящиеся методы.
1 1111 и
Рис. 2. График 110С-кривой для построенного классификатора
Рис. 3. Гистограмма процента задач тестового набора, ошибка прогнозирования на которых не превышает величины устанавливаемого порога, в зависимости от
порога
Для оценки точности прогнозирования аппроксимирующей функции проведен анализ относительной ошибки прогнозирования на тестовой выборке СЛАУ. Результаты представлены на рис. 3. Отметим, что для 26% всех оцениваемых методов решения СЛАУ ошибка прогнозирования составила менее 5%, что является хорошим результатом.
Заключение. Предложенный метод может быть обобщен не только для решения общей задачи выбора решателя СЛАУ и его параметров в библиотеке, но и для выбора решателя среди нескольких
библиотек и даже среди нескольких вычислительных платформ. Представленный метод автоматического выбора решателя СЛАУ в настоящее время находит реализацию в инструментальной системе, разрабатываемой на факультете ВМиК МГУ. Данная инструментальная система даст пользователям массивно-параллельных вычислительных платформ возможность по представленным ими разреженным СЛАУ получать рекомендацию, какой метод решения из заданного множества наиболее производителен и каковы значения его параметров.
СПИСОК ЛИТЕРАТУРЫ
1. Попова Н. Н., Воронов В.Ю., Джосан О.В., Медведев М. А. Сравнительный анализ эффективности параллельных вычислений с использованием современных параллельных математических библиотек на примере решения систем линейных уравнений // Труды Всероссийской научной конференции "Научный сервис в сети Интернет-2004". М.: Изд-во МГУ, 2004. С. 146-149.
2. Gupta A., George Т., Sarin V. An experimental evaluation of iterative solvers for large SPD systems of linear equations. IBM T.J. Watson Research Center. Tech. Rep. RC 24479. 2008.
3. Whaley R., Petitet A., Dongarra J. Automated empirical optimizations of software and the ATLAS project // Parallel Computing. 2001. 27. N 1. P. 3-35.
4. Vuduc R., Demmel J., Yelick K. OSKI: A library of automatically tuned sparse matrix kernels // J. Physics: Conference Series. 2005. 16. N 1. P. 521-530.
5. Xu S., Zhang J. A new data mining approach to predicting matrix condition numbers // Commun. Inform. Systems. 2004. 4. N 4. P. 325-340.
6. Вапник B.H. Восстановление зависимостей по эмпирическим данным. М.: Наука, 1979.
7. Cortes С., Vapnik V.N. Support-vector networks // Machine Learning. 1995. 20. P. 273-297.
8. S mola A., Scholkopf B. A tutorial on support vector regression // Statistics and Computing. 2004. 14. N 3. P. 199-222.
9. Fan R., Chen P., Lin C. Working set selection using second order information for training Support Vector Machines //J. Machine Learning Research. 2005. 6. P. 1889-1918.
10. Chen P., Fan R., Lin C. A study on SMO-type decomposition methods for support vector machines // IEEE Transactions on Neural Networks. 2006. 17. N 4. P. 893.
11. http://www.msc.anl.gov/petsc-2.
12. Eijkhout V., Fuentes E. A proposed standard for numerical metadata. Innovative Computing Laboratory. University of Tennessee. Tech. Rep. ICL-UT-03-02. 2003.
13. http://www.csie.ntu.edu.tw/ cjlin/libsvm.
14. Davis T. University of Florida sparse matrix collection // NA Digest. 1997. 97. N 23. P. 7.
15. Davis J., Goadrich M. The relationship between Precision-Recall and ROC curves // Proc. of the 23d Intern. Conf. on Machine Learning. N.Y., USA: ACM New York, 2006. P. 233-240.
Поступила в редакцию 10.12.08
