Animesh Datta
email: <adatta@ecn.purdue.edu

Swarup Bhunia
email: bhunias@ecn.purdue.edu

Saibal Mukhopadhyay

Nilanjan Banerjee
email: nbanerje@ecn.purdue.edu

Kaushik Roy
email: kaushik>@ecn.purdue.edu

Statistical Modeling of Pipeline Delay and Design of Pipeline under Process Variation to Enhance Yield in sub-100nm Technologies *

Operating frequency of a pipelined circuit is determined by the delay of the slowest pipeline stage. However, under statistical delay variation in sub-100nm technology regime, the slowest stage is not readily identifiable and the estimation of the pipeline yield with respect to a target delay is a challenging problem. We have proposed analytical models to estimate yield for a pipelined design based on delay distributions of individual pipe stages. Using the proposed models, we have shown that change in logic depth and imbalance between the stage delays can improve the yield of a pipeline. A statistical methodology has been developed to optimally design a pipeline circuit for enhancing yield. Optimization results show that, proper imbalance among the stage delays in a pipeline improves design yield by 9% for the same area and performance (and area reduction by about 8.4% under a yield constraint) over a balanced design.

Introduction

Increasing inter-die and intra-die variations in the process parameters, such as channel length, width, threshold voltage etc., result in large variation in the delay of logic circuits [START_REF] Bowman | Impact of Die-to-Die and Within-Die Parameter Fluctuations on the Maximum Clock Frequency Distribution for Gigascale Integration[END_REF]. Consequently, estimating circuit performance and designing high-performance circuits with high yield (probability that the design will meet certain delay target) under parameter variations have emerged as serious design challenges in sub-100nm regime [START_REF] Bowman | Impact of Die-to-Die and Within-Die Parameter Fluctuations on the Maximum Clock Frequency Distribution for Gigascale Integration[END_REF][START_REF] Jacobs | Gate Sizing Using a Statistical Delay Model[END_REF][START_REF] Borkar | Parameter Variations and Impact on Circuits and Microarchitecture[END_REF]. Statistical analysis of delay and techniques to enhance yield in combinational circuits have been proposed [START_REF] Jacobs | Gate Sizing Using a Statistical Delay Model[END_REF][START_REF] Choi | Novel Sizing Algorithm for Yield Improvement under Process Variation in Nanometer Technology[END_REF]. In the high-performance design, the throughput is primarily improved by pipelining the data and control paths [START_REF] Hennessy | Computer Architecture: A Quantitative Approach[END_REF]. In a synchronous pipelined circuit, the throughput is limited by the slowest pipe segment (segment with maximum delay) [START_REF] Hennessy | Computer Architecture: A Quantitative Approach[END_REF]. Under parameter variations, as the delays of all the stages vary considerably, the slowest stage is not readily identifiable. The variation in the stage delays thus result in variation in the overall pipeline delay (which determines the clock frequency and throughput). Thus, a statistical delay model is necessary to predict the delay distribution of the pipeline. Traditionally, the pipeline clock frequency has been enhanced by: a) increasing the number of pipeline stages, which, in essence, reduces the logic depth and hence, the delay of each stage [START_REF] Hennessy | Computer Architecture: A Quantitative Approach[END_REF]; and b) balancing the delay of the pipe stages, so that the maximum of stage delays are optimized [START_REF] Hennessy | Computer Architecture: A Quantitative Approach[END_REF]. However, it has been shown that if intra-die parameter variation is considered, reducing the logic depth increases the variability (defined as the ratio of standard deviation and mean) [START_REF] Borkar | Parameter Variations and Impact on Circuits and Microarchitecture[END_REF]. The effect of balancing the stage delays on the overall delay under parameter variation also needs to be analyzed. Thus, traditional deterministic approaches for maximizing the pipeline throughput need to be reinvestigated to understand their effect on pipeline yield under parameter variations. In this paper, we have developed statistical delay models and proposed a statistical design flow for enhancing the yield of a pipelined design considering inter-die and intra-die variations. In particular, in this paper, we have developed:

Analytical models to estimate the mean and the standard deviation of the overall pipeline delay from the individual stage delay distributions. Analytical models to predict yield of a pipelined design. Analysis of the effect of logic depth and imbalance in the stage delays and correlation among the different stage delays on the yield of a pipelined design.

A statistical design methodology to minimize the area (or power) of a pipelined design under a target yield and performance constraint. Our analysis shows that, under large intra-die process variations reducing the logic depth decreases the yield, which is in accordance with [START_REF] Borkar | Parameter Variations and Impact on Circuits and Microarchitecture[END_REF]. However, if inter-die process variation is dominant, the traditional approach of reducing the logic depth, results in better yield. It has also been shown that, balancing delays of different pipeline stages does not necessarily maximize the yield. Under parameter variations, proper imbalance among the stage delays result in yield improvement. We have come up with a simple heuristic (based on the area vs. delay trend of a circuit) to systematically incorporate proper imbalance among pipeline stages, so that yield can be improved with minimum area penalty. The heuristic is used to improve yield of an example 4 stage pipelined design by 9% from a balanced design (with stages independently optimized for equal stage delay). Based on the above observation, we have designed a fast gate-level sizing algorithm for the complete pipeline to minimize total area under a yield constraint. This optimization step reduces the area by 8.4% (while ensuring a 80% target yield) compared to a balanced design.

Yield estimation of a pipeline design under parameter variation

The delay of a pipeline is determined by the slowest pipe stage (Fig. 1(a)) [START_REF] Hennessy | Computer Architecture: A Quantitative Approach[END_REF]. The variation in process parameters results in a significant variation in the delay of each pipe stage and hence, the overall pipeline delay [START_REF] Bowman | Impact of Die-to-Die and Within-Die Parameter Fluctuations on the Maximum Clock Frequency Distribution for Gigascale Integration[END_REF][START_REF] Borkar | Parameter Variations and Impact on Circuits and Microarchitecture[END_REF].

Variation in pipeline delay and yield

The delay of a pipeline stage (SD) consists of the clock-to-Q delay of the latch (T C-Q ), propagation delay through the combinational logic (T comb ) and the setup time (T setup ) [START_REF] Hennessy | Computer Architecture: A Quantitative Approach[END_REF]. Both the inter-die and the intra-die distribution in process parameters result in the variation of the delay of the pipe stages. Inter-die distribution shifts the delay of each stage in the same direction (i.e. either all increases or all decreases). Due to intra-die distribution, the delay of different stages may shift in different directions. The random component of the intra-die distribution (e.g. V th variation due to Random Dopant Fluctuation [START_REF] Mahmoodi | Estimation of Delay Variations Due to Random-dopant Fluctuations in Nano-Scaled CMOS circuits[END_REF]) makes the delays of the different stages in a pipeline completely independent of each other (uncorrelated stage delays). On the other hand, the systematic variation in the parameters (e.g. spatially correlated W, L, T ox variations) makes the delays of the different stages correlated [START_REF] Bowman | Impact of Die-to-Die and Within-Die Parameter Fluctuations on the Maximum Clock Frequency Distribution for Gigascale Integration[END_REF]. Moreover, the stages can also be electrically correlated if the delay of the combinational logic of one stage affects the delay of the subsequent stage [START_REF] Mahmoodi | Estimation of Delay Variations Due to Random-dopant Fluctuations in Nano-Scaled CMOS circuits[END_REF].

Let us consider a pipeline consisting of N stages (Fig. 1). If SD i denotes the delay of i-th stage, then the overall pipeline delay (T P ) is the maximum of N individual stage delays and is given by [START_REF] Hennessy | Computer Architecture: A Quantitative Approach[END_REF]:

i i i P i C Q Comb Setup i N i N T T M a x S D M a x T T (1)
The stage delays can be represented as correlated Gaussian random variables (SD i ~N( i , i )) with mean µ i and standard deviation i [START_REF] Bowman | Impact of Die-to-Die and Within-Die Parameter Fluctuations on the Maximum Clock Frequency Distribution for Gigascale Integration[END_REF][START_REF] Choi | Novel Sizing Algorithm for Yield Improvement under Process Variation in Nanometer Technology[END_REF]. If we neglect the spatial and the electrical correlation, the stage delays can be considered as independent random variables. The overall pipeline delay T P given in (1) will also be a random variable with the mean T and the standard deviation T .

The probability (P D ) that the pipelined design will meet a specific delay requirement (say T TARGET ) is given by:

1,..., Pr{ } Pr{max } D P TARGET i TARGET i N P Yield T T SD T ( 2
)
We define this probability as a measure of the yield of the pipelined design. In this section we will discuss methods to analytically estimate the overall delay distribution of the pipeline (i.e. T , T ) and the yield (i.e. P D ).

Estimation of distribution of pipeline delay

The mean ( T ) and the standard deviation ( T ) of the pipeline delay (T P ) depend on the mean ( i ) and the standard deviation ( i ) of the each stage. A lower bound on T is given by (using the Jensen's inequality [START_REF] Ross | Useful Bounds on the Expected Maximum of Correlated Normal Variables[END_REF][START_REF] Clark | The Greatest of a Finite Set of Random Variables[END_REF]):

1,..., 1,..., [ ] [ ] { [ ]} T P i i i N i N E T E max SD max E SD (3)
The above equation states that the mean of the overall pipeline delay will be larger than the maximum of the mean delay of all the stages. To obtain an exact estimate of T and T we approximate T P as [START_REF] Clark | The Greatest of a Finite Set of Random Variables[END_REF]:

1 2 1 1 2 1 1 2 2
where, N n-1,n represents the normal approximation to max (SD n- 1 , SD n ). The mean ( n-1,n ) and the standard deviation ( n-1,n ) of N n-1,n can be approximated as [START_REF] Clark | The Greatest of a Finite Set of Random Variables[END_REF]:

1 1 2 2 2 2 2 1 1 1 2 2 2 1 1 1 1 , 1, 1, 1 2 1, 2 1 ; 2 .

Mean of N m

Std Deviation m m [START_REF] Borkar | Parameter Variations and Impact on Circuits and Microarchitecture[END_REF] where, represents the Cumulative Distribution Function (CDF) and represents the Probability Distribution Function (PDF) ( ( ) = (2 ) -1/2 exp(-2 /2)) of a standard normal ( = 0 and = 1) Gaussian variable [START_REF] Clark | The Greatest of a Finite Set of Random Variables[END_REF]. The correlation coefficient between SD n-1 and SD n is given by n-1,n . The correlation coefficient between N n,n-1 and SD n-2 (say SD n-2 , N n-1,n )) is given by [START_REF] Clark | The Greatest of a Finite Set of Random Variables[END_REF]:

2 , 1 2 , 1 2 , 2 , () ( )
n nn n n n n n n n SD N ( 6
) SD n-2 , N n-1,n
) is used to estimate the mean and the standard deviation of max(SD n-2 , N n,n-1 ). The above process is repeated by taking two variables at a time and finally T and T are estimated.

Estimation of yield

Using (2) the probability of delay failure (P D ) or yield can be estimated as [START_REF] Ross | Useful Bounds on the Expected Maximum of Correlated Normal Variables[END_REF]: The exact solution of ( 7) is possible by assuming the stage delays (SD i ) to be independent Gaussian random variables, as shown below [START_REF] Ross | Useful Bounds on the Expected Maximum of Correlated Normal Variables[END_REF]:

1,..., 1,..., Pr{max } Pr{ ( )} D i TARGET i TARGET i N i N P S D T S D T (7)
1 1,..., Pr{ ( )} n TARGET i D i TARGET i i N i T P S D T (8)
If the variables are correlated such a simplification is not possible. To estimate P D considering correlated SD i s, we approximate the overall pipeline delay (T P ) as a Gaussian random variable (with T and T estimated using ( 5)). Using this assumption P D is given by [START_REF] Ross | Useful Bounds on the Expected Maximum of Correlated Normal Variables[END_REF]:

Pr{ } TARGET T D D TARGET T T P T T (9)

Model verification

We have estimated the delay distributions in different pipeline structures (number of pipe stage logic depth in a stage). The models have been first verified by Monte-Carlo simulation of inverter chain pipelines with the transmission gate Masterslave flip-flops in 70nm BPTM [9]. SPICE Monte-Carlo simulation is first used to determine the mean and the standard deviation of the delay of each stage. The simulated i and i values for each stage are then fed into the proposed model to determine the distribution of the pipeline delay. It can be observed that, the delay distribution predicted by the proposed model closely matches the SPICE simulation result of a 5x8 pipeline for: (a) only random intra-die variations (i.e. stage delays are considered as independent) (Fig. 2 -I). The major source of error in the proposed modeling method is the assumption that the maximum of two Gaussian variables is also a Gaussian one (i.e. N n-1,n is Gaussian) [START_REF] Ross | Useful Bounds on the Expected Maximum of Correlated Normal Variables[END_REF][START_REF] Clark | The Greatest of a Finite Set of Random Variables[END_REF]. This assumption is valid only for two uncorrelated variables [START_REF] Ross | Useful Bounds on the Expected Maximum of Correlated Normal Variables[END_REF][START_REF] Clark | The Greatest of a Finite Set of Random Variables[END_REF]. Hence, the errors in the estimation of the mean and the standard deviation are expected to increase with the increase in the number of pipeline stages (Fig. 3(a)) and the correlation between the stage delays (Fig. 3(b)) [START_REF] Ross | Useful Bounds on the Expected Maximum of Correlated Normal Variables[END_REF][START_REF] Clark | The Greatest of a Finite Set of Random Variables[END_REF]. It can be observed that the increase in error in the standard deviation is more significant in both cases. However, in all cases the error in the standard deviation and the mean is less than 3% and 0.2%, respectively. The modeling error also depends on the ordering of the variables SD i in (2). It has been shown that the error is the minimum if the variables are ordered in increasing (or decreasing) sequence of their means [START_REF] Ross | Useful Bounds on the Expected Maximum of Correlated Normal Variables[END_REF]. We have used this ordering in our estimation to minimize the modeling error.

Table-I Modeling and simulation results of delay distribution and yield for different pipeline configurations

Monte

Estimation of the design space

Using the proposed models, we can estimate the design space for the mean and standard deviation of different stages of a pipeline that will satisfy a yield constraint. Assuming the overall delay of the pipeline to be Gaussian and using (3), the upper bound of the mean of each stage delay ( i ) is given by: However, this bound does not provide any estimate of the design space of i of each stage. A relaxed upper bound on the µ i and i for the i th stage can be obtained (Fig. 4) by assuming that any stage j i meets the yield requirement with probability 1 (since maximum value of (x) =1) as shown below:

1 i i D T A R G E T P T ( 11
)
Equation ( 11) states that, if the mean and/or the standard deviation of any stage falls outside this bound, no pipelined design with that stage can ever meet the target delay and yield.

A more stringent bound is obtained by assuming uncorrelated and equal stage delays and is given by:

1/ 1 S S N N TARGET i D i i TARGET D i T P P T (12)
where, N S is the number of pipeline stages. Depending on N s , it gives the set of values for µ i and i that can meet the target yield (Fig. 4 shows 2 such equality bounds for N s = n1, n2 with n1 < n2). Note that, there is a minimum bound on the µ i and i which depends on the minimum allowable logic depth and process specification (Fig. 4). Moreover, the µ i and i of a combinational circuit are related parameters and the relation determines the realizable design space for the µ i and i . For example, if we model each stage delay to be a chain of N L inverters then a simple relation between µ i and i is given by: 2 2 min min min min and

i L i L i i N N (13)
where, µ min and min are the mean and the standard deviations of an minimum sized inverter, respectively. Similarly for maximum sized inverter having parameters µ max and max , there will be another bound on realizable µ i and i (Fig. 4, curve marked Realizable upper bound). Hence, for equal critical stage delays a realizable design region bounded by a set of curves given by equation ( 12) to (13) is obtained (Fig. 4).

Observations on statistical pipeline design

Using the analytical models presented in section 2, we have analyzed the effect of logic depth and the unbalancing of the stage delays on the pipeline yield.

Trade-off between number of pipeline stages and logic depth

The delay of a pipelined design depends on the logic depth of the pipe stages. A reduction in the logic depth, which increases the number of stages in a pipeline, improves the operating frequency [START_REF] Hennessy | Computer Architecture: A Quantitative Approach[END_REF]. In this section we analyze the effect of logic depth and number of stages on the variability of a pipeline delay. A design with a lower variability has a higher probability of meeting a target delay (i.e. better yield) [START_REF] Borkar | Parameter Variations and Impact on Circuits and Microarchitecture[END_REF].

If only random intra-die variation is considered, the delays of the different gates in a combinational logic behave as independent variables. Under this condition, increasing logic depth reduces the variability ( / ) due to cancellation effect (i.e. delays of some gates increase while others decrease resulting in lower overall change) (Fig. 5(a)) [START_REF] Borkar | Parameter Variations and Impact on Circuits and Microarchitecture[END_REF]. However, correlation in the delays of the different gates (due to inter-die variation and spatial correlation) minimizes the cancellation effect. Hence, the variability becomes a weaker function of the logic depth (Fig. 5(a)). On the other hand, increasing the number of elements in the max function reduces its variability (Fig. 5(b), evaluated using a inverter chain pipeline with constant logic depth in each stage [START_REF] Ross | Useful Bounds on the Expected Maximum of Correlated Normal Variables[END_REF]. It can be further, observed that as the stage delays becomes more and more correlated, the sensitivity of the variability to the number of stages reduces (Fig. 5(b)). In order to understand the effect of logic depth (N L ), number of stages (N S ) and relative strength of inter-die and intra-die variations, we have estimated the variability of an inverter chain pipeline (in BPTM 70nm technology node) of 120 levels with different configurations (i.e. with different N L and N S such that N L x N S =120). In each case the delay distribution proposed by the analytical models closely match the Monte-Carlo simulation results from SPICE. With only intra-die variation, the effect of logic depth prevails over the effect of the max function. Consequently, increasing the number of stages (i.e. reducing the logic depth of each stage) increase the variability (Fig. 5(c)). On the other hand, as the strength of inter-die variation increases (i.e. stage delays become more correlated) the / ratio of each stage becomes a weaker function of its logic depth. Hence, the impact of max function prevails and variability of the overall pipeline delay reduces with an increase in the number of stages (Fig. 5(c)).

Perfectly balanced vs. unbalanced pipeline design

Traditionally, the pipeline stages are designed for equal delay to maximize the throughput [START_REF] Hennessy | Computer Architecture: A Quantitative Approach[END_REF]. However, incorporating imbalance among the pipeline stage delays can have a positive impact on the yield under process variation. This is because of the fact that a balanced pipeline has more number of critical paths than an unbalanced design that adversely affects the yield [START_REF] Borkar | Parameter Variations and Impact on Circuits and Microarchitecture[END_REF]. Imbalance can be incorporated in a balanced pipeline by using transistor sizing and/or logic re-structuring. However, the impact of imbalance on the overall pipeline delay needs to be estimated.

We have performed experiments with a 3-stage pipeline structure to understand the effects of imbalance on the pipeline delay. For example, Fig. 6 shows a three-stage pipelined ALU-Decoder circuit. The combinational logic of each stage is first optimized for minimum area (using [START_REF] Choi | Novel Sizing Algorithm for Yield Improvement under Process Variation in Nanometer Technology[END_REF]) for a specific target pipeline delay (179ps) and yield (80%). First, we kept the target delay and the target yield constant for all the three stages (i.e. neglecting any correlation, the yield target for each stage was kept at (0.80) 1/3 = 0.9283 using ( 12)). In the next step, we have introduced imbalance among the three stages (by transistor sizing) in such a way that the total area remains constant. The overall delay distributions and the yield of the balanced and unbalanced pipeline are shown in Fig. 7. Resizing the transistors reduces the mean but increases variance in pipeline delay (Fig. 7(a)). To understand the reason behind this yield improvement, let us consider the area vs. delay curves for each stage (Fig. 8). They are initially designed for equal delays indicated by line L1 in Fig. 8. This results in yield of Y 0 for each stage (pipeline yield = Y 0 3 ). The total area for this design is the sum of the stage areas (A1+A2+A3). Now, we introduce imbalance by reducing the area of stage 1 and 3 (by dA1 and dA3) increasing their delays to line L2. This reduces yields of stages 1 and 3, to Y 1 , Y 3 (both less than Y 0 ). However, this extra area (dA1+dA3) can be added to the stage 2, thereby reducing its delay to line L3. This improves the yield of stage 2 (i.e. Y 2 > Y 0 ). If, in this particular case,

3 1 2 3 0 Y Y Y
Y the overall pipeline yield improves. This trend has been observed in the 3 stage pipeline circuits as shown in Fig. 7(b). However, introducing excess imbalance in stage delays, we might get diminishing returns when pipeline performance is governed by the µ of the slowest stage (Fig. 7(b), worst case unbalancing results). Hence, it is necessary to appropriately apply imbalance among the stages. Using the above observations, we have developed a simple heuristic to provide imbalance among the stages as shown below: else increase in small area results in large improvement in delay endif endfor (14) In order to improve the yield of a design with minimum impact on area, the delays of the stages with (R i < 1) should be more effective to reduce (small area penalty). On the contrary, to reduce area (or power) with minimum penalty on yield, the area of the stages with (R i > 1) should more effective to reduce (small delay increase).

Pipeline design flow under statistical delay distribution

Using the models developed in section 2 and based on the observations in section 3, we propose a statistical design method for the complete pipelined circuit. The proposed method is targeted to optimize area (hence, power) while ensuring the yield. In a conventional pipeline design flow, individual stages are designed and optimized independently of others before they are glued together to form a pipeline.

Under statistical delay variation, a global optimization of the complete pipeline is necessary for two reasons:

Yield target (Y) for the complete pipeline may not be met if some of the stages fail to meet their yield target (Y 0 ). In that case, the delay of the stages meeting yield target Y 0 can be further improved to compensate for the stages failed to meet the yield.

Although the yield constraint is satisfied for each stage, the integrated pipelined circuit has opportunity for optimization with respect to the area/power (section 3.2).

Global optimization of the pipeline for target yield

The optimization problem for minimizing area of a pipeline under a given yield constraint can be formulated as below, Minimize L i and U i are the minimum and maximum size factors */ To solve the above problem, we use gate-level sizing algorithm for minimizing total area under yield constraint given in [START_REF] Choi | Novel Sizing Algorithm for Yield Improvement under Process Variation in Nanometer Technology[END_REF]. It is an iterative low-complexity algorithm based on Lagrangian Relaxation (LR) [START_REF] Choi | Novel Sizing Algorithm for Yield Improvement under Process Variation in Nanometer Technology[END_REF]. We have developed a global design optimization algorithm to solve the above problem efficiently (Fig. 9). Application of the proposed algorithm directly to the complete pipeline, where all the stages are sized simultaneously, is computationally very expensive. This is because of the fact that, a pipelined circuit can be very large and will take inordinately large run time and storage to converge [START_REF] Choi | Novel Sizing Algorithm for Yield Improvement under Process Variation in Nanometer Technology[END_REF]. The algorithm employs the principle of divide-and-conquer where we size one stage at a time in such a way that the target yield for the complete pipeline is satisfied while the total area is minimized. It should be noted that, for a pipelined design with a target yield Y with respect to a delay of T TARGET , the yield constraint to each stage is now set to Y (for T TARGET ). Moreover, statistical delay analysis is performed over the complete pipeline, although the sizing is done for only one stage. It helps to make the algorithm computationally efficient, since we avoid application of the sizing routine on all the stages simultaneously. It is worth noting that, improvement in the total area of the pipeline strongly depends on the ordering at which the stages are chosen for sizing. In our algorithm, the ordering is based on the position of each stage in their area vs. delay curve as explained in section 3.2 (14). This minimizes the total design area to meet the target pipeline yield during global optimization, since stages with lower R i are sized before stages with higher R i . For each stage, we apply the algorithm in [START_REF] Choi | Novel Sizing Algorithm for Yield Improvement under Process Variation in Nanometer Technology[END_REF], which starts with assigning an initial delay constraint A o (pipeline delay) to the stage and iteratively find optimal sizes with respect to it. At the end of each iteration, statistical timing analysis on the complete pipeline is performed. Depending on the new µ and of the pipeline delay, the pipeline target delay is modified to A' o (steps 4-7 in Fig. 9). Since logic gates of only one stage are being sized at a time, we perform full timing analysis on one stage only to estimate its mean and the standard deviation. The proposed analytical models (section 2.2) are then used to produce the resulting pipeline delay distribution. This incremental timing analysis further improves the computational efficiency of the algorithm. The global optimization algorithm proposed here is significantly faster and takes much less storage compared to the case where all the stages are sized simultaneously. The LR based sizing algorithm proposed in [START_REF] Choi | Novel Sizing Algorithm for Yield Improvement under Process Variation in Nanometer Technology[END_REF] has a computational complexity of O(n 2 ) where n is the number of logic gates to size. For m pipeline stages each having n gates the simultaneous sizing approach runs with a complexity of O(m 2 n 2 ) (with space complexity of O(mn)). The proposed algorithm improves the complexity to O(mn 2 ) (with space complexity of O(n)). The algorithm is applied on an example 4-stage pipelined circuit (stages are designed with ISCAS85 benchmark circuits) and the results were compared with the case where the logic stages were individually optimized using method in [START_REF] Choi | Novel Sizing Algorithm for Yield Improvement under Process Variation in Nanometer Technology[END_REF]. We observed that, preferentially incorporating imbalance among the stages using the proposed algorithm ensures target yield (Table-II) or minimizes area for a target yield (Table-III). The highlighted rows in the tables show the stages chosen by the algorithm for yield improvement with very small area increase. Similarly, the plain rows (Table-II and III) show stages chosen by the algorithm for area saving with small penalty in yield. Overall a 9% improvement in yield can be obtained with a small area penalty 2% (Table-II). On the other hand, about 8.4% area improvement can be obtained for the same yield (table-III).

Conclusions

We have investigated pipeline delay distribution under interand intra-die parameter variations. Analytical models for estimating yield of a pipelined circuit are presented. We have observed the impact of logic depth (or number of pipe stages) and imbalance among the stage delays on the variability of the pipeline delay. An efficient sizing algorithm for pipeline to minimize area under yield constraint is presented. Our investigations show that a statistical design of a complete pipeline (not only the individual stages) is effective to improve yield in presence of parameter variations.

= 1 Figure 1 :

 11 Figure 1: A simple 5-stage pipeline and its throughput for (a) static delay model, (b) statistical delay model

Figure 2 :

 2 Delay distribution of a 12-stage inverter chain pipeline (stage logic depth = 10) under process variation with (a) only random Intra-die variation, (b) only Inter-die variation, (c) Inter-and Intra-die variation with both random and systematic components

  (a)) (b) only inter-die variations (i.e. stage delays are perfectly correlated) (Fig.2(b)) and (c) both inter and intra-die distributions with spatial correlation (i.e. stage delays are partially correlated) (Fig.2(c)). Close match have also been observed for several other pipeline configurations (Table

Figure 4 :Figure 3 :Figure 5 :

 435 Figure 4: Range of permissible mean and standard deviation for each stage to meet a target yield

Figure 6 :

 6 Figure 6: A 3-stage pipeline with different stage delays

  change of area with delay for each stage R D For each stage if R reduction in large area results in small increase in delay .

  SD i stands for delay of the i-th stage *i is the size factor for a logic gate.

Figure 7 :Figure 8 :

 78 Figure 7: Effect of unbalancing on a) pipeline delay, b) yield

Figure 9 :

 9 Figure 9: Algorithm for the global optimization of pipeline for ensuring yield with minimum area overhead

Table - II Ensuring Y TARGET (80%) with small area penalty

 -

	Stage	Individually Optimized	Proposed Algorithm
	Logic	Area (%)	Yield (%)	Area (%)	Yield (%)
	c3540	47.4	86.3	47.3	86
	c2670	25.7	95	27.4	99.1
	c1980	20.4	95	20.7	95.6
	c432	6.5	95	6.6	99.2
	Pipeline:	100	73.9	102	80.5

Table - III Area reduction for a target yield (80%)

 -

	Stage	Individually Optimized	Proposed Algorithm
	Logic	Area (%)	Yield (%)	Area (%)	Yield (%)
	c3540	50	94	45	90.5
	c2670	23.2	95	21.2	99.1
	c1980	20.3	95	19.1	90.5
	c432	6.5	94.5	6.3	99.2
	Pipeline:	100	80.3	91.6	80.5

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

* The work is sponsored in part by Marco Gigascale Systems Research Center (GSRC),


S D Cotofana

J T J Van Eijndhoven

A M Molnos
email: molnos@natlab.research.philips.com

M J M Heijligers

Compositional Memory Systems for Multimedia Communicating Tasks

  L'archive ouverte pluridisciplinaire

Introduction

The system's high performance and predictability are the required characteristics for state-of-the-art embedded multimedia applications. To meet their performance requirements parallel processing of data on multi-processor architectures is needed. The low cost demands of this application domain make the use of general purpose architectures with several GHz clock frequencies impossible.

The applications' increase in size and complexity demand hardware platforms with large guaranteed memory bandwidth. The increase in bandwidth to off-chip memories is not growing as fast as the increase in speed of computation units. A possible approach to cope with this processor-memory gap is to use memory hierarchies (caches) [START_REF] Hennesy | Computer Architecture: A Quantitative Approach[END_REF]. For embedded systems caches induce undesired unpredictability because tasks can influence each other's performance by flushing each other's data out of the cache. In this way the system is not compositional so the overall performance cannot be predicted based on the system's parts performance and changing of components has design implications to the overall system causing longer time-to-market with economical implications.

In the real-time domain one could think of replacing caches with scratch pad memories, but in practice removing caches is not a good choice because they provide flexibility. They don't need to be resized or modified if the application changes. Detailed analysis regarding the memory requirements of multimedia applications' code and data at every change of standards or features takes much effort and negatively affects the time-to-market.

The importance of managing caches in general was identified already in the single processor with multiple threads domain. Several cache partitioning and management techniques are proposed for time-sharing environment [START_REF] Suh | Dynamic cache partitioning for simultaneous multithreading systems[END_REF], instruction controlled, compile time configured data cache [START_REF] Muller | Caches with compositional performance[END_REF], first level compiler-based partitioned cache [START_REF] Mueller | Compiler support for software-based cache partitioning[END_REF]. However, for the multiprocessor domain there is a lack of cache management methods. This article presents a new cache management technique applicable for multimedia applications when running on multiprocessor architectures that achieve systems' performance compositionality. The method allows tuning the memory system according to different optimization criteria. Our proposal is based on the observation that the following holds true: in a multiprocessor system, the levels of unified cache shared between processors are the most affected by the inter-task conflicts. Thus we propose to achieve compositionality by exclusively assigning sets of the last level of cache to tasks and to the communication buffers. An accurate analytical formulation and a practical approximation to find the cache partitioning ratio for throughput and power optimization are presented. The influence of task to processor assignment in this formulation is also discussed. This paper presents a continuation of [START_REF] Molnos | Compositional memory systems for data intensive applications[END_REF] where the case of non-communicating multimedia applications is discussed.

To evaluate the performance implications of our method we consider two applications (two instances of jpeg decoder plus one canny edge detection algorithm and an mpeg2 decoder) running on a multiprocessor CAKE platform [START_REF] Stravers | Homogeneous multiprocessing and the future of silicon design paradigms[END_REF] with 4 VLIW processors and 512KB of L2 cache. The simulations indicate that when the proposed paradigm is applied to the L2 cache the system becomes compositional and also the number of L2 misses is diminished up to 5 times for the first applications and 6.5 times for the second application.

The outline of this paper is as follows. Section 2 overviews the state of the art in the field of cache management for multitasking environments. Section 3 discuss the proposed method to induce compositionality together with the optimization possibilities and an approach to determine the static cache partitioning ratio. Section 4 and 5 present the experimental framework and results. Finally, in Section 6 conclusions are drawn.

Related work

In the single processor domain different hardware ( [START_REF] Suh | Dynamic cache partitioning for simultaneous multithreading systems[END_REF], [START_REF] Stone | Optimal partitioning of cache memory[END_REF], [START_REF] Kirk | Smart (strategic memory allocation for realtime) cache design[END_REF]), software ( [START_REF] Mueller | Compiler support for software-based cache partitioning[END_REF]) or mixed ( [START_REF] Muller | Caches with compositional performance[END_REF]) cache management methods were proposed. In [START_REF] Suh | Dynamic cache partitioning for simultaneous multithreading systems[END_REF] the problem of cache efficiency for simultaneous threads in a time-sharing environment is tackled using dynamic column caching, in a besteffort way. Based on their number of misses tasks are dynamically "stealing" each other cache ways, such that the overall number of misses is improved. In [START_REF] Stone | Optimal partitioning of cache memory[END_REF] the problem of optimal allocation of cache between two competing processes that minimizes the overall miss rate is presented. In [START_REF] Kirk | Smart (strategic memory allocation for realtime) cache design[END_REF] a method to divide a cache into partitions for each realtime task and a larger partition called the shared pool for the non-real-time tasks is described. In [START_REF] Mueller | Compiler support for software-based cache partitioning[END_REF] the cache is partitioned among tasks at compile and link time. In [START_REF] Muller | Caches with compositional performance[END_REF] a new data cache organization is proposed. The L1 direct mapped cache can be partitioned and configured compile time and controlled by specific cache instructions at run time, bringing considerable performance gain.

The main differences between our work and the previous research are that we tackle the case of applications running on platforms with unified set-associative cache shared between the processors. On such platforms, in our opinion without explicitly taking into account the inter-task communication, no compositionality can be achieved and performance cannot be guaranteed. None of the above approaches can be straightforwardly extended to our case. Column caching used [START_REF] Suh | Dynamic cache partitioning for simultaneous multithreading systems[END_REF] and [START_REF] Stone | Optimal partitioning of cache memory[END_REF] is based on ways allocation from every cache set to the tasks. This partitioning type severely restricts the granularity of cache allocation to the associativity of the cache. In [START_REF] Kirk | Smart (strategic memory allocation for realtime) cache design[END_REF] data access to the shared structures do not have predictable cache access. The method in [START_REF] Mueller | Compiler support for software-based cache partitioning[END_REF] targets the L1 cache, requires a specific compiler and it is rather inflexible. The method in [START_REF] Muller | Caches with compositional performance[END_REF] is limited to L1, direct mapped data caches for uniprocessors.

Cache partitioning for system compositionality

The platform model for the proposed method is a homogeneous on-chip multiprocessor having high bandwidth communication network with the partitionable shared unified on-chip cache. On this platform parallel tasks that communicate through the memory system are executed. Resource contention for the shared buses and the shared levels of cache are generating unpredictable inter-task conflicts. The communication network has high bandwidth so the resource contention there is low. Since the levels of cache private to each processor are usually small and task switching rate in multimedia application is typically low enough, the first levels of cache can be considered private to each task. The method to induce compositionality to the system is by allocating to the tasks and the communication buffers exclusive parts of the level of cache that is shared through the processors.

The communication buffers addresses in cache are accessed by multiple tasks. If every task has a private cache part, assigning one buffer for example in the consumer's cache will allow the producer to unpredictably evict some consumer cached data. If all the buffers are allocated in a shared cache pool they will still unpredictably evict each other data, depending on their addresses and the timing of the productions and the consumptions. In order to have no cache interactions between tasks the communication buffers need their own exclusive cache partitions.

Let's now assume that the cache allocated for every communication buffer is smaller than its size and the production and consumption happens in parallel. If the producer is slow enough not to flush former produced data out of the cache the consumer will have a hit. Depending on the producing/consuming rates the number of misses varies. So, in order to have predictability, for the communication buffers one of the following should be made:

-ensure that all accesses are hits (except cold misses).

-ensure that all accesses miss in the cache.

-predict the number of misses at the buffer level -this is rather difficult because requires cycle accurate producer and consumer rates which are difficult to estimate and usually not constant in a multiprocessor environment.

Optimization opportunities

In a compositional environment the performance of one task is independent of the other tasks and this is achieved by cache partitioning. Different optimization criteria, statically prioritizing tasks and guaranteeing performance are possible by tuning the cache partitioning ratio.

An application consisting of tasks and ¡ communication buffers can be represented as a graph ¢ £ ¥ § © , where

§ £ " $ & ( 0 1 2 2 2 5 6 1 8
is the set of nodes denoting the tasks and

£ A 0 E " $ 0 E & ( 0 1 2 2 2 5 6 1 8
is the set of edges denoting inter-task communication. An edge A 0 E is present in the application's graph if the task exchange data with the task E . We consider the case of homogeneous multiprocessors and the set of processors in the system is

I £ Q R " $ R & ( 0 1 2 2 2 V 6 1 8
. Typically a multimedia application executes for an infinite time in a periodic manner. Every task has assigned the following functions: X ¥ expressing its allocated cache size and a ¥ © X ¥ expressing its execution time for one application period.

In a multiprocessor platform the execution time of every task is specified by a set d £ f h " $ i q r s 0

h & ( 0 1 0 2 2 t q 8
of time slices when executes on different processors. Depending on the cache size X ¥ the overall execution time of task becomes:

a ¥ © X ¥ £ w E r y q a ¥ © Q X ¥ © X ¥
where a ¥ © Q X ¥ © X ¥ is the execution time of slice and Q X ¥ denotes the processor executing slice .

For the multimedia applications a good optimization criteria is the execution throughput. The execution's throughput can be defined as the number of application's complete executions in a time unit. Let's denote with ¥ Q R the time processor Q R is used for the completion of one application execution. The throughput is given by the

$ r 8 ¥ ¥ Q R . Let § R d §
be the subset of tasks that execute at least one of their slices on Q R , then the processor execution time becomes:

¥ Q R £ w i q r s E r y q f g h $ E 8 & a ¥ © Q R © X ¥ i l m i h n i o p q
where l m i h n and o p q are the amount of time when Q R does task switching respectively is idle.

In an environment which allows task migration and dynamic scheduling policies the task set § R executed on Q R and the slices d are execution dependent so ¥ Q R cannot be accurately computed. In order to have an exact analytical model to reason about the overall throughput when having a certain cache partitioning a static assigning of tasks to the processors is required. In this way, tasks assigned to the same processor are executed overall sequentially, independent of the scheduling decisions for that processor and the processor's overall execution time will be:

¥ Q R © § R £ w i q r s a ¥ © Q R © X ¥ i l m i h n i o p q
where a ¥ © Q R © X ¥ is the execution time of task on processor Q R .

To optimize the throughput, the task to processor assignment and the cache allocation should be such that

¥ ¥ Q R © § R
is minimized. The power consumed by the system and the execution time are other possible optimization criteria for the multimedia systems. The consumed power depends by the time and the memory traffic that the system needs to complete all its tasks. Optimizing the overall execution time (respectively the number of misses) gives the most power consumptions reduction:

t v w w i q r s a ¥ © X ¥ y
The method is applicable for systems with statically scheduling and allocation of tasks.

Proposed optimization method

Because in our experimental system task migration and dynamic scheduling are allowed, in the proposed method for finding the cache partitioning ratio the throughput and power consumed are improved by optimizing the overall application's number of L2 misses (so execution time).

The problem of minimization of the total number of cache misses is formulated as a (Mixed) Integer Linear problem. Let ¥ © X ¥ be the number of misses of task with X ¥ allocated cache sets. The objective is to find X ¥ values for every task such that the overall number of cache misses is minimized:

t v w w i q r s ¥ © X ¥ y
We denote the set of valid cache sizes with

z R " R & ( 0 1 0 2 2 2 0 { 6 1
and X ¥ } z R "

. Due to implementation reasons z R can be, for example, limited to powers of two. The number of misses of task with

z R cache sets ¥ R £ ¥ © z R
can be obtained by simulation or program analysis. In our model we use and average R over the R obtained out of different simulation of task having z R cache. We use a set

R " $ & ( 0 1 2 2 2 5 6 1 0 R & ( 0 1 0 2 2 2 0 { 6 1 8
of variables that specify if the allocated cache for task is z R or not. A task can have only one cache partition allocated, thus

{ 6 1 R & ( R £ .
As a result the ILP solver produces a R " assignment such that the overall number of misses is minimized. Using R " and R the number of misses for task can be written as:

¥ © X ¥ £ { 6 1 w R & ( R R
where the allocated cache size for every task is:

X ¥ £ { 6 1 w R & ( R z R
The sum of the cache sizes allocated to all tasks and communication buffers should be in the limits of the available cache size £ . The size of the cache assigned to the communication buffer A p is X ¥ A p } z R " and this results in the following constraint:

5 6 1 w & ( X ¥ i ¤ 6 1 w p & ( X ¥ A p ¥ £

Experimental framework

Application Model

To describe the applications we use Y-chart Applications Programmers Interface (YAPI) [START_REF] De Kock | Yapi: application modeling for signal processing systems[END_REF]. The model of computation in YAPI is based on Kahn Process Networks. Such a network consists of parallel tasks that communicate through (theoretically unbounded) FIFOs. Parallelism and communication are explicit and synchronization is done implicitly at read from empty FIFO (and write in full FIFO in the practical case where FIFOs are bounded). In practical video YAPI applications the "memory-active" entities of the system are tasks, FIFOs and frame buffers, so the cache is exclusively split among them. The application model is not restricted to YAPI descriptions and can be used with any multitasking programming model with shared memory primitives, e.g., POSIX threads.

The FIFOs access predictability is achieved by allocating them cache of the same size as the FIFO size.

For the frame buffers the production and consumption is done actually sequential (this is intrinsic for the YAPI application). The frames used for the prediction are accessed only after they are completely produced. This means that by allocating an exclusive cache partition for every frame buffer the compositionality of the system is preserved because the access to this buffer is sequential.

If all the tasks share the heap another factor that influences compositionality is dynamic memory allocation. Depending on which data was previously allocated (so depending on task scheduling and mapping) a tasks' data structure will have different address, resulting in different inside task In the experiments in Section 5 we assume that the memory allocation is done during the initialization period and the overall allocation order is always the same.

Hardware platform and operating system

In the experiments a practical instance of the CAKE multiprocessor architecture is used [START_REF] Stravers | Homogeneous multiprocessing and the future of silicon design paradigms[END_REF]. The CAKE platform consists of a homogeneous network of computing tiles on a chip. Each tile contains CPUs (Trimedia and/or MIPS cores), a router (for out of tile communication) and memory banks (Figure 1). The processors are connected to memory by a fast, high-bandwidth snooping interconnection network. The on-tile memory is actually used as a L2 cache, shared between tasks, facilitating a fast access to the main memory which is outside the chip. The addressing space is linear.

Allocating sets of the L2 cache is implemented by changing the conventional index part of an address to a new index. In order to implement this, the cache has to be able to relate memory accesses to tasks and communication buffers. Ideally, each access should be labelled with a task id or buffer id. In this way the cache can translate the address according to a table indexed by this id. In the CAKE platform, the task id is kept in a register and can therefore be used directly.

There are several ways to obtain an id for communication buffers. A buffer id register could be used. This would imply additional requirements for the compiler, which then need to keep that register up to date. Alternatively, a part of the address could be used to encode the buffer id. This reduces the usable address space and also requires adaptation of the compiler for handling shared static data structures. Nevertheless, for dynamic memory allocation the partitioning can be implemented relatively straightforward by providing a dedicated malloc for shared buffers. A third approach is to keep a table with intervals of shared memory. This table needs to be loaded by the operating system. Then for every access the cache can lookup if the address has an associated buffer id.

The third alternative is chosen as for our experiments we are mainly interested in the system level aspects (e.g., inducing the compositionality property, implication in number of misses) and we are not concerned for the implementation's performance issues The third approach more generic than the others because any address range can be placed in any place in the cache. This easily allows for other experiments, like for example separating tasks' instructions, static initialized variables (data) and static uninitialized variables (bss) in the cache or sharing some cache partitions.

We have adapted the operating system, such that it manages the necessary translation tables for the cache. For this, it offers primitives of cache allocation for tasks and for shared memory.

Experimental results

We evaluate the proposed technique using as workload two applications. The first application (15 tasks) consists of two jpeg decoders [START_REF] De Kock | Multiprocessor mapping of process networks: a jpeg decoding case study[END_REF] working on different pictures formats and one line based canny edge detection algorithm. The second application (13 tasks) is a mpeg2 video decoder [START_REF] Van Der Wolf | An mpeg-2 decoder case study as a driver for a system level design methodology[END_REF]. In both examples there is a run-time operating system that has an exclusive part allocated such that it does not interfere with the application's tasks.

The instance of the CAKE platform used has four Trimedia processors, 512KB, 4 ways associative L2. On the experimental platform the application and run time system static allocated data (data and bss) is shared between tasks so in order to have predictable access, with the same consideration as for communication buffers (Section 3), exclusive cache partitions are allocated for them as well. In the last row of Tables 1 and2 (corresponding to the two examples) the allocated cache sizes are presented.

We performed simulations in order to determine the number of misses depending of the allocated cache size (for z R £ ¢

) and we applied the method in Section 3.2. The

Figure 2. Shared vs. best partitioned cache for every task and communication buffer

obtained partitioning ratios are presented for the two examples in Table 1 and Table 2. To give a measure of the intertask cache interaction in Figure 2 the number of misses for the shared cache case and partitioned cache case are presented for all tasks and communication buffers. For both examples we can observe that with only few sets of exclusive cache assigned to static allocated data a major improvement in performance is obtained. With respect to performance cache partitioning has two effects: on one hand the misses due to inter-task conflicts are alleviated but on the other hand some tasks do not benefit from all the cache space, increasing their number of misses. The balance between the two effects can be seen in Figure 2. The results also indicate that for the first application the L2 miss rate is improved from 9.46% to 2.21% and as a result the number of cycles per instruction (CPI) of every processor is reduced with approximate 20% (from 1.4 cpi to 1.1 cpi). For the mpeg2 application the L2 miss rate is improved from 5.1% to 0.8% and as a result the number of CPI of every processor is reduced with approximate 4% (from 1.7-1.8 cpi to 1.6-1.7 cpi). For this application the reduction in CPI is not so large because the used mpeg2 implementation was not optimized for Trimedia so it is more L1 and processor bounded than L2 bounded. The mpeg2 was also simulated with 1MB of shared L2 cache and for that the L2 miss rate was 0.6% with a corresponding 1.7 cpi for every processor.

Figure 3. Expected-simulated performance comparison for every task

The graphs in Figure 3 presents the difference between the number of misses for every task ( \ 4

) expected using the model from Section 3.2 and the ones obtained by simulation of the determined best cache partitioning ratio for the application. These differences give a measure of compositionality and reflect the variation in number of L2 misses due to the neglected effects like task switching, L1 and bus contention. For both examples the largest difference for a task between the expected and simulated number of misses relative to the overall simulated number of misses is 2%. The results suggest that system behaves as expected and the compositionality property is induced.

Conclusions

For multimedia applications when running on embedded multiprocessor systems we tackle the problem of non compositionality due to inter task conflicts in the shared cache. In the proposed approach compositionality is induced by partitioning the unified level of cache shared between processors using a scheme that exclusively allocates sets of cache to the tasks and inter-task communication buffers. Different optimization criteria, statically prioritizing tasks and guaranteeing performance are possible by tuning the cache partitioning ratio. For homogeneous multiprocessors the problems of finding the cache partitioning ratio for optimizing the throughput and the power are analytically formulated and the influence of task to processor assignment is discussed. A practical approximation of these formulations is presented and it is validated using two applications: two jpegs decoders and a canny edge detection running in parallel and a parallel mpeg2 decoder. Applying the proposed method to a CAKE instance with 4 Trimedia processors and 512KB L2 cache compositionality was achieved. The largest difference for a task between the expected and simulated number of misses relative to the overall simulated number of misses is 2%. In terms of performance, compared with the shared cache case there were 5 times less L2 cache misses resulting in 20% reduction in cycles per instruction for the first example and 6.5 times less L2 cache misses resulting in a 4% cpi reduction for the second example.

Figure 1 .

 1 Figure 1. CAKE architecture -inside-tile view

Table 1 . L2 allocated sets for 2 jpegs & canny

 1

	tasks	FrontEnd1	IDCT1	Raster1	BackEnd1
	alloc. L2 sets	4	1	32	16
	tasks	FrontEnd2	IDCT2	Raster2	BackEnd2
	alloc. L2 sets	4	1	16	16
	tasks	Fr. canny	LowPass	HorizSobel	VertSobel
	alloc. L2 sets	4	16	8		16
	tasks	HorizNMS	VertNMS	MaxTreshold
	alloc. L2 sets	8	8	4
	data	appl data	appl bss	rt data	rt bss
	alloc. L2 sets	2	2	4		4
	tasks	input	vld	hdr	isiq	memMan
	alloc. L2 sets	2	4	16	8	1
	tasks	idct	add	decMV	predict	predictRD
	alloc. L2 sets	4	4	8	16	2
	tasks	writeMB	store	output
	alloc. L2 sets	8	2	1
	data	appl data	appl bss	rt data	rt bss
	alloc. L2 sets	4	1	8	1

Table 2 . L2 allocated sets to tasks for mpeg2

 2

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)


Sarah Wassermann
email: sarah.wassermann@inria.fr

John P Rula
email: john.rula@eecs.northwestern.edu

Fabián E Bustamante
email: fabianb@eecs.northwestern.edu

Pedro Casas
email: pedro.casas@ait.ac.at

Anycast on the Move: A Look at Mobile Anycast Performance

The appeal and clear operational and economic benefits of anycast to service providers have motivated a number of recent experimental studies on its potential performance impact for end users. For CDNs on mobile networks, in particular, anycast provides a simpler alternative to existing routing systems challenged by a growing, complex, and commonly opaque cellular infrastructure. This paper presents the first analysis of anycast performance for mobile users. In particular, our evaluation focuses on two distinct anycast services, both providing part of the DNS Root zone and together covering all major geographical regions. Our results show that mobile clients tend to be routed to suboptimal replicas in terms of geographical distance, more frequently while on a cellular connection than on WiFi, with a significant impact on latency. We find that this is not simply an issue of lacking better alternatives, and that the problem is not specific to particular geographic areas or autonomous systems. We close with a first analysis of the root causes of this phenomenon and describe some of the major classes of anycast anomalies revealed during our study, additionally including a systematic approach to automatically detect such anomalies without any sort of training or annotated measurements. We release our datasets to the networking community.

I. INTRODUCTION

The impressive growth in mobile devices and use has led to an unprecedented amount of cellular traffic. The number of mobile subscriptions has grown rapidly in just a few years, surpassing 7.8 billion in the third quarter of 2017 [START_REF] Ericsson | Ericsson Mobility Report -On the Pulse of the Networked Society[END_REF]. Today, users spend most of their time browsing on their mobile phones, more than on any other device. In the United States, for instance, the average smartphone user spends 2x-3x more hours (87 hours per month) on her mobile device than she does on desktop machines (34 hours per month) [START_REF] Ofcom | The Communications Market Report[END_REF].

Most of this content is delivered to end users by content delivery networks (CDNs). CDNs deploy servers around the world and redirect clients to nearby replicas to improve performance and reliability. The process of replica selectionthe mapping of each client to a close replica -is key to CDN performance. Most commonly, it has relied on a DNS-based approach pioneered by Akamai.

Traditional replica selection systems, however, are being challenged by the rapid growth, increased complexity, and common opacity of cellular infrastructure [START_REF] Rula | Behind the Curtain: Cellular DNS and Content Replica Selection[END_REF]. Anycast offers mobile CDNs an alternative approach.

This work has been carried out while Sarah Wassermann was an MSc. student at the University of Liège.

With IP anycast, services advertise a single IP address from many physical locations (anycast sites) and clients' requests are directed, based on BGP routing policies, to a "nearby" replica [START_REF] Abley | Operation of Anycast Services[END_REF]. The approach is being used for redirecting clients in a range of applications, from naming (e.g., root servers and top-level domain resolvers) to content delivery. By letting BGP control request routes, anycast routing obviates the need for fine-grained infrastructure or client information. BGP routing also offers some degree of robustness, adapting to changes in service or network availability, and allows for some policy control. The benefits of anycast to service providers have motivated a number of recent experimental studies on its potential performance impact for end users. Prior analyses have shown that anycast routing can be suboptimal [START_REF] Schmidt | Anycast Latency: How Many Sites Are Enough?[END_REF], unstable [START_REF] Bret | TCP over IP Anycast -Pipe Dream or Reality[END_REF], [START_REF] Hiebert | Determining the Cause and Frequency of Routing Instability with Anycast[END_REF], and seemingly chaotic [START_REF] Bellis | Researching F-root Anycast Placement Using RIPE Atlas[END_REF], [START_REF] Kuipers | Analysing the K-root Anycast Infrastructure[END_REF], as routing policies have not only technical motivations, but could be dictated by political or commercial reasons. Routing changes can silently shift traffic from one site to another with a consequent loss of shared state and potential performance impact [START_REF] Wei | Does Anycast Hang up on You?[END_REF]. Yet, despite the growing dominance of mobile Internet access, no previous studies have evaluated the effectiveness of anycast for mobile users, including its routing behavior and performance.

In this paper, we take a first look at anycast performance for mobile users. In particular, our evaluation focuses on two distinct anycast services, K-and F-Root DNS, each providing part the DNS Root zone and together covering all major geographical regions. Similarly to previous work [START_REF] Kuipers | Analysing the K-root Anycast Infrastructure[END_REF], [START_REF] Li | Longitudinal Analysis of Root Server Anycast Inefficiencies[END_REF], [START_REF] Liu | Two Days in the Life of the DNS Anycast Root Servers[END_REF], we use the geographic distance as our evaluation metric. Geographic distance is a useful metric to estimate RTT [START_REF] Huffaker | Distance Metrics in the Internet[END_REF], and IP-to-location mapping techniques such as GeoPing [START_REF] Padmanabhan | An Investigation of Geographic Mapping Techniques for Internet Hosts[END_REF] and Constraint-Based Geolocation (CBG) [START_REF] Gueye | Constraint-Based Geolocation of Internet Hosts[END_REF] rely on the correlation between network delay and geographic distance. In the specific case of wireless networks hosting users on the move, we consider geographic distance to be a more suitable proximity metric than delay. Indeed, delay is influenced by multiple factors which can vary greatly, even if the user does not move much. Note that, while we are interested in the impact and behavior of delay towards user-requested content, geographic distance provides a more stable ground truth to use as baseline for the comparison of measurements. We show that mobile clients tend to be routed to suboptimal replicas in terms of geographical distance, more frequently while on a cellular connection than on WiFi, with a significant impact on perceived service performance.

In the following section, we provide background on cellular networks and anycast routing and review prior work. We describe our methodology and datasets for measuring anycast routing for mobile end users in Section III, and present our findings in Section IV. We further investigate the causes of aberrant anycast behavior for mobile users in Section V and present a systematic approach to identify such anomalies in an unsupervised manner through clustering techniques in Section VI. Finally, we conclude in Section VII.

II. ANYCAST AND MOBILE USERS

Previous efforts have focused on techniques for characterizing anycast deployment [START_REF] Cicalese | Characterizing IPv4 Anycast Adoption and Deployment[END_REF], enumerating [START_REF] Cicalese | A Fistful of Pings: Accurate and Lightweight Anycast Enumeration and Geolocation[END_REF]- [START_REF] Madory | Who are the Anycasters?[END_REF], and geolocating servers based on latency measurements [START_REF] Cicalese | A Fistful of Pings: Accurate and Lightweight Anycast Enumeration and Geolocation[END_REF], [START_REF]Latency-based Anycast Geolocalization: Algorithms, Software and Datasets[END_REF].

Cicalese et al. [START_REF] Cicalese | Characterizing IPv4 Anycast Adoption and Deployment[END_REF] presented the first Internet-wide anycast census and showed that major players in the Internet ecosystem are relying on anycast, even though only a small part of the IPv4 space has so far been anycasted.

Calder et al. [START_REF] Calder | Analyzing the Performance of an Anycast CDN[END_REF] analyzed the performance of anycast in the Bing CDN. The authors showed that, although anycast performed well in general, about 20% of the clients were redirected to a suboptimal replica, with a negative impact on performance. Kuipers et al. [START_REF] Kuipers | Analysing the K-root Anycast Infrastructure[END_REF] reported a related study on the performance of the K-Root DNS service using the RIPE Atlas framework [START_REF]RIPE Atlas: a Global Internet Measurement Network[END_REF]. Like Calder et al., they observed that anycast was suboptimal for multiple probes with, for instance, 46% of them being redirected to a suboptimal DNS server, both in terms of latency and geographical distance. Li and Spring [START_REF] Li | Tier-1's break Anycast DNS[END_REF] claimed that the main causes for this behavior are Tier-1 ASes, which forward almost all requests systematically to the same replica, irrespectively of the user's location.

Addressing the flaws of anycast remains a major challenge. Schmidt et al. [START_REF] Schmidt | Anycast Latency: How Many Sites Are Enough?[END_REF] analyzed the impact of the number of anycast sites on the performance of this paradigm in the DNS infrastructure. The authors concluded that increasing the number of sites does not solve the suboptimal mapping observed nor, counterintuitively, reduces clients' latencies. In [START_REF] Bellis | Researching F-root Anycast Placement Using RIPE Atlas[END_REF], Bellis et al. described their use of RIPE Atlas to detect and address issues related to F-Root, such as the redirection of requests to replicas in different continents (e.g., from a probe in Europe to a server in Atlanta).

Zarifis et al. [START_REF] Zarifis | Diagnosing Path Inflation of Mobile Client Traffic[END_REF] studied the impact of the Internet topology and routing on the performance perceived by mobile users, and revealed that a significant fraction of Internet paths are inflated, with a non-negligible performance penalty. The authors investigated the root causes for route inflation and found the lack of carrier ingress points to be one of the main reasons.

Our work builds on approaches and insights from many of these past efforts. Yet, to the best of our knowledge, this is the first work focused on the performance of anycast in mobile networks.

III. METHODOLOGY AND DATASETS

This section describes our datasets and measurement methodology for analyzing anycast from mobile devices. We focus our analysis on two major root DNS services: F-Root and K-Root. We selected these services since they are both widely replicated, covering together all geographic areas (with approximately 60 servers for K-Root and nearly 140 for F-Root) and both with publicly available site locations and unicast IP addresses. The latter point is key to let us evaluate the performance of anycast routing relative to its "optimal" (in terms of unicast) site location.

We collected active measurements from geographically distributed clients on both cellular and WiFi networks from September 2016 to April 2017, using the ALICE engine 1 . ALICE conducts mobile network measurements by executing a small, self-contained experiment script, run approximately every hour. In each experiment, for each of the two root DNS services, clients launched ping and traceroute measurements towards the root server's anycast address, as well as towards five unicast addresses of the analyzed DNS service determined to be the geographically closest to the client. Target unicast addresses were selected based on the origin country of each client's IP address, as reported by whois.

Note that, besides our previous justification for choosing geographical distance over latency as selection criteria, we consider the closest unicast replicas in terms of geographical distance and not latency also due to practical concerns. Indeed, we cannot know the closest unicast replicas in terms of latency without flooding the user's connection with latency measurements towards all available unicast replicas before starting every new experiment, which would be impractical and highly invasive. Moreover, a previous study [START_REF] Kuipers | Analysing the K-root Anycast Infrastructure[END_REF] shows that geographically distant replicas generally result in poorer latency, further supporting our methodology.

We collected data from mobile devices while clients were connected to either WiFi or cellular networks. We thus divide our data into a set containing experiments launched from cellular networks (CELL) and another one including the experiments issued from WiFi (WIFI). As measurements were collected opportunistically, based on connection availability and resource usage, the number of experiments in each of these sets is not the same.

CELL dataset. CELL includes more than 20,000 experiments, issued from 151 different clients. Our cellular users were scattered across nearly 40 different countries, with 70% of them located in the United States, Greece, Brazil, and France. Furthermore, the analyzed clients were hosted in a range of major ASes, including AS 29247 (GR-Cosmote, Greece), AS 26599 (Vivo, Brazil), and AS 22394 (Verizon Wireless, US).

WIFI dataset. WIFI encompasses three times as many experiments as CELL, issued from 251 clients. Clients were located in nearly 50 different countries around the world with 70% of them in Greece and the United States. Some of the most active clients (i.e., those that have launched the most experiments) were hosted in AS 6799 (OTENET-GR, Greece), AS 9121 (TTNet, Turkey), and AS 7922 (COMCAST, US). A total of 125 clients launched experiments from both cellular and WiFi networks.

We are making both the CELL and WIFI datasets publicly available on GitHub2 .

IV. FINDINGS AND OBSERVATIONS

In this section, we present key observations from our analysis of mobile anycast performance. In particular, our results suggest that mobile clients are most of the time routed towards suboptimal replicas in terms of geographical distance (Section IV-A), more frequently while on a cellular connection than on WiFi, and that the additional distance traveled by requests has a significant impact on performance (Section IV-B). We analyze the scope of this phenomenon (Section IV-C) and the existence of better alternatives (Section IV-D).

A. The Travel Distance Problem

Our study is partially motivated by anecdotal evidence of long geographic distances between mobile clients and their assigned anycast DNS servers. To characterize these distances, given that anycast IP addresses cannot be geolocated [START_REF] Mcpherson | Architectural Considerations of IP Anycast[END_REF], we geolocate the penultimate hop on the path from a client to the anycast server and estimate its location using the Akamai EdgeScape service 3 . We use this and the client location to compute the geographic distance between a client and her assigned anycast server. For each experiment, the client recorded her anonymized geographic location to a 10 km 2 area. Figure 1 shows the distribution of these distances, in kilometers, for both WiFi and cellular users.

Figure 1(a) presents the travel distance between clients and their K-Root anycast servers (D client→anycast ). We see a significant difference between cell and WiFi: D client→anycast is smaller than 4,000 kilometers for approximately 75% of the experiments carried out on a WiFi network, whereas this holds for merely 50% of the experiments issued from cellular clients. While we can observe the same phenomenon for F-Root on Figure 1(b), the differences are not as striking as for the K-Root service. The difference in number and geographic distribution between F-and K-Root services explains the measured differences in the distances traveled by a client request on a cellular or WiFi connection. Still, clients on cellular networks travel farther to their assigned replicas (with potential latency implications).

B. Impact on Latency

Does the additional distance traveled by users' requests result in worse user performance? We evaluate the impact of the distance between clients and their replicas in terms of latency. Figure 2 shows the distribution of their latencies. Indeed, the figure reveals that latencies are higher for the measurements carried out on cellular networks than for the ones performed on WiFi. Again, the difference is more pronounced for K-Root (Figure 2(a)) than for F-root (Figure 2(b)). A comparison of both figures shows that, while the distribution of latencies is similar for K-Root and F-Root on WiFi, we see significantly worse performance for K-Root compared to F-Root on cellular networks: 90% of the experiments output a RTT lower than 200 milliseconds for F-Root, but only 70% for K-Root DNS.

Other factors beyond distance traveled, such as signal strength and radio access latency, could result in performance degradations for clients on cellular networks. To understand the importance of geographic distance on observed latency, we look at the correlation between distance and minimum RTT for cellular measurements. Figure 3, a scatter plot of this correlation, clearly shows that, while not the only factor, the distance between cellular clients and their assigned anycast servers is indeed an important aspect to take into account when it comes to performance optimization. We observed the same behavior for WiFi clients; the graph is not included due to space constraints.

C. Where Does the Travel Distance Problem Occur?

Is it possible that the observed issues are limited to some specific geographical regions or a few particularly mismanaged ASes? To explore this, we now focus our analysis on five ASes that provide the most measurements for clients in cellular networks. This leaves us with 2,310 measurements issued from 16 users for K-Root and 1,860 measurements retrieved from 15 mobile clients for F-Root DNS.

Figure 4 presents the distributions of the distances traveled to anycast servers for clients in these ASes. While we see large variability within and across ASes, particularly in USbased networks potentially covering larger geographic areas, there are no clear patterns to argue for the identified problem to be a region-or AS-specific issue.

As noted before, F-Root DNS replicas seem to be closer to their clients, regardless of the network. Indeed, for 80% of the launched measurements towards F-Root anycast servers, D client→anycast is less than approximately 4,300 kilometers, while this distance is larger than 8,500 kilometers for 20% of the measurements issued towards K-Root servers.

D. Do Cellular Clients Have a Closer Option?

Another possible explanation for the long distances would be that these cellular clients are simply located far away from all the available replicas. Further analysis shows that this is not necessarily the case. We compute the distance between cellular clients and both their assigned anycast server (D client→anycast ) and their geographically closest unicast replica (D client→unicast ). We compute the additional distance traveled by anycast as the difference between these two distances (i.e., D client→anycast -D client→unicast ), which we denote by δ client .

Figure 5 presents the additional distance traveled by anycast requests, δ client , for the measurements launched from the five ASes providing the most measurements for K-and F-Root DNS. We can easily infer from this graph that our clients are most of the time routed to a suboptimal replica and that this issue does not seem to be specific to one AS or region. Nevertheless, comparing Figures 5(a) and 5(b) leads us to conclude that anycast routing is significantly worse for K-Root than for F-Root, even though the mappings for F-Root are far from optimal. Particularly clear examples for the K-Root problem are clients residing in the United States, but being routed towards anycast servers located in London and Tehran. This is even more surprising as there are at least seven K-Root servers distributed across the United States, with four hosted on the East Coast and two on the West Coast. Finally, we observe that, in each of the examined ASes, some mobile clients are redirected to the nearest anycast F-Root server, while the ideal case never occurs for the K-Root service.

V. WHY TRAVELING SO FAR? A FIRST LOOK INTO CELLULAR ANOMALIES

In the following paragraphs, we explore some of the causes of anycast routing anomalies for cellular clients. In particular, we focus on the measurements with a value of δ client higher than 1,000 kilometers. This corresponds to more than 80% of the experiments for K-Root and to 70% of the experiments for F-Root. We refer to these measurements as suboptimal anycast measurements.

We first compare the lengths of the AS-paths leading from the cellular clients, on the one hand, towards the anycastassigned replicas and, on the other hand, towards the nearest unicast servers. We then explore whether the announcement policies of different root nodes -whether they are locally or globally announced -affects anycast quality. Finally, we identify three classes of anomalies including: (i) distant client packet gateways, (ii) poor anycast routing within Tier-1 networks, and (iii) improper routing within cellular networks.

AS-path comparison. We analyze the lengths of the AS paths between cellular clients and their assigned anycast server, and between those clients and the geographically closest unicast replica. We find that, when anycast servers are geographically farther away than unicast servers, they still tend to be closer in terms of number of traversed ASes as one would expect with anycast routing. The paths client → anycast server are shorter than the paths client → unicast replica in nearly 40% of the measurements towards K-Root and in more than 75% towards F-Root. While more than 80% of the paths to F-Root anycast servers have a length of at most six AS hops (more than 60%, considering K-Root), approximately 55% of unicast paths for both DNS services are longer than six hops.

Root announcement policies. We also analyze the node types of both DNS services. While all but one K-Root servers are global, there are only four global servers for F-Root (one in the Netherlands and three in the United States). We find that clients from Singapore, Australia, the United States, and the Dominican Republic are routed towards global nodes in the United States, with a significant latency penalty for most of them, despite the presence of local nodes in their geographical surroundings. Note that, in cases where both a local and global node reside in the same city, we are not able to distinguish between them from traceroutes alone.

Further investigation into anomalous anycast routing for our cellular clients reveals three classes of anycast performance problems. While by no means exhaustive of all problems encountered by cellular clients, these common problems appear to be either unique, or significantly more common, to cellular with respect to fixed-line clients.

A. Distant Packet Gateway

As previously reported in [START_REF] Xu | Cellular Data Network Infrastructure Characterization and Implication on Mobile Content Placement[END_REF], a cellular client's packet gateway (PGW) largely determines the client's network position and locality, since all client traffic is routed through that PGW. This implies that client packets first encounter Internet routing once they are beyond the PGW and thus, depending on the relative distance between a client's PGW and the assigned server, the wrong assignment of a client to a distant PGW can be the determinant factor of anycast routing performance. Knowing the location of a client's PGW can therefore greatly aid in diagnosing anomalous anycast routing. We investigate the consequences of this phenomenon in terms of latency. We find that for about 35% of the considered anycast measurements, the clients' request spends more than 50% of its time on the way towards the PGW. In particular, for 45% of the traceroutes, the clients' probe needs more than 70 milliseconds to reach the assigned packet gateway. These latencies are far from being negligible and suggest that the routing within the users' cellular network could be significantly enhanced.

o] vš >} š]}vW }•š}vU D W't >} š]}vW >}• vP o •U vÇ •š •š]v š]}v E OE •š hv] •š •š]v š]}v ,}‰ /W ^E ,}‰ /W ^E írñ WOE]À š OE •• • r r írñ WOE]À š OE •• • r r ò îìôXíôñXíòìXíôî òðòí >}• vP o •U ò îìôXíôñXíòìXíôî òðòí >}• vP o •U ó îìôXíôñXíòìXíôí òðòí >}• vP o •U ó îìôXíôñXíòìXíôí òðòí >}• vP o •U ô òðXíîñXîôXðí òðòí ^ v :}• U ô òðXíîñXîôXðñ òðòí ^ v :}• U õ òðXíîñXïìXîïì òðòí ^ v :}• U õ òðXíîñXïíXíñ òðòí ^ v :}• U
A striking example is depicted in Figure 6. These paths correspond to a client located in Boston (MA), but whose PGW is situated in Los Angeles (CA). While the closest unicast replica is hosted in St. George (UT), the client's anycast request is routed to London, with a large latency penalty of almost 100 milliseconds. In this scenario, the journey of the packet is not only costly in terms of kilometers and latency due to the suboptimal geolocation mapping between the client and an anycast server, but also because the packet needs to travel more than 4,000 kilometers to reach the PGW.

A more detailed analysis of our suboptimal anycast measurements shows that, in 40% of the cases, the concerned clients are located more than 1,000 kilometers away from their assigned PGW, and, for approximately 10%, this distance is larger than 3,000 kilometers. Unfortunately, even when the PGW is far away from a cellular client, it does not mean that the assigned anycast replica is near the packet gateway. Indeed, for the users presenting a distance to their PGW larger than 1,000 kilometers, their packet still has to travel more than 2,000 kilometers to reach its destination in more than 65% of these measurements.

We take a closer look at the top five ASes in terms of number of different clients suffering from suboptimal anycast measurements. Figure 7 depicts the distribution of the distances between the clients and their PGW for each investigated anycast case. Our results show that distant PGWs do not seem to be a systematic issue: while this distance varies widely in ASes 21928 (T-Mobile, US), 22394 (Verizon Wireless, US), and 26599 (Vivo, Brazil), it is relatively stable and small for ASes 12361 (Vodafone-Panafon, Greece) and 23243 (Comcel, Guatemala). Specifically, for all the considered measurements launched in ASes 12361 and 23243, the cellular clients present a distance to their PGW smaller than 700 kilometers, which is the case for less than 1% of the measurements seen in AS 26599. For AS 21928, we note that 60% of the measurements have been carried out while clients were very close to their assigned PGW (between 30 and 60 kilometers), while the remaining 40% have been launched when they were farther away from it (up to 3,000 kilometers).

B. Tier-1 Routing

We observe that many of the instances of poor anycast performance occur when anycast paths traverse Tier-1 transit networks. While Tier-1 routing problems are not solely specific to cellular networks [START_REF] Li | Tier-1's break Anycast DNS[END_REF], we found them to be much more prevalent on suboptimal cellular paths than on WiFi ones. We found that 64% of the paths leading towards the assigned anycast replica traverse at least one Tier-1 AS, while this is the case for 73% of the paths towards the geographically closest unicast server.

With large Tier-1 networks, we see clients often routed to the same anycast replica, regardless of where clients enter these networks. Many problematic cases we investigated appear to be caused by packets remaining in the transit network until routed to a distant destination. We find this behavior in Tier-1 networks with varying levels of consistency. For example, clients entering AT&T (AS 7018) split anycast destinations between sites in Reno (NV) and London. We also see clients routed through AboveNet (AS 6461) being consistently routed to K-Root sites in London. An observed example corresponds to a cellular client residing in Los Angeles and whose PGW is in the same city. Even though the geographically closest unicast server is in Reno, her packet is routed to a server situated in London. Analyzing the corresponding traceroutes reveals that more than 50% of the hops in the path leading to the anycast server lie in AboveNet, a well-known Tier-1 AS. As opposed to the anycast path, the path connecting the cellular user to her nearest unicast server exits the AS 6461 fairly quickly: less than 25% of the traceroute hops are in this Tier-1 AS. As a consequence, the latency towards the anycast server is much higher than the one towards the closest unicast replica. In this scenario, the seen traceroutes show a similar behavior as the ones in Figure 6.

Further analysis highlights that AboveNet is by far the most often encountered Tier-1 AS. Indeed, this Tier-1 is seen in 55% of the anycast traceroute paths traversing at least one of these ASes, while the second-most frequently observed one (AS 12956, Telxius) in only 24% of these traceroutes. When looking at the cases for which we have a δ client between 5,000 and 10,000 kilometers, AboveNet is still the most popular Tier-1 (appearing in 40% of the traceroutes with Tier-1 hops), but ASes such as 7018 (AT&T) and 2914 (NTT Communications) are traversed more frequently with respect to the cases with smaller values of δ client . In addition to that, we were curious about the impact of Tier-1s on latency. Our results reveal that anycast is suboptimal in terms of latency with respect to the closest unicast replica in more than 70% of the cases when the packet traverses at least one Tier-1 AS on its way to the assigned anycast server.

C. Improper Cellular Network Routing

We found cases where paths leaving towards the anycast and nearest unicast servers diverge immediately at the IP level after exiting a client's PGW, even though they both end up in the same AS for the next IP hop. We discovered that 17% of our analyzed cases suffer from improper routing within the cellular network. However, this issue does not seem to be tied to specific ISPs or regions, as we detected it in multiple ASes located across Europe and the United States.

Figure 8 illustrates this phenomenon. In the figure, the IP paths diverge right after they exit the client's PGW. These diverging paths remain in the operator's same AS (AS 27699), and even have the same next hop AS. In this scenario, the client's packet reaches the unicast replica significantly faster than the assigned anycast server; we observe an impressive latency difference of 151 milliseconds. While we cannot know what caused this exact instance, we saw for several operators multiple transit providers connected to, or very near, cellular network PGWs. As we noted in the previous class, the choice of transit can play a large role in the behavior of anycast routes. This is especially true if the provider is a Tier-1 network, as many commonly are for large cellular networks.

VI. TOWARDS AUTOMATIC CHARACTERIZATION OF ANYCAST CELLULAR ANOMALIES

We devote the last section to further investigate the aforementioned suboptimal anycast measurements (i.e., with δ client higher than 1,000 kilometers) and corresponding anomalies. We take an exploratory approach, and use machine-learning techniques to provide first steps in the automatic identification and characterization of anycast anomalies in cellular networks. Given the general lack of ground truth regarding the nature and root causes of the anomalies, we perform an unsupervised analysis, relying on clustering techniques.

We take two different approaches: firstly, we use the well-known DBSCAN clustering algorithm [START_REF] Ester | A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise[END_REF] to identify homogeneous groupings and outliers within the suboptimal measurements. Then, we apply a more advanced, DBSCANbased algorithm we have previously conceived in [START_REF] Casas | Unsupervised Network Intrusion Detection Systems: Detecting the Unknown without Knowledge[END_REF] to automatically spot out the most relevant anomalies within the measurements, without requiring any sort of training or labeled data.

A. DBSCAN-based Analysis

By applying DBSCAN to the suboptimal anycast measurements, we are able to identify four consistent clusters C i=1..4 (silhouette scores above 0.7) and a group of outliers, i.e., measurements not belonging to any of the four clusters. Following the results obtained in [START_REF] Casas | Unsupervised Network Intrusion Detection Systems: Detecting the Unknown without Knowledge[END_REF], we use auto-calibration for DBSCAN parameters ρ (distance search coefficient) and λ (minimum connected region size) [START_REF] Ester | A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise[END_REF], taking ρ = α × n (α = 0.01) and λ = β × d (β = 0.2), where n is the number of measurements and d the average distance between all the different pairs of measurements. Each measurement is described by a vector of 12 features: δ client , client to anycast RTT (ARTT) and distance (CAD), client to closest unicast RTT (URTT), client to PGW RTT (CPRTT) and distance (CPD), PGW to anycast RTT and distance, anycast IP-and AS-path length (IPPL and ASPL), number and fraction of Tier-1 hops in the anycast path.

Figure 9 presents an overall characterization of the results. According to Figures 9(b) and 9(c), anomalous or highly suboptimal anycast measurements as compared to unicast are located in clusters C 2 , C 3 and within the outliers, representing altogether more than 80% of all the suboptimal measurements. Indeed, despite the higher geographical distance of the anycast replicas, C 1 and C 4 measurements tend to have a similar performance between anycast and unicast, and both clusters consist of measurements with comparable latencies, trading closer PGWs with farther anycast replicas between clusters.

Previous anomaly types (e.g., distant PGWs and routing/path-inflation issues) are clearly observed in C 2 , C 3 , and the outliers. Outliers systematically correspond to poor anycast performance as compared to unicast, with large latencies, far away replicas with long paths (potentially linked to routing and path inflation issues), PGW selection at distant locations, and very high differences in terms of geographical distance to selected anycast versus unicast replicas. C 3 measurements are characterized by large δ client values and much higher latencies as compared to unicast replicas. Finally, C 2 measurements are more spread and correspond to mostly worse-than-unicast performance scenarios, also showing far located PGWs and anycast replicas, but resulting in RTTs below 250 milliseconds for more than 85% of the measurements.

B. Anomaly Detection with Sub-Space Clustering

To conclude, we explore the possibility of automatically spotting out the aforementioned anomalies to better support the root cause analysis process. In particular, we apply a fully unsupervised approach for anomaly detection (UAD) previously conceived in [START_REF] Casas | Unsupervised Network Intrusion Detection Systems: Detecting the Unknown without Knowledge[END_REF] to the set of suboptimal anycast measurements. In a nutshell, the UAD algorithm uses DB-SCAN to cluster the measurements by relying on different sub-spaces of the complete feature space, and accumulates a weighted distance between outliers and clusters on these sub-spaces to compute an abnormality score; the higher this score, the more different (i.e., anomalous) is the corresponding measurement from the rest. DBSCAN parameters are set through the same auto-calibration approach described before. We refer the reader to [START_REF] Casas | Unsupervised Network Intrusion Detection Systems: Detecting the Unknown without Knowledge[END_REF] for further details on the algorithm.

Figure 10 depicts the resulting abnormality scores obtained by UAD, sorted in descending order. The number of identified outliers corresponds to approximately 1% of the suboptimal anycast measurements. There are different regions clearly visible in the ranking, with knees and breaks showing different levels of abnormality. To shed some light on the detected anomalies, Table I reports the corresponding feature values for the top-10 ranked anomalies.

The first two anomalies are characterized by very large RTTs to both anycast and unicast replicas; anomalies with IDs 4, 6, 7, 8 correspond to similar measurements, characterized by a very high δ client . Interestingly, anomalies 5 and 9 correspond to some of the anomalies manually studied in Section V (cf. Figures 6 and8). Remaining anomalies 3 and 10 correspond to far located anycast replicas and reasonably close PGWs, suggesting anomalies linked to routing. In all cases, the performance of anycast in terms of latency is far worse than that of unicast, with a δ client higher than 2,000 kilometers.

VII. CONCLUSIONS

In this paper, we presented findings from the first look at anycast performance for mobile users. Using data collected from a crowdsourced platform, we showed that mobile clients are very frequently mapped to a geographically suboptimal anycast replica, with a significant impact on performance in terms of latency. We found that the long distances between clients and their assigned anycast server are not due to a lack of better, closer unicast replicas, and that this phenomenon is not bound to specific regions or particular ASes. Exploring root causes highlighted three classes of anycast anomalies, namely distant client packet gateways, poor anycast routing within Tier-1 networks, and improper routing out of cellular networks. We finally presented a clustering-based analysis of the suboptimal measurements, including a fully unsupervised, automatic approach to identify the most critical cases. In addition to that, we release our datasets to the networking community. In ongoing work, we are exploring some of these issues, including patterns of cellular specific anycast problems, their role on CDN request routing, and the impact of these pitfalls on clients' quality of experience.

Fig. 1 :

 1 Fig. 1: Travel distance from clients to anycast servers.

  Latency to F-Root.

Fig. 2 :

 2 Fig. 2: Latency from clients to anycast servers.

Fig. 3 :

 3 Fig. 3: Correlation between geographic distance and minimum RTT.

  D client→anycast observed for F-Root.

Fig. 4 :

 4 Fig. 4: D client→anycast in top 5 ASes (w.r.t. number of measurements).

Fig. 5 :

 5 Fig. 5: δ client in top 5 ASes (w.r.t. number of measurements).

Fig. 7 :

 7 Fig. 7: Distance client → assigned PGW in top 5 ASes (w.r.t. number of different clients).

Fig. 8 :

 8 Fig. 8: Improper cellular network routing.

  δ client . (b) anycastRT T -unicastRT T . (c) anycastRT T . IP path length. (e) client → PGW distance.(f) PGW → anycast distance.

Fig. 9 :

 9 Fig. 9: Characterization of the suboptimal anycast measurements using DBSCAN clustering. Highly suboptimal anycast measurements (as compared to unicast) are located in clusters C 2 , C 3 and within the outliers.

Fig. 10 :

 10 Fig. 10: Anomalies detected by sub-space clustering.

TABLE I :

 I Top 10 detected anomalies.

	ID	δ client (km)	ARTT (ms)	URTT (ms)	IPPL	CAD (km)	CPD (km)
	1	2,387	6,828	2,431	12	3,871	138
	2	5,022	1,638	1,428	12	8,040	1,457
	3	5,027	229	97	12	8,042	1,454
	4	13,097	273	95	14	14,582	1,084
	5	3,714	271	119	12	7,392	2,118
	6	13,073	250	77	14	14,567	1,099
	7	13,084	256	76	14	14,573	1,093
	8	13,063	261	90	14	14,559	1,108
	9	4,519	258	160	18	5,281	4,185
	10	9,597	322	138	13	10,316	848

http://aqualab.cs.northwestern.edu/projects/261-alice

https://github.com/SAWassermann/mobile-anycast

https://www.akamai.com/us/en/products/web-performance