Mingqiang Xue
email: xuem@i2r.a-star.edu.sg

Panagiotis Karras
email: karras@business.rutgers.edu

Chedy Raïssi
email: chedy.raissi@inria.fr

Panos Kalnis
email: panos.kalnis@kaust.edu.sa

Hung Keng Pung

Delineating social network data anonymization via random edge perturbation

Keywords: H.2.7 [Database Management]: Database Administration-Security, integrity, and protection, H.2.8 [Database Management]: Database applications-Data mining, K.4.1 [Computers and Society]: Public Policy Issues-Privacy privacy, social network, random perturbation, graph utility

INTRODUCTION

Today data owners (e.g., private corporations, public organizations) can store large volumes of digital human interactions. Examples include interactions among groups of players in an online game, for file-sharing in P2P networks, and among individuals in online social networks. In most cases, the network at hand can be represented by a graph G = (V, E), where vertices stand for individuals and edges for interactions or relations. Data owners wish to publish such graphs, e.g. for social studies and marketing.

Nevertheless, publication should not reveal sensitive personal information contained in the network. A naïve graph anonymization approach is to remove identifying attributes such as names, e-mails or IP addresses, and still publish the network in a way that allows all kind of analysis of its graph properties.

Unfortunately, this naïve anonymization approach does not guarantee privacy protection. As Hay et al. [START_REF] Hay | Resisting structural re-identification in anonymized social networks[END_REF] discerned, the structural characteristics of the published graph, combined with background knowledge, can expose individuals. For instance, assume an adversary has access to graph Gn in Figure 1(b), the naïvely anonymized version of the graph in Figure 1(a) and and knows that Suzanne has 3 friends, none of whom has more than 2 friends. With this knowledge, the adversary can identify Suzanne as node 3.

Backstrom et al. [START_REF] Backstrom | Wherefore art thou r3579x?: anonymized social networks, hidden patterns, and structural steganography[END_REF] have shown that a practical way for an adversary to gain structural knowledge is to embed a known subgraph (e.g., register dummy users in an online social network and link them to victims). Figure 1(c) shows our example graph including an embedded subgraph (nodes {11, 12, 13, 14, 15}). Malicious nodes were connected to Suzanne and Robert (e.g., by tricking them to respond to a bogus friendship invitation). An adversary capable to locate the inserted subgraph efficiently can also identify Suzanne and Robert and hence infer other attributes of theirs.

Let GA = (VA, EA) be the inserted subgraph and k = |VA| the number of vertices in GA. Typically, it suffices for k to be in the order of log |Go|, as the number of possible k-node subgraphs, hence the probability that the inserted subgraph is unique, grows exponentially with k 2 . Locating GA is an intractable problem [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF], hence impractical for large graphs. Therefore, [START_REF] Backstrom | Wherefore art thou r3579x?: anonymized social networks, hidden patterns, and structural steganography[END_REF] only considers subgraphs having a connected backbone (i.e., a Hamiltonian path of their k nodes), and describes an efficient walk-based attack: an adversary can locate GA by searching for k-node paths having the same degree sequence as GA (e.g., the subgraph of nodes {11 . . . 15} in Figure 1(c) has degree sequence [START_REF] Bonchi | Identity obfuscation in graphs through the information theoretic lens[END_REF][START_REF] Da | Characterization of complex networks: A survey of measurements[END_REF][START_REF] Cheng | k-isomorphism: privacy preserving network publication against structural attacks[END_REF][START_REF] Cormode | Anonymizing bipartite graph data using safe groupings[END_REF][START_REF] Bonchi | Identity obfuscation in graphs through the information theoretic lens[END_REF]).

In this paper, we study the effectiveness of random edge perturbation as a tool for preventing structural identification in an anony- mized graph. We assume that the original graph Go is first transformed to Gn by naïve anonymization, and then Gn is altered to Gp by random edge perturbation with perturbation probability µ. Figure 1(d) shows a perturbed graph Gp, in which edges [START_REF] Cormode | Anonymizing bipartite graph data using safe groupings[END_REF][START_REF] Da | Characterization of complex networks: A survey of measurements[END_REF], [START_REF] Hay | Resisting structural re-identification in anonymized social networks[END_REF][START_REF] Hay | Anonymizing social networks[END_REF] and [START_REF] Rizvi | Maintaining data privacy in association rule mining[END_REF][START_REF] Zhou | Preserving privacy in social networks against neighborhood attacks[END_REF] in Gn have been removed, whereas edges (13, 15) and (5, 10) have been added. Fittingly, the degree sequence of GA cannot be found in Gp; therefore the walk-based attack fails.

Remarkably, Ying and Wu [START_REF] Ying | Randomizing social networks: a spectrum preserving approach[END_REF] reject random edge perturbation as a method for privacy preservation on the grounds that it obfuscates important graph characteristics; Bonchi and Tassa arrive at a similar, albeit more positive result [START_REF] Bonchi | Identity obfuscation in graphs through the information theoretic lens[END_REF]. We inspect such claims by showing that estimation algorithms can accurately recover important graph properties (e.g., density, degree distribution, transitivity, modularity, etc) from the perturbed graph, not only for the complete original graph, but more crucially, even for subgraphs thereof -for which publishing such measures directly to the end-user would not be a solution. We apply this methodology on several graph metrics used in graph analysis [START_REF] Da | Characterization of complex networks: A survey of measurements[END_REF]. Our estimation algorithms are not applicable only on the presented metrics; in addition, we introduce a generic framework for estimating a class of utility metrics.

Nevertheless, we are wary that an adversary can employ a methodology similar to ours so as to launch sophisticated attacks. We demonstrate this potential by generalizing the walk-based attack of [START_REF] Backstrom | Wherefore art thou r3579x?: anonymized social networks, hidden patterns, and structural steganography[END_REF] to a probabilistic walk-based attack, which examines multiple possible degree sequences, each with its own probability of appearance. We show that, while being more computationally intensive than the deterministic variant, our probabilistic attack is practicable in perturbed graphs with higher success probability. Furthermore, we generalize our study to take into account any possible structural attack, even by extremely powerful adversaries having the computational power to enumerate all possible subgraphs of size k inside Gn. Our analysis delineates the probability of success in that case, for subgraph size k and perturbation probability µ.

To summarize, our contributions are outlined as follows: We show how important graph properties can be accurately recovered from a perturbed graph, even though this perturbed graph does not look like the original. We introduce a novel probabilistic attack, more powerful than the one [START_REF] Backstrom | Wherefore art thou r3579x?: anonymized social networks, hidden patterns, and structural steganography[END_REF]. We compute the probability of success of any structural attack, offering a yardstick for assessing the privacy risk entailed by the publication of perturbed graph data.

RELATED WORK

The idea that privacy can be preserved after randomization while data characteristics can be partially recovered, which we espouse, was debuted in the context of association rule mining [START_REF] Evfimievski | Privacy preserving mining of association rules[END_REF][START_REF] Rizvi | Maintaining data privacy in association rule mining[END_REF].

Backstrom et al. [START_REF] Backstrom | Wherefore art thou r3579x?: anonymized social networks, hidden patterns, and structural steganography[END_REF] were the first to present an attack demonstrating that naïve anonymization is insufficient to ensure privacy in graph data. They emphasize the differences between active attacks, where the adversary adds nodes and edges before publication, and passive attacks, where they do not, and propose the walk-based and cut-based attacks; the extent to which random edge perturba-tion can protect from this family of active attacks was not examined in [START_REF] Backstrom | Wherefore art thou r3579x?: anonymized social networks, hidden patterns, and structural steganography[END_REF]. Hay et al. [START_REF] Hay | Anonymizing social networks[END_REF] proposed a privacy protection technique based on random edge insertions and deletions, similar to the random edge perturbation that we study. However, the efficiency and utility preservation of their scheme were not studied. In [START_REF] Hay | Resisting structural re-identification in anonymized social networks[END_REF], Hay et al. formalize two adversary knowledge models based on subgraph queries and hub fingerprints, and propose an approach using structure aggregation and sampling to protect from such attacks; however this method does not preserve the complete graph structure. Ying and Wu [START_REF] Ying | Randomizing social networks: a spectrum preserving approach[END_REF] assert that random edge perturbation significantly degrades a graph's utility; this assertion is true if utility metrics are measured on the perturbed graph; however [START_REF] Ying | Randomizing social networks: a spectrum preserving approach[END_REF] does not attempt to recover graph properties by estimation algorithms.

Zhou and Pei [START_REF] Zhou | Preserving privacy in social networks against neighborhood attacks[END_REF] use edge addition (along with generalization of node labels) in order to ensure that, for every node in the graph, there exist at least k-1 other nodes that share isomorphic neighborhoods. Similarly, Liu and Terzi [START_REF] Liu | Towards identity anonymization on graphs[END_REF] bring a graph in a form such that, for each node, there exist at least k -1 other nodes having the same degree, while Cormode et al. [START_REF] Cormode | Anonymizing bipartite graph data using safe groupings[END_REF] study the anonymization problem for bipartite graphs and propose the notion of safe groupings. Zou et al. [START_REF] Zou | k-automorphism: A general framework for privacy preserving network publication[END_REF] propose the model of k-automorphism, which guarantees that the anonymized graph consists of k structurally indistinguishable subgraphs; Cheng et al. [START_REF] Cheng | k-isomorphism: privacy preserving network publication against structural attacks[END_REF] expand on this concept with k-isomorphism, in which a graph is published in the form of k identical disjoint subgraphs. Unfortunately, such syntactic anonymization methods can severely alter the nature of the published network so as to conform to the privacy guarantee they provide. For example, by k-isomorphism, a connected graph is converted to an assortment of k identical disjoint graphs; thus, the published graph is not an anonymized version of the whole original, but only of 1 k thereof. Besides, the computational complexity of these syntactic schemes renders them inapplicable on very large graphs.

Recently, Bonchi et al [START_REF] Bonchi | Identity obfuscation in graphs through the information theoretic lens[END_REF] used an entropy-based privacy metric to evaluate several graph anonymization techniques; they conclude that techniques based on random perturbation are comparable to syntactic methods. Our work builds upon these findings by inspecting the potential to estimate a graph's original properties and the viability of structural attacks after random perturbation.

UTILITY MEASURE ESTIMATION

Let Go = (V, E) be an undirected graph representing the original social network and let N = |V |. Without loss of generality, we assume that V is ordered and V [i] refers to the i th node in the set. Naïve anonymization replaces all labels in Go with random pseudonyms. The resulting graph is denoted by Gn, where new labels are given by L(V [i] ,Gn), i ∈ [1, N ]. As Gn and Go contain the same set of nodes and edges, Gn = (V, E).

Random edge perturbation yields a perturbed graph Gp = (V, Ep) from a naïvely anonymized graph Gn by adding or removing edges. Let µ ∈ [0, 1 2 ) be a user-defined perturbation probability, known to the legitimate data recipient, and (V [i], V [j]) denote the edge between nodes V [i] and V [j]. For every pair of nodes

V [i], V [j] ∈ V it holds that: Pr (V [i], V [j]) ∈ Ep | (V [i], V [j]) ∈ E = µ Pr (V [i], V [j]) ∈ Ep | (V [i], V [j]) ∈ E = µ
We assume that the adversary knows a subgraph GA = (VA, EA) in the original graph Go. GA contains k ordered nodes, with labels L(VA[i], GA) known to the adversary; the latter aims to identify GA in Gp, i.e., an ordered set Y of k nodes in Gp, such that L(Y [i], Go) = L(VA[i], GA), for 1 ≤ i ≤ k. We are interested to quantify the probability that the adversary will succeed.

Before we study the adversary's success probability, we study the effect of perturbation on different graph utility metrics. No single metric for measuring a graph's utility exists, as organizations may use a graph for different purposes. Still, there are a few widely accepted metrics [START_REF] Da | Characterization of complex networks: A survey of measurements[END_REF]. We study four popular metrics: density, degree distribution, transitivity, and modularity. While past research has asserted that such metrics rapidly degrade with random edge perturbation [START_REF] Ying | Randomizing social networks: a spectrum preserving approach[END_REF], we outline how, even when the metrics measured on the perturbed graph vary greatly from those of the original, we can still accurately estimate their original values.

Density Estimation

The density metric for a general graph is the ratio of the number of edges in the graph over the number of possible edges. It describes the average level of connectivity between nodes. Formally, the density value is defined as follows:

density = 2|E| |V |(|V | -1) (1)
The density metric is widely used in social network data analysis. Admittedly, due to perturbation, the number of edges in the perturbed graph will differ from that in the original graph. However, we can use the known perturbation probability to estimate the original density. Since the density value depends on the graph size and the number of edges, when the graph size is known, we only need to estimate the original number of edges. Let h = |E| be the (unknown) number of edges in Go, and hp = |Ep| the observed number of edges in the perturbed graph Gp. We can estimate the value of h via the maximum-likelihood estimator ĥ, defined as:

ĥ = arg max h (Pr(hp|h)) (2)
where Pr(hp|h) is the probability of observing hp edges in the perturbed graph when the original graph has h edges. Intuitively, the estimator ĥ is the number of edges in the original graph that is most likely to result in hp edges in the perturbed graph.

To estimate ĥ, we need to find a value h that maximizes Pr(hp|h). Alternatively, we can find the value of ĥ by deriving an equation with the following rationale: since each edge removal and addition during perturbation can be viewed as an independent Bernoulli trial, the value of h that maximizes Pr(hp|h) is the one that renders the expected number of edges in the perturbed graph equal to hp. Given that each existing edge survives with probability (1µ) and each non-existing one is created with probability µ, the equality in consideration is expressed as follows:

hp = ĥ • (1 -µ) + M -ĥ • µ (3)
where

M = N (N -1)

2

is the number of possible edges in the graph. The quantity on the right-hand side of the equation is the expected number of edges in Gp when the number of edges in Go is ĥ. Solving Equation 3 for ĥ, we get:

ĥ = hp -M • µ 1 -2µ (4)

Degree Distribution Estimation

The degree is a cardinal feature of a node. In social networks, degrees describe the number of friends a person has. In many realworld networks, degrees exhibit a power-law distribution.

Our estimation for degree distribution is similar to that for density estimation, as a node's original degree can be estimated via its degree in the perturbed graph. While there may be a high error in the estimation for an individual node, the total degree distribution can still be accurately estimated. We first focus on estimating the original degree of a single node. Let dp be a node's observed degree in Gp, and do our estimator of its degree in Go, respectively. As in our density estimation solution, we can estimate do by setting dp to be equal to the expected degree of the node at hand in Gp when its original degree is do. Therefore it is:

dp = do • (1 -µ) + (N -1 -do) • µ (5)
Solving the latter for do, we get:

do = dp -(N -1) • µ (1 -2µ) (6)
Let do[i] be the estimator of the i th node's degree. Then [ do [START_REF] Backstrom | Wherefore art thou r3579x?: anonymized social networks, hidden patterns, and structural steganography[END_REF], do [START_REF] Bonchi | Identity obfuscation in graphs through the information theoretic lens[END_REF], . . . , do [N ]] forms an estimation of the degree sequence of all nodes in Go, whence we derive an estimate of the distribution.

Transitivity Estimation

Transitivity measures the incidence of order-3 loops in a graph. In a social network, if v1 is connected to both v2 and v3, then there is a relatively high probability that v2 and v3 are also connected (i.e., a friend's friend is a friend). Most social networks present high transitivity. Formally, a network's transitivity T is defined as:

T = 3N∆ N3 (7)
In the above formula, N∆ is the number of triangles in the network, and N3 is the number of connected triplets, i.e., subgraphs of exactly 3 connected nodes (albeit not necessary triangular).

To compute a graph's transitivity, we have to count the number of triangles and the number of connected triplets therein. Such structures may be destroyed post-perturbation, yet we can still count the number of different triplet structures in Gp, and thereby estimate the number of triangles and connected triplets in Go.

Given any triplet of nodes in Go, the possible edge connections therein fall into the following four patterns:

• Pattern 1: They form a triangle (i.e., a triplet of three edges).

We denote the number of triangles in Go as To and their observed number in Gp as Tp.

• Pattern 2: They form a connected triplet with two edges.

Let Xo (Xp) be the original (observed) number of connected triplets in Go (Gp).

• Pattern 3: They form a disconnected triplet with only one edge. Again, let Io (Ip) be the original (observed) numbers of this pattern in Go (Gp). Given a triplet pattern in Go, we can calculate the probability that it will be transformed into one of the three other patterns in Gp. For instance, a triangle Go can become a triplet of either Pattern 2, Pattern 3 or Pattern 4, remain unchanged. Figure 2 depicts the calculations of the probabilities that a pattern in Go is converted to another pattern in Gp. We emphasize that these calare obtained under the simplistic assumption that each triplet conversion is independent of others. This assumption serves our objective of deriving an estimation, which is ultimately tested in our experiments. Thus, for example, the probability that a triangle in Go remains unchanged in Gp is (1-µ) 3 ; the probability that a triangle becomes a triplet of Pattern 2 (for which three possible cases as shown in Figure 2) is 3(1µ) 2 µ. Given the counts of triplets in different patterns (Tp, Xp, Ip and Dp) in Gp, and the pattern convertion probabilities in Figure 2, we can devise a system of linear equations with four unknown estimator variables To, Xo, Îo and Do, as we did for density and degree distribution estimation. Solving this system of equations, we get:

T p X p I p D p T o (1 µ) 3 (1 µ) 2 µ( 1 µ) 2 µ( 1 µ) 2 µ( 1 µ)µ 2 (1 µ)µ 2 (1 µ)µ 2 µ 3 X o (1 µ) 2 µ( 1 µ) 3 (1 µ)µ 2 (1 µ)µ 2 µ 3 (1 µ) 2 µ( 1 µ) 2 µ( 1 µ)µ 2 I o (1 µ)µ 2 (1 µ) 2 µµ 3 (1 µ) 2 µ( 1 µ)µ 2 (1 µ)µ 2 (1 µ) 3 (1 µ) 2 µ D o µ 3 (1 µ)µ 2 (1 µ)µ 2 (1 µ)µ 2 (1 µ) 2 µ( 1 µ) 2 µ( 1 µ) 2 µ( 1 µ) 3
To = 1 (2µ-1) 3 (Dpµ 3 -Ipµ 2 + Ipµ 3 -Tp +3µTp -3µ 2 Tp + µ 3 Tp + µXp -2µ 2 Xp + µ 3 Xp) (8) Xo = 1 (2µ-1) 3 (2Ipµ -3Dpµ 2 -4Ipµ 2 +3Dpµ 3 + 3Ipµ 3 + 3µTp -6µ 2 Tp + 3µ 3 Tp -Xp + 3µXp -5µ 2 Xp + 3µ 3 Xp) (9)
Then, given the maximum-likelihood estimators of triangles To and connected triplets Xo, we estimate the transitivity of Go as:

T = 3 To Xo (10)

Modularity Estimation

Social networks exhibit community structures. A characteristic of such a community is that intra-community links are more than inter-community links. The modularity is a metric measuring whether a partition of graph exhibits community properties. To compute modularity, the network is partitioned into a fixed number of communities. A symmetric matrix A is formed such that the elements A[i, i] (i.e., the diagonal of matrix A) are the fractions of links between the nodes in the same community i. other elements A[i, j] are the fractions of links between communities i and j. The network's modularity is defined as:

M = i A[i, i] - j A[i, j] 2 (11)
The modularity value depends on the entries in matrix A. To estimate modularity, we first create an estimator  for A. To do so, we use an auxiliary symmetric matrix, B, in which B[i, i] is the number of edges within community i and B[i, j] the number of edges between community i and community j. Bp and B refer to the matrix in the perturbed graph and the estimator for B, respectively. Once B is computed,  can be derived from B by dividing each entry by the estimated total number of edges. Since the entries in B are counts of edges between nodes, we can apply similar technique as in density estimation to estimate B. Bp is directly by counting the edges within and between partitions in the perturbed graph. The relation between B[i, i] and Bp[i, i] is:

B[i, i] = 2Bp[i, i] -(z 2 i + zi) • µ 2 -4µ (12)
In the above equation zi is the number of nodes in partition i. Assuming the number of nodes in community i (j) is zi (zj), the relation between B[i, j] and Bp[i, j] is:

B[i, j] = Bp[i, j] -zi • zj • µ 1 -µ (13)
Each entry in  is then computed by dividing the corresponding entry in B the estimated number of edges in the network.

A Generic Estimation Framework

While it is impossible to provide an exhaustive enumeration of estimation methods for every possible graph utility metric, in this section we provide a generic framework for the estimation a particular class of such metrics.

A common characteristic of the above four metrics is that the utility value relies on the counts of certain substructures (a subgraph) in the graph. For example, the density value relies on the count of disconnected pairs of nodes and the count of connected pairs of nodes. Similarly, the transitivity value relies on the count of connected triplets and the count of triangles in the graph. Generally, the class of the utility metrics whose values rely on the counts of substructures can be estimated using a generic framework. The generic framework for estimation is described as follows.

Let S be a set of smax substructures relevant to a particular utility metric, where Si refers to the i th substructure. For example, in transitivity estimation, smax = 2, and the two substructures are triangles and connected triplets. Si.c refers to the count of the i th substructure in the original graph. Therefore, the utility value is a function of the counts of different structures in S, i.e. f (S1.c, S2.c, . . . , Ss max .c). To estimate the utility value, we need to estimate the value of Si.c for all substructures in S. Let Si.s be the number of nodes in the i th substructure. For transitivity estimation, the number of nodes in the two substructures (triangles and triplets) are both 3. In the perturbed graph, we count the number of all substructures involving S1.s, S2.s, . . . , Ss max .s nodes. The count of all the substructures in the perturbed graph is denoted as p1, p2, . . ., pt max . For transitivity estimation, they are the counts of the four patterns involving three nodes. The Ŝi.c ∀i, are maximum likelihood estimations for the parameters Si.c, and can be derived by solving the following maximization problem: Ŝ1.c, Ŝ2.c, . . . , Ŝsmax .c= arg max S i .c ∀i (Pr(p1, p2, . . . , pt max |Si.c ∀i))

Given the estimation Ŝi.c ∀i, we can compute the estimated utility value using function f . Note that the counting of all substructures is only feasible for substructures involving a small number of nodes. For example, in density, we count substructures involving two nodes, and in transitivity substructures involving three nodes only. For utility metrics involving substructures of a larger number of nodes, we can use sampling (defining an upper limit for counting or execution time) to estimate their values.

While not all utility metrics can be effectively recovered, there exists an ineffective algorithm (impractical due to the computation cost) that follows standard procedures to estimate all the metrics. The algorithm is described as follows: Let G denote the set of all possible graphs on N number of nodes. For each possible graph G[i], there is a probability that G[i] = Go, which can be computed based on G[i], Gp and µ, and denoted by Pr(G[i] = Go|Gp). Consider a particular metric Z, we can measure its value vali(Z) on each G[i]. Lastly, the sum of vali(Z)•Pr(G[i] = Go|Gp), ∀i forms an estimation to the metric Z in the original graph. Although theoretically sound, this algorithm is impracticable as it requires the enumeration of all possible graphs on N nodes. A more efficient recovery of specific other metrics is left as an open problem.

The quality of our estimations is also a concern. To assess that quality, we can compute the standard error of the mean, StErr. For example, from Equation 4, we estimate the variance of ĥ as:

σ 2 ( ĥ) = σ 2 (hp) (1 -2µ) 2 (14)
Let R = 1 (R = 0) denote the event there exists (does not exist) an edge between a particular pair of nodes. Therefore,

Pr(R = = h M • (1 -µ) + (1 -h M ) • µ, and Pr(R = 0) = h M • µ + (1 - h M ) • (1 -µ). Therefore, σ 2 (hp) = M • Pr(R = 1) • Pr(R = 0).
Substituting in Equation 14, we get:

σ( ĥ) = M • 1 16 • ( 1 2 -µ) 2 - h M - 1 2 2 (15)

ATTACKING THE PERTURBED GRAPH

As we argued, by carefully designing estimation algorithms, many of the original properties can be accurately estimated. Nevertheless, the same potential can also be exploited by an adversary to launch more sophisticated attacks.

Principles of the probabilistic attack

The probabilistic attack starts out from the observation that an adversary can confidently estimate a degree interval (i.e., range) for each embedded malicious node. Using an approach similar to that of the walk-based attack [START_REF] Backstrom | Wherefore art thou r3579x?: anonymized social networks, hidden patterns, and structural steganography[END_REF], the can efficiently enumerate a list of candidate degree sequences will include, with high probability, the one that represents the embedded subgraph GA.

most cases, the adversary can arrive at a single degree sequence that represents with high probability the subgraph embedded in Go. This disposition is enabled by the fact that candidate degree sequences are filtered out in our attack methodology using two tests: the interval degree check and the error-tolerant edge check. These tests work differently from those of the walk-based attack [START_REF] Backstrom | Wherefore art thou r3579x?: anonymized social networks, hidden patterns, and structural steganography[END_REF]. Table 1: Probability that k-path in GA is preserved. We face a few challenges in demonstrating the feasibility of the attack: first, the prediction of degree ranges should be correct with high probability. Second, the length of the predicted interval for an adversary's node should be as small as possible, as a large interval may result in a large number of nodes passing the degree check, causing a severe penalty to the attack's time complexity. Last, in order for our attack to succeed, the perturbation may alter the attacker's subgraph but must not have destroyed the k-path (i.e., path of k nodes) therein. In Table 1 we show the probability that this k-path is preserved (i.e., all its edges are preserved after perturbation) for several combinations of different µ and k, computed as (1µ) k-1 under perturbation probability µ. As this probability is quite high, the adversary can reasonably assume that the embedded subgraph's k-path is preserved in Gp after perturbation.

Predicting the degree interval

Let do be the degree of a malicious node in Go and dp the degree of the same node after perturbation. In the following, we first compute Pr(dp|do), i.e., the probability that, given that the node's original node degree is do, its degree after perturbation is dp.

Let r be the number of neighbors eventually removed from the set of do's neighbors in Go, and a the number of new neighbors added, due to perturbation. Without loss of generality, suppose that the dp neighbors of the malicious node at hand are generated in two steps: first, r neighbors are disconnected and the number of remaining neighbors is dor; then, a = dp -(dor) nodes are connected and become neighbors, so the total number of neighbors is dor + a = dp. Then the following three inequalities hold: Thus, an adversary can efficiently compute Pr(dp|do). The posvalues of dp after perturbation ranges from 1 to N -1. Yet the distribution of these values is not uniform. For each embedded node, the adversary can select a small subset of dp values and an interval I representing the range of possible degrees for that node. We check inclusion in this interval as our degree check.

r ≥ 0, do -r ≥ 0, dp -(do -r) ≥ 0 (16)
The removal and addition of neighbors of an embedded node can be viewed as two independent Binomial processes. The expected values for r and a are E

[r] = do • µ and E[a] = (N -do -1) • µ, respectively. Pr(dp|do) is maximized for r = E[r] and a = E[a] (hence dp = E[dp]
). Then the chosen interval I for the embedded at hand is centered at dp = E[dp], with w other values to its left and right, where w is a small non-negative integer. Eventually, the predicted degree interval for a selected malicious node in Gp is

I = [E[dp] -w, E[dp] + w].
Let Pr(dp ∈ I) be the probability that the embedded node's degree is in I after perturbation. An effective attack is possible if the adversary can find a fine-tuned value of w such that Pr(dp ∈ I) is sufficiently large and yet the width of I is small. Let I V A [i] be the degree interval for embedded node VA[i], and Dp be the degree sequence of the embedded graph GA in Gp. Then the probability that all embedded nodes' degrees fall into their respective intervals is

k i=1 Pr(Dp[i] ∈ I V A [i]
). This calculation is made under the assumption that the events of different embedded nodes' degrees falling into their respective intervals are independent of each other, which helps us derive our estimate.

Table 2 shows the values of Pr(dp ∈ I) for selected values of µ and w with k = 12. We observe that, when the number of nodes in the network is 10, 000, the perturbation probability is 0.001 and w = 4, then the probability that a single embedded node falls into the interval I is close to 1. Moreover, the probability that all the embedded nodes' degrees fall into their respective intervals after perturbation is also close to 1 in the same configuration. We conclude that, with this configuration, the attacker is almost sure that all embedded nodes will pass the interval degree check. For example, in the graph in Figure 1(c), the adversary's degree sequence is [START_REF] Bonchi | Identity obfuscation in graphs through the information theoretic lens[END_REF][START_REF] Da | Characterization of complex networks: A survey of measurements[END_REF][START_REF] Cheng | k-isomorphism: privacy preserving network publication against structural attacks[END_REF][START_REF] Cormode | Anonymizing bipartite graph data using safe groupings[END_REF][START_REF] Bonchi | Identity obfuscation in graphs through the information theoretic lens[END_REF]. Yet in the perturbed graph Gp the degree of the node labeled 11 has become 1. Then, if a walk-based attack is launched, this node will not be detected. Still, with probabilistic attack, the adversary is able to estimate the degree interval for each embedded node. For example (after integer rounding) the estimated degree intervals can be [ [START_REF] Backstrom | Wherefore art thou r3579x?: anonymized social networks, hidden patterns, and structural steganography[END_REF][START_REF] Bonchi | Identity obfuscation in graphs through the information theoretic lens[END_REF], [START_REF] Cormode | Anonymizing bipartite graph data using safe groupings[END_REF][START_REF] Da | Characterization of complex networks: A survey of measurements[END_REF], [START_REF] Cheng | k-isomorphism: privacy preserving network publication against structural attacks[END_REF][START_REF] Cormode | Anonymizing bipartite graph data using safe groupings[END_REF], [START_REF] Cheng | k-isomorphism: privacy preserving network publication against structural attacks[END_REF][START_REF] Cormode | Anonymizing bipartite graph data using safe groupings[END_REF], [START_REF] Bonchi | Identity obfuscation in graphs through the information theoretic lens[END_REF][START_REF] Cheng | k-isomorphism: privacy preserving network publication against structural attacks[END_REF] ]. In this scenario node 11 can still be accepted by the adversary as a candidate embedded node, allowing for a successful attack.

µ w Pr(dp ∈ I) k i=1 Pr(Dp[i] ∈ I V A [i] ) µ = 0.

Attack Algorithm

Algorithm 1 describes our probabilistic attack. T is the tree of all candidate subgraphs, w is the width parameter for the degree intervals used in the interval degree check (chosen by the adversary as discussed) and m is the maximum number of errors that we allow in the employed error-tolerant edge check. The main block of the algorithm is a loop that continues until a k-path is found in T by the localSearch function (shown below), or both w and m reach their predefined maximum thresholds wmax and mmax.

Our probabilistic attack algorithm follows the pattern of a walkbased attack [START_REF] Backstrom | Wherefore art thou r3579x?: anonymized social networks, hidden patterns, and structural steganography[END_REF], with the key difference being on the two tests it performs (Line 3 in localSearch). To pass the interval degree check, a node's degree should fall in the predicted degree interval VA[level]. To pass the error-tolerant edge check, the number of errors in edge checks accumulated in the path from this node to the root should not be larger than m. At each iteration of the main loop, if a k-path is not found, we relax the searching condition by either increasing w or m. However, using large w and m enlarges search space. The maximum w and m values that can be used depend on the attacker's computational resources.

Building edges to target the victims

In order to compromise the victims' privacy, the adversary has to correctly identify the victims. However, due to perturbation, the links between the victim and the adversary's nodes may have changed which raises new challenges for identifying the victims. We propose a method that minimizes the impact of perturbation and establishes robust links against perturbation. Let the set of nodes that represent the victims in the network be VT = {τ1, τ2, . . . τq}. Sτ i ⊂ VA is the set of maliciously embedded nodes that are linked to victim τi, 1 ≤ i ≤ q. Our approach is as follows: we define two parameters and ρ2, where ρ1 defines the minimum size of Sτ i , and ρ2 defines the minimum number of different members between the two sets Sτ i and Sτ j , for i = j. Formally:

|Sτ i | ≥ ρ1 ∀i ∈ [1, q] |Sτ i \Sτ j | ≥ ρ2 ∀i, j ∈ [1, q] and i = j (17)
Moreover, we require that none of the adversary's nodes share common neighbors other than the nodes in VT and VA. To prove the robustness of the links between the victims and the adversary's nodes, we show analytically that the probabilities of the three events that affect the identification of the victims are negligible.

Claim 1:

The probability that Sτ i for any τi changes due to perturbation tends to be 0 when µ → 0.

PROOF. Sτ i is preserved for any τi if the edge relations between this τi and all the adversary's nodes are preserved (i.e, (1µ) k ). Therefore, the probability that

Sτ i changes is 1 -(1 -µ) k . When µ → 0, 1 -(1 -µ) k → 0.
Claim 2: The probability that there is a node v outside sets VA and VT such that Ev, the set of edges between v and the nodes from VA, is equal to Sτ i for the victim τi, decreases fast with the increase of ρ1.

PROOF. Let us consider a particular node v which already has an edge with a node in Sτ i . The probability that it forms new edges with all other nodes in Sτ i but not with the nodes in VA -Sτ i is at most

µ ρ 1 -1 • (1 -µ) k-ρ 1 +1 . Moreover, the total number of such possible nodes v in the network is N -k -q. Therefore, (N -k -q) • µ ρ 1 -1 • (1 -µ) k-ρ 1 +1
is the probability for the event in this claim. When µ = c N , this probability at most c ρ 1 -1 N ρ 1 -1 (by taking Nkq as N and (1µ) k-ρ 1 +1 as 1), which decreases fast with the increase of ρ1.

Claim 3: The probability that the set Sτ i of malicious nodes connected to victim τi becomes the same as the set Sτ j of malicious nodes connected to victim τj after perturbation decreases fast with the increase of ρ2.

PROOF. Let the number of non-common elements in Sτ i and Sτ j be xij. Similarly to the derivation of equation 27, Sτ i is con-verted to Sτ j by perturbation if and only if an xij number of edge additions and deletions occurs. Therefore, after perturbation, Pr(Sτ

i = Sτ j ) = (1 -µ) k-x ij • µ x ij . Since ρ2 ≤ xij, Pr(Sτ i = Sτ j ) ≤ µ ρ 2 ,
which decreases fast with the increase of ρ2.

A trivial algorithm that builds our robust links operates as follows: for each victim τi, repetitively select a subset of nodes in VA at random, until a subset that satisfies both ρ1 and ρ2 requirements is found; then link τi to all the nodes in this subset of VA.

Preventing the probabilistic attack

The adversary can successfully identify his subgraph based on the assumption that the k-path is not broken. However, the publisher can increase the perturbation probability so that, with high probability, the k-path is broken and therefore the probabilistic attack is infeasible. Let ε, the secure parameter, be the maximum probability that the k-path is preserved. Therefore,

(1 -µ) k-1 ≤ ε (18)
The following inequality gives the minimum µ to be used so as to prevent the probabilistic attack with probability no less than 1ε:

µ ≥ 1 -k-1 √ ε (19)

GENERIC STRUCTURAL ATTACK

We now provide a generic analysis of a structural attack's probability of success. We assume that any structural attack can be translated into an instance of the graph isomorphism problem and show that an adversary who can enumerate all permutations of k nodes in Gp, will also be able, with high probability, to detect the embedded nodes under a random perturbation scheme for low µ.

Let Yi be the i th permutation of k nodes in the network. We assume an extreme case where the adversary has infinite computational power and is able to enumerate all permutations of k nodes in the network:

Y = {Yi : 1 ≤ i ≤ P N k } where P N k = N ! (N -k)!
is the total number of permutations. The adversary will then choose a particular permutation Y ∈ Y as a candidate for VA and will assume that Y [i] is VA[i]. However, the adversary the following two challenges when choosing Y : first, the perturbation may have changed the embedded graph GA in such a way that GA cannot be detected in Gp. Second, even if the adversary succeeds in locating a permutation Y that matches subgraph GA in appearance, there is still a probability that Y is not the actual VA, but only a look-alike.

We define λY as the probability that the chosen Y is VA, given Gp. For the sake of conciseness, we use E Y A to denote the set of edges among the nodes in Y , after perturbation, in the event that this set has been created by perturbing the set EA, i.e., in the event that Y = VA. Then λY = Pr(E Y A |Ep). Let E Y p be the set of edges among the nodes in Y in the perturbed graph Ep. Then the following theorem holds: THEOREM 5.1. For a perturbed graph Gp of size N with a perturbation value µ, the probability that an adversary detects the embedded subgraph GA, for a permutation of k nodes Y , is:

λY = Pr(E Y p |E Y A ) P N k i=1 Pr(E Y i p |E Y i A ) (20)
PROOF. First, we rewrite the expression using Bayes' theorem,

Pr(E Y A |Ep) = Pr(Ep|E Y A ) • Pr(E Y A ) P N k i=1 Pr(Ep|E Y i A ) • Pr(E Y i A ) (21)
In the above equation, Pr(E Y i A ) is the prior probability of Y being the attacker's nodes, which is the same for all i. Therefore:

Pr(E Y A |Ep) = Pr(Ep|E Y A ) P N k i=1 Pr(Ep|E Y i A ) (22)
We first focus on the numerator Pr(Ep|E Y A ) in Equation 22. We split the set of edges in the perturbed graph into two sets: E Y p , the edges between the nodes in Y , only and E Y p , the other edges in Ep. By definition, E Y p =Ep-E Y p . E Y p and E Y p are independent, hence:

Pr(Ep) = Pr(E Y p ) • Pr(E Y p ) (23)
Inserting the conditional variable E Y A to Equation 23, we get:

Pr(Ep|E Y A ) = Pr(E Y p |E Y A ) • Pr(E Y p |E Y A ) (24)
Since

E Y p and E Y A are independent, Pr(E Y p |E Y A ) = Pr(E Y p ) = Pr(Ep)
Pr(E Y p ) (by Equation 23). Then Equation 24 becomes:

Pr(Ep|E Y A ) = Pr(E Y p |E Y A ) • Pr(Ep) Pr(E Y p ) (25)
Combining Equations 22 and 25, and using E Y i p for E Y p , we get:

Pr(E Y A |Ep) = Pr(E Y p |E Y A ) P N k i=1 Pr(E Y i p |E Y i A ) (26)
The intuition behind the derived equation for λY is the following: the numerator Pr(E Y p |E Y A ) is the likelihood of the particular chosen permutation Y to be VA. The denominator is the sum of the likelihoods of each permutation Yi in the network being VA. The ratio describes the probability of success of the particular selection Y being VA. Thus, the value of λY depends on the value of Pr(E Y p |E Y A ) and the sum of Pr(E Y i p |E Y i A ) for all i. In effect, the computation of the exact λY would require the enumeration of all permutation of k nodes in the graph. In

Pr(E Y p |E Y A ) = µ m • (1 -µ) M -m (27)

Estimation of λ Y

From a (computationally powerful) adversary's point of view, Equation (20) can be used to compute λY for each Y = Yi. The attacker will then assume that the set Y gives the maximal value of λY is the embedded subgraph VA in Gp. In particular, given that the denominator is constant for a given EA and Gp, the adversary will simply opt for the set Y that maximizes the numerator in Equation ( 20). The best-case scenario for an adversary is that m = 0 (i.e., E Y p and E Y A are identical). Otherwise, an adversary will choose the Y that gives the most similar subgraph to GA.

In the following we provide a simple method for estimating the value of λY . For any permutation of k nodes Yi, let

mY i = |E Y i A ∪ E Y i p -E Y i A ∩ E Y i p |.
We consider E Y i p as a random chosen from all possible edges among the nodes of Yi. Then the expected value of mY i is M 2 and the expected value of Pr(E

Y i p |E Y i A ) is µ M 2 • (1 -µ) M 2 .
Then our simple estimate of λY is:

λY = min 1, 1 P N k 1 -µ µ M 2 -m (28)
Thus, our estimate of λY depends on µ and m. Table 3 lists different λY estimates with respect to different µ and m combinations in a network of 10, 000 nodes with 10 malicious nodes. Intuitively, m can be seen as the number of altered edges from E Y A to E Y p . We note that, unless both µ and m are high, λY approaches 1.

We observe that tends to be large only when m is small. If the lowest value of m is kept large (e.g., greater than 10, in the last row of Table 3), then λY is kept low. In effect, to assess the potential for a successful attack, we study the distribution of m under random edge perturbation. This perturbation can be seen as a binomial process that adds and deletes edges with probability µ. Then the probability distribution of m is:

Pr(m ≤ ℓ) = ℓ m=0 M m (1 -µ) M -m µ m (29)
The last column of Table 3 shows instances of the probability distribution for m with the values of k and µ fixed. We observe that Pr(m = 0) equals 0.9955 and 0.9559 for µ = 0.0001 and µ = 0.001. This result shows that m tends to be rather small with high probability. Thus, we conclude that λY tends to be large for such perturbation probabilities. This is good news for an adversary who has the computational powers to enumerate all subgraphs and choose the one most similar to the embedded one.

EXPERIMENTAL EVALUATION

now present our experimental evaluation. First, we investigate the probability of success of the probabilistic attack under different values of µ, and measure its execution time. Next, we investigate the effect of perturbation on the graph properties.

All experiments ran on a 2.33GHz CPU, Windows-XP machine with 3.25GB RAM. We employ two real datasets: The Enron dataset 1 is the graph of email exchange among employees of Enron, having 4,644 accounts. Each account corresponds to a node and two accounts are linked if they have exchanged emails in both directions. The DBLP dataset 2 is a random subset of 20,000 authors from the

Assessing the Probabilistic Attack

In our first experiment, we assess the probability of success of our probabilistic attack as opposed to that of the classical walkbased attack [START_REF] Backstrom | Wherefore art thou r3579x?: anonymized social networks, hidden patterns, and structural steganography[END_REF]. We first test the walk-based attack on the Enron and DBLP data, measuring its success rate in trials of 200 separate attack runs, as a function of the perturbation probability µ. An attack run is considered to be successful, if the adversaries can detect the sequence of embedded nodes and re-identify at least one victim node. In [START_REF] Backstrom | Wherefore art thou r3579x?: anonymized social networks, hidden patterns, and structural steganography[END_REF], the suggested number of malicious nodes is Θ(log(N )) and the number of victim nodes is q = O(log 2 (N )). Following this suggestion, we vary k at values 20, 25, 30 for both graphs, with number of victims 100, 157 and 225, respectively.

Figures 4a,e show our results, which provide a glimpse of the probability that an adversary successfully identifies the embedded nodes in perturbed DBLP and Enron data using a walk-based attack. When µ is 0, all attacks are 100% successful. Still, already for rather small values of µ (10 -7 to 10 -6 ), the success rate drops drastically to very low values. In addition, the success rate is lower for larger under the same perturbation value µ; that is because, with larger k, the node degree sequence of the malicious nodes is more likely to be changed or the backbone to be broken, making the attack more likely to fail. In effect, the walk-based attack can be effectively prevented through random edge perturbation, with minimal impact on the graph's structure (as µ negligible).

On the other hand, Figures 4b-d,f-h show the success rate of the probabilistic attack on the same DBLP and Enron data, again in trials of 200 runs each. We show results for several values of the interval-width parameter w. As the search space of the attack algorithm grows with w, the success rate also rises with it. For µ ≃ 10 -4 , the probabilistic attack succeeds in almost 100% of the cases, in stark contrast to the walk-based one. Still, as µ grows further, the observed success rate swiftly drops for all values of w. As with the walk-based attack, the success rate falls as k grows.

Figures 3a,b show the execution time of our attack as a function of interval-width w and error-tolerance m, with µ = 10 -3 . The number of malicious nodes is k = 4 log(N ) and the number of victims q = log 2 (N ). The algorithm's search space grows with both w and m (Section 4.2), hence the execution time also ascends with them, yet remains lower than 3 minutes, rendering the attack rather feasible on reasonably-sized real-world data sets.

An adversary who successfully identifies the embedded subgraph inside the perturbed graph may yet not locate the target victims, as the edges between the embedded nodes and the victims may have been removed. Table 4 shows the measured percentage of victims that can be identified in a successful attack. The number of malicious nodes and the victims remain at k = 4 log(N ) and at ➭ (10 -4 ) Density ( 10-3 ) original perturbed estimated q = log 2 (N ). As the table shows, more than 91% of victims are identified when the attack succeeds.

✵ ✵✁ ✵ ✵✁✁ ✵ ✵✁✂ ✵ ✵✁✄ ✵ ✵✁☎ ✁ ✂ ✄ ☎ ✥ ✆ ✝ ✞ ✟ ✁✵ ♠ ✠✁✵ ✲✡ ✮ ❚☛☞✌✍✎✏✎✑✎✏✒ ♦☛✎✓✎✌☞✔ ♣✕☛✏✖☛✗✕✘ ✕✍✏✎❡☞✏✕✘ ✵ ✁ ✵ ✁✁ ✵ ✁✂ ✵ ✁✄ ✵ ✁☎ ✵ ✁✆ ✁ ✂ ✄ ☎ ✆ ✥ ✝ ✞ ✟ ✁✵ ♠ ✠✁✵
µ Enron DBLP 1×10 -4 95.2% 98.3% 2×10 -4 93.6% 98.3% 3×10 -4 92.7% 96.7% 4×10 -4 91.9% 94.2%

Table 4: Percentage of affected victims

The error-tolerance m also affects an attack's probability of success. Table 5 shows an instance of this effect: in a trial of 100 attacks with m = 0 on the Enron graph perturbed with µ = 0.04, there are 53 successes, 35 failures due to false predicted interval, and 6 due to broken path or edge check failures. Still, when we relax the requirement for passing the edge check test, we can increase the number of successes to 58 and 59. To sum up, the probabilistic attack is more effective than the walk-based one and feasible in terms of runtime. Still, both can be prevented under random perturbation with sufficiently large µ.

Assessing Utility Preservation

We use the perturbed data derived in previous experiments to evaluate the extent to which graph properties are preserved. Figures 5a-f show the density, transitivity and degree distribution for perturbed DBLP and Enron data. Each figure shows the original, perturbed, and estimated values. Density and transitivity values vary with the perturbation probability µ, while the degree distribution is given as a single snapshot for µ = 10 -3 . We observe that, while graph properties deviate significantly from the originals as µ grows, our derived estimates are resilient to perturbation and the original values well. Significantly, our estimation algorithms provide accurate approximations even under perturbation probabilities that render probabilistic attacks unviable.

Furthermore, we illustrate the benefit of our method in terms of measuring the properties of subgraphs of the original graph. This application of our technique offers an undeniable benefit, as publishing the graph cannot be substituted by directly publishing the measures of interest. We measure the same six structural properties over subgraphs of different sizes extracted from the original For degree distribution, we now present the sum of absolute errors between the perturbed/estimate and original distribution values for each degree. Our results reconfirm our findings even with subgraphs extracted from the parent graphs. Overall, our results corroborate our claim that a published graph allows accurate estimations of graph properties not only for the whole graph, but also for subgraphs thereof, where publishing the measures directly to end-users would not be a satisfactory solution.

Distance-based classification

We now attempt to perform a specific data mining task, distancebased node classification, over perturbed data.

Social networks often possess hubs, i.e., nodes with very high degree. A particular person's connectivity pattern to the hubs indicates that person's interests. For example, in a social election, a person's voting pattern may indicate its political views. Thus, node classification based on such patterns is useful. We consider a classification of nodes based on the distance between their hub connectivity pattern (HCP) and some target patterns (TPs). Given a set of hubs, a node's HCP is the subset of hubs that this node has connectivity to. Each TPi is a subset of the hubs defined by the analyst. Given a set of k hubs, HCP and TPi are k-dimensional binary vectors. The distance between HPC and a particular TPi is the edit distance between the two vectors. For each TPi, a group of nodes Gi is formed by assigning group membership to the nodes that have closest distances to TPi than to all other TPj(i = j). We aim to classify each node to the right group.

We use the Wiki graph. Hubs are chosen as the nodes that have top-10 degrees in the graph (ranging from 482 to 1,053). We extract a subset of 200 nodes for classification, each having at least 4 connections to the hubs. For instance, the 10-dimensional binary data (0 0 1 0 1 0 0 1 1 0) represent the HCP for a node that has connectivity to the 3rd, 5th, 8th and 9th hubs. We define four target patterns, TP1 =(0 0 0 0 0 0 0 0 0 0), TP2 =(0 1 0 1 0 1 0 1 0 1), TP3 =(1 1 1 1 1 1 1 1 1 1) and TP4 =(1 0 1 0 1 0 1 0 1 0), hence 4 classes of nodes. We assign IDs to the nodes so that nodes that are classified into the same group have consecutive IDs.

CONCLUSIONS

In this paper we have delineated the effectiveness of random edge perturbation as a tool for privacy-preserving publication of graph data in the face of structural attacks. We have shown that sophisticated attacks based on probabilistic heuristics are feasible when the perturbation probability is sufficiently low, but can be thwarted as that probability grows. Moreover, we conducted an indepth theoretical study of the probability of success for any structural attack. Our analysis offers a yardstick for assessing the privacy risk entailed by the publication of perturbed graph data. On the other hand, we developed methods that accurately estimate the properties of the original graph from the perturbed data. A thorough experimental study validates our analysis and explores the tradeoff between privacy and utility. Last, we have effectively demonstrated that data mining tasks such as distance-based node classification can be successfully performed even under random edge perturbation probabilities with which structural attacks typically fail.

  (a) Original graph Go. (b) naïvely anonymized Gn. (c) Gn with malicious nodes. (d) Perturbed graph Gp.

Figure 1 :

 1 Figure 1: Example of a social network graph.

Figure 2 :• Pattern 4 :

 24 Figure 2: Converting a pattern in Go to another. • Pattern 4: They are completely disconnected with no edges.Similarly, we denote the original (observed) numbers of this pattern in Go (Gp) as Do (Dp).

k = 10 k

 10 = 20 k = 30 µ = 0.0001 0.9991 0.9981 0.9971 µ = 0.001 0.9910 0.9812 0.9714 µ = 0.01 0.9135 0.8262 0.7472

  It follows that r is in the range of [max{0, do -dp}, do]. The r neighbors can be removed in do r ways, and new neighbors added in Ndo -1 dpdo + r ways. Hence, Pr(dp|do) is equal to: r do r N -do -1 dp -do +r •(1-µ) N -1-(dp-do+2r) µ dp-do+2r

Algorithm 1 : 2 T 4 T 5 foreach neighbor nb of curnode do 6 localSearch

 12456 The probabilistic attack Data: Gp, GA, µ, w as chosen, m = 0, k-path; Result: A k-path containing identifiers of nodes in VA; 1 while k-path not found and wmax,mmax unreached do =new Tree(); level = 0; 3 foreach V [i] in Gp do 4 localSearch(V [i], level, T .root()); passes Interval Degree and ET-Edge checks then .add(curnode, parent); (nb, level + +, curnode);

1 Table 3 :

 13 λY with k = 10, M = 45, N = 10, 000.

Figure 3 :

 3 Time vs. w (prob.) (b) Time vs. m (prob.) Efficiency of the probabilistic attack.

5 Figure 4 :

 54 Figure 4: Evaluation of probabilistic attack.

Figure 5 :

 5 Figure 5: Evaluation of utility preservation.

1 Figure 6 :

 16 Figure 6: Classification of nodes under perturbation DBLP and Enron graphs, and average our results over 10 different randomly extracted subgraphs. Figures5g-lshow our results as a function of subgraph size, with µ = 0.001. density and transitivity, we show again the original, perturbed, and estimated values, for each subgraph size. For degree distribution, we now present the sum of absolute errors between the perturbed/estimate and original distribution values for each degree. Our results reconfirm our findings even with subgraphs extracted from the parent graphs.Overall, our results corroborate our claim that a published graph allows accurate estimations of graph properties not only for the whole graph, but also for subgraphs thereof, where publishing the measures directly to end-users would not be a satisfactory solution.

  Figure 6(a) visualizes the original classification. Figures 6(b,c,d,e) show the classification obtained from the perturbed graph with increasing µ. Remarkably, the classification is faithful even with µ = 0.01, with which both walk-based attack and probabilistic attack fail (see Figure 4). Classification error becomes apparent with larger µ as Figures 6(d,e) show, while most nodes are still correctly classified.

Table 2 :

 2 Pr(dp ∈ I) with N = 10, 000 and do = 50.

	001 0 µ = 0.001 2	0.3670 0.9814	5.9643 × 10 -6 0.7983
	µ = 0.001 4	0.9994	0.9931
	µ = 0.01 µ = 0.01	0 4	0.1246 0.8488	1.3935 × 10 -11 0.1399
	µ = 0.01	8	0.9927	0.9158

  the following, we study the conditional probability Pr(E Y p |E Y A ) for a particular Y . Given that Y is the set of embedded nodes and E Y A is the set of edges among the nodes in Y before perturbation, Pr(E Y p |E Y A ) is the probability that the set of edges among the nodes in Y becomes E Y p after perturbation. With this fact in mind, the derivation of Pr(E Y p |E Y A ) becomes easy: let m be the number of non-common edges in E Y p and E Y A . The minimum value of m is 0, achieved when E Y A and E Y p are identical, and its maximum value is M =

	k 2 -k 2 , obtained when E Y A and E Y p are totally foreign to each other. Since each edge removal or addition occurs with probability µ, the
	probability that E Y A is converted to E Y p is:

  1 http://www.cs.cmu.edu/˜enron 2 http://dblp.uni-trier.de/xml/ DBLP bibliography. Each author corresponds to a node and two authors are linked if they are coauthors in at least one paper. The Wiki dataset 3 is a network Wikipedia encyclopedia writers around the world. It consists of 7,115 nodes and 103,689 edges.

Table 5 :

 5 Effect of m.

	Events	m=0 m=1 m=2
	Success	53	58	59
	False prediction	35	35	35
	Broken path	6	6	6
	Edge check fail	6	1	0

ACKNOWLEDGEMENTS

We thank Yufei Tao for sharing insights on this topic, and the anonymous reviewers for their constructive feedback. Mingqiang Xue is supported by I 2 R grant and the work was started when he was a NUS Ph.D student.


Olivier Devillers
email: olivier.devillers@sophia.inria.fr

Monique Teillaud
email: monique.teillaud@sophia.inria.fr

Perturbations and Vertex Removal in a 3D Delaunay Triangulation

Though Delaunay triangulations are very well known geometric data structures, the problem of the robust removal of a vertex in a three-dimensional Delaunay triangulation is still a problem in practice. We propose a simple method that allows to remove any vertex even when the points are in very degenerate configurations. The solution is available in cgal 1 .

Introduction

Let us first briefly recall that the Delaunay triangulation of a set S of points in R 3 is the partition of the convex hull of the points consisting of the tetrahedra whose circumscribing ball does not contain any point of S in its interior. This definition implicitly implies that the tetrahedra are not flat, otherwise their circumscribing ball is not defined. In fact, it is always possible to build a Delaunay triangulation of a given set of points without any flat tetrahedron, as will be noticed in § 4.2. When a vertex v is removed from the Delaunay triangulation DT (S) of S, the tetrahedra incident to v are removed, which creates a polyhedral hole, and the interior of this polyhedron H must be triangulated with Delaunay tetrahedra. Triangulating H is exactly the inverse operation of inserting v in DT (S \ {v}). After the removal, the Delaunay triangulation is the triangulation that would have been obtained if v had never been inserted. When degeneracies occur (at least five cospherical points in S), the Delaunay triangulation is not unique. In general, any of the possible triangulations of the set of cospherical points can be returned. The difficulty arises when there are at least four cocircular points p 1 , p 2 , p 3 , p 4 , i.e. points that are both cospherical and coplanar, on the boundary of H. Indeed, in this case, for any point p 5 on H, the five points p 1 , p 2 , p 3 , p 4 , p 5 are cospherical, and there are two possible triangulations of this set of points, corresponding to the two different choices for the diagonal of the convex polygon p 1 , p 2 , p 3 , p 4 which is a facet of H (or a subfacet in the case when there are more cocircular points). Depending on this choice, we will get a different triangulation of this facet. In fact, there is no choice: since the outside of H is already triangulated, the triangulations of all the facets of H are given and must be respected.

Note that a similar problem cannot occur in the planar case, since the hole created by the removal of a vertex is a polygon, so, even in the case when this polygon has collinear vertices, there is no choice to be made: the edges of the polygon are uniquely determined by the order of the collinear vertices on the line that contains them. Any triangulation of the interior of the polygon will have the same edges on the polygon.

The problem can be seen as a special instance of the following question:

Is it always possible to compute a Delaunay triangulation of a given polyhedron H in such a way that the triangulation of its facets is respected?

(1.1)
In this paper, we show that, though this general question has a negative answer, the more restrictive instance of the problem posed by the vertex removal from a 3D Delaunay triangulation constructed incrementally can be solved, even in degenerate situations, and we propose a simple method to achieve it.

Previous Work

2.1 Algorithms. The deletion of vertices in two dimensional triangulations has been widely studied from both theoretical and practical points of view. A degree k vertex v can be removed in optimal Θ(k) time [START_REF] Aggarwal | A linear-time algorithm for computing the Voronoi diagram of a convex polygon[END_REF], [START_REF] Chew | Building Voronoi diagrams for convex polygons in linear expected time[END_REF], but in practice people often use sub-optimal solutions in O(k 2 ) or O(k log k) time [START_REF] Devillers | On deletion in Delaunay triangulation[END_REF].

In higher dimensions, the problem has received less attention since the complexity of the result can reach Ω(k 2 ). Thus, a worst-case optimal solution consists in computing the Delaunay triangulation of the neighbors of v in DT (S) (i.e. the vertices of H), then removing from this small triangulation the tetrahedra that are outside H and finally "sewing" this triangulation in the hole of DT (S) created by the deletion of v. Though this approach leads to a worst-case asymptotic optimal algorithm, it must be noticed that not only the hole is triangulated, but the whole convex hull of this hole, so, many tetrahedra are constructed and then deleted, which is useless memory management work.

To avoid this sewing technique and triangulate only the hole created by the deletion of v, another approach is the generalization of the 2D ear algorithm [START_REF] Heller | Triangulation algorithms for adaptive terrain modeling[END_REF], [START_REF] Devillers | On deletion in Delaunay triangulation[END_REF]. The hole H is a simple polyhedron; the general idea of the algorithm is to add new tetrahedra one by one in a way that ensures that the hole stays simple; such a tetrahedron is called an ear. There are two kinds of ears: ears formed by two triangles on the boundary of the hole, that share an edge, and ears formed by three triangles incident to a vertex that has degree 3 on the boundary of the hole. At each step, there may be several ears satisfying the Delaunay criterion; Devillers proves that the ear for which the removed point has the smallest power with respect to the circumscribing sphere belongs to the Delaunay triangulation [START_REF] Devillers | On deletion in Delaunay triangulation[END_REF], and he proposes to choose this ear. By maintaining a priority queue of the candidate ears, this yields a theoretical complexity of O(f log k), where f is the number of tetrahedra to be created.

Perturbation Techniques.

A general approach to solve the degenerate cases is the use of symbolic perturbations [START_REF] Seidel | The nature and meaning of perturbations in geometric computing[END_REF]. The general idea is to make the problem dependent of a parameter ε such that:

• there exists ε 0 > 0 such that the parameterized problem is in general position for ε ∈ (0, ε 0 ],

• if the original problem is in general position, the solution of the parameterized problem tends to the solution of the original problem when ε goes to zero,

• if the original problem is not in general position, the solution of the parameterized problem tends to something relevant when ε goes to zero by positive values.

More precisely, if S is not in general position, the Delaunay triangulation is not uniquely defined and the aim of a perturbation technique is to select one of the possible Delaunay triangulations.

Different perturbation methods can be used in that case: either the input of the problem can be perturbed or the definition of the problem can be slightly changed.

When using a perturbation, some properties can be lost in the result: to check whether they are satisfied, we must -check that the parameterized solution satisfies the properties -check that the properties are still true at the limit.

Perturbing the input can have very serious drawbacks: if the points move with ε [2], then a non flat tetrahedron can become flat at the limit. By flat tetrahedron, we mean a tetrahedron whose four vertices are coplanar.

As explained in § 4, we want to avoid this.

A standard way of changing the problem definition is to add to each point a weight depending on ε which goes to zero when ε goes to zero and to compute the regular triangulation (defined as the dual of the power diagram) instead of the Delaunay triangulation [START_REF] Edelsbrunner | Simulation of simplicity: A technique to cope with degenerate cases in geometric algorithms[END_REF].

One important advantage is that the points do not move, thus if a tetrahedron belongs to the limit solution, the same tetrahedron exists in the solution of a non degenerate problem and thus this tetrahedron is non flat. But the drawback of this technique is that some properties of Delaunay triangulations satisfied for nondegenerate input may be not satisfied by the limit solution (e.g. [START_REF] Dillencourt | Realizability of Delaunay triangulations[END_REF]). Another feature of the technique is that to assign to each point a function of ε as weight, we need a oneto-one map between points and functions, which is usually achieved by indexing the points. This is both a drawback and an advantage: it is a drawback because the result depends of that indexing; it is an advantage because the indexing can be chosen freely to fit other purposes.

Contributions of this Paper

In this paper, we propose a simple method based on the algorithm using ears presented in § 2.1. Rather than using Devillers' algorithm [START_REF] Devillers | On deletion in Delaunay triangulation[END_REF], we choose a non-optimal variant that consists in maintaining dynamically a list of candidate ears and testing if their circumscribing ball is empty by applying the in sphere predicate to all the vertices of the hole. When an ear is found to belong to DT (S \ {v}), then the hole and the list of candidate ears are updated.

Devillers proves that if the Delaunay tetrahedra are added by increasing powers, then they are ears on the hole at the moment they are added. In our case, the tetrahedra are added in a different order. If we consider H * the hole at some stage of the algorithm and T the Delaunay tetrahedron of smallest power inside H * . If we consider H T the hole before the addition of T in Devillers algorithm then by Devillers' result T is an ear of H T , thus T is also an ear of H * since H * ⊂ H T , which proves that there is at least one Delaunay tetrahedron T which is an ear. So, in all cases, there will exist Delaunay ears in the list of candidate ears at each stage.

The number of candidates during the whole algorithm is O(f ), the test takes O(k) for each ear, thus the whole complexity is O(f k), where, as in § 2.1, k is the degree of the removed vertex and f is the number of created tetrahedra.

Since both f and k are small constants in general, this method will be quite efficient in practice. For example, the expected value of k when the input points follow a Poisson distribution is 96 70 π 2 +2 15.5 and the expected value of f is 72 35 π 2 20 [START_REF] Okabe | Spatial Tessellations: Concepts and Applications of Voronoi Diagrams[END_REF]. The reasons for choosing this sub-optimal algorithm are the following:

• The algorithms proposed by cgal are parameterized by a traits class defining the basic geometric operations needed by the algorithms. The user can provides the algorithms with his/her own traits class. Therefore we need to keep the requirements for the traits class minimal. Devillers' algorithm would require an additional predicate for comparing the powers of a point with respect to two spheres, that is not necessary for the construction of the Delaunay triangulation.

• There are two ways of implementing Devillers' algorithm:

-A predicate comparing the powers of a point with respect to two spheres defined by four vertices of the triangulation can be used. In this case, the same subexpressions are computed several times: the power of the removed point with the same sphere has to be compared with the power of the same point with several other spheres.

-Avoiding redundant computations is possible by computing and storing the values of the powers. This leads to major problems: the fact of computing values instead of computing only signs yields important precision issues; an exact type could be used, which would give a bad efficiency. Filtering constructions such as in look [START_REF] Funke | Look: A lazy objectoriented kernel for geometric computation[END_REF] could in principle be used to improve the running time, but combining such a mechanism while keeping the genericity of cgal is still an open problem; this mechanism is not available yet in cgal, which proposes only filtered predicates.

• In order to cope with degeneracies, Devillers' method would need a perturbation of the comparison of the power of a point with respect to different ears. These power tests have higher algebraic degree than the predicates that are really needed by the computation of the Delaunay triangulation, namely the orient and the in sphere predicates (see § 5). A perturbation of these power tests would have to be compatible with the perturbation of the in sphere predicate, which would be quite involved.

Thus, though this algorithm could be coded very efficiently as a stand-alone program, its integration in cgal or in any library providing genericity and robustness is quite problematic. It seems to be better to use only the predicates that are intrinsic to the Delaunay triangulation, namely the in sphere and the orient predicates, and to avoid any construction.

In the sequel, we focus on the perturbation used for solving degeneracies rather than on the algorithmic issues.

• The perturbation method presented in § 5 consists in adding a symbolic weight to each point. It does not create any flat tetrahedron. Allowing flat tetrahedra would have allowed us to use more trivial algorithms, but, as mentioned in § 4, we consider this as unacceptable both for theoretical and practical reasons.

• We also remark that the solution can be slightly modified in such a way that it does not depend on the insertion order of the points and defines the Delaunay triangulation uniquely even for degenerate configurations.

The method is implemented in cgal. As far as we know, cgal is the only publicly available software proposing a fully dynamic 3D Delaunay triangulation. v is removed, whatever choice will be made, retriangulating the hole H creates edges that cross edges of the tetrahedra outside the hole. At first sight, the example of Figure 1 seems to prove that avoiding such edge crossings is impossible. In fact, this case cannot be built by an incremental algorithm coded as follows. As mentioned, the only problematic case is when there are cospherical (and coplanar at the same time) points. When a point p is inserted, the set of tetrahedra conflicting with p (defined as the tetrahedra whose circumscribing ball contains p in its interior) is determined, and these tetrahedra are deleted from the triangulation. If the algorithm considers as non-conflicting the tetrahedra whose circumscribing ball has p on its boundary, and thus if it does not delete these tetrahedra, then we get a unique construction of a Delaunay triangulation. This triangulation is unique for a given order of the points, but it depends on the order of insertion of the points.

Moreover, this incremental construction does not create any flat tetrahedron: if the point p to be inserted is coplanar with a triangle t that is a facet of the Delaunay triangulation, then the two tetrahedra incident to t will have the same conflict status with respect to p, which means that either they will both stay, and t will still be their common facet in the updated triangulation, or they will both be deleted and t will disappear. Thus, p will not form a flat tetrahedron together with t.

Assuming that the configuration shown in Figure 1 has been constructed by this incremental algorithm leads to the conclusion that the order of insertion of the points is cyclic, which is of course impossible. So, such a configuration cannot occur in our case.

Our problem can then be more precisely rephrased as: assuming that the Delaunay triangulation is constructed with the incremental algorithm considering each last inserted point p as not conflicting with tetrahedra whose vertices are cospherical with p, find the Delaunay triangulation of the hole that would have been computed if the removed point had never been inserted and if the other points had been inserted in the same order.

4.3 Allowing Flat Tetrahedra. Going back to the example in Figure 1, the new dashed edge that would be created to triangulate the prism would cross an edge of a tetrahedron that is exterior to the hole. The four vertices of these two crossing edges can be seen as a flat tetrahedron that could be added to the triangulation in order to make it combinatorially valid.

In the same way, if there are k ≥ 4 cocircular vertices on some facet of the hole created when removing a vertex of the Delaunay triangulation, the Delaunay triangulation of the vertices of the hole is not uniquely defined, and retriangulating it may lead to two different triangulations of the same facet: one given by the interior tetrahedra, and one given by the exterior tetrahedra. It is well known that any triangulation in the plane can be obtained from any other triangulation having the same vertices by a sequence of O(k 2 ) flips of diagonals of convex quadrilaterals [START_REF] Hurtado | Flipping edges in triangulations[END_REF]. Each such edge flip in 2D can be seen in 3D as a flat tetrahedron formed by the four (coplanar) vertices of the quadrilateral. Adding the O(k 2 ) flat tetrahedra would yield a combinatorially valid triangulation. Thus, allowing flat tetrahedra would allow us to get a simple solution to the problem. However, we do not consider this as an acceptable solution for the cgal users. Indeed, these flat tetrahedra correspond neither with the definition of a Delaunay triangulation, since the circumscribing ball of a flat tetrahedron is not uniquely defined, nor to the usual intuition. Moreover, in practice, this would lead to a heavier user code: before applying geometric operations to a tetrahedron, such as computing its circumcenter, the user would have to check that the tetrahedron is not flat.

Since we know that a Delaunay triangulation without any flat tetrahedron always exists (it can be computed for instance using the incremental construction seen in § 4.2), our goal is to find a method that does not create any flat tetrahedron.

Perturbing the in sphere Predicate

The in sphere predicate is the basic predicate used when inserting or removing a point in a Delaunay triangulation. Let p 0 , p 1 , p 2 , p 3 be four non-coplanar points in R 3 .

in sphere (p 0 , p 1 , p 2 , p 3 , p 4 )

> 0 if p 4 is outside = 0 if p 4 is on the boundary of < 0 if p 4 is inside    the ball circumscribing p 0 , p 1 , p 2 , p 3 .
Let us notice straightaway that, if the triangulation admits no flat tetrahedron, then this predicate will be needed only when the points p 0 , p 1 , p 2 , p 3 are actually non coplanar, and that no degenerate case has to be considered here: when inserting a new point p 4 , it is only tested against tetrahedra that are already existing in the triangulation; this is also the case when filling the hole formed by the removal of another vertex, because we only consider ears of the hole formed by non coplanar facets, since the goal triangulation does not contain any flat tetrahedron.

The decision whether to create or delete a tetrahedron (p 0 , p 1 , p 2 , p 3 ) in the Delaunay triangulation is easy when all the in sphere tests against all points p 4 give > 0 and < 0 results. When 0 results appear, decisions must be made using other criteria.

Notice that coding the insertion as mentioned in § 4, by considering that the last inserted point is not conflicting with tetrahedra such that the result of the in sphere test is 0 (i.e. considering that the point lies outside their circumscribing sphere), can be seen as using an implicit symbolic perturbation. It will be clear later that this perturbation corresponds to the perturbation proposed in the sequel.

It is well known that the in sphere test can be computed in the following way: in sphere (p 0 , p 1 , p 2 , p 3 , p 4 ) = sign Det(p 0 , p 1 , p 2 , p 3 , p 4 ) orient(p 0 , p 1 , p 2 , p 3 ) where, if point p i has coordinates (x i , y i , z i ) for each i, The sign Det(p 0 , p 1 , p 2 , p 3 , p 4 ) predicate in R 3 can be seen as an orientation predicate in R 4 , if each point p = (x, y, z) of R 3 is projected onto a point π(p) = (x, y, z, t) on the unit paraboloid Π of R 4 with equation t = x 2 + y 2 + z 2 [START_REF] Edelsbrunner | Voronoi diagrams and arrangements[END_REF], [START_REF] Devillers | The space of spheres, a geometric tool to unify duality results on Voronoi diagrams[END_REF]. Five points of R 3 are cospherical if and only if their projections on Π lie in the same hyperplane of R 4 . The symbolic perturbation of the in sphere test consists in adding respectively some value to the fourth coordinate of π(p 0 ), π(p 1 ), π(p 2 ), π(p 3 ), π(p 4 ) so that these points are not in the same hyperplane any more in R 4 . Then the predicate answers positive or negative instead of zero.

orient(p 0 , p 1 , p 2 , p 3 ) = sign 1 1 1 1 x 0 x 1 x 2 x 3 y 0 y 1 y 2 y 3 z 0 z 1 z 2 z 3 and Det(p 0 , p 1 , p 2 , p 3 , p 4 ) = 1 1 1 1 1 x 0 x 1 x 2 x 3 x 4 y 0 y 1 y 2 y 3 y 4 z 0 z 1 z 2 z 3 z 4 x 2 0 +y 2 0 +z 2 0 x 2 1 +y 2 1 +z 2 1 x 2
Let S be the set of n points {p 0 , p 1 , . . . , p n-1 } of R 3 .

We add here ε n-i to the fourth coordinate of the point π(p i ) of R 4 corresponding to p i . The quantity each point is perturbed with depends on its index. The way the points are indexed will be discussed in § 5.1.

Det(p i , p j , p k , p l , p m ) is perturbed into Det ε (p i , p j , p k , p l , p m ) = 1 1 1 1 1 x i x j x k x l x m y i y j y k y l y m z i z j z k z l z m t i +ε n-i t j +ε n-j t k +ε n-k t l +ε n-l t m +ε n-m where t = x 2 + y 2 + z 2 for = i, j, k, l, m
Developing with respect to the last raw yields a polynomial in ε Det ε (p i , p j , p k , p l , p m ) = Det(p i , p j , p k , p l , p m ) + orient(p i , p j , p k , p l )ε n-m -orient(p i , p j , p k , p m )ε n-l + orient(p i , p j , p l , p m )ε n-k -orient(p i , p k , p l , p m )ε n-j + orient(p j , p k , p l , p m )ε n-i From now on, the five points p i , p j , p k , p l , p m will be assumed to be non coplanar. As noticed at the beginning of this section, this assumption makes no restriction in practice since this is always the case when in sphere is used. When p i , p j , p k , p l , p m are cospherical, the constant term Det(p i , p j , p k , p l , p m ) of this polynomial vanishes to zero. At most one of the other coefficients orient() can be zero, otherwise two subsets with four points would consist of four coplanar points: either the five points would be coplanar, which contradicts our hypothesis, or the three common points of the two subsets would be collinear, which is impossible since these three points are cospherical. So, the polynomial Det ε is not identically zero.

Let us examine the coefficients of Det ε in order to determine the sign of Det ε (p i , p j , p k , p l , p m ). If m = max{i, j, k, l, m} In fact, in the case when p i , p j , p k , p l are coplanar, the sign of in sphere (p i , p j , p k , p m , p l ) does not depend at all of the comparison between m and i, j, k, l. It depends only of the comparison between i, j, k, l. When four points are cocircular (coplanar and cospherical), the most perturbed point (the point with highest index) is considered to be out of the disk circumscribing the other three points in their common plane.

5.1 Indexing the Points. The current choice in cgal for indexing the points is to use their order of insertion in the Delaunay triangulation. This choice corresponds to the implicit perturbation used for the insertion: at each step, the last point is more perturbed than all the previous ones, so, when it is cospherical with the four vertices of an already existing tetrahedron, it is considered as being outside the circumscribing ball. Some users reported this choice as being annoying for their application [START_REF] Tamal | Detecting undersampling in surface reconstruction[END_REF]. In fact, any other way of indexing the points could be used. We are considering the possibility of changing the indexing in the future: we can choose for instance the lexicographical order on the points, or any other (such as the comparison of memory addresses). This change would be straightforward in the vertex removal since it is already coded using the perturbation scheme, but it would require a change in the insertion: the perturbation would have to be made explicit.

Choosing the lexicographical ordering of points would have the advantage of giving a unique definition of the Delaunay triangulation, even in degenerate cases. On the other hand, it may lead to a slower construction of the triangulation for some input. We may think of leaving the choice of the order to the user.

Weights and Regular Triangulation.

As mentioned in § 2.2, the perturbation of in sphere can also be seen as adding a weight ε n-i to the points and computing the regular triangulation of the weighted points. This is an analogy with the "sliver exudation" method [START_REF] Cheng | Sliver exudation[END_REF], [START_REF] Cheng | Quality meshing with weighted Delaunay refinement[END_REF] that consists in associating weights to points, chosen so that the almost flat tetrahedra that are unavoidable in a Delaunay triangulation disappear in the regular triangulation. In our case, the tetrahedra are not almost flat, but really flat, the weights are symbolic, and the triangulation we compute is a Delaunay triangulation without any flat tetrahedron.

Conclusion

We proposed a method for removing a vertex from a Delaunay triangulation, that works even in degenerate situations. This method uses a symbolic perturbation technique that does not produce any flat tetrahedron. As can be seen in Appendix, the code for perturbing the in sphere predicate is quite simple. The solution is implemented in cgal (releases 2.3 and 2.4).

Acknowledgments. The authors wish to thank Jeff Erickson for informal discussions at the very beginning of this research, Andreas Fabri for coding the ears algorithm of § 3 in the cgal 3D triangulation package, Sylvain Pion for his improvements of the efficiency of this package, cgal users, especially Tamal Dey and his students, for encouraging us in improving cgal, and the members of the Prisme group at INRIA for preliminary discussions on the problem.

1 Figure 1 :

 11 Figure 1: (1.1) has negative answer

  • if p i , p j , p k , p l are non coplanar, then DetWhen five points are cospherical, the point with highest index among the five points will always be considered as lying outside the sphere defined by the four other points, if they are non coplanar.• if p i , p j , p k , p l are coplanar, then, by hypothesis, p m does not lie in the same plane. The predicate in sphere (p i , p j , p k , p l , p m ) has no geometric meaning since p i , p j , p k , p l do not define a sphere. Let us assume that l = max{i, j, k, l} (the other cases are similar, they just reduce to permutations in the columns of the determinant Det ε ). As noticed above, p i , p j , p k , p m define a sphere and in sphere (p i , p j , p k , p m , p l ) has a geometric meaning.Det ε (p i , p j , p k , p m , p l ) = -Det ε (p i , p j , p k , p l , p m ) has the sign of orient(p i , p j , p k , p m ).

ε (p i , p j , p k , p l , p m ) has the sign of orient(p i , p j , p k , p l ), thus in sphere (p i , p j , p k , p l , p m ) > 0 and p m lies out of the ball defined by p i , p j , p k , p l .

// as seen in Section 5, at most one orientation test // is zero, so a result has been returned }