extfieldsofstudy:
Computer Science
Mathematics

text:
Sequential Change-point Detection for High-dimensional and non-Euclidean Data

In many applications, it is often of practical and scientific interest to detect anomaly events in a streaming sequence of high-dimensional or non-Euclidean observations. We study a non-parametric framework that utilizes nearest neighbor information among the observations to detect changes in an online setting. It can be applied to data in arbitrary dimension and non-Euclidean data as long as a similarity measure on the sample space can be defined. We consider new test statistics under this framework that can detect anomaly events more effectively than the existing test while keeping the false discovery rate controlled at a fixed level. Analytic formulas approximating the average run lengths of the new approaches are derived to make them fast applicable to modern datasets. Simulation studies are provided to support theoretical results. The proposed approach is illustrated with an analysis of the NYC taxi dataset.

Sequential Change-point Detection for High-dimensional and non-Euclidean Data Lynna Chu and Hao Chen
Abstract-In many applications, it is often of practical and scientific interest to detect anomaly events in a streaming sequence of high-dimensional or non-Euclidean observations. We study a non-parametric framework that utilizes nearest neighbor information among the observations to detect changes in an online setting. It can be applied to data in arbitrary dimension and non-Euclidean data as long as a similarity measure on the sample space can be defined. We consider new test statistics under this framework that can detect anomaly events more effectively than the existing test while keeping the false discovery rate controlled at a fixed level. Analytic formulas approximating the average run lengths of the new approaches are derived to make them fast applicable to modern datasets. Simulation studies are provided to support theoretical results. The proposed approach is illustrated with an analysis of the NYC taxi dataset.
Index Terms-Anomaly detection; Online change-point; Streaming data; Graph-based tests; Non-parametric.

I. INTRODUCTION
S EQUENTIAL change-point detection aims to detect abrupt anomalies, observations that deviate from regular behavior, in a streaming sequence of observations as quickly as possible, while controlling the number of false alarms. In many modern applications, the sequence of observations may consist of high-dimensional observations (where the dimensional of each observation is larger than the sample size) or non-Euclidean objects (for example, sequences of networks or images). Examples include fraud detection involving large amounts of customer transactions [1]; disease surveillance or medical monitoring using sequences of images or multiple diagnostic measures [2], [3]; and consumer-based data streams such as wearable and smart device monitoring, internet of things sensors, network security systems, cybersecurity, and other web applications [4]- [6]. In all these examples, the goal is to detect a change (if a change is present) as quickly as possible once it occurs, while also limiting the risk of false discovery.
The sequential change-point setting can be formulated as follows: let the observation at t be denoted as Y t , t = 1, 2, . . . , n, . . .. Here, t could be the time index or some other meaningful indices, Y t could be a vector, image, or network, and n is the index for the observation currently being observed. When there is no change-point, Y t 's are identically distributed The work was supported in part by the NSF awards DMS-1513653 and DMS-1848579. Lynna  from an unknown distribution, denoted as F 0 . If there is a change-point at τ , the observation after τ are from a different (unknown) distribution: where F 0 and F 1 are two different probability measures.
This formulation is very general. F 0 and F 1 are unknown and not specified. We do not impose any constraints on how the change happens here. For example, if the observation is a high-dimensional vector, the change may occur in a subset of (unknown) data streams and the subset may be of size one.
Our aim is to construct a stopping rule (denoted as T ) that will alert us to a change as quickly as possible, once it occurs, but keep the number of false discoveries at a fixed level. Moreover, in keeping with the spirit of the problem formulation, the stopping rule should be able to handle modern data types and should be relatively fast and easy to implement. To be precise, this means we would like to define T such that the detection delay, E τ (T − τ |T > τ ), is small, subject to a fixed average run length (formally E ∞ (T ) ≥ c, where c is a pre-specified large value) without making assumptions on the underlying sequence of observations. Here E τ denotes the expectation under the hypothesis that the true changepoint happens at τ and E ∞ denotes the expectation under the hypothesis of no change.

A. Related Works
When the data is univariate (or scalar), the sequential change-point detection has been studied extensively (see [7] and [8] for a review). For low-dimensional data, likelihood based methods have been explored which require knowledge or parameter estimates of the probability density functions (see for example: [9]- [12]). For high-dimensional data, the available methods are somewhat limited and may impose strict assumptions on the data. As an example, many existing methods have the assumption that the different data streams are independent [13]- [17]. Other works allow for more flexible application; for example kernel-based methods [18] and a modified sliding window algorithm [19]. Computationally efficient methods have also been proposed which combine summary statistics based on geometric entropy minimization (GEM) with the cumulative sum (CUSUM) algorithm [20], [21]. For network data, [22] proposed a sequential approach that acts by embedding each graph into a vector domain, where a conventional multivariate change-point detection procedure can be then applied. In general, for non-parametric methods applicable to high-dimensional and non-Euclidean data, theoretical analysis establishing false discovery control is very difficult to carry out.
Recently, [23] proposed a new non-parametric framework that utilizes nearest neighbor information to detect changes in an online setting. They also provided a general, analytical formula for false discovery control. This method can be applied to data in arbitrary dimensions (with no assumption that the different data streams are independent) and to non-Euclidean data. The author proposed to use the following stopping rule: where n 0 , n 1 , and L are pre-specified values, N 0 is the number of historical observations with no change-point, n > N 0 , and Z L|y (t, n) is a two-sample test statistic that tests whether {Y n−L+1 , . . . , Y t } and {Y t+1 , . . . , Y n } are from the same distribution. We refer to Z L|y (t, n) as the edge-count twosample test based on k-nearest neighbor (k-NN) in the following for simplicity. For more details of this test, please see Section II. The author also provided an analytic formula to compute b Z such that the average run length E ∞ (T Z (b Z )) is controlled at a pre-determined value. Simulation studies show that this method beats likelihood-based methods when the dimension is high.
Despite these nice properties, we find that the edge-count two-sample test on k-NN can have low power for some common types of changes when dimension is moderate to high, causing the stopping rule (1) to behave unexpectedly. To illustrate, consider a simple scenario where data are from a d-dimensional Gaussian distribution and there is a change at τ = 201: with Σ(i, j) = 0.3 |i−j| . We consider two types of changes: • Scenario 1 (only mean differs): ||µ|| 2 = ∆. • Scenario 2 (both mean and variance differ): ||µ|| 2 = ∆ and σ. Tables I presents the performance of T Z (b Z ) for both scenarios based on 1,000 simulation runs. Here k = 5, L = 200, n 0 = 25, and n 1 = 175. The table reports the fraction of trials (out of 1,000) to successfully detect the change within 30 (or 50) observations after the change occurs. The average detection delay (EDD) is estimated as the average elapsed time between when the change occurs (τ = 201) and when T Z (b Z ) detects a change. False alarms are not counted here. In each scenario, the threshold b Z is computed by formulas given in [23] such that E ∞ (T Z (b Z )) = 2000.
In Table I, we see that the method performs worse in Scenario 2 than Scenario 1 under two different comparison criteria: the fraction of trials that can be detected given a fixed time is smaller under Scenario 2 and the average detection delay (EDD) is larger under Scenario 2. However, common sense tells us that the additional change in variance should make the two distributions more different and the change in Scenario 2 easier to detect. This phenomenon is due to the curse-of-dimensionality (a more detailed explanation can be found in Section II-B) and results in the existing method having diminished power to detect general changes.

B. Our contributions
To address the problem of the stopping rule T Z (b Z ), we propose three new stopping rules: The definitions of S L|y (t, n), W L|y (t, n), and M L|y (t, n) are provided in Sections II-C, II-D and II-E, respectively. Under the same setup detailed in Table I, Table II shows that these new stopping rules are more successful in detecting the change quickly after it has occurred. They also have shorter detection delays than T Z (b Z ) under both above scenarios and all have shorter detection delays in Scenario 2 than that in Scenario 1. Further comparisons between the stopping rules can be found in Section IV; these demonstrate that the new stopping rule have improved power and detection delay compared to existing methods over a range of general scenarios. To construct these new stopping rules, we propose new twosample tests on k-NN. Specifically, we extend the generalized (S) /weighted (W ) /max-type (M ) edge-count test defined on an undirected similarity graph [24], [25] to the directed k-NN graph. The generalized edge-count and max-type tests on k-NN are well defined except for a particular construction of a k-NN graph (see Theorem 1). The detailed definitions of these stopping rules are given in Section II.
To make the new stopping rules useful for real-data applications, we provide analytical formulas to compute the stopping thresholds, b S , b W , and b M , such that the average run length (ARL) for each new stopping rule is controlled at a pre-determined value. This involves studying how the k-NN graph updates and obtaining expressions that explicitly characterize these dynamics, which lead to the development of new theoretical treatments. Specifically, for all the stopping rules, more accurate expressions of the graph updates are derived and the techniques used are an improvement over the approach utilized in [23]. We demonstrate that the analytical formulas for the stopping thresholds are reasonable to use and we further improve upon their accuracy for finite sample sizes by implementing a skewness correction technique on the thresholds.
In general, each test statistic has its own niche where it dominates. When interested in general change (for example, both mean and variance change), the stopping rules T S (b S ) and T M (b M ) are recommended. The stopping rule T M (b M ) has an advantage over T S (b S ) in that we can obtain more accurate analytical expression for the ARL for false discovery control. If the change of interest is in mean only, the stopping rule based on T W (b W ) is recommended. See Section IV for a comparison of their performance. In this paper, the types of changes we explore are confined to mean and/or variance change. However, the approach can be used to detect other changes in distribution, such as changes in covariance. For illustration, Table III shows the performance of the stopping rules under covariance change only. The data are again generated from d-dimensional Gaussian with d = 100, k = 5, L = 200, n 0 = 25, n 1 = 175, and there is change at τ = 201: with Σ 0 (i, j) = 0.3 |i−j| and Σ 1 (i, j) = 0.68 |i−j| . We see here again that the new stopping rules based on W , S, and M perform better than the stopping rule based on Z.
These new approaches are implemented in an R package gStream.
The organization of the rest of the paper is as follows. Section II discusses the new stopping rules in details. Section III studies the asymptotic properties of the proposed stopping rules and analytic ways to determine the thresholds. The performance of the new methods are further explored in Section IV and the new methods are illustrated on a real data set in Section V.

II. NEW TESTS
The test statistics in the stopping rules (2) -(4) trace from the offline version of the problem studied in [25] where the statistics were defined on an undirected similarity graph. For online detection, observations keep arriving and the similarity graph updates as a new observation arrives. Therefore one needs to understand the dynamics of the series of similarity graphs. [23] studied the directed nearest neighbor graphs as the dynamics of nearest neighbor graphs can be well understood. In this work, we continue to use the directed nearest neighbor graphs to construct improved statistics. The extension to other types of graphs is saved for future work.

A. Notation
First we define a random variable which indicates whether or not an observation Y i is among the rth nearest neighbor to another observation Y j among the observations in n L . Specifically, for any n > N 0 and i, j ∈ n L ∆ = {n − L + 1, . . . , n}, we let A (r) is the indicator function that takes value 1 if the event is true and 0 otherwise. In terms of graph construction, each observation points to its k nearest neighbors. For example, if A (r) to be the indicator function that Y j is one of the first k NNs of Y i among the observations in n L . We use y i 's to denote the realizations of Y i 's and let a + n L ,ij = k r=1 a (r) n L ,ij with a (r) n L ,ij = I(y j is the rth NN of y i among y n−L+1 , . . . , y n ). For any n, each t ∈ {n − L + 1, . . . , n} divides the data sequence into two groups: one group being the observations before t : {Y n−L+1 , . . . , Y t } (Group 1) and the other group being the observations after t : {Y t+1 , . . . , Y n } (Group 2). Define, b 0,ij (t, n L ) =I(n − L + 1 ≤ i ≤ t and t < j ≤ n or t < i ≤ n and n − L + 1 ≤ j ≤ t), b 1,ij (t, n L ) =I(n − L + 1 ≤ i ≤ t and n − L + 1 ≤ j ≤ t), b 2,ij (t, n L ) =I(t < i ≤ n and t < j ≤ n).
Then b 0,ij is the indicator function that Y i and Y j belong to different groups, b 1,ij is the indicator function that Y i and Y j both belong to Group 1, and b 2,ij is the indicator function that Y i and Y j both belong to Group 2.
We define our test statistics as follows:  B 0,ij (t, n L ), B 1,ij (t, n L ), and B 2,ij (t, n L ) are the random variable versions of b 0,ij (t, n L ), b 1,ij (t, n L ), and b 2,ij (t, n L ) such that the distribution of these random variables is defined to be the permutation distribution. The permutation distribution is the distribution induced by all L! possible permutations of observation indices among the observations in n L . The null hypothesis and independence assumption imply that the observation indices are exchangeable, and therefore under the permutation distribution, the true (unknown) distribution of the observations remains unchanged when the null hypothesis is true.
It is clear that R 0,L (t, n) is twice the number of edges in the k-NN graph connecting observations before t and after t, R 1,L (t, n) is twice the number of edges connecting observations prior to t, and R 2,L (t, n) is twice the number of edges that connect observations after t. The notation of the graph-based test quantities emphasizes their dependency on the graph which is constructed on the L most recent observations, with the most recent observation indexed at n. Figure 1 illustrates these test statistics constructed for different times t and n. In the top row of Figure 1, n = 30 and we construct the graph on the L = 20 most recent observations: y 11 , y 12 , . . . , y 30 . The edge-counts are calculated for different values of t ∈ {11, . . . , 30}. R 0,L (t, n) is twice the number of directed black edges, R 1,L (t, n) is twice the number of directed red edges, and R 2,L (t, n) is twice the number of directed blue edges. It is clear for a fixed n, the graph does not change: each graph in the first row of Figure 1 is the same. But as new observations continually arrive, the graph itself will update. For example when n = 40, we now construct the graph on observations y 21 , y 22 , . . . , y 40 (see the second row of Figure 1). For n = 40, the edge-counts are then calculated for different values of t ∈ {21, . . . , 40}.

B. Limitations of the method based on the edge-count test (Z)
The stopping rule T Z (b Z ) (1) is based on the edge-count two-sample test statistic Z L|y (t, n) [23]. To obtain Z L|y (t, n), first a directed k-NN graph is constructed based on a similarity measure (for example, Euclidean distance). Then the number of edges in the k-NN graph that connect observations before t (Y n−L+1 , . . . , Y t ) and after t (Y t+1 , . . . , Y n ) is counted (we refer to this as the between-sample edge-count). Please see Figure 1 for an illustration of how R 0, . Analytical expressions of E(R 0,L (t, n)) and Var(R 0,L (t, n)) can be obtained under the permutation distribution, which is defined as all L! possible rearrangements of the observation indices among n L , and are omitted here for brevity. A relatively low between-sample edge-count (or a large Z L|y (t, n)) indicates the observations before and after t are less mixed and this is evidence against the null hypothesis of no change. The intuition is that observations before t tend to find their nearest neighbor among other observations before t (and similarly for observations after t), which implies a distributional difference between the two groups of observations.
The rationale of a relatively small between sample edgecount (R 0,L (t, n))) holds particularly well under the scenario of mean change only and/or low-dimensional data. However, when the dimension of the data is moderate/high and the change in distribution is not only in mean, for example an additional variance change is present, this rationale breaks down. This is because in the high-dimensional setting observations that are similar (i.e. from the same distribution) are not necessarily close in distance. Consider two distributions that differ in both mean and variance and the dimension of the data is d = 100. The observations would be separated into two layers: observations with the smaller variance in the inner layer and observations with the larger variance in the outer layer. Since the volume of a d-dimensional space increases exponentially in d, for practical sample sizes it is not uncommon for observations in the outer layer to find themselves to be closer to observations in the inner layer compared to other observations in the outer layer. In this scenario the between-sample edge-count would be relatively large under the alternative, rendering Z L|y (t, n) ineffective. The new two-sample test statistics S L|y (t, n) and M L|y (t, n) address this curse-of-dimensionality problem.
Moreover, even for mean change only the stopping rule based on Z L|y (t, n) can still suffer from increased detection delay. This limitation of Z L|y (t, n) is due to a variance boosting problem under unequal sample sizes in the twosample test setting. For more details on the variance boosting issue in the two-sample test setting see [26]. In the sequential setting, this means that using Z L|y (t, n) can lead to increased detection delay since it may only be able to detect the change when it is near the middle of the sequence. To resolve this, we propose the weighted edge-count two-sample test W L|y (t, n).

C. The method based on the generalized edge-count test (S)
The generalized edge-count two-sample test at t ∈ {n−L+ 1, . . . , n} under k-NN can be defined as and Σ y (t, n) = Var((R 1,L (t, n), R 2,L (t, n)) T |y) such that E(·) and Var(·) denote the expectation and variance taken under the permutation distribution.
If a change-point τ > N 0 occurs in the sequence, we would expect S L|y (t, n) to be large when n > τ and t close to τ . The test statistic is defined in this way so that either direction of deviations in the number of within-sample edges from its null expectation would contribute to the test statistic. For example, under the location alternatives, we would expect both R 1,L (t, n) and R 2,L (t, n) to be larger than their null expectations, which would lead to a large S L|y (t, n). Under the scale alternatives, the group with the smaller variance would have a within-edge count larger than its null expectation and the group with a larger variance would have a within-edge count smaller than its null expectation (due to the curse-ofdimensionality), which would also lead to a large S L|y (t, n). Therefore, this test is powerful for both location and scale alternatives.
The following theorem ensures that the statistic is well defined except in very rare scenarios. These scenarios can be checked by calculating the node in-degree of each observation in the k-NN graph. A node's in-degree is the number of other observations that find that node to be its nearest neighbors. Theorem 1. For L ≥ 5, the generalized edge-count twosample test statistic under k-NN is well defined except for when all nodes have an in-degree of exactly k, i.e. d i = k ∀i, where d i = j∈n L a + ji . Based on Theorem 1, as long as each node in the graph does not have an in-degree of k and the graph is constructed on at least 5 observations, then Σ y (t, n) is invertible and S L|y (t, n) is well-defined. A proof of the above theorem is in the supplement.
This leads us to the stopping rule based on the generalized edge-count test under k-NN (2): D. The method based on the weighted edge-count test (W ) Following the same notations in Section II-C, for each t ∈ {n − L + 1, . . . , n}, the weighted edge-count two-sample test statistic under k-NN can be defined as Since it is more difficult for the sample with a smaller sample size to form an edge within the same sample, R 1,L (t, n) and R 2,L (t, n) are weighted by the inverse of their corresponding sample sizes. The test statistic defined in this way resolves the variance boosting problem described in Section II-B.
Relatively large values of R w,L (t, n) are evidence against the null hypothesis of no change. Let Var(R w,L (t, n)|y) .
Under the permutation distribution, analytical formulas for E(R w,L (t, n)) and var(R w,L (t, n)|y) can be calculated based on E(R 1,L (t, n)), E(R 2,L (t, n)), and Σ y (t, n) provided in Section II-C: The stopping rule based on the weighted edge-count test under k-NN is (3): E. The method based on the max-type edge-count test (M ) We can define the max-type test statistic under k-NN based on the following lemma: Lemma 1. The generalized edge-count two-sample test under k-NN can be expressed as The proof of this lemma is in the supplement. The analytical expressions for the expectation and variance of R diff,L (t, n) under the permutation null are: The variance of R diff,L (t, n)|y is well-defined as long as d i = k ∀ i, in other words as long as each node does not have an in-degree of k.
From the above lemma, S L|y (t, n) is the sum of squares of two uncorrelated quantities (these two quantities are further asymptotically independent; details given in Section III). Here, W L|y (t, n) is sensitive to location changes: when the change is in mean, W L|y (t, n) tends to be large. On the other hand, D L|y (t, n) is more sensitive to scale changes: when the change is in variance, |D L|y (t, n)| tends to be large. The sign of D L|y (t, n) depends on whether the distribution after the change has a larger spread or not. This leads to the following max-type edge-count two-sample test statistic under k-NN: M L|y (t, n) = max(|D L|y (t, n)|, W L|y (t, n)).
When there is a change in location and/or scale, depending on the signal of interest, it is useful to consider an extended version of the max-type edge-count two-sample test: where ξ ≥ 0. Different choices of ξ lead to different focuses of the alternatives. For example, if we are more interested in locational changes, we could choose a large ξ. On the other hand, setting ξ to be small would favor detecting scale changes. When ξ = 1, the test reduces to the plain max-type edge-count test. For more detailed discussion on how to select ξ, see Supplement H in [25] (under the offline change-point detection setting, but similar arguments apply to the online setting). The stopping rule based on the max-type edge-count test under k-NN is as follows : where n 1 = n − n 1 and n 0 = n − n 0 . This reduces to (4) when ξ = 1.

III. AVERAGE RUN LENGTH
Given the new stopping rules presented in Section II, we would like to determine the thresholds b S , b W , and b M ξ in an analytic way such that the false discovery rate is controlled at a pre-specified value. A common way to measure the false discovery rate under the online change-point detection is the average run length, i.e., the expected time to stop when there is no change-point, which we denote as In the comparisons in Section I-B (Table II), the thresholds were chosen such that the average run lengths are 2,000 based on simulation runs. This is doable when the underlying distribution of the sequence is known. However, in many applications, the distribution of the sequence is unknown. Furthermore, since new observations keep arriving, resampling based methods, such as permutation and bootstrap, are not appropriate here and even if they were, directly resampling could be very time consuming. Therefore, to make the method fast applicable, we seek to derive analytical expressions for the average run lengths. Given the non-parametric nature of the proposed method, we would not be able to get exact analytic formulas for the average run lengths under finite L. In the following, we first approach the problem asymptotically (Section III-A), and then make adjustments for finite samples (Section III-B).

A. Asymptotic results
To derive analytical expressions for the average run length, we must study the asymptotic distribution of the stopping rules. Since these stopping rules are composed of the random fields {S L|y (t, n L )}, {W L|y (t, n)}, and {M ξ,L|y (t, n)}, we study their asymptotic properties. To obtain the limiting distribution of these random fields, we only need to focus on {D L|y (t, n)} and {W L|y (t, n)}. We show that the limiting distribution of the random fields converge to independent twodimensional Gaussian random fields (Theorem 2). This proof utilizes Stein's method [27]. To fully specify the Gaussian processes, we must derive the covariance functions of the new processes (Theorem 3). Since the k-NN graph updates each time new observations arrive, we must study the dynamics of the k-NN series for the new test statistics. While [23] laid a framework for graph-based sequential detection, the techniques developed in [23] were not directly applicable to the stopping rules proposed in this manuscript and needed to be adapted. Integration of techniques from [25] were useful in this development. However, since the new stopping rules can cover more types of change, we found that a direct extension of these previous works were not sufficient in developing accurate analytical formulas approximating ARL. To push forward the theory, we carefully studied the dynamics of directed k-NN graphs and considered additional ways the graph can update. A more in-depth understanding of the graph updates led us to derive analytical expressions for these updates and incorporate them into our analytical formulas for the thresholds. A comparison of the improvement over a direct extension of previous works can be found in the supplement. Putting together Theorems 2 and 3 with results from [28], allow us to obtain analytical expressions for the average run lengths of the new stopping rules.
In the following, we provide a sketch of the key steps to obtain each result. We defer readers to the supplement for a more technical treatment. The subsequent results are derived under the following condition: Condition 1. There is a positive constant C, 1 ≤ C < ∞, depending only on k, such that In k-NN, each observation points to its first k NNs, so the out-degree of each observation is k. However the in-degree of each observation can vary. This condition states that the indegree of each observation is bounded. It is satisfied almost surely for multivariate data [29], [30]. Let Var(R diff,L (t, n)) , Var(R w,L (t, n)) , which replaces the conditional variances of the test statistics with the unconditional variances. See the supplement for more details.
Our proof that the limiting distributions converge to independent two-dimensional random fields depends on utilizing Stein's method [27]. The general idea of Stein's method is to show for a random variable W and a standard normal random variable Z, that for some family of functions h ∈ Lip(1), the following bound holds: where δ depends on the structure of W . The specific form of Stein's method used here requires A + n L ,ij to be locally dependent. However, even for different i, j, l, r, A + n L ,ij and A + n L ,lr are dependent due to the constraint that j∈n L A + n L ,ij = k for all i ∈ n L . Following [23], we relax these dependencies by considering a similar set of Bernoulli random variables {Ã + n L ,ij } i,j∈n L . We keep the following probabilities unchanged: P(Ã + n L ,ij = 1) = P(A + n L ,ij = 1), P(Ã + n L ,ij = 1,Ã + n L ,ji = 1) = P(A + n L ,ij = 1, A + n L ,ji = 1), P(Ã + n L ,ji = 1,Ã + n L ,li = 1) = P(A + n L ,ji = 1, A + n L ,li = 1), but relax the other dependencies such thatÃ + n L ,ij is independent of {Ã + n L ,il ,Ã + n L ,li } l =j , andÃ + n L ,ij andÃ + n L ,lr are independent when i, j, l, r are all different. ThenÃ + n L ,ij are only locally dependent and can be analyzed through the Stein's method [27].
We are now ready to present the main results.
converge to independent two-dimensional Gaussian random fields, which we denote as Here [x] denotes the largest integer smaller than or equal to x for any real number x.
The detailed proof for Theorem 2 is in the supplement. Based on Theorem 2, we can approximate E ∞ (T S (b S )), E ∞ (T W (b W )), and E ∞ (T M ξ (b M ξ )) by examining the asymptotic behavior of our stopping rules: where n 1 = n − n 1 and n 0 = n − n 0 . Our approximations involve the function ν(x) defined as This function is closely related to the Laplace transform of the overshoot over the boundary of a random walk. A simple approximation given in [31] is sufficient for numerical purposes: where Φ(·) is the cumulative distribution function of the standard normal distribution and φ(·) the density function of the standard normal distribution. According to [28], We now have approximations of the average run lengths for the new stopping rules. The only remaining unspecified quantities are the directional partial derivatives of the covariance functions of the Gaussian random fields. Their analytical expressions are derived in the following theorem.
Here, p k,∞ is the limiting expected number of mutual NNs a node has in k-NN, q k,∞ is the limiting expected number of nodes that share a NN with another node in k-NN, p ∞ (r, s) is the limiting expected number of mutual NNs shared between the rth and sth NNs, and similarly q ∞ (r, s) is the limiting expected number of nodes shared between the rth and sth NNs. Explicitly, p k,∞ = k r=1 k s=1 p ∞ (r, s), and q k, To derive these partial derivatives, we studied the dynamics of the k-NN series as new observations are added. It turns out that a few key quantities are enough to characterize the dynamics in the asymptotic domain. The proof of this theorem is in the supplement.

B. Finite L
We now consider the practical scenario where L is finite. Based on results in Section III-A, E ∞ (T S (b S )), E ∞ (T W (b W )), and E ∞ (T M ξ (b M )) can be approximated by with H L () and G L () are finite sample versions of H() and G(), respectively. In practice, when L is finite we use g W,2 (L, x) and g D,2 (L, x) in place of g W,2 (x) and g D,2 (x) in the above formulas, respectively, where Here, p k,L , p k+1,∞ and can be estimated in a data-driven way.

C. Skewness correction
Analytical approximations provided in Section III-B become less precise for finite L when n 0 is relatively small. This is mainly because the convergence of W L (t, n) and D L (t, n) to normal is slow if (n − t)/L is close to 0 or 1. This problem becomes more severe when dimension is high. To improve upon the analytic approximations for finite sample sizes, we perform skewness correction. We adopt a skewness correction approach discussed in [32] that does the correction up to different extents based on the amount of skewness at each value of t. In particular, we provide better approximations to the marginal probabilities P(W (u − x, w) ∈ b + du) and P(D (u − x, w) ∈ b + du). Following the method based on cumulant-generating functions and change of measure (details refer to [32]), we can approximate the marginal probability by where θ b is chosen such thatψ(θ b ) = b. By a third Taylor approximation, we get where γ L (t, n) := E P (Z L (t, n) 3 ) and explicit expressions are derived using combinatorial analysis. Skewness corrected thresholds are only obtained for W L|y (t, n) and D L|y (t, n), but not S L|y (t, n). This is because for S L|y (t, n) the integrand can easily be non-finite and the approach depends heavily on extrapolation. Therefore, the stopping rule based on M L|y (t, n) is often recommended because it can detect general changes but we can obtain more accurate stopping thresholds.

D. Checking accuracy of analytic formulas for the average run lengths
Here, we check the accuracy of the analytic formulas for the average run lengths. For all three new tests, we have analytic formulas based on asymptotic results (19), (20) and (21), and for the tests based on the weighted/max-type edge-count tests, we have analytic formulas after skewness correction, provided in the supplement. We compare the empirical ARL obtained from these analytic formulas (with skewness correction with applicable) to those obtained from 1, 000 Monte Carlo simulations. The analytical thresholds are obtained so that the average run length is 2, 000. We generated data from three different settings: multivariate normal with d = 10 (denoted by C1), multivariate t 5 with d = 100 (denoted by C2), and multivariate log-normal with d = 1000 (denoted by C3).
Results for different choices of n 0 are shown in Tables IV -VIII. We set n 1 = L − n 0 and k = 5. The asymptotic analytic results are denoted by 'A1' and the skewness corrected approximations are denoted by 'A2'. We see in general that the empirical ARLs obtained from the asymptotic approximations ('A1') are not very close, illustrating the need for skewness correction here. After skewness correction, the empirical ARLs are much closer to 2, 000. It is clear that the accuracy of the skewness corrected approximations depends on n 0 : in general, when n 0 = 40, the skewness corrected approximations do well across most dimensions.
We also investigate how the analytical threshold approximations perform compare to the Monte Carlo thresholds in the presence of change. In Tables V, VII, and IX we report the power, defined as the fraction of trials able to detect the change within 50 observations after the change occurs, and the average detection delay (reported in parenthesis) obtained using analytical formulas ('A1'), analytical formulas with skewness correction ('A2') and Monte Carlo thresholds ('MC') . The thresholds are obtained so that the ARL is set to be 2, 000. The amount of signal used correspond to the power setup described for Tables X, XII, and XI in Section IV, respectively. It is clear that in the presence of change, the analytical formulas (with skewness correction) and Monte Carlo thresholds all lead to similar results in power and detection delay.

IV. POWER ASSESSMENT A. Multivariate data
To examine the performance of the three new test statistics, we compare them to the existing approach in [23] (max t Z L|y (t, n)) and two parametric likelihood-based approaches: Hotelling's T 2 test when there is change in mean and the generalized likelihood ratio test when there is variance change (both these two-sample tests are adapted to the scan statistic setting). The simulation setup is as follows: there are N 0 = 200 historical observations and a change occurs at t = 400 (200 new observations after the start of the test). The observations are independent and follow a d-dimensional distribution. When there is a change in mean, the observations are shifted from 0 by amount ∆ in Euclidean distance. When the covariance matrix changes, to make the change less significant, only the first d/5 of the diagonal elements change with a multiple of σ, and the rest are unchanged. The amount of change is chosen so that the tests have moderate power to be comparable. For fair comparison, we use Monte Carlo simulations to determine the threshold for each of the test so that their average run lengths are all 2,000. Power is reported as the fraction of trials for which the change-point is detected within 50 observations after the change occurred. In the following, we use 'HT' to refer to the scan statistic over the Hotelling's T 2 statistic and use 'GLR' to refer to the scan statistic over the generalized likelihood ratio statistic. For d > 100, in order for HT and GLR to be applicable in higher dimensions, we treat each data stream as if it is independent so that the covariance matrix's inverse and determinant are well-defined.
When there is both mean and variance change, Table X shows the results under the Gaussian setting. When d ≤ 500, the GLR dominates in power. However, when the dimension increases, GLR is no longer able to retain competitive power compared to S and M . On the other hand, this setting is not well-suited for W which is meant to capture mean change only and its performance is the worst here.
To consider other distributions, we also compared the tests for multivariate log-normal data and multivariate t 5 data. The results for the log-normal data are shown in Table XI. Here there is a change in mean parameter only and ∆ is chosen such that the location change dominates. In this setting, W 's performance dominates. We see that all the new tests outperform Z and the parametric tests for d > 10.
The result for the multivariate t 5 data are shown in Table  XII. When there is a change in both mean and variance, Z is unable to outperform the new test statistics. Among the new test statistics, their performance depends on whether the mean or variance signal dominates. When the mean change is stronger (for example, when d = 100), W performs comparably well. However, when the variance change is stronger (for example, when d ≥ 1000), M and S dominate.
Based on the results of these tables, we see that the new graph-based methods perform well under various scenarios and have improved detection delay over the existing method in [23]. In general, if one is certain that the change is locational, the test based on W is recommended; while for more general changes, the tests based on S and M are recommended.

B. Network data
Non-Euclidean data, also referred to as object data, is simply data that does not lie in the Euclidean space. Examples include networks, images, shapes, and trees. Fundamental statistical tools involving vector space analysis are no longer applicable for object data, making it challenging to directly analyze these data types. To demonstrate the new test statistics' power on non-Euclidean data, we evaluate the proposed methods on a sequence of networks. We generate random networks using the configuration model, which is a specific method that allows us to generate random networks with a specified degree sequence (please see [33], [34] for an overview). For every node with degree k i , we create k half-edges (referred to as 'stubs'). The network is created by iteratively selecting two stubs uniformly at random and connecting them to form an edge. This is done until no stubs remain and will result in a network with a predefined degree sequence. For the graph-based results, we show the results for graphs constructed using two different similarity measures. Specifically, for a network at time t, we encode the network using an adjacency matrix A t with 1 for element (i, j) if node i and j are connected, and 0 otherwise. The similarity measures are: 1) Similarity 1: ||A t − A s || F , 2) Similarity 2: Under this setting, we compare the graph-based approach to another method designed specifically to detect changes in a stream of networks/graphs [22]. The approach proposed in [22] is quite general: it does not make explicit assumptions on the network/graph and can be applied to a stream of graphs of varying sizes. Their approach consists of two steps: (1) each graph G t is mapped to a vector y t through a prototype-based embedding, and (2) a change is then detected in the stream of vectors y t using any conventional multivariate changedetection procedure. Specifically, embedding is carried out by assessing the dissimilarity between graphs and a selection of prototypes. In the paper, they rely on the graph edit distance (GED), which count and weights the edit operations that are needed in order to make two input graphs equal and is applicable to graphs where the nodes are unidentified. To carry out their approach, we follow the implementation proposed in [22]; specific details can be found in Section IV of [22]. Since computing the GED is quite computationally expensive, we allow each network to have only 6 nodes. We set a change at τ = 101. Before and after the change, the total node degree remains the same, but the sequence of node degree changes. Before the change, each network is constructed such that nodes 1 and 4 have node degree of 1 and the remaining nodes have node degree 2. After the change, nodes 3 and 6 have node degree of 1 and the remaining nodes have node degree 2. Observe that in this setting Similarity 1 and 2 are equivalent and so we only report results for Similarity 1. The first 50 observations are treated as training observations for [22]. We set L = 50 for the graph-based approach. For fair comparison with the method in [22], we set the window size to be 50 and we implement their method using both GED and the Forbenius norm (Similarity 1) as the similarity measure between networks. The ARL is set to be 2,000. Power is defined as the number of trials (out of 100) where the change is detected within 50 observations. Among those trials where the change is detected, we also report the expected detection delay (EDD) and its standard deviation when applicable. We note that the approach in [22] depends on the choice of embedding and its hyper-parameters, which should be chosen carefully. A thorough discussion is provided in [22]. We carry out simulations under a variety of settings for the number of dimensions to embed (d) and the number of prototypes to select (nproto). Table XIII show the power performance of the graph-based methods compared to the approach in [22]. We can see that all of the graph-based statistics do better in terms of power and detection delay compared to [22]. For further exploration of the graph-based approach, we generate a sequence of networks such that each network has 20 nodes with a pre-specified degree sequence (see below for details). A change in the sequence of networks happens at τ = 101. The first 50 networks in the sequence are treated as historical observations. Power is reported as the fraction of trials for which the change-point is detected within 50 observations after the change occurred. Among those trials where the change is detected within the first 50 observations, the expected detection delay and its standard deviation is also reported (EDD).
The networks are generated under two different settings: 1) A fixed degree change in the network: before the change, half of the nodes have out-degree and in-degree 1 and half the nodes have out-degree and in-degree 3. After the change, 5 nodes have out-degree and in-degree k 1 and 5 nodes have out-degree and in-degree k 2 . The remaining nodes remain unchanged. 2) A random degree change in the network: before the change, half of the nodes have out-degree and in-degree 1 and half the nodes have out-degree and in-degree 3. After the change, half of the nodes have out-degree and indegree randomly selected from k 3 to k 4 and the remaining half of the nodes have out-degree and in-degree 3.
Tables XIV and XV report the power and expected detection delay of the graph-based test statistics for similarity measure 1 and 2, respectively. For similarity measure 1, the performance of S and M dominate in almost all settings, with the exception of when k 1 = 2, k 2 = 2 ; in general S and M do well with respect to both power and expected detection delay. For similarity measure 2, the performance of W and Z improve substantially, while the performance of S and M remain stable.

V. A REAL DATA APPLICATION
We compare the new approaches to the method in [23] using the yellow taxi trip records data. The data set is publicly available on the NYC Taxi & Limousine Commission (TLC) website. It provides information on the taxi pickup and dropoff date/times, longitude and latitude coordinates of pickup and drop-off locations, trip distances, fares, rate types, payments types, and driver-reported passenger counts.
Based on this data set, a natural question to ask is: Can we detect a change in traffic patterns during peak travel seasons? Here, we focus on those trips that began at John F. Kennedy International Airport and we look at two different time periods: the months of June through August and November through December in 2015. The dataset has been completely collected at the time of analysis. However, we treat it as if the data were being observed in order to illustrate how the proposed method works. For simplicity, the boundary of JFK airport was set to be 40.63 to 40.66 latitude and −73.80 to −73.77 longitude. For those trips that began with a pickup at JFK, we extract information on their longitude and latitude drop-off coordinates. Using longitude/latitude coordinates, we create a 30 by 30 grid of New York City and count the number of taxi drop-offs that fall within each cell, where each cell represents a longitude, latitude coordinate range. Then for each day, we have a 30 by 30 matrix such that each element represents the number of taxi drop-offs in each location. We apply the new approaches, as well as Z, to detect changes in the months of June through August 2015. We use data from the month of May as historical data. Applying the offline change-point detection method in [32] and [25] on the observations in May, we find there is no change-point in the first 30 days, so we set L = 30, n 0 = 5, and n 1 = L − n 0 . We denote A i to be the 30 by 30 matrix on day i and v i to be the vector form of A i , which is now 900 by 1. The L 1 norm is used to construct the k-NN graph representing similarity between days. Here, the new test statistics (W , M , and S) all report a stopping time of July 3 and July 4 whereas Z is unable to detect any anomaly event (Table XVI). The change-point triggering these stopping times is estimated to be June 29. To perform a sanity check, we plot a heatmap of the L 1 distance matrix used to the construct the k-NN graph (see Figure 2, left panel). Based on the heatmap, we can see there is a clear signal happening around Day 30, which corresponds with the results from the new test statistics.
To detect changes in November and December 2015, we use data from the months of September and October 2015 as historical data. Applying the offline change-point detection method in [32] and [25] on the observations in September and October, we find there is no change-point in the first 50 days. Therefore, we treat the first 50 observations from Sept. 1 -Oct. 20 as historical observations and we begin the test at Oct. 21. We set L = 50, n 0 = 8, and n 1 = L − n 0 . The stopping times based on the new test statistics report back dates that seem to be quite reasonable (see Table XVII). We see that multiple stopping times are caused by the same anomaly event. When the signal is large enough, the new test statistics (W , M , and S) and Z perform similarly: all are able to detect a change in travel pattern close to Thanksgiving and the Christmas holidays. Again to check our results, we plot a heatmap of the L 1 distance matrix used to the construct the k-NN graph. We can see that there is a clear signal starting roughly around Day 30 and again around Day 60 and Day 67, which matches the results reported from the test statistics. In comparison with the heatmap from the months of June through August, the signal from the summer months is much weaker and in that case Z is unable to detect any anomaly event.

VI. CONCLUSION
We propose new graph-based test statistics under k-NN for detecting change-points sequentially as data are generated. We study the asymptotic properties of the stopping rules based on the new test statistics, and derived the analytic formulas to approximate the average run lengths of the new stopping rules. To accommodate finite samples, skewness corrected approximations were also derived for the weighted and maxtype edge-count statistic under k-NN. The skewness-corrected versions give much more accurate approximations to the average run lengths and can be used reliably in practice. The performance of the proposed test statistics are examined under various common scenarios. Simulation studies reveal that the new test statistics have shorter detection delays for a wider range of alternatives and exhibit power gains for scale change when compared to parametric tests and the test statistic proposed in [23]. Specifically, simulation results show that the weighted-edge count statistic (W ) is useful at quickly detecting mean changes. When a change in variance is also of interest, the generalized edge-count statistic (S) and maxtype edge-count statistic (M ) are more effective in detecting changes and obtain faster detection. Together with the fact that skewness corrected average run length approximations can be obtained for the max-type edge-count statistic, the stopping rule T M ξ is recommended for sequential detection of general changes.


text:
A Bi-layered Parallel Training Architecture for Large-scale Convolutional Neural Networks

Benefitting from large-scale training datasets and the complex training network, Convolutional Neural Networks (CNNs) are widely applied in various fields with high accuracy. However, the training process of CNNs is very time-consuming, where large amounts of training samples and iterative operations are required to obtain high-quality weight parameters. In this paper, we focus on the time-consuming training process of large-scale CNNs and propose a Bi-layered Parallel Training (BPT-CNN) architecture in distributed computing environments. BPT-CNN consists of two main components: (a) an outer-layer parallel training for multiple CNN subnetworks on separate data subsets, and (b) an inner-layer parallel training for each subnetwork. In the outer-layer parallelism, we address critical issues of distributed and parallel computing, including data communication, synchronization, and workload balance. A heterogeneous-aware Incremental Data Partitioning and Allocation (IDPA) strategy is proposed, where large-scale training datasets are partitioned and allocated to the computing nodes in batches according to their computing power. To minimize the synchronization waiting during the global weight update process, an Asynchronous Global Weight Update (AGWU) strategy is proposed. In the inner-layer parallelism, we further accelerate the training process for each CNN subnetwork on each computer, where computation steps of convolutional layer and the local weight training are parallelized based on task-parallelism. We introduce task decomposition and scheduling strategies with the objectives of thread-level load balancing and minimum waiting time for critical paths. Extensive experimental results indicate that the proposed BPT-CNN effectively improves the training performance of CNNs while maintaining the accuracy.

INTRODUCTION
I N recent years, Deep Learning (DL) techniques have achieved promising results in various domains [1,2]. Convolutional Neural Network (CNN) algorithm is an important branch of DL. Benefitting from large-scale training datasets and the complex training network, CNN achieves high accuracy and is widely applied in various domains, such as image classification [3], speech recognition [4], and text processing [5]. However, the training process of CNN is very time-consuming, in which large amounts of training samples and iterative operations are required to obtain highquality weight parameters [6,7]. It is critical to accelerate the training process and improve the performance of CNN. Cloud computing, high-performance computing cluster, and supercomputing provides strong computing power for various applications [8][9][10]. Therefore, it is a critical issue that how to design an effective parallel CNN training model based on distributed computing clusters and address the challenges of data communication, synchronization, and workload balancing, while maintaining high performance and high accuracy.
Numerous enhancements were proposed to accelerate the CNN and DL algorithms by improving the execution environments, mainly focusing in two directions [8,[11][12][13][14][15][16]: (1) using multi-core CPUs or GPUs platform to provide high-speed processing capacities, and (2) using multiple distributed computers to increase the available computing power. The CPU/GPU based methods [8,12,13] can perform more arithmetic operations and are suitable for the training of modestly sized DL models. However, a known limitation of the these methods is that the training acceleration is small when the scale of training datasets or DL models exceeds the GPU memory capacity. To process largescale datasets and DL models, some distributed architecture based solutions were proposed, such as DistBelief [14], Caffe [15], and Tensorflow [16]. Each of these approaches achieves significant progress in different perspectives, with implicit own complications in engineering and implementation. Considerable improvements might be obtained by combining the advantages in both aspects and make use of the computing power of multiple machines in a distributed cluster and the high-performance CPUs or GPUs on each machine.
There exists multiple challenges in this regard. Firstly, the entire CNN model contains multiple CNN subnetwork models that are trained in parallel on different machines, which requires synchronization and integration operations. It is required to minimize the synchronization waiting problem between subnetwork models, and guarantee the accuracy of the integrated model. Moreover, high-quality CNN models often require large-scale training dataset and a large number of iterations. Hence, an effective parallel mechanism should be carefully designed to minimize the data communication overhead between different iteration steps, tasks in different threads/CPUs/GPUs, and distributed computers. Furthermore, considering the heterogeneity of distributed computing clusters, computing nodes might be equipped with different CPU or GPU structures and have different training speed. How to partition the training dataset into these computers and how many parallel training tasks are started on each computer to maximize the computing power and workload balance of each computer. Finally, how to design a scalable parallel mechanism according to the characteristics of the CNN network, which can be easily deployed on elastic computing clusters and meet the application requirements in different areas.
In this paper, we aim to address the above challenges and fully utilize the parallel computing capacity of computing clusters and multi-core CPU to accelerate the training process of large-scale CNNs. We propose a Bi-layered Parallel Training-CNN (BPT-CNN) architecture in distributed computing environments. The outer-layer parallelism is deployed in a distributed cluster to train data subsets in parallel with minimal data communication and maximal workload balance. The inner-layer parallelism is performed on each computer using multi-threaded platforms. Experiments on large-scale datasets indicate the advantages of the proposed BPT-CNN in terms of performance, data communication, workload balance, and scalability. The contributions of this paper are summarized as follows.
• In the outer-layer parallelism, an Incremental Data Partitioning and Allocation (IDPA) strategy is proposed to maximize the workload balance and minimize data communication among computers, where large-scale training datasets are partitioned and allocated to computers in batches according to their computing power. • An Asynchronous Global Weight Updating (AGWU) strategy is proposed to integrate CNN subnetwork models from different computers and to address the synchronization waiting problem during the global weight update process. • In the inner-layer parallelism, two time-consuming training steps of the CNN model are parallelized on each computer based on task-parallelism, including convolutional layer and local weight training process. • To achieve thread-level load balancing and critical paths waiting time minimization, we introduce task decomposition and scheduling strategies for CNN training tasks with multi-threaded parallelism. The rest of the paper is organized as follows. Section 2 reviews the related work. Section 3 presents the BPT-CNN architecture and the outer-layer parallelization process. Section 4 describes the inner-layer parallel training of BPT-CNN. Experimental results and evaluations are discussed in Section 5. Finally, Section 6 concludes the paper with a discussion of future work and research directions.

RELATED WORK
Previous works have proposed various hardware designs for CNNs and other deep learning algorithms acceleration [8,17]. FPGAs have been widely explored as hardware accelerators for CNNs because of their reconfigurability and energy efficiency [18]. In [11], a parallel solution of CNNs was designed on many-core architecture, in which the model is parallelized on a new platform of Intel Xeon Phi Coprocessor with OpenMP. Caffe [15] provides multimedia scientists and practitioners with a clean and modifiable framework for state-of-the-art deep learning algorithms. Chung et al. proposed a parallel Deep Neural Networks (DNNs) training approach for big data on the IBM Blue Gene/Q system, in which issues of regarding programming model and data-dependent imbalances were addressed [18]. In [19], a massively parallel coprocessor was designed as a meta-operator for CNNs, which consists of parallel 2D convolution primitives and programmable units.
To efficiently handle large-scale CNN and big data, outstanding distributed architecture based solutions were implemented in [14,20]. Adam [21] is an efficient and scalable deep learning training system, optimizing and balancing workload computation and communication through entire system co-design. An energy-efficient reconfigurable accelerator was presented in [22] for deep CNN. To minimize the energy cost of data movement for any CNN shape, a processing row stationary dataflow was introduced to reconfigure the computation mapping of a given shape. In [14], Dean et al. introduced a distributed system (termed DistBelief) for training large neural networks on massive amounts of data. DistBelief uses two complementary types of parallelism: distributed parallel between multiple models and in each model, respectively. In addition, an asynchronous Stochastic Gradient Descent (SGD) procedure was employed to support a large number of model replicas.
Tensorflow [16] is a popular framework for large-scale machine learning on heterogeneous distributed systems. The computation model of Tensorflow is based on dataflow graphs with mutable state, where the graph nodes can be distributed and executed in parallel on different workers, multi-core CPUs, and general-purpose GPUs. In addition, Tensorflow uses a declarative programming paradigm, and developers can focus on the symbolic definition and computation logic instead of the implementation details. However, there are some shortcomings in Tensorflow: (a) TensorFlow attempts to occupy all available GPU memory in the initial phase, which makes the machines deploying the Tensorflow program infeasible to share with other applications, and (b) many high-level operations and interfaces in TensorFlow are nested and chaotically packaged, making it difficult to customize programming. Hence, we study the parallelism idea of the DistBelief and Tensorflow approaches, and implement a bi-layered parallel training architecture for large-scale CNNs by combining the advantages of both distributed computing and CPU/GPU parallel computing.
Comparing with existing efforts, the proposed BPT-CNN architecture in this paper fully utilizes the parallel capacity of both the distributed cluster and multi-core CPU of individual machines. Benefitting from the proposed IDPA and AGWU strategies, we effectively improve the training performance of the CNN model and address the problems of data communication, synchronization, and workload balancing of distributed cluster. Moreover, according to task decomposition and scheduling strategies, BPT-CNN achieves the optimization objectives of thread-level load balancing and waiting time minimization of critical paths.

Convolutional Neural Networks
CNN model is one of the most representative network structures of DL technologies and has become one of the hot topics in various fields of science. The common architecture of a CNN network includes two components: a feature extractor and a fully-connected classifier. In a convolutional layer, each batch of the input dataset is analyzed to obtain different abstract features. Given an input X with scale (D x , H x , W x ), where D x , H x , and W x refer to the depth, height, and width of X. Assuming that a filter with scale (D f , H f , W f ) is used in the convolutional layer to extract a feature map A, then a i,j denotes the value of the jth column of the i-th row of the current feature map, as calculated in Eq. (1): where D f , H f , and W f are the depth, height, and width of the current filter, and f () is an activation function, such as sigmoid, tanh, or relu function [23].
Pooling layer is utilized on each feature map to reduce the feature dimensions. There exist various pooling methods, i.e., max pooling and mean pooling. Fully-connected layer is a classification layer of CNN, where all output features of convolutional or pooling layers are connected to all hidden neurons with weight parameters. An example of a CNN architecture is illustrated in Fig. 1. Massive training datasets and iterative training process guarantee the high precision of CNN. However, they become the performance bottleneck when the training network structure is complex and the computing power is insufficient. The characteristics and the corresponding limitations of CNNs are summarized as follows.
• Massive weight parameters: There are a large number of weight parameters in the network. If there are different connection structures among layers, the network will become more complex, which requires more computation time cost. • Large-scale training datasets and massive iterative operations for weight updating: High accuracy of CNNs requires high-quality weights, which depends on largescale training datasets and iterative training.

Bi-layered Parallel Training Architecture
To accelerate the training process of CNNs, we propose a bilayered parallel training architecture for large-scale CNNs. We describe the distributed computing environment and training process of the BPT-CNN architecture.

BPT-CNN Architecture
BPT-CNN architecture is composed of two main components: (a) an outer-layer parallel training for multiple CNN subnetworks on separate data subsets, and (b) an inner-layer parallel training for each subnetwork. The proposed BPT-CNN architecture is illustrated in Fig. 2. (1) Outer-layer parallel training. A data-parallelism strategy is adopted in the outer-layer parallel training process, where a large-scale dataset is split into multiple subsets and allocated to different computing nodes to be trained in parallel. At the parameter server, the global weight parameters of the entire CNN network are updated depending on the local weights from each training branch. The updated global weight parameters are shared to each machine for the next iterative training.
(2) Inner-layer parallel training. The inner layer adopts a task-parallelism strategy to further accelerate the training process of each CNN subnetwork on each computer. Two time-consuming computation tasks are parallelized, including convolutional layer and the local weight training process. Computation tasks on these processes are decomposed depending upon their logical and data dependence, and are executed with multi-threaded parallelism.

Distributed Computing Cluster for BPT-CNN
We construct a distributed computing cluster for the proposed BPT-CNN architecture to efficiently handle massive training datasets and large-scale CNN models. The distributed cluster mainly consists of a main server, several computing nodes with mult-core CPU, and a parameter server, as shown in Fig. 3.
The main server is responsible for CNN training task management as well as data partition and allocation. It copies the CNN training network and allocates them to each computing node. Meanwhile, the training dataset is split into a series of subsets and allocated to the corresponding computing nodes. During the parallel training, the main server monitors the training time costs on computing nodes and migrates datasets for the optimization objective of synchronization delay minimization.
On each computing node, samples in the subset are calculated by the corresponding CNN subnetwork, while the network weight parameters are trained as a local weight set. The training process on each computer is executed in parallel. In addition, each computer is equipped with a multi-core CPU platform. In an inner-layer parallel training, the training process of each CNN subnetwork is further parallelized using multi-threaded parallelism.
The parameter server collects the trained local weight parameters from each computing node and updates the global weight parameters. Then, the updated global weight set is re-allocated to each computing node for the next epoch of training.

Outer-layer Parallel Training of BPT-CNN
In BPT-CNN's outer-layer parallel training architecture, we address critical issues of distributed and parallel computing, including data communication, synchronization, and workload balancing. Firstly, considering the heterogeneity of the distributed computing cluster, we propose an incremental data partitioning and allocation strategy to maximize cluster's workload balancing and minimize data communication overhead. In addition, we propose an asynchronous global weight updating strategy to further minimize the synchronous wait in the global weight update process.

Incremental Data Partitioning and Allocation Strategy
Considering the heterogeneity of computing nodes and their different training speed, to maximize the workload balance of the distributed cluster and minimize the synchronization in global weight update process, we propose an Incremental Data Partitioning and Allocation (IDPA) strategy based on heterogeneous sensing. As there are no dependencies between training samples, they can be partitioned and allocated in batches instead of done at once, according to the computing power of the computing nodes. Assume that there are N training samples, m computing nodes in a distributed cluster, and K training iterations are required for the CNN model. Let A (A < K) be the number of batches of data partitioning, that is, the entire training dataset is incrementally partitioned and allocated in A times, and each time N A new samples are processed. Initially, we take the first N A samples as the training dataset in the first batch. Before the training iteration, we use the constant characteristics of the computing nodes to represent their heterogeneity, i.e., the CPU/GPU frequency is measured. Let µ j be the CPU/GPU frequency of computing node C j , and the number of samples that will be partitioned and allocated to C j is calculated as: After receiving the training samples, each computing node begins the first iteration of training. At the same time, we monitor the execution time of each computing node to complete the iteration and evaluate its actual computing power. Being of the opinion that there might be more applications from different employers executing on the compute nodes, although we can predict the computing power based on the computing node's CPU/GPU frequency, it is more accurate to evaluate its actual computing power by actual execution time. Therefore, after the first training iteration, we can partition the training dataset according to the actual computing power of the machines. Let T j be the execution time of computing node C j to train n (1) j samples in the current iteration, then we can get the average execution time of C j for a sample as t j = Tj n (1) j . We collect the execution time of the computing nodes in the current iteration and predict the execution time required by all computing nodes in the next iteration. Note that N A new samples will be partitioned and allocated in the a-th batch. Namely, there is a total of N A × a samples on the computing nodes. The average execution time of the computing nodes in the a-th training iteration is calculated as: where t = 1 m m j=1 t j is the average execution time for training a sample by any compute node. To minimize the synchronization latency among computing nodes during the global weight update process, we expect all nodes to complete each iteration as close as possible. Assume that n j samples on C j after the a-th batch partitioning and allocation, we can calculate the value of n j as: Accordingly, we can obtain the number of samples that C j can accept in the a-th batch allocation according to its actual computing power, as calculated as: Repeat this process A times until the entire training dataset is partitioned and allocated to the heterogeneous computing cluster. By considering the heterogeneity of computing nodes, each computing node receives a corresponding number of training samples based on its actual computing power. The total number of samples allocated to each computing node C j is denoted as n j , and n j = A a=1 n (a) j . The detail steps of the IDPA strategy are described in Algorithm 3.1. 1: if a = 1 then 2: for Cj in the a-th iteration; 11: 13: return ns.
Benefitting from the IDPA strategy, the training dataset is well partitioned and allocated to the computing nodes, allowing them to complete each iteration in same duration, achieving minimal synchronization delay and maximum workload balancing. Moreover, no data migration is required among compute nodes during the training process, thereby no unnecessary data communication overhead is incurred.
Recall that K iterations are required for the CNN model, that is, each sample has an average of K times to train the weight parameter set of the CNN network model. After the A iterations in the data partitioning process, each computing node continues to execute the remaining iterations on the n j samples. Since these samples are incrementally allocated to the computing nodes, the actual training times of samples in instead of N ×A. Therefore, we should recalculate the remaining iterations ∆K of the training process, as defined below: The total number of training iterations of the CNN model is To simplify the expression, we denote K as K in the remaining context.

Global Weights Updating Strategies
There are massive connections with different weight parameters among all layers in a CNN network. We define these weight parameters as a weight set. We need to collect the training results on each computing node to update the global weight set for the entire CNN network. In this section, we propose two global weight updating strategies for the CNN network. We respectively define the local weight of each CNN subnetwork and the global weight of the entire CNN network as follows.
Definition 1: Local weight set. The weight parameters among all training layers of a CNN training network are denoted as a weight set. At each computing node, the weight set of a CNN subnetwork is defined as the local weight set of the corresponding subnetwork. The local weight set is trained based on the related data subset. In a distributed computing cluster, there is a local weight set on each computing node, which is updated after training a sample.
Definition 2: Global weight set. The weight set of the entire CNN network is defined as the global weight set. We provide a parameter server for calculating the global weight set by combining parts or all of the local weight sets. The global weight set is aggregated by each local weight set and shared to all computing nodes for the next epoch of training.
(1) Synchronous global weight updating strategy. We propose a Synchronous Global Weight Updating (SGWU) strategy for BPT-CNN, where each computing node trains all the samples of the current subset and updates the local weight set for an iteration. The local weight sets trained by all computing nodes in the current iteration are gathered at the parameter server, where a new version of the global weight set is generated. The workflow of the SGWU strategy is illustrated in Fig. 4. Considering that different local weight sets are trained by the corresponding subsets on different computers, having different contributions for the global weight set. We verify the accuracy of each CNN subnetwork after completing an epoch of local iteration training and use it as the contribution of the local weight set. After all computers finish an epoch of local iteration training, the latest local weight set trained on each computer is aggregated to the parameter server to update the global weight set as a new version W (1) . The global weight set W (i) for the (i)-th epoch of iteration training is defined in Eq. (7): where W (i−1) j and Q (i−1) j are the local weight set and the corresponding accuracy of the CNN subnetwork on computer C j , which is obtained in the (i − 1)-th epoch of local iteration training.
In a distributed computing cluster, especially for one equipped with heterogeneous computers, although we use the IDPA strategy to maximize cluster's workload balancing, the SGWU strategy inevitably faces the synchronization problem during the global weight update process. Due to the different available computing capabilities, computers need different time costs to execute each training iteration. Let t (i) j be the execution duration for the (i)-th training iteration on computer C j . The waiting time for synchronization of the entire computing cluster is defined in Eq. (8): where K is the number of iteration training and m is the number of computing nodes.
To address the synchronization problem of SGWU, we propose an Asynchronous Global Weight Updating (AGWU) strategy. In AGWU, once a computing node completes a training iteration on the local samples, the updated local weight set is submitted to the parameter server to immediately generate a new version of the global weight set, without waiting for other computing nodes. Compared with SGWU, AGWU can effectively solve the synchronization waiting problem without increasing the communication overhead. The workflow of the AGWU strategy is shown in Fig. 5. Considering the heterogeneity of computing nodes, according to the IDPA strategy, each computing node may contain different scales of training subset. In addition, due to the different training speeds, each computing node may also submit its local weight set to the parameter server at different time points, and get different versions of the global weight set. For example, for a computing node C j with n j samples, we assume that C j train the local weight set W (k) j based on based on the version W k of the global weight set in the current iteration. During the current training iteration of C j , the global weight set has been updated from W (k) to W (i) by other computing nodes. In this case, the low speed computers train the local weight set based on the old version of the global weight set, while the high speed computers based on the newer version. Denote (W (k) j − W (k) ) as the increment between the submitted local weight set W (k) j and its base version global weight set W (k) . Assume that there is another local weight set W (i−1) j on C j and it is trained based on the version W (i−1) of the global weight set, where k < i − 1. It is easy to know that W (k) j has less impact than W (k) j in the process of updating W (i) . Hence, we can conclude that the local weight sets using the old version of the global weight set have less impact on the global weight updating than those using the new version of the global weight set. Therefore, we adopt a time attenuation factor γ to measure the impact of each local weight set to the current global weight update process. Denote γ (k) j as the time attenuation factor of the local weight set W (k) j submitted from C j , as calculated in Eq. (9): where (i − 1) is the latest version of the global weight set, and k is the version of the global weight set that used to train W (k) j . Since there is no dependence among the training subsets on different computing nodes, the global weight update process does not require the training results from all computing nodes at the same time. Once a local weight set is submitted, the current global weight set W (i−1) is immediately updated to a new version W (i) , without waiting for other computing nodes. The i-th version of the global weight set is updated as: where ∆W k→i j is the update component from W (k) j , and Q (k) j is the accuracy of the CNN subnetwork on computer C j , which is evaluated by the output of the current local iteration training.
After obtaining the updated global weight set W (i) , the parameter server shares W (i) to W j for the next iteration training. Subsequently submitted local weight sets from other computing nodes will update the global weight set based on the latest version. The steps of the AGWU strategy of BPT-CNN is described in Algorithm 3.2.
In comparison with the SGWU strategy, in the AGWU strategy, each computing node independently participates in the global weight update process, so there is no synchronization waiting problem in AGWU. Furthermore, from the perspective of the entire training process, the update of the global weight set depends on the training results of all  In BPT-CNN, data communication only incurs between each computing node and the parameter server for the global weight updating and sharing. In AGWU, to reduce the synchronization cost and data communication overhead from the perspective of computing nodes, after receiving a version of the global weight set, each computing node begin a training iteration. It does not participate in the global weight update and receive a new version before completing the current iteration. In both of SGWU and AGWU strategies, the global weight set is updated for every epoch of iteration training. Therefore, both strategies produce the same data communication overhead. Denote K as the number of CNN iteration training, data communication C SGW U in SGWU and C SGW U in AGWU between the parameter server and all computing nodes is calculated in Eq. (12): where c w is a unit communication cost for transmitting the global weight set between the parameter server and a computing node. For each update of the global weight set, there exist 2 iterations of data communication: (1) submitting the local weight set from a computing node to the parameter server, and (2) sharing the updated global weight set from the latter to the former.

INNER-LAYER PARALLEL TRAINING OF BPT-CNN
In the inner-layer parallel training of BPT-CNN, we further parallelize the training process for each CNN subnetwork on each computing node. Two time-consuming training steps are parallelized based on task-parallelism, including convolutional layer and the weight training process. In addition, we propose task decomposition and scheduling solutions to realize thread-level load balancing and critical paths waiting time minimization.

Parallelization of Convolutional Layer
In the training process of a CNN network, convolutional layers take more than 85.18% of the total training duration, but only train 5.32 -6.63% of the weight parameters [? ]. Fortunately, the matrix-parallel-based method provides an effective way of performing convolutional operations in parallel. We introduce the parallel mechanism of the convolutional operations into the inner-layer parallel training architecture of BPT-CNN. We use the data partitioning method of the input matrix in CNN and extract all convolution areas from the input matrix. Then, by sharing the filter matrix, all convolution areas are convoluted in parallel with the shared filter matrix. Given an input matrix X with the shape of (D x × H x × W x ), where D x , H x , and W x are the depth, height, and width of X. Providing a filter parameter matrix F with the shape of (D f × H f × W f ), a feature map A is generated via convolutional multiplication on X and F . Based on the scales of X and F , the shape of A is calculated as: where D a , H a , and W a are the depth, height, and width of A, respectively. Based on the scales of X, F , and A, we calculate the times K C of convolutional operations in the current convolutional layer, which will be executed in parallel. K C is calculated in Eq. (13): where S is the stride of the convolutional operation and P is the number of the zero padding, which means appending P laps elements around X with the value of 0.
To execute these K C operations in parallel, we need to identify the convolution areas of the input matrix X for each task. A convolution area X[r begin : r end , c begin : c end ] of X includes the begin and end rows and columns. In each convolutional operation task, an element-by-element multiplication is executed on X[r begin : r end , c begin : c end ] and F to generate the corresponding element a i,j of A. For each element a i,j in A, location indexes of the convolution area in X is calculated in Eq. (14): After obtaining location indexes of each convolution area, we extract the contents of different convolution areas and perform the related convolutional operations in parallel, without waiting for the end of the previous convolutional operations. These parallel convolutional operations on different areas access the input and filter matrices repeatedly and simultaneously from the same memory without updating the contents. Without data dependence among these tasks, different tasks can access different convolution areas in X simultaneously. An example of the parallel convolutional operation of each CNN subnetwork in BPT-CNN is illustrated in Fig. 6 and the steps of this process are described in Algorithm 4.1.
As defined in Eq. (13), the maximum parallelism degree of a convolutional layer is equal to the number of elements of the output feature map, which is computed according to the scale of X and F . The total execution duration T Conv of a convolutional layer is calculated in Eq. (15):  where |A| = H A × W A is the number of elements in A and T i is the execution duration of the i-th operation task.

Parallelization of Local Weight Training Process
To distinguish the weight set of the entire CNN network and that of each CNN subnetwork, we respectively define the global weight set and local weight sets in Section 3.3.2.
In this section, training process of the local weight set of each CNN subnetwork is parallelized on each computer. After obtaining the outputs of a CNN subnetwork, the error (loss function) of each layer is evaluated from the output layer to the first convolutional layer using the Back Propagation (BP) method. The Stochastic Gradient Descent (SGD) process [24,25] is involved in updating the weight parameters among all layers of the current CNN subnetwork. In the output layer, the square error of all neurons is taken as the objective function of weight training, as defined in Eq. (16): where E x denotes the loss function of the input x, and y i and y i are the label and the output of the neuron a i in the output layer, respectively. The error δ i of a i is the inverse of the partial derivative of the error of the input of a i , as calculated in Eq. (17): where x ji is the input of the neuron a i that connected with a j , that is, x ji is the output of a j . w ji is the weight of the connection between neurons a j and a i . Let δ l be the set of errors of neurons in the l-th layer L l . Based on δ l , the error set δ l−1 of neurons in L l−1 is calculated in Eq. (18): where W l−1 is the weight set of L l−1 and net l−1 is the weighted input of L l−1 , as defined as: where a l−2 is the output matrix of L l−2 , consisting of each element a l−2 i,j . An example of the calculation process of loss function between layers L l and L l−1 is shown in Fig. 7. We parallelize the process for the loss function calculation, where the errors of neurons in the same layer are computed in parallel. In the convolutional layer, each neuron in the output layer (a feature map) is connected to a part of neurons in the input layer (an input matrix). In such a case, the error calculation of neurons in the previous layer L l−1 depends on the results of a part of neurons in the next layer L l . Hence, we parallelize this process depending on neurons in L l−1 . An example of the loss function calculation parallelization is shown in Fig. 8. After obtaining the error set of neurons in L l , we calculate the error of each neuron in L l−1 . Let δ l−1 i,j be the error component of neuron a j in L l for a i in L l−1 , as defined as: where H f and W f are the height and width of the filter parameter matrix F between L l−1 and L l . Based on the error set of neurons, the weight parameters of F are computed subsequently. The gradient of each weight w i,j is calculated in parallel, as defined in Eq. (21): The gradient of the bias weight w b is computed in Eq. (22): Based on the gradient values, each weight w i,j is updated in Eq. (23): where η is the learning rate of the CNN network.

Implementation of Inner-layer Parallel Training
We implement the inner-layer parallel training of BPT-CNN on computing nodes equipped with multi-core CPUs. Based on the parallel models proposed in the previous section, computing tasks of these training phases are decomposed into several subtasks. The workflow of task decomposition for a CNN subnetwork is illustrated in Fig. 9. (1) Task priority marking. According to the logical and data dependence of the decomposed subtasks, a task Directed Acyclic Graph (DAG) is created. With the thread-level load balancing and completion time minimization as the optimization goal, the priorities of tasks in the task DAG are marked. We set a maximum value for the entrance task of the task DAG graph. Then, the priorities of tasks in each level are set according to the tasks' level. Specifically, upstream tasks' priorities are higher than that of downstream tasks, while tasks at the same level have the same priority.
Based on the priorities of tasks, we allocate these tasks to threads on the multi-core CPU platform using the priority task scheduling algorithm [26]. Based on the task priorities, tasks of the entire CNN training network are allocated to different threads on the different CPU cores. An example of the task scheduling of the CNN training network with multi-threaded parallelism is illustrated in Fig. 10.  8: find thread T h k from T Hs with minimal workload; 9: call Assignment A(P T s) ← (ti, T h k ); 10: remove ti from P T s; 11: return A(P T s).

Experimental Settings
All of the experiments are conducted on a distributed computing cluster built with 30 high-performance computing nodes, and each of them is equipped with Intel Xeon Nehalem EX CPU and 48 GB main memory, respectively. Each Nehalem-EX processor features up to 8 cores inside a single chip supporting 16 threads and 24MB of cache. Comparison experiments are conducted to evaluate the proposed BPT-CNN by comparing with Tensorflow CNN [16], DisBelief [14], and DC-CNN [23] algorithms, in terms of accuracy and performance evaluation. Large-scale public image datasets from ImageNet [3] with 14,197,122 samples are used in the experiments.

Accuracy Evaluation
We and 0.639, respectively. Because of the parallel training and global weight updating, BPT-CNN narrows the impact of local overfitting and obtains more stable and robust global network weights. As the epoch of iteration training increases, both of accuracy and AUC of BPT-CNN steadily increases. AUC of BPT-CNN is greater than that of Tensorflow by 5.91%, on average, 9.56% higher than that of Dis-Belief, and 10.09% higher than that of DC-CNN. Therefore, compared with Tensorflow, DisBelief, and DC-CNN, BPT-CNN does not reduce the accuracy of CNNs. Moreover, benefitting from the global weight updating strategy, BPT-CNN achieves more robustness than compared algorithms.

Execution Time of Comparison Algorithms
The execution time of these algorithms is compared using 100 training iterations in various configurations: different data sizes and computing cluster scales. The comparison of the average execution time of each algorithm in each case is shown in Fig. 12.
(a) Different data sizes (b) Different cluster scales As can be seen in Fig. 12 (a) and (b), the proposed BPT-CNN algorithm achieves higher performance than the compared algorithms in most of the cases. Benefitting from the data-parallelism strategy, when the data size increases, the volume of each partitioned subset on each computer is slightly increased, leading to a slight increase in the average workload of each computer. For example, when the number of training samples increases from 100,000 to 700,000, the execution time of BPT-CNN rises from 62.77s to 307.35s, while that of Tensorflow increases from 54.38s to 454.23s, and that of DC-CNN sharply increases from 91.21s to 929.74s. In addition, taking advantage of the IDPA strategy, the proposed BPT-CNN algorithm owns scalability over the compared algorithms. When the scale of the computing cluster expended, the execution time of BPT-CNN and Tensorflow is significantly reduced. Experimental results indicate that BPT-CNN achieves high performance and scalability in distributed computing clusters.

Execution Time Comparison for Fixed Accuracy
Considering the different training architectures of various comparison algorithms, we discuss how these algorithm trade off performance and accuracy with resource consumption. We discuss the training iterations required for each algorithm to achieve different accuracy, and then measure the execution time each algorithm takes under different computing resources. The comparison results are shown in Table 1 and Fig. 13.  From Table 1, all algorithms use similar iterations to achieve an accuracy of 0.650. However, to achieve higher accuracy, BPT-CNN requires fewer iterations than Tensorflow, DisBelief, and DC-CNN. For example, BPT-CNN requires 42 iterations to achieve an accuracy of 0.750, while Tensorflow uses 64 and DisBelief uses 85, and DC-CNN requires up to 147. In addition, to achieve an accuracy of 0.750, we compare the actual execution times of each algorithm under different numbers of computing nodes and CPU cores, as shown in Fig. 13 (a) and (b). When the scale of the computing cluster and CPU cores expended, the execution time of BPT-CNN and Tensorflow is significantly reduced. In contrast, the execution time of DisBelief and DC-CNN algorithms is increased when the cluster scale reaches a certain amount (e.g., [25][26][27][28][29][30][31][32][33][34][35], which is caused by the more data communication among the increasing machines. Experimental results indicate that BPT-CNN achieves higher accuracy and performance than other algorithms using the same computing resource. Moreover, when the scale of computing nodes and CPU cores increases, the performance benefits of BPT-CNN is more noticeable.

Execution Time of BPT-CNN with Different Strategies
We further evaluate the performance of the proposed BPT-CNN under different global weight update and data partitioning strategies. To evaluate the effectiveness of the IDPA strategy, we perform the same work using the Uniform Data Partitioning and Allocation (UDPA) strategy, where the training dataset is uniformly partitioned into m partitions and allocated into the m computers. Comparison experiments are conducted in terms of data size, computing cluster scale, CNN network scale, and thread size. The average execution time of BPT-CNN with different strategies is presented in Fig. 14.
In Fig 14 (a), 7 different scales of CNN network are constructed in the experiments, as described in Table 2.
Here "layers(Conv)" and "filters(Conv)" denote the number of the convolutional layer and that of filters at each layer, respectively. "layers(FC)" and "neurons(FC)" denote the number of layers in the fully-connected layers and number of neurons in each layer, respectively.  By comparing strategies AGWU and SGWU, the execution time of BPT-CNN using AGWU is obviously lower than SGWU in most cases. In AGWU, because of the asynchronous update of the global weight set, each computing node uses the minimum time to wait for the global weight update and trains almost continuously. In addition, by comparing data partitioning strategies IDPA and UDPA, benefitting from the incrementa data partitioning, the workload of computing nodes stays well balanced, which further shortens the waiting time among different machines. As shown in each case in Fig. 14, the execution time of BPT-CNN with IDPA is significantly lower than that with UDPA strategy. Hence, BPT-CNN using AGWU+IDPA strategies exhibits the most efficient performance against other cases. Moreover, with the increase of data size or CNN network scale, the execution time of BPT-CNN using AGWU+IDPA maintains a slow rise. When the computing cluster scale and the number of threads on each machine increases, the benefits of AGWU+IDPA are more noticeable. Taking advantage of the IDPA and AGWU strategies, BPT-CNN achieves significant strength in terms of performance.

Data Communication and Workload Balancing
We evaluate the proposed BPT-CNN architecture in the view of data communication overhead and workload balancing by comparing with Tensorflow, DisBeilef, and DC-CNN algorithms. 600,000 training samples are used in the experiments, and the number of computing nodes increases from 5 to 35 in each case. Experiment results of data communication and workload balancing are shown in Fig. 15. (a) Impact of data size (b) Impact of cluster scale It is clear from Fig. 15 (a) and (b) that, in most cases, BPT-CNN owns significant workload balancing and lower data communication costs than other algorithms. Due to the use of the IDPA strategy in BPT-CNN, there is only communication overhead between the computing nodes for transmitting local/global weight parameters, and no training sample migration is required. Hence, as the number of computing nodes increases from 5 to 35, the communication overhead of BPT-CNN slowly increases from 2.35 MB to 11.44 MB. In contrast, due to dynamic resource scheduling, Tensorflow generates communication overhead from 2.73 MB for 5 computers to 45.23 MB between 35 computers. Moreover, to achieve workload balancing, DisBelief and DC-CNN use data migration operations during training, which results in heavy communication overhead between computers.
We compare the workload balance of each algorithm under different scales of the computing cluster, as shown in Fig. 15 (b). Our BPT-CNN architecture considers the heterogeneity of compute nodes and allocates corresponding workloads based on the actual computing power of each compute node. Hence, as the scale of the cluster increases, BPT-CNN achieves a stable workload balance, keeping between 0.89 and 0.80. In contrast, without heterogeneityaware data allocation, the workload of other comparison algorithms is not as balanced as BPT-CNN. The unbalanced workload further leads to long waiting time for synchronization and more execution time for the entire CNN network. Experimental results demonstrate that BPT-CNN significantly improves the workload balance of the distributed computing cluster with acceptable communication overhead.

CONCLUSIONS
This paper presented a bi-layered parallel training architecture to accelerate the training process of large-scale CNNs.
In the outer-layer parallel training, the performance of the entire CNN network is significantly improved based on data-parallelism optimization, where the issues of data communication, workload balance, and synchronization, are well addressed. In the inner-layer parallelism, the training process of each CNN subnetwork is further accelerated using task-parallelism optimization. Extensive experimental results on large-scale datasets indicate that the proposed BPT-CNN effectively improves the training performance of CNNs in distributed computing clusters with minimum data communication and synchronization waiting.
For future work, we will further concentrate on scalable CNN models and the parallelization of deep learning algorithms on high-performance computers. In addition, development of deep learning algorithms specific applications is also an interesting topic, such as scalable CNNs for images and LSTMs for time series.