The Use of MPI and OpenMP Technologies for Subsequence Similarity Search in Very Large Time Series on Computer Cluster System with Nodes Based on the Intel Xeon Phi Knights Landing Many-core Processor
Yana Kraeva Mikhail Zymbler

	
УДК 004.272.25, 004.421, 004.032.24










M[1]>m#1
N@m0pt@












СОВМЕСТНОЕ ИСПОЛЬЗОВАНИЕ ТЕХНОЛОГИЙ MPI и OpenMP 

ДЛЯ ПАРАЛЛЕЛЬНОГО ПОИСКА ПОХОЖИХ ПОДПОСЛЕДОВАТЕЛЬНОСТЕЙ 

В СВЕРХБОЛЬШИХ ВРЕМЕННЫХ РЯДАХ НА ВЫЧИСЛИТЕЛЬНОМ КЛАСТЕРЕ 

С УЗЛАМИ НА БАЗЕ МНОГОЯДЕРНЫХ ПРОЦЕССОРОВ 

INTEL XEON PHI KNIGHTS LANDING
    Я.А. Краева[Южно-Уральский государственный университет (национальный исследовательский университет), кафедра системного программирования, просп. Ленина, 76, 454080, Челябинск; магистрант, e-mail: mailto:kraevaya@susu.rukraevaya@susu.ru], М.Л. Цымблер[Южно-Уральский государственный университет (национальный исследовательский университет), кафедра системного программирования, просп. Ленина, 76, 454080, Челябинск; к.ф.-м.н., доцент, e-mail: mailto:mzym@susu.rumzym@susu.ru]
    Received: date / Accepted: date
==============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





В настоящее время поиск похожих подпоследовательностей требуется в широком спектре приложений интеллектуального анализа временных рядов: моделирование климата, финансовые прогнозы, медицинские исследования и др. В большинстве указанных приложений при поиске используется мера схожести Dynamic Time Warping (DTW), поскольку на сегодня научное сообщество признает меру DTW одной из лучших для большинства предметных областей. Мера DTW имеет квадратичную вычислительную сложность относительно длины искомой подпоследовательности, в силу чего разработан ряд параллельных алгоритмов ее вычисления на устройствах FPGA и многоядерных ускорителях с архитектурами GPU и Intel MIC. В данной статье предлагается новый параллельный алгоритм для поиска похожих подпоследовательностей в сверхбольших временных рядах на кластерных системах с узлами на базе многоядерных процессоров Intel Xeon Phi поколения Knights Landing (KNL). Вычисления распараллеливаются на двух уровнях: на уровне всех узлов кластера — с помощью технологии MPI, в рамках одного узла кластера — с помощью технологии OpenMP. Алгоритм предполагает использование дополнительных структур данных и избыточных вычислений, позволяющих эффективно задействовать возможности векторизации вычислений на процессорных системах Phi KNL. Эксперименты, проведенные на синтетических и реальных наборах данных, показали хорошую масштабируемость алгоритма.





Я.А. Краева, М.Л. ЦымблерСовместное использование технологий MPI и OpenMP для параллельного поиска…



§ ВВЕДЕНИЕ.



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

На сегодня мера DTW (Dynamic Time Warping, динамическая трансформация временной шкалы) <cit.> признается научным сообществом лучшей мерой схожести временных рядов для многих предметных областей <cit.>, поскольку позволяет адекватно измерять схожесть временных рядов, которые имеют различные длины либо сдвиг, растяжение или сжатие друг относительно друга. Поскольку мера DTW имеет квадратичную вычислительную сложность относительно длины поискового запроса, были предложены методы снижения вычислительной нагрузки: индексирование <cit.>, заблаговременный отказ от избыточных вычислений, повторное использование результатов ранее выполненных вычислений <cit.> и др. Последовательный алгоритм UCR-DTW <cit.>, интегрирующий множество указанных техник, является на сегодня, вероятно, самым быстрым последовательным алгоритмом поиска похожих подпоследовательностей <cit.>. В то же время предложены параллельные алгоритмы поиска похожих подпоследовательностей для устройств FPGA <cit.>, многоядерных ускорителей архитектур GPU <cit.> и Intel Many Integrated Core (MIC) <cit.>.

Рассматриваемые в данной статье многоядерные системы архитектуры Intel MIC являются конкурентоспособной альтернативой более распространенным системам NVIDIA GPU и FPGA. Системы Intel MIC основаны на распространенной архитектуре Intel x86 и поддерживают соответствующие модели и инструменты параллельного программирования. Устройства MIC обеспечивают большое количество энергоэффективных вычислительных ядер с малой (относительно обычных многоядерных процессоров Intel) тактовой частотой, обладающих высокой пропускной способностью локальной памяти и 512-битными векторными регистрами. В настоящее время компания Intel предлагает два поколения MIC-устройств, имеющих базовое имя Xeon Phi: сопроцессор Knights Corner (KNC) <cit.> с 57–61 ядрами и самостоятельная процессорная система Knights Landing (KNL) <cit.> с 64–72 ядрами соответственно. Системы MIC, как правило, дают наибольшую производительность в таких приложениях, где имеет место большой объем данных, вовлеченных в вычисления (десятки миллионов элементов и более), которые могут быть векторизованы компилятором <cit.>. Под векторизацией понимается способность компилятора заменить несколько скалярных операций в теле цикла с фиксированным количеством повторений в одну векторную операцию <cit.>. 

Параллельный поиск похожих подпоследовательностей в сверхбольших временных рядах (порядка сотен миллиардов точек) на кластерных системах рассмотрен в работах <cit.>. В статье <cit.> предложен способ 
распараллеливания последовательного алгоритма поиска похожих подпоследовательностей UCR-DTW для вычислительного кластера на основе использования фреймворка Apache Spark. В рамках данного фреймворка приложение запускается как задача, координируемая мастер-узлом вычислительного кластера, на котором установлен соответствующий драйвер. Алгоритм предполагает фрагментацию временного ряда, однако количество фрагментов не совпадает
с количеством вычислительных узлов кластерной системы. Фрагменты сохраняются в виде отдельных файлов, доступных всем узлам в силу использования распределенной файловой системы HDFS (Hadoop Distributed File System). Каждый фрагмент обрабатывается отдельной нитью, реализующей последовательный алгоритм UCR-DTW. Количество процессов, запускаемых на узле кластера, равно количеству ядер процессора на данном узле. В экспериментах разработанный алгоритм опережает UCR-DTW в 2–5 раз, однако представленные результаты его запусков на кластере из 2–6 узлов не позволяют говорить о хорошей масштабируемости решения.

В предыдущих исследованиях авторов данной статьи <cit.> разработан параллельный алгоритм поиска похожих подпоследовательностей на вычислительном кластере с узлами на базе сопроцессоров Phi KNC. Предложенная вычислительная схема <<CPU+Phi>> предполагает фрагментацию ряда по узлам кластера и следующее разделение функций между процессором и KNC внутри узла. KNC выполняет вычисление меры DTW для подпоследовательностей, загруженных в его память процессором, и выбирает подпоследовательность, максимально похожую на образец поиска. Процессор исключает из обработки заведомо непохожие на поисковый запрос подпоследовательности на основе вычисления нижних оценок формирует пакеты подпоследовательностей для выгрузки данных на сопроцессор. Предложенная схема вычислений обеспечивает высокую производительность, однако слабо задействует возможности векторизации вычислений Phi KNC. Масштабируемость, близкая к линейной, достигается алгоритмом в случае, когда поисковый запрос имеет длину от 4000 точек, что относительно редко встречается на практике. Помимо этого, данная схема не может быть применена для случая кластерной системы на базе узлов с системами следующего поколения Phi, поскольку KNL является независимым многоядерным устройством, запускающим приложения только в нативном режиме. В работе <cit.> авторами предложен подход к эффективному поиску похожих подпоследовательностей в оперативной памяти процессора Phi KNL. 

В данной статье предлагается новый параллельный алгоритм PhiBestMatch для поиска похожих подпоследовательностей в сверхбольших временных рядах на кластерных системах с узлами на базе многоядерных процессоров Intel Xeon Phi поколения Knights Landing. Вычисления распараллеливаются на двух уровнях: на уровне всех узлов кластера — с помощью технологии MPI, в рамках одного узла кластера — с помощью технологии OpenMP. Алгоритм предполагает использование дополнительных структур данных и избыточных вычислений, позволяющих эффективно задействовать возможности векторизации вычислений на процессорных системах Phi KNL.

Статья организована следующим образом. В разделе <ref> приводится формальная постановка задачи и дано краткое описание последовательного алгоритма UCR-DTW, используемого в качестве базиса нового параллельного алгоритма. Раздел <ref> содержит описание алгоритма PhiBestMatch. В разделе <ref> описаны вычислительные эксперименты, исследующие эффективность предложенного алгоритма. Заключение резюмирует результаты, полученные в рамках исследования.



§ ПОСТАНОВКА ЗАДАЧИ.





 §.§ Формальные определения и нотация.




Дадим определения используемых терминов в соответствии с работой <cit.>. Временной ряд (time series) T — это  хронологически упорядоченная последовательность вещественных значений: T=(t_1, t_2,…,t_m), t_i ∈ℝ. Число m обозначается | T | и называется длиной временного ряда.



Мерой схожести (similarity measure) 𝒟 между двумя временными рядами X и Y называется вещественная неотрицательная функция, вычисляющая расстояние между данными рядами: 𝒟(X,Y) ⩾ 0.


Пусть даны два временных ряда X = (x_1, x_2,…,x_m) и Y = (y_1, y_2,…,y_m). Тогда динамической трансформацией временной шкалы (Dynamic Time Warping, DTW) называется мера схожести DTW(X,Y), вычисляемая следующим образом <cit.>:


    DTW(X,Y) = d(m,m), 
    
    d(i,j) = (x_i - y_j)^2 + mind(i-1,j)
    
    d(i, j-1)
    
    d(i-1,j-1)
    , 
    
    d(0,0) = 0, d(i,0) = d(0,j) = ∞, 1⩽ i ⩽ m, 1⩽ j ⩽ m.



Мера DTW позволяет сравнивать ряды, которые смещены вдоль временной оси, сжаты, растянуты или имеют разные длины. Далее нами будет рассматриваться случай, когда временные ряды имеют одинаковые длины, поскольку это упрощает изложение и не ограничивает общность <cit.>.

В формуле (<ref>) матрица (d_ij) ∈ℝ^m × m выражает соответствие между точками сравниваемых временных рядов и называется матрицей трансформации.


Путь трансформации (warping path) P представляет собой последовательность элементов матрицы трансформации, которая определяет соответствие между временными рядами Q и C. Пусть элемент p_t пути трансформации P определяется как p_t=(i,j)_t, тогда 


    P = p_1,p_2,…,p_t,…,p_T,  m⩽ T ⩽ 2n-1.



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


Объем вычислений при подсчете меры схожести DTW может быть уменьшен за счет огрубления схожести. Для этого на путь трансформации могут налагаться дополнительные ограничения. Одним из наиболее часто применяемых ограничений является т.н. полоса Сако—Чиба <cit.>, не позволяющая пути трансформации отклоняться более чем на r элементов от диагонали матрицы трансформации. 

Подпоследовательностью (subsequence) T_i,n временного ряда T называется непрерывное подмножество T, состоящее из n элементов и начинающееся с позиции i: T_i,n = (t_i, t_i+1,…,t_i+n-1),  1 ⩽ i ⩽ m-n+1,  k ≪ m. 



Определим решаемую в данной работе задачу поиска самой похожей подпоследовательности ряда в смысле меры DTW. Пусть имеется длинный временной ряд T и заданный пользователем существенно более короткий временной ряд Q (называемый поисковым запросом (query)), где m=| T |≫| Q | = n. Тогда самой похожей подпоследовательностью (best match) назовем такую подпоследовательность T_i,n, форма которой максимально похожа на запрос Q в смысле меры DTW:


    ∃ i ∀ k   DTW(Q, T_i,n) ⩽ DTW(Q, T_k,n), 1 ⩽ i,k ⩽ m-n+1.



Далее для краткости изложения алгоритма подпоследовательность T_i,n будем называть кандидатом (на максимальную схожесть с запросом Q) и обозначать как C.



 §.§ Последовательный алгоритм.



В настоящее время алгоритм UCR-DTW <cit.> является одним из самых быстрых алгоритмов поиска похожих подпоследовательностей, интегрируя множество техник, ускоряющих поиск. Поскольку предлагаемый параллельный алгоритм использует UCR-DTW в качестве основы, ниже кратко описаны его основные идеи.

Использование квадрата Евклидова расстояния. Евклидово расстояние (ED) между двумя временными рядами Q и C, где | Q | = | C |, определяется следующим образом:

    ED(Q,C)=√(∑_i=1^n(q_i-c_i)^2).


Вычисление квадратного корня в ED (<ref>) и DTW (<ref>) можно опустить, т.к. это не изменит относительного ранжирования подпоследовательностей, похожих на запрос, поскольку обе функции монотонные и вогнутые.

Z-нормализация. Перед вычислением меры схожести запрос и подпоследовательность временного ряда необходимо подвергнуть z-нормализации <cit.>. Z-нормализацией временного ряда T называется временной ряд T̂=(t̂_1,…,t̂_m), элементы которого вычисляются следующим образом:


    t̂_̂î = t_i-μ/σ,  1⩽ i ⩽ m;  μ=1/m∑_i=1^mt_i;   σ^2=1/m∑_i=1^mt_i^2-μ^2.
	

Z-нормализация позволяет сравнивать формы рядов, которые отличны по амплитуде. После нормализации среднее арифметическое временного ряда приблизительно равно 0, а среднеквадратичное отклонение близко к 1. 

Каскадное применение нижних границ схожести. Нижняя граница схожести (lower bound, LB) представляет собой функцию, вычислительная сложность которой меньше вычислительной сложности меры DTW. Нижняя граница используется для отбрасывания кандидатов, заведомо непохожих на запрос, без вычисления меры DTW <cit.>. 

Обозначим за bsf (best-so-far) лучшую текущую нижнюю оценку схожести текущей подпоследовательности T_i,n и поискового запроса Q. Если нижняя граница для кандидата превышает порог bsf, то значение меры DTW для данного кандидата также превысит bsf, и кандидат заведомо непохож на запрос. В процессе сканирования временного ряда алгоритм UCR-DTW пытается улучшить (уменьшить) значение bsf. Значение bsf инициализируется как +∞ и на i-м шаге поиска вычисляется следующим образом: 


    bsf_(i)= min(bsf_(i-1), {[                         +∞ , LB(Q, T_i,n) > bsf_(i-1);              DTW(Q, T_i,n)                , otherwise ]
    ).


Алгоритм UCR-DTW использует следующие нижние границы схожести: LB_KimFL <cit.>, LB_KeoghEC и LB_KeoghEQ <cit.>, которые применяются каскадным образом. 

Нижняя граница LB_KimFL представляет собой Евклидово расстояние между первой и последней парами точек Q и C: 

    LB_KimFL(Q,C) = ED(q̂_1,ĉ_1) + ED(q̂_n, ĉ_n).
	

Нижняя граница LB_KeoghEC показывает минимальную схожесть между оболочкой запроса E (envelope) и кандидатом Ĉ и вычисляется следующим образом:


    LB_KeoghEC(Q, C)=∑_i=1^n
    
    (ĉ_i-u_i)^2,     if ĉ_i>u_i
    
    (ĉ_i-ℓ_i)^2,     if ĉ_i<ℓ_i
    
    0,     otherwise.


В формуле (<ref>) последовательности U = (u_1,..,u_n) и L = (ℓ_1,..,ℓ_n) обозначают верхнюю (upper) и нижнюю (lower) границы оболочки запроса Q, которые вычисляются по формуле (<ref>):


    u_i=max(q̂_i-r,…,q̂_i+r),
    ℓ_i=min(q̂_i-r,…,q̂_i+r).


Нижняя граница LB_KeoghEQ представляет собой Евклидово расстояние между запросом Q и оболочкой кандидата C, то есть по сравнению с LB_KeoghEC роли запроса и кандидата меняются местами:


    LB_KeoghEQ(Q,C) := LB_KeoghEC(C,Q).




Выполнение алгоритма UCR-DTW происходит следующим образом. Сначала нормализуется запрос и вычисляется его оболочка, порог bsf инициализируется значением ∞. Далее алгоритм сканирует временной ряд от начала до конца, применяя к текущей подпоследовательности каскад нижних границ схожести. Если подпоследовательность не была отброшена, то вычисляется расстояние DTW. Далее значение bsf заменяется на вычисленное значение меры DTW, если последнее меньше лучшей текущей оценки схожести. Таким образом, по окончании сканирования ряда алгоритм находит самую похожую подпоследовательность.



§ ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ PHIBESTMATCH.



В данном разделе представлен новый параллельный алгоритм поиска самой похожей подпоследовательности временного ряда PhiBestMatch для кластерных систем с узлами на базе многоядерных процессорных систем архитектуры Intel MIC. Распараллеливание алгоритма основано на следующих основных принципах: параллелизм по данным, выравнивание данных в памяти, векторизация вычислений.

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

Для обработки фрагмента на каждом узле запускается вычислительный процесс, и все процессы используют один и тот же алгоритм. В процессе обработки процессы обмениваются данными для сокращения объема вычислений. Обмены данными между процессами реализуются с помощью технологии MPI <cit.>. Параллелизм по данным внутри узла кластерной системы реализуется следующим образом. В рамках вычислительного процесса на ядрах процессора узла запускаются (fork) нити, разделяющие оперативную память узла. Нити осуществляют параллельную обработку фрагмента на базе технологии OpenMP <cit.>. 

Векторизация циклов является одним из ключевых условий достижения высокой производительности вычислительных программ на параллельных архитектурах <cit.>. Векторизация циклов заключается в преобразовании компилятором последовательности скалярных операторов из тела цикла в один векторный оператор. Таким образом, чтобы повысить производительность поиска похожих подпоследовательностей на Phi KNL, необходимо организовывать вычисления таким образом, чтобы увеличить в алгоритме, насколько это возможно, количество векторизуемых циклов.

Эффективность векторизации циклов, однако, может существенно снижена в силу не выровненного доступа к данным в оперативной памяти, который порождает эффект разделения цикла (loop peeling) <cit.>. Если начальный адрес массива не выровнен на ширину векторного регистра (т.е. количество элементов, которые могут быть загружены в векторный регистр), то компилятор разбивает цикл на три части. Первая часть итераций, которые обращаются к памяти с начального адреса до первого выровненного адреса, и третья часть итераций с последнего выровненного адреса до конечного адреса векторизуется отдельно.

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



 §.§ Фрагментация временного ряда.



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

Пусть N = | T |-n +1 = m-n+1 — количество подпоследовательностей, которые необходимо обработать, F — количество фрагментов, T^(k) — k-й фрагмент временного ряда T=(t_1, t_2, …, t_m), где 0 ⩽ k ⩽ F-1. Тогда фрагмент T^(k) определяется как подпоследовательность T_start,   len, где 


    start = k ·⌊NF⌋ + 1 
    
    len = ⌊NF⌋ + (N mod F) + n - 1,     k=F-1 
    ⌊NF⌋ + n - 1,     otherwise.
 

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



 §.§ Компоновка данных.



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

Выравнивание данных выполняется следующим образом.

Пусть обработка подпоследовательности T_i,n ряда T осуществляется с использованием векторного регистра, вмещающего w вещественных чисел. Если длина подпоследовательности не кратна w, то подпоследовательность дополняется фиктивными нулевыми элементами. Обозначим количество фиктивных элементов за pad = w - (n mod w), тогда выровненная подпоследовательность T̃_i,n определяется следующим образом:


    T̃_i,n = 
    t_i, t_i+1,…,t_i+n-1, 0, 0,…,0_pad,     if  n mod w > 0
     
    t_i, t_i+1,…,t_i+n-1,     otherwise.



Согласно определению (<ref>), ∀ Q, C: DTW(Q,C) = DTW(Q̃,C̃). В дальнейшем изложении мы предполагаем, что запрос и подпоследовательности-кандидаты выровнены в памяти, и для упрощения используем обозначения и Q и C соответственно. Выравнивание подпоследовательностей позволяет избежать накладных расходов на разбивание циклов на несколько итераций (loop peeling).

Все (выровненные) подпоследовательности временного ряда сохраняются в виде матрицы подпоследовательностей для обеспечения векторизации вычислений.

Матрица подпоследовательностей S_T^n ∈ℝ^N×(n+pad) определяется следующим образом:


    S_T^n(i,j):=t̃_i+j-1.
 



Обозначим количество нижних границ схожести, используемых в алгоритме поиска самой похожей подпоследовательности, за lb_max (lb_max⩾ 1), и за LB_1, LB_2, …, LB_lb_max обозначим эти границы, перечисленные в порядке их вычисления в каскаде применения нижних оценок. Тогда L_T^n ∈ℝ^N × lb_max, матрица нижних границ схожести всех подпоследовательностей длины n временного ряда T с поисковым запросом Q, имеет следующий вид:


    L_T^n(i,j):=LB_j(T_i,n, Q).




Карта схожести представляет собой вектор-столбец B_T^n ∈𝔹^N, который для каждой подпоследовательности длины n ряда T хранит конъюнкцию применения всех нижних границ схожести этой подпоследовательности с текущим значением порога (bsf):


    B_T^n(i):=⋀_j=1^lb_max (L_T^n(i,j) < bsf).
 


Введем матрицу для хранения подпоследовательностей-кандидатов, то есть тех подпоследовательностей из матрицы S_T^n, которые не были отброшены при применении нижней границы схожести. Данная матрица будет обрабатываться параллельно для вычисления значений меры DTW между кандидатами и запросом. Далее минимальное из вычисленных значение меры схожести DTW используется в качестве порога bsf.

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

Обозначим число нитей, используемых параллельным алгоритмом, за p (p⩾1). Введем размер сегмента s∈ℕ (s ⩽⌈N/p⌉) — количество строк матрицы кандидатов, обрабатываемых одной нитью. Тогда матрица кандидатов C_T ^ n ∈ℝ ^ (s· p) × (n + pad) определяется следующим образом:


    C_T^n(i,·):=S_T^n(k,·): B_T^n(i)=.
 


Несмотря на то, что вычисление матрицы кандидатов является распараллеливаемой операцией, векторизация операторов соответствующего цикла затруднена, поскольку в соответствии с определением (<ref>) при вычислении меры схожести DTW имеет место зависимость по данным.



 §.§ Вычислительная схема алгоритма.



Предлагаемая в данной работе параллельная реализация поиска похожих подпоследовательностей временного ряда представлена в алг. <ref> и на рис. <ref>. 
Алгоритм выполняется следующим образом. 

	

Для инициализации вычислений алгоритм определяет номер текущего вычислительного процесса myrank с помощью функций библиотеки MPI. В дальнейшем в рамках алгоритма каждый вычислительный процесс с номером myrank обрабатывает матрицу подпоследовательностей S^n_T^(myrank) фрагмента T^(myrank) исходного временного ряда T. Переменная bsf инициализируется значением меры схожести DTW между поисковым запросом и случайной подпоследовательностью из данного фрагмента временного ряда. 



Далее выполняется подготовка данных для вычислений. Формируется матрица выровненных подпоследовательностей. В цикле, выполняемом по строкам матрицы подпоследовательностей, каждая подпоследовательность z-нормализуется, и выполняются вычисление строк матрицы нижних границ схожести. Следует обратить внимание, что, строго говоря, предварительное вычисление матрицы нижних границ схожести представляет собой избыточные накладные расходы. В этом состоит ключевое отличие от алгоритма UCR-DTW, где нижние границы схожести вычисляются каскадом (следующая нижняя граница вычисляется только в случае, если с помощью предыдущей границы  не выявлено, что текущий кандидат не является заведомо непохожим на запрос). Тем не менее, эти предвычисления обоснованы тем, что выполняются однократно, и могут быть реализованы с помощью распараллеливаемых и векторизуемых компилятором циклов.

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

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

Для распараллеливания описанных выше операций матрица нижних границ схожести разбивается на сегменты (по числу используемых алгоритмом нитей). По построению матрица нижних границ схожести имеет существенно большее количество строк, чем матрица кандидатов. Соответственно, после заполнения и вычисления матрицы кандидатов, и обновления порога bsf каждая нить должна продолжить сканирование своего сегмента вплоть до его исчерпания. Для хранения номера последнего обработанного кандидата в сегменте вводится индексный массив Pos ∈ℕ^p, где 




    pos_i :=  k: p · (i-1)+1 ⩽ k ⩽⌈Ni · p⌉ ∧ ∀ j, 1 ⩽ j ⩽ lb_max, LB_T^n(k,j)<bsf.
 

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

Чтобы вывести индекс самой похожей на запрос подпоследовательности, введен индексный массив Idx ∈ℕ^s · p, который предназначен для хранения позиции подпоследовательности во временном ряде и определяется следующим образом:


    idx_i := k: 1 ⩽ k ⩽ N ∧∃ S_T^n(i,·) ⇔∃ T_i,n⇔ k = (i-1) · n+1.
 

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

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

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



§ ВЫЧИСЛИТЕЛЬНЫЕ ЭКСПЕРИМЕНТЫ.





 §.§ Аппаратная платформа.



Для исследования эффективности разработанного алгоритма были проведены вычислительные эксперименты. В качестве аппаратной платформы экспериментов использовались узлы кластерной системы <<Торнадо ЮУрГУ>> <cit.>, характеристики которых приведены в tab:HardwarePhiBestMatch.



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

Исследование на одном вычислительном узле кластерной системы позволяет определить, насколько эффективно реализованы параллельные вычисления, выполняемые многоядерной процессорной системой Intel Xeon Phi. Для такого исследования использовалась упрощенная версия алг. <ref> <cit.>, в которой ряд отождествлен с одним фрагментом, вызовы коммуникационной библиотеки MPI не задействуются, и цикл обработки подпоследовательностей ряда продолжается до исчерпания подпоследовательностей-кандидатов. 

Во всех экспериментах размер сегмента матрицы кандидатов — количество строк матрицы кандидатов, обрабатываемых одной нитью в рамках одного вычислительного узла кластерной системы (см. формулу <ref>) — имеет значение s=100.



 §.§ Эффективность алгоритма на одном узле кластера.


Исследование эффективности алгоритма на одном вычислительном узле кластера <<Торнадо ЮУрГУ>> производилось на наборах данных, которые представлены в tab:DatasetsSingleNodePhiBestMatch.



Временной ряд Random Walk получен искусственно на основе модели случайных блужданий <cit.>. Проведение экспериментов по исследованию свойств алгоритмов поиска подпоследовательностей на базе меры схожести DTW с использованием подобных временных рядов является общепринятой практикой (см., например, работы <cit.>). 

Временной ряд EPG (Electrical Penetration Graph, график электрического проникновения) использован в экспериментах в работе <cit.>. EPG представляет собой набор сигналов, используемых энтомологами для сравнения поведения исследуемого насекомого с поведением цикадок macrosteles quadrilineatus, которые являются переносчиками болезней растений и ежегодно наносят сельскому хозяйству США ущерб более чем на 2 млн. долларов <cit.>.

В экспериментах исследовались производительность и масштабируемость алгоритма PhiBestMatch. Под производительностью понимается время работы алгоритма без учета времени загрузки данных в память и выдачи результата. Масштабируемость параллельного алгоритма означает его способность адекватно адаптироваться к увеличению параллельно работающих вычислительных элементов (процессов, процессоров, нитей и др.) и характеризуется ускорением и параллельной эффективностью, которые определяются следующим образом <cit.>. Ускорение и параллельная эффективность параллельного алгоритма, запускаемого на k нитях, вычисляются как s(k)=t_1t_k и e(k)=s(k)k
соответственно, где t_1 и t_k — время работы алгоритма на одной и k
нитях соответственно.

В экспериментах рассматривались вышеупомянутые показатели в зависимости от изменения параметра r (ширина полосы Сако—Чиба), значения которого брались в долях длины поискового запроса n. 

Время работы алгоритма PhiBestMatch на многоядерных платформах сравнивалось с временем работы  алгоритма UCR-DTW <cit.>.

Результаты экспериментов по исследованию производительности алгоритма представлены на рис. <ref>. Можно видеть, что PhiBestMatch работает до 5 раз быстрее, чем алгоритм UCR-DTW. Вместе с тем видно, на производительность алгоритма PhiBestMatch на платформах двухпроцессорного узла Intel Xeon и многоядерной системы Intel Xeon Phi влияют следующие два параметра: ширина полосы Сако—Чиба r и длина поискового запроса n.



При малых значениях этих параметров (примерно 0 < r ⩽ 0.5n и n<512) алгоритм PhiBestMatch на платформе двухпроцессорного узла Intel Xeon работает несколько быстрее или примерно с тем же быстродействием, что и на платформе многоядерной системы Intel Xeon Phi. При бо́льших значениях данных параметров (0.5n<r ⩽ n и n ⩾ 512) алгоритм работает быстрее на платформе Intel Xeon Phi. Это означает, что заложенные в алгоритм при проектировании возможности векторизации наилучшим образом проявляются при увеличении общего объема вычислений. 

Поисковые запросы с длиной n ⩾ 512, равно как и значение параметра r=1 требуются на практике в ряде приложений, требующих при определении схожести подпоследовательностей как можно более высокую точность, например, в медицине при исследовании ЭКГ <cit.>, в 
энтомологии <cit.>, в астрономии <cit.> и др.  





Результаты экспериментов по исследованию масштабируемости алгоритма представлены на fig:PhiBestMatch-SingleNode-Scalability-RandomWalk и <ref>. Можно видеть, что PhiBestMatch демонстрирует близкое к линейному ускорение и параллельную эффективность, близкую к 100%, если количество нитей, на которых запущен алгоритм, совпадает с количеством физических ядер системы Intel Xeon Phi.

При увеличении количества нитей, запускаемых на одном физическом ядре системы, ускорение становится сублинейным, равно как наблюдается и падение параллельной эффективности. При этом наилучшие показатели ускорения и параллельной эффективности ожидаемо наблюдаются при значениях параметра r 0.8 и 1 от длины поискового запроса n, обеспечивающих алгоритму наибольшую вычислительную нагрузку. Например, при задействовании 240 нитей при r=0.8n обработка ряда Random Walk выполняется с ускорением 120 и параллельной эффективностью 50%; при r=n обработка ряда EPG — с ускорением 130 и параллельной эффективностью 52%. 

Полученные результаты позволяют сделать заключение о хорошей масштабируемости разработанного алгоритма и эффективном использовании им возможностей векторизации вычислений на многоядерной системе Intel Xeon Phi в рамках одного вычислительного узла кластерной системы. Указанные свойства проявляются при значениях параметров r (ширина полосы Сако—Чиба) и n (длина поискового запроса), обеспечивающих алгоритму наибольшую вычислительную нагрузку: 0.8n ⩽ r ⩽ n и n ⩾ 512 соответственно. 



 §.§ Эффективность алгоритма на кластерной системе.




В исследовании эффективности алгоритма на кластерной системе в целом было задействовано от 16 до 128 вычислительных узлов суперкомпьютера <<Торнадо ЮУрГУ>> (см. tab:HardwarePhiBestMatch). Исследование производилось на наборах данных, которые представлены в tab:DatasetsAllNodesPhiBestMatch.



В экспериментах исследовалось ускорение масштабируемости (scaled speedup) параллельного алгоритма, которое определяется как ускорение, демонстрируемое алгоритмом при линейном увеличении объема данных и количества используемых вычислительных узлов <cit.> и вычисляется следующим образом: s_scaled=p · mt_p(p · m), где p — количество задействованных вычислительных узлов, m — объем исходных данных, t_p(p · m) — время выполнения алгоритма на p узлах при обработке исходных данных, имеющих объем p · m. 

При этом значение параметра r было зафиксировано как r=n и варьировалась длина поискового запроса.

Результаты экспериментов по исследованию ускорения масштабируемости алгоритма представлены на fig:PhiBestMatch-AlleNodes-ScaledSpeedup-RandomWalk и <ref>.



Можно видеть, что PhiBestMatch демонстрирует  ускорение масштабируемости, близкое к линейному, как для синтетических, так и для реальных данных. При этом поиск подпоследовательности большей длины показывает более высокое ускорение масштабируемости, поскольку это обеспечивает больший объем вычислений в рамках одного вычислительного узла кластерной системы.



Полученные результаты позволяют сделать заключение о хорошей масштабируемости разработанного алгоритма при работе на кластерной системе с вычислительными узлами на базе многоядерных процессоров Intel Xeon Phi. Указанные свойства проявляются при значениях параметров r (ширина полосы Сако—Чиба) и n (длина поискового запроса), обеспечивающих алгоритму наибольшую вычислительную нагрузку: 0.8n ⩽ r ⩽ n и n ⩾ 512 соответственно.



§ ЗАКЛЮЧЕНИЕ.


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

Предложен новый параллельный алгоритм для поиска похожих подпоследовательностей временного ряда на кластерной системе с вычислительными узлами на базе многоядерных процессоров Intel Xeon Phi поколения Knights Landing, названный PhiBestMatch. Алгоритм предполагает двухуровневое распараллеливание вычислений: на уровне всех узлов кластерной системы используется технология MPI, в рамках одного вычислительного узла кластера используется технология OpenMP. Для эффективного использования возможностей векторизации вычислений процессоров Phi KNL в рамках одного вычислительного узла используются дополнительные матричные структуры данных и избыточные вычисления. Проведены вычислительные эксперименты, исследующие быстродействие и масштабируемость алгоритма на синтетических и реальных наборах данных как на кластерной системе в целом, так и в рамках одного вычислительного узла кластера. Результаты  экспериментов показали хорошую масштабируемость алгоритма PhiBestMatch на одном узле и на кластерной системе в целом.



Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 17-07-00463), Правительства РФ в соответствии с Постановлением № 211 от 16.03.2013 (соглашение № 02.A03.21.0011) и Министерства образования и науки РФ (государственное задание 2.7905.2017/8.9).


	
	
	
	
	VoevodinV02
	Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. СПб.:
	БХВ-Петербург, 2002.
	
	
	
	
	
	
	
	DBLP:conf/sigmod/AthitsosPPKG08
	Athitsos V., Papapetrou P., Potamias M. et al. Approximate embedding-based subsequence matching of time series // Proc. of the ACM SIGMOD Int. Conf. on Management of Data (Vancouver, Canada, June 10–12, 2008). New York: ACM, 2008. 365–378.  DOI: https://dx.doi.org/10.1145/1376616.137665610.1145/1376616.1376656.
	
	bacon
	Bacon D.F., Graham S.L., Sharp O.J. Compiler transformations for high-performance computing // ACM Comput. Surv. 1994. 26, N 4. 345–420. DOI: https://dx.doi.org/10.1145/197405.19740610.1145/197405.197406.
	
	DBLP:conf/kdd/BerndtC94
	Berndt D.J., Clifford J. Using dynamic time warping to find patterns in time series // Proc. of the 1994 AAAI Workshop on Knowledge Discovery in Databases (Seattle, Washington, July, 1994). AAAI Press, 1994. 359–370.
	
	chrysos
	Chrysos G. Intel Xeon Phi Coprocessor (codename Knights Corner) // 2012 IEEE Hot Chips 24th Symposium (HCS) (Cupertino, CA, USA, August 27–29, 2012). IEEE, 2012. 1–31. DOI: https://dx.doi.org/10.1109/HOTCHIPS.2012.747648710.1109/HOTCHIPS.2012.7476487.
	
	DBLP:journals/pvldb/DingTSWK08
	Ding H., Trajcevski G., Scheuermann P., Wang X., Keogh E. Querying and mining of time series data: experimental comparison of representations and distance measures // Proc. of the VLDB Endowment. 2008. 1, N 2. 1542–1552.
	
	fu
	Fu A.W., Keogh E.J., Lau L.Y.H., Ratanamahatana C.A., Wong R.C. Scaling and time warping in time series querying // VLDB. 2008. 17, N 4. 899–921. DOI: https://doi.org/10.1007/s00778-006-0040-z10.1007/s00778-006-0040-z.
	
	DBLP:conf/pvm/Gropp12
	Gropp W. MPI 3 and beyond: why MPI is successful and what challenges it faces // Proc. of the 19th European MPI Users' Group Meeting. Lecture Notes in Computer Science. Vol. 7490. Berlin: Springer, 2012. 1–9. DOI: https://dx.doi.org/10.1007/978-3-642-33518-1_110.1007/978-3-642-33518-1_1.
	
	keogh
	Keogh E.J., Ratanamahatana C.A. Exact indexing of dynamic time warping // Knowl. Inf. Syst. 2005. 7, N 3. 358–386. DOI: https://doi.org/10.1007/s10115-004-0154-910.1007/s10115-004-0154-9.
	
	
	
	
	KostenetskyS16
	Kostenetskiy P.S., Safonov A.Y. SUSU supercomputer resources // Proc. of the 10th Annual Int. Scientific Conf. on Parallel Computing Technologies (PCT 2016). CEUR Workshop Proceedings. Vol. 1576. 2016. 561–573.
	
	KraevaZ18
	Kraeva Ya., Zymbler M. An Efficient Subsequence Similarity Search on Modern Intel Many-core Processors for Data Intensive Applications // Аналитика и управление данными в областях с интенсивным использованием данных: Сборник научных трудов XX Международной конференции DAMDID / RCDL'2018 (9–12 октября 2018 г., Москва, Россия). С. 116–124.
	
	DBLP:books/bc/KumarGGK94
	Kumar V., Grama A., Gupta A., Karypis G. Introduction to parallel
	computing. Benjamin/Cummings, 1994.
	
	lim
	Lim S., Park H., Kim S. Using multiple indexes for efficient subsequence matching in time-series databases // Proc. of the 11th Int. Conf. on Database Systems for Advanced Applications. Lecture Notes in Computer Science. Vol. 3882. Berlin: Springer, 2006. 65–79. DOI: https://doi.org/10.1007/11733836_710.1007/11733836_7.
	
	DBLP:conf/sc/Mattson06
	Mattson T. Introduction to OpenMP // Proc. of the 2006 ACM/IEEE Conf. on Supercomputing (Tampa, FL, USA, November 11–17, 2006). New York: ACM Press, 2006. 209.
	DOI: https://dx.doi.org/10.1145/1188455.118867310.1145/1188455.1188673.
	
	DBLP:journals/cmpb/MazandaraniM18
	Mazandarani F.N., Mohebbi M. Wide complex tachycardia discrimination using dynamic time warping of ECG beats // Computer Methods and Programs in Biomedicine. 2018. 164. 238–249. DOI: https://dx.doi.org/10.1016/j.cmpb.2018.04.00910.1016/j.cmpb.2018.04.009.
	
	MovchanZ16
	Movchan A.V., Zymbler M.L. Parallel implementation of searching the most similar subsequence in time series for computer systems with distributed memory // Proc. of the 10th Annual Int. Scientific Conf. on Parallel Computing Technologies (PCT 2016). CEUR Workshop Proceedings. Vol. 1576. 2016. 615–628. 
	
	Pearson1905
	Pearson K. The problem of the random walk // Nature. 1905. 72, N 1865. 294. DOI: https://dx.doi.org/10.1038/072342a010.1038/072342a0.
	
	rakthanmanon
	Rakthanmanon T., Campana B.J.L., Mueen A., Batista G.E.A.P.A., Westover M.B., et al. Searching and mining trillions of time series subsequences under dynamic time warping // Proc. of the 18th ACM SIGKDD Int. Conf. on Knowledge Discovery and Data Mining (Beijing, China, August 12–16, 2012). New York: ACM, 2012. 262–270. DOI: http://doi.acm.org/10.1145/2339530.233957610.1145/2339530.2339576.
	
	Rebbapragada2009	
	Rebbapragada U., Protopapas P., Brodley C.E., Alcock C. Finding anomalous periodic time series // Machine Learning. 2009. 74, N 3. 281–313. DOI: https://dx.doi.org/10.1007/s10994-008-5093-310.1007/s10994-008-5093-3.	
	
	sakoe
	Sakoe H., Chiba S. Readings in speech recognition. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc., 1990. 159–165.
	
	sakurai
	Sakurai Y., Faloutsos C., Yamamuro M. Stream monitoring under the time warping distance // Proc. of the 23rd Int. Conf. on Data Engineering (Istanbul, Turkey, April 15–20, 2007). Washington, DC, USA: IEEE Computer Society, 2007. 1046–1055. DOI: https://dx.doi.org/10.1109/ICDE.2007.36896310.1109/ICDE.2007.368963.	
	
	sart
	Sart D., Mueen A., Najjar W.A., Keogh E.J., Niennattrakul V. Accelerating dynamic time warping subsequence search with GPUs and FPGAs // Proc. of the 2010 IEEE Int. Conf. on Data Mining (Sydney, Australia, December 14–17, 2010). Washington, DC, USA: IEEE Computer Society, 2010. 1001–1006. DOI:https://doi.org/10.1109/ICDM.2010.2110.1109/ICDM.2010.21.
	
	sodani
	Sodani A. Knights Landing (KNL): 2nd generation Intel Xeon Phi processor // 2015 IEEE Hot Chips 27th Symposium (HCS) (Cupertino, CA, USA, August 22–25, 2015). IEEE, 2015. 1–24.	
	
	SokolinskayaS16
	Sokolinskaya I.,  Sokolinsky L. Revised pursuit algorithm for solving non-stationary linear programming problems on modern computing clusters with manycore accelerators // Communications in Computer and Information Science. 2016. 687. 212–223. DOI: https://dx.doi.org/10.1007/978-3-319-55669-7_1710.1007/978-3-319-55669-7_17.
	
	DBLP:conf/icacci/ShabibNNDPSATS15
	Shabib A., Narang A., Niddodi C.P. et al. Parallelization of searching and mining time series data using Dynamic Time Warping // 2015 Int. Conf. on Advances in Computing, Communications and Informatics (Kochi, India, August 10–13, 2015). IEEE, 2015. 343–348. DOI: https://dx.doi.org/10.1109/ICACCI.2015.727563310.1109/ICACCI.2015.7275633.
	
	srikanthan
	Srikanthan S., Kumar A., Gupta R. Implementing the dynamic time warping algorithm in multithreaded environments for real time and unsupervised pattern discovery // 2011 2nd Int. Conf. on Computer and Communication Technology (Allahabad, India, September 15–17, 2011). IEEE, 2015. 394–398. DOI: https://doi.org/10.1109/ICCCT.2011.607511110.1109/ICCCT.2011.6075111.
	
	takahashi
	Takahashi N., Yoshihisa T., Sakurai Y., Kanazawa M. A parallelized data stream processing system using dynamic time warping distance // 2009 Int. Conf. on Complex, Intelligent and Software Intensive Systems (Fukuoka, Japan, March 16–19, 2009). IEEE, 2009. 1100–1105. DOI: https://doi.org/10.1109/CISIS.2009.7710.1109/CISIS.2009.77.
	
	tarango
	Tarango J., Keogh E.J., Brisk P. Instruction set extensions for dynamic time warping // Proc. of the Int. Conf. on Hardware/Software Codesign and System Synthesis (Montreal, QC, Canada, September 29–October 4, 2013). IEEE, 2013. 18:1–18:10. DOI: https://doi.org/10.1109/CODES-ISSS.2013.665900510.1109/CODES-ISSS.2013.6659005.
	
	wang
	Wang Z., Huang S., Wang L., Li H., Wang Y., et al. Accelerating subsequence similarity search based on dynamic time warping distance with FPGA // Proc. of the ACM/SIGDA Int. Symposium on Field Programmable Gate Arrays (Monterey, CA, USA, February 11–13, 2013). New York: ACM, 2013. 53–62. DOI: http://doi.acm.org/10.1145/2435264.243527710.1145/2435264.2435277.
	
	zhang
	Zhang Y., Adl K., Glass J.R. Fast spoken query detection using
	lower-bound dynamic time warping on graphical processing units // 2012 IEEE Int. Conf. on Acoustics, Speech and Signal Proc. (Kyoto, Japan, March 25–30, 2012). IEEE, 2012. 5173–5176. DOI: https://doi.org/10.1109/ICASSP.2012.628908510.1109/ICASSP.2012.6289085.
	
	zymbler
	Zymbler M. Best-match time series subsequence search on the Intel Many Integrated Core architecture // Proc. of the 19th East European Conf. on Advances in Databases and Information Systems (Poitiers, France, September 8–11, 2015). Lecture Notes in Computer Science. 9282. Heidelberg: Springer, 2015. 275–286. DOI: https://doi.org/10.1007/978-3-319-23135-8_1910.1007/978-3-319-23135-8_19.
	
	


[t]3.7cmПоступила в редакцию
04.12.2018





 








	^1 South Ural State University, Department of System Programming; prospekt Lenina 76, Chelyabinsk, 454080, Russia; Master student, e-mail:  mailto:kraevaya@susu.rukraevaya@susu.ru

	^2 South Ural State University, Department of System Programming; prospekt Lenina 76, Chelyabinsk, 454080, Russia; Cand. Sci., Associate Professor, e-mail: mailto:mzym@susu.rumzym@susu.ru





December 4, 2018

The Use of MPI and OpenMP Technologies for Subsequence Similarity Search 
 in Very Large Time Series on Computer Cluster System 
with Nodes Based on the Intel Xeon Phi Knights Landing Many-core Processor
    Ya.A. Kraeva^1 and M.L. Zymbler^2
    Received: date / Accepted: date
===========================================================================================================================================================================================================




	Nowadays, subsequence similarity search is required in a wide range of time series mining applications: climate modeling, financial forecasts, medical research, etc. In most of these applications, the Dynamic Time Warping (DTW) similarity measure is used since DTW is empirically confirmed as one of the best similarity measure for most subject domains. Since the DTW measure has a quadratic computational complexity w.r.t. the length of query subsequence, a number of parallel algorithms for various many-core architectures have been developed, namely FPGA, GPU, and Intel MIC. In this article, we propose a new parallel algorithm for subsequence similarity search in very large time series on computer cluster systems with nodes based on Intel Xeon Phi Knights Landing (KNL) many-core processors. Computations are parallelized on two levels as follows: through MPI at the level of all cluster nodes, and through OpenMP within one cluster node. The algorithm involves additional data structures and redundant computations, which make it possible to effectively use the capabilities of vector computations on Phi KNL. Experimental evaluation of the algorithm on real-world and synthetic datasets shows that the proposed algorithm is highly scalable.
	
	time series, similarity search, Dynamic Time Warping, parallel algorithm, OpenMP, Intel Xeon Phi, Knights Landing, data layout, vectorization.






	
	
	
	
	VoevodinV02en
	V. V. Voevodin, and Vl. V. Voevodin, The Parallel Computing. (BHV-Petersburg, St. Petersburg, 2002).
	
	
	
	
	
	
	
	DBLP:conf/sigmod/AthitsosPPKG08en
	V. Athitsos, P. Papapetrou, M. Potamias, et al., “Approximate Embedding-Based Subsequence Matching of Time Series,” in Proc. of the ACM SIGMOD Int. Conf. on Management of Data. Vancouver, Canada, June 10–12, 2008 (ACM, New York, 2008), pp. 365–378.  doi: https://dx.doi.org/10.1145/1376616.137665610.1145/1376616.1376656.
	
	baconen
	D. F. Bacon, S. L. Graham, and O. J. Sharp, “Compiler transformations for high-performance computing,” ACM Comput. Surv. 26 (4), 345–420 (1994).  doi: https://dx.doi.org/10.1145/197405.19740610.1145/197405.197406.
	
	DBLP:conf/kdd/BerndtC94en
	D. J. Berndt, and J. Clifford, “Using Dynamic Time Warping to Find Patterns in Time Series,” in Proc. of the 1994 AAAI Workshop on Knowledge Discovery in Databases, Seattle, Washington, July, 1994 (AAAI Press, 1994), pp. 359–370.
	
	chrysosen
	G. Chrysos, “Intel Xeon Phi Coprocessor (Codename Knights Corner),” in 2012 IEEE Hot Chips 24th Symposium (HCS), Cupertino, CA, USA, August 27–29, 2012 (IEEE, 2012), pp. 1–31. doi: https://dx.doi.org/10.1109/HOTCHIPS.2012.747648710.1109/HOTCHIPS.2012.7476487.
	
	DBLP:journals/pvldb/DingTSWK08en
	H. Ding, G. Trajcevski, P. Scheuermann, X. Wang, and E. Keogh, “Querying and Mining of Time Series Data: Experimental Comparison of Representations and Distance Measures,” Proc. of the VLDB Endowment. 1 (2), 1542–1552 (2008).
	
	fuen
	A. W. Fu, E. J. Keogh, L. Y. H. Lau, C. A. Ratanamahatana, and R. C. Wong, “Scaling and Time Warping in Time Series Querying,” VLDB. 17 (4), 899–921 (2008).  doi: https://doi.org/10.1007/s00778-006-0040-z10.1007/s00778-006-0040-z.
	
	DBLP:conf/pvm/Gropp12en
	W. Gropp, “MPI 3 and Beyond: Why MPI is Successful and What Challenges It Faces,” in Proc. of the 19th European MPI Users' Group Meeting. Lecture Notes in Computer Science (Springer, Berlin, 2012), Vol. 7490, pp. 1–9. doi: https://dx.doi.org/10.1007/978-3-642-33518-1_110.1007/978-3-642-33518-1_1.
	
	keoghen
	E. J. Keogh, and C. A. Ratanamahatana, “Exact Indexing of Dynamic Time Warping,” Knowl. Inf. Syst. 7 (3). pp. 358–386. (2005). doi: https://doi.org/10.1007/s10115-004-0154-910.1007/s10115-004-0154-9.
	
	
	
	
	KostenetskyS16en
	P. S. Kostenetskiy, and A. Y. Safonov, “SUSU Supercomputer Resources,” in Proc. of the 10th Int. Conf. on Parallel Computing Technologies (PCT 2016), Arkhangelsk, Russia, March 29–31, 2016. CEUR Workshop Proc. Vol. 1576, pp. 561–573.
	
	KraevaZ18en
	Ya. Kraeva, and M. Zymbler “An Efficient Subsequence Similarity Search on Modern Intel Many-core Processors for Data Intensive Applications,” in Proc. of the 20th Int. Conf. on Data Analytics and Management in Data Intensive Domains, DAMDID/RCDL 2018, Moscow, Russia, October 9–12, 2018, pp. 116–124.
	
	DBLP:books/bc/KumarGGK94en
	V. Kumar, A. Grama, A. Gupta, and G. Karypis. Introduction to Parallel Computing (Benjamin/Cummings, 1994).
	
	limen
	S. Lim, H. Park, and S. Kim, “Using Multiple Indexes for Efficient Subsequence Matching in Time-Series Databases,” in Proc. of the 11th Int. Conf. on Database Systems for Advanced Applications. Lecture Notes in Computer Science (Springer, Berlin, 2006), Vol. 3882, pp. 65–79. doi: https://doi.org/10.1007/11733836_710.1007/11733836_7.
	
	DBLP:conf/sc/Mattson06en
	T. Mattson, “Introduction to OpenMP,” in Proc. of the 2006 ACM/IEEE Conf. on Supercomputing, Tampa, FL, USA, November 11–17, 2006 (ACM Press, New York, 2006), pp. 209.
	doi: https://dx.doi.org/10.1145/1188455.118867310.1145/1188455.1188673.
	
	DBLP:journals/cmpb/MazandaraniM18en
	F. N. Mazandarani, and M. Mohebbi, “Wide Complex Tachycardia Discrimination Using Dynamic Time Warping of ECG Beats,” Computer Methods and Programs in Biomedicine. 164, 238–249 (2018). doi: https://dx.doi.org/10.1016/j.cmpb.2018.04.00910.1016/j.cmpb.2018.04.009.
	
	MovchanZ16en
	A. V. Movchan, and M. L. Zymbler, “Parallel Implementation of Searching the Most Similar Subsequence in Time Series for Computer Systems with Distributed Memory,” in Proc. of the 10th Annual International Scientific Conference on Parallel Computing Technologies (PCT 2016), Arkhangelsk, Russia, March 29–31, 2016. CEUR Workshop Proc. Vol. 1576, pp. 615–628. 
	
	Pearson1905en
	K. Pearson, “The Problem of the Random Walk,” Nature. 72 (1865). 294. (1905). doi: https://dx.doi.org/10.1038/072342a010.1038/072342a0.
	
	rakthanmanonen
	T. Rakthanmanon, B. J. L. Campana, A. Mueen, G. E. A. P. A. Batista, M. B. Westover, et al., “Searching and Mining Trillions of Time Series subsequences under Dynamic Time Warping,” in Proc. of the 18th ACM SIGKDD Int. Conf. on Knowledge Discovery and Data Mining, Beijing, China, August 12–16, 2012 (ACM, New York, 2012), pp. 262–270. doi: http://doi.acm.org/10.1145/2339530.233957610.1145/2339530.2339576.
	
	Rebbapragada2009en	
	U. Rebbapragada, P. Protopapas, C. E. Brodley, and C. Alcock, “Finding Anomalous Periodic Time Series,” Machine Learning 74 (3), 281–313 (2009). doi: https://dx.doi.org/10.1007/s10994-008-5093-310.1007/s10994-008-5093-3.	
	
	sakoeen
	H. Sakoe, and S. Chiba, Readings in Speech Recognition. (San Francisco, CA, USA, Morgan Kaufmann Publishers Inc., 1990), pp. 159–165.
	
	sakuraien
	Y. Sakurai, C. Faloutsos, and M. Yamamuro, “Stream Monitoring under the Time Warping Distance,” in Proc. of the 23rd Int. Conf. on Data Engineering, Istanbul, Turkey, April 15–20, 2007 (IEEE Computer Society, Washington, DC, USA, 2007). pp. 1046–1055. doi: https://dx.doi.org/10.1109/ICDE.2007.36896310.1109/ICDE.2007.368963.	
	
	sarten
	D. Sart, A. Mueen, W. A. Najjar, E. J. Keogh, and V. Niennattrakul, “Accelerating Dynamic Time Warping Subsequence Search with GPUs and FPGAs,” in Proc. of the 2010 IEEE Int. Conf. on Data Mining, Sydney, Australia, December 14–17, 2010 (IEEE Computer Society, Washington, DC, USA, 2010), pp. 1001–1006. doi: https://doi.org/10.1109/ICDM.2010.2110.1109/ICDM.2010.21.
	
	sodani
	A. Sodani, “Knights Landing (KNL): 2nd generation Intel Xeon Phi processor,” in 2015 IEEE Hot Chips 27th Symposium (HCS), Cupertino, CA, USA, August 22–25, 2015 (IEEE, 2015), pp. 1–24.	
	
	SokolinskayaS16en
	I. Sokolinskaya, and L. Sokolinsky, “Revised Pursuit Algorithm for Solving Non-Stationary Linear Programming Problems on Modern Computing Clusters with Manycore Accelerators,” Communications in Computer and Information Science. 687, 212–223 (2016). doi: https://dx.doi.org/10.1007/978-3-319-55669-7_1710.1007/978-3-319-55669-7_17.
	
	DBLP:conf/icacci/ShabibNNDPSATS15en
	A. Shabib, A. Narang, C. P. Niddodi, et al., “Parallelization of Searching and Mining Time Series Data using Dynamic Time Warping,” in 2015 Int. Conf. on Advances in Computing, Communications and Informatics, Kochi, India, August 10–13, 2015 (IEEE, 2015), pp. 343–348. doi: https://dx.doi.org/10.1109/ICACCI.2015.727563310.1109/ICACCI.2015.7275633.
	
	srikanthanen
	S. Srikanthan, A. Kumar, and R. Gupta, “Implementing the Dynamic Time Warping Algorithm in Multithreaded Environments for Real Time and Unsupervised Pattern Discovery,” in 2011 2nd International Conference on Computer and Communication Technology, Allahabad, India, September 15–17, 2011 (IEEE, 2015), pp. 394–398. doi: https://doi.org/10.1109/ICCCT.2011.607511110.1109/ICCCT.2011.6075111.
	
	takahashien
	N. Takahashi, T. Yoshihisa, Y. Sakurai, and M. Kanazawa, “A Parallelized Data Stream Processing System using Dynamic Time Warping Distance,” in 2009 Int. Conf. on Complex, Intelligent and Software Intensive Systems, Fukuoka, Japan, March 16–19, 2009 (IEEE, 2009), pp. 1100–1105. doi: https://doi.org/10.1109/CISIS.2009.7710.1109/CISIS.2009.77.
	
	tarango
	J. Tarango, E. J. Keogh, and P. Brisk, “Instruction Set Extensions for Dynamic Time Warping” in Proc. of the Int. Conf. on Hardware/Software Codesign and System Synthesis, Montreal, QC, Canada, September 29–October 4, 2013 (IEEE, 2013), pp. 18:1–18:10. doi: https://doi.org/10.1109/CODES-ISSS.2013.665900510.1109/CODES-ISSS.2013.6659005.
	
	wang
	Z. Wang, S. Huang, L. Wang, H. Li, Y. Wang, et al., “Accelerating Subsequence Similarity Search Based on Dynamic Time Warping Distance with FPGA” in Proc. of the ACM/SIGDA International Symposium on Field Programmable Gate Arrays, Monterey, CA, USA, February 11–13, 2013 (ACM, New York, 2013), pp. 53–62. doi: http://doi.acm.org/10.1145/2435264.243527710.1145/2435264.2435277.
	
	zhang
	Y. Zhang, K. Adl, and J. R. Glass, “Fast Spoken Query Detection using
	Lower-Bound Dynamic Time Warping on Graphical Processing Units,” in 2012 IEEE Int. Conf. on Acoustics, Speech and Signal Proc., Kyoto, Japan, March 25–30, 2012 (IEEE, 2012), pp. 5173–5176. doi: https://doi.org/10.1109/ICASSP.2012.628908510.1109/ICASSP.2012.6289085.
	
	zymbler
	M. Zymbler, “Best-Match Time Series Subsequence Search on the Intel Many Integrated Core Architecture” in Proc. of the 19th East European Conf. on Advances in Databases and Information Systems, Poitiers, France, September 8–11, 2015. Lecture Notes in Computer Science (Springer, Heidelberg, 2015), Vol. 9282, pp. 275–286. doi: https://doi.org/10.1007/978-3-319-23135-8_1910.1007/978-3-319-23135-8_19.
	


