{'fraction_non_alphanumeric': 0.039985133335398146, 'fraction_numerical': 0.01235791495028959, 'mean_word_length': 3.8766047424860295, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On role extraction for digraphs via neighbourhood pattern similarity

We analyse the recovery of different roles in a network modelled by a directed graph, based on the so-called Neighbourhood Pattern Similarity approach. Our analysis uses results from random matrix theory to show that when assuming the graph is generated as a particular Stochastic Block Model with Bernoulli probability distributions for the different blocks, then the recovery is asymptotically correct when the graph has a sufficiently large dimension. Under these assumptions there is a sufficient gap between the dominant and dominated eigenvalues of the similarity matrix, which guarantees the asymptotic correct identification of the number of different roles. We also comment on the connections with the literature on Stochastic Block Models, including the case of probabilities of order log(n)/n where n is the graph size. We provide numerical experiments to assess the effectiveness of the method when applied to practical networks of finite size.


Introduction
The analysis of large graphs frequently assumes that there is an underlying structure in the graph that allows us to represent it in a simpler manner. A typical example of this is the detection of communities, which are groups of nodes that have most of their connections with other nodes of the same group, and few connections with nodes of other groups. Various measures and algorithms have been developed to identify community structures [60] and many applications have also been found for these model structures [36,4,59,58]. Yet, many graph structures cannot be modelled using communities: for example, arrowhead and tree graph structures, which appear in overlapping communities, human protein-protein interaction networks, and food and web networks [4,48,67]. These more general types of network structures can be modelled as role structures, and the process of finding them is called the role extraction problem, or block modelling [11,12,51,71,72]. The role extraction problem is a generalization of the community detection problem and it determines a representation of a network by a smaller structured graph, where nodes are grouped together based upon their interactions with nodes in either the same group or in different groups called roles. If no a priori information is available, one needs to verify all possible group and role assignments in order to determine the best role structure for the data, which leads to an NP-hard problem [21,12] for both the community detection problem and the more general role extraction problem.
There are many algorithms proposed for community detection, both for directed and undirected graphs [33,60,43,26,63,49], but they often do not state any conclusive results about the exact recovery of communities, because they make no statistical assumption about the underlying model of the graph. On the other hand, if one assumes that the adjacency matrix of the graph is a sample of a random matrix that follows certain rules, then the problem of recovering the correct underlying block structure may become tractable. The Stochastic Block Model (SBM) is precisely such a model: the interactions between all nodes of a particular group with all nodes of another group follow exactly the same distribution [40]. There is a considerable literature on SBM [37,51,76], including variants that address diagonal scaling of the SBM [68].
A class of algorithms that has been largely employed in the past years for such purpose are the so-called spectral methods [47,75,10,34,35]. Broadly speaking, a spectral method first performs an eigendecomposition of a symmetric matrix encoding the properties of the graph. Then the community membership is inferred by applying a clustering algorithm, typically K-means, to the rows of the matrix formed by the first few leading eigenvectors. Spectral clustering is easier to implement and computationally less demanding than many other methods, which amount to computationally intractable combinatorial searches. From a theoretical standpoint, spectral clustering has been shown to enjoy good theoretical properties in stochastic block models [24,42,18]. In the computer science literature, spectral clustering is also a standard procedure for graph partitioning and for solving the planted partition model, a special case of the SBM [45].
As their first step requires the eigendecomposition of a symmetric matrix, spectral methods are commonly applied to undirected graphs. Moreover, when they do consider directed graphs, their analysis does not include the recovery of the underlying block structure [62].
In this paper, we will show that a particular method, using the so-called Neighbourhood Pattern Similarity (NPS) matrices [22,54], allows us to give a positive answer to the following question: Can we recover asymptotically the block structure of a general directed graph with a stochastic block model structure? A NPS matrix is a real symmetric positive semi-definite matrix, also for directed graphs, and therefore has real eigenvalues and eigenvectors. We then show that for sufficiently large graphs, the gap between dominant and dominated eigenvalues allows a convergent recovery of which nodes are associated with the different roles in the model. The nearest results available in the literature are the successful extraction of the correct roles in SBM for the community detection problem of undirected graphs [2,51], and the use of spectral clustering for the directed role extraction problem [69], in which a different type of stochastic block model is used. The present paper extends the asymptotic analysis of the general role modelling problem to specific symmetrizations of the standard SBM model for directed graphs, the NPS matrices, for which a correctness result is still missing in the literature. This is in particular one of the few existing results of correctness for the spectral clustering algorithm applied to directed graph with a SBM structure. Furthermore, our results can be seen as covering a whole class of methods, in the sense that our asymptotic analysis applies to all the admissible values of the scaling factor β and all NPS matrices S k (including both any finite value of k and the limit S = lim k→∞ S k , which is the NPS matrix); see Subsection 2.3 for the definitions of β and S k . While in this paper we focus on the theoretical analysis of the method, determining optimal values of k and β for a practical implementation of the algorithm to analyse actual graphs is an interesting possible subject of future research.
In Section 2 we go over several preliminaries related to graphs, random matrices, stochastic block models and role modelling. Section 3 then yields the spectral bounds for the NPS matrix associated with the graph (as well as for the matrices S k whose limit is the NPS matrix) and in Section 4 we describe the asymptotic behaviour of the clustering error. In Section 5 we give a few numerical experiments illustrating our theoretical analysis, and we conclude with a few final remarks in Section 6. Several technical proofs are moved to the Appendices for the sake of readability.

Graph theory and role extraction
An unweighted directed graph, or digraph, G = (V, E) is an ordered pair of sets where the elements of V = [n] are called vertices, or nodes, while the elements of E ⊆ V × V are called edges, or links. A walk of length on the digraph G from i 1 to i +1 is a sequence of vertices of the form i 1 , i 2 , . . . , i +1 such that for all j = 1, . . . , , (i j , i j+1 ) ∈ E. G is said to be strongly connected if for all i, j ∈ V there is a path on G from i to j.
The adjacency matrix of an unweighted digraph is defined as In particular, A is a non-negative matrix and A ∈ {0, 1} n×n . It is well known that A is irreducible if and only if G is strongly connected; in that case, the Perron-Frobenius spectral theory for irreducible non-negative matrices applies. Manifestly, there is a bijection between adjacency matrices and digraphs. Moreover, two digraphs are isomorphic (i.e. they coincide up to a relabelling of the vertices) if and only if their adjacency matrices are permutation similar. Given a graph G with adjacency matrix A, the problem of role extraction consists in finding a positive integer r ≤ n and an assignment function ξ : [n] → [r] such that A can be well approximated by an ideal adjacency matrix E such that E ij only depends on ξ(i) and ξ(j). Equivalently, if π is any permutation that reorders the nodes such that nodes of the same group are adjacent to each other, and P is the corresponding permutation matrix, then P T AP is approximately constant in the r blocks induced by the assignment. One can then associate with it a so-called ideal adjacency matrix A id as illustrated in Figure 1 : for the blocks of A where 1 dominates, put all elements of the corresponding block in A id equal to 1 and for the blocks where 0 dominates, put them all equal to 0 in A id . In doing so, the nodes in each block of A id are regularly equivalent [20], i.e. they have the same parents and the same children [22,54]. The above approximation problem for P T AP can thus be viewed as finding a nearby regularly equivalent graph to a given graph. Figure 1: Associating a regularly equivalent graph A id to the permuted graph P T AP The role extraction problem can be also generalized into finding clusters of nodes C 1 , . . . , C q , or equivalently an assignment function, such that for a given node x in cluster C a , the number of edges between x and the cluster C b depends only on a and b. In the next section we will find that this problem is better formalized by the Stochastic Block Model.

Stochastic Block Model
Let (Ω, F, P) be a probability triple and consider the space of random variables Ω → {0, 1} a×b . A random digraph G(ω), ω ∈ Ω, is a graph whose adjacency matrix A(ω) is one such random variable. We denote the expectation of a random matrix A by E[A]. We construct a random unweighted digraph as follows: 1. The nodes are partitioned in q clusters of nodes, C 1 , C 2 , . . . , C q , of size m 1 n, m 2 n, . . . , m q n respectively.
2. There is an edge between a node in cluster C a to a node in cluster where θ a,b depend only on a and b and max a,b θ a,b = 1.
Since p a,b are probabilities, necessarily f (n) = O(1), but from classical information theory we know that exact recovery for the clusters requires [1,57] and it is also more restrictive than the sufficient condition for clustering detection [56]. The adjacency matrix A n of such a random graph is an mn × mn random matrix, where m = q i=1 m i . Suppose that m i may vary with n, but for every n, i we have 0 < m min ≤ m i ≤ m max where m min , m max , q, θ a,b do not depend on n. As a consequence m may vary, but it is always bounded between absolute constants qm min ≤ m ≤ qm max . Suppose i ∈ C a , j ∈ C b : then, A ij is distributed as a Bernoulli variable centered on {0, 1} with P(1) = p ab . In this section, we assume that the nodes of the same cluster are adjacent to each other, in order to simplify the notation. This does not affect the generality of our results.
Denoting by 1 k ∈ R k the vector of all ones and by M n = E[A n ], then then Z n := Z n D −1 / √ n has orthogonal columns and the nonzero singular values of M n are those of nf (n)DΥD. We have in particular that The above scenario is what arises in the theory of Stochastic Block Model, but in most references the matrix Υ is taken symmetric. In the following sections we will analyse the model described above, together with a spectral method designed to extract the clusters, which will be called roles, through the use of a similarity matrix S. For this reason, we report here a result we will need in our arguments about the matrix Y n := A n − M n .
almost surely.
Since the entries of when n is big enough. In what follows, we will bound the norm of here is not sharp. In fact, both Theorem 4.1 in [39] and the experiments we will present suggest that the result holds with the tighter constant 1 + √ 2, following the classical bound on the Marchenko-Pastur distribution. Since there is no such result in literature, we formulate it here as a conjecture. Conjecture 1. Let Z N be N ×N random matrices with independent, mean zero, and uniformly bounded entries. Suppose that σ 2 bounds the second moments of all entries, independently on N . If we call almost surely. Moreover, if every entry has the same second moment σ 2 , the bound is attained.
From now on, when we say "for any n big enough" or a similar formulation, we always implicitly mean that the result holds almost surely.

Role extraction via the similarity matrix S
In [22,54], it was proposed to solve the problem of role extraction for a digraph with adjacency matrix A by means of a Neighbourhood Pattern Similarity matrix S, which is defined as the limit of the sequence of SPD matrices where the operator (depending on the matrix parameter W ) Γ W is defined as It was shown in [22,54] that the sequence is convergent if and only if It was also shown there that element (i, j) of the matrix S k is the weighted sum of the walks of length up to k between nodes i and j, and that this can exploited to find nodes that should be associated with the same role. Note that S 1 is a known symmetrization for direct graphs called "Bibliometric Symmetrization" [63]. Throughout this document, the parameter β in (5) is always assumed to satisfy β 2 Γ A < 1, which is sufficient for the sequence (S k ) k∈N to converge. Here and thereafter, we measure the norm of the operator Γ W induced by the spectral norm of its matrix argument. More concretely, X .
Lemma 1. The norm of the linear matrix mapping Γ W : From the previous result, whose proof can be found in the appendices, Γ W ≤ 2 W 2 , so that it is easy to compute a good enough β with very low computational effort. In fact, we can always choose, for example, β 2 = 1/(4 A 2 ) and obtain that necessarily β 2 Γ A ≤ 1/2.
Consider the Stochastic Block Model described in Section 2.2. From now on, we drop for the sake of notational simplicity the suffixes emphasizing the dependence on the size n, so, for example, we simply write A, M, Z, Y for A n , M n , Z n , Y n . Given the random adjacency matrix A, suppose that f (n)Υ is a minimal role matrix, defined as follows.
a deterministic block matrix, and the following result shows that it is possible to recover the original clusters by analysing any of the matrices T k generated as the S k but replacing A with M .
the orthogonal matrix in the reduced SVD (or, equivalently, reduced eigendecomposition) of T k , it follows that the set of the vectors of R r that are a row of V k has precisely q distinct elements. Moreover, the q original clusters of the graph coincide with the partition of [mn] into the q subsets associated with the row indices that correspond to each distinct vector that is a row of V k .
As a consequence, it is enough to perform an eigendecomposition of T k , extract the reduced orthogonal matrix V k and then identify the repeated q rows to recover the clustering. A natural thought is to try and apply the same method to the random symmetric matrix S k generated by the recurrence (5), but some issues arise.
• T k has rank r ≤ q, while S k is with high probability full rank, so we need a way to determine the truncation parameter r for the SVD.
• In the truncated eigendecomposition of S k , the orthogonal matrix U k has usually distinct rows. In order to retrieve the clusters, we thus need to estimate the number of roles q and perform a K-means algorithm on the rows.
There is method to do this, commonly referred to as Spectral Clustering of the matrix S k . A detailed description is given in [54] and a concise description as pseudocode is given below: • Inputs: adjacency matrix A, number of roles q, scaling factor β, integer k.
• Output: a partitioning of the nodes of the graph into q clusters.
• Procedure: -Compute the matrix X 1 whose columns are the q dominant singular vectors of A A T ; -For h = 2, .., k compute the matrix X h whose columns are the q dominant singular vectors of Y h = βAX h−1 βA T X h−1 X 1 ; -Apply the K-means algorithm to the rows of the matrix X k .
The sparse singular value decomposition of [A, A T ] can be computed using the Lanczos bidiagonalization procedure [38] and its complexity is O(µq 2 ) because each matrix vector multiplication requires exactly 2µ flops, where µ is the number of edges in the graph, i.e., the number of nonzero entries of the mn × mn matrix A. For the same reason, the construction of the matrix Y h requires exactly 2(µ + mn)q flops. The singular value decomposition of the economy size singular value decomposition of the dense mn × 3q matrix Y h , requires O(mnq 2 ) + O(q 3 ) flops [13]. Altogether, we thus have a complexity of O(kq(µ + mn + mnq + q 2 )) to compute the low rank factor X k , which scales well with µ. The subsequent clustering of the rows of X k is then constrained to a q-dimensional space, and requires on average O(mnq 2 ) flops per iteration of the K-means algorithm [53].
In the next sections, we show that the matrices S k sport a clear gap between the eigenvalues λ r (S k ) and λ r+1 (S k ) that lets us identify the rank r with high probability for big n. Moreover, when the matrix f (n)[Υ Υ T ] is full-rank, so that f (n)Υ is minimal and r = q, we estimate the clustering relative error for the K-means algorithm on S k , and show that it is proportional to (nf (n)) −1 .

Spectral Bounds
We now consider the recurrence relation using the expected value M rather than A since this yields a good approximation for the S k matrices. We denote these matrices as T k and their recurrence is thus given by Note that in (6) the matrix parameter in the operator Γ is set to M = E[A], in contrast with (5) where it was set to A. Again, the parameter β 2 is chosen such that β 2 Γ M < 1, which is required for the sequence T k to converge to T = lim k→∞ T k . In order to choose an appropriate β we need an estimation of Γ M depending only on the matrix A.
Lemma 2. Let δ 2 = 4mnf (n). For n large enough, it holds Using the last result and Lemma 1, we find that for n large enough, β 2 Γ M ≤ 3β 2 A 2 and β 2 Γ A ≤ 2β 2 A 2 so from now on, we always suppose that β 2 ≤ 1/6 A 2 and consequently It was pointed out in [54] that the matrices S k and T k are all positive semi-definite, and that both sequences are ordered in the Loewner ordering : Moreover, as shown in Theorem 2, if T k were available then we would be able to recover exactly the original clustering that generated the random directed graph. Since we can only work on S k , that is an approximation of T k , it is essential to analyse the proximity between the two matrices more accurately. This will let us study how well the properties of T k transfer to S k and how effective is a spectral clustering algorithm applied to S k .
where the last inequality holds also for S − T .

Spectral Gap
From Theorem 2 we know that T k has rank r, so it stands with reason to expect that its approximation S k has r dominant singular values (which we refer to as the "signal") and mn − r small singular values (which we refer to as the "noise"). Here we report estimations for the eigenvalues of S k and T k and then derive bounds on the respective gaps. for every k ≥ 1 and n big enough.
for every k ≥ 1 and n big enough.
The gaps λ r (S k ) − λ r+1 (S k ) and λ r (S k )/λ r+1 (S k ) between signal and noise, are expected to be large enough to allow for the correct truncation for the SVD of S k , and a correct assignment of the different nodes in each "role", as we will show in the next section. This separation becomes more pronounced when the dimensions of the matrix and its subgroups increase, as we can see by applying Lemma 3 and Theorem 4.
• For the absolute gap, that is order of magnitudes greater than the following absolute gaps, since • For the relative gap, that is order of magnitudes greater than the previous relative gaps, since As a consequence, a comparison of the gaps between signal and noise with the other gaps is a clear indicator of the right rank r with which one should perform the truncated SVD in the algorithm. This holds also for the limit matrix S. We can note that all the estimates get worse as σ r ([Υ Υ T ]) gets close to zero. This has to be expected since it is harder to compute the rank of an almost singular matrix. In fact, for example, in the case where all the probabilities θ a,b are close to each other, it is harder to distinguish between different groups and clusters.

Dominant Subspaces
In this subsection, we study the dominant subspace of a real symmetric matrix S k (i.e. the invariant subspace associated with the largest r eigenvalues) and argue why, for sufficiently large n, it allows role extraction. Classically, distances between subspaces are measured via the concept of principal angles [44]: a multidimensional generalization of the acute angle between the unit vectors u, v, i.e., 0 ≤ θ(u, v) := arccos |u T v| ≤ π/2. More generally, if U, V are subspaces whose orthonormal bases are given, respectively, as the columns of the matrices U, V , then the min{dim U, dim V} largest singular values of U T V are the sines of the principal angles between U ⊥ (orthogonal complement of U) and V. Just as in the one-dimensional case, the principal angles between two subspaces are all zero if and only if the subspaces coincide, and more generally the smaller the principal angles the closer the subspaces.
To set up notation, fix k ∈ N, and let E and F be the dominant subspaces of dimension r for S k and T k respectively. By classical results in geometry and linear algebra [30,41], the r largest singular values of the matrix are the sines of the principal angles between the dominant subspaces of T k and that of S k , where Π E and Π F are the orthogonal projection matrices on the relative subspaces. Hence, the spectral norm of sin Θ measures how well the dominant subspace of the similarity matrix S k approximates the one of the ideal graph. We rely on Davis-Kahan's sine theta theorem [30], in the form given by [41,Theorem 5.3]. Call S k the best r-rank approximation of S k . Since the r-th eigenvalue of T k is larger than the (r + 1)-th eigenvalue of S k (which is 0), the assumptions of [41,Theorem 5.3] apply and thus Remark 1. We could apply [41,Theorem 5.3] reverting the roles of S k and T k , obtaining In this case, though, we prefer to deal with the deterministic quantity λ r (T k ) instead of the aleatory λ r (S k ), even if the estimation gets worse by a constant factor 2.
Using the results of the previous section, in turn this yields for sufficiently large n Therefore, we can state Corollary 1. Asymptotically as n → ∞, the principal angles between the dominant subspaces of S k and T k tend to 0 at least as fast as δ −1 .

Clustering Error
In the previous sections we have estimated how close the matrix S k is to the deterministic matrix T k and how this influences their spectral properties and their dominant subspaces. Here we show that the same estimates can be used to bound the clustering error of the proposed method on S k , under the technical hypothesis r = q that is, the matrix f (n)[Υ Υ T ] is full rank. Note that Υ is still allowed to be singular.
Recall that the model is generated by the clusters C 1 , . . . , C q , where C i has cardinality n i := m i n. Suppose that T 1 , . . . , T q are the resulting clusters from the algorithm operated on the similarity matrix S k . Define the misclassification error f as where is the symmetric difference of sets defined as the elements belonging to exactly one of the two sets, or equivalently A B := (A \ B) ∪ (B \ A). S q is the q-th symmetric group, that contains all the permutations on q elements. f is thus a measure of the maximum rate of misclassified points over all clusters, up to the assignment of the correct clusters C i to the T i derived by the algorithm. In the appendix, we give a proof for the following bound on f . Theorem 5. There exists an absolute constant C such that asymptotically in n Note that the error goes asymptotically to zero as long as nf (n) → ∞, which is exactly condition (1).

Remark 2.
The proof of the Theorem follows the same steps as [69] and [46]. In particular, in the former we find a similar algorithm applied directly on the adjacency matrix A instead of S k , but the analysis is limited to the case where f (n)Υ has full rank, while we work under the more general condition that f (n)[Υ Υ T ] is full rank.
Observe moreover that all the results of Section 3 hold without any assumption on Υ, so we still have all the spectral bounds and the convergence of the dominant subspace of S k to the one of T k also in the general case.
Yet, for the algorithm to make sense, we need f (n)Υ to be a minimal role matrix as defined in Definition 1. Moreover, since r < q, it is necessary to apply the K-means algorithm on S k for K = r, r + 1, r + 2, . . . and look at the error in order to find the optimal number of clusters.

Numerical examples
In this section we illustrate the theoretical results of the paper using an example generated according to the rules of a Stochastic Block Model where only two different probabilities are used, namely p and 1 − p. We chose q = 3, m 1 = m 2 = m 3 = 10 and hence m = 30, and We then ran simulations for matrices A with n = 10 : 50.
In Figure 2 we took β = 0 which means that the sequences S k and T k are constant after one step, and hence that S = S 1 and T = T 1 . The q dominant eigenvalues of S 1 are the circles in each plot (because of the structure of Υ there are two repeated ones). The full lines are their estimates obtained from the rank q matrix T 1 , and clearly they are very accurate estimates as indicated in Theorem 3. The squares correspond to the "noise" eigenvalue λ r+1 (S 1 ) and the dashed line is its estimate (3 + √ 8)p(1 − p)mn according to Conjecture 1 and λ r+1 ( It is clear from these plots that this is also a very good estimate and that the ratio λ r (S 1 )/λ r+1 (S 1 ) grows like O(n). Moreover, the plots show that the gap |λ r (S 1 ) − λ r+1 (S 1 )| shrinks with p getting closer to 0.5, which is expected since for p = 0.5 the rank of Υ drops to 1. This means that for p getting closer to 0.5, one has to require larger dimensions of the graph in order to recover an accurate enough grouping. In Figure 3, we performed the same experiment, but now with β chosen such that β 2 Γ A (·) ≈ 1 2 , which guarantees convergence of the method. In order to reduce the complexity of the method, we computed S 10 and T 10 rather than the limits S and T , since in 10 steps we should have reasonably good estimates of these limits. We can see from the plots that one has to wait for larger values of n to reach a sufficiently large gap |λ r (S 10 ) − λ r+1 (S 10 )| than for |λ r (S 1 ) − λ r+1 (S 1 )| in Figure 2.
In Figure 4, we computed the misclassification errorf of the clustering associated to the matrix Υ for p = 0.6. Using the same parameters m i , n as before we show an averagedf over 60000 instances for the clusters extracted from S 1 and S 10 , where we took β = (2 [A A T ] 2 ) −1 . For comparison, we also plot the function 3/(10n + 24) and note that it fits well both plots, thus confirming the bound O(1/n) predicted by Theorem 5.

Concluding remarks
In this paper, we showed that the Neighbourhood Pattern Similarity matrices S k of a directed graph with adjacency matrix A have spectra that are well separated into two groups of eigenvalues, provided that the graph is sufficiently large and that it is generated according to the Stochastic Block Model with blocks where all elements in each block follow a Bernoulli distribution with the same probability O(f (n)) where nf (n) → ∞.
The large eigenvalues are then associated with the nonzero eigenvalues of the expected value E[A], which is a low rank matrix, and the small eigenvalues are associated with the mean and variance of the random distribution used in the Stochastic Block Model. Moreover, the gap between the "large" eigenvalues and the "small" ones This analysis was primarily based on the recurrences defining the matrices S and T and on the fact that the underlying adjacency matrix is generated according to a Stochastic Block Model. It is likely that our results can be extended for other types of distributions and that weighted graphs can also be dealt with, but our analysis here was limited to unweighted adjacency matrices for directed graphs.
We point out that the same analysis could in principle be conducted in the sparse limit case f (n) = O(1/n), but since most of the results are formulated asymptotically in nf (n) one has to explicitly compute all the implicit multiplicative constants. A technical work of refinement is also needed on each proposition to obtain the best constants and thus meaningful results. For these reasons, we postpone the limit sparse case analysis to future work. Note that for all X Thus we have Γ W ≤ W W T 2 ≤ 2 W 2 , whose first bound is satisfied for X = I since

A.2 Proof of Lemma 2
For any matrix X, if we rewrite we readily see that Using (3) and (4), and recalling that δ 2 = 4mnf (n), Then the desired bound follows easily. As for the last bound, note that δ 2 ∼ nf (n) and in virtue of (1), (2) and (4) we have If we call C the constant Υ m min /4qm max , then

B.1 Proof of Theorem 3
Denoting the increments ∆ S i+1 := S i+1 − S i and ∆ T i+1 := T i+1 − T i we obtain Observing that Hence, by induction on k ≥ 1, it is not difficult to obtain the upper bound , which finally yields the bound In virtue of (7), we can let k go to ∞ and find that and the same holds for M , thus the desired bound follows.

C Clustering Error
Theorem 6. [73, Theorem 1] Call M the set of nm × q matrices that have only q distinct rows. Let E, F be mn × q matrices, where F ∈ M, whose rows µ 1 , . . . , µ q identify the clusters C i . Call n i := |C i | and Suppose there exists ρ ≥ 100 such that µ i − µ j ≥ ρ(∆ i + ∆ j ) for any i = j. Let G ∈ M be a 10-approximation of the K-means algorithm on E, that is and call ν 1 , . . . , ν k the rows of G. Partition the indices 1, . . . , mn into q clusters T i according to G and E as in T r := {i : E i,: − ν r ≤ E i,: − ν s ∀s}. We have that there exists a permutation π and an absolute constant C such that |C r T π(r) | ≤ Cn r /ρ 2 for every r.

C.1 Proof of Theorem 5
In order to analyse the method, we need first to better characterize the eigenvalue decomposition (EVD) of T k . In fact, from Theorem 2, we know that there exists a full-rank PSD matrix T k such that T k = Z T k Z T . Recall now that D = diag( √ m 1 , . . . , √ m q ), and that Z = ZD −1 / √ n has orthonormal columns. If D T k D = U k Σ k U T k is its EVD, then T k = Z T k Z T = n ZD T k D Z T = n ZU k Σ k U T k Z T , where W k := ZU k = ZD −1 U k / √ n has orthogonal columns, so that W k nΣ k W T k is the q-reduced EVD of T k . Note that W k ∈ M since its rows coincide with the ones of D −1 U k / √ n, that is a full rank q × q matrix. For the same reason, we have that W k Q ∈ M for every orthogonal q × q matrix Q. Moreover, the clustering induced by W k and W k Q are the same, and coincide with the original clustering C 1 , . . . , C q . It follows that if W k , the orthogonal matrix in the q-truncated SVD of S k , is close to W k Q for even one matrix Q, then it has good chance to generate a good clustering. Here we report two results formalizing the concept.
By Lemma 4 and (9) there exists a k × k orthogonal matrix Q such that where sin(Θ) are the sines of the principal angles between the subspaces W k and W k , thus If we call µ 1 , . . . , µ q the distinct rows of W k Q = ZD −1 U k Q/ √ n, then they are in the form u i / √ nm i where u i are the rows of U k Q, that is an orthogonal matrix, so .
By Corollary 1, sin(Θ) = O(δ −1 ), so ρ > 100 for n big enough. The K-means algorithm applied to the matrix W k outputs the clusters T 1 , . . . , T q and Theorem 6 assures us that there is an absolute constant C for which We can finally conclude that by (10) and incorporating all the absolute constants into C,