Introduction {#s1}
============

The theory of network science has significantly improved our understanding of complex systems. In the last fifteen years, an explosive growth of interests and activities on the structure and dynamics of complex networks is appearing. One of the most favorable but challenging tasks in network science is cluster analysis, which is aimed at revealing possible partitions of a network into subsets of nodes (clusters, or communities) [@pone.0082964-Lancichinetti1]. Markov chains are frequently used as analytic models in the quantitative evaluations of stochastic systems. Examples of their use may be found in diverse areas such as computer, biological, physical and social sciences as well as in business, economics and engineering [@pone.0082964-Albert1]--[@pone.0082964-Brner1]. The fundamental property characterizing the model, referred to as the Markov property, is that given the present, past and future transitions of the system are independent of each other [@pone.0082964-Shi1], [@pone.0082964-Meila1]. The information that is often sought from this model is either the transient or stationary probability of the system being in a given state. When the number of states is small, it is relatively easy to obtain the transient and stationary solutions allowing the prediction of the system behavior. However, as models become more complex the process of obtaining these solutions becomes much more difficult. There is also a wide class of situations, where the modeler does not need information about each state of the system but about classes of states only. This leads to the consideration of a new process, to be called the aggregated or lumped, whose states are the state classes of the original Markov chain. The new stochastic process need not to be Markovian. In order to be able to utilize all the power of the Markov chain theory, it is important to be able to claim that for a given initial distribution the aggregated process has the Markov property.

In a previous paper [@pone.0082964-Weinan1], a *k*-means approach is proposed to partition the networks based on optimal prediction theory proposed by Chorin and coworkers [@pone.0082964-Chorin1], [@pone.0082964-Chorin2]. The basic idea is to associate the network with the random walker Markovian dynamics [@pone.0082964-Lovsz1], then introduce a metric on the space of Markov chains (stochastic matrices), and optimally reduce the chain under this metric. The final minimization problem is solved by an analogy to the traditional *k*-means algorithm [@pone.0082964-Devijver1], [@pone.0082964-Hastie1] in clustering analysis. This approach is motivated by the diffusion maps [@pone.0082964-Lafon1] and MNCut algorithms in imaging science [@pone.0082964-Shi1].

To give a coarse definition about the study of complex networks from the viewpoints of applied mathematics, it is about the research of dynamical systems on graphs. The graph structure may be fixed, or time-varying; the dynamical system may be deterministic, or stochastic. Since these networks are typically very complex, it is of great interest to see whether they can be reduced to much simpler systems [@pone.0082964-Shi1]--[@pone.0082964-Weinan1], [@pone.0082964-Lafon1]--[@pone.0082964-Hofman1]. And in a broader aspect, it is also closely related to the model reduction theory of differential equations [@pone.0082964-Schilders1], [@pone.0082964-Danon1]. The concept of lumpability on hard partitions is a useful tool to analysis the dynamic of the network based on the a coarse grain view. Lumped Markov chain model of the random walker (i.e., a reduced-order Markov chain in which the clusters of the original network become nodes) which is easily derived from the original (highorder) Markov chain model. This notion of proximity of 'nodes' in the lumped markov chain reflects the intrinsic geometry of the meta-node set in terms of connectivity of the cluster in a diffusion process.

In the current paper, we first analyze the property of lumped markov chain and propose two novel ways of inferring the lumped markov transition matrix. Some useful results are proposed based on the analysis of the properties of the lumped markov process. Furthermore, we construct two algorithms --- the steepest descent method with projection (SDP) and the reduced conjugate gradient method with projection (CGP) --- from minimizing the objective function under the generalized framework in this paper. According to two choices of projection operators P1, P2, we obtain the formulations --- SDP1, SDP2, CGP1, CGP2 --- which have been applied to two artificial networks, including the ad hoc network, as well as real-world network, the karate club network. The proposed algorithms are easy to be implemented with reasonable computational effort and the final results do make sense in the considered models. It is demonstrated by these experiments that the algorithms can always perform successfully during the learning process and lead to a good clustering result.

Materials and Methods {#s2}
=====================

We will start with a brief review of markov random walks on complex network [@pone.0082964-Chorin1], [@pone.0082964-Chorin2]. Let be a network with *N* nodes and *M* edges, where *S* is the nodes set, is the weight matrix and is the weight for the edge connecting the nodes *x* and *y*. A simple example of the weight matrix is given by the adjacency matrix: or 1, depending whether *x* and *y* are connected. We can relate this network to a discrete-time Markov chain with stochastic matrix *p* with entries given bywhere is the degree of the node *x* [@pone.0082964-Lovsz1]--[@pone.0082964-Lafon1], [@pone.0082964-Chung1]. This Markov chain has stationary distributionand it satisfies the detailed balance condition

We take a partition of *S* as with if . Our aim is to aggregate the nodes in each cluster in order to give the exact expression of lumped markov chain. To do so, first, we regard each set in the state space as corresponding to the nodes of a -nodes network , where , and the weight on the edge that connects and is defined aswhere the sum involves all the transition probabilities between . From the detailed balance condition (4), it can be verified that . By setting , one can define the lumped Markov chain on graph with stationary distribution and lumped transition probabilitiesIt can be easily shown that is a stochastic matrix on the state space and satisfies a detailed balance condition with respect to , i.e.This construction of a lumped markov chain is shown in [Figure 1](#pone-0082964-g001){ref-type="fig"}.

![Lumped markov process of a network.\
For a given partition on a network , we define a lumped network by aggregating all nodes belonging to a subset into a meta-node. The new weights are computed via weight averaging the transition probabilities between nodes and --, and the lumped Markov chain with transition probabilities can also be obtained.](pone.0082964.g001){#pone-0082964-g001}

To give a more clear probability expression, we denote such probability function as to represent the probability which the node *x* belongs to the *k*-th cluster with. Naturally we need the assumption thatfor all .

Thus, the lumped markov transition matrix can be rewritten based on

Another approach of deriving the lumped markov transition matrix is based on lifted transition matrix. Any lumped stochastic matrix can be naturally lifted to the space of stochastic matrices on the original state space *S* viawhere

The basic idea in [@pone.0082964-Ma1], [@pone.0082964-Ma2] is to introduce a metric, also called the Hilbert-Schmidt norm, in the space of stochastic matrices. The optimal partition and the corresponding reduced Markov chain is found by minimizingwheresubject to the constraint [Eq.(7)](#pone.0082964.e042){ref-type="disp-formula"}.

To minimize the objective function *J* in [Eq.(12)](#pone.0082964.e051){ref-type="disp-formula"}, we definewhich has the similar form as in hard clustering. Then is indeed a stochastic matrix because for all *k*. Furthermore we have that satisfies the detailed balance condition with respect to , that isWith the above background, we have the following basic results.

The expressions of and according to [Eq.(12)](#pone.0082964.e051){ref-type="disp-formula"} are given by where is a matrix and is a matrix with entriesThe diagonal matrices , are and respectively, with entries where and are both Kronecker delta symbols.

The proof of Lemma 1 can be found in [Appendix_S1](#pone.0082964.s001){ref-type="supplementary-material"}.pdf. [Eq.(15a)](#pone.0082964.e061){ref-type="disp-formula"} and [Eq.(15b)](#pone.0082964.e062){ref-type="disp-formula"} give the partial derivatives of the objective function in [Eq.(12)](#pone.0082964.e051){ref-type="disp-formula"}, which are the critical points of constructing gradient methods.

Based upon this formulation, we can find the optimal for any fixed partition . From the optimality conditionwe can obtain [Eq.(8)](#pone.0082964.e045){ref-type="disp-formula"}. It indicates that when the partition is known, the minimizer of [Eq.(11)](#pone.0082964.e050){ref-type="disp-formula"} is unique, which can be given by [Eq.(8)](#pone.0082964.e045){ref-type="disp-formula"}, and the corresponding is the stochastic matrix in the class [Eq.(9)](#pone.0082964.e047){ref-type="disp-formula"} which provides the best rank- approximation of the original one under the metric [Eq.(11)](#pone.0082964.e050){ref-type="disp-formula"}.

SDP {#s2a}
---

In this section, based on the lumped markov transition matrix and corresponding results derived before, we devote to find ways to uncover partitions. An obvious choice to solve the constrained optimization [Eq.(11)](#pone.0082964.e050){ref-type="disp-formula"} is the steepest descent method [@pone.0082964-Hastie1]. However, the components of and may become negative and non-normalized during the descent procedure for our problems, which makes the probabilistic interpretation useless. To ensure the nonnegativity and normalization conditions for and , we add a projection step after each renewal. This leads to the following Steepest Descent method with Projections(SDP).

1.  Set up the initial state as the indicator matrix for each node in the network with the -means algorithm in [@pone.0082964-Weinan1], .

2.  Perform the following iteration until : Here is some type of projection operator which maps a real vector into a probability vector, is the learning rate and is a prescribed tolerance.

3.  The final gives the fuzzy partition for each node.

Two choices of the projection operator are used in our computations, but the final results seem to be insensitive to them. Now suppose , and when , we make projection as any of the following two choices

P1 Direct projection to the boundary.

When , we set ; otherwise we set .

P2 Iterative projections.

At first make projection to the hyperplane , then check each component of the projected . If , we take and make projection again to a dimension reduced hyperplane . Repeat the projection procedure to a lower and lower dimensional hyperplane until no component is negative.

The application of projection operators for SD can simply solve the original constrained optimization [Eq.(11)](#pone.0082964.e050){ref-type="disp-formula"}, while to solve it exactly may involve much more complexity. Note that with this SDP method, we can not guarantee that and are the exact minimizer of the original problem [Eq.(11)](#pone.0082964.e050){ref-type="disp-formula"} with non-negative and normalization constraints, but we take them as an approximate minimizer. The numerical results show that this strategy works fine for many examples.

The learning rate was usually chosen that started from a reasonable initial value and then reduced to zero with the iteration number in such a way that , andThe typical example of such case is , where is a positive constant [@pone.0082964-Ma3]. Another choice is to fix the learning rate as a positive constant [@pone.0082964-Ma1], [@pone.0082964-Ma2], which we utilize here, since the initial partition is good enough that the objective function [Eq.(12)](#pone.0082964.e051){ref-type="disp-formula"} descends much more slowly when the learning rate becomes smaller, while larger values of cause blow up.

The number of the iteration steps is difficult to be estimated, which may depend on the structure of the network itself, the choice of the initial values, etc. It usually converges fast for well-clustered networks and may converge slowly for diffusive networks.

Now let us estimate the computational cost in each iteration. In the iteration step for , all of the matrices are of order and full according to [Eq.(15a)](#pone.0082964.e061){ref-type="disp-formula"}. It is easy to find that the computation of costs , and the computation of costs . Note that the stochastic matrix is sparse with entries, so the computation for costs , where represents the number of edges, which is usually assumed in realistic networks. So finally, we obtain the cost in the step for is . The cost for is also according to [Eq.(15b)](#pone.0082964.e062){ref-type="disp-formula"}, since is involved in the equations.

CGP {#s2b}
---

Another choice is to minimize the objective function using a simplified formulation of traditional conjugate gradient method, which is frequently used in machine learning [@pone.0082964-Qian1]. It can be also regard as the above steepest descent method with a non-zero momentum term, which leads to the following reduced Conjugate Gradient method with Projections(CGP).

1.  Set up the initial state as the indicator matrix for each node in the network with the -means algorithm in [@pone.0082964-Lafon1], .

2.  Perform the following iteration until : Here is some type of projection operator which maps a real vector into a probability vector, is the learning rates and is a prescribed tolerance.

3.  The final gives the fuzzy partition for each node.

We again note that this is just a reduced form of conjugate gradient method, and it is demonstrated by simulation experiments that such method performs more efficiently than SD, just like the superiority traditional conjugate gradient has. The learning rate and are still chosen as constants by experience due to the same reason mentioned above. The computational cost in each iteration of CGP is the same as SDP, which is also for both and . Associating the two projections described above with SDP and CGP respectively, we refer to the derived algorithms: SDP1, SDP2, CGP1, CGP2, as the fuzzy partitioning algorithms for networks.

Results {#s3}
=======

In this section, simulation experiments on artificial network, the ad hoc network with 128 nodes, is carried out to demonstrate the performance of the proposed algorithms, via comparing the clustering results with some priori quantities. Moreover, the algorithms is applied to real-world network, the social interactions between members of a karate club.

Ad hoc network with 128 nodes {#s3a}
-----------------------------

The first example is the ad hoc network with 128 nodes. The ad hoc network is a benchmark problem used in many papers [@pone.0082964-Weinan1], [@pone.0082964-Newman3], [@pone.0082964-Newman4], [@pone.0082964-Danon1]. It has a known partition and is constructed as follows. Suppose we choose nodes, split them into 4 clusters with 32 nodes each. Assume that pairs of nodes belonging to the same clusters are linked with probability , and pairs belonging to different clusters with probability . These values are chosen so that the average node degree is fixed at . In other words, and are related asWe will denote .

To test on a less diffusive network, we take and generate the network according to [Eq.(1)](#pone.0082964.e006){ref-type="disp-formula"}. This network has a fuzzy clustering structure that some nodes should have immediate weights belonging to different clusters. We set the parameters by , and the learning rates in this model computation. is obtained after initialization [@pone.0082964-Hastie1]. The numerical results are shown in [Table 1](#pone-0082964-t001){ref-type="table"}. Here we compare with an interesting quantity, the degree fraction , which is defined aswhere is the number of nodes that are connected with and lie in cluster . Thus we have . With this definition, means the fraction of the edges connected with the node in the -th cluster. Note that this is not the same as the clustering probability, even though it is an interesting quantity to be compared with. We expect that the degree fraction [Eq.(2)](#pone.0082964.e008){ref-type="disp-formula"} is close to our result for network though generally this assertion needs to be justified or disconfirmed theoretically. To verify this fact, we define the mean and maximal -error of for error comparing. [Table 1](#pone-0082964-t001){ref-type="table"} shows that the deviation between these two is about . Obviously CG algorithm improves the convergence rate of SD. The projection P1 has the smallest , while the projection P2 reaches a better minimum which indicates a more accurate result. In [Figure 2](#pone-0082964-g002){ref-type="fig"} we plot the probability distribution function (pdf) of and . We observe that the shape of the pdf for or are almost the same. Note that all the \'s have a lower peak centered at about 0.85, which corresponds to the nodes in this cluster, and a higher peak centered at about 0.05, which corresponds to the other nodes outside of this cluster. The case for is similar but with the lower peak centered at about 0.7 and the higher peak centered at about 0.1. We note here that the center corresponds to the choice of the parameters . If we classify the nodes according to the majority rule, i.e. if for a given node then we set , we obtain the 4-cluster partition exactly for this model. This also verifies the accuracy of our algorithms, but fuzzy algorithms give more detailed information for each node.

![The pdf of and for the given ad hoc network with 128 nodes and .\
The solid lines and dashed lines represent the pdf of and respectively. In each figure, the lower peak corresponds to the nodes in this cluster, and the higher peak corresponds to the other nodes outside of this cluster.](pone.0082964.g002){#pone-0082964-g002}

10.1371/journal.pone.0082964.t001

###### .

![](pone.0082964.t001){#pone-0082964-t001-1}

               Iterations                    
  ----------- ------------ -------- -------- --------
  SDP1&CGP1      110&55     5.8794   0.1191   0.2389
  SDP2&CGP2      119&54     5.8757   0.1207   0.2374

The iterations, the value of the objective function and the mean and maximum -error of defined in [Eq.(3)](#pone.0082964.e009){ref-type="disp-formula"} for algorithms.

Karate club network {#s3b}
-------------------

This network was constructed by Wayne Zachary after he observed social interactions between members of a karate club at an American university [@pone.0082964-Zachary1]. Soon after, a dispute arose between the clubs administrator and main teacher and the club split into two smaller clubs. It has been used in several papers to test the algorithms for finding clusters in networks [@pone.0082964-Weinan1], [@pone.0082964-Girvan1]--[@pone.0082964-Newman5].

There are 34 nodes in karate club network, where each node represents one member in the club. In Zachary\'s original partition, each node belongs to only one sub-club after splitting. We label it as red or yellow color in the figures to show its attribute in the graph representation. From the viewpoint of the fuzzy clustering, the attribute of each node is no longer an indicator function but rather a discrete probability distribution. In our following notations, the association probability and means the probability of each node belonging to red or yellow colored cluster respectively.

We set the parameters by . Here is obtained after initialization [@pone.0082964-Weinan1]. The numerical results are shown in [Table 2](#pone-0082964-t002){ref-type="table"}. It can be demonstrate again that CGP is more efficient, and P2 can reach a smaller value of . [Figure 3](#pone-0082964-g003){ref-type="fig"} shows the convergence history during 5--30 iterations for the methods. It can be obviously seen that CGP, which decreased the objective function faster than SDP, performs more efficiently. The final association probabilities are presented in [Table 3](#pone-0082964-t003){ref-type="table"}, where and are the probability of belonging to the red or yellow colored group shown in [Figure 4](#pone-0082964-g004){ref-type="fig"} respectively. Comparing or between P1 and P2, we find that almost all the errors are less than , but the association probability is quite different from the 0--1 distributions obtained in the -means algorithm. Now let us compare the association probability and obtained by our methods with the original partition result obtained by Zachary. In [@pone.0082964-Zachary1], Zachary gave the partition and . If we classify the nodes according to the majority rule, i.e., if then we set , otherwise we set , we obtain the same partition as Zachary\'s (see [Figure 4(a)](#pone-0082964-g004){ref-type="fig"}). We note that this hard partition deduced by fuzzy partition is more reasonable than the result of *k*-means(see [Figure 4(b)](#pone-0082964-g004){ref-type="fig"}), since node 10 is classified correctly this time.

![The convergence history of the objective function during 5--30 iterations of the karate club network.](pone.0082964.g003){#pone-0082964-g003}

![The visualization of the partitioning results of the karete club network.\
(a)Hard partition with thresholding operation according to the node\'s maximal weight. (b)Hard partition by *k*-means in [@pone.0082964-Lafon1]; (c)Fuzzy partition by CGP2.](pone.0082964.g004){#pone-0082964-g004}

10.1371/journal.pone.0082964.t002

###### of the karate club.

![](pone.0082964.t002){#pone-0082964-t002-2}

               Iterations  
  ----------- ------------ --------
  SDP1&CGP1      153&55     4.0394
  SDP2&CGP2      173&53     4.0389

The iterations and minimized objective function values for the algorithms of the karate club network.

10.1371/journal.pone.0082964.t003

###### The association probability.

![](pone.0082964.t003){#pone-0082964-t003-3}

               Nodes     1       2       3       4     5   6   7     8       9       10     11   12
  ----------- ------- ------- ------- ------- ------- --- --- --- ------- -------- ------- ---- ----
  SDP1&CGP1            0.039   0.079   0.44      0     0   0   0   0.004   0.675    0.764   0    0
                       0.96    0.92    0.56      1     1   1   1   0.996   0.3254   0.236   1    1
  SDP2&CGP2            0.052   0.098   0.452   0.018   0   0   0   0.03    0.678    0.766   0    0
                       0.948   0.902   0.548   0.982   1   1   1   0.97    0.3222   0.234   1    1

               Nodes   13    14     15   16   17   18   19    20     21   22   23   24
  ----------- ------- ---- ------- ---- ---- ---- ---- ---- ------- ---- ---- ---- ----
  SDP1&CGP1            0    0.227   1    1    0    0    1    0.303   1    0    1    1
                       1    0.773   0    0    1    1    0    0.697   0    1    0    0
  SDP2&CGP2            0    0.245   1    1    0    0    1    0.317   1    0    1    1
                       1    0.755   0    0    1    1    0    0.683   0    1    0    0

               Nodes   25   26   27    28      29     30    31      32     33    34       
  ----------- ------- ---- ---- ---- ------- ------- ---- ------- ------- ---- ------- -- --
  SDP1&CGP1            1    1    1    0.965   0.858   1    0.734   0.91    1    0.963     
                       0    0    0    0.035   0.142   0    0.266   0.09    0    0.037     
  SDP2&CGP2            1    1    1    0.963   0.856   1    0.737   0.907   1    0.96      
                       0    0    0    0.037   0.144   0    0.263   0.093   0    0.04      

The association probability that each node belongs to different clusters. and means the probability belonging to red or yellow colored cluster in [Figure 4](#pone-0082964-g004){ref-type="fig"} respectively.

However, we have more detailed information in fact. From [Table 3](#pone-0082964-t003){ref-type="table"}, we find for nodes , which lie at the boundary of the yellow colored group; and for nodes , which mostly lie at the boundary of the red colored group. The others belong to the red and yellow colored groups with nonzero probability, especially the nodes have more diffusive probability and they play the role of transition nodes between the red and yellow colored groups. We can visualize the data more transparently with the color vector [Eq.(23)](#pone.0082964.e166){ref-type="disp-formula"} for each node. We can conclude from [Figure 4(c)](#pone-0082964-g004){ref-type="fig"} that how much probability each of the 34 members stands by both parts with. Members and have an obvious attitude on following their leader, i.e. the club administrator or the main teacher. Others such as hold neutralism that they can support either leader according their weights.

Discussion {#s4}
==========

We address the expression of lumped markov transition matrix for networks with two novel method in this paper. This can also be considered as a generalization of markov random walk dynamic in statistics for the networks. We successfully constructed the steepest descent method with projection (SDP) and the reduced conjugate gradient method with projection (CGP). They are derived to search for the local minimum of the objective function in [Eq.(12)](#pone.0082964.e051){ref-type="disp-formula"} under the fuzzy clustering framework, which is extended from a deterministic framework for network partition based on the optimal prediction of a random walker Markovian dynamics [@pone.0082964-Weinan1]. The simulation experiments have shown that the algorithms can efficiently determine the fuzzy partition matrix. Partitioning the network with thresholding operation according to the node\'s maximal weight can give a more reasonable clustering result than the previous *k*-means algorithm [@pone.0082964-Weinan1]. We use two datasets(Ad hoc network with 128 nodes and Karate club network) to validate algorithms and achieve good results Numerical results show that our algorithms with two different projections produce similar results, while the CGP2 algorithm has better efficiency and accuracy. Moreover, the algorithms succeed in real-world learning tasks.

Supporting Information {#s5}
======================

###### 

**The ad hoc network is a benchmark problem used in many papers** [@pone.0082964-Chorin1] **,** [@pone.0082964-Chorin2] **,** [@pone.0082964-Hastie1] **,** [@pone.0082964-Newman3] **.** It has a known partition and is constructed as follows. Suppose we choose *N* = 128 nodes, split them into 4 clusters with 32 nodes each.

(TXT)

###### 

Click here for additional data file.

###### 

**The file karate.txt contains the network of friendships between the 34 members of a karate club at a US university, as described by Wayne Zachary in 1977.** If you use these data in your work, please cite W. W. Zachary, An information flow model for conflict and fission in small groups, Journal of Anthropological Research 33, 452--473 (1977).

(TXT)

###### 

Click here for additional data file.

###### 

**The Proof of Lemma 1.**

(PDF)

###### 

Click here for additional data file.

[^1]: **Competing Interests:**The authors have declared that no competing interests exist.

[^2]: Conceived and designed the experiments: FJ XJ. Performed the experiments: FJ XJ. Analyzed the data: WJ QJ. Contributed reagents/materials/analysis tools: FJ XJ. Wrote the paper: FJ XJ WJ QJ. Ad hoc network with 128 nodes: WJ. Karate club network: OJ. Wrote Introduction and Result: FJ. Wrote Discussion and Materials and Methods: XJ. Wrote Author Summary: WJ Wrote Discussion: QJ.
