BILCO: An Efficient Algorithm for Joint Alignment of Time Series
Authors: Xuelong Mi
Mengfan Wang
Alex Bo-Yuan Chen
Jing-Xuan Lim
Yizhi Wang
Misha Ahrens
Guoqiang Yu
Venue: NeurIPS 2022 Accept (2022)
Keywords: Joint alignment
graphical time warping
bidirectional-pushing strategy
linear component operation
BILCO

Abstract: Multiple time series data occur in many real applications and the alignment among them is usually a fundamental step of data analysis. Frequently, these multiple time series are inter-dependent, which provides extra information for the alignment task and this information cannot be well utilized in the conventional pairwise alignment methods. Recently, the joint alignment was modeled as a max-flow problem, in which both the profile similarity between the aligned time series and the distance between adjacent warping functions are jointly optimized. However, despite the new model having elegant mathematical formulation and superior alignment accuracy, the long computation time and large memory usage, due to the use of the existing general-purpose max-flow algorithms, limit significantly its well-deserved wide use. In this report, we present BIdirectional pushing with Linear Component Operations (BILCO), a novel algorithm that solves the joint alignment max-flow problems efficiently and exactly. We develop the strategy of linear component operations that integrates dynamic programming technique and the push-relabel approach. This strategy is motivated by the fact that the joint alignment max-flow problem is a generalization of dynamic time warping (DTW) and numerous individual DTW problems are embedded. Further, a bidirectional-pushing strategy is proposed to introduce prior knowledge and reduce unnecessary computation, by leveraging another fact that good initialization can be easily computed for the joint alignment max-flow problem. We demonstrate the efficiency of BILCO using both synthetic and real experiments. Tested on thousands of datasets under various simulated scenarios and in three distinct application categories, BILCO consistently achieves at least 10 and averagely 20-folds increase in speed, and uses at most 1/8 and averagely 1/10 memory compared with the best existing max-flow method. Our source code can be found at https://github.com/yu-lab-vt/BILCO.

BILCO: An Efficient Algorithm for Joint Alignment of Time Series
Xuelong Mi', MMngfan Wang', Alex Bo-Yuan Chen', Jing-Xuan Lim?, Yizhi Wang', Misha Ahrens", Guoqiang Yu²
1 Dept. of Electrical and Computer Engineering, Virginia Tech
Howard Hughes Medical Institute, Janelia Research Campus
{mix118, mengfan yzwang yugh@vt. edu 2 2 limj2, ahrensm anelia. hhmi org
Abstract
Multiple time series data occur in many real applications and the alignment among them is usually fundamental step of data analysis. Frequently, these multiple time series are inter-dependent, which provides extra information for the alignment task and this information cannot be well utilized in the conventional pairwise alignment methods. Recently, the joint alignment was modeled as a max-flow problem, which both the profile similarity between the aligned time series and the distance between adjacent warping functions are jointly optimized. However, despite the model having elegant mathematical formulation and superior alignment accuracy, the long computation time and large memory usage, due to the use of the existing ceneral-puuposs max- flowwlow algorithms, limit significantly its well-deserved wide use. In this report, we present Bidirectional pushing with Linear Component Operations (BILCO), novel algorithm that solves the joint alignment max-flow problems efficiently and exactly. We develop the strategy of linear component operations that integrates dynamic programming technique and the push-relabel approach. This strategy is motivated by the fact that the joint alignment max- flow problem a generalization of dynamic time warping (DTW) and numerous individual DTW problems are embedded. Further, hidirectional-pushigg strategy is proposed introduce prior knowledge and reduce unnecessary computation, by leveraging another fact that good initialization can be easily computed for the joint alignment max-flow problem. We demonstrate the efficiency BILCO using both synthetic and real experiments. Tested thousands of datasets under various simulated scenarios and in three distinct application categories, BILCO consistently achieves at least 10 and averagely 20-folds increase in speed, and uses at most 1/8 and averagely 1/10 memory compared with the best existing max-flow method.
Our source code can be found at
Introduction
Time series data appear naturally in wide range of fields, such as motion capture, speech recognition, and bioinformatics. As sequences are usually not directly comparable due to possible delay and distortion, alignment among them is prerequisite step before further analysis. The most popular alignment technique, dynamic time warping (DTW) (3), finds an optimal match between two sequences linear time using dynamic programming (DP). However, DTW and its variants 2 only compare a single pair of time series, while real data often contain structural information. For example, in time-lapse mieroscopy data, the pixel intensities are recorded over time in a 2D grid, and adjacent pixels tend to have similar temporal patterns. Thus, joint alignment of more than a single pair is expected to achieve better performance by taking advantage of dependency from the structure.
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

[+]
Single pair alignment
[d] Joint alignment
Xr
Binorolar stereo [eeftiitee
p0
graph for tst and
Patrwine aligmment
(b)
destinatioo
Joint algnment
DTW graph
GTW graph
subgraphs Cross edges
Degth
Fearnie 10
Frame 18
Framie
John slagranes
Aomer Lime
graph/GTW subgraph
(a-c) The alignment between sequences x and y. DTW graph of alignment shown in (b) and its dual graph in (c). In (b), if warping path crosses (k1, k2). 20[c]] and y(ky] are matched. Each warping path in (b) is dual to a corresponding cut (c). (d) GTW graph construction for time series arranged on 2D grid. Details and Pa (top pancl), their corresponding cuts segmented the dual graph into the sink side (red) and source side (blue) (bottom panel). dist P..: P.) (the green shadow in the top panel) is proportional to the number of cross edges linking source side and sink side. (e) application of pairwise alignment and joint alignment calculating signal propagation. (f) application to extracting depth information (data from (6)).
Figure I:
shown for two adjacent pairs (m, n.) Given warping functions Pm
For long time, it was not clear how to incorporate the structural information in a principled way. Various heuristic tricks were developed [3 19I. In a recent theoretical breakthrough, the joint alignment problem was elegantly modeled as max-flow problem of graphical time warping (GTW) graph by optimizing both the time series pairwise similarity and the distance between warping functions (22). GTW was shown to have superior alignment accuracy many applications, such as brain activity analysis (21) and liquid chromatography- mass spectrometry (LC-MS) proteomics (23, where warping functions represent the delay of the real signal at each time point relative to the reference time series and adjacent time series should propagate similarly. The binocular stereo vision formulation in [ is also a special case of GTW, where warping functions denote the depth information contained in two views and it is natural to assume the depths are close in adjacent positions. By adding penalties on warping function dissimilarity, GTW utilizes the structural information and could obtain better alignment performance, the long computation time and large memory usage severely limit its potential broad applications, because GTW was solved by existing general purpose max-flow algorithms and the graph is often huge. Indeed, for a typical dataset with 5000 pairs of time series, if each time series contains 200 time points, all popular methods, such incremental breadth- first search (IBFS) [8] Hochbaum's 's pseudoflow (HPF) 10, Boykov-Kolmogord (BK) max-flow [4] [2], and highest-label push-relabel (HIPR)
shown in Fig Fig (e) However,
5, cost hours and more than 100 gigabytes memory.
In this paper, we identify two important properties of the joint alignment max-flor problem and show that they can be leveraged to design novel algorithm with improvements in both speed and memory efficiency. Specifically, first, joint alignment is a generalization of pairwise alignment and numerous individual DTW problems are embedded, as shown in Fig [d)()). If the dependency is neglected, it is reduced to multiple independent DTW problems, which can be solved in linear time through DP. Although the dependency makes the problem more complex, the property of the DTW problem can still be utilized. Second, a coarse approximate solution to joint alignment can be readily estimated in many applications. Such prior knowledge can be incorporated accelerate the max-flow algorithms.
Taking advantage of these two properties, we develop Bldirectional pushing with I. Linear Component Operations (BILCO), a novel algorithm that solves the joint alignment mlx-flow problem exactly and

---

more efficiently. The algorithm consists of two major strategies: (a) By utilizing the first property. we propose Excess pushing with Linear Component Operations (ELCO) that integrates DP and the push-relahel approach [9] Each component is defined as small subgraph including connected nodes and their related edges. The components are automatically and adaptively determined. We show that all operations on each component can be achieved in linear time by combining DP and subgraph properties. With linear component operations, ELCO can achieve higher efficiency than the generic push-relabel approach. (b) By leveraging the second property, we design bidirectional-puse strategy to utilize prior knowledge as initialization. The strategy simplifies the original problem into two smaller sub-pree sub-problems with opposite pushing directions. In each sub-problem, ELCO is applied, but the separation two pushing directions dramatically reduces unnecessary computation.
BILCO has the same theoretical time complexity as the best popular methods such as HIPR, but it provides a significant empirical efficiency boost without sacrificing the accuracy in the task of joint alignment. Its effectiveness is evidenced through both synthetic and real experiments. Compared with IBFS, HPF, BK, and HIPR on thousands of datasets under various simulated scenarios and three real applications in distinct categories, BILCO improves the speed by 10-50 folds and costs averagely 1/10 memory relative to the best peer method.
2 Problem formulation
The joint alignment problem can be formulated 22
min Pa.
cost( +
Σ
dist(Pm my Pr))
(1)
(min)
where Pa denotes the warping path for the nthh time series pair, N is the total number of time series pairs that are jointly aligned, cost(Px) ) is the alignment cost of the пın time series pair, dist(Pm. Pa) the warping path distance defined by the area of the region bounded by Pm. and Pn- Neib is the set of pair indices (m, n.) representing the adjacent time series, and is the hyperparameter to balance the alignment cost term and the distance term. For example, for 2D grid time-series data, Neib may include the pair of neighboring pixels, and K represents a prior similarity between the pixels.
As shown in Fig Equation ①1) can be converted flow network and solved by finding the min-eut of it. The constructed graph, GTW graph, consists of N GTW subgraphs (G" = (V", Fr)] = 1, V.....N. and cross edges Feross with capacity 5 (c)). The edges within GTW subgraph are called Ewithin- For convenience, in the following content, we refer to "GTW subgraph" as "subgraph". Each subgraph G" is dual DTW graph, which represents the warping between pair of time series. The cut within a subgraph is dual to a warping path of the time series pair [I] (Fig (c)) and thus the min-cut of one DTW graph can be solved in linear time through DP by finding the shortest path. The cross edges constrain the difference between cuts in neighboring subgraphs, corresponding to the distance term in Equation (T) (Fig dd). To ensure the monotonicity and continuity of warping paths, the capacities of reverse edges in each subgraph are set infinite [22].
It is worth noting that GTW framework is flexible, broadly applicable, and be used building block in solving many problems while utilizing structural information. The assumed neighbor structure and similarity strength be application-specfii different pairs of warping functions, GTW can set different hyperparameters k. For simplicity, in this paper we use the same hyperparameter x for all warping function pairs in the following.
user-des,gned. For example, for
3 Method
BILCO contains two major parts, ELCO and bidirectional pushing, as discussed in Sections 3.1 and respectively. Based on the initialization of warping functions, the bidirectional pushing strategy converts a max-flow problem into two sub-problems and we solve each of them with ELCO. The integration leads to BILCCO (Section 3.3), followed by the analysis of time complexity and memory usage in Section 3.4 and 3.5 All lemmas and theorems are proved in the supplementary.
3.2
3

---

(6)
Aeiter
water
water
Sink
G2
(c)
(e)
Inactive
Discharge through tubes
Active sgain,
New cut
Active
component
G2
G2
Subgraph/wwterr tank
cross edges/tube
Fartttio/outt
Passable adges
Blocked tube O Unblocked tube
An analogy joint alignment max-flow problem. Water, tube, tank, and partition represent excess, cross edge, subgraph, and cut, respectively. (a) In the initial stage, the source has pushed all water it water tanks. (b) "Drain" operation: water in the tanks is drained sink. (c) "Discharge" operation: water flows from an active component to neighbor component in another water tank, where water needs concentrate at the connecting tube first. (d) After (c), the previously inactive component in G2 becomes active and water can be drained again. Due to the concentration step in (c), (e) One intermediate stage of tank, where multiple components are segmented in a subgraph and some tubes are blocked.
Figure 2:
there one new partition in G that further segments it into three components.
3.1 Excess pushing with linear component operations (ELCO)
Since GTW graph contains many subgraphs dual DTW graphs and each DTW graph can be linearly solved by DP, we hope the DTW graph and DP algorithm can be exploited in our approach. While the cross edges between subgraphs prohibit the direct use DP, we found that the component as defined in Subsection 3.1.1 is the maximum subgraph that can allow the use of DP. By combining the newly designed graph conversion strategy (Fig and DP, we establish that the operations on components, "Drain" and "Discharge" (defined in Subsections 3.1.1). 1), can indeed be linearly implemented as in Subsections 3.1.2 and 3.1.3. To guarantee global optimality, we design new labeling function as in Subsection 3.1.4 1.4 and use that to guide our component operations.
3.1.1 The basic unit operation component, not node subgraph
To exploit the property of each subgraph, ELCO borrows the idea from the generic push-relabe algorithm, whose operations are localized. ELCO allows the existence of excess, which is the surplus between the entering flow and outgoing flow of one node. The nodes that carry excess are regarded as "active"". Based on the location of flow exchanges, "Drain" pushes excess to the sink through Ewithin. which is within- within-subbraphh operation. While "Discharge" is cross-subgraph operation that pushes excess neighbor subgraphs through Ecrozs. By alternatively executing these two operations, more and more excess can be sent to sink (Fig(2). The spread of excess leads to many edges being saturated, resulting multiple cuts the same subgraph. These cuts segment the subgraph into different components Ee)) Formally, we define one component as a subset GTW subgraph bounded by two adjacent cuts. ELCO uses component, instead of node or subgraph, as the basic operation unit, because it is the largest possible unit to use DP to push flows.
define two operations, "Drain" and "Discharge".
Lemma 1: two nodes are in the same component, then there is at least one path linking them.
By Lemma if one node in component is active, the excess component. Thus, the whole component can be seen as one unit. On the contrary, cuts block the flow across components in the same subgraph, it is inapplicable to take a whole subgraph as one unit.
spread anywhere in the same
Using the component as the basic operation unit, ELCO can be described as follows: Initially, each subgraph is one active component. "Drain" operation sends maximum excess directly to the sink through Ewithin- A new cut will be generated, which blocks some excess in segmented components near the source. Then it is the turn of "Discharge" to seek opportunities to move excess across subgraph toward components that can execute "Drain".. By executing these two component operations alternatively, the global max-flow can be achieved. An illustrative example is shown

---

(b)
(c)
(d)
(o)
sink
sink
lat b
tar ci
source
SOUICE
6+- Entering flow
Ourgoing flow
Current cut
Extra cost for current out
0+C
Figure 3: An example of graph conversion to incorporate flows Ecross into Emithhnn (a) Subgraph with flows on Eeross in the first column. (b) Entering flow and outgoing flow are like extra edges linking to source and sink. (c) Purple cut would segment the subgraph into source side and sink side, where extra edges lead to counting two more cuts the cost. (d) Costs are added to all paths/cuts below the edge linking to the source and are added to all paths/cuts above the edge linking to the sink. (c) The first column of the converted planar DTW graph. (f) The converted GTW subgraph.
in the supplementary. In the following, we will show "Drain" and "Discharge" operations can be implemented in linear time complexity.
3.1.2 Linear time implementation of "Drain" component operation
Regarding "Drain" operation, we only need to calculate the min-cut of one active component that links to the sink, then all the excess above the new cut is drained and the part below the new cut is identified new active component (Fig (b))). As mentioned, the min-cut can be quickly solved by DP on its dual graph. However, with flows on Epross: DP cannot be applied directly. The graph is not planar anymore (ab)) and its dual graph does not exist. Here we design linear-time graph conversion strategy as shown in Fig 3see (see algorithms in the supplementary), which can incorporate known flows Eerose into Ewithin and obtain equivalent planar graph where DP can still be applied. Combining linear conversion strategy with DP, the "Drain" component operation can be implemented in linear time.
Lemma 2: The corresponding cut values/path costs before and after graph conversion are the same.
3.1.3 Linear time implementation of "Discharge" component operation
To discharge component, try push excess out each node much possible to make sure the component has sent all the excess it can since excess can flow anywhere in the same component.
Lemma 3:
The amount of maximum possible excess on node is the min-cut value in residual graph by taking v as sink.
Working on the modified DTW graph. DP can still be applied to calculate such excess by finding the difference between the shortest path below v and the shortest path of the whole graph. The difference can represent the min-cut value mentioned in Lemma 3. To implement linear-complexity component discharge operation, here we use another layer of DP that reuses the stored dynamic matrices that save the shortest distance recursively. The detailed algorithm of "Discharge" is given in the supplementary.
3.1.4 A new labeling function guide excess approaching the sink across components
We design a new labeling function to guide excess approaching across components from high label to low label. Different from the in generic aush-relabel [9], our labeling function implies the distance from one component to sink rather than from the node. Since the cut blocks the connection of components in the same subgraph, only Ecross is counted. Derived from such distance, we define our labeling function d + to be valid if for all residual edges,
(v, w) E Ewithin d(u) S d(w) +1 (v, w) E Ecross
d(v) S d(w)
(2a)
(2b)
With such definitions, we have,
Lemma 4: All nodes in the same component have the same label.
Lemma 5: The new labeling function ( consistent with generic push-relabel labeling function treating each component one unit.
5

---

Algorithm ELCO
R2 G², d(Ra) 0, cul Drain (R:). Split (Ri, cut) for -,22.-, 2,... N while active region exists do
Initialization
Pick active region with highest label, R2 d( R.) = 0 then
Highest-label selection rule & sink component First type operation and update I- No direct edge linking to sink
cut Drain Split (Ri, cut)
else
if JR that R; + R2 and d(R;) d(R,) cut = Discharge( (R,). Split (R,, cut)
then
>
Second type operation and update
else
Relabel(R,)
> Relabel
Gap-relahel heuristic
Relabel may result label gap
end if end if end while
Algorithm 2 Relabel( R), RE G
AminCross = min (d(())) + Rg. R2 & G*}
D Lowest reachable label in neighbor subgraphs Lowest reachable label in same subgraph
AminWithin min R2) [ Rz. R2 G*)
if dminnross < AmminWthie then
d(R) UminCross +
else
A Consider the validity
d(R) unnnW
R = R R R RR,RR, RE RE G", G", E G") into R
Merge
Relabel(R)
end if'
Lemma 4 suggests that we can set the same label to one component and all the nodes in it. And lemma implies that ELCO can be seen as an alternative push-relabel approach that uses component as the operation unit, and all the statements in the generic push-relabel approach [9] hold on the component level. Since the validity of the labeling function is maintained in ELCO (see the proof in supplementary), the optimal solution can be achieved in polynomial component operations.
The ELCO is outlined in Algo and Algo Details and illustrative example can be found in the supplementary. R is the symbol of component. R2 + R2 if for o E R4 and w E R2 that (v, w) has positive residual capacity. Highest-label selection rule is applied, and both global-relabel heuristic and gap-relabel heuristic in 9 are utilized for acceleration.
3.2 Bidirectional-pushing strategy
Generic push-relabel methods, and ELCO alike, attempts to push as much as excess to the sink. If part of the excess fails to reach the sink due to the limited edge capacity, they need to be absorbed back by the source, which may take a long time and incur high computational cost. Since this part of excess can neither make impact on the max-flow nor change the final cut, such computation is redundant. ELCO reduces part the redundancy within each component through linear component operations, but it cannot solve the problem across components, especially with large (large capacities of Ecross). Under large K, Ewithin are easier to saturated and components are more likely to be segmented into smaller size. Thus, the redundancy is hard to be reduced if we focus each component alone.
Consider an extreme case in which all nodes are divided into the source side Vy according the cut. All excess can then reach the sink and all push operations are not redundant. In the other extreme, all nodes belong to the source side Vg. Most excess cannot reach the sink but is still pushed, leading to redundant computation. Reversely, rather than pushing excess from the source, better strategy pushing deficit from the sink. Deficit means input flow less than output flow. Pushing deficit is just pushing excess in graph with source and sink switched. Since all nodes are divided into Vss. all operations for pushing deficit in the opposite direction are not redundant. In conclusion, pushing excess in the sink side VT or pushing deficit in the source side Vg will not result in invalid

---

[e] Initialization
[b) Push excess
[c] Residuel jgraph
(d) Push deficit
V,
Vi, V4
Co
G
G
Ve
Cannot more Vt.
deficit
Figure 4: Bidirectional pushing. (a) Initial cut Co segments graph into Vs. and VTa (b) Replace nodes in Vs. by source, pushing excess leads to the new cut C1 which segments the new sink side Vi; out. (c) The residual graph of (b), while Vs, cannot push more before replacing nodes. (d) Replace nodes in Vr, by sink, the optimal cut of the whole problem is found by pushing deficit.
computation. If prior knowledge is available to guide us find an initialization of warping paths, we only need to push excess on the sink side and deficit on source side. After that, little computation would be wasted. This idea motivates us to design the bidirectional-pushin strategy (Fig(4)
Initialization: Estimate initial cut Co for GTW graph. The corresponding source side and sink side are denoted as Vs, and VT, respectively.
Push excess: Replace all the nodes in Vss by source, solve the max-flow problem by pushing excess. The new sink side segmented by the new min-cut C₁ is denoted as Vr,
Push deficit: Replace all the nodes in VT, by sink, solve the max-flow problem by pushing deficit. The obtained min-cut C2 is the min-cut of the original GTW graph.
The strategy is guaranteed to achieve optimal solution by the following statements:
Lemma Assume the min-cut segments the nodes V into source side Vs and sink side Vr, then replacing the nodes in Vss or VT by source or sink does not impact the min-cut.
Lemma 7: Assume Vr represents the nodes of sink side segmented by the real min-cut of original GTW graph, then Vr, S Vr-
Theorem 1: The obtained min-cut in bidirectional strategy is the min-cut of the original graph.
In joint alignment problem, since the final tend be similar because larger brings larger similarity penalty, a good initialization is easy to be found.
3.3 Bidirectional pushing with linear component operations
We propose BILCO as the top-level algorithm integrating ELCO and the bidirectional-pusieee strategy, its source code can be found at component is usually large so that ELCO can work efficiently. While under large k, the cuts in different subgraphs are similar and the initial solution can be easily estimated. Therefore, the bidirectional-pusee -pushing strategy works well. By setting initialization each subgraph the optimal warping path of the averaged subgraph, which is also the solution of joint alignment with infinite K, BIICCO could take both advantages and has a good performance under any k. We use such initialization in our experiments and get good performance. Therefore, we recommend it as the default setting for BILCO, least for similar applications.
Under small k, the
3.4 Complexity
Bidirectional-pushing strategy can convert the original problem into two smaller sub-problems, which may help accelerate any max-flow methods. For example, in the best case, the initialization is accurate and divides the problem into two half-size sub-problems. Then the speed for one O(VV /[EI) method may be accelerated nearly times since both |V| and [E] are halved. Even with a bad initialization, the strategy can only increase the operation number without impacting the complexity bound.
ELCO is somehow like the integration of DTW (with complexity e((Vl)) and HIPR (with complexity D E (5). On one hand, Ecross makes ELCO more complex than DTW, thus it has lower
7

---

N 2600, =0.2
[b)
= 2500,
N 2500. T
10
10
10
10
10
20
30
40
50
10
20
30
50
secconos
(d)
T 50, e =.22
T=50,
0.2
10
10
BILCO
HIPR
10P
HPP
BI-HPR
10
6q((/N
S0/TUN)
(a)-(d) compare the rumning time and memory usage under different graph size. (e) compares the running time of max-flow methods under different K, including BI-HIPR and ELCO.
Figure 5:
bound (2(VV). On the other hand, using linear component operations to push excess, ELCO has fewer operation units and less complexity than HIPR. Even in the worst case where each component contains only one node, EL.CO is still equivalent to HIPR with the same upper bound.
Therefore, the worst complexity of BILCO O(IVV
3.5 Memory usage
BILCO can also improve memory efficiency by utilizing the structure of GTW graph. In general max-flow methods, the majority of memory is used to record the relationship between nodes and edges. With known graph structure, BILCO stores nodes and edges in coordinate order that the relationship can be derived from their positions. The expense for storing components is insignificant in real applications since the component number is usually neglectable compared with large V and EEE. Table I compares the memory usage of BILCO and peer methods based on their implementation.
Table 1: Memory usage (bytes) among BILCO, IBFS, HPF, BK, and HIPR.
BILCO
IBFS
HPF
BK
HIPR
8 8
E
48 + F]
156|V| 8
48 V| + 64 F]
64 112
Experiments
In this section, we compare the efficiency of BILCO with four peer methods IBFS, HPF, BK, and HIPR under synthetic simulation and real applications. All experiments were conducted in MATLAB with Intel(R) Xeon(R) Gold 6140@2.30Hz 128GB memory, Windows 10 64-bit, and Mierosoft VC++ compiler. No GPU is used. All methods implemented in C/C++ with MATLAB wrapper.
4.1 Synthetic data
Following a similar scheme in 122, we simulated an image signal propagation dataset with varying pixels (N) and frames (2). 4-connected neighbors are used and thus the GTW graph has around [ 2 nodes and EE 7NT edges. A bell-shaped signal propagated from the center of the image to the boundary. Gaussian noise was added so that the signal-to-novl mimics real scenario. To make hyperparameter comparable, normalized the synthetic data by dividing the standard deviation of the noise. We tested 20 instances for each combination of N, T, and
ratio 10dB, which
Fig (a-d) compares BILCO and peer methods under different graph sizes, where BILCO least 10 times faster than any peer method and costs only 1 1/10 memory. The impact of hyperparameter is demonstrated in Fig 5 (e), where HIPR combined with our iinrrectional-usshin strategy (BI-HIPR) and ELCO method are also compared. As shown, the running time of the majority of the methods

---

Table Efficiency comparison of different methods on real data.
Problem
BILCO
IBFS
HPF
BK
HIPR
Time
Time/(Time of BILCO)
Name
IVI
Memory
Memoryy(Memory of BILCO)
Calculate signal propagation in imaging data (first generated by last two from
20 and 21
586s
Out
Out
Out of
Glia Ca
1.1 10"
3.8 10"
25.1GB
Memory
Memory
Memory
Memory
Glia Ca?+
7.0 107
2.3 10
7s
x50.0
x109.4
x348.0
x385.6
1.5GB
x12.1
x21.1
x
x18.6
32s
x17.8
x19.5
x480.2
x76.3
5 rat astro ATP
7.9 107
2.7 10
1.6GB
x
x22.4
x
x19.8
26s
x20.2
x18.4
x665.7
x74.1
Exvivo Ca
6.5 107
2.2 10
1.4GB
x12.1
x21.1
x
x18.6
Extract depth information binocular stereo (data from [16]
60655
Out
Out
Out of
Out of
penduluml
2.9x 10
7.2
59.2GB
Memory
Memory
Memory
Memory
7316s
Out
Out of
Out of
Out of
artroom2
2.9x10
7.2x10
58.9GB
Memory
Memory
Memory
Memory
35s
10.0
x 113.6
x32.5
x72.6
penduluml*
7.7x107
1,910
1.9GB
x7.9
x16.7
x7.9
x14.1
33s
x9.8
x x12.1
X x
x58,9
artroom2*
1.7x10²
109
L.8GB
x8.3
x16.7
x8.3
x14.8
Signature identification (data from 14)
1,7 x 10
5.9 10*
204s
x15.2
X 19.5
x34.8
X 1556.3
Real signatures
6.3 106
107
0.14GB
x11.4
x20.7
x11.4
x21.4
1.7 x 10*
8.8 10*
564s
X20.9
X x16.6
x77.9
X 145.4
Forgeries
1.6 10
107
0.53GB
x8.3
x14.7
x8.3
x14.5
grows as becomes larger, because larger capacities of Ecross give more freedom to exchange flow and make the graph more complex. Among all methods, BILCO shows the best efficiency under any k. By comparing ELCO with BILCO, or HIPR with BI-HIPR, we find that strategy speeds up both methods, especially under large s. It's because the default initialization in Section 3.3 gets closer to the optimal solution when k gets larger, which effectively removes unnecessary computation. Fig(5c) also shows that the two strategies, ELCO and bidirectional pushing are complementary and synergistic. When is small, ELCO dominates the acceleration. When is large, the bidirectional pushing makes bigger impact. Interestingly, the two strategies synergistically improve the speed for the large k. For example, under the largest K, both BI-HIPR and ELCO are nearly times faster than HIPR, suggesting speedup of 16 folds for independent effects, while BILCO shows around 100-fold speedup compared to HIPR.
bidirectional-pushing
4.2
Real data
Here compare our BILCO with four peer methods in three distinct application categories: caleu- lating signal propagation (1), extracting depth information [I], and signature identification [13). Since all these max-flow methods would give the same results, here we only compare running time and memory usage, as shown in Table. where in the name represents the spatially downsampled data. BILCO is used as the baseline and show how many times others cost. The first application is similar to the one in Section 4.1 and we won't 't go into detail here. In the second application, the same row between two images, while the
the depth can be derived from the misalignment depth distributions in adjacent rows are assumed to be similar. To avoid unnecessary computation while allowing relatively large disparity, we set window size with 1/5 sequence length in this experiment. In the third application, 20 persons' signatures [1] are identified where totally 500 real signatures and 500 forgeries are compared with the reference signature. Smoothness imposed on the feature sequences of the same signature, and in the table we show the summation of running time

---

and maximum memory usage for two signature categories. The detailed preprocessing steps and experiment settings are given in the supplementary. As shown, BILCO always performs best for all these joint alignment applications. It shows around 10 to 50 times speed improvement and only costs nearly 1 110 memory compared with the best peer method. Especially for the application of estimating signal propagation, BILCO performs averagely 29 times (range from 18 to 50) faster than other methods. Considering these applications are diverse and the dependency structures are various from 2-connected neighbors for depth, 4-connected neighbors for propagation. to all-connected neighbors for signature, this forcefully demonstrates the superiority and broad applicability of BILCO.
Conclusion and discussion
Compared with DTW. joint alignment utilizes the dependency in the data and could obtain better performance, however, at the expense of speed and memory. In this paper, developed BILCO algorithm to minimize such an expense by exploiting the special properties of the problem. It can efficiently solve joint alignment max-flow problems and get exactly the global optimal solution without sacrificing accuracy. Besides theoretical analysis, we demonstrated the efficiency of BILCO through both synthetic experiments and real applications, where BILCO showed around 10 to 50 times speed improvement and only cost averagely 1/10 memory compared with the best one among peer methods IBFS, HPF, BK, and HIPR. The testing was conducted on wide range of scenarios, suggesting the broad applicability of BILCO. We expect our work could facilitate greatly the application of joint alignment and help obtain better performance in different fields.
We observed that the bidirectional pusshiin strategy can not only be regarded as a general push-relabel acceleration trick but also work synergistically with linear component operations for BILCO. Because of the bidirectional strategy, the graph is separated so that the less unnecessary excess or deficit needs to be pushed, the fewer saturated edges and cuts appear, and then the larger components will be. With larger components, the number of operations can be reduced greatly, which makes BILCO even faster x100x) than the multiplication of BI-HIPR (4x) and ELCO (4x) when compared HIPR as mentioned in Section [.
Our BILCO algorithm is specialized approach to specific nlx-flow problem. While inventing max-flow algorithms has a long history, all the existing methods fail to solve our problem efficiently, no matter whether they are based augmenting path 44 [2] 8) localized operations in push- relabel 9 31. the combination of the two 77. or integrating graph decomposition [8 17.. To the best of our knowledge,
are the first of integrating DP and the generic push-relabe method to solve max-flow problems. Although there is a general feeling in the field that initialization does not help much for max-flow algorithms, initialization and is particularly powerful when coupled with linear component operations. We hope efforts here will be not only useful for the application of joint alignment but also inspiring for the methodology development of specific network flow problems. We have general sense that while it hard to invent better generic algorithms for classical problems, there are lot of opportunities for more specialized problems and indeed great demands due to the super large scale of recent problems.
showed bidirectional pushing is good strategy leverage
Acknowledgments and Disclosure of F Funding
We thank the anonymous reviewers for their helpful suggestions. This work was supported by grants NIH R01MH1 U19NS123719, and NSF 1750931.
References
[1] RK Ahujia, Thomas L Magnanti, and James B Orlin. Network flows: Theory, algorithms and applications. New Jersey: Rentice-Hall, 1993.
[2] Ghazi Al-Naymat, Sanjay Chawla, and Javid Taheri. Sparsedtw: A novel approach speed up dynamic time warping. arXiv preprint Ar0i: 2012.
[3] Donald Berndt and James Clifford, Using dynamic time warping to find patterns in time series. In KDD workshop, volume 10, pages 370. Seattle, WA, USA:, 1994.
10

---

[4] Yuri Boykov and Vladimir Kolmogorov. An experimental comparison of min-cut/max-Die flow algorithms for energy minimization in vision. IEEE transactions on pattern analysis and machine intelligence. 26(9):1 2004.
[5] Boris V Cherkassky and Andrew V Goldberg. On implementing push-relabel method for the maximum flow problem. In International Conference on Integer Programming and Combinatorial Optimization, pages 157-171. -priner, 1995.
[6] Andrew Delong and Yuri Boykov. A scalable graph-cut algorithm for nd grids. In 2008 IEEE Conference on Computer Vision and Pattern Recognition, pages 1-8. IEEE, 2008.
[7] Andrew V Goldberg, Sagi Hed, Haim Kaplan, Pushmeet Kohli, Robert E Tarjan, and Renato F Werneck. Faster and more dynamic maximum flow by incremental breadth-firs, search. In Algorithms- 2015, pages 619-630. Springer, 2015.
[8] Andrew V Goldberg. Sagi Hed, Haim Kaplan, Robert E Tarjan, and Renato F Werneck. Maximum flows by incremental breadth-fine search. In European Symposium on Algorithms, pages 457-468. Springer. 2011.
[9] Andrew V Goldberg and Robert Tarjan. A new approach to the maximum- -low problem. Journal of the ACM (JACM), 35(4))921-940, 1988.
[10] Dorit S Hochbaum. The pscudoflow algorithm: A new algorithm for the maximum-flov problem. Operations research, 56(4))992-1009, 2008.
[11] Hiroshi Ishikawa and Davi Geiger. Occlusions. discontinuities, and epipolar lines stereo. In European conference on compuder vision, pages 232-248. Springer, 1998.
[12] Pushmeet Kohli and Philip Torr Dynamic graph cuts for efficient inference in markov random fields. IEEE transactions on pattern analysis and machine intelligence, 29(12):2079-2085 2007.
[13] Manabu Okawa. Template matching using time-series averaging and dtw with dependent warping for online signature verification. IEEE Access, 151010-810199 2019.
[14] Javier Ortega-Garcia, J Fierrez-Aguilar, D Simon, Gonzalez, Marcos Faundez-Zanuy, V Espinosa, A Satue, Hermaez, J-J Igarza, Vivaracho, et al. Mcyt baseline corpus: a bimodal biometric database. IEE Proceedings Vision, Image and Signal Processing, 150(6):395-401, 2003.
[15] Thomas Prätzlich, Jonathan Driedger, and Meinard Müller. Memory-restricted multiscale dynamic time warping. In 2016 IEFF International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 569-572 IEEE, 2016.
[16] Daniel Scharstein, Heiko Hirschmüller, York Kitajima, Greg Krathwohl, Nera Nešić, Xi Wang, and Porter Westling. High-resolution sterco datasets with subpixel-accurate ground truth. In German conference on pattern recognition, pages 31 42. Springer, 2014.
[17] Alexander Shekhovtsov and Václav Hlaváč. A distributed mincut/maxflow algorithm combining path augmentation and push-relabel. International journal of computer vision, 104(3):35-342,, 2013.
[18] Petter Strandmark and Fredrik Kahl. Parallel and distributed graph cuts by dual decomposition. In 2010 IEF. Computer Society Conference on Computer Vision and Pottern Recognition, pages 2085-2092. IERF, 2010.
[19] Gineke A Ten Holt, Marcel JT Reinders, and Emile A Hendriks. Multi-dimensional dynamic time warping for gesture recognition. In Thirteenth annual conference of the Advanced School for Computing and Imaging, volume 300, page 2007.
[20] Yinxue Wang, Guilai Shi, David Miller, Yizhi Wang, Congchao Wang, Gerard Broussard, Yue Wang, Lin Tian, and Guoqiang Yu. Automated functional analysis of astrocytes from chronic time-lapse calcium imaging data. Frontiers in neuroinformatics, 11:48, 2017.
[21] Yizhi Wang, Nicole V DelRasso, Trisha V Vaidyanathan, Michelle K Cahill, Michael F. Reitman, Silvia Pittoo,, Xuelong Mi, Guoqiang Yu, and Kira E. Poskanzer. Accurate quantification astrocyte and neuro- transmitter fluorescence dynamics for single-cell and population-level physiology. Nature neuroscience. 22(11)119661144 2019.
[22] Yizhi Wang. David J Miller, Kira Poskanzer, Yue Wang, Lin Tian, and Guoqiang Yu. Graphical time warping for joint alignment of multiple curves. In Proceedings of the 30th International Conference on Neural Information Processing Systems, pages 6655-3663, 2016.

---

[23] Chiung- Ting Wu, Yizhi Wang. Yinxue Wang, Timothy Ebbels, Ibrahim Karaman, Gonçalo Graça, Rui Pinto, David M Herrington, Yue Wang. and Guoqiang Yu. Targeted realigmment of lc-ms profiles by neighbor-wise canpound-speciffie graphical time warping with misalignment detection, Bioinformatics. 36(9):2862-23) 2871, 2020.
Checklist
1. For all authors...
(a) Do main claims made the abstract and introduction accurately reflect the paper's contribu- tions and scope? [Yes]
(b) Did you describe the limitations your work? [Yes] See Section
(c) Did you discuss any potential negative socictal impacts your work? [No] Not applicable
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] See Section 2
(b) Did you include complete proofs of all thcoretical results? [Yes] All given in supplementary duc the space limitation.
3. If you ran experiments..
(a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as URL)? [Yes) See supplemental material
(b) you specify the training details (e.g.. data splits, hyperparameters, how they were chosen)?
[Yes] See the default hyperparameter setting in supplementary
(c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [Yes] See Fig.5
(d) Did you include the total amount of compute and the type of resources used (c.g., type GPUs, internal cluster, or cloud provider)? [Yes] See Section 4
4. If you using existing assets (e.g., code, data, models) curating/relcasing new assets...
(a) your work uses existing assets, did you cite the creators? [Yes]
(b) Did you mention the license the assets? [No] All referenced methods are public. Part of the data are obtained from the coauthor, part of the data are public and cited.
(c) Did you include any new assets either in the supplemental material or as a URL? [Yes] Our code of new method is given through in supplemental material.
(d) Did you discuss whether and how consent was obtained from people whose data you're us- ing/curating? [No] Part of the data are obtained from the coauthor, part of the data are public and cited.
(e) Did you discuss whether the data you using/curating contains personally identifiable informa- tion or offensive content? [No] The data does contain personally identifiable information offensive content.
5. If you used crowd sourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable?
[N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
12


Reviews:
[{"review_id": "l54AnfsoI_g", "reviewer": "Reviewer_ycT1", "summary": "This paper investigates the joint alignment of multiple time series. Compared to the existing max-flow method on GTW, the proposed method makes use of two properties of the joint alignment max-flow problem, i.e., joint alignment is a generalization of pairwise alignment, and a coarse approximate solution to joint alignment can be readily estimated, for significantly reducing computational time and space.", "questions": "1. Is it reasonable to include the comparison of the alignment performance of the compared methods?\n2. Are there any trade-off between the alignment performance and the efficiency of the proposed method?\n", "limitations": "The authors addressed the potential negative societal impact of the work well.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n1. The observation of two important properties of the joint alignment max-flow problem.\n2. The corresponding designs to the two properties, i.e., Excess pushing with Linear Component Operations (ELCO) that integrates DP and the push-relabel approach, and a bidirectional-pushing strategy to utilize prior knowledge as initialization, are reasonable for improving efficiency.\n3. The complexity analysis provides good insights in comparison.\n4. The experiments on both synthetic and real datasets indicate the efficiency of the proposed method.\n\nWeakness:\n1. The experiments focus on comparing time and space efficiency of different methods. Alignment accuracy could be compared to indicate whether there is a trade-off between effectiveness and efficiency of the proposed method.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "VVLk4pQHDq8", "reviewer": "Reviewer_dysU", "summary": "This paper improves the computation speed and memory cost of GTW [Wang, et al, 2016] method, for the multiple time-series joint alignment problem. Two major enhancements are introduced. The first is excess pushing with linear component operations, which treats components (the number is much less than nodes) as the basic unit. The second is a bi-directional pushing (push excess and push deficit), where pairwise DTW initialization is used.\nComprehensive experiments on both synthetic and real world datasets are conducted to evaluate the speed and memory improvement of the proposed method.\n", "questions": "See above and I am willing to hear more opinions.\n", "limitations": "NA", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n1. The joint alignment problem is of great importance and interest to the community. The GTW method converts the alignment problem into a max-flow problem that sheds the light on solving it optimally via techniques in graph theory. This paper further pushes this direction to speed up the process and reduce the memory cost significantly. This makes the GTW approach feasible in real world applications.\n2. The paper has illustrative figures to demonstrate the approach.\n3. Time and space complexity analysis is provided.\n4. Several max-flow algorithms are included in the experiments part. The proposed method performs consistently much better than those baselines.\n\nWeakness:\n1. The paper lacks some introduction / background on graph techniques, making it slightly difficult to follow especially when defining ‘drain’, ‘discharge’ and ‘component’. With an illustrative example that demonstrates how these operations are done, would significantly improve the readability of the paper.\n2. The most critical problem lies in the definition of GTW. In GTW, the joint alignment problem is converted to a slightly different problem: constraint the similarity of pairwise warping functions / paths. The warping function could vary a lot, but the effect could be very similar. Let me use edit distance to illustrate (slightly different than DTW, but shares similar DP spirits, so I use this as an easy example). Given a reference string, there is one string that needs prepending a letter to match the reference (this is one warping). There is another string that needs appending a letter to match the reference (another warping). These two warping / editing functions are very different, but their integral is the same (both accumulate to only one insertion edit), so their effects are also similar as well. Therefore, some external metric, like the integral of the warping function, their difference should be considered as the penalty, instead of the function itself. GTW makes this definition for ease of computation via max-flow algorithms, but lacks the practical meaning as it is often hard to interpret the warping functions’ difference, and hard to apply to practical tasks like finding the centroids of time series. I think this formulation reduces the interests of GTW based approaches to the community.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "OiGhEFCpvnh", "reviewer": "Reviewer_JqmR", "summary": "This paper proposes an algorithm BILCO for solving graphical time warping, an alignment method for multiple time series data. The authors realized a fast and memory-saving algorithm by focusing on the special characteristics of DTW graphs. Experimental results, including results for real data, show that their two proposed methods, ELCO and bidirectional pushing, dramatically reduce computation and memory consumption.", "questions": "- Are there any other applications of the proposed method other than graphical time warping?\n\ntypo\n- l.170: source side $V_T$", "limitations": "This paper proposes an algorithm for solving time series data alignment, which does not involve negative social impact.", "rating": 7, "confidence": 1, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "- The authors propose an algorithm for joint alignment that is much faster (10x~) and memory-saving than other methods.\n- Experimentally, the proposed method is applied to three real datasets that require multiple alignment, and shows impressive speedup and memory-saving performance.\n\n\n- The proposed method is not for general multiple time series alignment problem, but is an algorithm for solving a narrower problem, graphical time warping, in which warping paths in given neighbor relationships should be close, as written in Eq.(1). The application of the proposed method to general multiple alignment problem and other problems is not discussed.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept