Jason Wyse

Nial Friel

Pierre Latouche

Inferring structure in bipartite networks using the latent blockmodel and exact ICL

Keywords: Bipartite network, network clustering, bi-clustering, blockmodel, model selection, integrated complete likelihood

come

Inferring structure in bipartite networks using the latent blockmodel and exact ICL

Introduction

Bipartite networks are those containing two types of nodes, say types A and B. Nodes of type A may be linked to nodes of type B, and vice versa, but links between two nodes of the same type are not considered. There are many real life networks that can be naturally viewed in this way. Take for example relational networks where a user rates a movie. A user is a member of node type A and node type B represents the movies. One may ask a number of quantitative questions in such a situation. Can users be grouped by the types of movies they watch and rate? How many substantive genres of movies are defined by users?

Bipartite or two-mode networks have seen much attention in the social networks and machine learning literature, see for example [START_REF] Borgatti | Network analysis of 2-mode data[END_REF], [START_REF] Doreian | Generalized blockmodeling of two-mode network data[END_REF], [START_REF] Brusco | A tabu search heuristic for deterministic two-mode blockmodelling of binary network matrices[END_REF], [START_REF] Brusco | A variable neighbourhood search method for a two-mode blockmodelling problem in social network analysis[END_REF], [START_REF] Doreian | Partitioning large signed two-mode networks: Problems and prospects[END_REF]. [START_REF] Rohe | Co-clustering directed graphs to discover asymmetries and directional communities[END_REF] discuss how their stochastic co-blockmodel may be extended to a bipartite setting. The reason for this high level of interest is due to their wide applicability and the fact that it is often 1 very natural and fruitful to model interactions between node sets in this way. Clustering or partitioning the node sets simultaneously can reveal structure and give considerable insight into the entities in the network. Such insights may be allusive to the more classical network measures [START_REF] Wasserman | Social Network Analysis: Methods and Applications[END_REF], some of which have been adapted from the classical literature to the bipartite or two mode situation (for example, the clustering coefficients of [START_REF] Opsahl | Triadic closure in two-mode networks: Redefining the global and local clustering coefficients[END_REF]). For this reason much attention is focused on clustering or grouping the node sets in tandem; this practice is referred to using many terms in the literature: bi-clustering, co-clustering, block-clustering, two-mode blockmodelling and others. One of the pioneering papers of this area was that of [START_REF] Hartigan | Direct Clustering of a Data Matrix[END_REF]. [START_REF] Flynn | Consistent Biclustering[END_REF] examine asymptotic theoretical guarantees for bi-clustering.

Approaches to clustering or blockmodelling of two-mode networks fall into two classes. The first of these is often called deterministic, whereby the clustering of the network (or adjacency matrix) is obtained by minimizing an objective function which measures discrepancy from an ideal block structure. Examples of this are the work of [START_REF] Doreian | Generalized blockmodeling of two-mode network data[END_REF], [START_REF] Brusco | Inducing a blockmodel structure on two-mode data using seriation procedures[END_REF], [START_REF] Brusco | A tabu search heuristic for deterministic two-mode blockmodelling of binary network matrices[END_REF], [START_REF] Brusco | A variable neighbourhood search method for a two-mode blockmodelling problem in social network analysis[END_REF] and [START_REF] Doreian | Partitioning large signed two-mode networks: Problems and prospects[END_REF]. The second type of approach is stochastic. In stochastic blockmodelling procedures, one assumes that the probability of links between the node sets in the network can be modelled by a parameterized distribution. These parameters are usually estimated (learned) and then used as a representative embodiment of the true network linking behaviour. The stochastic approach may also be referred to as model-based; that is, a statistical model is used for links in the network. Examples include the work of [START_REF] Govaert | Simultaneous clustering of rows and columns[END_REF], [START_REF] Govaert | Comparison of the mixture and the classification maximum likelihood in cluster analysis when data are binary[END_REF], [START_REF] Govaert | Clustering with block mixture models[END_REF], [START_REF] Govaert | An em algorithm for the block mixture model[END_REF], [START_REF] Govaert | Block clustering with Bernoulli mixture models: Comparison of different approaches[END_REF], [START_REF] Rohe | Co-clustering directed graphs to discover asymmetries and directional communities[END_REF], [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] and [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF].

This paper is concerned with the stochastic approach to blockmodelling of two-mode networks. In particular, we take the latent blockmodel (LBM) developed in a series of papers by Gérard Govaert and Mohamed Nadif;[START_REF] Govaert | Simultaneous clustering of rows and columns[END_REF], [START_REF] Govaert | Comparison of the mixture and the classification maximum likelihood in cluster analysis when data are binary[END_REF], [START_REF] Govaert | Clustering with block mixture models[END_REF], [START_REF] Govaert | An em algorithm for the block mixture model[END_REF], [START_REF] Govaert | Block clustering with Bernoulli mixture models: Comparison of different approaches[END_REF]. We consider this model as applied in the context of bipartite networks. This is a desirable model, as it provides a model-based clustering of both node sets and has richer modelling capability than only absence/presence data for ties between nodes should this information be observed. The LBM is based around an intuitive generative structure as outlined in Section 2.2. We point out that the LBM is a different model than the stochastic blockmodel (SBM), as in, for example, [START_REF] Nowicki | Estimation and Prediction for Stochastic Blockstructures[END_REF]. The LBM operates on items and objects as opposed to the SBM which focusses on modelling interactions between items.

In the LBM, the posterior distribution over the latent label vectors, given the model parameters and observed data, cannot be factorized due to conditional dependency. Therefore standard optimization techniques such as the expectation maximization (EM) algorithm cannot be used directly for clustering. To tackle this issue, approximation methods like variational EM [START_REF] Govaert | Block clustering with Bernoulli mixture models: Comparison of different approaches[END_REF] or stochastic EM [START_REF] Keribin | Estimation d'un modèle à blocs latents par l'algorithme SEM[END_REF]) have been proposed. Moreover, in practice, the numbers of clusters in each node set have to be estimated. This has led to treatements such as the one by [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] who use Markov chain Monte Carlo (MCMC) to do inference for the number of clusters and the members of the nodes to the groups. We also refer to the recent work of [START_REF] Keribin | Model selection for the binary latent block model[END_REF] and [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF] who relied on model selection criteria to estimate the number of clusters.

Unlike [START_REF] Govaert | Block clustering with Bernoulli mixture models: Comparison of different approaches[END_REF], [START_REF] Keribin | Model selection for the binary latent block model[END_REF], [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF], our ap-proach allows the number of clusters in both node sets to be estimated while simultaneously partitioning the nodes. This is based on a clustering criterion termed the exact integrated complete likelihood (ICL), and a method to search over partitions of the nodes. The main ideas of using the exact ICL come from [START_REF] Côme | Model selection and clustering in stochastic block models based on the exact integrated complete data likelihood[END_REF] who use this in estimation of the SBM of [START_REF] Nowicki | Estimation and Prediction for Stochastic Blockstructures[END_REF]. Our work can be seen as somewhat complementary to the work in [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF], but there are many advantages to using the framework presented here, similar to the SBM in [START_REF] Côme | Model selection and clustering in stochastic block models based on the exact integrated complete data likelihood[END_REF]. Firstly, it is more scalable than the MCMC approach and secondly we do not have to worry about the mixing rates of the MCMC algorithm in larger settings. One drawback of our new approach is that it does not provide a joint posterior distribution for the number of clusters in both node sets. The remainder of the paper is organised as follows. Section 2 introduces ideas of blockmodels for bipartite networks (Section 2.1) and the LBM of [START_REF] Govaert | Block clustering with Bernoulli mixture models: Comparison of different approaches[END_REF] (Section 2.2) and discusses its relevance for this context. Section 2.3 reviews existing estimation techniques for the LBM and discusses their advantages and limitations. In Section 3 we introduce the exact ICL along with a greedy algorithm for inference purposes. We also draw strong parallels with the approach of [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] in this section, and highlight the sensible and pragmatic nature of the algorithm. The computational complexity of the algorithm is discussed in Section 3.4. In Section 3.6 we show how large computational savings can be made over a naive implementation of this algorithm. Section 4 compares the proposed approach with the competing MCMC algorithm through a simulation study. Section 5 applies the approach to a number of datasets. We conclude with a discussion.

2 Bipartite networks and the latent blockmodel

Bipartite networks

Bipartite networks consist of possible ties between members of two different node sets. Let the node sets be A and B. There are N nodes in A, A = {a 1 , . . . , a N } and M in B, B = {b 1 , . . . , b M }. We use the terms node type and set interchangeably. Node i from A may have a tie with node j from B, this tie, if it exists, is undirected. Alternatively, a tie may not exist. Either way, this information on the linking attribute between the two nodes is contained in the observation y ij . We record all the observed information about the network in an N × M adjacency matrix

Y =    y 11 . . . y 1M . . . . . . . . . y N 1 . . . y N M   
having a row for each node from A and a column for each in B. The problem of primary focus for us is the following. Can we group or partition the A nodes and the B nodes simultaneously to reveal subgroups or subsets of the A nodes that have linking attributes of similar nature to subgroups of the B nodes? We use the term "similar nature" here to highlight that we approach this problem by modelling properties of the linking attribute which could be binary in nature, categorical (> 2 categories), a count or an observation with a continuous value.

Suppose a grouping of the node sets into subsets like that mentioned above does exist. The nodes in A are partitioned into K and those in B into G disjoint subsets

A = K k=1 A k , A k ∩ A l = ∅, ∀k = l and B = G g=1 B g , B g ∩ B h = ∅, ∀g = h.
Then the nodes in set A k will be seen as having similar natured linking attributes to their respective nodes in set B g . We suppose now that we model linking attributes using some parametric distribution p(y ij |θ) for the y ij . In order to capture differences in the nature of linking attributes between the different subsets, we allow this distribution to have a different valued parameter θ kg for each subset pairing A k and B g , so that if node i is in A k and node j is in B g , the density y ij is p(y ij |θ kg ).

Latent blockmodel

The LBM was introduced in [START_REF] Govaert | Block clustering with Bernoulli mixture models: Comparison of different approaches[END_REF] as a means to provide a concise summary of a large data matrix. In order to index the different possible partitionings of A and B into K and G disjoint subsets (clusters) we introduce respective label vectors c and w. These are such that c i = k if node i ∈ A k and similarly w j = g if j ∈ B g . Using the model for linking attributes, the probability of observing the adjacency Y can be written down, if we assume we know the partitioning of the nodes, i.e. that know c and w

p(Y |c, w, Θ, K, G) = K k=1 G g=1 i:c i =k j:w j =g p(y ij |θ kg ).
Here Θ denotes the collection of θ kg .

Of course, in practice c and w will not be known. Attach to each clustering of nodes c, w a probability p(c, w|Φ, K, G), depending on some hyperparameters Φ. This allows the density of the adjacency matrix Y to be written down as a mixture model. Let C K and W G be collections containing all possible clusterings of A into K groups and B into G groups respectively. Then

p(Y |Θ, Φ, K, G) = (c,w)∈C K ×W G p(c, w|Φ, K, G)p(Y |c, w, Θ, K, G).
(1)

Here the mixture is over all clusterings of the nodes into K and G groups. The p(c, w|Φ, K, G) terms represent probabilities of particular clusterings being generated parameterized by Φ. Suppose that, a priori, no information exists on the joint clustering of nodes in sets A and B. Then it is reasonable to assume that p(c, w|Φ, K, G) = p(c|Φ, K)p(w|Φ, G).

Having made this assumption, LBM assumes further that there are weights associated with each subset A k and B g , such that Pr{c i = k|K} = ω k and Pr{w j = g|G} = ρ g . This defines a multinomial distribution for the node labels, where the weights sum to unity:

K k=1 ω k = 1, G g=1 ρ g = 1.
The parameters Φ represent the weight vectors ω, ρ. Thus we can write

p(c, w|Φ, K, G) = K k=1 ω N k k G g=1 ρ Mg g where N k = |A k | and M g = |B g |.
Effectively the LBM defines a probability distribution over clustering of the node sets A and B into K and G disjoint subsets. Examining equation (1) one sees that the sum is over K N G M terms, which is clearly intractable for even moderate N and M . In practice, we are interested in finding optimal or near optimal clusterings, corresponding to large values of p(c, w|Φ, K, G)p(Y |c, w, Θ, K, G) in (1) which contribute the most to the sum. It is thus usual to rephrase the problem into one which optimises p(c, w|Φ, K, G)p(Y |c, w, Θ, K, G) jointly over the clusterings (c, w) and parameters Φ, Θ in some way. Something to note here is that the particular cluster configurations are conditional on the value of K and G being known, something which will rarely be the case in real applications.

Estimation of the LBM

Estimation of the LBM can be performed either through hybrid variational and stochastic EM type algorithms [START_REF] Govaert | Block clustering with Bernoulli mixture models: Comparison of different approaches[END_REF][START_REF] Keribin | Estimation d'un modèle à blocs latents par l'algorithme SEM[END_REF] or through Bayesian estimation [START_REF] Van Dijk | A Bayesian approach to two-mode clustering[END_REF], [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF]). We refer the reader to Section 2.1.1 of [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] for a review of the former. Here we focus on Bayesian estimation as it is most relevant for what follows. The one remark to be made about EM type estimation of the LBM however is that it conditions on the values of K and G. As of yet, there is no widely accepted information criterion for choosing the best values of K and G as outlined in Section 2.1.2 of [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF]. [START_REF] Keribin | Model selection for the binary latent block model[END_REF], [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF] relied on a ICL criterion for model selection purposes, while [START_REF] Van Dijk | A Bayesian approach to two-mode clustering[END_REF] proposed a Gibbs sampler for Bayesian estimation of the latent block model along with an AIC-3 [START_REF] Bozdogan | Mixture-model cluster analysis using model selection criteria and a new informational measure of complexity[END_REF] criterion.

Wyse & Friel (2012) circumvented the problem of choosing K and G by including these as unknowns in their Bayesian formulation of LBM. They begin by taking priors on the unknowns in the model, namely, Θ, Φ, K and G and writing down the full posterior given the data

π(K, G, c, w, Φ, Θ|Y ) ∝ p(c, w|Φ, K, G)p(Y |c, w, Θ, K, G) ×π(Θ|K, G)π(Φ|K, G)π(K)π(G) (2) = p(c, w|Φ, K, G)π(Φ|K, G) ×p(Y |c, w, Θ, K, G)π(Θ|K, G)π(K)π(G)
where the last three distributions on the right hand side are priors for Θ, K, G. Note here the assumption that, a priori, the number of node clusters K and G of types A and B are independent. If other prior information exists beforehand, this prior assumption can be replaced by one which represents this information. Integrating both sides of the above proportionality relation with respect to Φ and Θ returns the joint marginal distribution of K, G, c, w given the observed network Y :

π(K, G, c, w|Y ) ∝ p(c, w|Φ, K, G)π(Φ|K, G) dΦ × p(Y |c, w, Θ, K, G)π(Θ|K, G) dΘ ×π(K) π(G) = π(c, w|K, G)π(Y |c, w, K, G)π(K)π(G).
The key observation in [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] is that the quantities

π(c, w|K, G) = p(c, w|Φ, K, G)π(Φ|K, G) dΦ π(Y |c, w, K, G) = p(Y |c, w, Θ, K, G)π(Θ|K, G) dΘ
can be obtained exactly (analytically) using relatively standard prior assumptions. For example π(c, w|K, G) can be computed exactly by assuming independent Dirichlet priors on the weights vectors

π(Φ|K, G) = Dir K (ω; α 0 , . . . , α 0 ) × Dir G (ρ; β 0 , . . . , β 0 )
where Dir p represents the density of the Dirichlet distribution on the p -1 dimensional simplex. In this case

π(c, w|K, G) = Γ{α 0 K} Γ{α 0 } K K k=1 Γ{N k + α 0 } Γ{N + α 0 K} × Γ{β 0 G} Γ{β 0 } G G g=1 Γ{M g + β 0 } Γ{M + β 0 G}
where Γ(•) is the gamma function. See [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] for further details. Furthermore π(Y |c, w, K, G) can be computed exactly if we make the following assumptions. Assume that the prior for Θ can be expressed as independent priors for the θ kg ,

π(Θ|K, G) = K k=1 G g=1 π(θ kg ). We can compute π(Y |c, w, K, G) exactly if we can compute Λ kg = π(θ kg ) i:c i =k j:w j =g p(y ij |θ kg ) dθ kg
exactly. This happens when the prior π(θ kg ) is fully conjugate to p(y|θ kg ). There are many widely used and standard situations where this is the case as we discuss in Section 3.2. In our notation we suppress dependence on the values of the hyperparameters, say α 0 , β 0 and the hyperparameters of π(θ kg ) and similar for brevity. This is to be understood.

With these quantities available exactly, the joint marginal posterior of clusterings and number of clusters of the nodes is given by

π(K, G, c, w|Y ) ∝ π(c, w|K, G) K k=1 G g=1 Λ kg π(K)π(G).
(3)

Notice that this posterior is defined over a large discrete model space. If we allow a maximum of K max and G max subsets for nodes in A and B, then the size of the support of the posterior is Kmax k=1 Gmax g=1 k N g M . [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] developed an MCMC algorithm that generates samples from (3). The algorithm moves are involved and we refer to [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] for exact details. Node labels are sampled using a Gibbs step and there are more sophisticated moves for adding and removing clusters (that is, changing K and G). The main idea is that iteratively the chain will sample high probability configurations (K, G, c, w), so that in the sampler output one gets a joint posterior distribution for K and G with corresponding label vectors. There are upper bounds K max and G max placed on the number of groups the nodes may be partitioned into. Usually this will be conservatively large. A drawback of this algorithm in general is that mixing of the Markov chain may disimprove with increasing M and/or N , resulting in infrequent jumps to models with different K and G from the current state. So for high dimensional problems, one can observe most of the empirical (approximate) posterior mass centred on one or two (K, G) combinations. However, the work of [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] is a step forward in estimation of LBMs as it is an automatic way to perform inference for K and G also, while clustering the nodes, for which no other such approach existed before.

3 Exact ICL and greedy ICL algorithm for bipartite networks

ICL and exact ICL

The integrated completed likelihood (ICL) was first introduced in [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF], as a model selection criterion, in the context of Gaussian mixture models. The rationale for using ICL is that in the finite mixture model context (similar to our situation here), analysis is often carried out using a latent label vector which it is difficult to integrate from the model. This latent vector is often termed the allocation vector [START_REF] Nobile | Bayesian finite mixtures with an unknown number of components: The allocation sampler[END_REF] and it provides a clustering of the data points to component densities. For this reason [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] argue that the evidence for a particular clustering should be taken into account when determining the number of mixture components and thus they focus on the integrated completed data likelihood. So instead of marginalising these labels from the model, ICL includes them as part of the information criterion. The number of components in the mixture (in our case values of K, G) which gives the largest ICL is the most supported by the data. For full details see [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF] have used ICL in the context of the latent blockmodel to choose the number of clusters, however, it is not the same approach as we adopt here as they do not use the available form of the exact ICL. Instead, their analysis employs a penalty term in the ICL criterion, which we avoid by computing ICL analytically. Using the notation introduced earlier, the ICL can be written as

ICL = log p(c, w|Φ, K, G)p(Y |c, w, K, G)π(Θ|K, G)π(Φ|K, G) dΘ dΦ .
Typically this quantity cannot be computed exactly and so ICL is usually approximated using a high probability configuration and a penalty term. In Section 2.3 we encountered a similar expression to ICL in the joint collapsed posterior of [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF], and we see that we can write

ICL = log p(c, w|Φ, K, G)p(Y |c, w, K, G)π(Θ|K, G)π(Φ|K, G) dΘ dΦ = log p(c, w|Φ, K, G)π(Φ|K, G)dΦ × p(Y |c, w, K, G)π(Θ|K, G)dΘ = log p(c, w|Φ, K, G)π(Φ|K, G)dΦ + log p(Y |c, w, K, G)π(Θ|K, G)dΘ = log π(c, w|K, G) + log π(Y |c, w, K, G).
The conditions for the ICL to be available exactly are the same as those given in Section 2.3. If this is the case we term this quantity exact ICL. Now suppose we can compute the exact ICL so that

ICL = log π(c, w|K, G) + K k=1 G g=1 log Λ kg .
Note the similarity between this expression and the log of the right hand side of (3). In fact, the only difference is that in (3) we also have two log prior terms for K and G. Thus, finding the highest ICL is almost identical to finding regions of the support of (3) with high posterior mass (bar the prior terms). This is elucidated more in Section 3.5.

Linking attribute models

We now turn attention to the types of linking attribute models for which exact ICL can be computed. The requirement for the exact ICL to be available is that

Λ kg = π(θ kg ) i:c i =k j:w j =g p(y ij |θ kg ) dθ kg
can be computed analytically. Λ kg can be referred to as the marginal likelihood for block (k, g), as we pick up nodes from cluster k of A and cluster g of set B and average their likelihood over the prior for the parameter of the linking attributed between the two node subsets. We now give some examples of very commonly used models where the marginal likelihood can be computed exactly.

Absent/present linking attribute

The most commonly observed networks will be those where edges are represented by binary indicators when there is a link between a node in set A and a node in set B or not. For example, students being members of university clubs or not. For i ∈ A k and j ∈ B g the most natural way to model this is using a Bernoulli random variable with probability of a tie equal to θ kg p(y ij |θ kg ) = θ

y ij kg (1 -θ kg ) 1-y ij .
A conjugate prior for this data distribution is a beta prior

θ kg ∼ Beta(η, η)
with hyperparameter κ. Then it can be easily shown that

Λ kg = Γ{2η} Γ{η}Γ{η} Γ{N 1 kg + η}Γ{N k M g -N 1 kg + η} Γ{N k M g + 2η}
where N 1 kg = i:c i =k j:w j =g I(y ij = 1).

Multinomial links with Dirichlet prior

This is a generalization of the previous model for more than two categories. Linking attributes which are categorical arising from C categories could naturally be modelled using the multinomial distribution. This could be useful say if the college years of the students in the university clubs above were also recorded e.g. junior freshman, senior freshman, junior sophister, senior sophister. The distribution is

p(y ij |θ kg ) = C l=1 θ l kg I(y ij =l)
where θ l kg is the probability y ij takes category l. Taking a symmetric Dirichlet prior with parameter ζ the block integrated likelihood can be computed exactly as

Λ kg = Γ(ζC) Γ(ζ) C C l=1 Γ(N l kg + ζ) Γ(N k M g + ζC)
where N l kg = i:c i =k j:w j =g I(y ij = l) is the number taking on category l in block (k, g).

Poisson model for count links with Gamma prior

Count data may be modelled using a Poisson distribution with rate θ kg . This could be the number of emails exchanged between students in the university clubs. Taking a Gamma(δ, γ) prior on θ kg leads to

Λ kg = γ δ Γ(δ) Γ(S kg + δ) (N k M g + γ) S kg +δ   i:c i =k j:w j =g y ij !   -1
where S kg = i:c i =k j:w j =g y ij is the sum of the counts in a block.

Gaussian model with Gaussian-Gamma prior

For continuous data, if appropriate one can assume a Gaussian distribution with mean µ kg and precision τ kg for observations in a block. This may be the number of hours spent by the university students on club activities. For example, one would expect a sports club to require greater time commitment than, say, a car appreciation club. Taking a N(ξ, {κτ kg } -1 ) prior for µ kg conditional on τ kg and a Gamma(γ/2, δ/2) prior on τ kg

Λ kg = π -N k Mg/2 κ 1/2 δ γ/2 (N k M g + κ) 1/2 Γ({N k M g + γ}/2) Γ(γ/2) SS kg + κξ 2 - (S kg + κξ) 2 N k M g + κ + δ -(N k Mg+γ)/2
.

Here S kg = i:c i =k j:w j =g y ij and SS kg = i:c i =k j:w j =g y 2 ij .

A greedy search strategy

Based on the observations about the ICL in Section 3.1 it makes sense to try to find some strategy to optimise the exact ICL. There are in essence four unknowns in our problem K, G, c, w. From Section 3.1 the largest value of the exact ICL (assuming it can be computed), gives the clustering of the nodes into subsets which is most preferable. The idea used by [START_REF] Côme | Model selection and clustering in stochastic block models based on the exact integrated complete data likelihood[END_REF] for the stochastic blockmodel was to optimise the ICL criterion using a greedy search over labels and the number of node clusters. We adopt a similar approach here. First define the ICL using any instance of numbers of clusters and labels

ICL(K, G, c, w) = log π(c, w|K, G) + K k=1 G g=1 log Λ kg .
This frames the exact ICL as an objective function in parameters (K, G, c, w) which we wish to optimise. Of course, the values of K and G place constraints on the possible values of c and w. The exact ICL is optimised iteratively by cycling through the following smaller optimisation techniques repeatedly until no further increases in it can be obtained.

Update node labels in set A

To update the node labels for set A, we firstly shuffle the order of nodes in the update. If node i is currently in A k (cluster k ), that is c i = k , compute the change in exact ICL when moving node i to A l for all l = k . First of all, suppose that currently |A k | > 1, then the change in exact ICL is given by

∆ k →l = ICL(K, G, c * , w) -ICL(K, G, c, w)
where c * is such that c * r = c r , r = i and c * i = l. Clearly ∆ k →k = 0. If all of the ∆ k →l < 0, do not move node i from subset k . Otherwise, move i to the node cluster A l * where l * = arg max l ∆ k →l . This process is illustrated in a flow diagram in Figure 1 which summarizes the primary parts of the update.

If |A k | = 1, then our algorithm assumes that moving node i from A k causes that subset to vanish, so that K, the number of subsets is reduced by 1. This is similar to component absorption in the case of algorithms for finite mixture models. In this case, the computation of the change in exact ICL must be modified accordingly to

∆ k →l = ICL(K -1, G, c * , w) -ICL(K, G, c, w).
This updating process is repeated for each node i in A.

Update node labels in set B

The node labels in set B are updated in an analogous way to those in set A, so we omit a description for brevity.

Algorithm initialization and termination

The algorithm is initialized randomly, with a random assignment of nodes to the subsets A k , B g , k = 1 . . . , K, g = 1, . . . , G. Initially we assume a large number of clusters of nodes. To draw analogy with the MCMC sampler described in Section 2.3 we set K = K max and G = G max at the beginning. Nodes in sets A and B are then processed, and the algorithm terminates when no improvement (or further increase) in the exact ICL can be obtained. Clearly, the algorithm will cause convergence of the ICL to a local maximum, which may not necessarily be the global maximum. Thus, the greedy algorithm is run a number of times (this could be done in parallel), and the run given the highest exact ICL value used. As suggested by one reviewer, in practice, using the algorithm of [START_REF] Rohe | Co-clustering directed graphs to discover asymmetries and directional communities[END_REF] (for the bipartite setting) to provide an initial clustering before invoking the greedy search algorithm could well lead to faster convergence of the algorithm and quicker run times. We have found this can be a useful strategy (see Section 5.1). The code for the algorithm of [START_REF] Rohe | Co-clustering directed graphs to discover asymmetries and directional communities[END_REF] is available from (http://pages.stat.wisc.edu/∼karlrohe/homepage/Welcome.html). We have included this algorithm as an initializer option in the implementation of our algorithm for the absent/present linking attribute.

Using a merging move at termination

At termination of the algorithm, we can attempt to merge clusters if this increases the exact ICL further. Merging the clusters is more of a "mass node" move than the one node updates used in the greedy search. To merge two A node clusters, we compute ICL(K -1, G, c, w) where if attempting to merge clusters k and k , c is such that all nodes with label k in c become label k nodes. This necessitates G new block marginal likelihood calculations to get

ICL(K-1, G, c, w)-ICL(K, G, c, w) = log π(c, w|K -1, G) π(c, w|K, G) + G g=1 log Λkg -log(Λ kg Λ k g )
where Λkg is obtained by merging the sufficient statistics for blocks (k, g) and (k , g) and recomputing the marginal likelihood for these new statistics. If this difference is greater than 0 we merge the clusters k and k , otherwise no merge is performed. All pairwise merges are considered after the termination of the greedy exact ICL algorithm. shows gives an overview of each stage of the process. The bottom portion outlines the computations carried out.

Compute

Computational complexity of the greedy algorithm

Quantifying exactly the computational complexity of the greedy algorithm is not possible, as it will be problem dependent. However, we can analyse its various components and give some guidelines and comparison with the MCMC algorithm of [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF]. Assume that the average cost of computing one of the Λ kg terms is Q Λ (this clearly depends on the model chosen for the linking attributes). To update node i in cluster A k we must compute K -1 terms of the form

∆ k →l = ICL(K, G, c * , w) -ICL(K, G, c, w) = log π(c * , w|K, G) + K k=1 G g=1 log Λ * kg -log π(c, w|K, G) + K k=1 G g=1 log Λ kg assuming that |A k | > 1,
where the Λ * kg refer to the label vector c * . One can see that Λ * kg = Λ kg for all but the affected subsets k and l, so that cancellation results in

∆ k →l = log π(c * , w|K, G) + G g=1 log Λ (-i) k g + log Λ (+i) lg -log π(c, w|K, G) + G g=1 (log Λ k g + log Λ lg ) . (4)
Where now we have substituted the Λ * kg for more meaningful and specific notation ("+i" meaning add i to the block, "-i" meaning take i out of the block). For all of the linking attribute models mentioned in Section 3.2, one can store the sufficient statistics of each (k, g) block to compute the Λ kg , meaning that Λ * kg can be computed by a slight modification of these statistics. Also we can store the Λ kg from the most recent state of the algorithm. The term log π(c * , w|K, G) -log π(c, w|K, G) is equal to

log Γ{N k -1 + α 0 }Γ{N l + 1 + α 0 } Γ{N k + α 0 }Γ{N l + α 0 } .
We assume that the cost of computing this quantity is dominated by the computation of the Λ kg . This is reasonable for all of the models in Section 3.2. This means that (4) can be computed with O(2Q Λ G) operations (for the case

|A k | = 1 it can be shown that O(Q Λ G)
operations are required for the calculation). So updating the label of node i is O(2Q Λ KG).

One can show that updating nodes in set B has the same cost, giving O(2N M Q Λ KG) overall per sweep of the greedy algorithm. If we compare to [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF], the total cost of the Gibbs update for nodes for T it iterations is fixed at O(2T it N M Q Λ KG). Usually we would expect convergence of the greedy algorithm in much fewer than T it iterations. Also there are other moves in [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF], which are costly but whose complexity is difficult to quantify exactly. In short, we expect our approach to be much more scalable than the MCMC based inference. Another point to note is that the performance in terms of mixing of the algorithm of [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] usually declines as the dimensions N and/or M increase, so that although we may be expending much more computing power, we may not necessarily be gaining much more information about the clustering. The MCMC algorithm also has limitations in size, in that to run it on larger datasets may not even be feasible from a time consideration.

Why the greedy search is both sensible and pragmatic

Consider a Gibbs sampler for updating the labels of the nodes from set A for example similar to that used in [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF]. Using (3) the full conditional distribution of the new label for node i with

c i = k , is π(c * i = l|c -i , w, Y, K, G) ∝ Γ{N k -1 + α 0 }Γ{N l + 1 + α 0 } G g=1 Λ (-i) k g Λ (+i) lg
for l = k where c -i denotes c without the i th entry and Λ

(-i)
k g denotes the quantity computed by removing the statistics associated with node i. Also

π(c * i = k |c -i , w, Y, K, G) ∝ Γ{N k + α 0 }Γ{N l + α 0 } G g=1 Λ k g Λ lg reflects no change in label. Notice that π(c * i = l|c -i , w, Y, K, G) π(c * i = k |c -i , w, Y, K, G) ∝ Γ{N k -1 + α 0 }Γ{N l + 1 + α 0 } Γ{N k + α 0 }Γ{N l + α 0 } G g=1 Λ (-i) k g Λ (+i) lg G g=1 Λ k g Λ lg = exp{∆ k →l } from (4). Using this relation, it is straightforward to show that π(c * i = l|c -i , w, Y, K, G) = exp{∆ k →l } K k=1 exp{∆ k →k } . ( 5
)
This implies that changing the label of node i to the label that maximizes the change in exact ICL (or alternatively not changing the node label if none of the changes in exact ICL are positive) corresponds to maximizing the full conditional distribution of the label in the Gibbs sampler. So instead of stochastically sampling a label from the full conditional, we deterministically choose the label that maximizes it in our greedy exact ICL search. This is similar to the iterated conditional modes algorithm of [START_REF] Besag | On the statistical analysis of dirty pictures (with discussion)[END_REF]. For this reason we argue that using the exact ICL as an objective function is sensible. The pragmatic nature of the algorithm comes from the reasoning that if there is a strong cohesion to subsets (strong clustering and good separation of node subsets), then the label that will be sampled from the full conditional will usually be the one which corresponds to the largest value of the change in exact ICL. If our objective is to get the optimal partitioning of the node sets, then this appears an appealing approach.

Computational savings for the greedy ICL algorithm

There are two main ways we can reduce the computational complexity of the greedy algorithm discussed in Section 3.4. The first of these concerns the correspondence between maximizing exact ICL and the full conditional distribution of labels as shown in Section 3.5. The second is by exploiting the fact that often observed bipartite networks are sparse.

Greedy search pruning

One might expect that after a few sweeps of the greedy algorithm, some degree of degeneracy will appear in the full conditional distributions of the labels. In such situations one will observe that π(c * i = l|c -i , w, Y, K, G) will be quite small compared with the full conditional for some other label, say π(c * i = k |c -i , w, Y, K, G). Suppose that currently c i = k. Compute ∆ k→l for all l = k, where ∆ k→k = 0. Let k = arg max r ∆ k→r . Consider labels l which are very unlikely for node i compared to k

π(c * i = l|c -i , w, Y, K, G) π(c * i = l|c -i , w, Y, K, G) + π(c * i = k |c -i , w, Y, K, G) < ε
for some small number ε. Using (5) this may be written

exp{∆ k→l } exp{∆ k→l } + exp{∆ k→k } < ε.
This can be rearranged to give

∆ k→k -∆ k→l > T (ε) (6)
where T (ε) is a threshold depending on the value of ε. To reduce our search space in future iterations of the greedy search, we ignore computing the exact ICL for label changes to l when ( 6) is satisfied. In this paper we take T (ε) = 150. This appeared to give satisfactory results without excluding important parts of the search space and introducing error. It can be seen that removal of these poor configurations from the search space will reduce the O(4Q Λ KG) cost of optimizing the labels at each step. The procedure has been described here for nodes from set A but applies analogously to those in set B.

Sparse representations

When dealing with observed bipartite networks one often encounters much sparsity. By sparsity we mean many observed non-ties. For this reason, it can be useful to use a sparse representation of the adjacency matrix Y , denoted Y s , and use this to perform the exact ICL calculations outlined earlier. Only the non-zero valued ties of Y are stored in Y s . For example, suppose we would like to compute the ICL difference given in (4) for Y s stored in sparse form. We denote a non-tie generically by "0". Suppose there are m i non-zero entries in the row corresponding to node i (from A) of Y s in columns j i 1 , . . . , j i m i . Let the indexes of the non-zero columns which correspond to nodes from B and are members of B g be stored in the set J i g = {l : w j l = g} and denote the nodes that have corresponding zero columns and which are members of B g by J i g = {l : w l = g}\J i g . For each g compute

υ i : nz g = m i l=1 I(w j i l = g)
which gives the number of non-zero entries in B g which in row i of the adjacency matrix, corresponding to node i. Then the number of zero entries in columns indexed by B g for row i is

υ i : z g = υ g -υ i : nz g . Then Λ kg = π(θ kg ) i:c i =k j:w j =g p(y ij |θ kg ) dθ kg = π(θ kg ) i:c i =k   j∈J i g p(y ij |θ kg ) j∈ J i g p(y ij |θ kg )   dθ kg = π(θ kg ) i:c i =k   p("0"|θ kg ) υ i : z g j∈J i g p(y s ij |θ kg )   dθ kg . (7)
Thus by storing the number of zeros, and carrying out the calculation for these terms collectively, we can save many evaluations. In the experiments below, we found that exploiting the sparsity of some problems gave a considerable speed up of the algorithm.

Algorithm types

The greedy search algorithm can be split into four types outlined in Table 1, algorithms A2 and A3 using the pruning ideas of Section 3.6.1, and A1 and A3 using the sparse representations of Section 3.6.2. Observing the algorithm types, A0 is the naive implementation and we would expect this to be the slowest. On the other end of the scale, we would expect A3 to be fastest. Between A1 and A2, we would expect A1 to be faster than A2 when the matrix is quite sparse, and A2 to be faster than A1 when there is strong separation of nodes into subsets (strong cohesion) i.e. very strong blocking. For small to moderate sizes of adjacency matrices, differences in run times may not be apparent, but for larger and sparser adjacency matrices we would expect to notice differences. When updating labels of A nodes or B nodes, these are processed in a random order. This introduces some amount of stochasticity into the ultimate destination of the greedy search. This will be investigated in due course in the examples below. For the examples we take α 0 and β 0 both equal to 1 throughout. For the pruning algorithms we allow five full sweeps of the data before pruning begins.

Simulation study

To investigate the performance of the greedy search approaches, we compared them with the MCMC approach of [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] on simulated data.

Study set-up

The data was generated from the model outlined in Section 2, with K = 5 and G = 5. The label distributions were chosen as ω k = 1/5, k = 1, . . . , 5 and ρ g = 1/5, g = 1, . . . , 5. The node sets were chosen to be of size N = 100 and M = 100, giving a 100 × 100 adjacency matrix. Using the absent/present linking attribute outlined in Section 3.2, the ties were generated so as to give different levels of overlap in the blocks. There are five blocks of varying clarity in the matrix following

θ kg = 1 -q if k = g q otherwise
where we let q vary from 0.0125 to 0.5 in steps of 0.0125, giving 40 values of q in total. As q gets closer to 0.5 the clustering task becomes more difficult, as the probability of a tie becomes constant across all blocks. Twenty datasets were generated at each value of q.

Comparing clusterings

As a means of comparing different clusterings of the data, we extend the normalized mutual information measure introduced in Vinh et al. ( 2010) and used by, for example [START_REF] Côme | Model selection and clustering in stochastic block models based on the exact integrated complete data likelihood[END_REF]. [START_REF] Vinh | Information theoretic measures for clusterings comparison: Variants, properties, normalization and correction for chance[END_REF] provide extensive and rigorous justification for using this such a measure for clustering comparison. We define a combined measure for the two node sets (rows and columns of the adjacency) to account for the fact that we are effectively doing two clustering tasks simultaneously. Consider first node set A. Denote the true clustering (from the simulation) by the labels c t and the estimated one (from the algorithm in question) by c e . The mutual information between the two clusterings is

I A (c e , c t ) = K k,l
P kl log P kl P e k P t l where

P kl = 1 N N i,j I(c e i = k, c t j = l), P e k = 1 N N i=1 I(c e i = k), P t l = 1 N N i=1 I(c t i = l).
The mutual information measures how much is learned about the true clustering if the estimated one is known, and vice versa. We normalize this quantity when the clusterings can have a different number of clusters (as in our case)

N I A (c e , c t ) = I(c e , c t ) max(H(c e ), H(c t )) with H A (c) = -K k=1 P k log P k and P k = 1 N N i=1 I(c i = k).
We compute the same quantity analogously for node set B, and add the two together to give an overall measure of mutual information between estimated and true clusterings as N I A (c e , c t ) + N I B (w e , w t ) which has a maximum value of 2 with high agreement and minimum 0 with little or no agreement.

Different approaches

The three different approaches taken were • a run of the collapsed MCMC algorithm of [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] of 25,000 iterations taking the first 5,000 as a burn in

• the basic greedy search algorithm A0

• the pruned greedy search algorithm A2.

The reason the sparse forms (A1 and A3) of the greedy search algorithms are not used here is the dimension of the adjacency matrix; a 100 × 100 matrix runs quickly enough using algorithms A0 and A2 that exploiting any possibly sparsity is not of huge benefit here. For the MCMC algorithm, the estimated clustering was computed by first finding the MAP estimates of K and G from the approximated posterior, and then finding the MAP of the labels conditioning on these values. It should be noted here that the MCMC takes considerably longer to run than either of the greedy algorithms, which are about 2,000 times faster in this instance. The pruned algorithm is usually about 10 times faster than the basic greedy algorithm. Figure 2 shows the average value of the normalized mutual information over the twenty datasets for each value of q for each of the algorithms. The MCMC algorithm is shown in blue, the greedy algorithm in red and the pruned greedy algorithm in green. It can be seen that as the value of q approaches 0.5, the greedy algorithms slightly outperform the MCMC algorithm, with the MCMC outperforming the greedy algorithm around q = 0.35. However, the observation most of note here is that the greedy algorithm gives essentially the same information as the longer MCMC run in a very small fraction of the time.

Examples

Congressional voting data

First of all the congressional voting data was analysed by treating abstains as nay's. The data matrix is 435 × 16. Here the node sets to be grouped are the congressmen (435) and the issues which they vote on ( 16). Ten instances of the algorithms in affects the results, one would expect this to me more pronounced in the pruned versions A2 and A3, since restarting here corresponds to resetting the search for each row/column to the set of all possible clusters as opposed to the pruned set.

To allow for direct comparison of the four algorithms, the random number generator was seeded the same-this also allows a direct time comparison for the same search progression as well as the maximum exact ICL reached. Table 2 shows the average runtime and maximum exact ICL reached for each algorithm over the ten runs. It can be seen that the sparse versions of the algorithm give a reduction in runtime with no noticeable reduction due to pruning.

The experiment was repeated, this time with different random seeds for each algorithm. Figure 3 shows the re-ordered data matrices for the four highest ICL over ten runs of the four algorithms. Here the effect of the stochasticity on the outcome of the greedy search is noted. The four algorithms converge to different numbers of clusters, with different resulting maximum exact ICL (Table 3). Visually, all four algorithms appear to give a sensible result despite the fact that there are some differences. This highlights the possibility of many local maxima in general applications and the need to run the algorithm a number of times to give it the best chance to find a close to global maximum.

To take a closer look at the effect of starting value and the randomness introduced into the greedy search by processing nodes in a random fashion when updating the labels, we took 100 runs of each algorithm. Histograms of the maximum ICL obtained are shown in Figure 4. From this it can be seen that all four algorithms most often get to a maximum exact ICL around -3560.00 and less frequently obtaining values larger than this. As mentioned in Section 3.3 and suggested by one reviewer, the optimal ICL could potentially be improved by providing a good starting value for the algorithm. In the absence of such an initialization, we'd recommend running the algorithm a number of times (10 or more) and taking the output of the run that gives the maximum ICL over these. This is similar to the general approach when the EM algorithm is fitted to finite mixture models; the run which results in the maximum value of the likelihood is taken as the output run.

As a rough comparison to the MCMC based algorithm of [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF], their analysis took about thirty minutes on this network, gleaning much the same information that we have obtained here in a fraction of the time. The MCMC analysis here was coded in parallel with the implementation for computing the exact ICL and is more efficient than the original implementation in [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF], hence the difference in reported run time. The mixing of their algorithm, even for this moderately sized data was slow for the congressman dimension with only about 1.5% of proposed cluster additions/deletions being accepted. In terms of scalability, MCMC is not an option when dealing with adjacency matrices of any large size. The fact that both the MCMC approach and the exact ICL approach are based on almost identical posteriors indicates that the greedy search is the only viable option of the two for larger applications.

Instead of using a random initialization for the clustering, it is possible to use the algorithm of [START_REF] Rohe | Co-clustering directed graphs to discover asymmetries and directional communities[END_REF]. Using the default (from the code mentioned in Section 3.3) algorithm configuration and an initial number of ten row and column clusters (chosen so that it is possible to use the algorithm of [START_REF] Rohe | Co-clustering directed graphs to discover asymmetries and directional communities[END_REF]) gave a maximum ICL of -3581.034 over ten runs with K = 4 row clusters and G = 9 column clusters. The lower number of clusters is most likely due to the lower number of initializing clusters. Also, using a different regularization parameter than the default value used (-1) could have led to a different number of clusters or value of the ICL (we suspect that the cluster initialization was already in a local maximum with respect to the ICL). The average run time was 0.37 seconds, which is faster than any of the other algorithms. This is, of course, impacted by the initial number of clusters assumed, but suggests this as a good candidate for initializing the algorithm.

Movie-Lens 100k

The Movie-Lens dataset gives 100,000 ratings of movies by users. In total, there are 943 users and 1,682 movies. Many users do not rate particular movies so the network is quite sparse (93%) and it is useful to use the approach mentioned in Section 3.6.2. The ratings (linking attributes) are assumed Poisson distributed, due to the ordinal nature of the ratings going from 1-5. This is a sensible model choice in this context due to both the ordinal and discrete nature of the ratings. The clustering task will cluster groups of individuals giving a similar rating to specific groups of movies. One could expect that this would glean more information than treating the 1-5 ratings as purely categorical since the Poisson model captures in some ways the ordinal features of the data. A missing rating is just given the value 0 when assuming the ratings are Poisson. The data is shown in Figure 5. One run of each of the algorithms was carried out each with two restarts with results in of clusters to those that do not involve pruning. This could be due to the effect of the value of the difference in exact ICL considered for the pruning. A value of this threshold which is too high could introduce some degeneracy into the search space, with labels settling in sub-optimal groups. The result for this dataset can vary depending on the pruning threshold chosen. However, pruning does give a considerable reduction in computing time which is favourable from a scalability viewpoint. This could be seen as an option when the data size is such that a run of the full algorithm is just too expensive.

Discussion

We have presented an approach to grouping or clustering subsets of nodes in bipartite networks which share similar linking properties. There are two main advantages to our work.

The first is that we model linking patterns between nodes in the two node sets in a statistical way. The second is that we provide a inference technique for the number of clusters in the node sets by appealing to a known information criterion due to [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. The fact that the exact ICL may be computed due to the tractability assumptions we make, means that it can be iteratively updated for each proposed new labelling of a node in either node set. The approach presented here can be considered favourable to the MCMC techniques of [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] for a few reasons. If only an optimal block structure is required our algorithm provides this in a more scalable framework than MCMC. Also, in summarizing MCMC output, often the maximum a posteriori (MAP) estimate is used as an overall summary, since it is not possible to average clusterings with a difference number of clusters. The concept of this MAP estimator is in somewhat equivalent to the configuration of labels and clusters that maximizes the ICL. Probably the most important reason is that MCMC is simply not a viable option for large bipartite networks due to the declining mixing performance observed and the long run times required in such cases. The simulation study shows that we gain as much information from the greedy algorithm as we would from a longer and more time consuming MCMC run, which is encouraging. This being said, there is always a risk that the ICL greedy algorithm we discuss could get trapped in a local maximum. Due to the fast run time, we can circumvent this somewhat in practice by running, say 10 instances of the greedy algorithm and taking the one that gives the highest ICL (indicating the best clustering). The greedy algorithm (especially the sparse versions) scale well enough to make this feasible. The risk that this will happen to the MCMC algorithm is less due to its stochastic nature. ICL does still possess many of the benefits of MCMC, the strongest of which is that it allows for automatic choice of the number of clusters. Considering the complexity of the task it performs (model selection and clustering), the ICL algorithm does appear to scale quite well, running on the Movie-Lens 100k data in less than 3 minutes, and the congressional voting data in less than a second. For many of the applications of bi-clustering in the literature, this is competitive. In addition, we provide code (written in C with an R wrapper) which is available from the first author's webpage https://sites.google.com/site/jsnwyse/ (a Makefile is provided to compile on Unix-alikes).

In Section 3.6 we showed how sparsity (a common feature of observed networks) can be exploited to provide a faster version of the greedy search exact ICL algorithm. Although not possible to provide exact quantification for the order of speed up (in terms of N and M ) we observed considerable computational gain in exploiting sparsity. The applicability of these ideas may stretch well beyond implementation of the LBM. For example, it may be possible to use the same kind of representation for SBMs.

We believe that the latent blockmodel is an appealing model for bipartite networks, although to our knowledge, no authors have appeared to explicitly discuss it in the context of network modelling. One drawback of the LBM in its original form was the lack of any principled information criterion (such as Bayesian Information Criterion) for choosing the number of clusters. However, the work of [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF] and the ideas presented in this paper overcome this issue. For this reason, the LBM can now be seen as a rich model for which there is capacity to model group structure in bipartite networks. The approach we take also means it is scalable for larger networks where sparsity can be exploited.

Figure 2 :

 2 Figure 2: Comparison of the normalized mutual information for the three algorithms. MCMC is shown in blue, greedy algorithm in red and the pruned greedy algorithm in green.

Figure 3 :

 3 Figure 3: Re-ordered data matrix of maximum ICL configuration output over the ten runs of each algorithm.

Figure 4 :Figure 5 :

 45 Figure 4: Histograms of maximum ICL for 100 runs of each of algorithms A0-A4.Movie-Lens data

  Flowchart showing the process for updating the label of a node from set A. Here the c

	Current labels and current ICL when node i in cluster k	ICL differences for change of label of node i from k to l, ∆ k →l	Identify the best label up-date for node i based on ICL	Move node i updating the label
		c * (1) , w giving ∆ k →1
			Take cluster l
			that gives max
		c * (2) , w giving ∆ k →2	∆ k →l and move
	c, w current labels		node i to this	c, w updated labels
		(K) , w giving ∆ k →K c * . . . . . .	cluster (leave put if max ∆ k →l is 0)
	Figure 1:

* are subscripted with the candidate cluster for clarity. The top portion (boldface)

Table 1 :

 1 Different types of greedy search algorithm.

	Algorithm Pruning Sparse form
	A0	No	No
	A1	No	Yes
	A2	Yes	No
	A3	Yes	Yes

Table 2 :

 2 Table 1 were run, each with two restarts and the maximum ICL recorded from the ten runs. If restarting the algorithm Results for ten runs of each of the four algorithms on Congressional voting data with same random seed.

	Algorithm maximum ICL average time (sec) (K, G)
	A0	-3543.062	0.62	(6,12)
	A1	-3543.062	0.56	(6,12)
	A2	-3543.062	0.62	(6,12)
	A3	-3543.062	0.56	(6,12)
	Algorithm maximum ICL average time (sec) (K, G)
	A0	-3546.968	0.63	(6,11)
	A1	-3543.812	0.58	(6,10)
	A2	-3537.503	0.65	(6,11)
	A3	-3537.550	0.64	(6,11)

Table 3 :

 3 Results for ten runs of each of the four algorithms on Congressional voting data with different random seed.

Table 4 .

 4 The reordered data is shown for algorithms A1 and A3 in Figure5.A glance at these results reveals that the pruning based algorithms give different number

	Algorithm maximum ICL run time (sec) (K, G)
	A0	-646597.7	1098.529	(55,64)
	A1	-646597.7	368.032	(55,64)
	A2	-646268.2	525.818	(56,62)
	A3	-646268.2	160.162	(56,62)

Table 4 :

 4 Results for one run of each of the four algorithms on Movie-Lens data with same random seed.

Acknowledgements

Part of this work was carried out while the first author worked at the Insight Centre for Data Analytics at University College Dublin. The Insight Centre for Data Analytics is supported by Science Foundation Ireland under Grant Number SFI/12/RC/2289. Nial Friel's research was also supported by a Science Foundation Ireland grant: 12/IP/1424.


Feng Zhang

André Nicolle

Lili Xing

Stephen J Klippenstein

Recombination of Aromatic Radicals with Molecular Oxygen

Keywords: hydrocarbon peroxy radicals, aromatic radicals, recombination reactions, quantum chemical calculations, variable reaction coordinate transition state theory

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Aromatic compounds contribute greatly to the emission of polycyclic aromatic hydrocarbons (PAH) and soot and are key precursors to the formation of tropospheric ozone and secondary organic aerosols [1,2]. As such, an accurate understanding of their oxidation kinetics is crucial to the modeling of both atmospheric and combustion chemistry [3][4][5][START_REF] Calvert | The Mechanisms of Atmospheric Oxidation of Aromatic Hydrocarbons[END_REF]. In the low temperature oxidation environment, ROO chemistry is of direct relevance to the chain reactions that ultimately lead to autoignition [4,[START_REF] Calvert | The Mechanisms of Atmospheric Oxidation of Aromatic Hydrocarbons[END_REF]. The recombination reaction between hydrocarbon radicals and oxygen molecule is a significant source of ROO, triggering the subsequent ROO chemistry.

Combustion modeling studies on the oxidation of aromatic hydrocarbon fuels such as benzene, toluene, and naphthalene (see, e.g, [START_REF] Yuan | [END_REF][8][9]) have generally estimated the kinetic data for the fuel radical R + O2 reaction from experience or via extrapolation of relatively low-temperature experimental measurements [10][11][12]. For phenyl, measurements of Yu et al. [10] over the 297-473 K temperature range suggest that the O2 recombination rate constant decreases with temperature, while Schaugg et al.'s measurements showed a positive temperature dependence at slightly higher temperatures (418-815 K) [START_REF] Schaugg | Rate Coefficients for the Fast Reactions of Phenyl Radicals with NO and O2 from Mass Spectrometric Measurements of Phenyl Decays, Transport phenomena in combustion[END_REF]. For benzyl, both Hoyermann et al. [START_REF] Hoyermann | [END_REF] and Nelson et al. [15] observed no temperature dependence over the 243-373 K range, while Fenter et al. indicate a significant decline with increasing temperature at analogous temperatures (298-398 K) [11]. Only Park et al. have measured the rate constant for 2-naphthyl, finding a modest decrease over the 299-444 K temperature range. These discrepancies in the observed temperature behavior add to the uncertainty in any attempted extrapolations to combustion temperatures.

There is also discordancy in recent theoretical predictions of the temperature dependence for aromatic radical + O2 recombination kinetics. The canonical variational transition state theory (VTST) calculations of da Silva et al. [16] predict that the rate coefficient for phenyl + O2 should decrease by about a factor of 4 over the 300-2000 K temperature range. In contrast, the very recent variable reaction coordinate transition state theory (VRC-TST) [17] calculations of Kislov et al. [18] predict that for phenyl, 1-and 2-naphthyl radicals the O2 recombination rate constant should rise over the 1500-2500 K range, in qualitative accord with the results of Schaugg et al. [START_REF] Schaugg | Rate Coefficients for the Fast Reactions of Phenyl Radicals with NO and O2 from Mass Spectrometric Measurements of Phenyl Decays, Transport phenomena in combustion[END_REF]18].

In this work, the O2 recombination reaction was studied for four representative aromatic hydrocarbon radicals: phenyl, benzyl, 1-and 2-naphthyl radicals. The radical characteristics of triplet O2 result in a barrierless minimum energy path for its addition to radicals [19,20]. Prior applications of the VRC-TST approach, such as the recent application to the closely related C2H3 + O2 addition reaction [19], have amply demonstrated its ability to accurately predict the rate coefficients for such barrierless radical-radical reactions.

Hence, the present analysis employs VRC-TST as in the work of Kislov et al. [18], but proceeds beyond it in a variety of ways. Their approximate representation of the transitional mode interaction potential is replaced here with direct CASPT2 determinations. Furthermore, accurate estimates of the minimum energy path (MEP) energies are obtained through multireference based evaluations of the spin-splitting combined with coupled-cluster based evaluations for the high spin quartet state, as in the C2H3 + O2 study.

Methodology

Kinetic theory

The key assumption in VRC-TST is the separation of the "conserved" and "translational" modes in the transition state region of the potential. The conserved modes, which correspond to the vibrational modes of the reacting fragments, are treated as quantized harmonic oscillators. The transitional modes correspond to the remaining modes, which describe the rotational and translational motions of the two reacting fragments. The essential advantage of the VRC-TST approach is its quantitative treatment of the full coupling and anharmonicity of the transitional modes through Monte Carlo integration over the coordinate space degrees of freedom in the classical phase space representation for the reactive flux for a variable set of reaction coordinates.

The variability of the reaction coordinate in VRC-TST is achieved by choosing a set of pivot points whose locations and separations determine the transition state dividing surface. The rate constant is then minimized with respect to both the pivot point locations and separation distances.

For the four R + O2 reactions studied in this work, two sets of dividing surfaces were used. The first set employs center-of-mass (CoM) pivot points, which are appropriate for large separations, while the second set employs orbital centered pivot points, which are appropriate for shorter separations.

A separate 6-point-mesh with a variable grid spacing of 1.0-4.0 au was used for the larger separations (9.0-21.0 au) where the CoM pivot points are used. A 12-point-mesh with finer grids (0.3-1.0 au) was used for the shorter separations (3.8-9.0 au) with orbital centered pivot points. For the orbital centered pivot points, their displacement from the bonding atoms was varied from 0.5 to 1.5 au, with a step size of 0.5 au. With the preferred E,J-resolved VRC-TST, which was employed here, the reaction coordinate minimizations are performed for each energy, E, and total angular momentum J. A dynamical correction (0.85) was used for the final rate coefficients.

Quantum chemical calculations

The implementation of the VRC-TST approach requires accurate interfragment interaction energies for ~10 4 configurations for each reaction. The program VaReCoF [21], which was used in the present work, allows for the evaluation of these interactions on-the-fly, with direct calls to ab initio electronic structure methods. The complete active space self-consistent field theory with second-order perturbation (CASPT2) [22] method combined with the cc-pVDZ basis set efficiently yields interaction energies of suitable accuracy [19,20]. Here, a 7-electron 5-orbital (7e,5o) active space was chosen to optimize the MEP for all four R + O2 reactions. This space correlates with the two pairs of O-O π and π* orbitals, and the radical orbital of R. An ionization potential-electronic affinity (IPEA) shift of 0.25, was introduced into the CASPT2 calculations [23].

The CASPT2/cc-pVDZ energy computed for each sampled configuration is corrected by an orientationally independent factor [ΔVtot (RCO)], where RCO is the bond length for the incipient CO bond. Since the internal coordinates of the fragments are fixed during the VRC-TST evaluations, this correction factor includes both a relaxation correction and a higher level energy correction. The correction factor was obtained by multi-step composite computations: first, the relaxation correction (ΔVrelax) was defined as the energy difference between the MEP built by rigid and relaxed scanning at the CASPT2(7e,5o)/cc-pVDZ level; second, along the MEP, CASPT2 calculations were performed with the cc-pVTZ and cc-pVQZ basis sets, and the energies were extrapolated to the complete basis set limit (CBS) [24]. A basis set correction (ΔVbasis) is defined as the CASPT2/CBS to CASPT2/cc-pVDZ energy difference; finally, a method correction (ΔVmethod) was obtained from CCSD(T)-F12/VDZ-F12 [25,26], MP2-F12/VDZ-F12 [27] and MP2-F12/VTZ-F12 evaluations for the high-spin quartet state:

∆ / / / / (E1)
Here the superscript "4" denotes an energy calculated for the quartet state. For phenyl + O2, the smaller size allowed for additional calculations with larger basis sets, i.e., cc-pVTZ and cc-pVQZ. In summary, the total correction term (ΔVtot) includes three parts, i.e., a geometry relaxation correction, a basis set correction, and a method correction.

For the present method correction, the interaction energy is approximated as the sum of the CASPT2 calculated quartet-doublet splitting subtracted from the CCSD(T) calculated quartet energy.

Note that direct application of CCSD(T) to the doublet state MEP does not yield accurate energies due to the multi-reference nature of the doublet wavefunction in the transition state region. The method correction is effective because the CCSD(T) method accurately treats the high-spin quartet state, while multi-reference treatments of the spin-splitting are more consistent and reliable than corresponding treatments of the absolute interaction energy.

For phenyl and benzyl, we also explored the difference between CASPT2 and multireference configuration interaction MEPs, with the latter performed with (MRCI+Q) and without (MRCI) the Davidson correction. For benzyl, we also explored a CASPT2(13e,11o) active space (which includes the additional (6e,6o) π-space) in order to improve our treatment of the resonance stabilization of the radical.

Reaction enthalpies for benzyl-and naphthyl-+ O2 reactions were computed at the CCSD(T)-F12/VDZ-F12, MP2-F12/VDZ-F12, and MP2-F12/VTZ-F12 levels. The final single point energies were obtained from the expression

/ / / (E2)
The geometries and zero-point energies for these reaction enthalpy evaluations were obtained with the B2PLYP-D3 density functional method [28], employing the aug-cc-pVDZ basis for C6H5OO and C7H7OO, and the cc-pVDZ basis for C10H7OO. These density functional theory (DFT) calculations were performed with the Gaussian program [START_REF] Frisch | Gaussian 09, Revision B.01[END_REF]. All of the present CASPT2, CCSD(T)-F12, and MP2-F12 calculations were performed with the Molpro program package [START_REF] Werner | [END_REF].

Results and Discussions

Structures and reaction enthalpies

The B2PLYP-D3 optimized structures and reaction enthalpies (in brackets) for the ROO complexes arising from the addition of O2 to phenyl, benzyl, 1-and 2-naphthyl radicals are illustrated in Fig. 1. The Cartesian coordinates for all optimized structures are provided in the Supplemental Material. There are three distinct ROO adducts for the benzyl radical arising from its various resonance configurations. The C7H7O2 radical with O2 added on the methyl branch [i.e., the C7H7O2-1 radical shown in Fig. 1 kcal/mol more stable than 1-C10H7OO-L. In contrast, the -OO orientation has very little effect on the structures or enthalpies for the O2 addition to form 2-C10H7OO, which is similar to the symmetric C6H5OO structure.

The one-dimensional energy profiles required to evaluate ΔVtot were calculated for both 1-C10H7OO-R and 1-C10H7OO-L. The correction energy used in the VRC-TST calculations for 1-C10H7 + O2 is the mean value of that from the 1-C10H7OO-R and 1-C10H7OO-L MEPs. For 2-C10H7OO formation only the -R configuration was considered as the energy correction for 2-C10H7OO-L is expected to be nearly identical.

Phenyl + O 2

Figure 2(a) illustrates the basis set and method dependence of the calculated MEP for the phenyl + O2 association reaction over RCO distances of 2 to 4 Å, which roughly spans the range of separations for the variational transition state. The potential energy curves were evaluated for geometries along the CASPT2(7,5)/cc-pVDZ relaxed MEP. The CBS correction is seen to decrease the attractiveness at long-range while increasing it at shorter RCO. The M1 method correction (from Eq. (E1)) has only a minor effect on the CASPT2/CBS MEP. Although not shown here, employing an expression similar to Eq. (E1) but using cc-pVTZ and cc-pVQZ bases, yields a method correction, M2, that is essentially identical to M1.

Figure 2(b) reveals the effect of the method correction on the calculated CBS limit MEP for CASPT2, MRCI, and MRCI+Q calculations. The closer agreement between the three solid lines indicates that the method correction provides a more consistent description of the MEP. After the method correction, the maximum difference between the MRCI+Q and CASPT2 MEP energies is only 25%.

Fig. 2 MEP energies for the phenyl + O2 reaction from different corrections. PT2 refers to CASPT2, CI refers to MRCI, while M1 represents the method correction evaluated from Eq. (E1).

Figure 3 illustrates the present VRC-TST predictions for the phenyl + O2 high pressure recombination rate coefficient. The black solid and dash lines illustrate our two best results. They include the CBS, relaxation, and M2 method correction and differ only in whether the CBS and method corrections are evaluated with CASPT2 or MRCI+Q calculations. The difference between the two provides some indication of the uncertainty in our predictions. At room temperature, the nearly factor of two difference is due to the strong sensitivity of the predicted rate to the energetics at RCO ~2.5 Å. Above 1000 K the two results are in much improved agreement, differing by less than a factor of 1.35. In comparing with experiment we take the average of these two predictions as our best prediction.

The blue line with circles illustrates the results obtained when all correction terms are ignored.

At room temperature this reference result is nearly 3 times larger than the best results. But by 1000 K the reference and fully corrected results are within about 20%. The uncertainty in our predictions is smaller when these results agree, because the error arising from the 1-dimensional assumptions implicit in the corrections should then be less important. When only the basis set and relaxation correction are included (pink solid line) the rate constants are overestimated by up to 30% compared with those including the method correction.

Fig. 3 VRC-TST predicted recombination rate constants for phenyl + O2.

The computed rate constants show a relatively modest temperature dependence, with first a negative T dependence for T < 900 K, and then a positive T dependence for higher T. This minimum in the rate constant is indicative of the delicate balance between entropic and enthalpic effects. For both T ranges, the transition state is moving to shorter separations with increasing temperature and entropic changes correlate with a negative T dependence, while enthalpic changes correlate with a positive T dependence. At low temperature, the transition state lies at large separation, where the MEP has a small slope and the decreasing entropy dominates the T dependence. At higher T, the transition state lies at shorter separation where the MEP has a larger slope and the enthalpic changes dominate the T dependence.

Our best prediction for phenyl + O2 is compared with previously reported experimental and theoretical rate constants in Fig. 4. The VRC-TST predictions are seen to be in good agreement with Yu et al.'s measurements with the CRD technique, and even accurately reproduce the observed negative temperature dependence [10]. In contrast, there is a clear discrepancy between the present VRC-TST predictions and the observation by Schaugg et al. in their fast flow reactor, which shows a positive T dependence over the T range of 418-815 K at 0.7 mbar [START_REF] Schaugg | Rate Coefficients for the Fast Reactions of Phenyl Radicals with NO and O2 from Mass Spectrometric Measurements of Phenyl Decays, Transport phenomena in combustion[END_REF]. Recently, Tanaka et al. also studied this association reaction with the CRD method, measuring the rate coefficient at 296 K and 30 Torr [START_REF] Tanaka | [END_REF], which is in good agreement with our predictions and with the measurements of Yu et al. [10]. The VTST results of da Silva et al. [16], which are based on a G3B3 scaled O3LYP/6-31G(d) PES, are in reasonable agreement with the present predictions, differing by a factor of 2 or less.

However, they do show a stronger negative temperature dependence and fail to predict the slightly positive temperature dependence at high temperature. Very recently, Kislov et al. employed VRC-TST to compute the rate constant at high temperatures (1500-2500 K) [18,[START_REF] Klippenstein | Variflex, version 1.0[END_REF]. Their results, which are based on an approximate representation of the transitional mode potential from fits to CASPT2(19,14)//CASSCF(9,9)/aug-cc-pVDZ evaluations, are ~3.5 times greater than ours and show a stronger temperature dependence.

Fig. 4 Comparison between our best prediction and previously reported rate constants for phenyl + O2 recombination.

1-and 2-naphthyl + O 2

Figure 5 provides a comparison of the MEP for O2 adding to 1-and 2-naphthyl. The potential energy has been corrected by all three terms, with Eq. (E1) employed for the method correction.

1-C10H7OO-L and 1-C10H7OO-R are distinguished by the orientation of the -OO group, as indicated by Fig. 1. As expected, the potential energy for 1-C10H7OO-R is very similar to that for 2-C10H7OO.

In contrast, there are considerably greater steric interactions with nearby H atoms for 1-C10H7OO-L, which is up to 2.5 kcal/mol more repulsive than 1-C10H7OO-R in the RCO range of 2-4 Å. The present predictions for the high pressure recombination rates for naphthyl + O2 are illustrated in Fig. 6 together with previously reported data, including for phenyl + O2. The present calculations suggest that the rate constant for adding to 1-naphthyl is about 1.5 times lower than that for adding to phenyl. This is expected on the basis of the increased steric repulsion for 1-naphthyl. In contrast, the predicted rate constant for adding to 2-naphthyl is somewhat greater than that for adding to phenyl. The predicted modest decrease with increasing temperature is again in reasonable accord with the experiments of Park et al. [10,12]. However, the experimental data suggest that O2 addition to 2-naphthyl should actually be slower than to phenyl (by about a factor of 2) for temperatures in the 299-444 K range. The reason for this discordance is not clear. It is perhaps worth noting that the theoretical analysis of Kislov also predicts a modest increase in the rate constant from phenyl to 2-naphthyl. However, their rate predictions are again about 3 times higher than the present ones.

Fig. 6 Comparison between VRC-TST rate constants and previous studies for C6H5 + O2, 1-C10H7 + O2, and 2-C10H7 + O2.

Benzyl + O 2

The present kinetic analysis for benzyl radical focuses on the formation of the C7H7O2-1 radical since additions to the radical site for the two other resonance structures of C7H7 are endothermic (cf., Fig. 1). Figure 7 illustrates the MEP for formation of C7H7O2-1 in a trans conformation in Cs symmetry at various levels. Notably, due to the resonance stabilization in the C7H7 radical, which must be broken before the radical-O2 bond is formed, the MEP is significantly less attractive than it is for phenyl. Indeed, for most methods there is even a saddle point at about RCO = 2.3 Å.

The MEP curves from the CASPT2(7e,5o), CASPT2(13e,11o), and MRCI+Q calculations are again more consistent after the method correction is included. However, there is still considerable variation in the predicted MEP potential curves even after including the method correction. For example, the CASPT2(7e, 5o) and CASPT2(13e, 11o) predicted saddle points differ by 0.9 kcal/mol.

The MRCI+Q MEP is intermediate in value, and might be expected to best handle the effects of resonance stabilization of the radical. The present MRCI+Q based predictions, which are expected to provide the best results, are in reasonably satisfactory agreement with the experimental data. They show a very modest minimum in the rate constant near 400 K and the decline from room temperature to 400 K is in reasonable accord with experiment. Near room temperature the CASPT2(7e,5o) and CASPT2(13e,11o) predictions are about a factor of two higher and lower, respectively. With increasing temperature the discrepancy between the three predictions decreases somewhat. The discrepancy between these predictions roughly correlates with the overall uncertainty in our best prediction for this reaction.

Conclusions

The high pressure recombination kinetics for O2 adding to a series of prototypical aromatic hydrocarbon radicals was studied with direct CASPT2 based VRC-TST theory. The predictions for phenyl, benzyl, 1-and 2-naphthyl radicals over the 300 to 2500 K temperature range are accurately reproduced by the modified Arrhenius expressions reported in the Supplementary material. The predicted kinetics for C6H5, 1-C10H7, and 2-C10H7 are fairly similar, with the rates being lowest for 1-C10H7 due to a modest steric effect. In contrast, the resonantly-stabilized character of the benzyl radical leads to much lower reaction enthalpies and rate coefficients.

For the O2 addition to C6H5, 1-C10H7, and 2-C10H7, the calculations predict a modest negative temperature dependence near room temperature, then show a minimum near 600 K, and finally rise by up to a factor of 2 by 2500 K. The predicted temperature dependence for the rate constants is in excellent agreement with many of the experimental observations, which are generally limited to temperatures of 500 K and lower. The predicted magnitudes are in satisfactory agreement with experiment, with maximum discrepancies within 50%. The resonantly-stabilized benzyl radical is more difficult to treat accurately, especially near room temperature where the predictions are sensitive to the predicted saddle point energy. Nevertheless, our best predictions for benzyl are still in reasonable agreement with experiment. The present systematic kinetic studies for the oxidation of these four typical aromatic radicals enhances our understanding of the structure-activity relationship for such reactions, and should help to improve the rate rules for R + O2 reactions.

  (b1)] is ~19.1 kcal/mol (in the trans conformation) more stable than C7H7 + O2, while the C7H7O2-2 and C7H7O2-3 radicals [cf. Fig.1(b2) and (b3)] are 5-10 kcal/mol endothermic relative to C7H7 + O2. Only the trans conformation was considered in the calculations of the one-dimensional corrections for the C7H7O2-1 radical.The orientation of the -OO group in the C10H7OO complexes has two forms, denoted as C10H7OO-R (right), and C10H7OO-L (left) in Fig.1 (c) and (d). The steric repulsion arising from the interaction of the -OO group with the nearest H atom is more significant in 1-C10H7OO than in 2-C10H7OO. Indeed, 1-C10H7OO-R maintains Cs symmetry, while the OOCC dihedral in 1-C10H7OO-L is ~50 degrees. Even with this torsional rotation in 1-C10H7OO-L, the distance from the O to the nearest H is still 0.08 Å shorter than in 1-C10H7OO-R, implying greater steric hindrance in the former configuration. According to the computed reaction enthalpies, 1-C10H7OO-R is 2.2

Fig. 1

 1 Fig. 1 Optimized ROO complexes for (a) phenyl, (b) benzyl, (c) 1-naphthyl, and (d) 2-naphthyl radicals.

Fig. 5

 5 Fig. 5 MEP for 1-C10H7 + O2 and 2-C10H7 + O2.

Fig. 7

 7 Fig. 7 Potential energy curves for benzyl + O2 at various levels. PT2(75) refers to the CASPT2(7e,5o) method, PT2(1311) refers to the CASPT2(13e,11o) method, and (CI+Q) refers to MRCI+Q(7e,5o) energies.

Figure 8 Fig. 8

 88 Figure 8 presents a comparison of our VRC-TST predictions for the benzyl + O2 high pressure

Acknowledgements

This material is based in part on work at Argonne supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences under Contract No. DE-AC02-06CH11357. This work is also financially supported by National Natural Science Foundation of China under grants 21303173 and 51376170. We gratefully acknowledge numerous constructive discussions with Yuri Georgievskii and Lawrence B.

Harding.