Anne Benoit

Aurélien Cavelan

Valentin Le Fèvre

Yves Robert

Keywords: resilience, failures, checkpoint, checkpointing period, replication, heterogeneous platforms résilience, faute, checkpoint, période de checkpoint, réplication, plates-formes hétérogènes

In this paper, we design and analyze strategies to replicate the execution of an application on two different platforms subject to failures, using checkpointing on a shared stable storage. We derive the optimal pattern size W for a periodic checkpointing strategy where both platforms concurrently try and execute W units of work before checkpointing. The first platform that completes its pattern takes a checkpoint, and the other platform interrupts its execution to synchronize from that checkpoint. We compare this strategy to a simpler on-failure checkpointing strategy, where a checkpoint is taken by one platform only whenever the other platform encounters a failure. We use first or second-order approximations to compute overheads and optimal pattern sizes, and show through extensive simulations that these models are very accurate. The simulations show the usefulness of a secondary platform to reduce execution time, even when the platforms have relatively different speeds: in average, over a wide range of scenarios, the overhead is reduced by 30%. The simulations also demonstrate that the periodic checkpointing strategy is globally more efficient, unless platform speeds are quite close.

Optimal checkpointing period with replicated execution on heterogeneous platforms 1 Introduction

One of the most important challenges faced by large-scale computing systems is the frequent occurence of failures (a.k.a. fails-top errors) [START_REF] Cappello | Toward exascale resilience: 2014 update[END_REF][START_REF] Snir | Addressing failures in exascale computing[END_REF]. Platform sizes have become so large that failures are likely to strike during the execution of an application. Consider the mean time between failures µ (usually denoted as MTBF) of a platform with p processors: µ decreases linearly with p, since µ = µ ind p , where µ ind is the MTBF of each individual component (see Proposition 1.2 in [START_REF] Hérault | Fault-Tolerance Techniques for High-Performance Computing, Computer Communications and Networks[END_REF]). For instance, with µ ind = 10 years and p = 10 5 , we have µ ≈ 50 minutes, and it goes down to a failure every 5 minutes for p = 10 6 .

The classical technique to deal with failures is to use a checkpoint-restart mechanism: the state of the application is periodically checkpointed on stable storage, and when a failure occurs, we can recover from the last valid checkpoint and resume the execution, rather than starting again from scratch. Checkpointing policies have been widely studied, see [START_REF] Hérault | Fault-Tolerance Techniques for High-Performance Computing, Computer Communications and Networks[END_REF] for a survey of various protocols and the derivation of the Young/Daly formula [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF][START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF] for the optimal checkpointing period. Recent advances include multi-level approaches, or the use of SSD or NVRAM as secondary storage [START_REF] Cappello | Toward exascale resilience: 2014 update[END_REF].

Another technique that has been advocated to deal with failures is process replication, where each process in a parallel MPI (Message Passing Interface) application is duplicated to increase the mean-time to interruption. More precisely, each processor of the platform is paired with a replica so that the execution can continue whenever one of them is struck by a failure. Given the high rate of failures expected in current systems, process replication is usually combined with a periodic checkpointing mechanism, as proposed in [START_REF] Schroeder | Understanding Failures in Petascale Computers[END_REF][START_REF] Zheng | Reliability-aware scalability models for high performance computing[END_REF][START_REF] Engelmann | The case for modular redundancy in largescale highh performance computing systems[END_REF] for HPC platforms, and in [START_REF] Leblanc | Volpexmpi: An mpi library for execution of parallel applications on volatile nodes[END_REF][START_REF] Yi | Using Replication and Checkpointing for Reliable Task Management in Computational Grids[END_REF] for grid computing. These approaches use process replication: each processor of the platform is paired with a replica so that the execution can continue whenever one is struck by a failure.

Another approach introduced in [START_REF] Casanova | Using group replication for resilience on exascale systems[END_REF] is group replication, a technique that can be used whenever process replication is not available. Group replication is agnostic to the parallel programming model, and thus views the application as an unmodified black box. Group replication consists in executing multiple application instances concurrently. For example, two distinct p-process application instances could be executed on a 2p-processor platform. Once an instance saves a checkpoint, the other instance can use this checkpoint immediately to "jump ahead" in its execution. Hence, group replication is more efficient than the mere independent execution of two instances: each time one instance successfully completes a given "chunk of work", the other instance immediately benefits from this success.

In this work, we extend group replication to the case of two different computing platforms executing concurrently and cooperating to the success of a given application. To the best of our knowledge, this scenario has not been explored yet. The two platforms share a set of remote disks, used as a stable storage for checkpointing. Typically, these platforms would be clusters, which may have different number of processors, and hence different MTBFs and execution speeds. Our goal is to determine the best way to have both platforms cooperate so that the execution time of the application is minimized. We design and analyze two strategies: 1. A periodic checkpointing strategy, where both platforms checkpoint periodically once they have executed a chunk of work of size W . Both platforms synchronize through the shared storage as soon as one of them has completed the execution of a chunk (at the time of the checkpoint). We provide a thorough analysis to express the overhead given the checkpointing period W , and we derive the size of the optimal pattern.

2. An on-failure checkpointing strategy, where each platform progresses at its own speed, and checkpoints only when a failure occurs on the other platform. Hence, when a failure occurs on one of the platforms (say platform A), the other one (platform B) checkpoints, and platform A gets a copy of this checkpoint to restart its execution at this point. Intuitively, if both platforms have the same speed, we will never roll back with this strategy, unless a failure occurs during checkpoint. We compare both strategies through extensive simulations, and show the gain (opr the absence thereof) compared to using a single platform. We also assess the accuracy of the model and of our first or second-order approximations.

The rest of the paper is organized as follows. We introduce the execution model in Section 2, and derive the optimal pattern for the periodic checkpointing strategy in Section 3. The analysis for the checkpoint-on-failure strategy is given in Section 4. Section 5 is devoted to the experimental evaluation. Finally, we provide concluding remarks and directions for future work in Section 6.

Model

We consider a black-box application and replicate its execution on two different computing platforms P 1 and P 2 . The platforms may well be heterogeneous, with different processor numbers, different MTBF values and different execution speeds. Both platforms use the same stable storage system. A typical instance is the case of two clusters that share a set of storage disks. We assume that both executions can synchronize through checkpointing. Checkpoint time is C on either platform, and this includes the time to update the state of the application on the other platform. We make no further hypothesis: The checkpointing protocol can be single-level or multi-level, and the update of the application state on the other platform can take place either through the network or via the file system. Recovery time is R, independently of which platform has taken the last checkpoint.

We partition the execution of the application into periodic patterns, i.e., computational units that repeat over time. Each pattern includes W units of work (we also say a chunk of size W ) and ends with a checkpoint. With a single platform, the optimal pattern length is wellknown and obeys the Young/Daly formula [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF][START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF]. With two platforms executing concurrently, both platforms execute the pattern concurrently, and repeat until success. Once a platform succeeds, the other one stops executing and synchronizes on checkpoint. Computing the optimal pattern length turns out a challenging problem in this case.

We assume that failures independently strike the platforms with an Exponential distribution. Platform P 1 has failure rate λ 1 , which means its MTBF (Mean Time Between Failures) is µ 1 = 1 λ 1 . Similarly, P 2 has failure rate λ 2 , and MTBF µ 2 = 1 λ 2 . We let σ 1 be the execution speed of the application on platform P 1 , and σ 2 be the speed on P 2 . We assume that P 1 is the fast platform, so that σ 1 ≥ σ 2 .

The expected execution time of the pattern is E(P ): we have to take expectations, as the computation time is not deterministic. Letting T 1 = W σ 1 , we note that E(P ) > T 1 + C, the failure-free time on the fast platform. An optimal pattern is defined as the one minimizing the ratio E(P )

T 1 , or equivalently the ratio H(P ) = E(P ) T 1 -1. This latter ratio is the relative overhead paid for executing the pattern. The smaller this overhead, the faster the progress of the execution. For the theoretical analysis, we assume that checkpoint and recovery are failure-free, because this assumption does not modify the dominant terms of the overhead (see Section 3.3 for details), but for the simulations, we do account for failures striking anytime. Finally, to be able to write Taylor expansions, we also let λ be the global failure rate and write λ 1 = α 1 λ and λ 2 = α 2 λ, with α 1 + α 2 = 1.

Optimal pattern

In this section, we show how to derive the optimal pattern length. The derivation is quite technical, and the reader may want to skip the proofs.

Expected execution time

Consider a pattern P of size W , and let E(P ) denote the expected execution time of the pattern. Because we assume that checkpoints are failure-free, we have E(P ) = E(W ) + C, where E(W ) is the expected time to execute a chunk of size W .

We start with some background on well-known results on E(W ) with a single platform P 1 , before moving on to our problem with two platforms. With a single platform P 1 , let T 1 = W σ 1 and p 1 = 1 -e -λ 1 T 1 be the probability of a failure on P 1 while attempting to execute the chunk of size W . We can write

E(W ) = (1 -p 1 )T 1 + p 1 (E lost + R + E(W )).
The first term corresponds to a successful execution, while the second term accounts for a failure striking during execution, with expected time lost E lost , recovery time R and calling E(W ) recursively to restart from scratch. We know from [START_REF] Hérault | Fault-Tolerance Techniques for High-Performance Computing, Computer Communications and Networks[END_REF] that

E lost = 1 λ 1 -T 1 e λ 1 T 1 -1
, and after simplification we get E(W ) = ( 1 λ 1 + R)(e λ 1 T 1 -1) (see [START_REF] Hérault | Fault-Tolerance Techniques for High-Performance Computing, Computer Communications and Networks[END_REF]for details). We aim at minimizing the pattern overhead H(P ) = E(P )

T 1 -1. To get a first-order approximation, we assume that λ 1 W is small so that we can expand p 1 = 1 -e λ 1 T 1 into

p 1 = λ 1 W σ 1 + 1 2 λ 1 W σ 1 2 + o λ 1 W σ 1 2 .
We then derive that

H(P ) = Cσ 1 W + λ 1 W 2σ 1 + o( √ λ 1 ). The first two-terms show that W opt = Θ(λ -1/2 1
) and we retrieve the Young/Daly formula W opt = σ 1 2C λ 1 . For the optimal pattern, we have

H opt = √ 2Cλ 1 + o( √ λ 1
). Equipped with these results for a single platform, we can now tackle the problem with two platforms. We will need a second-order approximation of the form

H(P ) = Cσ 1 W + β λ W σ 1 + γ λ W σ 1 2 + δλ + o (λW ) 2 ,
where λ = λ 1 + λ 2 is the total failure rate, and β, γ and δ are constants that we derive below. With a single platform, we had β = 1 2 . With two platforms, we obtain a complicated expression for β, whose value will always be nonnegative. If β is strictly positive and above a reasonable threshold, we will proceed as above and be satisfied with the first-order approximation that gives W opt = σ 1 C βλ = Θ(λ -1/2 ). However, if β is zero or close to zero, we will need to resort to the second-order expansion to derive an accurate approximation of W opt . In particular, when P 1 and P 2 are same-speed platforms, we will find that β = 0, γ > 0 and W opt = Θ(λ -2/3 ).

As above, let E(W ) be the expected time to execute a chunk of size W with both platforms. Let T 1 = W σ 1 and p 1 = 1 -e -λ 1 T 1 as before. We write:

E(W ) = ∞ i=0 p i 1 (1 -p 1 )E i ,
where E i denotes the expected time to execute W successfully, knowing that there were i failures on P 1 before P 1 executes the chunk successfully. We point out that we condition E i on the number of failures on P 1 , independently on what is happening on P 2 . In other words, we let P 1 execute until success, but we do account for the fact that P 2 may have completed before P 1 when computing E i . Similarly, letting T 2 = W σ 2 and p 2 = 1 -e -λ 2 T 2 be the probability of a failure on P 2 , we write

E i = ∞ j=0 p j 2 (1 -p 2 )E i,j ,
where E i,j denotes the expected execution time of the pattern, knowing there were i failures on P 1 and j failures on P 2 before both platforms execute successfully.

Theorem 1. The expected execution time of a pattern E(P ), whose execution is replicated on two platforms P 1 and P 2 , is

E(P ) = (1 -p 1 )T 1 + p 1 (1 -p 1 -p 2 )E 1,0 + p 1 p 2 E 1,1 + p 2 1 E 2,0 + C + o(λ 2 W 3 )
, where p 1 and p 2 denote the probability of having a failure during the execution of the pattern on P 1 or P 2 , respectively.

Proof. We aim at deriving a second-order approximation of E. When i = 0, there is no failure on P 1 , and the execution takes time T 1 , regardless of the number of failures on P 2 . Therefore:

E(W ) = (1 -p 1 )T 1 + p 1 (1 -p 1 )E 1 + p 2 1 (1 -p 1 )E 2 + p 3 1 E
3+ where E 3+ denote the expected time to execute the pattern successfully, knowing that there were at least 3 failures on P 1 . The rationale to introduce E 3+ is that it represents a lowerorder term that can be neglected. Indeed, we have E 3+ ≤ 3W σ 1 + 3R + E. To see that: for each failure, we lose at most T 1 + R; in the worst case, the three failures strike right before the checkpoint on P 1 , and for each of them we lose the entire pattern time T 1 plus the recovery R (and P 2 has not completed execution yet). Then we re-execute the application once more, which is accounted for in re-introducing E. Similarly, we can write E 1 and E 2 as follows:

E 1 = (1 -p 2 )E 1,0 + p 2 (1 -p 2 )E 1,1 + p 2 2 E 1,2+ E 2 = (1 -p 2 )E 2,0 + p 2 E 2,1+ , where E 1,2+ ≤ W σ 1 + R + E 1 and E 2,1+ ≤ 2W σ 1 + 2R + E 2 .
Then, we use Taylor series to approximate p 1 and p 2 to α 1 λ σ 1 W + O(λ 2 W 2 ) and α 2 λ σ 2 W + o(λ 2 W 2 ), respectively. Solving for E 1 and E 2 , we derive that:

E 1 ≤ 1 1 -p 2 2 (1 -p 2 )E 1,0 + p 2 (1 -p 2 )E 1,1 + p 2 2 1W σ 1 + 1R E 2 ≤ 1 1 -p 2 (1 -p 2 )E 2,0 + p 2 2W σ 1 + 2R . Time Time C 1 T1 C 1 C 1 T2 C 1 (P1) (P2)
Figure 1: I 0 -no failure on P 1 (there can be failures on P 2 ); P 1 always finishes first.

Note that

1 1-p 2 2 = 1 + O(λ 2 W 2 ) and that 1 1-p 2 = 1 + O(λW ). Altogether, p 2 2 W σ 1 + R = O(λ 2 W 3 ) and p 2 2W σ 1 + 2R = O(λW 2
). Therefore, we derive that:

E 1 = (1 -p 2 )E 1,0 + p 2 E 1,1 + O(λ 2 W 3 ) E 2 = E 2,0 + O(λW 2 ) .
Then, putting E 1 , E 2 , and E 3+ back into E and solving for E, we obtain:

E(W ) ≤ 1 1 -p 3 (1 -p 1 )T 1 + p 1 (1 -p 1 ) (1 -p 2 )E 1,0 + p 2 E 1,1 + O(λ 2 W 3 ) + p 2 1 (1 -p 1 ) E 2,0 + O(λW 2 ) + p 3 1 3W σ 1 + 3R .
Finally, note that 1 1-p 3 = 1 + O(λ 3 W 3 ) and p 3 1 3W

σ 1 + 3R = O(λ 3 W 4 )
. Therefore, keeping second-order terms only, we obtain

E(W ) = (1 -p 1 )T 1 (I 0 ) + p 1 (1 -p 1 -p 2 )E 1,0 (I 1 ) + p 1 p 2 E 1,1 (I 2 ) + p 2 1 E 2,0 (I 3 ) + O(λ 3 W 4 ) ,
where I 0 , I 1 , I 2 and I 3 denote the four possible outcomes of the execution (up to two failures), with their associated probability. Finally, plugging E(W ) back into E(P ) = E(W ) + C, we retrieve the equation of Theorem 1.

Computing I 0 (Figure 1). Let I 0 denote the expected execution time associated with having no failures on P 1 . With probability (1 -p 1 ), P 1 finishes faster than P 2 in T 1 time, and we can write:

I 0 = (1 -p 1 )T 1 .
Using Taylor expansions to approximate

p 1 to λ 1 T 1 + λ 2 1 T 2 1 2 + o(λ 2 T 2 1
), we can write: Computing I 1 (Figure 2). Let I 1 denote the expected execution time when having exactly one failure on P 1 . Letting X ∼ exp(λ 1 ) denote the failure inter-arrival time, we have:

I 0 = 1 -λ 1 T 1 - λ 2 1 T 2 1 2 + o(λ 2 T 2 1 ) T 1 = T 1 -λ 1 T 2 1 - λ 2 T 3 1 2 + o(λ 2 T 3 1 ) . RR n°9055 Time Time C 1 t1 R T1 C 1 C 1 T2 C 1 (P1) (P2) (a) Time Time C 1 t1 R T1 C 1 C 1 T2 C 1 (P1) (P2) (b)
I 1 = p 1 (1 -p 1 -p 2 ) ∞ 0 P(X = t|X ≤ T 1 ) min(t + R + T 1 , T 2 )dt = p 1 (1 -p 1 -p 2 ) 1 P(X ≤ T 1 ) T 1 0 P(X = t) min(t + R + T 1 , T 2 )dt .
By definition, P(X ≤ T 1 ) = p 1 and P(X = t) = λ 1 e -λ 1 t , therefore:

I 1 = (1 -p 1 -p 2 ) T 1 0 λ 1 e -λ 1 t min(t + R + T 1 , T 2 )dt .
Note that min(t + R + T 1 , T 2 ) is in order of O(W ). Using Taylor series to approximate p 1 to

λ 1 T 1 + o(λW ), p 2 to λ 2 T 2 + o(λW ), e -λ 1 t to 1 -λ 1 t + o(λt)
and keeping second-order terms only, we can get:

I 1 = λ 1 (1 -λ 1 T 1 -λ 2 T 2 ) T 1 0 (1 -λ 1 T 1 ) min(t + R + T 1 , T 2 )dt + o(λ 2 W 3 ) .
The minimum depends on which platform finishes first. We know that t

+ R + T 1 ≤ T 2 ⇐⇒ t ≤ T 2 -T 1 -R
, so that we break the integral into two parts to address both cases, as follows:

I 1 = λ 1 (1 -λ 1 T 1 -λ 2 T 2 ) T 2 -T 1 -R 0 (1 -λ 1 t)(t + R + T 1 )dt + T 1 T 2 -T 1 -R (1 -λ 1 t)T 2 dt + o(λ 2 W 3 ) ,
where T 2 -T 1 -R must be both positive and less that T 1 . Finally, let r 1 = max(min(T 2 -T 1 -R, T 1 ), 0), and we can write:

I 1 = λ 1 (1 -λ 1 T 1 -λ 2 T 2 ) r 1 0 (1 -λ 1 t)(t + R + T 1 )dt + T 1 r 1 (1 -λ 1 t)T 2 dt + o(λ 2 W 3 ) . RR n°9055 Time Time C 1 t1 R T1 C 1 C 1 t2 R T2 C 1 (P1) (P2) (a) Time Time C 1 t1 R T1 C 1 C 1 t2 R T2 C 1 (P1) (P2) (b)
Figure 3: I 2 -there is one failure on P 1 and one failure on P 2 ; depending on the failure arrival times t 1 and t 2 , P 1 finishes either first (a) or last (b).

Finally, note that I 1 depends on the value of r 1 as follows:

r 1 =    T 2 -T 1 -R, if 0 ≤ T 2 -T 1 -R ≤ T 1 T 1 , if T 2 -T 1 -R > T 1 0, otherwise.
Assuming R is small in front T 1 and T 2 , we derive:

r 1 = T 2 -T 1 -R, if 1 ≤ σ 1 σ 2 ≤ 2 T 1 , if 2 < σ 1 σ 2 .
Computing I 2 (Figure 3). Let I 2 denote the expected execution time when having one failure in P 1 and one failure in P 2 . Let X 1 ∼ exp(λ 1 ) and X 2 ∼ exp(λ 2 ) denote the failure arrival time in P 1 and P 2 , respectively. We can write:

I 2 = p 1 p 2 ∞ 0 P(X 1 = t 1 |X 1 ≤ T 1 )P(X 2 = t 2 |X 2 ≤ T 2 ) min(t 1 + R + T 1 , t 2 + R + T 2 )dt 1 dt 2 = p 1 p 2 P(X 1 ≤ T 1 )P(X 2 ≤ T 2 ) T 1 0 T 2 0 P(X 1 = t 1 )P(X 2 = t 2 ) min(t 1 + R + T 1 , t 2 + R + T 2 )dt 1 dt 2 .
Again, we have P(X 1 ≤ T 1 ) = p 1 and P(X 2 ≤ T 2 ) = p 2 , as well as P(X 1 = t 1 ) = λ 1 e -λ 1 t 1 and P(X 2 = t 2 ) = λ 2 e -λ 2 t 2 . Therefore, we can write:

I 2 = T 1 0 T 2 0 λ 1 e -λ 1 t 1 λ 2 e -λ 2 t 2 min(t 1 +R+T 1 , t 2 +R+T 2 )dt 1 dt 2 .
Using Taylor series to approximate the Exponential terms to 1+o(1) and keeping second-order terms only, we can get: As before, platform P 2 finishes faster ⇐⇒

I 2 = λ 1 λ 2 T 1 0 T 2 0 min(t 1 +R+T 1 , t 2 +R+T 2 )dt 1 dt 2 + o(λ 2 W 3 ) . RR n°9055 Time Time C 1 t1 R t2 R T1 C 1 C 1 T2 C 1 (P1) (P2) (a) Time Time C 1 t1 R t2 R T1 C 1 C 1 T2 C 1 (P1) (P2) (b)
t 2 + R + T 2 ≤ t 1 + R + T 1 ⇐⇒ t 2 ≤ t 1 + T 1 -T 2 .
Therefore, we can can break the second integral into two parts, and we get:

I 2 = λ 1 λ 2 T 1 0 e -λ 1 t 1 t 1 +T 1 -T 2 0 (t 2 + R + T 2 )dt 2 + T 2 t 1 +T 1 -T 2 (t 1 + R + T 1 )dt 2 dt 1 + o(λ 2 W 3 ) ,
where t 1 +T 1 -T 2 must be both positive and less than T 2 . We find that

t 1 +T 1 -T 2 ≤ T 2 ⇐⇒ t 1 ≤ 2T 2 -T 1
, which is always true, since t 1 is comprised between 0 and T 1 . There remains one condition, and we find that

t 1 + T 1 -T 2 ≤ 0 ⇐⇒ t 1 ≥ T 2 -T 1
, so that we can break the first integral into parts. Let r 2 = min(T 2 -T 1 , T 1 ) (note that T 2 -T 2 is always positive), and we can write:

I 2 = λ 1 λ 2 r 2 0 T 2 0 (t 1 + R + T 1 )dt 2 dt 1 + T 1 r 2 t 1 +T 1 -T 2 0 (t 2 + R + T 2 )dt 2 + T 2 t 1 +T 1 -T 2 (t 1 + R + T 1 )dt 2 dt 1 + o(λ 2 W 3 ) .
Finally, note that, similarly to I 1 , I 2 depends on the value of r 2 , which can be either:

r 2 = T 2 -T 1 , if T 1 -T 2 ≤ T 1 T 1 , otherwise.
Simplifying, we find that:

r 2 = T 2 -T 1 , if 1 ≤ σ 1 σ 2 ≤ 2 T 1 , if 2 < σ 1 σ 2 .
Computing I 3 (Figure 4). Let I 3 denote the expected execution time when having two failures on P 1 . Let X 1 ∼ exp(λ 1 ) and X 2 ∼ exp(λ 1 ) denote the failure arrival time of each failure relative to their execution. We can write:

I 3 = p 1 p 1 ∞ 0 P(X 1 = t 1 |X 1 ≤ T 1 ) ∞ 0 P(X 2 = t 2 |X 2 ≤ T 1 ) min(t 1 + t 2 + 2R + T 1 , T 2 )dt 1 dt 2 = p 1 p 1 P(X 1 ≤ T 1 )P(X 2 ≤ T 1 ) T 1 0 P(X 1 = t 1 ) T 1 0 P(X 2 = t 2 ) min(t 1 + t 2 + 2R + T 1 , T 2 )dt 1 dt 2 ,
where P(X 1 ≤ T 1 ) = P(X 2 ≤ T 1 ) = p 1 , and P(X 1 = t 1 ) = λ 1 e -λ 1 t 1 and P(X 2 = t 2 ) = λ 1 e -λ 1 t 2 . Therefore, we can get:

I 3 = T 1 0 λ 1 e -λ 1 t 1 T 1 0 λ 1 e -λ 1 t 2 min(t 1 + t 2 + 2R + T 1 , T 2 )dt 1 dt 2 = T 1 0 λ 1 e -λ 1 t 1 t 1 +R+T 1 t 1 +R λ 1 e -λ 1 t 2 min(t 2 + R + T 1 , T 2 )dt 1 dt 2 .
Using Taylor series to approximate e -λ 1 t to 1 + o(1) and keeping second-order terms only, we can get:

I 3 = λ 2 1 T 1 0 t 1 +R+T 1 t 1 +R min(t 2 + R + T 1 , T 2 )dt 1 dt 2 + o(λ 2 W 3 ) .
We find that

t 2 + R + T 1 ≤ T 2 ⇐⇒ t 2 ≤ T 2 -T 1 -R,
and we break the second integral into two parts, as follows:

I 3 = λ 2 1 T 1 0 T 2 -T 1 -R t 1 +R (t 2 + R + T 1 )dt 2 + t 1 +R+T 1 T 2 -T 1 -R T 2 dt 2 dt 1 + o(λ 2 W 3 ) ,
where T 2 -T 1 -R must be both greater than t 1 + R and less than

t 1 + R + T 1 . Again, we find that T 2 -T 1 -R ≤ t 1 + R + T 1 ⇐⇒ t 1 ≤ T 2 -2T 1 -2R
, and we can break the first integral into two parts again, as follows:

I 3 = λ 2 1 T 2 -2T 1 -2R 0 t 1 +R+T 1 t 1 +R (t 2 + R + T 1 )dt 2 dt 1 + T 1 T 2 -2T 1 -2R T 2 -T 1 -R t 1 +R (t 2 + R + T 1 )dt 2 + t 1 +R+T 1 T 2 -T 1 -R T 2 dt 2 dt 1 + o(λ 2 W 3 ) ,
where T 2 -2T 1 -2R must be both positive and less than T 1 . Let r 3,1 = max(min(T 2 -2T 1 -2R, T 1 ), 0). Then, thanks to the last step, we know that the condition

T 2 -T 1 -R ≤ t 1 +R+T 1 is always verified, but there remains t 1 + R ≤ T 2 -T 1 -R ⇐⇒ t 1 ≤ T 2 -T 1 -2R.
Therefore, we can break the second term into parts, as follows:

I 3 = λ 2 1 r 3,1 0 t 1 +R+T 1 t 1 +R (t 2 + R + T 1 )dt 2 dt 1 + T 2 -T 1 -2R r 3,1 T 2 -T 1 -R t 1 +R (t 2 + R + T 1 )dt 2 + t 1 +R+T 1 T 2 -T 1 -R T 2 dt 2 dt 1 + T 1 T 2 -T 1 -2R t 1 +R+T 1 t 1 +R T 2 dt 2 dt 1 + o(λ 2 W 3 ) ,
where T 2 -T 1 -2R must be both greater than r 3,1 , and less than T 1 . Let r 3,2 = max(min(T 2 -T 1 -2R, T 1 ), 0), and we can write:

I 3 = λ 2 1 r 3,1 0 t 1 +R+T 1 t 1 +R (t 2 + R + T 1 )dt 2 dt 1 + r 3,2 r 3,1 T 2 -T 1 -R t 1 +R (t 2 + R + T 1 )dt 2 + t 1 +R+T 1 T 2 -T 1 -R T 2 dt 2 dt 1 + T 1 r 3,2 t 1 +R+T 1 t 1 +R T 2 dt 2 dt 1 + o(λ 2 W 3 ) .
Finally, and similarly to I 1 and I 2 before, I 3 depends on the value of r 3,1 and r 3,2 , which are as follows:

r 3,1 =    T 2 -2T 1 -2R, if 0 ≤ T 2 -2T 1 -2R ≤ T 1 T 1 , if T 2 -2T 1 -2R > T 1 0, otherwise. r 3,2 =    T 2 -T 1 -2R, if 0 ≤ T 2 -T 1 -2R ≤ T 1 T 1 , if T 2 -T 1 -2R > T 1 0, otherwise.
Asymptotically, the constant R is small in front of T 1 and T 2 , which tend to infinity when λ tends to zero, and we derive:

r 3,1 =    0, if 1 ≤ σ 1 σ 2 < 2. T 2 -2T 1 -2R, if 2 ≤ σ 1 σ 2 ≤ 3 T 1 , if 3 ≤ σ 1 σ 2 r 3,2 = T 2 -T 1 -2R, if 1 ≤ σ 1 σ 2 ≤ 2 T 1 , if 2 < σ 1 σ 2 .

Expected overhead

Theorem 2. The expected overhead of a pattern H(P ), whose execution is replicated on two independent platforms P 1 and P 2 is

H(P ) = Cσ 1 W + β λ W σ 1 + γ λ W σ 1 2 + δλ + o (λW ) 2 , (1)
where λ 1 = α 1 λ and λ 2 = α 2 λ with α 1 + α 2 = 1. The values of the constants β, γ and δ are provided by the following case analysis:

Case 1: 1 ≤ σ 1 σ 2 ≤ 2. β = α 1 2 -σ 2 1 + 4σ 1 σ 2 -3σ 2 2 σ 2 2 , γ = α 2 1 2 σ 2 1 -3σ 1 σ 2 + 2σ 2 2 σ 2 2 + α 1 α 2 3 2σ 3 1 -9σ 2 1 σ 2 + 12σ 1 σ 2 2 -4σ 3 2 σ 3 2 , δ = R σ 1 -σ 2 σ 2 .
Case 2: 2 ≤ σ 1 σ 2 < 3.

β = α 1 2 γ = α 2 1 6 σ 3 1 -9σ 2 1 σ 2 + 27σ 1 σ 2 2 -26σ 3 2 σ 3 2 δ = α 1 R . Case 3: 3 ≤ σ 1 σ 2 . β = α 1 2 γ = α 2 1 δ = α 1 R .
The optimal checkpointing period W opt can be obtained by solving the following third-degree equation numerically:

∂H(P ) ∂W = - Cσ 1 W 2 + β λ σ 1 + 2γ λW σ 2 1 = 0 . (2)
Proof. Let H(P ) = E(P ) T 1 -1. We can write:

H(P ) = σ 1 C W + σ 1 I 0 + I 1 + I 2 + I 3 W -1 + o(λ 2 W 2 ) .
Then, computing I 0 , I 1 , I 2 and I 3 according to the values of r 1 , r 2 , r 3,1 and r 3,2 presented in Section 3.1, we need to consider three cases, depending upon the ratio σ 1 σ 2 . The values of β, γ and δ for each of the three cases are reported above. Finally, in order to get the corresponding optimal period, we need to solve ∂H(P ) ∂W = 0, which amounts to solve Equation (2). This can be done numerically.

For cases 2 and 3 (where σ 1 ≥ 2σ 2 ), we have β = α 1 2 . If we use the first-order approximation, we neglect the last two terms with γ and δ in H(P ). Then we obtain W opt = σ 1 C βλ , a similar formula as with a single platform. We experimentally check the accuracy of the first-order approximation in Section 5.

On the contrary for case 1 (where

σ 1 ≥ 2σ 2 ), we have β = α 1 2 ( σ 1 σ 2 -1)(3 -σ 1 σ 2 ) ≥ 0 but β = 0 ⇐⇒ σ 1 = σ 2 .
We can still use the first-order approximation when β is not too close to 0. For same-speed platforms, we need to use the second-order approximation: Theorem 3. For same-speed platforms (σ 2 = σ 1 ), the expected overhead is

H(P ) = Cσ 1 W + α 1 α 2 λ 2 W 2 3σ 2 1 + o(λ 2 W 2 ) . (3)
and the associated optimal checkpointing period is

W opt = σ 1 3 3C 2α 1 α 2 λ 2 . ( 4
)
Proof. With two same-speed platforms, we have σ 2 = σ 1 . This corresponds to case 1 with β = δ = 0 and γ = α 1 α 2 3 , hence we retrieve Equation (3). Then, differentiating and solving for W , we obtain Equation (4).

It is striking to note that W opt = Θ(λ -2/3 ) for same speed platforms, instead of W opt = Θ(λ -1/2 ). Finally, with two identical platforms (α 1 = α 2 = 1 2 and λ = 2λ 1 ), we obtain

W opt = σ 1 3 3C 2λ 2 1 .

Failures in checkpoints and recoveries

So far, we have assumed that failures do not strike during checkpoints and recoveries. In this section, we show how to handle failures during these operations, and that the approximations derived in the preceding section remain valid as long as the platform MTBF µ = 1/λ is large in front of the other resilience parameters.

Let E(R) and E(C) denote the expected time to perform a recovery and a checkpoint, respectively. The probability of a failure occurring during a process of length L on platform P i is given by p

L i = 1 -e -λ i L σ i .
If a failure strikes during the recovery, we lose E lost R time due to the failure, and we account for the time to try again by calling E(R) recursively. If a failure strikes during the checkpoint, we lose E lost C time due to the failure, we account for the recovery time E(R), and we then need to re-execute the entire pattern, which is accounted for by calling E(W ) and E(C), recursively. Altogether, we have:

E(R) = p R i E lost R + E(R) + (1 -p R i )R , E(C) = p C i E lost C + E(R) + E(W ) + E(C) + (1 -p C i )C .
We know from [START_REF] Hérault | Fault-Tolerance Techniques for High-Performance Computing, Computer Communications and Networks[END_REF] that

E lost L = 1 λ i -L e λ i L -1 .
Solving the above equations and simplifying, we derive that:

E(R) = e λ i R -1 λ i , (5)
E(C) = e λ i C -1 λ 1 + (e λ i C -1) (E(R) + E(W )) . (6)
Now, recall from our previous analysis in Section 3 that the optimal pattern length satisfies W opt = Θ(λ -1/2 ) and that H opt (P ) = Θ(λ 1/2 ). Hence, in an optimized pattern, we will have E(W ) ≤ E(P ) = Wopt σ i (1 + H opt (P )) = Θ(λ -1/2 ). Then, using Taylor expansions to approximate Equations 5 and 6, we can derive the following results:

E(R) = R + O(λ) , E(C) = C + O( √ λ) .
This suggests that the expected costs to perform checkpoints and recoveries are dominated by their original costs, under the assumption of a large MTBF. Intuitively, this is due to the small probability of encountering a failure during these operations. Thus, in Section 3, replacing R and C by their expected values does not affect the expected execution time of the pattern, neither in the first-order nor second-order approximations.

On-failure checkpointing

In this section, we present another strategy. Contrarily to the approach of Section 3, the work is not divided into periodic patterns. We only checkpoint when a failure strikes either platform. More precisely, when a failure f strikes one platform, we use the other platform to checkpoint the work, so that both platforms can resume their execution from this checkpoint, in a synchronized fashion. This scheme is exposed to the risk of having a second failure f striking the other platform during its checkpoint, which would cause to roll-back and reexecute from the previous checkpoint (which was taken right after the failure preceding f , which may be a long time ago). Such a risk can be neglected in most practical settings. As before, we will assume that failures do not strike during checkpoints.

Intuitively, this checkpoint-on-failure strategy is appealing, because we checkpoint a minimum number of times. And when a failure strikes the slow platform P 2 , we do not roll-back. However, when a failure strikes the fast platform P 1 , we have to roll-back to the state of P 2 . Altogether, we expect this strategy to work better when platform speeds are close. We will experimentally assess the checkpoint-on-failure strategy in Section 5.

Expected execution time

Let E(A) denote the expected time needed to execute the application successfully, and let T base = W base σ 1 denote the total execution time of the application on the fast platform P 1 , without any resilience mechanism nor failures. Here W base denotes the total amount of work of the application.

Theorem 4. The expected execution time of the application is

E(A) = T base + T base µ C + α 1 µ σ 1 -σ 2 σ 1 . ( 7
)
where µ = 1 λ is the MTBF. Proof. We first consider the case of two identical platforms, i.e. σ 1 = σ 2 and λ 1 = λ 2 = λ 2 . In this case, as soon as a failure occurs on either platform, the other one immediately checkpoints, and both platforms synchronize on this checkpoint, before resuming execution. In other words, the execution never rolls back, and no work is ever lost. Now, in order to compute the expected execution time, we need to account for the time needed to execute the entire application T base , as well as the time lost due to failures. When a failure occurs, we only need to account for the time C to checkpoint and synchronize. In addition, we can estimate the expected number of failures as T base µ in average, and we write:

E(A) = T base + T base µ C .
This is fine for two identical platforms. However, when failure rates and speeds differ, there are two cases: (i) a failure strikes the fast platform P 1 . Then platform P 2 checkpoints, but because it is slower than P 1 , P 1 needs to rollback and we lose the extra amount of work that P 1 has computed since the last failure and synchronization; (ii) a failure strikes the slow platform P 2 . Then platform P 1 checkpoints, and because it is faster, P 2 will roll-forward instead, catching up with the execution of P 1 .

Assuming failures are Exponentially distributed, and given that a failure (from either platform) strikes during the execution of the segment, the probability that the failure belongs to a particular platform is proportional to the failure rate of that platform [START_REF] Mitzenmacher | Probability and Computing: Randomized Algorithms and Probabilistic Analysis[END_REF], i.e. the probability that the failure belongs to P 1 and P 2 are λ 1 λ = α 1 and λ 2 λ = α 2 , respectively. In order to compute the expected execution time, we first need to account for T base , which is the time to execute the application once, without failures. Then, when a failure strikes, either it strikes P 2 , with probability α 2 , and we only lose the time to checkpoint C; or it strikes P 1 , with probability α 1 , and we lose the difference between the amount of work executed on P 1 and P 2 since the last synchronization. In average, the last synchronization was when the last failure occurred, that is µ time-steps ago. During that time, P 1 and P 2 have executed µσ 1 and µσ 2 units of work, respectively, and we have lost µ σ 1 -σ 2 µ due to the failure. Altogether, we can write:

E(A) = T base + T base µ C + α 1 µ σ 1 -σ 2 σ 1 .

Expected overhead

Theorem 5. The expected overhead is

H(A) = C µ + α 1 σ 1 -σ 2 σ 1 . ( 8
)
Proof. Let H(A) = E(A) T base -1. We write:

H(A) = 1 µ C + α 1 µ σ 1 -σ 2 σ 1 .
Then, simplifying, we obtain Equation (8).

Experimental evaluation

In this section, we conduct a set of simulations, whose goal is three-fold: (i) assess the accuracy of the proposed models; (ii) compare the performance of the two replication strategies in different scenarios; and (iii) evaluate the performance improvement of the approach over classical periodic checkpointing with a single platform.

Simulation setup

This section describes the parameters used for the simulations. First, we set R = C in all cases. Indeed, the recovery time and checkpointing time are equivalent to a read (recovery) and a write (checkpoint) operation, and they take approximately the same amount of time.

Then, we set the other parameters according to real behaviors on today's supercomputers.

Because the typical failure rate for the most powerful Top500 platforms [START_REF]Top500 Supercomputer Sites[END_REF] is around 1 or 2 failures per day, we choose µ 1 = 50, 000s ≈ 14h and µ 2 = 100, 000s ≈ 28h. The speeds were set using the Rmax value (maximum performance achieved when executing LINPACK) in PFlops of Top500 platforms (list of November 2016). We always set σ 1 = 17.6 (units in Petaflops, corresponding to the Titan platform), and we build four different cases aiming at having different σ 1 σ 2 ratios: σ 2 can be either 14.0 (Cori), 10.5 (K computer), 8.1 (Trinity) or 5.1 (Theta). We also have two possible configurations for the checkpointing (and recovery) time: a small checkpoint of 60 seconds and a large checkpoint of 1800 seconds. Overall, the parameters used by default for each platform are summarized in Table 1.

For each experiment, we setup the simulator with the resilience parameters λ 1 , λ 2 , C and R, and we compute the optimal pattern length W opt , which is obtained by solving Equation 1 numerically. The total amount of work in the simulation is fixed to be 1000W opt , and each simulation is repeated 1000 times. All the figures report the optimal overhead H opt as a function of some parameter. The solid lines are simulation results: green for the fastest machine alone (with Young/Daly period), blue for the periodic checkpoint strategy, red for the on-failure checkpoint strrtegy. The dashed lines are model predictions: blue for the periodic checkpoint strategy, red for the on-failure checkpoint strategy. The simulator is publicly available at http://perso.ens-lyon.fr/aurelien.cavelan/replication-ftxs.zip.

Accuracy of the models

In this section, we study the accuracy of the models and we assess the usefulness of the second-order approximation by comparing the results obtained with both first and secondorder formulas. We take the fastest machine Titan and let its speed σ 1 vary, while keeping all other parameters fixed. Hence we always have µ 1 = 50, 000s and four possible second platforms (Cori, K-computer, Trinity, Theta) whose parameters are given in Table 1.

Figure 5 presents the evolution of the overhead as a function of σ 1 varying from σ 2 to 5σ 2 , and using a checkpointing time of 60s (left), and 1800s (right). We observe that the model matches very well the results of the simulations: the maximum relative error is 5% with C = 1800s, and is within 0.2% with C = 60s. The latter result is expected: we do not account for failures during checkpoints t in the analysis, hence the approximation gets less accurate as checkpoint time increases.

For each value of σ 1 varying from σ 2 to 5σ 2 , we set β, γ and δ in Equation 1, according to the ratio σ 1 σ 2 , which shows the accuracy of the formula in all three cases. Finally, we note that the overhead increases with larger speeds σ 1 , but the expected throughput (time per unit of Regarding on-failure checkpointing, we observe that the precision of the formula quickly degrades with larger σ 1 , because it does not take into account failures that can occur during the re-execution work, which corresponds to the factor µ( σ 1 -σ 2 σ 1 ) in Equation 8. Note that this factor grows when σ 1 increases (or when σ 2 decreases), and it is not surprising to find that the overhead is always underestimated when the two speeds are quite different.

Next in Figure 6, we compare the simulated and theoretical overheads obtained with the first and second-order approximations. Note that the plot colors have a different meaning in this figure. The difference is small when using small checkpoint time (left), but when the two speeds get close and the checkpoint cost is high (right), the first-order approximation collapses and the theoretical overhead increases dramatically (H opt = 0.5). This is because the coefficient in O(λW ) tends to 0, and the first-order approximation used to get W opt is not valid anymore. However, we show that using the second-order approximation, (i.e. considering additional terms in O(λ 2 W 2 )) still yields good results (H opt = 0.128).

Comparison of the two strategies

In this section, we compare the overhead with the two strategies against that with a single platform. Coming back to Figure 5, we make two observations. First, when the ratio between σ 1 and σ 2 is large (greater than 2 with a small checkpoint C, somewhat higher when C increases), using a periodic pattern with replication is the same as using the fast platform only: the slow platform is not useful. Second, when this ratio between between σ 1 and σ 2 increases, the on-failure checkpointing strategy becomes worst than using the fast platform alone, especially with small checkpoint costs (left). This can be explained as follows: we wait for a failure on the slow platform to checkpoint the work done by the fast platform. But given the value of µ 2 , the slow platform is struck less frequently than the fast one, hence we often lose a lot of work in expectation (remember we lose µ(σ 1 -σ 2 ) units of work when a failure strikes on P 1 ). Figures 7 to 9 show the evolution of the overhead when parameters µ 1 , µ 2 and C, R vary. Overall, we observe again that the work lost when a failure occurs on P 1 is important with the on-failure checkpointing strategy, whose overhead strongly depends upon on the second platform used. For instance, the overhead for µ 1 = 10, 000s and C = 60s goes from 0.236 (using Cori) to 1.81 (using Theta), whereas the overhead of the periodic checkpointing remains small (between 0.074 and 0.125). This observation is confirmed by Figure 8, where the overhead increases when the number of faults actually decreases on the slow platform! We see the benefits of using replication when looking at Figure 7. When µ 1 becomes small (10, 000s, or 8.6 failures per day), the overhead with a single platform (green) increases a lot, while the overhead with the periodic strategy (blue) increases only a little, even when the second platform is twice slower than the first one. For instance we have an overhead of 1.36 for P 1 alone when C = 1800s, whereas we get 0.894 when using P 1 in conjunction with Trinity, i.e. a reduction of 34%. However, when the second platform gets too slow, the improvement brought by the use of P 2 is only meaningful when the checkpointing cost is large: on Figure 9, we get 15% of improvement if C ≥ 10s with Cori, if C ≥ 760s with K, if C ≥ 4460s with Trinity, and more than 5000s with Theta.

Figure 10 presents the case of same-speed platforms. In this case, for all parameter choices (C, R, µ 1 , µ 2 ), it is interesting to see that on-failure checkpointing is the best strategy, while it was less efficient than periodic checkpointing in almost all the other scenarios that we considered. This can be explained by the fact that there is no work lost at all with this strategy, except when there is a failure during a checkpoint.

Summary

We summarize simulation results as follows:

• The model is very accurate, as long as the resilience parameters remain reasonably small. • On-failure checkpointing is generally less efficient than periodic checkpointing, except when the speeds of the two platforms are equal (σ 2 = σ 1 ). • If P 2 is really too slow compared to P 1 (σ 2 < σ 1

2 ) or if the checkpointing cost is small, there is little reason to use a second platform.

• In all other cases ( σ 1 2 ≤ σ 2 < σ 1 ), the periodic checkpointing strategy reduces the overhead by 30% in average, and up to 90% in some particular cases.

Conclusion

This work has addressed group replication for a black-box application executing on two heterogeneous platforms. We designed and thoroughly analyzed two strategies, periodic checkpointing and on-failure checkpointing. For periodic checkpointing, we have been able to analytically derive the best pattern length, using either first-order or second-order approximations. These results nicely extend the Young/Daly formula. Simulations show that the model is quite accurate. As expected, when the platform speeds have different orders of magnitude, it is better to use only the fast platform. However, periodic checkpointing is useful for a wide range of speeds, and generally more efficient than on-failure checkpointing. The latter strategy is to be preferred only when the platform speeds are close.

Future work will be devoted to extending replication with heterogeneous platforms to deal with more complex applications, such as scientific workflows arranged as linear chains or fork-join graphs. Another interesting direction is to study the bi-criteria problem with energy consumption as a second metric, in addition to total execution time, in order to better assess the cost of replication.

Figure 2 :

 2 Figure 2: I 1 -there is one failure on P 1 ; depending on the failure arrival time t 1 , P 1 finishes either first (a) or last (b).

Figure 4 :

 4 Figure 4: I 3 -there are two failures on P 1 ; depending on the failure arrival times t 1 and t 2 , P 1 finishes either first (a) or last (b).

Figure 5 :Figure 6 :

 56 Figure 5: Evolution of overhead when σ 1 varies with C = R = 60s on the left and C = R = 1800s on the right.

Figure 7 :Figure 8 :

 78 Figure 7: Evolution of overhead when µ 1 varies with C = R = 60s on the left and C = R = 1800s on the right.

Figure 9 :

 9 Figure 9: Evolution of overhead when C and R vary.

Figure 10 :

 10 Figure 10: Evolution of overhead when parameters vary, using two same-speed platforms.

Table 1 :

 1 Summary of parameters used for simulations for each platform.

	Name	Titan Cori K computer Trinity Theta
	Speed (PFlops)	17.6	14.0	10.5	8.1	5.1
	MTBF (s)	50,000		100,000

RR n°9055

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399


Xinyu Wang

Julien Ah-Pine

Jerome Darmont

Jérôme Darmont

A New Test of Cluster Hypothesis Using a Scalable Similarity-Based Agglomerative Hierarchical Clustering Framework

Keywords: hypothèse de cluster1, classification ascendante hiérarchique2, efficacité3 cluster hypothesis1, agglomerative hierarchical clustering2, effectiveness3

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Introduced in [START_REF] Jardine | The use of hierarchic clustering in information retrieval[END_REF], the Cluster Hypothesis states that "the association between documents conveys information about the relevance of documents to requests", i.e., similar documents are expected to be relevant to the same query. In previous research works, this hypothesis has been extensively tested using four agglomerative hierarchical clustering (AHC) methods, they are single link, complete link, average link and Ward's method. However, in terms of retrieval effectiveness, results of these tests are not consistent concerning a certain clustering method. For example, [START_REF] Griffiths | Hierarchic agglomerative clustering methods for automatic document classification[END_REF] claims that "average link gave the best results" in his test, while [START_REF] Willett | Recent trends in hierarchic document clustering : a critical review[END_REF] concludes that "compete link is probably the most effective method", and [START_REF] Griffiths | Using interdocument similarity information in document retrieval systems[END_REF] states that "Ward's method was found to give the best overall results". Besides, computing efficiency is usually out of scope in these works, as complexity of conventional AHC is up to O(N 3 ), efficiency is an important factor to consider in practice.

In this paper, we extend our previous work [START_REF] Ah-Pine | Similarity Based Hierarchical Clustering with an Application to Text Collections[END_REF] to a new test of cluster hypothesis by applying a scalable similarity-based AHC framework. Different from our previous work, here we focus on the retrieval effectiveness and efficiency using AHC methods in the context of cluster hypothesis. To serve this objective, our experimental setting and evaluation measure vary. The contributions of this work include : (I) optimal cluster search is applied to compare the retrieval effectiveness of seven AHC methods ; (II) computing efficiency is addressed, and the impact of improving efficiency to effectiveness is discussed. To our knowledge, our work is the first test on cluster hypothesis that addresses the efficiency issue and extends the scope to seven AHC methods.

The outline of this paper is organized as follows : Section 2 introduces the state of the art. The scalable similarity-based AHC framework that we proposed in (Ah-Pine et Wang, 2016) is described in Section 3. Section 4 details our experiment and presents obtained results. Conclusion and future work is in Section 5.

State of the Art

Cluster Hypothesis Tests and Optimal Cluster Search

Research works on cluster hypothesis tests are largely inspired by the overlap test [START_REF] Jardine | The use of hierarchic clustering in information retrieval[END_REF], the nearest neighbor test [START_REF] Voorhees | The cluster hypothesis revisited[END_REF] and the density test [START_REF] El-Hamdouchi | Techniques for the measurement of clustering tendency in document retrieval systems[END_REF]. The objective of these tests is to verify whether or not the cluster hypothesis characterizes a document collection. [START_REF] Jardine | The use of hierarchic clustering in information retrieval[END_REF]) conduct the test by measuring the overlap between distributions of similarities of RR pairs and RN pairs of documents (R :relevant, N :non-relevant). A collection with a low overlap value is believed to cluster strongly for a set of queries, and thus better separated from the non-relevant documents. Conclusion of this test claims that cluster search has potential to outperform the inverted file system. Ho-wever, this test is only experimented on one dataset. Besides, as in a collection most documents are non-relevant, RN pairs dominate RR pairs, result based on the number of RN and RR pairs is biased. [START_REF] Voorhees | The cluster hypothesis revisited[END_REF] questions the overlap test and provides an alternative. This work assumes that if a document is relevant to a query, then the document's k nearest neighbors contain relevant documents as well. With k set between [0, 5] and four tested datasets, this test concludes that cluster hypothesis holds for one collection but not for another. Limitation of this test is that it is dependent on k, setting k to different values may alter the conclusion. (El-Hamdouchi et [START_REF] El-Hamdouchi | Techniques for the measurement of clustering tendency in document retrieval systems[END_REF] apply the density test to eliminate impact of external parameter on the experiment. This test is empirically demonstrated to be more correlated than the overlap test and the nearest-neighbor test, in proving that cluster-based retrieval outperforms document-based retrieval.

Apart from verifying cluster hypothesis, other tests are carried out to compare retrieval effectiveness using different clustering techniques, among which, four hierarchical clustering methods, single link, complete link, average link and Ward's method are widely experimented [START_REF] Griffiths | Hierarchic agglomerative clustering methods for automatic document classification[END_REF] [START_REF] Willett | Recent trends in hierarchic document clustering : a critical review[END_REF] [START_REF] Griffiths | Using interdocument similarity information in document retrieval systems[END_REF]. As stated in Section 1, these works result in different conclusions. [START_REF] Tombros | The effectiveness of query-specific hierarchic clustering in information retrieval[END_REF] agrees with the conclusion in [START_REF] Griffiths | Hierarchic agglomerative clustering methods for automatic document classification[END_REF] that average link is found to be the most effective, however unlike other works, it computes retrieval effectiveness with partial but not the entire dendrogam of document clusters.

Another branch of tests focus on comparison of cluster-based searching strategies using hierarchical clustering, such as bottom-up search and top-down search [START_REF] Van Rijsbergen | Document clustering : An evaluation of some experiments with the Cranfield 1400 collection[END_REF] [START_REF] Willett | Recent trends in hierarchic document clustering : a critical review[END_REF]. Both searching strategies require to compute the similarity between a query and a cluster representative, mostly cluster centroid is used as the representative of a cluster. These works suggest that bottom-up search results in higher retrieval effectiveness than top-down search. Different from the two searching strategies, optimal cluster search [START_REF] Jardine | The use of hierarchic clustering in information retrieval[END_REF] does not involve actual matching between query and cluster centroids. It searches for the optimal cluster for a given query by scanning the whole hierarchies of the dendrogram. The advantage of this searching strategy is that, it directly concerns the internal connections of hierarchies when computing the retrieval effectiveness for a query, and it eliminates any bias brought from external sources to the document hierarchy.

Agglomerative Hierarchical Clustering and Lance-Williams Formula

Compared to flat clustering, hierarchical clustering is more informative as it thoroughly reveals internal connections among clusters. This property is favorable in many text mining tasks. Hierarchical clustering family has two principle members, the agglomerative (AHC) and the divisive (DHC). DHC is a top-bottom process that begins with splitting a given dataset and results in clusters of which each is a data instance. AHC works in the opposite way. For DHC, if the dataset is large enough, finding an optimal split can be NP-hard. However, the worst complexity of AHC is O(N 3 ).

There are many works that improve the performance of AHC [START_REF] Xu | Survey of clustering algorithms[END_REF]. In this paper, our scope is the conventional AHC and its methods.

Given a dataset of N instances, the general procedure of conventional AHC is shown as in Algorithm 1, the output is a binary tree-like structure, named dendrogram. In the Vector Space Model and when TF-IDF weighting system is applied, documents can be projected into the feature space where the projection is sphere-like. In this case Euclidean distance is a suitable choice to determine the dissimilarity of two document vectors, D(x, y). In AHC, an essential step is to compute D(C ij , C k ) after cluster C i and cluster C j are merged. There are seven methods to determine D(C ij , C k ) by using either clustering centroids or graphic representations of clusters. [START_REF] Lance | A general theory of classificatory sorting strategies II. Clustering systems[END_REF] formulate the seven conventional clustering methods in a unified framework and name it the Lance-Williams (LW) formula. Equation 1and Table 1 display the formula and its parameter values. Thanks to their work, one can simply switch parameter values to compute D(C ij , C k ) using a certain method.

D(C ij , C k ) =α i D(C i , C k ) + α j D(C j , C k ) + βD(C i , C j ) + γ|D(C i , C k ) -D(C j , C k )| [1]

A Scalable Similarity-based AHC Framework

In (Ah-Pine et Wang, 2016) we establish a new expression of the LW formula by replacing distances with similarities, we use "Sim-AHC" to refer this approach. Sim-AHC is mathematically and experimentally proved to be equivalent to its counterpart. Input of Sim-AHC is the pairwise similarity matrix S of a normalized dataset, for example, S can be a cosine similarity matrix. Having all values in S between [0, 1] permits us to sparsify S by a threshold value, so that less memory is required to store S and less time is required to compute on it. This is an advantageous property of Sim-

Methods α i α j β γ Single 1/2 1/2 0 -1/2 Complete 1/2 1/2 0 1/2 Average |Ci| |Ci|+|Cj | |Cj | |Ci|+|Cj | 0 0 Mcquitty 1/2 1/2 0 0 Centroid |Ci| |Ci|+|Cj | |Cj | |Ci|+|Cj | - |Ci||Cj | (|Ci|+|Cj |) 2 0 Median 1/2 1/2 -1/4 0 Ward |Ci|+|C k | |Ci|+|Cj |+|C k | |Cj |+|C k | |Ci|+|Cj |+|C k | - |C k | |Ci|+|Cj |+|C k | 0 Table 1.
Lance-Williams formula : methods and parameter values AHC. However, this property cannot be found in the conventional AHC that uses distances, because in this setting, zero and close-to-zero values signify high similarity of two data instances, which is of the most interest in clustering, thus these values have to be stored for computation instead of being ignored. This is a major difference between Sim-AHC and the conventional AHC. More importantly, Sim-AHC is capable to return the same or better clustering results even when S is largely sparsified. This characteristic offers Sim-AHC computing efficiency of substantially reduced memory usage and running time. To this extend, Sim-AHC is superior to the conventional AHC. To our knowledge, Sim-AHC is the first framework that complements the conventional AHC methods formulated by the LW formula, with advantageous characteristics.

As we know, the squared form of the Euclidean distance can be expressed by a linear combination of inner products, D 2 Euclidean =< x, x > + < y, y > -2 < x, y >, where x and y are two data vectors and <, > represents the inner product. With this connection, we assume that the similarity between x and y is defined by the inner product of their normalized form, shown as Equation 2. As S(x, x) = S(y, y) = 1 holds, the corresponding dissimilarity D(x, y) can be expressed as in Equation 3. Sim-AHC amounts to working with -1 2 D(x, y) instead of D(x, y). Due to the length limitation of this paper, we do not detail the mathematical deduction of Sim-AHC here. Interested readers can refer to (Ah-Pine et Wang, 2016) for more details. Algorithm 2 illustrates the computing procedure of Sim-AHC.

S(x, y)

=< x ||x|| , y ||y|| > [2] D(x, y) = || x ||x|| - y ||y|| || 2 = S(x, x) + S(y, y) -2S(x, y) = 2[1 -S(x, y)] [3]
Algorithm 2: Computing procedure of Sim-AHC Data: pairwise similarity matrix S of a normalized dataset with N instances Initialize a dendrogram of N leaves with null height values; while number of clusters > 1 do 1.

(C i , C j ) = argmax[S(C x , C y ) -1 2 (S(C x , C x ) + S(C y , C y ))],
i.e. search for the maximal similarity in S and the pair of clusters (C i , C j ); 2. Merge C i and C j into C ij and add a corresponding parent node in the dendrogram with height value

[S(C i , C j ) -1 2 (S(C i , C i ) + S(C j , C j ))]; 3. Update S by computing similarity of S(C ij , C k ) and the self similarity S(C ij , C ij ) : S(C ij , C k ) = α i S(C i , C k ) + α j S(C j , C k ) + βS(C i , C j ) -γ|S(C i , C k ) -S(C j , C k )| S(C ij , C ij ) = δ i S(C i , C i ) + δ j S(C j , C j ) end Result: a dendrogram of 2N -1 nodes
Sim-AHC keeps the values of parameters α i , α j , β and γ unchanged as in Table 1. To guarantee the equivalence for each individual clustering method, the values of newly added parameters δ i and δ j are set differently : δ i = δ j = 1 4 for median method ;

δ i = |Ci| 2 (|Ci|+|Cj |) 2 , δ j = |Cj | 2
(|Ci|+|Cj |) 2 for centroid method ; and for the other five methods, values of δ i and δ j can be determined freely as long as their sum is 1.

Experiments and Results

Our experiments apply Sim-AHC to verify : (I) if retrieval effectiveness comprises when efficiency is improved by sparsifying the similarity matrix, and (II) which clustering method outperforms the others in this setting. Four classic medium-sized datasets described in Table 2 are used in our tests. These datasets have been experimented in previous tests [START_REF] Voorhees | The cluster hypothesis revisited[END_REF] [START_REF] Tombros | The effectiveness of query-specific hierarchic clustering in information retrieval[END_REF], they contain a complete query set and a relevance judgment file, which allow us to evaluate our results. In preprocessing, each dataset is converted into a document-term matrix using the TF-IDF weighting scheme. Stop words are removed and terms are stemmed by Porter Stemmer. A rough feature selection is applied by removing terms that appear in less than 0.2% and more than 95% documents in each dataset, and each document vector is normalized by L 2 -norm.

We apply E-measure [START_REF] Jardine | The use of hierarchic clustering in information retrieval[END_REF] to evaluate the retrieval effectiveness. Expressed as E = (β 2 +1)P R β 2 P +R , smaller E value indicates better retrieval effectiveness. P and R represent precision and recall, respectively. β is the parameter that balances the importance between precision and recall, it takes values of 0.5, 1.0 and 2.0. Given a clustering method, Sim-AHC takes a preprocessed document-term Table 2. Description of datasets matrix as input, computes its pairwise cosine similarity matrix S and outputs a dendrogram. In order to test the impact of improving efficiency on retrieval effectiveness, we choose five threshold values at 10th, 25th, 50th, 75th and 90th percentiles of the distribution of similarities in S. Table 3 shows the threshold values at these percentiles for tested datasets.

Dataset/ Percentile 10% 25% 50% 75% 90% MED 0.0047 0.0094 0.0188 0.0351 0.0618 CACM 0.0111 0.0215 0.0396 0.0692 0.1240 CISI 0.0100 0.0200 0.0381 0.0659 0.1024 LISA 0.0092 0.0180 0.0353 0.0626 0.0981 Table 3. Threshold values at 10th, 25th, 50th, 75th and 90th percentiles

In each dataset, for each clustering method, we apply Sim-AHC and sparsify the cosine similarity matrix by the threshold values. After a dendrogram is obtained, we flatten it and apply optimal cluster search for each query. Given a query Q, E value at each level of the flattened dendrogram is computed. The optimal cluster is the one that has the minimal E value, and this E value is recorded for Q. In order to reflect the overall retrieval effectiveness of a clustering method, optimal E values for all queries in a dataset are averaged.

Comparison of seven clustering methods. Table 4 displays the averaged optimal E values for each clustering method using three β values obtained from the full-sized dendrograms of each dataset, no filtering strategy is applied to sparsify the cosine similarity matrices. Values highlighted in bold are column-wise minimums, indicating the best retrieval effectiveness among seven clustering methods. This result partially assents to previous finding that complete link and average link outperform single link and Ward's method. However, this only holds true for CACM and CISI datasets. It is obvious that Ward's method performs the best for LISA dataset. Besides, Mcquitty outperforms the other methods in some cases. This experiment demonstrates that complete link and average link do not always dominate the other methods as reported in previous works. Impact of sparsifying similarity matrix on retrieval effectiveness. We are interested to test if cluster retrieval effectiveness comprises when the computing efficiency of Sim-AHC is improved by sparsifying the similarity matrix S. In this experiment, our baseline is the absolute running time (in seconds) and memory usage when a fullsized similarity matrix S is used as input. We then record the relative running time and relative memory usage of Sim-AHC when S is sparsified by a certain threshold value. Results of this experiment are displayed in Figure 1, where x-axis corresponds to percentiles. The dotted line with circle signs and the solid line with triangle signs respectively represent the relative memory usage and the relative running time. The dashed line with plus sign, and dashed line with cross signs and the dashed line with square signs plot the E values at β = 0.5, 1.0 and 2.0, respectively. A surprising finding is that retrieval effectiveness does not comprise when S becomes more and more sparsified, as we can observe in each sub-figure, that the optimal E values are mostly invariant to the changes of threshold values. This demonstrates that even when S is largely sparsified, when memory usage and running time are significantly reduced, retrieval effectiveness can be guaranteed. This discovery also implies that Sim-AHC is a scalable algorithm, and it is favorable in tasks where efficiency is demanded.

Conclusion and Future Work

In this paper, we conduct a new cluster hypothesis test. The innovative features of this test are that (I) seven AHC methods unified by the LW formula are experimented, (II) a scalable similarity-based AHC framework is applied, the impact of improving efficiency on retrieval effectiveness is discussed, and (III) optimal cluster search is employed, the retrieval effectiveness is addressed based on the optimal E values. Our results reveal that complete link and average link are not always the best-performing methods in terms of retrieval effectiveness. We also discover that when applying filte-ring strategy in Sim-AHC, retrieval effectiveness does not comprise when the similarity matrix becomes more and more sparse, that is, computing efficiency can be achieved when retrieval effectiveness is guaranteed. This property inspires us to conduct future tests of Sim-AHC on larger datasets. Currently, we are working on an implementation of a distributed version of Sim-AHC supported by Apache Spark engine. We anticipate that future tests using this implementation would help us reveal more interesting facts of cluster hypothesis in the scale of big data environment.

  0.5 β = 1 β = 2 β = 0.5 β = 1 β = 2 Single 0

Figure 1 .

 1 Figure 1. Restuls of retrieval effectiveness, relative memory usage and relative running time when sparsifying similarity matrix

  Algorithm 1: General procedure of AHC Data: pairwise distance matrix D of a dataset Initialize a dendrogram of N leaves with null height values; while number of clusters > 1 do 1. (C i , C j ) = argminD(C x , C y ), i.e. search for the minimal distance in D and the pair of clusters (C i , C j ); 2. Merge C i and C j into C ij and add a corresponding parent node in the dendrogram with height value D(C i , C j ); 3. Compute distances between C ij and other clusters C k , and update D accordingly.

end Result: a dendrogram of 2N -1 nodes

Table 4 .

 4 Comparison of retrieval effectiveness for seven clustering methods

		.791	0.876 0.869	0.743	0.819 0.846
	Complete	0.755	0.849 0.857	0.714	0.793 0.807
	Average	0.771	0.855 0.851	0.745	0.816 0.830
	Mcquitty	0.768	0.851 0.850	0.742	0.812 0.828
	Centroid	0.847	0.913 0.885	0.759	0.837 0.864
	Median	0.813	0.894 0.881	0.756	0.834 0.863
	Ward	0.772	0.857 0.855	0.730	0.805 0.822
	Dataset		CISI		LISA
	Method β = 0.5 β = 1 β = 2 β = 0.5 β = 1 β = 2
	Single	0.821	0.881 0.849	0.795	0.859 0.880
	Complete	0.801	0.856 0.819	0.725	0.776 0.776
	Average	0.799	0.853 0.817	0.735	0.795 0.799
	Mcquitty	0.799	0.857 0.825	0.732	0.792 0.796
	Centroid	0.901	0.922 0.873	0.871	0.916 0.934
	Median	0.857	0.897 0.857	0.842	0.891 0.906
	Ward	0.801	0.858 0.823	0.710	0.765 0.764

Acknowledgment

This work is supported by the French national project REQUEST PIA/FSN.

Bibliographie

MED

CISI CACM LISA Single q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Complete q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Average q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

Mcquitty q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Centroid q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Median q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Ward q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0


Xinyu Wang
email: xinyu.wang@univ-lyon2.fr

Julien Ah-Pine
email: julien.ah-pine@univ-lyon2.fr

Jérôme Darmont
email: jerome.darmont@univ-lyon2.fr

SHCoClust, a Scalable Similarity-based Hierarchical Co-clustering Method and its Application to Textual Collections

Keywords:

research documents, whether they are published or not. The documents may come

I. INTRODUCTION

As an important branch of unsupervised methods, clustering techniques always succeed in drawing interests of researchers in various domains. Clustering data in different manners reveals different information conveyed by the data. For example, hierarchical clustering groups data instances in an iterative fashion and outputs a binary-tree like structure called dendrogram, which explicitly displays the internal connections among sub-clusters. However, the process of generating a dendrogram is computationally costly. Differently, co-clustering is a technique of subspace clustering, instead of grouping only data instances, it performs clustering in both data space and feature space. In a resulted co-cluster, a partition of data instances is described by a partition of features.

Interested in discovering hierarchical structures inside coclusters and in attenuating the influence of high complexity of hierarchical clustering, in this paper we propose a similaritybased hierarchical co-clustering algorithm, SHCoClust. It is built upon a bipartite spectral co-clustering method [1], and a recent approach of the agglomerative hierarchical clustering [2], which allows sparsifying the proximity matrix to reduce computational cost. SHCoClust is capable to organize coclusters in a form of dendrogram, and to discover cluster hierarchies inside co-clusters. We believe this property is valuable in understanding the content of data. However, none of the methods that contributes to establish SHCoClust is able to do so. Furthermore, with inner product at base, SHCoClust can employ kernel functions, which offers it the capability to process non-linearly separable datasets more effectively. To our knowledge, SHCoClust is the first method that is able to reveal cluster hierarchies of and inside co-clusters, and to apply kernel functions in co-clustering tasks.

Through this paper, we begin with introducing the methods that we use to establish SHCoClust in Section II. After detailing our approach in Section III, we illustrate the experiments and results in Section IV. Section V concludes our current work and presents further works.

II. STATE OF THE ART

A. The Bipartite Spectral Co-clustering

Since introduced in [3], co-clustering has been widely applied in the domains of text mining and bio-informatics. Capable to simultaneously partition data instances and their features, co-clustering is favored in searching for the most representative terms for a set of closely grouped documents, or in searching for the most expressive conditions for a set of highly related genes. Details of many published works on co-clustering can be found in [4], [5], [6].

Based on bipartite spectral graph partitioning, [1] models a collection of documents as an undirected bipartite graph, in which documents and terms are vertices and TF-IDF weights are edges. Discovering co-clusters of documents and terms is to find the best cuts in this graph, eventually resulting in K subgraphs, in which vertices are highly associated and are of equally-sized. To solve this NP-hard discrete optimization problem, an effective heuristic method, the spectral graph partitioning, is applied to offer a real relaxation. The proposed solutions are the second largest left and right singular vectors of the diagonal-normalized document-term matrix. For convenience, we refer this approach as BSC (bipartite spectral co-clustering).

B. The Agglomerative Hierarchical Clustering

There are two principle techniques of hierarchical clustering, the divisive (DHC) and the agglomerative (AHC). Given a dataset D of n instances and m features, DHC begins with splitting D into two sets iteratively and results in clusters of which each is a data instance. AHC performs in the opposite way, it starts with merging two individual instances and eventually outputs a cluster of size n. They output a binary tree-like structure of 2n -1 nodes, named dendrogram. Comparing these two techniques in terms of complexity, AHC is practically preferred as it is of O(N 3 ) in the worst case. However, finding an optimal split in DHC can be NP-hard, if n is large. [7] surveys many works that improve the performance of AHC. In the scope of this paper, we focus on the conventional AHC and its methods. The computing procedure of AHC is shown in Algorithm 1, where d(, ) is a distance metric, and each data instance is initially considered as a cluster, represented by C.

Algorithm 1 Conventional AHC Require: dataset D n×m . ∀x ∈ [1, . . . , n], C x ← D x 1: distance matrix D ← ∅ 2: for i = 2 to n do 3: for j = 1 to i -1 do 4: D ij ← d(D i , D j ) 5:
end for 6: end for 7: while num clusters> 1 do 8:

(i, j) ← arg minD x,y 9:

C ij ← C i ∪ C j 10: for all C k ∈ D such that C k = C i , C k = C j do 11: D update ← ----d(C ij , C k ) 12:
end for 13: end while

In the setting of Vector Space Model and the TF-IDF weighting scheme, documents can have a sphere-like projection in the feature space, in such case the Euclidean distance is a suitable choice to compute the distance of two document vectors. In AHC, an essential step, shown as line 11 in Algorithm 1, is to determine d(C ij , C k ). There are seven methods to do so by using either the centroids or the graphic representatives of clusters. The Lance-Williams (LW) formula is a framework that unifies the seven conventional clustering methods [8]. Equation 1 and Table I

(C ij , C k ) =α i d(C i , C k ) + α j d(C j , C k ) + βd(C i , C j ) +γ|d(C i , C k ) -d(C j , C k )|.
(1)

i α j β γ Single 1/2 1/2 0 -1/2 Complete 1/2 1/2 0 1/2 Average |C i | |C i |+|C j | |C j | |C i |+|C j | 0 0 Weighted 1/2 1/2 0 0 Centroid |C i | |C i |+|C j | |C j | |C i |+|C j | - |C i ||C j | (|C i |+|C j |) 2 0 Median 1/2 1/2 -1/4 0 Ward |C i |+|C k | |C i |+|C j |+|C k | |C j |+|C k | |C i |+|C j |+|C k | - |C k | |C i |+|C j |+|C k | 0 C. Sim AHC, the Similarity-based AHC A new expression of the LW formula is proposed in [2],
where distances are replaced by similarities, we refer this approach as Sim AHC. It is mathematically and experimentally proved to be equivalent to the original LW formula. Input to Sim AHC is the pairwise similarity matrix S of a normalized dataset (that each row has norm 1). As ∀s ∈ S, 1 ≥ s ≥ 0, S can be sparsified by a threshold value, so that less memory is demanded to store S and less time is required to compute on it. This is an advantageous property in terms of efficiency. However, this property cannot be found in the conventional AHC. Because when distances are used, zero and close-to-zero values signify high similarities, which are of the most interest in clustering, thus they have to be stored for computation instead of being ignored. The similarity updating process between a newly-merged cluster C ij and any other cluster C k in Sim AHC consists of two steps, updating the cross-cluster similarity and updating the cluster's self similarity:

s(C ij , C k ) = α i s(C i , C k ) + α j s(C j , C k ) + βs(C i , C j ) -γ|s(C i , C k ) -s(C j , C k )|. (2) s(C ij , C ij ) = δ i s(C i , C i ) + δ j s(C j , C j ). (3)
Values of parameters α i , α j , β and γ in Sim AHC are kept the same as in Table I. To guarantee the equivalence for each individual clustering method, the values of newly added parameters δ i and δ j are set accordingly: for median method, δ i = δ j = 1 4 ; for centroid method,

δ i = |Ci| 2 (|Ci|+|Cj |) 2 , δ j = |Cj | 2
(|Ci|+|Cj |) 2 ; and for the other five methods, values of δ i and δ j can be determined freely as long as their sum is 1.

III. SHCOCLUST

Co-clusters returned by BSC are of the same level, like clusters in K-means, they are flat. A dendrogram returned by AHC or by Sim AHC displays how data instances are merged step by step, nevertheless, the connections among data features are unknown. Discovering hierarchies of coclusters, and exploring connections of sub-clusters inside coclusters surely provide us a better understanding of our data.

For such purpose, we propose SHCoClust, a similarity-based hierarchical co-clustering method, detailed in Algorithm 2.

Algorithm 2 SHCoClust

Require: dataset D n×m 1: R n×n (i, i) ← m j=1 D ij , C m×m (j, j) ← n i D ij 2: D norm ← R -1 2 DC -1 2 3: k ← ⌈log 2 K⌉ + 1 4: U n×k , S k×k , V T k×m ← SV D(D norm , k) 5: if S 11 ≥ • • • ≥ S kk then 6: U ′ n×(k-1) ← U n×[2,...,k] 7: V ′ m×(k-1) ← V m×[2,...,k] 8: else if S kk ≥ • • • ≥ S 11 then 9: U ′ n×(k-1) ← U n×[1,...,k-1]
10:

V ′ m×(k-1) ← V m×[1,...,k-1] 11: end if 12: Z (n+m)×(k-1) ← R n×n U ′ n×(k-1) C m×m V ′ n×(k-1) 13: similarity matrix S ← ∅, ∀x ∈ [1, . . . , n + m], C x ← Z x 14: Z ′ ← norm(Z) such k-1 c=1 Z ′2 rc = 1, ∀r ∈ [1, . . . , n + m] 15: for i = 1 to n + m do 16: for j = 1 to i do 17: S ij ← s(Z ′ i , Z ′ j ) 18:
end for 19: end for 20: while num clusters> 1 do 21:

(i, j) ← arg maxS x,y 22: Shown as line 14 in Algorithm 2, Z is normalized so that if the similarity function s(, ) in line 17 takes a form of inner product, values in matrix S are cosine similarities. The diagonal entries of S are constantly 1. This nature allows s(, ) to be extended to a kernel function, K, such that K(x, y) = φ(x), φ(y) , where φ : I → F is a mapping from I to F, with dim(I) ≪ dim(F), and x, y ∈ D. To generalize K, we apply ∀x, y ∈ D: s(x, y) = K(x, y)/ K(x, x)K(y, y). In our experiments, linear kernel and Gaussian kernel are tested. Having ∀s ∈ S, 1 ≥ s ≥ 0, given a threshold τ ∈ (0, 1) on S, we define that s = 0 if s ≤ τ to sparsify S. Consequently, less memory is required to store S and thus less time is needed to compute the while loop of lines 20-26 in Algorithm 2.

C ij ← C i ∪ C j 23: for all C k ∈ Z ′ such that C k = C i and C k = C j do

IV. EXPERIMENTS

A. Datasets, Preprocessing and Evaluation Metrics

Three well-known corpora that have been widely tested in text clustering benchmarks are used in our experiments. These collections are Reuters-215781 [9], SMART [1] and 20NewsGroups (20NG)2 [10]. Table II details six experimented datasets that are sampled from the three collections. In the rest of the paper, theses datasets are referred by their indexes shown in column "Ind.". Based on the bag-of-words assumption and the Space Vector Model, each dataset is preprocessed into a doc-term matrix, whose rows are document vectors and columns are terms. A few steps are involved in preprocessing, in detail, stop words are removed, terms that have document frequency higher than 20% and lower than 1% are removed, and TF-IDF weighting strategy is applied. "nb.docs" and "nb.terms" in Table II are the numbers of documents and of terms after preprocessing, "nb.K" is the groundtruth number of clusters. "ARPACK" is used as solver in the function of Singular Value Decomposition that forms Z. ARI (Adjusted Rand Index) [11] and NMI (Normalized Mutual Information) [12] are used to evaluate clustering quality.

B. Comparisons of Clustering Quality

Our first experiment consists of four tests to compare clustering quality among conventional AHC, BSC, SHCoClust with and without sparsification. These tests are:

1) SHCoClust without sparsification v.s. conventional AHC 2) SHCoClust with sparsification v.s. without 3) SHCoClust without sparsification v.s. BSC 4) SHCoClust with sparsification v.s. BSC In tests that concern sparsification, threshold values τ at percentile ranks {10, 25, 50, 75, 90}% of similarities in S are used to sparsify S. In this experiment, we apply linear kernel to obtain S, which is in fact filled with cosine similarities. The results of this experiment are illustrated in Figure 1. From top to down, each row maps to one test that is marked by its index. Each column lists the corresponding results for a tested dataset. Graph in each cell is a bar chart that is highlighted by a vertical line at x=0. The y-axis in each graph is labeled by the abbreviations of seven clustering methods. Bars in these graphs indicate the "difference" of values for ARI (the red, or the dark bars in a gray-scale printed paper) and for NMI (the blue, or the light ones) using a clustering method, obtained by subtraction of two tested algorithms (A v.s. B, A -B).

In the case of sparsification, we chose the highest value of ARI or NMI obtained in each clustering method through τ s, then subtract it with the ARI or NMI value obtained by the approach in comparison.

In Figure 1, we can discovery a few interesting findings:

• Compared to conventional AHC, clustering quality is largely improved using SHCoClust for all clustering methods, except for single link and the Ward method. In datasets of S0 and S1, the increase is tremendous, ARI and NMI are raised up to 0.8 and 0.7 at maximum. to BSC for all clustering methods, except for single link. Though improvements can be found in R5, R7 and NG8, other datasets do not display any. However, when sparsification is applied in SHCoClust, some noticeable amelioration can be observed. Comparing graphs of test 3) and of test 4), bars generally exhibit a left-to-right drift, diminishing their heights in the left side of the vertical line at x=0, and growing their heights in the right side. This drift is very apparent in R7 and in NG20 datasets, where ARI and NMI bars are noticeably pulled towards the right side. For other datasets, single link seems to be the least affected. To summarize, SHCoClust significantly improves clustering quality compared to the conventional AHC for most clustering methods. Without sparsification it obtains close results to BSC. However, when sparsification is applied, SHCoClust demonstrates some obvious improvements over BSC.

C. Extension to Kernel Functions

In our second experiment, we extend the similarity function s(, ) in SHCoClust to Gaussian kernel, given by K(x, y) = exp(-γ xy 2 ) for x, y ∈ D, γ = 1/dim(I) by default3 . Three representative datasets R5, S0 and NG8 are experimented and are compared to their results obtained by linear kernel. The input kernel matrix is sparsified in the same fashion as in the previous experiment. For each threshold value, ARI, NMI, relative memory usage and relative running time are recorded. Figure 2 illustrates the results of this experiment. Indexed by the clustering methods in column and the names of datasets in row, each graph in Figure 2 contains four lines: the solid line with "+" sign denotes the NMI values, the dashed line with triangle sign is for ARI, while the solid line with "o" sign and the solid line with "x" sign represent the relative memory usage and the relative running time, respectively. The x-axis corresponds to the percentile ranks that define the threshold values τ (not shown). At each τ ∈ (0, 1), exact memory usage and running time are recorded, and are divided by those at τ = 0, i.e., when no thresholding is applied, to obtain the relative memory usage and the relative running time.

In Figure 2 we can see that, as the percentile rank increases, memory usage and running time decrease correspondingly, however, ARI and NMI tend to preserve their values at some level until the percentile rank approaches closely to 1. More precisely, before the percentile rank crosses 75%, in most cases, we can obtain ARI and NMI values as high as (or higher than) using the full-sized input. In some other cases, ARI and NMI even boost after percentile rank is over 75%, like in S0 (of linear kernel using average link). Among seven clustering methods, single link is the most peculiar, its ARI and NMI are almost invariant to the effect of sparsification. Exceptionally, in R5 of linear kernel and Gaussian kernel, ARI of single link boosts when percentile rank = 90%, at which point memory usage and running time are largely reduced. Comparing the two kernel functions, overall similar behavior can be observed regarding the four curves. In some cases, Gaussian kernel returns higher metric values than linear kernel, and vice versa. Among the three experimented datasets, results of S0 are globally better, while R5 and NG8 display limited difference in their results.

Principle conclusions we can draw from this experiment are that (1) SHCoClust is capable to guarantee the clustering quality even when its input is significantly sparsified, requiring considerably reduced memory usage and running time. τ at percentile rank of 75% is a heuristically good threshold in our experiment (90% for single link). (2) SHCoClust can be easily employed with a kernel function, though only linear and Gaussian kernels are experimented, we believe that its compatibility with kernel functions enables it to handle nonlinearly separable datasets more effectively in other tasks.

D. Discussion of Complexity and Scalability

In terms of storage and time complexity, computation with the input of a full-sized pairwise similarity matrix S of N instances and features corresponds to storage complexity of O(N 2 ) and time complexity of O(N 3 ) (lines 15-25 in Algorithm 2). Sparsifying S with τ results in M non-zero values stored for computation, M < N 2 , eventually reduces storage complexity to O(M ) and time complexity to O(N M ). The linear relationship between the storage and time complexity is demonstrated by the lines of relative memory usage and running time in Figure 2. Let τ * define a threshold value, at which clustering quality is preserved as high as possible, meanwhile S is sparsified to have M as small as possible. As stated previously, τ at percentile rank of 75% is the τ * for most cases in our experiment. From this table we can see that all ARIs* are higher than those of BSC, implying that better clustering quality is achieved by SHCoClust. Moreover, ARIs* are obtained with substantially reduced memory usage and running time. Except for R5 dataset, in which median method only gains 8% in time and 25% in memory, other datasets obtain up to 86% time gain and on average 75% memory gain. This proves that SHCoClust is good at economizing computing resources, meanwhile preserving clustering effectiveness. In other words, if given limited computing resources, SHCoClust is able to process relatively large datasets. This reflects its scalability.

E. Exploration of Hierarchical Co-clusters and Structure

We use a sampled dataset D sam from the SMART collection to illustrate the hierarchies of co-clusters and the connections of sub-clusters inside a co-cluster. D sam contains 30 documents (10 from each class) and 51 terms (after preprocessing). With the aid of (c), we are able to know that in (b) the bottomright co-cluster is in fact closer to the top-left one and farther from the middle-situated co-cluster. (d) presents the content of the left-most sub-dendrogram in (c), i.e., the bottom-right co-cluster in (b). Keeping the same structure as it is in the dendrogram, (d) allows us to know which documents and terms co-exit and how they are connected in their co-cluster.

V. CONCLUSIONS AND PERSPECTIVES

In this paper, we proposed a similarity-based hierarchical co-clustering method, SHCoClust. We illustrated its properties, examined its clustering quality, efficiency, discussed its complexity, scalability and visualized its output. We believe that it fills a blank in literature with its unique characteristics.

A future improvement that we intend to do for this method

  d

  whileThe above procedure of SHCoClust outputs a dendrogram of n + m leaves and n + m -1 internal nodes. If input D n×m is a doc-term matrix, the dendrogram grows by merging a pair of doc-doc, term-term or doc-term in each iteration. Cutting the dendrogram by a desired number of clusters K results in flattened co-clusters, each contains a sub-dendrogram of documents and terms. In Section IV-E, visualization of a sampled dataset is presented.

Fig. 1 .

 1 Fig. 1. Comparisons of clustering quality among conventional AHC, BSC, SHCoClust with and without sparsification

  (a) and (b) of Figure3exhibit the doc-term matrix of D sam before and after being processed by SHCoClust. In (b), we can clearly see that three co-clusters of different size lie along the diagonal. These co-clusters are obtained by the dendrogram cut at a height, shown in (c). Different from the dendrogram obtained by a conventional AHC method, the dendrogram of SHCoClust has negative heights, thus it looks like it grows downwards. In (c), a leave node of the dendrogram is marked by the number of members in it. The three sub-dendrograms above the cutting line (from left to right) map to the bottomright, the top-left and the middle-situated co-clusters in (b).

  display the formula and its parameter values. |C x | is the number of data instances in cluster C x . With this formula, one can simply switch parameter values to compute d(C ij , C k ) for a given method.

  Table III exhibits the conservation of computing resources measured by memory gain and time gain at τ * , at which highest ARIs are obtained (marked by *).

Distribution 1.0, the ApteMod version.

We used the same dataset as in http://qwone.com/ ∼ jason/20Newsgroups/.

This default setting is used in popular SVM packages. Besides, when γ is low, Gaussian kernel provides close-to-one values and higher similarities between pairs of points.