\documentclass[article,shortnames,nojss]{jss} \usepackage[english]{babel} \usepackage{array,thumbpdf} %% need no \usepackage{Sweave.sty} \SweaveOpts{engine=R, eps=FALSE} \DefineVerbatimEnvironment{Soutput}{Verbatim}{fontshape=n} %\VignetteIndexEntry{Computing and Visualizing Dynamic Time Warping Alignments in R: The dtw Package} %\VignetteDepends{dtw,proxy} %\VignetteKeywords{timeseries, alignment, dynamic programming, dynamic time warping} %\VignettePackage{dtw} \newcommand{\R}{\proglang{R}} \newcommand{\dtw}{\pkg{dtw}} \newtheorem{example}{Example} \author{Toni Giorgino\\National Research Council of Italy} \Plainauthor{Toni Giorgino} \title{Computing and Visualizing Dynamic Time Warping Alignments in \proglang{R}: The \pkg{dtw} Package} \Plaintitle{Computing and Visualizing Dynamic Time Warping Alignments in R: The dtw Package} \Shorttitle{\pkg{dtw}: Computing and Visualizing Dynamic Time Warping Alignments in \proglang{R}} \Volume{31} \Issue{7} \Month{August} \Year{2009} \Submitdate{2008-09-30} \Acceptdate{2009-06-26} \Abstract{ This introduction to the \proglang{R} package \pkg{dtw} is a (slightly) modified version of \cite{Giorgino:2009}, published in the \emph{Journal of Statistical Software}. Dynamic time warping is a popular technique for comparing time series, providing both a distance measure that is insensitive to local compression and stretches and the warping which optimally deforms one of the two input series onto the other. A variety of algorithms and constraints have been discussed in the literature. The \pkg{dtw} package provides an unification of them; it allows \R\ users to compute time series alignments mixing freely a variety of continuity constraints, restriction windows, endpoints, local distance definitions, and so on. The package also provides functions for visualizing alignments and constraints using several classic diagram types. } \Keywords{timeseries, alignment, dynamic programming, dynamic time warping} \Address{ Toni Giorgino\\ Institute of Biophysics (IBF-CNR)\\ National Research Council of Italy\\ Department of Biosciences, University of Milan\\ Via Celoria 26\\ I-20133, Milan, Italy\\ E-mail: \email{toni.giorgino@cnr.it} } \begin{document} <>= options(prompt= "R> ", continue = "+ ") @ \section{Introduction} Dynamic time warping (DTW) is the name of a class of algorithms for comparing series of values with each other. The rationale behind DTW is, given two time series, to stretch or compress them locally in order to make one resemble the other as much as possible. The distance between the two is computed, after stretching, by summing the distances of individual aligned elements (Figure~\ref{fig:example}). DTW algorithms have been proposed around 1970 in the context of speech recognition, to account for differences in speaking rates between speakers and utterances. The technique is useful, for example, when one is willing to find a low distance score between the sound signals corresponding to utterances \emph{now} and \emph{nooow} respectively, insensitive to the prolonged duration of the \emph{o} sound. Various types of DTW algorithms differ for the input feature space, the local distance assumed, presence of local and global constraints on the alignment, and so on. This freedom makes DTW a very flexible alignment approach. As said, it has been popularized in the '70s, when it was mainly applied to isolated word recognition \citep{Velichko1970,Sakoe1971}; since then, it has been employed for clustering and classification in countless domains: electro-cardiogram analysis \citep{Huang2002,dtw-ecg,Tuzcu2005}, clustering of gene expression profiles~\citep{Aach:2001:Bioinformatics:11395426,Hermans:2007:Bioinformatics:17237107}, biometrics~\citep{Faundez-Zanuy2007,Rath2003}, process monitoring~\citep{Gollmer1996}, just to name a few. The term \emph{time} series may even be misleading, because the warped dimensions can be other than time, e.g., an angle for shape recognition~\citep{Kartikeyan1989,WeiKX06,Tak2007}. This paper describes the \pkg{dtw} package for the \R\ statistical software \citep{R}, which provides a comprehensive solution for the computation and visualization of DTW alignments, whose theory is outlined in Section~\ref{sec:definition-algorithm}. The stable package version is available in source and binary form on the Comprehensive \proglang{R} Archive Network (CRAN), while development versions are hosted at \proglang{R}-Forge \citep{DtwCran}. The package allows users to select among a wide variety of algorithms and constraints described in the literature, simply passing arguments to a single function, \code{dtw}, introduced in Section~\ref{sec:computing-alignments}. The following sections discuss how the user can customize the classic constraints of the algorithm (local slope, endpoints, windowing). Finally, Section \ref{sec:plotting-results} describes how the alignment results can be conveniently plotted in a variety of ways. \section{Definition of the algorithm} \label{sec:definition-algorithm} In the following exposition, the choice of symbols will follow roughly the one in Chapter 4 of \citet{Rabiner1993}, which is an excellent reference on the subject. We assume that we want to compare two time series: a \emph{test}, or \emph{query}, $X=(x_1, \dots, x_N)$; and a \emph{reference} $Y=(y_1, \dots, y_M)$. For clarity, in the following we shall reserve the symbol $i=1 \dots N$ for indexing the elements in $X$ and $j=1 \dots M$ for those in $Y$. We also assume that a non-negative, \emph{local dissimilarity} function $f$ is defined between any pair of elements $x_i$ and $y_j$, with the shortcut: % \begin{equation} d(i,j) = f(x_i,y_j) \ge 0 \label{eq:cross-dist} \end{equation} Note that $d$, the cross-distance matrix between vectors $X$ and $Y$, is the only input to the DTW algorithm: elements $x_i$ and $y_j$ only enter the computation through the arguments of $f$. Therefore, the following discussion applies, with no loss of generality, to cases when $X$ and $Y$ are single- or multi-variate, continuous, nominal, or mixed, as long as $f(\cdot,\cdot)$ is suitably defined. While the most common choice is to assume the Euclidean distance, different definitions \citep[e.g., those provided by the \pkg{proxy} package,][]{proxy} may be useful as well. % Since $f$ needs not obey the triangle inequality, it is preferably % called ``dissimilarity'' rather than ``distance''. %, e.g., \cite{Gower1971} At the core of the technique lies the \emph{warping curve} $\phi(k), \, k = 1 \dots T $: % \begin{eqnarray*} \phi(k) & = & \Big( \phi_x(k),\phi_y(k) \Big) \quad \mbox{with}\\ \phi_x(k) & \in & \{ 1 \dots N \}, \\ \phi_y(k) & \in & \{ 1 \dots M \} \end{eqnarray*} The warping \emph{functions} $\phi_x$ and $\phi_y$ remap the time indices of $X$ and $Y$ respectively. Given $\phi$, we compute the average accumulated distortion between the warped time series $X$ and $Y$: % $$ d_\phi(X,Y) = \sum_{k=1}^T d(\phi_x(k),\phi_y(k)) \, m_\phi(k) / M_\phi $$ % where $m_\phi(k)$ is a per-step weighting coefficient and $M_\phi$ is the corresponding normalization constant, which ensures that the accumulated distortions are comparable along different paths. To ensure reasonable warps, \emph{constraints} are usually imposed on $\phi$. For example, \emph{monotonicity} is imposed to preserve their time ordering and avoid meaningless loops: % \begin{eqnarray*} \phi_x(k+1) \ge \phi_x(k) \\ \phi_y(k+1) \ge \phi_y(k) \\ \end{eqnarray*} The idea underlying DTW is to find the optimal alignment $\phi$ such that \begin{equation} D(X,Y) = \min_\phi d_\phi(X,Y) \label{e:dtw} \end{equation} In other words, one picks the deformation of the time axes of $X$ and $Y$ which brings the two time series as close as possible to each other. Remarkably, despite of the large search space, Equation~\ref{e:dtw} can be computed in $O(N \cdot M)$ time using dynamic programming. Interested readers can find the derivation of the solution in several references, e.g., \citet[Section~II.F]{Myers1980}. The output of the DTW algorithm is quite rich: first of all, the value of the function $D(X,Y)$ (minimum global dissimilarity, or ``DTW distance'') can be assumed as the stretch-insensitive measure of the ``inherent difference'' between two given time series. This distance has a straightforward application in hierarchical clustering and classification (e.g., with $k$-NN classifiers). The shape of the warping curve $\phi$ itself provides information about which point matches which, i.e., the pairwise correspondences of time points can be easily inspected. Once obtained, the warping function can be applied to any time series, allowing one to inspect post-hoc aligned signals or measure time distortions. \begin{figure}[t] \centering <>= library("dtw") data("aami3a") ref <- window(aami3a,start=0,end=2) test <- window(aami3a,start=2.7,end=5) plot(dtw(test,ref,k=TRUE),type="two",off=1,match.lty=2,match.indices=20) @ \caption{Aligning two time series: a simple example. A query (solid, left axis) and a reference (dashed, right axis) ECG time series, excerpted from \code{aami3a}.} \label{fig:example} \end{figure} \begin{figure} \centering \includegraphics[width=.95\textwidth]{figures/object-diagram} \caption{Block diagram of the main functions (rounded boxes) and object classes (rectangles) of the \dtw\ package and their relationship. Arrows going in functions are arguments, outgoing are returned values.} \label{fig:object-diagram} \end{figure} \section{Computing alignments} \label{sec:computing-alignments} After loading the \dtw\ package, alignments can be computed invoking the \code{dtw} function. In its simplest incarnation, the function takes two vector arguments representing the input time series, performs the minimization, and returns an object of class \code{dtw} encapsulating all the alignment information. The elements of the result can be retrieved with the usual \code{$} notation (see list in Table~\ref{tab:attributes}). By default, the \code{dtw} function computes a global alignment, with no windowing, a symmetric local continuity constraint, and the Euclidean local distance. In particular, the symmetric continuity constraint implies that arbitrary time compressions and expansions are allowed, and that all elements must be matched: % \begin{eqnarray} |\phi_x(k+1)-\phi_x(k)| & \le & 1, \nonumber \\ |\phi_y(k+1)-\phi_y(k)| & \le & 1. \label{eq:continuity} \end{eqnarray} % Several alternative forms for local continuity constraints however exist; Section~\ref{sec:step-patterns-local} discusses them and how they can be selected. Computing global alignments means that the time series' heads and tails are constrained to match each other. In other words, the following endpoint constraints are imposed: % \begin{eqnarray} & \phi_x(1)=\phi_y(1)=1; & \label{eq:startpoint} \\ & \phi_x(T)=N; \quad \phi_y(T)=M. & \label{eq:endpoint} \end{eqnarray} % As we shall see in Section~\ref{sec:unconstr-endp-pref}, one or both of these constraints can be relaxed in order to compute partial time series matches. % The function computes the alignment and returns the result in a list % (instance of class \code{dtw}). As the example below shows, using \code{dtw} is straightforward. \begin{example} The \code{aami3a} time series included in the package contains a reference electrocardiogram from the PhysioBank dataset \citep{Goldberger2000}. We extract two non-overlapping windows from it, and compute their optimal alignment with the default \code{dtw} settings. The two extracted segments are shown in Figure~\ref{fig:example}. <>= library("dtw") data("aami3a") ref <- window(aami3a,start=0,end=2) test <- window(aami3a,start=2.7,end=5) alignment <- dtw(test,ref) alignment$distance @ \end{example} In \code{alignment\$distance} one finds the cumulative (unnormalized) cost for the alignment, i.e., $M_\phi \cdot d_\phi(X,Y)$, while the warping functions $\phi_x$ and $\phi_y$ are found in components \code{\$index1} and \code{\$index2}. They are two integer vectors of the same length, listing the matching indices in the query and reference time series, respectively. Although the two vectors may be readily plotted with commands built-in in \R, in Section~\ref{sec:plotting-results} we shall show how this is accomplished even more easily with dedicated plotting commands. \begin{table} \centering \begin{tabular}{>{\ttfamily}ll>{\itshape}c} \hline \rmfamily {Element} & {Description} & \upshape {Remark} \\ \hline query & Query time series, if given & k,p \\ reference & Reference time series, if given & k,p \\ localCostMatrix & Local distance matrix & \\ stepPattern & Step pattern instance used & \\ N & Length of the query time series & \\ M & Length of the reference time series & \\ call & Function call & \\ \hline distance & Unnormalized minimum cumulative distance & \\ normalizedDistance & Normalized cumulative distance & n \\ index1 & Warping function $\phi_x(k)$ for the query & d \\ index2 & Warping function $\phi_y(k)$ for the reference & d \\ costMatrix & Computed cumulative cost matrix & k \\ directionMatrix & Transition chosen at each alignment point & k \\ jmin & End-point, for partial alignments & \\ \hline \end{tabular} \caption{Elements in a \code{dtw} result object, as of package version~1.13-1. Remarks: \emph{k} -- requires \code{keep = TRUE}; \emph{n} -- requires that the chosen step pattern is normalizable; \emph{d} -- unless \code{distance.only = TRUE}; \emph{p} -- unless a local distance matrix is used as input . } \label{tab:attributes} \end{table} \subsection{Step patterns and local slope constraints} \label{sec:step-patterns-local} The elegant formulation of the DTW alignment problem lends itself to a variety of extensions. For example, one usually wants to limit the number of consecutive elements which are ``skipped'' in either time series, i.e., are left unmatched (skipping elements is often entirely disallowed by a continuity constraint, indeed). Alignments are in general achieved by \emph{duplicating} elements, i.e., one lets a single time point in $X$ match multiple (consecutive) elements in $Y$, or vice-versa. How many repeated elements can be matched consecutively, or how many can be skipped, put limits on the local slope of the warping curve. This property can be controlled by a very flexible scheme called \emph{step patterns}. Step patterns list sets of allowed transitions between matched pairs, and the corresponding weights. In other words, step patterns specify the admissible values for $\phi(k+1)$ given $\phi(k), \phi(k-1)$, and so on. It may be useful to remark that in DTW there is no additive penalty for duplicating or skipping elements, as with other alignment algorithms like Smith-Waterman's or Levenshtein's. While most authors stay with the simplest recursion types, users of \dtw\ have access to almost the full range of step patterns defined in the literature, with no programming burden. Step patterns are represented by objects of class \code{stepPattern}. Patterns are selected by passing an appropriate instance to the \code{step.pattern} argument of the \code{dtw} function call, as in the following example: \begin{example} Compute the same alignment of the previous example, assuming the well-known \emph{asymmetric} pattern (bottom left of Figure~\ref{fig:step-patterns}). <>= alignment <- dtw(test,ref,step.pattern=asymmetric) alignment$distance @ \end{example} When printed, \code{stepPattern} objects display the corresponding DTW recursion in human-readable form; they can also be plotted via the \code{plot()} overloaded method, and transposed via \code{t()}. Transposing a pattern means that the role of the query and reference are interchanged. \begin{example} Display the recursion formula for the properly symmetric continuity constraint (top right in Figure~\ref{fig:step-patterns}), also known as the {\em symmetric} $P=0$ from \citet{Sakoe1978}. % <>= symmetric2 @ \end{example} The so-called \emph{symmetric} recursion printed above allows an unlimited number of elements of the query to be matched to a single element of the reference, and vice-versa; in other words, there is no limit in the amount of time expansion or compression allowed at any point. Once the package is loaded, an instance representing this recursion is pre-defined under the name \code{symmetric2}, which is also the default. For this recursion, the average cost per-step is computed by dividing the cumulative distance by $N+M$, where $N$ is the length of the query sequence and $M$ is the length of the reference. Other step patterns require different normalization formulas, as we shall see in Section~\ref{sec:normalization}. The normalized distance, if defined, is stored in the \code{$normalizedDistance} component of the result. Several step patterns have been discussed in the literature. A classic paper by \citet{Sakoe1978} classifies them according to two properties: their symmetry (symmetric/asymmetric), and the bounds imposed on the slope expressed through a parameter $P$. The eight step patterns shown in \citet[Table I]{Sakoe1978} are pre-defined in \dtw, with names \code{symmetricP1}, \code{asymmetricP05}, and so on.\footnote{Sakoe's ``asymmetric'' property should not be confused with the \code{asymmetric} step pattern.} All of them are normalizable. \citet[Chapter~4]{Rabiner1993} introduced a different classification with three attributes: local continuity constraint type (in Roman numerals, I to VII); slope weighting (Latin letters \emph{a} to \emph{d}); and their being ``smoothed'' or not (boolean). All of Rabiner's step patterns are available through the three-argument function \code{rabinerJuangStepPattern()}. Slope weighting types \emph{c} and \emph{d} are normalizable. Yet another classification (now obsoleted by the previous one) follows \citet{Myers1980}; it is similar to Rabiner's, except that only four types (I)-(IV) are defined. The corresponding instances are named like \code{typeIId}, \code{typeIIIc} and so on. The latter is a good approximation to the slope-limited pattern proposed by \citet{Itakura1975}. It should be noted that, in general, step patterns impose a lower and/or an upper bound to the \emph{local} slope of the alignment. In other words, they limit the maximum amount of time stretch and compression allowed at any point of the alignment. For example, the \code{asymmetric} step pattern limits time expansion to a factor of two; it would therefore be impossible to completely align a query with a reference more than twice as long. The Rabiner-Juang type IV, instead, generates the so-called Itakura parallelogram \citep{Itakura1975}. \begin{figure} \centering <>= plot(symmetric1) @ \includegraphics*[width=0.48\textwidth]{dtw-plot-symmetric1-step-pattern} % ---- <>= plot(symmetric2) @ \includegraphics[width=0.48\textwidth]{dtw-plot-symmetric2-step-pattern} \\ <>= plot(asymmetric) @ \includegraphics[width=0.48\textwidth]{dtw-plot-asymmetric-step-pattern} % ---- <>= plot(rabinerJuangStepPattern(4,"c",TRUE)) @ \includegraphics[width=0.48\textwidth]{dtw-plot-rj4c-step-pattern} \caption{Four well-known step patterns. Left to right, top to bottom: \code{symmetric1}, \code{symmetric2}, \code{asymmetric}, \code{rabinerJuangStepPattern(4, "c", TRUE)} (i.e., Rabiner-Juang's type IV with slope weighting \emph{c}). Numbers on transitions indicate the multiplicative weight $m_\phi$ for the local distance $d(i,j)$, should the corresponding step be followed. See \code{?stepPattern} for the full list.} \label{fig:step-patterns} \end{figure} \subsection{Normalization} \label{sec:normalization} Computing the average per-step distance along the warping curve is especially important in two cases: (1) when comparing alignments between time series of different lengths, to decide the best match (e.g., for classification); and (2) when performing partial matches. Each step patterns requires a different normalization function. In \dtw, objects of class \code{stepPattern} ``know'' the proper normalization formula they require, called ``normalization hint''. Not all step patterns are normalizable. For example, the \emph{quasi-symmetric} recursion (top left in Figure \ref{fig:step-patterns}, instance name \code{symmetric1}), like the symmetric one, does not restrict the slope; however, it favors ``diagonal'' steps over stair-stepping paths. Since similar warping curves have different weights, a path-independent per-step alignment cost is not defined. As a consequence, only the cumulative distance value is available in the \code{$distance} component of the result. Table~\ref{tab:normalization} lists the currently supported normalization rules, grouped by weighting type. ``R-J'' types follow Rabiner's aforementioned slope weight classification (cf.\ Figure~7 in~\cite{Myers1980}), while ``S-C'' refers to the symmetric/asymmetric categories introduced by Sakoe-Chiba (cf.\ Table I in~\cite{Sakoe1978}). Type (c$'$) includes asymmetric patterns similar to R-J type (c) with test and reference interchanged. Step patterns of this category have been used e.g., by \citet{Mori2006} (namely, instance \code{mori2006}), and by \citet{Oka1998}. Values $n$ and $m$ are the number of elements actually matched in the query and reference, respectively, i.e., % \begin{eqnarray*} n &=& \phi_x(T)-\phi_x(1)+1 \\ m &=& \phi_y(T)-\phi_y(1)+1. \end{eqnarray*} % For global alignments, they are equal to the lengths of the input time series, i.e., $n=N$ and $m=M$. \begin{table} \centering \begin{tabular}{ccccc} % \hline {Slope weighting} & {Symmetry} & {Formula} & {Hint} \\ {(Rabiner-Juang)} & {(Sakoe-Chiba)} & & \\ \hline R-J type (a) & --- & --- & \code{NA} \\ R-J type (b) & --- & --- & \code{NA} \\ R-J type (c) & S-C asymmetric & $n$ & \code{N} \\ R-J type (d) & S-C symmetric & $n+m$ & \code{N + M} \\ Type (c$'$) & S-C asymmetric & $m$ & \code{M} \\ \multicolumn{2}{c}{(Others)} & --- & \code{NA} % \hline \end{tabular} \caption{Normalization functions for well-known weighting types. } \label{tab:normalization} \end{table} \subsection{A note on indexing conventions and axes} \label{sec:note-index-conv} As a general rule, we stick to the convention that the first argument and indices refer to the query time series, and the second to the reference. This implies that when we print alignment-related matrices such as $d(i,j)$, the query index grows row-wise towards the bottom. The reader should not confuse this layout with plot axes, where the query and reference are usually arranged along the abscissa and ordinate respectively (Figure~\ref{fig:plotlike}). \begin{figure} \centering \includegraphics[width=.9\textwidth]{figures/plot-conventions} \caption{Different conventions for displaying distance matrices: plot-like (left) and matrix (right) arrangements. } \label{fig:plotlike} \end{figure} \subsection{Windowing and global constraints} A \emph{global} constraint, or \emph{window}, explicitly forbids warping curves to enter some region of the $(i,j)$ plane. A global constraint translates an a-priori knowledge about the fact that the time distortion is limited. For example, the well-known Sakoe-Chiba band \citep{Sakoe1978} enforces the additional constraint % $$| \phi_x(k)-\phi_y(k)| \le T_0$$ % where $T_0$ is the maximum allowable absolute time deviation between two matched elements. Intuitively, the constraint creates an allowed band of fixed width about the main diagonal of the alignment plane (Figure~\ref{fig:sakoechiba}). In \dtw, windowing constraints are enabled through the \code{window.type} argument of the \code{dtw} call; it can be either a character string (e.g., \code{"sakoechiba"}) or a function, specifying the shape of the allowed window. The window size, $T_0$, is passed through the \code{window.size} argument. It should be remarked that the Sakoe-Chiba band works well when $N \sim M$, but is inappropriate when the lengths of the two inputs differ significantly. In particular, when $|N-M|>T_0$, the $(N,M)$ endpoint lies outside of the band, thus violating (\ref{eq:endpoint}), and therefore no solution exists for a global alignment. In general, when enforcing global and/or local constraints one should take care that they are compatible with each other and with the time series' lengths. The \code{slantedBandWindow} creates a band centered around the jagged line segment which joins element $(1,1)$ to element $(N,M)$, and is $T_0$ elements wide along the first axis. In other words, the ``diagonal'' goes from one corner to the other of the possibly rectangular cost matrix, therefore having a slope of $M/N$, not 1. Arbitrary windows can be specified by supplying a function that takes (at least) the $i$ and $j$ integer arguments, and returns a boolean value specifying whether that match is allowed or not. Arguments unused in \code{dtw} are passed to the windowing function. More detail on windowing functions, including how to plot them (as in Figure~\ref{fig:sakoechiba}), are given in \code{?dtwWindowingFunctions}. % <>= dtwWindow.plot(sakoeChibaWindow, window.size=2,reference=17, query=13) @ \begin{figure} \centering \includegraphics{dtw-plot-sakoechiba-window} \caption{The allowed region for the warping curve under the \code{"sakoechiba"} global constraint. In this degenerate case, element $(N,M)$ at the upper-right is outside of the band, so the endpoint constraint (\ref{eq:endpoint}) can't be satisfied: this band is too narrow to be compatible with any global alignment. }\label{fig:sakoechiba} \end{figure} \begin{table} \centering \begin{tabular}{lcl} % \hline \hline {Attribute} & {Section in R-J} & {Documentation} \\ \hline Local continuity constraints & 4.7.2.3 & \code{?stepPattern} \\ Global path constraints & 4.7.2.4 & \code{?dtwWindowingFunctions} \\ Slope weighting & 4.7.2.5 & \code{?stepPattern} \\ % \hline \hline \end{tabular} \caption{Where to find documentation for various DTW parameters. R-J is the book by \citet{Rabiner1993}. } \label{tab:help} \end{table} \subsection{Unconstrained endpoints: Prefix and subsequence matches} \label{sec:unconstr-endp-pref} \nopagebreak Normally, the DTW distance is understood as a \emph{global} alignment, i.e., subject to the conditions (\ref{eq:startpoint}) and (\ref{eq:endpoint}). In certain applications, \emph{partial} matches are useful, and one consequently relaxes one or both the constraints. Interested readers can find a review of partial matching algorithms, applications, and the interplay with normalization, in a paper by \cite{Tormene2008}. The open-end (OE-DTW) or prefix-matching algorithm is defined by relaxing the end-point constraint (\ref{eq:endpoint}). The algorithm therefore returns the \emph{prefix} of the reference which best matches the query; it has been used e.g., by \citet{Sakoe1979} and \citet{Mori2006}. Open-end matching is achieved, in principle, by constructing several incomplete versions $Y^{(p)}$ of the reference $Y$, each truncated at the index $p=1, \dots, M$. One then computes their corresponding DTW distances from $X$, and picks the best match: % \begin{eqnarray} Y^{(p)} &=& (y_1, \dots, y_p) \nonumber \\ D_{OE}(X,Y) &=& \min_{1 \le j \le M} D(X,Y^{(j)}) \label{eq:dtwoe} \end{eqnarray} It is worthwhile noting that in Equation~\ref{eq:dtwoe} one compares alignments of different lengths with each other. This is only meaningful if the \emph{average} per-step distances are considered, and it is therefore essential to use normalizable step patterns. In \dtw, open-end alignment of time series is achieved simply setting the \code{open.end = TRUE} parameter. The element \code{$jmin} in the result will hold the size of the prefix matched, i.e., the $j$ that minimizes (\ref{eq:dtwoe}). This is the index of the last element matched in the reference, which in turn equals $\phi_y(T)$. Accordingly, when \code{oc} is a partial alignment result, \code{oc$jmin} equals \code{max(oc$index2)}. Subsequence match, also called ``unconstrained'' or open-begin-end (OBE-DTW), is achieved relaxing \emph{both} the start-point constraint in (\ref{eq:startpoint}) and the end-point constraint of (\ref{eq:endpoint}). Intuitively, subsequence matching discovers the contiguous part of the template which best matches the \emph{whole} query (Figure~\ref{fig:partial}). This form is used e.g., by \cite{Rabiner1978} as algorithm UE2-1; by \cite{Sakurai2007}; and others. Formally, we define $Y^{(p,q)}$ as the subsequence of $Y$ including elements $y_j$ with $p \le j \le q$; the OBE-DTW problem is to minimize the cumulative distance over the initial and final reference indices $p$ and $q$ simultaneously: % \begin{eqnarray} Y^{(p,q)} &=& (y_p, \dots, y_q) \nonumber \\ D_{OBE}(X,Y) &=& \min_{1 \le p \le q \le M} D(X,Y^{(p,q)}) \label{eq:dtwobe} \end{eqnarray} In \dtw, OBE-DTW alignments are computed setting both the \code{open.begin} and \code{open.end} arguments to \code{TRUE}.\footnote{The reader may wonder about the last remaining combination: constraining the ends of the time series but not their heads. This is also allowed, but not really necessary: a partial alignment with the query sequence reversed achieves the same effect.} Open-begin alignments are supported for step patterns with $n$-type normalization only. \begin{figure} \centering <>= idx<-seq(0,6.28,len=100) query<-sin(idx)+runif(100)/10 reference<-cos(idx) alignmentOBE <- dtw(query[44:88],reference, keep=TRUE,step=asymmetric, open.end=TRUE,open.begin=TRUE) plot(alignmentOBE,type="two",off=1) @ \caption{A partial alignment with both endpoints free: the whole query (solid line) is warped into a contiguous subsequence of the template (dashed). Partial alignments are computed setting the \code{open.begin = TRUE} and/or \code{open.end = TRUE} arguments. Code available in \code{example(dtw)}.} \label{fig:partial} \end{figure} \subsection{Dealing with multivariate time series} \label{sec:deal-with-mult} As seen in Section~\ref{sec:definition-algorithm}, the actual time series values only enter the DTW algorithm through their cross-distance matrix, i.e., matrix $d$ in Equation~\ref{eq:cross-dist}. In practice, the choice for the local distance function $f$ influences how ``strongly'' the alignment will avoid mismatching regions. A variety of dissimilarity functions are available, e.g., the Euclidean, squared Euclidean, Manhattan, Gower coefficient, and many others. Each of them takes two arguments, which may be multi-variate; for example, in speech recognition it is customary to align high-dimensional ``frames'', whose components are short-time spectral coefficients, or similar quantities. The \pkg{proxy} package maintains a database of distance definitions which can be used in cross-distance computations (see \code{summary(pr_DB)}). To deal with multivariate time series, one provides \code{dtw} with two matrices $X_{ic}$ and $Y_{jc}$, rather than two vectors as for the single-variate case. The time indices $i=1 \dots N$ and $j=1 \dots M$ are arranged along rows, while multivariate dimensions $c = 1 \dots C$ are arranged in columns. Multivariate time series objects (\R\ class \code{mts}) can be used, as well. By default, \code{dtw} assumes an Euclidean local distance, i.e., $$ d(i,j)^2 = \sum_{c=1}^C (X_{ic}-Y_{jc})^2 $$ Other distance functions can be selected through the \code{dist.method} character argument, which is forwarded to the \code{proxy::dist} function. Alternatively, the user can compute the local distance matrix by herself, and supply that instead of the time series. This is achieved invoking \code{dtw} with \emph{one} matrix argument rather than two. According to the convention in Section~\ref{sec:note-index-conv}, matrix rows (first index) are understood as the index in the query sequence, and matrix columns as the index in the reference. The following examples show both the two- and the one-argument call styles. \begin{example} Align two synthetic bivariate time series (a query with 10 time points, and a reference of length 5), assuming the Manhattan local distance for element pairs: <>= query <- cbind(1:10,1) ref <- cbind(11:15,2) dtw(query,ref,dist.method="Manhattan")$distance @ or equivalently, pre-computing \code{proxy::dist} <>= cxdist <- proxy::dist(query,ref,method="Manhattan") dtw(cxdist)$distance @ \end{example} The optimal alignment of {\em strings}, i.e., discrete values or {\em categorial} data, rather than of sequences of real values, is an ubiquitous task in bioinformatics, closely related to DTW. Functions to perform alignments of strings, like Levenshtein's distance or the Needleman-Wunsch algorithm, are implemented, e.g., in package \pkg{cba} (cf.\ function \code{sdists}), package \pkg{TraMineR} \citep{TraMineR}, and in \pkg{Bioconductor} \citep{Bioconductor}. The DTW algorithm can be used to align strings, as well, by defining a local distance function over all the pairs of symbols that can be formed in the alphabets considered. A cross-distance matrix can then be obtained from this function, and used in the single-argument call of \code{dtw}, as in the last example. \subsection{Computing several alignments at once} A common situation, arising e.g., in clustering and classification of time series, is to compute several DTW distances between pairs of elements in a database. For instance, we assume to have a database $\{ Q_k \}$ of $K$ time series, each of which is a single-variate vector $Q_k = ( q_{k1}, \dots, q_{kN} )$. We also assume that all $Q_k$ have the same length $N$. The database can therefore be encoded as a $K \times N$ matrix, $\{Q_k\}=q_{ki}$. To compute DTW alignments between all possible couples, we iterate two indices $h$ and $k$ over the rows of $q$, thus obtaining a $K \times K$ self-dissimilarity matrix: % $$\Lambda_{kh}=D(Q_k,Q_h)$$ % \noindent where $D$ is the dynamic time warping distance of Equation~\ref{e:dtw}. It should be noted that, since the DTW distance is not in general symmetric, $\Lambda$ will not be, either. As seen above, the \code{dtw} function leverages the \pkg{proxy} package for building the local distance matrix. There is, however, another important relationship between the two packages. Since dynamic time warping itself \emph{is} a dissimilarity function between vectors (understood as time series), \pkg{dtw} registers itself as a distance function in the database of distances \code{pr_DB}. This feature makes it straightforward to compute many-versus-many alignments. Once the package is loaded, if $q_{ki}$ is given as a matrix \code{q} (again, time series are arranged in rows), $\Lambda$ can be straightforwardly computed through \code{proxy::dist(q, q, method = "DTW")}.\footnote{We use the two-argument form because the single-argument form would return a symmetric distance matrix, while DTW is not -- in general -- symmetric.} Computing $K \times K$ alignments in this way is faster than iterating over the plain \code{dtw} call, because only the numeric value of the distance is computed, while the construction of the actual warping path is bypassed. Equally useful is the cross-distance form of \code{dist}, which can be used, for instance, to classify $K$ query time series at once (arranged in a matrix \code{q}) against a given database of $L$ templates (matrix \code{p}). Once the $K \times L$ cross-distance matrix is computed by \code{dist(q, p, method = "DTW")}, one could take column-wise minima to discover the template best matching each of the $K$ query elements. The dissimilarity matrices resulting from calls to \code{dist} can also be fed into the usual clustering functions. \subsection{Minimal variance matching} \citet{Latecki2007} proposed the minimal variance matching (MVM) algorithm to align a given test to a reference, allowing arbitrary portions of the template to be skipped (Figure~\ref{f:mvmsp}, left-hand side). Interestingly, MVM may be considered a special case of DTW with a large step pattern (Figure~\ref{f:mvmsp}, right-hand side). Like the asymmetric step pattern, the MVM algorithm constrains each query element to match exactly one time point in the template, so the degenerate empty-to-empty match is not allowed. Vice-versa, an arbitrary number of template elements can remain unmatched. Matched reference points must be in strict increasing order, although arbitrary gaps are allowed in the sequence. This translates into the continuity constraints % \begin{eqnarray*} \phi_x(k) &=& k, \qquad \mbox{with} \qquad k=1, \dots, N \\ \phi_y(k) &<& \phi_y(k+1) \end{eqnarray*} Due to the strictly-increasing requirement of the reference index, the minimum local slope of the warping curve is unity, and time compression is therefore ruled out. MVM matching is therefore undefined if the query is longer than the reference, and only one alignment exists if their lengths are equal. In \dtw, step pattern instances implementing MVM can be built through the function\linebreak \code{mvmStepPattern()}. The \code{elasticity} integer argument limits the maximum length of a contiguous subsequence which can be skipped in the reference time series. If no limit is desired, elasticity should be made at least as large as the reference length. \begin{figure} \centering % plot(dtw(reference,rep(query,2),keep=TRUE,step=mvmStepPattern(100),p=TRUE),type="three") \includegraphics[width=.7\textwidth]{figures/alignMVM} \includegraphics[width=.27\textwidth]{figures/mvmSP} \caption{Left: MVM alignment of a truncated query; ``jumps'' may occur in the middle of the query, e.g., around time 175. Note that since time compression is not allowed, the valley at query index 150 does not match index 125 in the reference. Right: a type (c) DTW step pattern equivalent to the MVM search algorithm.} \label{f:mvmsp} \end{figure} \subsection{A worked-out exercise} \label{sec:worked-out-exercise} As an additional example, we solve Exercise 4.7 in \citet[page 226]{Rabiner1993}. The first question is to find the best path through a $6 \times 6$ local distance matrix that is explicitly given in the text. We enter the given matrix as follows: <>= lm <- matrix(nrow = 6, ncol = 6, byrow = TRUE, c( 1, 1, 2, 2, 3, 3, 1, 1, 1, 2, 2, 2, 3, 1, 2, 2, 3, 3, 3, 1, 2, 1, 1, 2, 3, 2, 1, 2, 1, 2, 3, 3, 3, 2, 1, 2 )) @ Note that we entered \code{lm} according to matrix conventions: the first subscript indexes elements in the query time series, while the second indexes the reference (right hand side of Figure~\ref{fig:plotlike}). Therefore, the query is understood to grow row-wise towards the bottom, while the reference goes column-wise towards the right. Conversely, the exercise text displays the table in plot-like conventions ($i_x$ growing rightwards and $i_y$ upwards). The exercise requires a specific step pattern which is readily identified as the \code{asymmetric} recursion (bottom left of Figure~\ref{fig:step-patterns}). To solve the problem and find the best \emph{global} path going through the grid, we invoke \code{dtw} with one matrix argument. From the result, we extract the cost matrix (see next section) and the normalized distance, 7/6, thus solving the problem: <>= alignment <- dtw(lm,step=asymmetric,keep=TRUE) alignment$costMatrix alignment$normalizedDistance @ A follow-up question of the same exercise requires the optimal alignment between the whole test and any \emph{prefix} of the reference; in other words, we lift the constraint on the reference end-point. This is achieved with the setting \code{open.end = TRUE}. As above, we invoke \code{dtw} in its matrix form; the optimal end-point (that is, $\phi_y(T)$) is found in component \code{$jmin} of the result. <>= alignmentOE <- dtw(lm,step=asymmetric,keep=TRUE,open.end=TRUE) alignmentOE$jmin alignmentOE$normalizedDistance @ The exercise is thus solved. \subsection{Retrieving cost matrices} If the parameter \code{keep.internals} is \code{TRUE}, the local distance matrix and the cumulative cost matrix are preserved after the calculation, and stored in the result elements \code{$localCostMatrix} and \code{$costMatrix}, respectively. The matrices can be printed and manipulated as usual. Figure~\ref{fig:exercise-figure}, for example, shows how to display them along with the alignment path superimposed. In the left panel one can hand-check the result of the exercise solved in the previous section: no better warping path exists that passes by $(1,1)$ through $(6,6)$ under the constraints of the \code{asymmetric} pattern, and the cumulative cost is indeed found in the upper-right element of the right panel. An analogous plot style, suitable for larger alignments, is presented in Section~\ref{sec:plott-cost-dens}. % <>= lcm <- alignment$localCostMatrix image(x=1:nrow(lcm),y=1:ncol(lcm),lcm) text(row(lcm),col(lcm),label=lcm) lines(alignment$index1,alignment$index2) @ <>= ccm <- alignment$costMatrix image(x=1:nrow(ccm),y=1:ncol(ccm),ccm) text(row(ccm),col(ccm),label=ccm) lines(alignment$index1,alignment$index2) @ \begin{figure} \centering <>= <> @ <>= <> @ \includegraphics[width=0.48\textwidth]{dtw-costmatrix-figure-plot-left} \includegraphics[width=0.48\textwidth]{dtw-costmatrix-figure-plot-right} \caption{The local distance matrix for the exercise in Section~\ref{sec:worked-out-exercise} (left) and the corresponding cumulative distance matrix (right) under the \code{asymmetric} step pattern. The optimal alignment path is superimposed.} \label{fig:exercise-figure} \end{figure} \section{Displaying alignments} \label{sec:plotting-results} Plots are valuable tools to inspect time series pairs together with their alignments. Several plotting styles are available in \dtw\ for producing publication-quality figures. They are achieved through the \code{plot} method overloaded for objects of type \code{dtw}. The first two plot styles discussed below, namely two- and three-way, are used to inspect the time series themselves along with their alignment. For this reason, the two actual input time series must be available, and they have to be single-variate. For convenience, time series are retrieved from the alignment object itself, if available (elements \code{$query} and \code{$reference}); they can be overridden via the \code{xts} and \code{yts} plot arguments. A \emph{density} plot style is also available to display the relative costs of different warping curves. It builds on the global cost matrix only, and does not therefore require the knowledge of the two original time series. \subsection{Two-way plotting} One intuitive alignment visualization style places both time series in the same plane, and connects the matching point pairs with segments (see e.g., Figure~\ref{fig:example}). This plot style is selected passing the argument \code{type = "twoway"} to the plot call. The optional numeric argument \code{offset} can be used to visually separate query and reference time series; if set, the scales for the reference time series are displayed on the secondary (right-hand) vertical axis. Other options affecting the visual appearance are explained in the \code{?dtwPlotTwoWay} manual page. \begin{example} Generate the plot in Figure~\ref{fig:example}. <>= <> @ \end{example} \subsection{Three-way plotting} Another effective layout to display alignments places the query time series horizontally in a small lower panel, the reference time series vertically on the left; a larger inner panel holds the warping curve (Figure~\ref{f:mvmsp}). In this way, matching points can be recovered by tracing indices on the query time series, moving upwards until the warping curve is met, and then moving leftwards to discover the index of the reference matched. The advantage of this method is that the warping curve is directly exposed, so it becomes easy to visualize the impact of the local and global constraints. Three-way plots can be achieved with the \code{type = "threeway"} argument to the \code{plot} call. Match lines can also be explicitly visualized, as shown in Figure~\ref{fig:three-way}; documentation for all options is available in \code{?dtwPlotThreeWay}. \begin{figure} \centering <>= ## A noisy sine wave as query ## A cosine is for reference; sin and cos are offset by 25 samples idx<-seq(0,6.28,len=100) query<-sin(idx)+runif(100)/10 reference<-cos(idx) dtw(query,reference,step=asymmetric,keep=TRUE)->alignment ## Beware of the reference's y axis, may be confusing ## Equivalent to plot(alignment,type="three") #dtwPlotThreeWay(alignment) ## Highlight matches of chosen QUERY indices. We will do some index ## arithmetics to recover the corresponding indices along the warping ## curve hq <- (0:8)/8 hq <- round(hq*100) # indices in query for pi/4 .. 7/4 pi hw <- (alignment$index1 %in% hq) # where are they on the w. curve? hi <- (1:length(alignment$index1))[hw] # get the indices of TRUE elems dtwPlotThreeWay(alignment,match.indices=hi) @ \caption{Three-way plot of the \code{asymmetric} alignment between a noisy sine and a cosine in $[0,2\pi]$, with visual guide lines drawn every $\pi/4$. Code available in \code{example(dtwPlotThreeWay)}.} \label{fig:three-way} \end{figure} \subsection{Displaying the cost density} \label{sec:plott-cost-dens} A third plot style displays the ``cost density'' of alignments, and can therefore be useful to inspect qualitatively how much ``slack'' is present around the optimal alignment (Figure~\ref{fig:density}). By default, the \emph{cumulative} distance distribution is displayed as a density distribution with contours superimposed. The result is a terse display of the cumulative distance matrix, analogous to the right panel of Figure~\ref{fig:exercise-figure}. For normalizable step patterns, a plot of the \emph{per-step} density can also be requested with the \code{normalize = TRUE} argument. Since the density plot is based on the cumulative distance matrix, \code{keep = TRUE} is required in the \code{dtw} function call. Density plotting is selected with the \code{type = "density"} parameter to the plot call, and documented in \code{?dtwPlotDensity}. \begin{figure} \centering <>= dtw(query,reference,keep=TRUE,step=asymmetric)->alignment dtwPlotDensity(alignment,normalize=TRUE) @ \caption{Cost density plot: average per-step cost density of the sine-cosine global alignment. The local constraints and weighting are chosen from the \code{asymmetric} step pattern. Code in \code{example(dtwPlotDensity)}.} \label{fig:density} \end{figure} \section*{Computational details} The computing kernel of the \code{dtw} function is written in the \proglang{C} programming language for efficiency. Alignment computations are quite fast, as long as the cross-distance matrices fit in the machine's RAM. A standard quad-core Linux x86-64 PC with 4~GB of RAM and 4~GB of swap computes (using only one core) an unconstrained alignment of $100 \times 100$ time points in 7~ms, $6000 \times 6000$ points in under 10~s, and $8000 \times 8000$ points (close to the virtual memory limit) in 10~minutes. Larger problems may be addressed by approximate strategies, e.g., computing a preliminary alignment between downsampled time series \citep{Salvador2004}; indexing \citep{Keogh2005}; or breaking one of the sequences into chunks and then iterating subsequence matches. \section*{Acknowledgments} I am grateful to P.\ Tormene for testing the package extensively. Thanks to M.\ J.\ Harvey and the anonymous reviewers for their valuable suggestions. \bibliography{dtw} \end{document}