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

TL;DR

This paper presents a scalable parallel algorithm for subsequence similarity search in large time series using MPI and OpenMP on Intel Xeon Phi Knights Landing clusters, optimizing vector computations for improved performance.

Contribution

It introduces a novel parallel algorithm leveraging MPI and OpenMP for efficient subsequence search on KNL clusters, incorporating data structures and redundant computations for better vectorization.

Findings

Algorithm is highly scalable on real-world datasets

Effective utilization of vector computations on KNL processors

Parallelization across cluster nodes improves performance

Abstract

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 TimeWarping (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 it is highly scalable.

Tables

Table 1. Таблица 1: Аппаратная платформа экспериментов

Характеристика	Хост	Сопроцессор
Модель, Intel Xeon	X5680	Phi (KNC), SE10X
Количество физических ядер	2$\times$6	61
Гиперпоточность	2$\times$	4$\times$
Количество логических ядер	24	244
Частота, ГГц	3.33	1.1
Размер VPU, бит	128	512
Пиковая производительность, TFLOPS	0.371	1.076

Table 2. Таблица 2: Наборы данных для экспериментов на одном узле кластера

Набор данных	Вид	$\left|T\right|=m$	$\left|Q\right|=n$
Random Walk	Синтетический	${10}^6$	128
EPG	Реальный	$2.5\cdot {10}^5$	360

Table 3. Таблица 3: Наборы данных для экспериментов на кластерной системе

Набор данных	Вид	$\left|T\right|=m$	$\left|Q\right|=n$
Random Walk	Синтетический	$12.8\cdot {10}^7$	128, 512, 1024
ECG	Реальный	$12.8\cdot {10}^7$	432, 512, 1024

Equations

$$\begin{array}{c}DTW(X,Y)=d(m,m),\\ d(i,j)=(x_i-y_j{)}^2+\min \{\begin{array}{l}d(i-1,j)\\ d(i,j-1)\\ d(i-1,j-1)\end{array},\\ d(0,0)=0,d(i,0)=d(0,j)=\infty ,1⩽i⩽m,1⩽j⩽m.\end{array}$$

$$P=p_1,p_2,\dots ,p_t,\dots ,p_T,m⩽T⩽2n-1.$$

$$\exists i\forall kDTW(Q,T_{i,n})⩽DTW(Q,T_{k,n}),1⩽i,k⩽m-n+1.$$

$$ED(Q,C)=\sqrt{\sum _{i=1}^n(q_i-c_i{)}^2}.$$

$$\begin{array}{c}\widehat{t_i}=\frac{t_i-\mu }{\sigma },1⩽i⩽m;\mu =\frac{1}{m}\sum _{i=1}^m{t}_i;{\sigma }^2=\frac{1}{m}\sum _{i=1}^m{t}_i^2-{\mu }^2.\end{array}$$

$$bsf_{(i)}=\operatorname{min}(bsf_{(i-1)},\left\{\begin{array}[]{l l}+\infty&,LB(Q,T_{i,n})>bsf_{(i-1)}\\ DTW(Q,T_{i,n})&,otherwise\end{array})\right.$$

$$L{B}_{Kim}FL(Q,C)=ED({\hat{q}}_1,{\hat{c}}_1)+ED({\hat{q}}_n,{\hat{c}}_n).$$

$$L{B}_{Keogh}EC(Q,C)=\sum _{i=1}^n\{\begin{array}{ll}({\hat{c}}_i-u_i{)}^2,& if{\hat{c}}_i>u_i\\ ({\hat{c}}_i-{\ell }_i{)}^2,& if{\hat{c}}_i<{\ell }_i\\ 0,& otherwise\end{array}.$$

$$\begin{array}{c}{u}_i=\max ({\hat{q}}_{i-r},\dots ,{\hat{q}}_{i+r}),\\ {\ell }_i=\min ({\hat{q}}_{i-r},\dots ,{\hat{q}}_{i+r}).\end{array}$$

$$L{B}_{Keogh}EQ(Q,C):=L{B}_{Keogh}EC(C,Q).$$

$$\begin{array}{c}start=k\cdot \lfloor \frac{N}{F}\rfloor +1\\ len=\{\begin{array}{ll}\lfloor \frac{N}{F}\rfloor +(NmodF)+n-1,& k=F-1\\ \lfloor \frac{N}{F}\rfloor +n-1,& otherwise.\end{array}\end{array}$$

$${\tilde{T}}_{i,n}=\{\begin{array}{ll}{t}_i,t_{i+1},\dots ,t_{i+n-1},\underset{pad}{\underbrace{0,0,\dots ,0}},& ifnmodw>0\\ t_i,t_{i+1},\dots ,t_{i+n-1},& otherwise.\end{array}$$

$$S_T^n(i,j):={\tilde{t}}_{i+j-1}.$$

$$L_T^n(i,j):=L{B}_j(T_{i,n},Q).$$

$$B_T^n(i):=\underset{j=1}{\overset{l{b}_{max}}{\bigwedge }}(L_T^n(i,j)<bsf).$$

$$C_T^n(i,\cdot ):=S_T^n(k,\cdot ):B_T^n(i)=\text{TRUE}.$$

$$\begin{array}{c}po{s}_i:=k:p\cdot (i-1)+1⩽k⩽\lceil \frac{N}{i\cdot p}\rceil \text{leavevmode}\wedge \forall j,1⩽j⩽l{b}_{max},L{B}_T^n(k,j)<bsf.\end{array}$$

$$id{x}_i:=k:1⩽k⩽N\wedge \exists S_T^n(i,\cdot )\iff \exists T_{i,n}\iff k=(i-1)\cdot n+1.$$

Full text

1 Введение.

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

На сегодня мера DTW (Dynamic Time Warping, динамическая трансформация временной шкалы) [4] признается научным сообществом лучшей мерой схожести временных рядов для многих предметных областей [6], поскольку позволяет адекватно измерять схожесть временных рядов, которые имеют различные длины либо сдвиг, растяжение или сжатие друг относительно друга. Поскольку мера DTW имеет квадратичную вычислительную сложность относительно длины поискового запроса, были предложены методы снижения вычислительной нагрузки: индексирование [9], заблаговременный отказ от избыточных вычислений, повторное использование результатов ранее выполненных вычислений [21] и др. Последовательный алгоритм UCR-DTW [18], интегрирующий множество указанных техник, является на сегодня, вероятно, самым быстрым последовательным алгоритмом поиска похожих подпоследовательностей [22]. В то же время предложены параллельные алгоритмы поиска похожих подпоследовательностей для устройств FPGA [29, 30], многоядерных ускорителей архитектур GPU [22, 30] и Intel Many Integrated Core (MIC) [31].

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

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

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

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

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

2 Постановка задачи.

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

Дадим определения используемых терминов в соответствии с работой [18]. Временной ряд (time series) $T$ — это хронологически упорядоченная последовательность вещественных значений: $T=(t_{1},t_{2},\dots,t_{m})$, $t_{i}\in\mathbb{R}$. Число $m$ обозначается $\lvert T\rvert$ и называется длиной временного ряда.

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

Пусть даны два временных ряда $X=(x_{1},x_{2},\dots,x_{m})$ и $Y=(y_{1},y_{2},\dots,y_{m})$. Тогда динамической трансформацией временной шкалы (Dynamic Time Warping, DTW) называется мера схожести $DTW(X,Y)$, вычисляемая следующим образом [4]:

[TABLE]

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

В формуле (1) матрица $(d_{ij})\in\mathbb{R}^{m\times m}$ выражает соответствие между точками сравниваемых временных рядов и называется матрицей трансформации.

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

[TABLE]

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

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

Подпоследовательностью (subsequence) $T_{i,n}$ временного ряда $T$ называется непрерывное подмножество $T$, состоящее из $n$ элементов и начинающееся с позиции $i$: $T_{i,n}=(t_{i},t_{i+1},\dots,t_{i+n-1}),\;1\leqslant i\leqslant m-n+1,\;k\ll m$.

Определим решаемую в данной работе задачу поиска самой похожей подпоследовательности ряда в смысле меры DTW. Пусть имеется длинный временной ряд $T$ и заданный пользователем существенно более короткий временной ряд $Q$ (называемый поисковым запросом (query)), где $m=\lvert T\rvert\gg\lvert Q\rvert=n$. Тогда самой похожей подпоследовательностью (best match) назовем такую подпоследовательность $T_{i,n}$, форма которой максимально похожа на запрос $Q$ в смысле меры DTW:

[TABLE]

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

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

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

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

[TABLE]

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

Z-нормализация. Перед вычислением меры схожести запрос и подпоследовательность временного ряда необходимо подвергнуть z-нормализации [28]. Z-нормализацией временного ряда $T$ называется временной ряд $\hat{T}=(\hat{t}_{1},\dots,\hat{t}_{m})$, элементы которого вычисляются следующим образом:

[TABLE]

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

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

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

[TABLE]

Алгоритм UCR-DTW использует следующие нижние границы схожести: $LB_{Kim}FL$ [18], $LB_{Keogh}EC$ и $LB_{Keogh}EQ$ [9], которые применяются каскадным образом.

Нижняя граница $LB_{Kim}FL$ представляет собой Евклидово расстояние между первой и последней парами точек $Q$ и $C$:

[TABLE]

Нижняя граница $LB_{Keogh}EC$ показывает минимальную схожесть между оболочкой запроса $E$ (envelope) и кандидатом $\hat{C}$ и вычисляется следующим образом:

[TABLE]

В формуле (8) последовательности $U=(u_{1},..,u_{n})$ и $L=(\ell_{1},..,\ell_{n})$ обозначают верхнюю (upper) и нижнюю (lower) границы оболочки запроса $Q$, которые вычисляются по формуле (9):

[TABLE]

Нижняя граница $LB_{Keogh}EQ$ представляет собой Евклидово расстояние между запросом $Q$ и оболочкой кандидата $C$, то есть по сравнению с $LB_{Keogh}EC$ роли запроса и кандидата меняются местами:

[TABLE]

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

3 Параллельный алгоритм PhiBestMatch.

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

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

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

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

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

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

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

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

Пусть $N=\lvert T\rvert-n+1=m-n+1$ — количество подпоследовательностей, которые необходимо обработать, $F$ — количество фрагментов, $T^{(k)}$ — $k$-й фрагмент временного ряда $T=(t_{1},t_{2},\dots,t_{m})$, где $0\leqslant k\leqslant F-1$. Тогда фрагмент $T^{(k)}$ определяется как подпоследовательность $T_{start,\,len}$, где

[TABLE]

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

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

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

Выравнивание данных выполняется следующим образом. Пусть обработка подпоследовательности $T_{i,n}$ ряда $T$ осуществляется с использованием векторного регистра, вмещающего $w$ вещественных чисел. Если длина подпоследовательности не кратна $w$, то подпоследовательность дополняется фиктивными нулевыми элементами. Обозначим количество фиктивных элементов за $pad=w-(n\;mod\;w)$, тогда выровненная подпоследовательность $\tilde{T}_{i,n}$ определяется следующим образом:

[TABLE]

Согласно определению (1), $\forall Q,\,C:DTW(Q,C)=DTW(\tilde{Q},\tilde{C})$. В дальнейшем изложении мы предполагаем, что запрос и подпоследовательности-кандидаты выровнены в памяти, и для упрощения используем обозначения и $Q$ и $C$ соответственно. Выравнивание подпоследовательностей позволяет избежать накладных расходов на разбивание циклов на несколько итераций (loop peeling).

Все (выровненные) подпоследовательности временного ряда сохраняются в виде матрицы подпоследовательностей для обеспечения векторизации вычислений. Матрица подпоследовательностей $S_{T}^{n}\in\mathbb{R}^{N\times(n+pad)}$ определяется следующим образом:

[TABLE]

Обозначим количество нижних границ схожести, используемых в алгоритме поиска самой похожей подпоследовательности, за $lb_{max}$ ($lb_{max}\geqslant 1$), и за $LB_{1}$, $LB_{2}$, $\dots$, $LB_{lb_{max}}$ обозначим эти границы, перечисленные в порядке их вычисления в каскаде применения нижних оценок. Тогда $L_{T}^{n}\in\mathbb{R}^{N\times lb_{max}}$, матрица нижних границ схожести всех подпоследовательностей длины $n$ временного ряда $T$ с поисковым запросом $Q$, имеет следующий вид:

[TABLE]

Карта схожести представляет собой вектор-столбец $B_{T}^{n}\in\mathbb{B}^{N}$, который для каждой подпоследовательности длины $n$ ряда $T$ хранит конъюнкцию применения всех нижних границ схожести этой подпоследовательности с текущим значением порога (bsf):

[TABLE]

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

Параллельная обработка матрицы кандидатов использует разбиение строк этой матрицы на сегменты, обрабатываемые отдельными нитями одного узла кластера. Обозначим число нитей, используемых параллельным алгоритмом, за $p\;(p\geqslant 1)$. Введем размер сегмента $s\in\mathbb{N}$ ($s\leqslant\lceil\frac{N}{p}\rceil$) — количество строк матрицы кандидатов, обрабатываемых одной нитью. Тогда матрица кандидатов $C_{T}^{n}\in\mathbb{R}^{(s\cdot p)\times(n+pad)}$ определяется следующим образом:

[TABLE]

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

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

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

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

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

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

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

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

[TABLE]

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

Чтобы вывести индекс самой похожей на запрос подпоследовательности, введен индексный массив $Idx\in\mathbb{N}^{s\cdot p}$, который предназначен для хранения позиции подпоследовательности во временном ряде и определяется следующим образом:

[TABLE]

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

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

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

4 Вычислительные эксперименты.

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

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

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

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

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

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

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

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

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

В экспериментах исследовались производительность и масштабируемость алгоритма PhiBestMatch. Под производительностью понимается время работы алгоритма без учета времени загрузки данных в память и выдачи результата. Масштабируемость параллельного алгоритма означает его способность адекватно адаптироваться к увеличению параллельно работающих вычислительных элементов (процессов, процессоров, нитей и др.) и характеризуется ускорением и параллельной эффективностью, которые определяются следующим образом [1]. Ускорение и параллельная эффективность параллельного алгоритма, запускаемого на $k$ нитях, вычисляются как $s(k)=\tfrac{t_{1}}{t_{k}}$ и $e(k)=\tfrac{s(k)}{k}$ соответственно, где $t_{1}$ и $t_{k}$ — время работы алгоритма на одной и $k$ нитях соответственно.

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

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

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

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

Поисковые запросы с длиной $n\geqslant 512$, равно как и значение параметра $r=1$ требуются на практике в ряде приложений, требующих при определении схожести подпоследовательностей как можно более высокую точность, например, в медицине при исследовании ЭКГ [15], в энтомологии [2], в астрономии [19] и др.

Результаты экспериментов по исследованию масштабируемости алгоритма представлены на \figreffig:PhiBestMatch-SingleNode-Scalability-RandomWalk и 4. Можно видеть, что 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\leqslant r\leqslant n$ и $n\geqslant 512$ соответственно.

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

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

В экспериментах исследовалось ускорение масштабируемости (scaled speedup) параллельного алгоритма, которое определяется как ускорение, демонстрируемое алгоритмом при линейном увеличении объема данных и количества используемых вычислительных узлов [12] и вычисляется следующим образом: $s_{scaled}=\tfrac{p\cdot m}{t_{p(p\cdot m)}}$, где $p$ — количество задействованных вычислительных узлов, $m$ — объем исходных данных, $t_{p(p\cdot m)}$ — время выполнения алгоритма на $p$ узлах при обработке исходных данных, имеющих объем $p\cdot m$.

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

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

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

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

5 Заключение.

В работе рассмотрена проблема поиска похожих подпоследовательностей в сверхбольших временных рядах (сотни миллиардов точек) на основе использования меры схожести 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).

{biblio}

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

• [1]

Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. СПб.: БХВ-Петербург, 2002.

• [2]

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: 10.1145/1376616.1376656.

• [3]

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: 10.1145/197405.197406.

• [4]

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.

• [5]

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: 10.1109/HOTCHIPS.2012.7476487.

• [6]

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.

• [7]

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: 10.1007/s00778-006-0040-z.

• [8]

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: 10.1007/978-3-642-33518-1_1.

• [9]

Keogh E.J., Ratanamahatana C.A. Exact indexing of dynamic time warping // Knowl. Inf. Syst. 2005. 7, N 3. 358–386. DOI: 10.1007/s10115-004-0154-9.

• [10]

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.

• [11]

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.

• [12]

Kumar V., Grama A., Gupta A., Karypis G. Introduction to parallel computing. Benjamin/Cummings, 1994.

• [13]

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: 10.1007/11733836_7.

• [14]

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: 10.1145/1188455.1188673.

• [15]

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: 10.1016/j.cmpb.2018.04.009.

• [16]

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.

• [17]

Pearson K. The problem of the random walk // Nature. 1905. 72, N 1865. 294. DOI: 10.1038/072342a0.

• [18]

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: 10.1145/2339530.2339576.

• [19] Rebbapragada U., Protopapas P., Brodley C.E., Alcock C. Finding anomalous periodic time series // Machine Learning. 2009. 74, N 3. 281–313. DOI: 10.1007/s10994-008-5093-3.
• [20]

Sakoe H., Chiba S. Readings in speech recognition. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc., 1990. 159–165.

• [21]

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: 10.1109/ICDE.2007.368963.

• [22]

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:10.1109/ICDM.2010.21.

• [23]

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.

• [24]

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: 10.1007/978-3-319-55669-7_17.

• [25]

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: 10.1109/ICACCI.2015.7275633.

• [26]

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: 10.1109/ICCCT.2011.6075111.

• [27]

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: 10.1109/CISIS.2009.77.

• [28]

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: 10.1109/CODES-ISSS.2013.6659005.

• [29]

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: 10.1145/2435264.2435277.

• [30]

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: 10.1109/ICASSP.2012.6289085.

• [31]

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: 10.1007/978-3-319-23135-8_19.

Поступила в редакцию

04.12.2018

\classify 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. Kraeva1 and M.L. Zymbler2 \addressen1 South Ural State University, Department of System Programming; prospekt Lenina 76, Chelyabinsk, 454080, Russia; Master student, e-mail: kraevaya@susu.ru

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

{abstracten} 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. \keywordsentime series, similarity search, Dynamic Time Warping, parallel algorithm, OpenMP, Intel Xeon Phi, Knights Landing, data layout, vectorization.

{biblio_lat}

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

• [1]

V. V. Voevodin, and Vl. V. Voevodin, The Parallel Computing. (BHV-Petersburg, St. Petersburg, 2002).

• [2]

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: 10.1145/1376616.1376656.

• [3]

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: 10.1145/197405.197406.

• [4]

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.

• [5]

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: 10.1109/HOTCHIPS.2012.7476487.

• [6]

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).

• [7]

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: 10.1007/s00778-006-0040-z.

• [8]

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: 10.1007/978-3-642-33518-1_1.

• [9]

E. J. Keogh, and C. A. Ratanamahatana, ‘‘Exact Indexing of Dynamic Time Warping,’’ Knowl. Inf. Syst. 7 (3). pp. 358–386. (2005). doi: 10.1007/s10115-004-0154-9.

• [10]

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.

• [11]

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.

• [12]

V. Kumar, A. Grama, A. Gupta, and G. Karypis. Introduction to Parallel Computing (Benjamin/Cummings, 1994).

• [13]

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: 10.1007/11733836_7.

• [14]

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: 10.1145/1188455.1188673.

• [15]

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: 10.1016/j.cmpb.2018.04.009.

• [16]

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.

• [17]

K. Pearson, ‘‘The Problem of the Random Walk,’’ Nature. 72 (1865). 294. (1905). doi: 10.1038/072342a0.

• [18]

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: 10.1145/2339530.2339576.

• [19] U. Rebbapragada, P. Protopapas, C. E. Brodley, and C. Alcock, ‘‘Finding Anomalous Periodic Time Series,’’ Machine Learning 74 (3), 281–313 (2009). doi: 10.1007/s10994-008-5093-3.
• [20]

H. Sakoe, and S. Chiba, Readings in Speech Recognition. (San Francisco, CA, USA, Morgan Kaufmann Publishers Inc., 1990), pp. 159–165.

• [21]

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: 10.1109/ICDE.2007.368963.

• [22]

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: 10.1109/ICDM.2010.21.

• [23]

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.

• [24]

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: 10.1007/978-3-319-55669-7_17.

• [25]

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: 10.1109/ICACCI.2015.7275633.

• [26]

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: 10.1109/ICCCT.2011.6075111.

• [27]

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: 10.1109/CISIS.2009.77.

• [28]

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: 10.1109/CODES-ISSS.2013.6659005.

• [29]

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: 10.1145/2435264.2435277.

• [30]

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: 10.1109/ICASSP.2012.6289085.

• [31]

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: 10.1007/978-3-319-23135-8_19.