SparseDTW: A Novel Approach to Speed up Dynamic Time Warping

1.1 MB·English

by  Ghazi Al-Naymat

SparseDTW: A Novel Approach to Speed up Dynamic Time Warping
GhaziAl-Naymat1,∗ SanjayChawla2 JavidTaheri2
1 SchoolofComputerScienceandEngineering
TheUniversityofNewSouthWales
Sydney,NSW2052,Australia
Email:[email protected]
2
1
0 2 SchoolofInformationTechnologies
2 TheUniversityofSydney,Australia

n Email:{chawla,javidt}@it.usyd.edu.au
a
J

3
Abstract
1

] We present a new space-efficient approach, (SparseDTW), to compute the Dynamic Time Warping (DTW) distance be-
B
tween two time series that always yields the optimal result. This is in contrast to other known approaches which typically
D
sacrifice optimality to attain space efficiency. The main idea behind our approach is to dynamically exploit the existence
. ofsimilarityand/orcorrelationbetween thetimeseries. The more thesimilaritybetweenthetimeseriesthe lessspacere-
s
c quired to compute the DTW between them. To the best of our knowledge, all other techniques to speedup DTW, impose
[
aprioriconstraintsanddonotexploitsimilaritycharacteristicsthatmaybepresentinthedata. Weconductexperimentsand

1 demonstratethatSparseDTWoutperformspreviousapproaches.
v
9
6 Keywords: Timeseries,Similaritymeasures,Dynamictimewarping,Datamining
9
2
1 Introduction
.
1
0
2 Dynamic time warping (DTW) uses the dynamic programming paradigm to compute the alignment between two time
1 series. An alignment “warps” one time series onto another and can be used as a basis to determine the similarity between
v: thetimeseries. DTW hassimilaritiestosequencealignmentinbioinformaticsandcomputationallinguisticsexceptthatthe
i matchingprocessinsequencealignmentandwarpinghavetosatisfyadifferentsetofconstraintsandthereisnogapcondition
X
inwarping. DTW firstbecamepopularinthespeechrecognitioncommunity(Sakoe&Chiba1978)whereithasbeenused
r todetermineifthetwospeechwave-formsrepresentthesameunderlyingspokenphrase. Sincethenithasbeenadoptedin
a
manyotherdiverseareasandhasbecomethesimilaritymetricofchoiceintimeseriesanalysis(Keogh&Pazzani2000).
Likeinsequencealignment,thestandardDTW algorithmhasO(mn)spacecomplexitywheremandnarethelengthsof
thetwosequencesbeingaligned. Thislimitsthepracticalityofthealgorithmintodays“datarichenvironment”wherelong
sequencesareoftenthenormratherthantheexception. Forexample,considertwotimeserieswhichrepresentstockprices
at one second granularity. A typical stock is traded for at least eight hours on the stock exchange and that corresponds to
alengthof8×60×60 = 28800. Tocomputethesimilarity, DTW wouldhavetostoreamatrixwithatleast800million
entries!
Figure1(a)showsanexampleofanalignment(warping)betweentwosequencesSandQ. Itisclearthatthereareseveral
possiblealignmentsbutthechallengeistoselecttheonewhichhastheminimaloveralldistance.Thealignmenthastosatisfy
severalconstraintswhichwewillelaborateoninSection3.
Salvador&Chan(2007)haveprovidedasuccinctcategorizationofdifferenttechniquesthathavebeenusedtospeedup
DTW:
∗TheworkwasdonewhiletheauthoratSchoolofInformationTechnologies,TheUniversityofSydney,Australia.
1
• Constraints: By adding additional constraints the search space of possible alignments can be reduced. Two well
knownexemplarsofthisapproacharetheSakoe&Chiba(1978)andtheItakura(1975)constraintswhichlimithow
farthealignmentcandeviatefromthediagonal. Whiletheseapproachesprovideareliefinthespacecomplexity,they
donotguaranteetheoptimalityofthealignment.
• DataAbstraction: Inthisapproach,thewarpingpathiscomputedatalowerresolutionofthedataandthenmapped
backtotheoriginalresolution(Salvador&Chan2007). Again,optimalityofthealignmentisnotguaranteed.
• Indexing:Keogh&Ratanamahatana(2004),Sakuraietal.(2005),andLemire(2009)proposedanindexingapproach,
whichdoesnotdirectlyspeedupDTWbutlimitsthenumberofDTWcomputations. Forexample,supposethereexists
a database D of time series sequences and a query sequence q. We want to retrieve all sequences d ∈ D such that
DTW(q,d)<(cid:15). TheninsteadofcheckingqagainsteachandeverysequenceinD,aneasytocalculatelowerbound
functionLBFisfirstappliedbetweenqandD. Theargumentworksasfollows:
1. Byconstruction,LBF(q,d)<DTW(q,d).
2. Therefore,ifLBF(q,d)>(cid:15)thenDTW(q,d)>(cid:15)andDTW(q,d)doesnothavetobecomputed.
1.1 Main Contribution
Themaininsightbehindourproposedapproach,SparseDTW,istodynamicallyexploitthepossibleexistenceofinherent
similarity and correlation between the two time series whose DTW is being computed. This is the motivation behind the
Sakoe-ChibabandandtheItakuraParellelogrambutourapproachhasthreedistinctadvantages:
1. Bands in SparseDTW evolve dynamically and are, on average, much smaller than the traditional approaches. We
alwaysrepresentthewarpingmatrixusingsparsematrices,whichleadstobetteraveragespacecomplexitycompared
tootherapproaches(Figure9).
2. SparseDTW always yields the optimal warping path since we never have to set apriori constraints independently of
thedata. Forexample, inthetraditionalbandedapproaches, asub-optimalpathwillresultifallthepossibleoptimal
warpingpathshavetocrossthebands.
3. SinceSparseDTW yieldsanoptimalalignment,itcaneasilybeusedinconjunctionwithlowerboundapproaches.
1.2 Paper Outline
Therestofthepaperisorganizedasfollows: Section2describesrelatedworkonDTW.TheDTW algorithmisdescribed
inSection3. InSection4, wegiveanoverviewofthetechniquesusedtospeedupDTWbyaddingconstraints. Section5
reviews the Divide and Conquer approach for DTW which is guaranteed to take up O(m + n) space and O(mn) time.
Furthermore,weprovideanexamplewhichclearlyshowsthatthedivideandconquerapproachfailstoarriveattheoptimal
DTW result. The SparseDTW algorithm is introduced with a detailed example in Section 6. We analyze and discuss our
resultsinSection7,followedbyourconclusionsinSection8.
2 RelatedWork
DTWwasfirstintroducedinthedataminingcommunityinthecontextofminingtimeseries(Berndt&Clifford1994).
Since it is a flexible measure for time series similarity it is used extensively for ECGs (Electrocardiograms) (Caiani et al.
1998),speechprocessing(Rabiner&Juang1993),androbotics(Schmilletal.1999). ItisimportanttoknowthatDTW isa
measurenotametric,becauseDTW doesnotsatisfythetriangularinequality.
SeveraltechniqueshavebeenintroducedtospeedupDTW and/orreducethespaceoverhead(Hirschberg1975,Yietal.
1998,Kimetal.2001,Keogh&Ratanamahatana2004,Lemire2009).
Divideandconquer(DC)heuristicproposedbyHirschberg(1975); thatisadynamicprogrammingalgorithmthatfinds
the least cost sequence alignment between two strings in linear space and quadratic time. The algorithm was first used in
speechrecognitionareatosolvetheLongestCommonSubsequence(LCSS).Howeveraswewillshowwiththehelpofan
example,DCdoesnotguaranteetheoptimalityoftheDTW distance.
6 Q
D 1 2 2 1 0 1 1 2 1 2
5 3 4 5 6 10 19 23 27 28 32 33
4 13 88 9 15 26 28 32 31 37 36
4 5 29 17 17 25 40 42 44 40 47 45
3 33 18 18 21 30 34 38 39 43 44
S S
3 3 37 19 19 22 30 34 38 39 43 44
2 38 19 19 20 24 25 26 26 27 27
2
3 42 20 20 23 29 28 26 27 30 28
4 51 24 24 29 39 37 37 31 36 32
1
Q 2 52 24 24 25 29 30 31 31 32 32
3 56 25 25 28 34 33 34 32 35 33
0
(a) ThealignmentofmeasurementsformeasuringtheDTWdistancebetween(b) ThewarpingmatrixDproducedbyDTW;highlightedcellsconstitutethe
thetwosequencesSandQ. optimalwarpingpath.
Figure1.IllustrationofDTW.
Sakoe & Chiba (1978) speed up the DTW by constraining the warping path to lie within a band around the diagonal.
However,iftheoptimalpathcrossestheband,theresultwillnotbeoptimal.
Keogh&Ratanamahatana(2004)andLemire(2009)introducedefficientlowerboundsthatreducethenumberofDTW
computationsinatimeseriesdatabasecontext. However,theselowerboundsdonotreducethespacecomplexityoftheDTW
computation,whichistheobjectiveofourwork.
Sakurai et al. (2005) presented FTW, a search method for DTW; it adds no global constraints on DTW. Their method
designed based on a lower bounding distance measure that approximates the DTW distance. Therefore, it minimizes the
numberofDTWcomputationsbutdoesnotincreasethespeedtheDTWitself.
Salvador & Chan (2007) introduced an approximation algorithm for DTW called FastDTW. Their algorithm begins by
usingDTWinverylowresolution,andprogressestoahigherresolutionlinearlyinspaceandtime. FastDTWisperformedin
threesteps: coarseningshrinksthetimeseriesintoasmallertimeseries;thetimeseriesisprojectedbyfindingtheminimum
distance (warping path) in the lower resolution; and the warping path is an initial step for higher resolutions. The authors
refinedthewarpingpathusinglocaladjustment. FastDTW isanapproximationalgorithm,andthusthereisnoguaranteeit
willalwaysfindtheoptimalpath.Itrequiresthecoarseningsteptoberunseveraltimestoproducemanydifferentresolutions
ofthetimeseries. TheFastDTW approachdependsonaradiusparameterasaconstraintontheoptimalpath;however,our
techniquedoesnotplaceanyconstrainwhilecalculatingtheDTW distance.
DTW has been used in data streaming problems. Capitani & Ciaccia (2007) proposed a new technique, Stream-DTW
(STDW).ThismeasureisalowerboundoftheDTW.Theirmethodusesaslidingwindowofsize512. Theyincorporateda
bandconstraint,forcingthepathtostaywithinthebandfrontiers,asin(Sakoe&Chiba1978).
AlltheabovealgorithmswereproposedeithertospeedupDTW,byreducingitsspaceandtimecomplexity,orreducing
thenumberofDTW computations. Interestingly, theapproachofexploitingthesimilaritybetweenpoints(correlation)has
never,tothebestofourknowledge,beenusedinfindingtheoptimalitybetweentwotimeseries. SparseDTW considersthe
correlationbetweendatapoints,thatallowsustouseasparsematrixtostorethewarpingmatrixinsteadofafullmatrix. We
donotbelievethattheideaofsparsematrixhasbeenconsideredpreviouslytoreducetherequiredspace.
3 DynamicTimeWarping(DTW)
DTWisadynamicprogrammingtechniqueusedformeasuringthesimilaritybetweenanytwotimeserieswitharbitrary
lengths. ThissectiongivesanoverviewofDTW andhowitiscalculated. Thefollowingtwotimeseries(Equations1and2)
Algorithm1DTW:ThestandardDTWalgorithm.
Input: S: Sequenceoflengthn,Q: Sequenceoflengthm.
Output: DTWdistance.
1: InitializeD(i,1)⇐iδforeachi
2: InitializeD(1,j)⇐jδforeachj
3: forallisuchthat2≤i≤ndo
4: forallj suchthat2≤j ≤mdo
5: UseEquation3tocomputeD(i,j)
6: endfor
7: endfor
8: return D(n,m)
willbeusedinourexplanations.
S = s ,s ,s ,···,s ,···,s (1)
1 2 3 i n
Q = q ,q ,q ,···,q ,···,q (2)
1 2 3 j m
WherenandmrepresentthelengthoftimeseriesS andQ,respectively. iandj arethepointindicesinthetimeseries.
DTWisatimeseriesassociationalgorithmthatwasoriginallyusedinspeechrecognition(Sakoe&Chiba1978).Itrelates
twotimeseriesoffeaturevectorsbywarpingthetimeaxisofoneseriesontoanother.
As a dynamic programming technique, it divides the problem into several sub-problems, each of which contribute in
calculatingthedistancecumulatively. Equation3showstherecursionthatgovernsthecomputationsis:

D(i−1,j)

D(i,j)=d(i,j)+min D(i−1,j−1) (3)
 D(i,j−1).
ThefirststageintheDTW algorithmistofillalocaldistancematrixd. Thatmatrixhasn×melementswhichrepresent
the Euclidean distance between every two points in the time series (i.e., distance matrix). In the second stage, it fills the
warpingmatrixD (Figure1(b))onthebasisofEquation3. Lines1to7inAlgorithm1illustratetheprocessoffillingthe
warpingmatrix. Werefertothecostbetweentheithandthejthelementsasδasmentionedinline1and2.
Afterfillingthewarpingmatrix,thefinalstagefortheDTW istoreporttheoptimalwarpingpathandtheDTW distance.
Warping path is a set of adjacent matrix elements that identify the mapping between S and Q. It represents the path that
minimizestheoveralldistancebetweenS andQ. ThetotalnumberofelementsinthewarpingpathisK,whereK denotes
thenormalizingfactorandithasthefollowingattributes:
W =w ,w ,...,w
1 2 K
max(|S|,|Q|)≤K <(|S|+|Q|)
Everywarpingpathmustsatisfythefollowingconstraints(Keogh&Ratanamahatana2004,Salvador&Chan2007,Sakoe
&Chiba1978):
1. Monotonicity: AnytwoadjacentelementsofthewarpingpathW,w = (w ,w )andw = (w(cid:48),w(cid:48)),followthe
k i j k−1 i j
inequalities,w −w(cid:48) ≥0andw −w(cid:48) ≥0. Thisconstrainguaranteesthatthewarpingpathwillnotrollbackonitself.
i i j j
Thatis,bothindexesiandj eitherstaythesameorincrease(theyneverdecrease).
2. Continuity: Any two adjacent elements of the warping path W,w = (w ,w ) and w = (w(cid:48),w(cid:48)), follow the
k i j k+1 i j
inequalities,w −w(cid:48) ≤ 1andw −w(cid:48) ≤ 1. Thisconstraintguaranteesthatthewarpingpathadvancesonestepata
i i j j
time. Thatis,bothindexesiandj canonlyincreasebyatmost1oneachstepalongthepath.
3. Boundary: The warping path starts from the top left corner w = (1,1) and ends at the bottom right corner w =
1 k
(n,m). Thisconstraintguaranteesthatthewarpingpathcontainsallpointsofbothtimeseries.
Figure2.Globalconstraint(SakoeChibaBand),whichlimitsthewarpingscope. Thediagonalgreen
areascorrespondtothewarpingscopes.
Althoughtherearealargenumberofwarpingpathsthatsatisfyalloftheaboveconstraints,DTW isdesignedtofindthe
onethatminimizesthewarpingcost(distance). Figures1(a)and1(b)demonstrateanexampleofhowtwotimeseries(Sand
Q)arewarpedandthewaytheirdistanceiscalculated. Thecircledcellsshowtheoptimalwarpingpath,whichcrossesthe
gridfromthetopleftcornertothebottomrightcorner. TheDTW distancebetweenthetwotimeseriesiscalculatedbased
onthisoptimalwarpingpathusingthefollowingequation:
(cid:113)
 (cid:80)Kk=1Wk
DTW(S,Q)=min (4)
K

TheK inthedenominatorisusedtonormalizedifferentwarpingpathswithdifferentlengths.
SincetheDTW hastopotentiallyexamineeverycellinthewarpingmatrix,itsspaceandtimecomplexityisO(nm).
4 GlobalConstraint(BandDTW)
ThereareseveralmethodsthataddglobalconstraintsonDTWtoincreaseitsspeedbylimitinghowfarthewarpingpath
maystrayfromthediagonalofthewarpingmatrix(Tappert&Das1978,Berndt&Clifford1994,Myersetal.1980). Inthis
paperweuseSakoe-ChibaBand(henceforth,werefertoitasBandDTW)Sakoe&Chiba(1978)whencomparingwithour
proposedalgorithm(Figure2). BandDTWusedtospeeduptheDTW byaddingconstraintswhichforcethewarpingpathto
liewithinabandaroundthediagonal;iftheoptimalpathcrossestheband,theDTWdistancewillnotbeoptimal.
5 DivideandConquerTechnique(DC)
Intheprevioussection,wehaveshownhowtocomputetheoptimalalignmentusingthestandardDTWtechniquebetween
twotimeseries. InthissectionwewillshowanothertechniquethatusesaDivideandConquerheuristic,henceforthwerefer
toitas(DC),proposedbyHirschberg(1975). DCisadynamicprogrammingalgorithmusedtofindtheleastcostsequence
alignment between two strings. The algorithm was first introduced to solve the Longest Common Subsequence (LCSS)
(Hirschberg1975). Algorithm2givesahighleveldescriptionoftheDCalgorithm. Likeinthestandardsequencealignment,
theDC algorithmhasO(mn)timecomplexitybutO(m+n)spacecomplexity,wheremandnarethelengthsofthetwo
Algorithm2DC:DivideandConquertechnique.
Input: S: Sequenceoflengthn,Q: Sequenceoflengthm.
Output: DTWdistance.
1: Divide-Conquer-Alignment(S,Q)
2: n⇐|S|
3: m⇐|Q|
4: Mid⇐(cid:100)m/2(cid:101)
5: ifn≤2orm≤2then
6: ComputeoptimalalignmentusingstandardDTW
7: else
8: f ⇐ForwardsSpaceEfficientAlign(S,Q[1:Mid])
9: g ⇐BackwardsSpaceEfficientAlign(S,Q[Mid:m])
10: q ⇐indexthatminimizingf(q,Mid)+g(q,Mid)
11: Add(q,Mid)toglobalarrayP
12: Divide-Conquer-Alignment(S[1:q],Q[1:Mid])
13: Divide-Conquer-Alignment(S[q:n],Q[Mid:m])
14: endif
15: return P
sequencesbeingaligned. WewillbeusingAlgorithm2alongwithFigure3toexplainhowDC works. Intheexamplewe
usetwosequencesS = [3,4,5,3,3]andQ = [1,2,2,1,0]todeterminetheoptimalalignmentbetweenthem. Thereisonly
oneoptimalalignmentforthisexample(Figure3(e)), whereshadedcellsaretheoptimalwarpingpath. TheDC algorithm
worksasfollows:
1. ItfindsthemiddlepointinQwhichisMid = |Q|/2, (Figure3(a)). Thishelpstofindthesplitpointwhichdivides
the warping matrix into two parts (sub-problems). A forward space efficiency function (Line 8) uses S and the first
cut of Q = [1,2,2], then a backward step (Line 9) uses S and Q = [2,1,0] (Figure 3(a)). Then by adding the
lastcolumnfromtheforwardandbackwardstepstogetherandfindingtheindexoftheminimumvalue,theresultant
column indicates the row index that will be used along with the middle point to locate the split point (shaded cell
in Figure 3(a)). Thus, the first split point is D(4,3). At this stage of the algorithm, there are two sub-problems; the
alignmentofS =[3,4,5,3]withQ=[1,2,2]andofS =[3,3]withQ=[2,1,0].
2. DCisrecursivealgorithm,eachcallsplitstheproblemintotwoothersub-problemsifbothsequencesareoflength>2,
otherwiseitcallsthestandardDTW tofindtheoptimalpathforthatparticularsub-problem. Intheexample,thefirst
sub-problemwillbefedtoLine12whichwillfindanothersplitpoint,becausebothinputsequencesareoflength>2.
Figure3(b)showshowthenewsplitpointisfound. Figure3(c)showsthetwosplitpoints(shadedcells)whichyield
tohavesub-problemsofsequencesoflength≤2. InthiscaseDTW willbeusedtofindtheoptimalalignmentforeach
sub-problem.
3. TheDCalgorithmfindsthefinalalignmentbyconcatenatingtheresultsfromeachcallofthestandardDTW.
TheexampleinFigure3clarifiesthattheDCalgorithmdoesnotgivetheoptimalwarpingpath. Figures3(d)and(e)show
thepathsobtainedbytheDCandDTW algorithms,respectively.
DC does not yield the optimal path as it goes into infinite recursion because of how it calculates the middle point. DC
calculatesthemiddlepointasfollows:
Therearetwoscenarios: first,whenthemiddlepoint(Algorithm2Line4)isfloored(Mid = (cid:98)m/2(cid:99))andsecondwhen
itisroundedup(Mid=(cid:100)m/2(cid:101)). Thefirstscenariocausesinfiniterecursion,sincethesplitfromthepreviousstepgivesthe
samesub-sequences(i.e.,thealgorithmkeepsfindingthesamesplitpoint). ThesecondscenarioisshowninFigures3(a-d),
whichclearlyconfirmsthatthefinaloptimalpathisnotthesameastheoneretrievedbythestandardDTW 1.ThefinalDTW
distanceisdifferentaswell. TheshadedcellsinFigures3(d)and(e)showthatbothwarpingpathsaredifferent.
1Itshouldbenotedthatourexamplehasonlyoneoptimalpaththatgivestheoptimaldistance.
Middle point Middle point
Q Q Q Q
1 2 2 2 1 0 1 2 2 2
3 4 5 6 33 27 42 68 3 3 4 5 20 15 15 3
S S
S 4 13 8 9 35 26 38 59 4 4 13 8 22 14 14 4
5 29 17 17 39 22 29 43 5 5 29 17 27 10 10 5
3 33 18 18 32 14 13 18 3 3 33 18 20 2 1 3
3 37 19 19 33 14 13 9 3 Forward cut Backward cut
Forward cut Backward cut
(b) Finding second split point.
(a) Finding first split point.
Q Q Q
1 2 2 1 0 1 2 2 1 0 D 1 2 2 1 0
3 3 4 5 6 3 4 5 6 10 19
4 4 9 9 4 13 8 9 15 26
S S S
5 5 18 18 5 29 17 17 25 40
3 3 19 19 23 32 3 33 18 18 21 30
3 3 20 23 32 3 37 19 19 22 30
(c) Split points (sub-problems).
(d) Final optimal path using DC. (e) Final optimal path using standard DTW.
Figure3.AnexampletoshowthedifferencebetweenthestandardDTWandtheDCalgorithm.
6 SparseDynamicProgrammingApproach
In this section, we outline the main principles we use in SparseDTW and follow up with an illustrated example along
with the SparseDTW pseudo-code. We exploit the following facts in order to reduce space usage while avoiding any re-
computations:
1. Quantizingtheinputtimeseriestoexploitthesimilaritybetweenthepointsinthetwotimeseries.
2. Using a sparse matrix of size k, where k = n×m in the worst case. However, if the two sequences are similar, k
<<n×m.
3. Thewarpingmatrixiscalculatedusingdynamicprogrammingandsparsematrixindexing.
6.1 Key Concepts
Inthissectionweintroducethekeyconceptsusedinouralgorithm.
Definition1(SparseMatrixSM) isamatrixthatispopulatedlargelywithzeros.Itallowsthetechniquestotakeadvantage
ofthelargenumberofzeroelements. Figure4(a)showstheSM initialstate. SM islinearlyindexed,Thelittlenumbers,in
thetopleftcornerofSM’scells,representthecellindex. Forexample,theindicesofthecellsSM(1,1)andSM(5,5)are
1and25,respectively.
Definition2(LowerNeighbors(LowerNeighbors)) acellc ∈ SM hasthreelowerneighborswhicharethecellsofthe
indices(c−1),(c−n),and(c−(n+1))(wherenisthenumberofrowsinSM). Forexample,thelowerneighborsofcell
SM(12)areSM(6),SM(7)andSM(11)(Figure4(a)).
Q Q Q
SM 1 2 2 1 0 SM 1 2 2 1 0 SM 1 2 2 1 0
3 1 B 6 B 11B 16B 21B 3 4 B B 4 9 3 4 B B 4 9
4 2 B 7 B 12B 17B 22B 4 9 4 4 9 16 4 13 4 4 9 16
S S S
5 3 B 8 B 13B 18B 23B 5 16 9 9 16 B 5 29 9 9 16 25
4 9 14 19 24
3 B B B B B 3 4 B B 4 9 3 33 B B 4 9
3 5 B 10B 15B 20B 25B 3 4 B B 4 9 3 37 1 1 4 9
(a) SM initially blocked [B]. (b)SM after unblocking the optimal cells. (c) Unblocking upper neighbor
Values are Euclidean distances. (Shaded cell).
Q Q Q
SM 1 2 2 1 0 SM 1 2 2 1 0 D 1 2 2 1 0
3 4 B B 4 13 3 4 B B 4 13 3 4 5 6 10 19
4 13 8 12 13 20 4 13 8 12 13 20 4 13 8 9 15 26
S S
S 5 29 17 17 28 38 5 29 17 17 28 38 5 29 17 17 25 40
3 33 B B 21 30 3 33 B B 21 30 3 33 18 18 21 30
3 37 34 35 25 30 3 37 34 35 25 30 3 37 19 19 22 30
(d) Constructing SM. (e) Final optimal path using SparseDTW. (f) Final optimal path using standard DTW.
Figure4.AnexampleoftheSparseDTWalgorithmandthemethodoffindingtheoptimalpath.
Definition3(UpperNeighbors(UpperNeighbors)) acellc ∈ SM hasthreeupperneighborswhicharethecellsofthe
indices(c+1),(c+n),and(c+n+1)(wherenisthenumberofrowsinSM). Forexample,theupperneighborsofcell
SM(12)areSM(13),SM(17)andSM(18)(Figure4(a)).
Definition4(BlockedCell(B)) a cell c ∈ SM is blocked if its value is zero. The letter (B) refers to the blocked cells
(Figure4(a)).
Definition5(Unblocking) Givenacellc∈SM,ifSM(c)’supperneighbors(SM(c+1),SM(c+n),andSM(c+n+1))
areblocked,theywillbeunblocked. UnblockingisperformedbycalculatingtheEucDistforthesecellsandaddingthemto
SM. In other words, adding the distances to these cells means changing their state from blocked (B) into unblocked. For
example,SM(10)isablockedupperneighborofSM(5),inthiscaseSM(10)needstobeunblocked(Figure4(c)).
6.2 SparseDTW Algorithm
Algorithm3takesRes,theresolutionparameterasaninputthatdeterminesthenumberofbinsas 2 . Reswillhaveno
Res
impactontheoptimality. Wenowpresentanexampleofouralgorithmtoillustratesomeofthehighlightsofourapproach:
Westartwithtwosequences:
S =[3,4,5,3,3]andQ=[1,2,2,1,0].
InLine1,wefirstquantizethesequencesintotherange[0,1]usingEquation5:
Sk−min(Sk)
QuantizedSeqk = i . (5)
i max(Sk)−min(Sk)
WhereSk denotestheithelementofthekthtimeseries. Thisyieldsthefollowingsequences:
i
S(cid:48) =[0,0.5,1.0,0.0,0.0]andQ(cid:48) =[0.5,1.0,1.0,0.5,0]
InLines4to7wecreateoverlappingbins,governedbytwoparameters: bin-widthandtheoverlappingwidth(whichwe
refertoastheresolution). Itisimportanttonotethatthesetwoparametersdonotaffecttheoptimalityofthealignmentbut
do have an affect on the amount of space utilized. For this particular example, the bin-width is 0.5. We thus have 4 bins
whichareshowninTable1.
BinNumber Bin Indices Indices
(B ) Bounds ofS(cid:48) ofQ(cid:48)
k
1 0.0-0.5 1,2,4,5 1,4,5
2 0.25-0.75 2 1,4
3 0.5-1.0 2,3 1,2,3,4
4 0.75-1.25 3 2,3
Table1.Binsbounds,whereB isthekth bin.
k
Ourintuitionisthatpointsinsequenceswithsimilarprofileswillbemappedtootherpointsinthesamebinorneighboring
bins. Inwhichcasethenon-defaultentriesofthesparsematrixcanbeusedtocomputethewarpingpath. Otherwise,default
entriesofthematrixwillhavetobe“opened”,reducingthesparsityofthematrixbutneversacrificingtheoptimalalignment.
In Lines 3 to 13, the sparse warping matrix SM is constructed using the equation below. SM2 is a matrix that has
generallyfewnon-zero(or“interesting”)entries. Itcanberepresentedinmuchlessthann×mspace,wherenandmare
thelengthsofthetimeseriesS andQ,respectively.
(cid:26)
EucDist(S(i),Q(j)) ifS(i)andQ(j)∈B
SM(i,j)= k (6)
B otherwise
We assume that SM is linearly ordered and the default value of SM cells are zeros. That means the cells initially are
Blocked (B) (Figure 4(a)). Figure 4(a) shows the linear order of the SM matrix, where the little numbers on the top left
cornerofeachcellrepresenttheindexofthecells. InLine6and7,wefindtheindexofeachquantizedvaluethatfallsinthe
binbounds(Table1column2,3and4). TheInequality7isusedinLine6and7tofindtheindicesofthedefaultentriesof
theSM.
LowerBound≤QuantizedSeqk ≤UpperBound. (7)
i
Where LowerBound and UpperBound are the bin bounds and QuantizedSeqk represents the quantized time series
i
whichcanbecalculatedusingEquation5.
Lines8to12areusedtoinitializetheSM. ThatisbyjoiningallindicesinidxS andidxQtoopencorrespondingcells
inSM. Afterunblocking(opening)thecellsthatreflectthesimilaritybetweenpointsinbothsequences,theSM entriesare
showninFigure4(b).
Lines14to22areusedtocalculatethewarpingcost.InLine15,wefindthewarpingcostforeachopencellc∈SM (cell
cisthenumberfromthelinearorderofSM’scells)byfindingtheminimumofthecostsofitslowerneighbors,whichare
[c−1,c−n,c−(n+1)](blackarrowsinFigure4(d)showthelowerneighborsofeveryopencell).Thiscostisthenaddedto
thelocaldistanceofcellc(Line17). TheabovestepissimilartoDTW,however,wemayhavetoopennewcellsiftheupper
neighborsatagivenlocalcellc∈SM areblocked. Theindicesoftheupperneighborsare[c+1,c+n,c+n+1],where
n is the length of sequence S (i.e., number of rows in SM). Lines 18 to 21 are used to check always the upper neighbors
ofc ∈ SM. Thisisperformedasfollows: ifthe|UpperNeighbors| = 0foraparticularcell, itsupperneighborswillbe
unblocked. This is very useful when the algorithm traverses SM in reverse to find the final optimal path. In other words,
unblockingallowsthepathtobeconnected. Forexample,thecellSM(5)hasoneupperneighborthatiscellSM(10)which
isblocked(Figure4(b)),thereforethiscellwillbeunblockedbycalculatingtheEucDist(S(5),Q(2)). Thevaluewillbeaddto
theSM whichmeansthatcellSM(10)isnowanentryinSM (Figure4(c)). AlthoughunblockingaddscellstoSM which
meansthenumberofopencellswillincrease,buttheoverlappinginthebinsboundariesallowstheSM’sunblockedcellsto
beconnectedmostlythatmeanslessnumberofunblockingoperations. Figure4(d)showsthefinalentriesoftheSM after
calculatingthewarpingcostofallopencells.
2IftheEuclideandistance(EucDist)betweenS(i)andQ(j)iszero,thenSM(i,j) = −1,todistinguishbetweenablockedcellandanycellthat
representszerodistance.
Lines 23 to 32 return the warping path. hop initially represents the linear index for the (m,n) entry of SM, that is the
bottomrightcornerofSM inFigure4(e).Startingfromhop=n×mwechoosetheneighbors[hop−n,hop−1,hop−(n+
1)]withminimumwarpingcostandproceedrecursivelyuntilwereachthefirstentryofSM,namelySM(1)orhop=1. It
isinterestingthatwhilecalculatingthewarpingpathweonlyhavetolookattheopencells,whichmaybefewerinnumber
than3. Thispotentiallyreducestheoveralltimecomplexity.
Figure 4(e) demonstrates an example of how the two time series (S and Q) are warped and the way their distance is
calculatedusingSparseDTW.Thefilledcellsshowtheoptimalwarpingpath,whichcrossesthegridfromthetopleftcorner
tothebottomrightcorner. ThedistancebetweenthetwotimeseriesiscalculatedusingEquation4. Figure4(f)showsthe
standardDTW wherethefilledcellsaretheoptimalwarpingpath. Itisclearthatbothtechniquesgivetheoptimalwarping
pathwhichwillyieldtheoptimaldistance.
Algorithm3SparseDTW:Sparsedynamicprogrammingtechnique.
Input: S: Timeseriesoflengthn,Q: Timeseriesoflengthm,andRes.
Output: OptimalwarpingpathandSparseDTW distance.
1: [S(cid:48),Q(cid:48)]⇐Quantize(S,Q)
2: LowerBound⇐0,UpperBound⇐Res
3: forall0≤LowerBound≤1− Res do
2
4: IdxS ⇐find(LowerBound≤S(cid:48) ≤UpperBound)
5: IdxQ⇐find(LowerBound≤Q(cid:48) ≤UpperBound)
6: LowerBound⇐LowerBound+ Res
2
7: UpperBound⇐LowerBound+Res
8: forallidxi ∈IdxS do
9: forallidxj ∈IdxQdo
10: AddEucDist(idxi,idxj)toSM {WhenEucDist(idxi,idxj)=0,SM(i,j)=−1.}
11: endfor
12: endfor
13: endfor
{Note: SM islinearlyindexed.}
14: forallc∈SM do
15: LowerNeighbors⇐{(c−1),(c−n),(c−(n+1))}
16: minCost⇐min(SM(LowerNeighbors)){SM(LowerNeighbors)=-1meanscost=0.}
17: SM(c)⇐SM(c)+minCost
18: UpperNeighbors⇐{(c+1),(c+n),(c+n+1)}
19: if|UpperNeighbors|==0then
20: SM ∪EucDist(UpperNeighbors)
21: endif
22: endfor
23: WarpingPath⇐Φ
24: hop⇐SM(n×m){LastindexinSM.}
25: WarpingPath∪hop
26: whilehop(cid:54)=SM(1)do
27: LowerNeighbors⇐{(hop−1),(hop−n),(hop−(n+1))}
28: [minCost,index]⇐min[Cost([LowerNeighbors])
29: hop⇐index
30: WarpingPath∪hop
31: endwhile
32: WarpingPath∪SM(1)
33: return WarpingPath,SM(n×m)