{'fraction_non_alphanumeric': 0.03823601270461764, 'fraction_numerical': 0.01606401172733936, 'mean_word_length': 4.241382990075766, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Skewless Network Clock Synchronization Without Discontinuity: Convergence and Performance

This paper examines synchronization of computer clocks connected via a data network and proposes a skewless algorithm to synchronize them. Unlike existing solutions, which either estimate and compensate the frequency difference (skew) among clocks or introduce offset corrections that can generate jitter and possibly even backward jumps, our solution achieves synchronization without these problems. We first analyze the convergence property of the algorithm and provide explicit necessary and sufficient conditions on the parameters to guarantee synchronization. We then study the effect of noisy measurements (jitter) and frequency drift (wander) on the offsets and synchronization frequency, and further optimize the parameter values to minimize their variance. Our study reveals a few insights, for example, we show that our algorithm can converge even in the presence of timing loops and noise, provided that there is a well defined leader. This marks a clear contrast with current standards such as NTP and PTP, where timing loops are specifically avoided. Furthermore, timing loops can even be beneficial in our scheme as it is demonstrated that highly connected subnetworks can collectively outperform individual clients when the time source has large jitter. The results are supported by experiments running on a cluster of IBM BladeCenter servers with Linux.


I. INTRODUCTION
Keeping consistent time among different nodes in a network is a fundamental requirement of many distributed applications. Nodes' internal clocks are usually not accurate enough and tend to drift apart from each other over time, generating inconsistent time values. Network clock synchronization allows these devices to correct their clocks to match a global reference of time, such as the Universal Coordinated Time (UTC), by performing time Preliminary version of this paper appears in IEEE ICNP 2013 [1]. arXiv:1405.6477v2 [math.OC] 28 Jul 2014 measurements through a network. For example, for the Internet, network clock synchronization has been an important subject of research and several different protocols have been proposed [2]- [8]. These protocols are used in various applications with diverse precision requirements such as banking, communications, traffic monitoring and security.
In modern wireless cellular networks, for instance, time-sharing protocols need an accuracy of several microseconds to guarantee the efficient use of channel capacity. Another example is the recently announced Google Spanner [9], a globally-distributed database, which depends on globally-synchronized clocks within at most several milliseconds drifts.
The current de facto standard for IP networks is the Network Time Protocol (NTP) proposed by David Mills [2].
It is a low-cost, purely software-based solution whose accuracy mostly ranges from hundreds of microseconds to several milliseconds. On the other hand, IEEE 1588 (PTP) [4] gives superior performance by achieving submicrosecond or even nanosecond accuracy. However, it is relatively expensive as it requires special hardware support to achieve those accuracy levels and may not be fully compatible with legacy cluster systems.
Newer synchronization protocols have been proposed with the objective of balancing between accuracy and cost.
For example, IBM Coordinated Cluster Time (CCT) [10] is able to provide better performance than NTP without additional hardware. Its success is based on a skew estimation mechanism [11] that progressively adapts the clock frequency without offset corrections. Another solution that achieves this objective is the RADclock [5], [8] which decouples skew compensation from offset corrections by decomposing the clock into a high performance difference clock for measuring time differences and a less precise absolute clock that provides UTC time.
There are two major difficulties that make the network clock synchronization problem challenging. Firstly, the frequency of hardware clocks is sensitive to temperature and is constantly varying. Secondly, the latency introduced by the OS and network congestion delay results in errors in the time measurements which can be propagated through the network. Thus, most protocols introduce different ways of estimating the frequency mismatch (skew) [11], [12] and measuring the time difference (offset) [13], [14] while maintaining a simple network topology [2], [4]. This leads in particular to extensive literature on skew estimation [12], [15]- [17] which suggests that explicit skew estimation is necessary for clock synchronization. This paper takes a different approach and shows that using skew estimation is unnecessary. We provide a simple algorithm that is able to compensate the clock skew without any explicit estimation of it. Our algorithm only uses current offset information and an exponential average of the past offsets. Thus, it neither needs to store long offset history nor perform expensive computations on them. The solution provided in this paper achieves microsecond level accuracy without requiring any special hardware. Since we do not explicitly estimate the skew, the implementation is simpler and more robust to noise than IBM CCT, and does not introduce offset corrections, which avoids the need of decomposing the clock into several components to reduce jitter as in RADclock.
By looking at the synchronization problem from a new angle, this paper also provides several new insights. For example, a common practice in the clock synchronization community is to avoid timing loops in the network [2, p. 6] [4, p. 16, s. 6.2]. This is because it is thought that timing loops can introduce instability as stated in [2]: "Drawing from the experience of the telephone industry, which learned such lessons at considerable cost, the subnet topology... must never be allowed to form a loop." Even though for some parameter values loops can produce instability, we show that a set of proper parameters can guarantee convergence even in the presence of loops. Furthermore, we experimentally demonstrate in Section VI that timing loops among clients can actually help reduce the jitter of the synchronization error and is therefore desirable.
The rest of the paper is organized as follows. In Section II we provide some background on how clocks are actually implemented in computers and how different protocols discipline them. Section III motivates and describes our algorithm together with an intuitive explanation of why it works. In Section IV, we analyze the convergence property of the algorithm and determine the set of parameter values and connectivity patterns under which synchronization is guaranteed. The parameter values that guarantee synchronization depend on the network topology, but there exists a subset of them that is independent of topology and therefore of great practical interest. The effect of noisy measurement and wander is studied in Section V, together with an optimization procedure that finds optimal parameter values. Experimental results evaluating the performance of the algorithm are presented in Section VI.
We conclude in Section VII.

II. COMPUTER CLOCKS AND SYNCHRONIZATION
Most computer architectures keep their own estimate of time using a counter that is periodically increased by either hardware or kernel's interrupt service routines (ISRs). On Linux platforms for instance, there are usually several different clock devices that can be selected as the clock source by changing the clocksource kernel parameter.
One particular counter that has recently been used by several clock synchronization protocols [5], [10] is the Time Stamp Counter (TSC) that counts the number of CPU cycles since the last restart of the system. For example, in the IBM BladeCenter LS21 servers, the TSC is a 64-bit counter that increments every δ o = 0.416ns since the CPU Based on this counter, each server builds its own estimate x i (t) of the global time reference, UTC, denoted here by t. For example, if c i (t) denotes the counter's value of computer i at time t, then x i (t) can be computed using where x o is the estimate of the time when the server was turned on (t 0 ).
Thus, synchronizing computer clocks implies correcting x i (t) in order to match t, i.e. enforcing x i (t) = t. There are two difficulties on this estimation process. Firstly, the initial time t 0 in which the counter starts its unknown.
Secondly, the counter updating period δ i (δ i ≈ δ 0 ) is usually unknown with enough precision and therefore presents a skew r i = xi(t)−xi(t0) t−t0 = δ 0 δi . This is illustrated in Figure 1a where x i (t) not only increases at a different rate than t, but also starts from a value different from t 0 , represented by x o i . In practice, c i (t) can be approximated by a real value since the time between increments is extremely small (0.416ns) and the maximum count register value so large (2 64 − 1) that it would take more than 200 years to reach.
Therefore, x i (t) can be described by the linear map of the global reference t, i.e.
where s o i is an additional skew correction implemented to compensate the skew. Equation (2) also shows that if one can set s o i = 1/r i and x o i = t 0 , then we obtain a perfectly synchronized clock with x i (t) = t. The main problem is that not only neither t 0 nor r i can be explicitly estimated, but also r i varies with time as shown in Figure 2a. Thus, current protocols periodically update s o i and x o i in order to keep track of the changes of r i . These updates are made using the offset between the current estimate x i (t) and the global time t, i.e.
Offset and relative skew measurements , and the relative frequency error that is computed using two offset measurements separated by τ seconds, i.e.
Figure 1b provides an illustration of these measurements. In most protocols (see e.g. [2], [5], [10]) (3) goes through an additional filtering process to reduce the estimation noise. Here we will use f err i (t k ) to denote either the measurement obtained using (3) or a filtered version of it. To understand the differences between current protocols, we first rewrite the evolution of x i (t) based only on the time instants t k in which the clock corrections are performed. We allow the skew correction s o i to vary over time, i.e. s i (t k ), and write x i (t k+1 ) as a function of x i (t k ). Thus, we obtain where τ = t k+1 − t k is the time elapsed between adaptations; also known as poll interval [2]. The values u x i (t k ) and u s i (t k ) represent two different types of corrections that a given protocol chooses to do at time t k and are usually implemented within the interval (t k , t k+1 ). u x i (t k ) is usually referred to as offset correction and u s i (t k ) as skew correction. 1 See Figure 2b for an illustration of their effect on the linux time.
Remark 1: One of the implicit assumptions of the model (4) is that we require every server to update their clocks simultaneously at time instances {t k }. This may seem unrealistic since its implementation would require sharing a common time reference which is the whole purpose of the algorithm. However, the analysis presented in Section IV can be extended for pseudo-synchronous implementations as proposed in [18] where each node measures the offset with their neighbors and updates whenever x i (t) = kT .
We now proceed to summarize the different types of adaptations implemented by current protocols. To simplify the comparison, we assume that each server can connect directly to the source of UTC time (t). This assumption will be dropped in Section III after we describe our solution. The main differences between current protocols lies on whether they use offset corrections, skew corrections, or both, and whether they update using offset values D x i (t k ), relative frequency errors f err i (t k ), or both.

A. Offset corrections
These corrections consist in keeping the skew fixed and periodically introducing time changes of size u where κ 1 , κ 2 > 0. They are used by NTPv3 [19] and NTPv4 [2] respectively under ordinary conditions. These protocols have in general a slow initialization period as shown in Figure 3a. This is because the algorithm must first obtain a very accurate estimate of the initial frequency error f err i (t 0 ). Furthermore, these updates usually generate non-smooth time evolutions as in Figures 3b and 4a, and should be done carefully since they might introduce backward jumps (x i (t k+1 ) < x i (t k )), which can be problematic for some applications.

B. Skew corrections
Another alternative that avoids using steep changes in time is proposed by the IBM CCT solution [10]. This alternative does not introduce any offset correction, i.e. u x i (t k ) = 0, and updates the skew . The behavior of this algorithm is shown in Figure 4b. In [20] it was shown for a slightly modified version of ) the algorithm can achieve synchronization for very diverse network architectures. 1 These corrections can be implemented in Linux OS using the adjtimex() interface to update the system clock or by maintaining a virtual version of x i (t) and directly applying the corrections to it, as in IBM CCT [10] and RADclock [5]. The latter gives more control on how the corrections are implemented since it does not depend on kernel's routines. However, the estimation of f err i (t k ) is nontrivial as it is constantly changing with subsequent updates of s i (t k ) and it usually involves sophisticated computations [11], [12].

C. Skew and offset corrections
This type of corrections allow dependence on only offset information D x i (t k ) as input to u x i (t k ) and u s i (t k ). For instance, in [6] was proposed. This option allows the system to achieve synchronization without any skew estimation. But the cost of achieving it, is introducing offset corrections in x i (t) as shown in Figure 4c. Therefore, it suffers from the same problems discussed in II-A.
Another alternative that falls in into this category is the RADclock [5]. In this solution the offset correction u x i (t k ) is an exponential average of the past offsets and the skew compensation u s i (t k ) is a filtered version of f err i (t k ). The exponential average of offsets and filter stage in f err i (t k ) allows this solution to mitigate the jumps and become more robust to jitter. However, it does not necessarily prevent backward jumps unless the offset corrections are smaller than the precision of the clock.

III. CONTINUOUS SKEWLESS SYNCHRONIZATION
We now present an algorithm that overcomes the limitations of the solutions described in Section II. In other words, our solution has the following two properties: 1) Continuity: The protocol does not introduce steep changes on the time value, i.e. u x i (t k ) ≡ 0. 2) Skew independence: The protocol does not use skew information f err i (t k ) as input. A solution with these properties will therefore prevent unnecessary offset corrections that produce jitter and will be more robust to noise by avoiding skew estimation. After describing and motivating our algorithm, we show how the updating rule can be implemented in the context of a network environment.
The motivation behind the proposed solution comes from trying to compensate the problem that arises when one tries to naively impose properties 1) and 2), i.e. using Figure 5 shows that this type of clock corrections is unstable; the offset D x i (t k ) of the slave clock oscillates with an exponentially increasing amplitude.
The oscillations in Figure 5 arise due to the fundamental limitations of using offset to update frequency. This is better seen in the continuous time version of the system (4) with (5), i.e.
If we consider the offset D x i = t − x i (t) as the system state, then we havė . This is analogous to a spring mass system without friction. Thus, it has two purely imaginary eigenvalues that generate sustained oscillations; see [7], [21] for similar examples. 2 One way to damp these oscillations in the 2 In the discrete time system the oscillations increase in amplitude since there is a delay between the time the offset is measured t k and the time the update is made t k+1 which makes the system unstable.
spring-mass case is by adding friction. This implies adding a term that includes a frequency mismatch f err i (t) in our system, which is equivalent to the protocols of Section II-B, and therefore undesired.
However, there are other ways to damp these oscillations using passivity-based techniques from control theory [22]. The basic idea is to introduce an additional state y i that generates the desired friction to damp the oscillations.
Inspired by [22], we consider the exponentially weighted moving average of the offset and update x i (t k ) and s i (t k ) using: Figure 5 shows how the proposed strategy is able to compensate the oscillations without needing to estimate the value of f err i (t k ). The stability of the algorithm will depend on how κ 1 , κ 2 and p are chosen. A detailed specification of these values is given in Section IV-B.
Finally, since we are interested in studying the effect of timing loops, we move away from the client-server configuration implicitly assumed in Section II and allow mutual or cyclic interactions among nodes. The interactions between different nodes is described by a graph G(V, E), where V represents the set of n nodes (i ∈ V ) and E the set of directed edges ij; ij ∈ E means node i can measure its offset with respect to j, D . Within this context, a natural extension of (6)-(7) is to substitute D x i (t k ) with the weighted average of i's neighbors offsets. Thus, we propose the following algorithm to update the clocks in the network.
Algorithm 1 (Alg1): For each computer node i in the network, perform the following actions: -Compute the time offsets (D x ij (t k )) from i to every neighbor j at time t k . -Update the skew s i (t k+1 ) and the moving average y i (t k+1 ) at time t k+1 according to: where N i represents the set of neighbors of i and the weights α ij are positive.
Equation (8) can be interpreted as a discrete-time second-order consensus algorithm with an additional smoothing in which, besides using position information (time estimates x i (t k )), we use a smoothed version of the position errors (y i (t k )) to control speed (s i (t k )). Consensus algorithms have been a subject of intensive research since the seminar work of Jadbabaie et al. [23], see e.g. [22] and references therein. In particular, application of consensus ideas to computer clock synchronization can be found in [6] (second order consensus) and [18] (first order consensus).
Thus, the analysis presented in this paper also contributes to this rich literature by characterizing convergence of discrete-time consensus algorithms.
When using our algorithm, many servers can affect the final frequency of the system. Thus, when the system synchronizes, we have r * and x * are possibly different from their ideal values 1 and t 0 . Their final values depend on the initial condition of all different clocks as well as the topology, which we assume to be a connected graph in this paper.
Differences with RADclock: Although (8) seems to be similar to RADclock [5], there are some key differences that affect their behavior.
1) Even though both solutions used an exponentially weighted offset estimate, our filtering (8c) does not depend on the estimated offset error as in [5]. Moreover, while RADclock uses it to make offset corrections (changing u x i (t k )), we use our weighted offset measurement y i (t) to make skew correction (changing u s i (t k )). Therefore, neither the measurement itself nor its use are the same.
2) RADclock explicitly uses offset measurements to introduce correction on the offset (u x i (t k )) and an estimation of the skew to compensate it (u s i (t k )). Our algorithm only compensates the skew by using the last measured offsets D x ij (t k ) and our filtered offset measurement y i (t k ). Thus, we have neither explicit estimation of the skew nor explicit compensation of the offset, which makes synchronization rather unintuitive.
Notation: We use 0 m×n (1 m×n ) to denote the matrices of all zeros (ones) within R m×n and 0 n (1 n ) to denote the column vectors of appropriate dimensions. I n ∈ R n×n represents the identity matrix. Given a matrix

IV. CONVERGENCE ANALYSIS
We now analyze the asymptotic behavior of system (8) and provide a necessary and sufficient condition on the parameter values that guarantee its convergence to (9). Throughout this section we shall assume that the internal skew r i of each clock is constant and that the offset measurements D x ij (t k ) can be obtained without incurring in any error. These assumptions will be relaxed in Section V.
The key insight of our analysis comes from decomposing the system (12) into two complementary systems that keep track of two different physical properties. In particular, we will use the scalars to track the average behavior of the system and to track how each individual clock deviates from the collective mean.
left eigenvector of the zero eigenvalue of the Laplacian matrix L ∈ R n associated with G(V, E), i.e.
ξi ri . While in general ξ may not be unique, it becomes unique when G(E, V ) is connected.
It will be more convenient sometimes to use a vector form representation of (8) given by where whereÂ The convergence analysis of this section is done in two stages. First, we provide necessary and sufficient conditions for synchronization in terms of the eigenvalues of A (Section IV-A) and then use Hermite-Biehler Theorem [24] to relate these eigenvalues with the parameter values that can be directly used in practice (Section IV-B). All the proof details are included in the appendix for interested readers.

A. Asymptotic Behavior
We start by studying the asymptotic behavior of (12). That is, we are interested in finding under what conditions the series of elements {x i (t k )} converge to (9) as t k goes to infinity.
We will show that we can study (12) by looking at the evolution of (13). In particular we will show that (9) is equivalent to Consider the Jordan normal form [25] of where J = blockdiag(J l ) l∈I(A) , ζ i and η i are the right and left generalized eigenvectors of A such that The following lemmas will allow us to connect the behavior of (13) with (12). Furthermore, µ l (A) are the roots of where ν l = µ l (τ LR) and satisfies and (η 2 , η 1 ) respectively, associated with µ(A) = 1 and the eigenvectors ζ 3 and η 3 associated with µ(A) = 1 − p are given by Moreover, given ζ l = [x T l s T l y T l ] T , with l > 3, then the following conditions must be satisfied The proof of Lemmas 1 and 2 can be found in Appendices A-A and A-B.
Lemmas 1 and 2 also provide further information of the structure of J in (16). That is, J = blockdiag(Ĵ 1 ,Ĵ 2 ) Moreover, direct application of (21) shows that N = P [blockdiag(0 3×3 , I 3(n−1) )]P −1 which implies that A and N have the same eigenvectors and thereforê Similarly, it is easy to see that the Jordan normal form ofÃ is given bỹ Therefore, under the conditions of Lemmas 1 and 2 δz k andz k capture the behavior of a complementary set of eigenvalues of A. We are now ready to state our main convergence result.
Moreover, whenever the system synchronizes, we have The proof of Theorem 1 can be found in Appendix B. Theorem 1 provides an analytical tool to understand the influence of the different nodes of the graph in the final offset x * and frequency r * . For example, suppose that we know that node 1 has perfect knowledge of its own frequency (r 1 ) and the UTC time at t = t 0 (x 1 (t 0 ) = t 0 ), and configure the network such that node 1 is the unique leader like the top node in Figures 6a and 6c. It is easy to show that ξ 1 = 1 and ξ i = 0 ∀i = 1. Then, using (25a)-(25b) and definition of γ we can see that γ = r 1 and However, since node 1 knows r 1 and t 0 , it can choose x 1 (t 0 ) = t 0 , s 1 (t 0 ) = 1 r1 and y 1 (t 0 ) = 0. Thus, we obtain x * = t 0 and r * = 1 which implies by (9) that every node in the network will end up with x i (t) = t. In other words, Theorem 1 allows us to understand how the information propagates and how we can guarantee that every server will converge to the desired time. Notice that the initial condition used for server 1 is equivalent to assuming that server 1 is a reliable source of UTC like an atomic clock for instance.

B. Necessary and sufficient conditions for synchronization
We now provide necessary and sufficient conditions in terms of explicit parameter values (κ 1 , κ 2 ,τ and p) for Theorem 1 to hold. We will restrict our attention to graphs that have Laplacian matrices with real eigenvalues. This includes for example trees (Figure 6a), symmetric graphs with α ij = α ji (Figure 6b) and symmetric graphs with a leader (Figure 6c).
The proof consists on studying the Schur stability of g l (λ) and has several steps. We first perform a change of variable that maps the unit circle onto the left half-plane. This transforms the problem of studying the Schur stability into a Hurwitz stability problem which is solved using Hermite-Biehler Theorem which says: Given the polynomial P (s) = a n s n +...+a 0 , let P r (ω) and P i (ω) be the real and imaginary part of P (jω), i.e. P (jω) = P r (ω)+jP i (ω).
Then P (s) is a Hurwitz polynomial if and only if 1) a n a n−1 > 0 and 2) The zeros of P r (ω) and P i (ω) are all simple and real and interlace as ω runs from −∞ to +∞.  The proof of Theorem 2 can be found in Appendix C. Note that although µ max depends on r i which is in general unknown, it is easy to show that µ l (LR) ≤r max µ l (L) wherer max is an upper bound of the maximum rate deviation r i . Furthermore, using Greshgorin's circle theorem, it is easy to show that µ max (L) ≤ 2α max := 2 max i α ii .
Therefore, if we set convergence is guaranteed for every connected graph with real eigenvalues.

V. NETWORK DELAYS AND CLOCK WANDER
In the previous section we showed the in the absence of network delays and clock wander, the system was able to achieve synchronization on a wide variety of communication topologies. In other words, we assumed the internal clock skew r i was fixed and that each computer could measure its offset with a neighbor D x ij (t k ) without incurring in any error. We now study the behavior of our system when such assumptions are no longer true. We will model both, network delays and clock drifts using noise processes.
Network Delays: Since our algorithm does not perform skew estimation, the network errors only affect the offset measurements D x ij (t k ) in (8). We use g w ij w ij (t k ) to denote the error incurred in estimating the offset between nodes i and j at time t k , i.e. we replace D (8). This can be produced for instance by a congested connection between the two different nodes or due to path delay asymmetries [5]. We assume that w ij (t k ) has stationary mean E [w ij (t k )] =w ij ∀t k and unit variance E[(w ij (t k ) −w ij ) 2 ] = 1 and use g w ij to weigh the different connections.
Clock Wander: We model the clock wander as a stochastic input to the clock skew adaptation. That is, instead of (4b), s i (t k ) changes according to where d i (t k ) is a random variable with zero mean E[d i (t k )] = 0 and unit variance E[d i (t k ) 2 ] = 1 and g d i a positive scalar use to model clock heterogeneity. Equation (27) can be derived from a linear approximation of (4) with a time varying internal skew r i (t k ) driven by an auto regressive process [12], i.e. r i (t k+1 ) = (1−q i )r i +q i r i (t k )+g d i d i (t k ). We omit the details of the derivation due to space constraints.
This motivates the study of the stochastic process where B − G = min{B G , 0 n×m } and B G being the incidence matrix of G(V, E) 3 and w k = [w ij (t k )] T . The matrix C maps the system state z k to the performance metric v k and will be specified in Section VI.
One interesting difference between network delays (w ij (t k )) and wander (d i (t k )) is that in order to obtain good performance the algorithm should reject the noise from network delays w ij (t k ), but compensate the skew fast enough to follow d i (t k ).
In the remaining of this section, we first study the effect of biased network noise (w ij = 0) in the asymptotic frequency of the system and time offsets. In particular, we show that for arbitrarily distributed noise with stationary mean, the system's frequency tends to constantly drift unless there is a well defined leader in the topology. We then proceed to study how the parameters and network topology affect the systems performance, which is represented by the output signal v k of the stochastic process.

A. Frequency Drift and Time Offsets
Here, we concentrate on studying the evolution of the first moment of the stochastic process (28). That is, we want to understand how E[z k ] evolves as k → +∞. This is equivalent to study (28) in the case when the noise input e k is constant e k =ē = [d TwT ] T , whered = 0 n by definition.
Therefore, we will focus on understanding the effect of a constant inputē on (28). Again, we will usez k to understand howē affects the collective behavior of the clocks and δz k to understand how each individual clock deviates from the collective. The next theorem summarizes the effect of non-zero mean error on the collective behavior of the system.

Theorem 3 (Frequency Drift):
In the presence of noise and under the condition of Theorem 1 the collective frequencys(t k ) will constantly drift away from its mean with probability one (in the set of possiblew), unless the graph G(V, E) has a unique leader 4 . Whenever G(V, E) does have a leader, the mean frequency r * is given by (25b).
The proof of Theorem 3 can be found in Appendix D.
Remark 2: Theorem 3 provides a precise characterization of how network delays and loops can produce instabilities. Similar results can be obtained for any protocol that controls clock speeds based on neighbors information. Therefore, this theorem shows that current industry practice is conservative and allows us to explore a wider set of topologies with timing loops (provided that such loops avoid the leader).
We now show how the deviations δz k are affected byē.
where L † is the pseudo inverse of L and δw : w. The proof can be found in Appendix D.

B. H 2 Performance Optimization
We now proceed to study the effect of noisy measurements and wander on the output standard deviation of the system (||v k || 2 ) when the input e k is white noise (E[e k e T l ] = I m+n δ(l − k)). In other words, we seek to find parameter values that minimize Since in practice we want to avoid any frequency drift introduced by the noise we will consider only topologies with a well defined leader. Thus, all the randomness of the system is concentrated in δz k and we limit to study the stochastic process This optimization problem is standard in the control theory community and it can be shown to be equivalent to min X,κ1,κ2,p,αij where f (κ 1 , κ 2 , p, α ij ) := trace[XBB T ],Â is a function of (κ 1 , κ 2 , p, α ij ) and ρ * < 1. The constrain (30b) has been added in order to maintain the stability ofÂ.
While it is not in general easy to find the global minimum of (30) there has been intensive research in designing optimization algorithms that find local minimums of the H 2 norm of continuous time [26] and discrete time [27] systems. In this work we solve (30) using a discrete-time version of the package Hifoo [26], [28] known as Hifood [29]. Several adaptations were needed to use Hifood to solve (30). The details of these changes are documented in the Appendix E.
The output of this optimization problem will be used in Experiment 6 of next section to demonstrate that the standard belief that "clock precision degrades with the number of hops" is not necessarily true. To test our solution and analysis, we implement an asynchronous version of Algorithm 1 (Alg1) in C using the IBM CCT solution as our code base. Each server issues a thread to handle the connection with each neighbor.
Every τ seconds (using OS time) each client takes offset measurements with its assigned neighbor and reports it to the main thread. Similarly, the main thread wakes up every τ seconds and gathers the offset information from all the connections and performs the update described in (8). We do not perform any explicit filtering of offset values, besides discarding spurious offsets larger than 500ms in comparison with previous measurement. 5 Our program reads the TSC counter directly using the rdtsc assembly instruction to minimize reading latencies and maintains a virtual clock that can be directly updated. The list of neighbors is read from a configuration file and whenever there is no neighbor, the program follows the local Linux clock. Finally, offset measurements are taken using an improved ping pong mechanism proposed in [10].
We run our skewless protocol in a cluster of IBM BladeCenter LS21 servers with two AMD Opteron processors of 2.40GHz, and 16GB of memory. As shown in Figure 7, the servers serv1-serv10 are used to run the protocol.
The offset measurements are taken through a Gigabit Ethernet switch.
Server serv0 is used as a common reference. It runs the same program that implements Alg1, but without skew adaptations, to measure the offset between itself and the other servers (serv1-serv10). These measurements are obtained through a 10Gbps Cisco 4x InfiniBand Switch to minimize network latencies. Since the offset measurements performed by serv0 are done at different instances for different servers, we use linear interpolation to compensate this error. To compute the offset between two servers, say serv1 and serv2 (x 1 (t) − x 2 (t)), we subtract offset measurements obtained form serv0, i.e. (x 1 (t) − x 0 (t)) − (x 2 (t) − x 0 (t)). Finally, we also eliminate spurious measurements that generate offsets bigger than 1ms as these are clearly due to network or os latencies.
We use this testbed to validate the analysis in sections IV and V. First, we illustrate the effect of different parameters and analyze the effect of the network configuration on convergence (Experiment 1). Then we present a series of configurations that demonstrate how connectivity between clients is useful in reducing the jitter of a noisy clock source (Experiment 2). We compare the performance of our protocol with respect to NTP version 4 (Experiment 3) and IBM CCT (Experiment 4). Finally, we verify the presence frequency drift in the absence of a leader (Experiment 5), and study the interplay between network delays, wander and parameter values (Experiment 6).
The output performance signal v k will be the vector of offset differences between the leader 1 and every other .., n}. We will use a normalized version of it, referred here as mean relative deviation √ S n , as a performance metric. To make these comparison fair among different servers we correct our performance value by the empirical mean deviation to compensate biases due to path asymmetries. In other words, where < · > amounts to the sample average. We will also use the 99% Confidence Interval CI 99 and the maximum offset (CI 100 ) as metrics of accuracy. For example, if CI 99 = 10µs, then the 99% of the offset samples will be within 10µs of the leader.
Unless explicitly stated, the default parameter values are The scalar c is a commit or gain factor that will allow us to compensate the effect of τ . Notice that by definition of α ij , α ii = c for every node that is not the leader.
For fixed polling interval τ , the stability of the system depends on the value of µ max , which is determined by the underlying network topology and the values of α ij .

Experiment 1 (Convergence):
We first consider the client server configuration described in Figure 8a with a time step τ = 1s. In this configuration µ max ≈ c = 0.7 and therefore condition (iii) becomes τ < 1.2717s. Figure 9a shows the offset between serv1 (the leader) and serv2 (the client) in microseconds. There we can see how serv2 gradually updates s 2 (t k ) until the offset becomes negligible.  the servers. This is in fact true provided that network is a directed tree as in Figure 6a. The intuition behind this fact is that in a tree, each client connects only to one server. Thus, those connected to the leader will synchronize first and then subsequent layers will follow.
However, once loops appear in the network, there is no longer a clear dependency since two given nodes can mutually get information from each other. This type of dependency might make the algorithm unstable. Figure 9b shows an experiment with the same configuration as Figure 9a in which serv2 synchronizes with serv1 until a third server (serv3) appears after 60s. At that moment the system is reconfigured to have the topology of Figure 8b introducing a timing loop between serv2 and serv3. This timing loop makes the system unstable.
The instability arises since after serv3 starts, the new topology has µ max ≈ 1.5c = 1.05. Thus, the time step condition (iii) becomes τ < 847.8ms which is no longer satisfied by τ = 1s.
This may be solved for the new topology (Figure 8b) by using any τ smaller than 847.8ms. However, if we want a set of parameters that is independent of the topology, we can use (26) and notice that α max = c andr max ≈ 1.
We choose τ = 500ms < 890.2 2α max ms = 890.2 2c ms = 635.9ms. Figure 10 shows how now serv2 and serv3 can synchronize with serv1 after introducing this change.  We run Alg1 on 10 servers (serv1 through serv10). The connection setup is described in Figure 11. Every node is directly connected unidirectionally to the leader (serv1) and bidirectionally to 2K additional neighbors.

K=0 K=2
Notice that the arrow only shows a dependency relationship, the ping pong mechanism sends packets in both direction of the physical communication. We used a value of Jitter max = 10ms. Since the error was introduced in both directions of the ping pong, this is equivalent to a standard deviation of 6.05ms.  Figure 12 illustrates the relative offset between the two extreme cases; The star topology (K = 0) is shown in Figure 12a, and the complete subgraph (K = 4) is shown in Figure 12b.
The worst case offset for K = 0 is CI 100 = 5.1ms which is on the order of the standard deviation of the jitter.
However, when K = 4 we obtain a worst case offset of CI 100 = 690.8µs, an order of magnitude improvement. The change on the mean relative deviation √ S n as the connectivity among clients increases is studied in Figure   13. The results presented show that even without any offset filtering mechanism the network itself is able to perform a distributed filtering that achieves an improvement of up to a factor of 6.26 in √ S n , or equivalently a noise reduction of almost 8dB.

Experiment 3 (Comparison with NTPv4):
We now perform a thorough comparison between our protocol (Alg1) and NTPv4. We used a one hop configuration using serv1 as leader running an NTPv4 server and Alg1, and serv9 and serv10 as clients, connected only to serv1, running NTPv4 and Alg1 respectively.
In order to make a fair comparison, we need both algorithms to use the same polling interval. Thus, we fix τ = 16sec. This can be done for NTP by setting the parameters minpoll and maxpoll to 4 (2 4 = 16secs).
The remainder parameter values for Alg1 were obtained with our optimization framework (using g w ij = 100 and g d ij = 1e − 3) and are given by p = 1.98, κ 1 = 1.388 and κ 2 = 1.374.
(34) Figure 14a shows the time differences between the clients running NTPv4 and Alg1 (serv9 and serv10) , and the leader (serv1) over a period of 60 hours. It can be seen that Alg1 is able to track serv1's clock keeping an offset smaller than 5µs for most of the time while NTPv4 incurs in larger offsets during the same period of time. This difference is produced by the fact that Alg1 is able to react more rapidly to frequency changes while NTPv4 incurs in more offset corrections that generate larger jitter.  The mean offsets of Alg1 and NTPv4 are −0.48µs and 9.00µs. This difference in mean is mainly due to an asymmetric path on the measurements from serv9 to serv1. After compensating this bias, Alg1 achieves a performance of √ S n = 1.3µs, CI 99 = 4.9µs and a maximum offset of CI 100 = 20.6µs, while NTPv4 obtains √ S n = 6.4µs, CI 99 = 74.5µs and a maximum offset of CI 100 = 1.4ms. Thus, not only Alg1 achieves a reduction of √ S n by a factor of 5.0 (−7dB) with respect to NTPv4, but it also obtains smaller confidence intervals and maximum offset values.
A more detailed and comprehensive analysis is presented in Figure 15 where we plot the Cumulative Distribution Function (CDF) and Probability Density Function of the samples. The improvement of Alg1 with respect to NTPv4 is again clearly seen here. Finally, we investigate the speed of convergence. Starting from both clients synchronized to server serv1, we introduce a 25ms offset. Figure 14b shows how Alg1 is able to converge to a 20µs range within one hour while NTPv4 needs 4.5hours to achieve the same synchronization precision. In Figure 16a we present the mean relative deviation √ S n for two clients connected directly to the leader as the jitter is increased from Jitter max = 10µs to Jitter max = 160µs, doubling Jitter max each time, with a granularity in the random generator of 1µs. The worst case offset is shown in Figure 16b. Each data point is computed using a sample run of 250 seconds.
Our algorithm consistently outperforms IBM CCT in terms of both √ S n and worst case offset. The performance improvement is due to two reasons. Firstly, the noise filter used by the IBM CCT algorithm is tailored for noise distributions that are mostly concentrated close to zero with sporadic large errors. However, it does not work properly in cases where the distribution is more homogeneous as in this case. Secondly, by choosing δκ = κ 1 −κ 2 = 0.002 1 the protocol becomes very robust to offset errors.

Experiment 5 (Frequency drift without leader):
We now proceed to experimentally verify that without leader, the system tends to constantly drift the frequency. Our analysis predicts that even the minor bias in the offset measurements will produce this effect. To verify this phenomenon, we use the network topology in Figure 8b with τ = 0.5s and wait for the system to converge.
After 1000s the timing process of serv1 is turned off. Figure 17 shows how the offsets of serv2 and serv3 start to grow in a parabolic trajectory characteristic of a constant acceleration, i.e. constant drift. After 6600s serv1 is restarted and the system quickly recovers synchronization. A second order fit of the faulty trajectory was perform   the optimal parameter values depend on the different noise conditions within the network described in Figure 18.
We consider three different noise scenarios in which we either add jitter between server serv1 and servers serv2 and serv3, and/or add wander on severs serv2-serv7. In all the cases we used τ = 0.5s and make offset measurements through the InfiniBand switch to minimize the any additional source of noise.  We used different values of g w ij and g d i to differentiate the noise conditions in the optimization scheme. The large jitter scenario is represented in by g d i = 1e − 3 ∀i, g w 21 = g w 31 = 100 and g w ij = 1 otherwise. The large wander scenario is represented by g d i = 1e − 1 ∀i and g w ij = 1. Finally, the large jitter and wander scenario is represented using g d i = 1e − 1 ∀i, g w 21 = g w 31 = 100 and g w ij = 1 otherwise. The output parameter values for all three cases are present also in Figure 18.   Figure 19 shows the standard deviation of the offset between servers serv2-serv7 and serv1 in the three experi-mental scenarios and for the three different set of parameters shown in Figure 18. It can be seen that although the configuration tuned for jitter performs very well in cases with large jitter, it performs quite poorly in scenarios with large wander. Similarly, the configuration tuned for wander does not perform well in high jitter scenarios. However, the configuration tuned for jitter and wander is able to provide acceptable performance in all three experimental scenarios. Thus, we experimentally demonstrate a fundamental tradeoff between jitter and wander.
Finally, it is interesting to notice that due to the fact the optimization is solved using v i (t k ) = x i (t k ) − x 1 (t k ) as performance metric, the choice of parameters does not degrades the performance of each clock with the hop count.

VII. CONCLUSION
This paper presents a clock synchronization protocol that is able to synchronize networked nodes without explicit estimation of the clock skews and steep corrections on the time. Our solution is guaranteed to converge even in the presence of timing loops which allow different clients to share timing information and even collectively outperform individual clients when the time source has large jitter. The system is robust to noisy measurements and wander provided that the topology has a well defined leader, and we can optimize the parameter values to minimize noise variance. We implemented our solution on a cluster of IBM BladeCenter servers and empirically verified our predictions and our protocol's supremacy over some existing solutions.
Further evaluation of our protocol is needed. In particular, we are interested in comparing our solution with other protocols such as PTP and RADclock, as well as studying its robustness under stressed scenarios. Another interesting direction is to devise a distributed algorithm, exploiting our optimization framework, that can adapt the parameter values depending on the network condition.
3) =⇒ 1): The algorithm achieves synchronization whenever (9) holds. Then, it follows from (12) and (9) that asymptotically the system behaves according to Thus, since P := [ζ 1 ... ζ 3n ] is invertible, its columns ζ l are linearly independent. Therefore, if the system synchronizes for arbitrary initial condition, then it must be the case that the effect of the remaining modes µ l (A) vanishes, which can only happen if for every µ l (A) = 1, |µ l (A)| < 1 and the multiplicity of µ l (A) = 1 is two, i.e.
Finally, to obtain (25) we use a similar computation to the one of Lemma 1 to show thatP andP −1 in (24) are given byP Thus, sinceÃ k =PĴ k

1P
−1 a direct computation shows that Therefore the multiplicity of µ(A) = 1 is two and by (18) these are the roots of g n (λ) = (λ − 1) 2 (λ − 1 + p), which corresponds to the case ν n = 0. Thus, to satisfy Theorem 1 we need to show that the remaining eigenvalues are strictly in the unit circle. This is true for the remaining root of g n (λ) if and only if (i).
For the remaining g l (λ), this implies that are Schur polynomials. Thus, we will show that g l (λ) is a Schur polynomial if and only if (i)-(iii) hold. We drop the subindex l for the rest of the proof.
We first transform the Schur stability problem into a Hurwitz stability problem. Consider the change of variable λ = s+1 s−1 . Then |λ| < 1 if and only if R[s] < 0. Now, since ν > 0 by (18), let We will apply Hermite-Biehler Theorem to P (s), but first let us express what 1) and 2) of the Theorem mean here.
Therefore, since κ 1 = κ 2 ,s(t k ) constantly drifts unless Finally, there are two different scenarios in which (42) can be satisfied.
That is −ξ T B − G diag[α ij g w ij ]w = ξ 1 0 = 0 2) G does not have a well defined root: Thus, there are at least two nodes with ξ i = 0 andw is such that However, 2) is only satisfied by a set of values ofw with zero measure. Thus, there should be a unique leader for synchronization.

APPENDIX E
H 2 OPTIMIZATION USING HIFOOD The software package Hifood [29] does not solve (30) directly. Instead, it solves:  Using these definitions it is straight forward to verify that (A 1 + B 2 KC 2 ) =Â, B 1 + B 2 KD 21 =B and C 1 =Ĉ.
The main difficulty in solving (30) instead of (48) is that our controller K is a nonlinear function of the parameters K(κ 1 , κ 2 , p, α) and cannot be readily obtained using (48). Furthermore, the main source of nonlinearity comes from the products κ 1 diag[α ij ] and pdiag[α ij ]. This structure is not currently supported by traditional software distributions, which usually only support sparsity patterns, and therefore needs to be implemented.
Fortunately, Hifood only uses gradient information in their implementation of BGS and gradient bundle stages.
Thus, to implement discrete time H 2 optimization we generated a new Matlab subroutine that evaluated the H 2 norm f as well as its gradients.
The evaluation of the gradient is performed in three stages using the chain rule. We first compute the gradients of f with respect toĀ,B andC which are given by where Y is the solution to Y =ĀYĀ T +BB T .
Once ∇Āf , ∇Bf and ∇Cf are computed we can use the subroutines of hifood to compute ∂Ā ∂K , ∂B ∂K and ∂C ∂K . Finally, we obtain and similarly for other parameters.