Bi-SCORE for Weighted Bipartite Networks with Application in Knowledge Source Discovery

TL;DR

This paper introduces Bi-SCORE, a spectral clustering method tailored for weighted bipartite networks, effectively uncovering community structures in citation networks to reveal key research areas in statistics.

Contribution

The paper presents Bi-SCORE, a novel spectral clustering algorithm with theoretical guarantees, specifically designed for weighted bipartite networks with degree heterogeneity, applied to citation data.

Findings

Bi-SCORE successfully detects six distinct research communities.

The method outperforms existing community detection techniques.

Application reveals meaningful insights into statistical knowledge flow.

Abstract

Community detection in citation networks offers a powerful approach to understanding knowledge flow and identifying core research areas within academic disciplines. This study focuses on knowledge source discovery in statistics by analyzing a weighted bipartite journal citation network constructed from 16,119 articles published in eight core journals from 2001 to 2023. To capture the inherent asymmetry of citation behavior, we explicitly preserve the bipartite structure of the network, distinguishing between citing and cited journals. For this task, we propose Bi-SCORE (Bipartite Spectral Clustering on Ratios-of-Eigenvectors), a computationally efficient and initialization-free spectral method designed for community detection in weighted bipartite networks with degree heterogeneity. We establish rigorous theoretical guarantees for the performance of Bi-SCORE under the weighted bipartite degree-corrected stochastic block model. Furthermore, simulation studies demonstrate its robustness across varying levels of sparsity and degree heterogeneity, where it outperforms existing methods. When applied to the real-world citation network, Bi-SCORE uncovers a six-community structure corresponding to key research areas in statistics, including applied statistics, methodology, theory, computation, and econometrics. These findings provide valuable insights into the intricate citation patterns and knowledge flow among statistical journals.

Tables

Table 1. Table 1: An example of the basic information for publications, including the title, authors, publisher, publication year, citation count, and reference list.

Variable	Detail
Title	Using Heteroscedasticity to Identify and Estimate Mismeasured and Endogenous Regressor Models
Authors	Lewbel, A
Publisher	Journal of Business & Economic Statistics
Publication Year	2012
Keywords	Endogeneity, Heteroscedastic errors, Identification, Measurement error, Partly linear model, Simultaneous system
Abstract	This article proposes a new method of obtaining identification in mismeasured regressor models, triangular systems, and simultaneous equation systems. The method may be used in applications where other sources of identification, such as instrumental variables or repeated measurements …
Citation Count	1,315
Reference List	- Efficient estimation of models with conditional moment restrictions containing unknown functions (Econometrica) - The evidence on credit constraints in post-secondary schooling (Economic Journal) - Two-step GMM estimation of the errors-in-variables model using high-order moments (Econometric Theory) … (and more)

Table 2. Table 2: Number of publications, average citation, average reference, and self-citation rate for core journals from 2001 to 2023.

Journal Name	Number of Publications	Average Citation	Average Reference	Self-Citation Rate
AOAS	1,793	37.05	42.38	1.68%
AOS	2,627	79.84	36.19	16.24%
Biometrika	1,861	45.46	25.68	10.82%
JRSS-B	1,290	113.62	37.84	5.99%
JASA	3,758	65.51	37.96	10.17%
JCGS	1,679	34.27	33.76	3.36%
JBES	1,451	49.97	38.34	3.89%
SC	1,660	34.47	36.45	2.61%

Table 3. Table 3: Eight representative journals. Journal name, number of publications and JCR category are presented.

Journal Name	Number of Publications	JCR Category
Annals of Applied Statistics (AOAS)	1677	STATISTICS & PROBABILITY
Annals of Statistics (AOS)	2464	STATISTICS & PROBABILITY
Biometrika	1804	BIOLOGY MATHEMATICAL & COMPUTATIONAL BIOLOGY STATISTICS & PROBABILITY
Journal of the Royal Statistical Society Series B–statistical Methodology (JRSS-B)	1129	STATISTICS & PROBABILITY
Journal of the American Statistical Association (JASA)	3185	STATISTICS & PROBABILITY
Journal of Computational and Graphical Statistics (JCGS)	1614	STATISTICS & PROBABILITY
Journal of Business & Economic Statistics (JBES)	1271	ECONOMICS SOCIAL SCIENCES, MATHEMATICAL METHODS STATISTICS & PROBABILITY
Statistics and Computing (SC)	1613	COMPUTER SCIENCE, THEORY & METHODS STATISTICS & PROBABILITY

Table 4. Table 4: Top 10 references for core journals . Journal name abbreviation is presented with its weight in parentheses.

Journal Name	Top 10 References	Journal Name	Top 10 References
AOAS	JASA (4189)	JASA	JASA (12294)
	JRSS-B (2125)		AOS (8653)
	AOS (1751)		Biometrika (5396)
	Biometrika (1686)		JRSS-B (5122)
	Biometrics (1638)		Biometrics (3086)
	AOAS (1191)		Econometrica (1600)
	STAT MED (857)		STAT MED (1514)
	Bioinformatics (797)		JOE (1392)
	PNAS (794)		STAT SINICA (1382)
	Nature (712)		JMLR (1183)
AOS	AOS (14471)	JCGS	JASA (4366)
	JASA (4922)		AOS (2795)
	JRSS-B (2993)		JRSS-B (2549)
	Biometrika (2859)		Biometrika (1994)
	JMA (1380)		JCGS (1829)
	Bernoulli (1279)		Biometrics (839)
	TIT (1253)		CSDA (755)
	AMS (1180)		JMLR (740)
	JMLR (1152)		SC (704)
	STAT SINICA (1060)		STAT SINICA (650)
Biometrika	Biometrika (4998)	JBES	JOE (4411)
	JASA (4368)		Econometrica (3929)
	AOS (4176)		JASA (2340)
	JRSS-B (2518)		JBES (1894)
	Biometrics (1442)		AOS (1650)
	STAT MED (688)		ECONOMET THEOR (1008)
	STAT SINICA (643)		JOF (987)
	Econometrica (571)		AER (971)
	STAT SCI (506)		JAE (916)
	JMA (480)		RES (876)
JRSS-B	AOS (3618)	SC	JASA (3273)
	JASA (3593)		JRSS-B (2501)
	JRSS-B (2542)		AOS (2317)
	Biometrika (2320)		Biometrika (1869)
	Biometrics (831)		SC (1535)
	Econometrica (522)		CSDA (1163)
	STAT SCI (493)		JCGS (1036)
	JMLR (491)		Biometrics (736)
	STAT SINICA (472)		JMLR (727)
	ARXIV (420)		ARXIV (658)
\botrule

Table 5. Table 5: Journals in the “Natural Science” subcommunity. Journal name and weight are presented, where the weight refers to the total citation counts from the core journals.

Journal Name	Weight
Monthly Notices of the Royal Astronomical Society (MNRAS)	145
Methods in Ecology and Evolution (MEE)	83
Astronomy and Astrophysics (A&A)	80
Atmospheric Environment	79
PLoS Medicine (PLoS MED)	62
Journal of Hydrology	61
Analytical Chemistry	60
Demography	59
Astronomical Journal (AJ)	58
Journal of Neurophysiology	57
Philosophical Transactions of the Royal Society B	49
BMC Genomics	45
Remote Sensing of Environment (RSE)	45
Biophysical Journal	44
Climate Dynamics	44
Journal of Evolutionary Economics	44
PLoS Biology (PLoS Biol)	44
Journal of Theoretical Biology	43
Statistical Methodology	43
Blood	42
Clinical Cancer Research	42
Journal of Neuroscience Methods	42
Environmental Science and Technology	41
Nature Reviews Neuroscience	41
Cancer Cell	40

Equations

$$\mathbb{E}(A_{ij})={\theta }_i{\gamma }_j{B}_{c_i^r{c}_j^c},$$

$$\Omega =\mathbb{E}(A)=\Theta ZB{W}^{\top }\Gamma ,$$

$${\hat{R}}_{ik}^r=\{\begin{array}{ll}{\tau }_n,& \text{if }\frac{{\hat{U}}_{i(k+1)}}{{\hat{U}}_{i1}}>{\tau }_n\\ \frac{{\hat{U}}_{k+1}(i)}{{\hat{U}}_1(i)},& \text{if }\left|\frac{{\hat{U}}_{i(k+1)}}{{\hat{U}}_{i1}}\right|\le {\tau }_n\\ -{\tau }_n,& \text{if }\frac{{\hat{U}}_{i(k+1)}}{{\hat{U}}_{i1}}<-{\tau }_n\end{array},{\hat{R}}_{jk}^c=\{\begin{array}{ll}{\tau }_m,& \text{if }\frac{{\hat{V}}_{j(k+1)}}{{\hat{V}}_{j1}}>{\tau }_m\\ \frac{{\hat{V}}_{k+1}(i)}{{\hat{V}}_1(i)},& \text{if }\left|\frac{{\hat{V}}_{j(k+1)}}{{\hat{V}}_{j1}}\right|\le {\tau }_m.\\ -{\tau }_m,& \text{if }\frac{{\hat{V}}_{j(k+1)}}{{\hat{V}}_{j1}}<-{\tau }_m\end{array}$$

$${\hat{X}}^r=\underset{X^r\in {\mathcal{X}}_{n,\kappa -1,K}}{\text{mbox}argmin}{\Vert X^r-{\hat{R}}^r\Vert }_F^2,{\hat{X}}^c=\underset{X^c\in {\mathcal{X}}_{m,\kappa -1,L}}{\text{mbox}argmin}{\Vert X^c-{\hat{R}}^c\Vert }_F^2,$$

$$B\text{ is non-singular},0\le B_{kl}\le 1\text{for}k=1,\dots ,K,\text{and}l=1,\dots ,L,$$

$$B{B}^{\top }\text{and}{B}^{\top }B\text{are non-negative and irreducible.}$$

$${\theta }_i^{(k)}=\{\begin{array}{ll}{\theta }_i& \text{if }{c}_i^r=k\\ 0& \text{if }{c}_i^r\ne k\end{array}\text{and}{\gamma }_j^{(l)}=\{\begin{array}{ll}{\gamma }_j& \text{if }{c}_j^c=l\\ 0& \text{if }{c}_j^c\ne l.\end{array}$$

$$\mathcal{Z}=\max ({\theta }_{\max },{\gamma }_{\max })\max (\Vert \theta {\Vert }_1,\Vert \gamma {\Vert }_1),$$

$${\theta }_{\min }>0,{\gamma }_{\min }>0,$$

$$\Vert {\theta }^{(k_1)}\Vert \asymp \Vert {\theta }^{(k_2)}\Vert ,\Vert {\gamma }^{(l_1)}\Vert \asymp \Vert {\gamma }^{(l_2)}\Vert \text{for}1\le k_1,k_2\le K,\text{and}1\le l_1,l_2\le L,$$

$$\underset{n,m\to \infty }{\lim }\frac{\log (nm)\mathcal{Z}+{\log }^2(nm)}{{\theta }_{\min }{\gamma }_{\min }\Vert \theta {\Vert }_1\Vert \gamma {\Vert }_1}=0.$$

$$\Vert {\theta }^{(k)}\Vert \asymp \Vert \theta \Vert \text{and}\Vert {\gamma }^{(l)}\Vert \asymp \Vert \gamma \Vert \text{for}1\le k\le K\text{and}1\le l\le L.$$

$$\underset{n,m\to \infty }{\lim }\frac{\log (nm)\mathcal{Z}+{\log }^2(nm)}{\Vert \theta {\Vert }^2\Vert \gamma {\Vert }^2}\le \underset{n,m\to \infty }{\lim }\frac{\log (nm)\mathcal{Z}+{\log }^2(nm)}{{\theta }_{\min }{\gamma }_{\min }\Vert \theta {\Vert }_1\Vert \gamma {\Vert }_1}=0$$

$${\sigma }_k(\Omega )=\{\begin{array}{ll}\Vert \theta \Vert \Vert \gamma \Vert {\sigma }_k(S)& \text{if }1\le k\le \kappa ,\\ 0& \text{if }k>\kappa ,\end{array}$$

$$U_{i\cdot }=\frac{{\theta }_i}{\Vert {\theta }^{(c_i^r)}\Vert }{Y}_{c_i^r\cdot }\text{for}1\le i\le n\text{and}{V}_{j\cdot }=\frac{{\gamma }_j}{\Vert {\gamma }^{(c_j^c)}\Vert }{H}_{c_j^c\cdot }\text{for}1\le j\le m,$$

$$U_{\cdot k}=\sum _{k^{\prime }=1}^K\frac{{\theta }^{(k^{\prime })}}{\Vert {\theta }^{(k^{\prime })}\Vert }{Y}_{k^{\prime }k}\text{for}1\le k\le \kappa ,\text{and}{V}_{\cdot l}=\sum _{l^{\prime }=1}^L\frac{{\gamma }^{(l^{\prime })}}{\Vert {\gamma }^{(l^{\prime })}\Vert }{H}_{l^{\prime }l}\text{for}1\le l\le \kappa .$$

$$\Vert U_{i\cdot }\Vert \asymp \frac{{\theta }_i}{\Vert \theta \Vert }\text{for}1\le i\le n,\text{and}\Vert V_{j\cdot }\Vert \asymp \frac{{\gamma }_j}{\Vert \gamma \Vert }\text{for}1\le j\le m.$$

$$\Vert {\hat{U}}_{\cdot 1}-U_{\cdot 1}{C}_U\Vert \le C\frac{\sqrt{\log (nm)\mathcal{Z}}+\log (nm)}{\Vert \theta \Vert \Vert \gamma \Vert },$$

$$\Vert {\hat{U}}_{\cdot 2\sim \kappa }-U_{\cdot 2\sim \kappa }{O}_U{\Vert }_F\le C\frac{\sqrt{\log (nm)\mathcal{Z}}+\log (nm)}{\Vert \theta \Vert \Vert \gamma \Vert },$$

$$\Vert {\hat{V}}_{\cdot 1}-V_{\cdot 1}{C}_V\Vert \le C\frac{\sqrt{\log (nm)\mathcal{Z}}+\log (nm)}{\Vert \theta \Vert \Vert \gamma \Vert },$$

$$\Vert {\hat{V}}_{\cdot 2\sim \kappa }-V_{\cdot 2\sim \kappa }{O}_V{\Vert }_F\le C\frac{\sqrt{\log (nm)\mathcal{Z}}+\log (nm)}{\Vert \theta \Vert \Vert \gamma \Vert },$$

$$R_{i\cdot }^r=\frac{(U_{\cdot 2\sim \kappa }{O}_U{)}_{i\cdot }}{C_U{U}_{i1}}\text{and}{R}_{j\cdot }^c=\frac{(V_{\cdot 2\sim \kappa }{O}_V{)}_{j\cdot }}{C_V{V}_{j1}}.$$

$$R_{i\cdot }^r=\frac{(U_{\cdot 2\sim \kappa }{O}_U{)}_{i\cdot }}{C_U{U}_{i1}}=\frac{U_{i2\sim \kappa }{O}_U}{C_U{U}_{i1}}=\frac{\frac{{\theta }_i}{\Vert {\theta }^{(c_i^r)}\Vert }{Y}_{c_i^{r}2\sim \kappa }{O}_U}{C_U\frac{{\theta }_i}{\Vert {\theta }^{(c_i^r)}\Vert }{Y}_{c_i^{r}1}}=\frac{Y_{c_i^{r}2\sim \kappa }{O}_U}{C_U{Y}_{c_i^{r}1}},$$

$$\Vert R_{i_1\cdot }^r-R_{i_2\cdot }^r\Vert \ge 2,\text{if}{c}_{i_1}^r\ne c_{i_2}^r,\text{and}\Vert R_{i_1\cdot }^r-R_{i_2\cdot }^r\Vert =0,\text{if}{c}_{i_1}^r=c_{i_2}^r;$$

$$\Vert R_{j_1\cdot }^c-R_{j_2\cdot }^c\Vert \ge 2,\text{if}{c}_{j_1}^c\ne c_{j_2}^c,\text{and}\Vert R_{j_1\cdot }^c-R_{j_2\cdot }^c\Vert =0,\text{if}{c}_{j_1}^c=c_{j_2}^c.$$

$$\Vert {\hat{R}}^r-R^r{\Vert }_F^2\le \frac{C{\tau }_n^2{\left\{\sqrt{\log (nm)\mathcal{Z}}+\log (nm)\right\}}^2}{{\theta }_{\min }^2\Vert \gamma {\Vert }^2},$$

$$\Vert {\hat{R}}^c-R^c{\Vert }_F^2\le \frac{C{\tau }_m^2{\left\{\sqrt{\log (nm)\mathcal{Z}}+\log (nm)\right\}}^2}{{\gamma }_{\min }^2\Vert \theta {\Vert }^2},$$

$$\Vert {\hat{X}}^r-R^r{\Vert }_F^2\le \frac{C{\tau }_n^2{\left\{\sqrt{\log (nm)\mathcal{Z}}+\log (nm)\right\}}^2}{{\theta }_{\min }^2\Vert \gamma {\Vert }^2},$$

$$\Vert {\hat{X}}^c-R^c{\Vert }_F^2\le \frac{C{\tau }_m^2{\left\{\sqrt{\log (nm)\mathcal{Z}}+\log (nm)\right\}}^2}{{\gamma }_{\min }^2\Vert \theta {\Vert }^2}.$$

$$|{\mathcal{V}}^r\setminus {\mathcal{W}}^r|<\frac{C{\tau }_n^2{\left\{\sqrt{\log (nm)\mathcal{Z}}+\log (nm)\right\}}^2}{{\theta }_{\min }^2\Vert \gamma {\Vert }^2},$$

$$|{\mathcal{V}}^c\setminus {\mathcal{W}}^c|<\frac{C{\tau }_m^2{\left\{\sqrt{\log (nm)\mathcal{Z}}+\log (nm)\right\}}^2}{{\gamma }_{\min }^2\Vert \theta {\Vert }^2}.$$

$$\underset{n,m\to \infty }{\lim }\frac{|{\mathcal{V}}^r\setminus {\mathcal{W}}^r|}{n}<\underset{n,m\to \infty }{\lim }\frac{C{\log }^2(n)\log (nm)}{n}=0,$$

$$\underset{n,m\to \infty }{\lim }\frac{|{\mathcal{V}}^c\setminus {\mathcal{W}}^c|}{m}<\underset{n,m\to \infty }{\lim }\frac{C{\log }^2(m)\log (nm)}{m}=0.$$

$$B=\left[\begin{array}{ccc}1.0& 0.1& 0.2\\ 0.3& 0.9& 0.1\end{array}\right].$$

$$B=\left[\begin{array}{ccc}0.9& 0.1& 0.2\\ 0.2& 1.0& 0.1\end{array}\right].$$

$$B=\left[\begin{array}{ccc}1.0& 0.2& 0.1\\ 0.2& 0.9& 0.3\end{array}\right].$$

Full text

1 Introduction

Journal citation networks are powerful tools for studying the structure and evolution of science knowledge [31]. In these networks, nodes represent journals and edges represent citation relationships. Specifically, if a paper published in journal $A$ cites a paper published in journal $B$, there is a directed edge from $A$ to $B$. This results in a directed and weighted network, where the weight reflects the citation counts between journals. Journal citation networks have been widely applied in various disciplines, such as economics [35, 29], education [64], medicine [22, 24], neuroscience [20, 59], technology management [39], e-commerce [27], demography [45], materials science [5] and others. Specifically, journal citation networks can be used to identify influential journals [49, 54], uncover knowledge dissemination patterns [12, 53], and measure inter-disciplinarity [23, 41].

In particular, journal citation networks play a crucial role in knowledge source discovery by revealing key journals, researchers, methods, and theories within a field or interdisciplinary research [43, 69]. Existing techniques for identifying knowledge sources using journal citation networks fall into three main categories. The first is the network centrality analysis, which assesses the importance of journals by calculating measures such as degree, closeness, and betweenness centrality [40, 15]. The second is intermediary analysis, which traces citation paths to identify influential journals that serve as bridges [42]. The third is community detection, which clusters journals to uncover topical or disciplinary structures [14, 10, 66]. For instance, [26] identify research clusters based on direct citations, co-citation, and bibliographic coupling to generate a scientific knowledge map in library and information science. Later, [13] perform a comprehensive bibliometric analysis of publications in the business field, employing a label propagation algorithm to identify research themes. However, to the best of our knowledge, although the journal citation network plays a crucial role, its application in the field of statistics is relatively rare.

In this study, we focus on the field of statistics and collect articles from eight core journals (listed in Table 3) that span 2001 to 2023. The data are retrieved from the Web of Science (WOS, www.webofknowledge.com), including 16,119 articles and their corresponding 179,654 cited references. Using this dataset, we construct a bipartite (i.e., two-mode) journal citation network according to the citation relationships, where one mode comprises the eight core journals and the other the cited journals. This bipartite construction explicitly distinguishes the citing and cited roles, thereby preserving the inherent asymmetry of the citation behavior. This representation is particularly suitable for the discovery of the source of knowledge, which can be obscured in traditional one-mode projections [30]. Although prior studies mainly rely on one-mode networks that treat all journals as structurally equivalent [28, 17], the bipartite perspective offers a more accurate structural representation. This formulation choice forms one of our key contributions. Similar bipartite network formulations have been adopted in domains such as recommendation systems [46], ecological networks [55], and financial transaction analysis [16].

To reveal the structural patterns of citation behaviors between journals, community detection can be performed on such bipartite journal citation networks. Existing approaches to bipartite community detection generally fall into three categories. The first is the projection method. It simply converts the bipartite network into one-mode projections and then applies traditional community detection techniques to the one-mode networks. Although this approach is simple to implement, it often results in information loss in the bipartite structure [7, 51, 4]. The second is the heuristic method, which identifies communities using techniques such as matrix decomposition, motif-based clustering, or discrete optimization, without relying on explicit generative models. These methods emphasize computational efficiency and scalability and have been widely applied in computer science domains such as information retrieval, recommendation systems, and web mining [70, 62]. The third is the model-based method, which assumes that the network is generated from an underlying probabilistic model [52, 67]. A widely adopted model in this category is the Bipartite Stochastic Block Model (BiSBM), which assumes that edge probabilities depend solely on the community memberships of the two types of nodes [38]. However, the original BiSBM does not account for degree heterogeneity, which is an empirically common feature of real-world networks [36]. To address this limitation, recent extensions, such as degree-corrected BiSBM, introduce node-specific parameters to accommodate degree variability, improving both empirical performance and theoretical guarantees [34, 8, 9, 72].

Despite these advances, efficient and accurate community detection in large-scale and degree heterogeneous bipartite networks remains a challenge [11]. Traditional inference techniques such as pseudo-likelihood estimation [3, 61] and variational methods [44, 72] are often computationally intensive or sensitive to initialization [71]. As a scalable and initialization-free alternative, spectral clustering relies on the leading eigenvectors of the adjacency or Laplacian matrices. Several recent studies have explored spectral methods for community detection under BiSBM assumptions. For instance, [56] propose a distribution-free model and a degree-corrected extension for weighted bipartite networks. Then they introduce bipartite spectral clustering (BiSC) and normalized BiSC for community detection. These methods are further extended to overlapping bipartite networks in Qing and Wang [57]. Although these approaches represent important progress in model-based community detection, their theoretical guarantees often depend on the boundedness of deviations between adjacency entries and their expectations, a condition that may not hold in many weighted networks. Moreover, although normalized BiSC attempts to address degree heterogeneity, it can perform poorly in sparse settings. The limitation can also be observed in normalized principal component analysis for unipartite spectral clustering [34].

To overcome these challenges, we propose Bipartite Spectral Clustering on Ratios-of-Eigenvectors (Bi-SCORE), a new spectral methodology designed for bipartite degree-corrected block models. Bi-SCORE generalizes the widely used SCORE framework [34] in two key directions: it extends SCORE from unipartite to bipartite networks and from unweighted to weighted networks, thereby largely broadening its applicability. Specifically, Bi-SCORE utilizes the leading singular vectors of the bipartite adjacency matrix and applies a ratio transformation to mitigate the influence of nodal degrees. Unlike iterative optimization-based methods, Bi-SCORE is initialization-free and computationally efficient. Theoretically, Bi-SCORE offers strong consistency guarantees under mild regularity conditions, even in sparse and highly heterogeneous bipartite networks. Practically, we demonstrate the effectiveness of Bi-SCORE by applying it to a bipartite journal citation network in the field of statistics. The results reveal a clear and interpretable six-community structure, where each community corresponds to a meaningful subfield of statistical research, and the citation patterns among core journals further validate the domain-specific focus and cross-community knowledge flow captured by Bi-SCORE.

The remainder of this paper is organized as follows. In Section 2, we construct the journal citation network and provide a descriptive analysis. Section 3 proposes the Bi-SCORE method and presents its theoretical properties. Extensive simulation studies are conducted in Section 4 to evaluate the finite sample performance of our newly proposed methodology. The empirical results of the real data are shown in Section 5. Conclusions and discussions are presented in Section 6. All technical proofs are left in the Appendix.

2 Journal Citation Network

2.1 Data and Descriptive Analysis

This study collects publications from 2001-01-01 to 2023-12-31 in eight representative journals in the field of statistics from the WOS. The journal name, number of publications, and JCR category are listed in Table 3. These eight journals cover the core aspects of statistical science, i.e., statistical methodology, mathematical theory, applied statistics, and computation. Specifically, Biometrika, Journal of Business & Economic Statistics (JBES), Journal of the American Statistical Association (JASA), and Journal of the Royal Statistical Society Series B–Statistical Methodology (JRSS-B) aim to publish high-quality articles on statistical methodology. For instance, JBES focuses on methodological innovations in economics, business, and finance, including the development and improvement of statistical methods, the adaptation of methods from other fields, and advances in computation. Regarding mathematical theory, the Annals of Statistics (AOS) focuses on publishing research papers on the theory of mathematical statistics and applied/interdisciplinary statistics. In the aspect of applied statistics, the Annals of Applied Statistics (AOAS) aims to provide a timely and unified forum for all areas of applied statistics. Furthermore, the Journal of Computational and Graphical Statistics (JCGS) and Statistics and Computing (SC) specialize in computation, with JCGS having a particular focus on the latest techniques of graphical methods in statistics and data analysis. We refer to these eight journals as the core journals. In addition, we thoroughly retrieve the basic information for each publication, including the title, authors, publisher, publication year, keywords, abstract, citation count, and reference list. A concrete example is provided in Table 1. Specifically, the citation count represents the number of times each article has been cited by other works indexed by the Web of Science until 2023. This dataset differs primarily from the data previously collected by our team [25] in that it includes more detailed information about the journals in which the references are published.

To better understand the performance of the eight core journals, we define the following three metrics: average citation, average reference, and self-citation rate. Specifically, the average citation is defined as the total citation counts of all publications in a specific journal over the total number of publications in the same journal. This measurement reflects the influence of publications in a journal. A higher average citation indicates that the research outputs of the journal have greater academic value and impact [47]. The average reference is defined as the total reference counts of all publications in a specific journal over the total number of publications in the same journal. A higher average reference may suggest that publications in a journal rely more heavily on and integrate research from other fields [2]. In addition, the self-citation rate is defined as the total self-citation counts of all publications in a specific journal over the total citation counts of all publications in the same journal, multiplied by 100%. It measures the proportion of citation counts that a journal receives for its own previously published work. Research on the self-citation rate is important in the study of knowledge source discovery, as it often reflects structural differences among journals, such as journal size, publisher country, and subject area [58].

The results are shown in Table 2. Among the eight core journals, JRSS-B has the highest average citation at 113.62, highlighting its influence in the field of statistical methodology. In addition, AOS and JASA also have high average citations, at 79.84 and 65.51, respectively, indicating that these two journals hold a prominent position in statistical theory. Most journals fall within the 30 to 40 range in terms of average reference, suggesting a moderate level of reliance on external sources in their research. There is an obvious variation in the self-citation rate among the eight core journals. Specifically, AOS, Biometrika, and JASA have a relatively higher self-citation rate, all exceeding 10%, while AOAS and SC have a relatively lower self-citation rate, below 3%. This difference may be attributed to the fact that AOS, Biometrika, and JASA concentrate on theoretical and methodological aspects of statistics, and publications in these journals often build upon prior work in the same journal. As a result, publications in these journals are more likely to cite earlier research from the same journal. In contrast, AOAS and SC focus on applied statistics and statistical computing, fields that are more interdisciplinary in nature, thus requiring the citation of research from a broader range of disciplines. The second column of Table 2 represents the total number of publications in the eight core journals. It can be seen that JASA and AOS have higher volume of publication compared to other core journals, demonstrating stronger academic influence. These findings emphasize that the eight core journals have a high academic reputation and authoritative status in statistics and related fields.

2.2 The Journal Citation Network

To investigate the knowledge source of the eight core journals, we construct a bipartite journal citation network. Specifically, we denote the eight core journals as one set of nodes, indexed by $i=1,\dots,8$. Meanwhile, the journals cited by the core journals are denoted as the other set of nodes, indexed by $j=1,\dots,m$, where $m$ is the total number of journals including the core journals. As a result, the adjacency matrix of the journal citation network can be defined as $A=(A_{ij})\in\mathbb{R}^{8\times m}$, where $A_{ij}$ is a non-negative integer representing the frequency with which the core journal $i$ cites papers published in journal $j$. To focus more on key sources of knowledge and ensure a more concentrated analysis, we retain only journals with $A_{ij}\geq 40$, resulting in $m=334$. Figure 1 visualizes the journal citation network, where the core journals are colored orange. The thickness of the edges represents the citation counts between journals, and we only display about 25% names of all journals randomly, for the purpose of image clarity. It shows a clear community structure around the core journals, implying the need for community detection. For instance, the journals surrounding JBES are concentrated in the fields of economics and finance (e.g., Journal of Econometrics, Quantitative Economics), focusing primarily on statistical methods related to economic models. Similarly, the journals surrounding AOS are primarily associated with mathematical theory (e.g., Annals of Mathematics, Probability Surveys), reflecting its emphasis on probability theory, statistical inference, and other aspects of mathematical statistics. The above analysis provides a preliminary view of the community structure of the journal citation network.

Table 4 further lists the top 10 reference journals of the eight core journals. It can be seen that JASA, AOS, Biometrika and JRSS-B rank among almost the top five in terms of references in all core journals. This may be due to the high quality of their articles, making them an important knowledge source for other core journals during the observed period, as well as its comprehensive nature in statistical research and applications. In addition, Table 4 reflects the reference preferences of the core journals, which is in line with the journal scope as presented in Section 2.1. For instance, AOAS tends to cite articles related to applied statistics (e.g., Biometrics, Statistics in Medicine), while also referencing general science journals such as Proceedings of the National Academy of Sciences of the United States of America (PNAS) and Nature. This reflects its interdisciplinary orientation and broad scope of application. AOS primarily cites methodological and theoretical journals (e.g., Bernoulli, IEEE Transactions on Information Theory (TIT), Annals of Mathematical Statistics (AMS)), indicating its emphasis on foundational statistical and mathematical research. Biometrika tends to cite journals related to biostatistics and biomedical research (e.g., Biometrics, Statistics in Medicine). This indicates a concentrated emphasis on methodological development for biological and medical applications. JRSS-B and JASA cite journals across multiple statistical domains, including biostatistics (e.g., Biometrics), econometrics (e.g., Econometrica), and machine learning (e.g., Journal of Machine Learning Research (JMLR)). This reflects their integrative scope, including theoretical, applied, and computational statistics, as well as other interdisciplinary research. JCGS and SC tend to cite journals in computational statistics, such as Computational Statistics & Data Analysis (CSDA) and JMLR. This reflects their focus on statistical computing, data analysis, and the development of algorithmic methodologies. JBES cites more studies related to economics and finance (e.g., Journal of Econometrics (JOE), Econometrica). This reflects its emphasis on the use of statistical methods to address problems in economic and financial research.

In addition to characterizing reference patterns of each eight core journals, we then shift our attention to the journals being cited, evaluating their impact and importance within the citation network. In-degree is widely used to measure this importance, reflecting how frequently a journal is cited by others. Specifically, we calculate the weighted in-degree for each cited journal $j$ as $C_{w}(j)=\sum_{i=1}^{8}A_{ij}$, where $A_{ij}$ is the weight of the edge from core journal $i$ to journal $j$ [1]. The detailed definition of $A_{ij}$ is provided in Section 3.1. Figure 2 presents the distribution of the weighted in-degree. The histogram indicates a highly right-skewed distribution, where the majority of nodes exhibit low-weighted in-degree (between 0 and 1,500), while a few nodes possess extremely high centrality values. Notably, AOS, JASA, JRSS-B, and Biometrika exhibit the highest weighted in-degree. This phenomenon is common in journal citation networks because the right-skewed distribution is often associated with the preferential attachment mechanism, in which newly published articles tend to cite already well-established, high-impact journals [63]. Furthermore, structural centralization in academic disciplines reinforces this imbalance, making core journals crucial in academic dissemination and knowledge accumulation.

3 Methodology

In this section, we first introduce the Weighted Bipartite Degree-Corrected Stochastic Block Model (DCBM), which serves as the underlying generative model for the networks under study. Subsequently, we present Bi-SCORE, our novel spectral clustering algorithm designed to effectively identify community structures within these networks. Finally, we establish the theoretical properties of Bi-SCORE, providing rigorous guarantees on its performance.

3.1 Weighted Bipartite DCBM

The stochastic block model (SBM) is a generative model designed for networks with community structures. In its most basic form, it applies to undirected, unweighted, and single-mode networks, where each node is assigned to one of the $K$ communities. The probability of forming an edge between two nodes is entirely dependent on their community memberships. However, the basic SBM ignores the inherent degree heterogeneity observed in real-world networks. To address this limitation, [36] introduce the degree-corrected stochastic block model (DCBM), which accounts for degree heterogeneity, making it more suitable for analyzing empirical networks. Recently, [72] extend the DCBM to weighted bipartite networks. Consider a network composed of $n$ row nodes, denoted by $\mathcal{V}^{r}=\{1,\dots,n\}$, and $m$ column nodes, denoted by $\mathcal{V}^{c}=\{1,\dots,m\}$. Let $A=(A_{ij})\in\mathbb{R}^{n\times m}$ be the adjacency matrix, where $A_{ij}$ denotes a non-negative integer weight between the $i$-th row node and the $j$-th column node. Assume that the $n$ row nodes are partitioned into $K$ communities, with their community labels given by $c^{r}=(c^{r}_{1},\dots,c^{r}_{n})^{\top}\in\{1,\dots,K\}^{n}$. Similarly, the $m$ column nodes are classified into $L$ communities, with labels $c^{c}=(c^{c}_{1},\dots,c^{c}_{m})^{\top}\in\{1,\dots,L\}^{m}$. To capture degree heterogeneity, we further incorporate two sets of parameters $\theta=(\theta_{1},\dots,\theta_{n})^{\top}\in\mathbb{R}^{n}$ for row nodes and $\gamma=(\gamma_{1},\dots,\gamma_{m})^{\top}\in\mathbb{R}^{m}$ for column nodes. We call them the heterogeneity parameter. These parameters model the individual propensities of nodes to form edges. Conditional on the community labels, the edge weights $A_{ij},\ i=1,\dots,n,\ j=1,\dots,m$ are assumed to be independent Poisson random variables with expected values defined as

[TABLE]

where $B=(B_{c^{r}_{i}c^{c}_{j}})\in\mathbb{R}^{K\times L}$ is a full rank matrix that quantifies the interaction strengths between row and column communities. The community memberships can be represented by indicator matrices: $Z=(Z_{ik})\in\mathbb{R}^{n\times K}$ for row nodes, where $Z_{ik}=1$ if $c^{r}_{i}=k$ and $Z_{ik}=0$ otherwise; and $W=(W_{jl})\in\mathbb{R}^{m\times L}$ for column nodes, defined analogously. Throughout this paper, each row (column) node belongs to only one row (column) community. In other words, each row of $Z$ or $W$ has only one 1. Additionally, we assume that $\text{rank}(Z)=K$ and $\text{rank}(W)=L$, which ensures that every row and column community is non-empty (i.e., contains at least one node). In matrix form, the expected adjacency matrix $\Omega=\mathbb{E}(A)$ can be expressed as

[TABLE]

where $\Theta$ is an $n\times n$ diagonal matrix with $\Theta_{ii}=\theta_{i}$ and $\Gamma$ is an $m\times m$ diagonal matrix with $\Gamma_{jj}=\gamma_{j}$.

3.2 Algorithm

Before introducing the algorithm, we first define some notation. For a matrix $Q$, we use $Q_{i\cdot}$ to denote the $i$-th row of the matrix $Q$, and $Q_{\cdot j}$ to denote the $j$-th column of the matrix $Q$. In addition, we use $Q_{ij}$ to denote the $(i,j)$th entry in the matrix $Q$. For integer $0<j_{1}\leq j_{2}$, let $Q_{\cdot j_{1}\sim j_{2}}$ denote the matrix that is formed by extracting the $j_{1}$-th to $j_{2}$-th columns of the matrix $Q$, and $Q_{ij_{1}\sim j_{2}}$ denote the vector that is formed by extracting the $i$-th row and the $j_{1}$-th to $j_{2}$-th columns of the matrix $Q$. For a vector $u$ and fixed $q>0$, $\|u\|_{q}$ denote its $\ell_{q}$-norm. We remove the subscript if $q=2$. For a matrix $Q$, $\|Q\|$ denotes the spectral norm, and $\|Q\|_{F}$ denotes the Frobenius norm. We let $\|Q\|_{\min}$ denote the smallest singular value of the matrix $Q$. Let $\sigma_{k}(Q)$ denote the $k$-th largest singular value of matrix $Q$, and $\lambda_{k}(Q)$ denote the $k$-th largest eigenvalue of matrix $Q$ ordered by magnitude. For two positive sequences $\{a_{n}\}_{n=1}^{\infty}$ and $\{b_{n}\}_{n=1}^{\infty}$, we say $a_{n}\asymp b_{n}$ if there exists a positive constant $C$ such that $b_{n}/C\leq a_{n}\leq Cb_{n}$ for sufficiently large $n$, i.e., $a_{n}$ and $b_{n}$ are of the same order. For a set $\mathcal{V}$, $|\mathcal{V}|$ denotes its cardinality.

In order to detect meaningful communities in weighted bipartite networks, we develop an algorithm called Bi-SCORE. Specifically, Bi-SCORE utilizes spectral clustering techniques for bipartite networks, based on singular value decomposition (SVD). It takes as input a bipartite adjacency matrix $A\in\mathbb{R}^{n\times m}$ and requires the number of clusters $K$ and $L$ to be pre-specified. The algorithm then outputs the community detection results for both the row and column nodes. The complete details of the Bi-SCORE algorithm are shown in Algorithm 1.

Remark 1**.**

The main goal of the Bi-SCORE algorithm is to find community labels for row nodes and column nodes. We do not directly estimate the heterogeneity parameters ($\theta$ and $\gamma$) or the community interaction strength matrix ($B$). This keeps the algorithm simpler and faster because we avoid estimating many parameters. If the estimation of $\theta$, $\gamma$, and $B$ is required, their values can be further estimated using methods such as moment estimation, after the community partitions have been obtained through the Bi-SCORE algorithm.

3.3 Theoretical Property

Before obtaining the bounds of the misclustered nodes for Bi-SCORE, we first introduce some assumptions and propositions. We begin with an assumption regarding the matrix $B$.

Assumption 1**.**

The matrix $B$ satisfies

[TABLE]

This assumption restricts the connection strengths $B_{kl}$ to lie within the interval $[0,1]$, where the upper bound of 1 is primarily for notational simplicity. In practice, this condition can be relaxed as long as the entries of $B$ remain bounded. The irreducibility and non-negativity of $BB^{\top}$ and $B^{\top}B$ ensure that the leading left and right singular vectors of $B$ (associated with the largest singular value) are non-zero. This is crucial for the Bi-SCORE algorithm, as it guarantees non-vanishing denominators in spectral analysis, thereby enabling reliable community detection. This assumption is similar to Condition (2.7) of [34] and Assumption 1 of [65].

To introduce our assumptions about the heterogeneity parameters, we first define some notation. For each $k=1,\dots,K$ and $l=1,\dots,L$, define $\theta^{(k)}\in\mathbb{R}^{n}$ and $\gamma^{(l)}\in\mathbb{R}^{m}$ by

[TABLE]

Next, define the following quantities: $\theta_{\min}\equiv\min_{1\leq i\leq n}\theta_{i}$, $\theta_{\max}\equiv\max_{1\leq i\leq n}\theta_{i}$, $\gamma_{\min}\equiv\min_{1\leq j\leq m}\gamma_{j}$, and $\gamma_{\max}\equiv\max_{1\leq j\leq m}\gamma_{j}$. We also define

[TABLE]

a quantity that frequently appears in our analysis. In this work, we allow the heterogeneity parameters $\theta$ and $\gamma$ to scale with the network sizes $n$ and $m$.

Assumption 2**.**

The heterogeneity parameters $\theta$ and $\gamma$ satisfy

[TABLE]

Condition (5) ensures that every node has a strictly positive probability of forming connections. The requirement in (6) enforces a form of balance between communities: the total connectivity within each row (or column) community is comparable to that of any other, preventing disproportionate dominance. Finally, condition (7) implies that an accurate community recovery requires the total connection strength, represented by $\|\theta\|_{1}$ and $\|\gamma\|_{1}$, to grow faster than $\log(nm)$. Similar assumptions can be found in Conditions (2.6), (2.9), (2.13) of [34] and Assumption 2 of [65]. We then present a few properties that follow directly from Assumption 2. As a direct consequence of (6), we have:

[TABLE]

Moreover, since $\theta_{\min}\|\theta\|_{1}\leq\|\theta\|^{2}$ and $\gamma_{\min}\|\gamma\|_{1}\leq\|\gamma\|^{2}$, condition (7) implies:

[TABLE]

Next, we analyze the singular vector matrix of the expected adjacency matrix $\Omega$ defined in (2), which captures the essential structure for clustering. We then show that the properties of $\Omega$ closely approximate those of the observed random adjacency matrix $A$. To proceed, we define $S\equiv\Psi_{\theta}B\Psi_{\gamma}^{\top}$, where $\Psi_{\theta}$ is a $K\times K$ diagonal matrix with entries $(\Psi_{\theta})_{kk}=\|\theta^{(k)}\|/\|\theta\|$, and $\Psi_{\gamma}$ is an $L\times L$ diagonal matrix with entries $(\Psi_{\gamma})_{ll}=\|\gamma^{(l)}\|/\|\gamma\|$. It captures the core community-level interactions scaled by the relative sizes of the heterogeneous parameters.

Proposition 1**.**

Let $\Omega=U\Lambda_{\Omega}V^{\top}$ be the compact singular value decomposition of the expected adjacency matrix $\Omega$. Then, the singular values of $\Omega$ are given by

[TABLE]

where $S\equiv\Psi_{\theta}B\Psi_{\gamma}^{\top}$. Let $S=Y\Lambda_{S}H^{\top}$ denote the compact singular value decomposition of $S$. The singular vectors of $\Omega$ can be expressed in row form as

[TABLE]

and in column form as

[TABLE]

Furthermore, under Assumption 2, we have

[TABLE]

The proof of Proposition 1 can be found in Appendix A.1. Proposition 1 provides a precise characterization of the spectral structure of the expected adjacency matrix $\Omega$. Specifically, it states that $\Omega$ has rank $\kappa$, where $\kappa$ is the rank of the matrix $B$. As a result, $\Omega$ admits a compact singular value decomposition $\Omega=U\Lambda_{\Omega}V^{\top}$, where $U\in\mathbb{R}^{n\times\kappa}$, $V\in\mathbb{R}^{m\times\kappa}$, and $\Lambda_{\Omega}\in\mathbb{R}^{\kappa\times\kappa}$ contains the nonzero singular values. Moreover, the singular vectors of $\Omega$ inherit a structured form that aligns with the community memberships. In particular, each row of $U$ (respectively, $V$) lies in the direction of the left (respectively, right) singular vector of the corresponding community from the decomposition of $S$, scaled by the heterogeneity parameter $\theta_{i}$ (or $\gamma_{j}$). This implies that nodes within the same community have vector representations that are parallel but differ in magnitude. The column-wise form of the singular vectors further reveals that each column of $U$ and $V$ can be expressed as a linear combination of community-specific vectors, weighted by normalized heterogeneity. Finally, (13) illustrates that the spectral embedding magnitudes are governed by the connection tendencies of nodes.

The next proposition establishes bounds on the deviation between the singular vectors of the random adjacency matrix $A$ and those of its expectation $\Omega$.

Proposition 2**.**

Let $\widehat{U}_{\cdot 1},\dots,\widehat{U}_{\cdot\kappa}$ and $\widehat{V}_{\cdot 1},\dots,\widehat{V}_{\cdot\kappa}$ denote the top $\kappa$ left and right singular vectors of $A$, respectively, and let $U_{\cdot 1},\dots,U_{\cdot\kappa}$ and $V_{\cdot 1},\dots,V_{\cdot\kappa}$ be the corresponding singular vectors of $\Omega$. Under Assumptions 1 and 2, there exist constants $C_{U}$ and $C_{V}$ with absolute value 1, as well as two orthogonal $(\kappa-1)\times(\kappa-1)$ matrices $O_{U}$ and $O_{V}$, such that for sufficiently large $n$ and $m$, with probability at least $1-1/n-1/m$, the following bounds hold

[TABLE]

where $\mathcal{Z}$ is defined in (4).

The proof of Proposition 2 can be found in Appendix A.2. Proposition 2 quantifies how closely the leading singular vectors of the observed adjacency matrix $A$ approximate those of its expectation $\Omega$. Specifically, the leading singular vectors (the first columns of $U$ and $V$) are determined only up to a sign, which is captured by the constants $C_{U}$ and $C_{V}$ of modulus one. For the remaining $\kappa-1$ singular vectors, the ambiguity is up to an orthogonal rotation, captured by the matrices $O_{U}$ and $O_{V}$.

Based on the expression in (11) and the approximation result from Proposition 2, we are now ready to explain the rationale behind removing the influence of node heterogeneity from the singular vectors in Algorithm 1. The core idea is to normalize each row of the singular vector matrix by its first entry. Specifically, for each $i=1,\dots,n$ and $j=1,\dots,m$, we define the ratio matrices $R^{r}$ and $R^{c}$ as

[TABLE]

In other words, each row of $U$ is divided by its first entry, and the second through $\kappa$th columns are collected to form the ratio matrix $R^{r}$. The construction of $R^{c}$ follows analogously for the right singular vectors. To see the effect of this transformation, observe that

[TABLE]

where the third step follows from (12). As shown in (14), the ratio eliminates the heterogeneity parameter $\theta_{i}$, and the resulting row vector in $R^{r}$ depends only on the community label $c_{i}^{r}$. Therefore, all nodes belonging to the same row community share the same row in $R^{r}$. A similar argument holds for $R^{c}$ and the column communities. This analysis highlights the key role of the ratio transformation in Algorithm 1. To understand how the ratio matrices enable clear separation between communities, we present the following proposition.

Proposition 3**.**

Let $R^{r}$ and $R^{c}$ be the ratio matrices constructed from the singular vectors of the expected adjacency matrix $\Omega$. Then, under Assumptions 1 and 2, for all $1\leq i_{1}<i_{2}\leq n$ and $1\leq j_{1}<j_{2}\leq m$, the following properties hold:

[TABLE]

The proof of Proposition 3 can be found in Appendix A.3. Proposition 3 indicates that if two nodes $i_{1}$ and $i_{2}$ (or $j_{1}$ and $j_{2}$) belong to different communities (i.e., $c^{r}_{i_{1}}\neq c^{r}_{i_{2}}$ or $c^{c}_{j_{1}}\neq c^{c}_{j_{2}}$), then their rows in $R^{r}$ or $R^{c}$ will be sufficiently distinct. On the other hand, if two nodes belong to the same community (i.e., $c^{r}_{i_{1}}=c^{r}_{i_{2}}$ or $c^{c}_{j_{1}}=c^{c}_{j_{2}}$), their rows will be identical. This result suggests that the ratio matrices effectively separate nodes from different communities while clustering those from the same community together. Therefore, the structure captured in the ratio matrices $R^{r}$ and $R^{c}$ directly reflects the underlying community structure, which is crucial for the community detection task in Algorithm 1. To assess how well the empirical ratio matrices approximate their population counterparts, we present the following result.

Proposition 4**.**

For sufficient large $n$ and $m$, and under Assumptions 1 and 2, with probability at least $1-1/n-1/m$, the following inequalities hold:

[TABLE]

where $\tau_{n}$ and $\tau_{m}$ are threshold values defined in Algorithm 1.

The proof of Proposition 4 can be found in Appendix A.4. In the following proposition, we provide a bound on the difference between the empirical clustering matrices $\widehat{X}^{r}$ and $\widehat{X}^{c}$ (derived from Algorithm 1) and $R^{r}$ and $R^{c}$. Specifically, $\widehat{X}^{r}$ (or $\widehat{X}^{c}$) is the matrix with exactly $K$ (or $L$) different rows (or columns), and it is the one closest to the ratio matrix $\widehat{R}^{r}$ (or $\widehat{R}^{c}$) in terms of the Frobenius norm.

Proposition 5**.**

For sufficiently large $n$ and $m$, under Assumptions 1 and 2, with probability at least $1-1/n-1/m$, the following inequalities hold:

[TABLE]

The proof of Proposition 5 can be found in Appendix A.5. Proposition 5 suggests that the empirical clustering matrices will closely match the true community structure, which is essential for effective clustering in the Bi-SCORE algorithm.

Before presenting the main theorem for the Bi-SCORE algorithm, we first introduce some notation. Recall that $\mathcal{V}^{r}$ and $\mathcal{V}^{c}$ are the sets of all row and column nodes, respectively. We define the sets of “well-behaved” nodes as follows: $\mathcal{W}^{r}\equiv\{1\leq i\leq n:\|\widehat{X}_{i\cdot}^{r}-R_{i\cdot}^{r}\|\leq 1/2\}$ and $\mathcal{W}^{c}\equiv\{1\leq j\leq m:\|\widehat{X}_{j\cdot}^{c}-R_{j\cdot}^{c}\|\leq 1/2\}$. This specific threshold of $1/2$ is chosen to ensure that the estimated values $\widehat{X}^{r}$ and $\widehat{X}^{c}$ for nodes from different communities remain sufficiently separated, which is crucial for accurate clustering. While any constant within the interval $(0,1)$ could serve as a valid threshold, we adopt $1/2$ in line with the practice established in [65]. By this definition, the sets $\mathcal{V}^{r}\setminus\mathcal{W}^{r}$ and $\mathcal{V}^{c}\setminus\mathcal{W}^{c}$ thus consist of nodes that are “ill-behaved”. Let $n_{k}$ denote the number of row nodes in community $k$, for $k=1,\dots,K$, and $m_{l}$ denote the number of column nodes in community $l$, for $l=1,\dots,L$. The following theorem provides upper bounds on the number of misclustered row and column nodes under mild conditions.

Theorem 1**.**

Under Assumptions 1 and 2, suppose that $\lvert\mathcal{V}^{r}\setminus\mathcal{W}^{r}\rvert<\min\{n_{1},\dots,n_{K}\}$ and $\lvert\mathcal{V}^{c}\setminus\mathcal{W}^{c}\rvert<\min\{m_{1},\dots,m_{L}\}$. Then, all nodes in $\mathcal{W}^{r}$ and $\mathcal{W}^{c}$ are correctly clustered by the Bi-SCORE algorithm. Moreover, for $n$ and $m$ large enough, with probability at least $1-1/n-1/m$,

[TABLE]

The proof of Theorem 1 can be found in Appendix A.6. In Theorem 1, the assumption that the number of misclustered nodes is less than the size of the smallest community guarantees that at least one node of each community is correctly clustered. A similar assumption is made in [34] and [65]. In the special case where the row nodes coincide with the column nodes, the network is simplified to a one-mode graph. When $\mathcal{Z}$ is of the order larger than $\log(nm)$, the bound on the number of misclustered nodes aligns with the result presented in Theorem 1 of [65]. Next, consider a simple case where $\theta$ and $\gamma$ are bounded by constants, i.e., $0<\alpha\leq\min(\theta_{\min},\gamma_{\min})\leq\max(\theta_{\max},\gamma_{\max})\leq\beta<\infty$ with $\alpha$ and $\beta$ being positive constants and $\alpha\leq\beta$. For sufficiently large $n$ and $m$, we have $\mathcal{Z}\asymp\max(n,m)$, which dominates $\sqrt{\log(nm)}$. In this scenario, it follows that: $\mathcal{Z}/(\theta_{\min}^{2}\|\gamma\|^{2})\leq\beta^{2}\alpha^{-4}$ and $\mathcal{Z}/(\gamma_{\min}^{2}\|\theta\|^{2})\leq\beta^{2}\alpha^{-4}$. When $\tau_{n}=\log(n)$ and $\tau_{m}=\log(m)$ with $n\asymp m$, the misclustering rates satisfy:

[TABLE]

4 Simulation

4.1 Model Setup

To evaluate the finite sample performance of our newly proposed methodology, we conduct a number of simulation studies in this section. We adopt the weighted bipartite DCBM described in Section 3.1 to generate bipartite networks. First, we assign each of the $n$ row nodes and the $m$ column nodes to one of the $K$ and $L$ communities, respectively. Specifically, we independently sample the community label of each row node from the discrete uniform distribution over $\{1,\dots,K\}$, and similarly for each column node over $\{1,\dots,L\}$. The resulting community memberships are recorded in matrices $Z\in\mathbb{R}^{n\times K}$ and $W\in\mathbb{R}^{m\times L}$. Second, we specify the community-level interaction matrix $B\in\mathbb{R}^{K\times L}$, which controls the baseline connectivity strength between communities. The details of generating $B$ are described in each scenario. Third, we incorporate two degree parameters as $\theta_{i}=\sqrt{\rho}\cdot a_{i}$ for the row nodes and $\gamma_{j}=\sqrt{\rho}\cdot b_{j}$ for the column nodes, where $a_{i},b_{j}\sim\mathrm{Uniform}(0.5,1)$ with $i=1,\dots,n$ and $j=1,\dots,m$. The parameter $\rho\in(0,1]$ controls the overall strength of degree heterogeneity in the network. Specifically, higher values of $\rho$ lead to greater variability in nodal degrees, whereas smaller values result in more homogeneous degree distributions. As a result, the adjacency matrix $A\in\mathbb{R}^{n\times m}$ is generated with independent Poisson random variables whose means are specified by equation (1).

Scenario 1: Balanced Sample Size

In this scenario, we investigate the clustering performance under balanced sample size. Specifically, the combinations $(n,m)$ are chosen as $(500,525)$, $(1000,1050)$, $(1500,1575)$, $(2000,2100)$, $(2500,2625)$, and $(3000,3150)$. In addition, the degree heterogeneity parameter is set to $\rho=0.2$, and the number of communities is set to $K=2$ and $L=3$. Finally, the interaction matrix $B\in\mathbb{R}^{2\times 3}$ is specified as

[TABLE]

Scenario 2: Unbalanced Sample Size

In this scenario, we evaluate the effect of an unbalanced sample size between the two types of nodes. In particular, the combinations $(n,m)$ are chosen as $(50,1500)$, $(100,3000)$, $(150,4500)$, $(200,6000)$, $(250,7500)$, and $(300,9000)$. In addition, we set the degree heterogeneity parameter to $\rho=0.9$, and the number of communities is chosen the same as that in Scenario 1. Finally, the interaction matrix $B\in\mathbb{R}^{2\times 3}$ is specified as

[TABLE]

This setting mimics some bipartite networks in the real world, such as in academic collaboration networks, where one type of node is typically much smaller in number than the other [50, 73]. Specifically, the journal citation network constructed in our study also exhibits the same phenomenon.

Scenario 3: Degree Heterogeneity

In this scenario, we examine the influence of degree heterogeneity on clustering performance. We fix the number of row and column nodes at $(n,m)$ = $(1000,1050)$. In addition, the degree heterogeneity parameter $\rho$ is varied in $\{0.1,0.2,\ldots,1.0\}$. As $\rho$ increases, nodal degrees become more heterogeneous, raising greater challenges for community detection [18, 37]. Finally, the number of communities remains the same as in Scenario 1, and the interaction matrix $B$ is specified as

[TABLE]

Scenario 4: Number of Communities

In this scenario, we investigate the effect of varying the number of communities on clustering performance. We fix the number of nodes in the row and column at $(n,m)=(1000,1050)$, and set the degree heterogeneity parameter to $\rho=0.9$. Specifically, we consider three combinations of community numbers $(K=2,L=3)$, $(K=3,L=5)$, and $(K=3,L=7)$. The corresponding community interaction matrices $B_{1}$, $B_{2}$, and $B_{3}$ are given as

[TABLE]

4.2 Performance Measurement and Results

To assess the performance of different clustering methods, we conduct simulation experiments with $T=1{,}000$ replications for each scenario. In the $t$th replication, the clustering performance is evaluated using the error rate (ErrorRate) and the adjusted Rand index (ARI) [33], which are defined as follows.

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

In the above equations, $n_{kl}^{r(t)}$ denotes the number of row nodes that belong to cluster $k$ in the true partition and are assigned to cluster $l$ in the predicted partition in the $t$th replication, where $1\leq k,l\leq K$. The terms $n_{k\cdot}^{r(t)}=\sum_{l}n_{kl}^{r(t)}$ and $n_{\cdot l}^{r(t)}=\sum_{k}n_{kl}^{r(t)}$ represent the marginal totals of the row and column, respectively. The quantities $n_{kl}^{c(t)}$, $n_{k\cdot}^{c(t)}$, and $n_{\cdot l}^{c(t)}$ for column nodes are defined analogously, with $1\leq k,l\leq L$. Specifically, the adjusted Rand index is bounded above by $1$ and equals [math] under random labeling, which higher values indicate better agreement [60]. The final error rate and adjusted Rand index are taken as the averages over the $1{,}000$ replications. For each setting, we first generate the bipartite network according to the model described in Section 3.1, and extract the giant component if necessary. Subsequently, we apply Bi-SCORE, nBiSC [56], and spectral clustering [19]. The results are summarized in Figure 3.

Figure 3 displays the error rate curves in four scenarios. In general, Bi-SCORE performs better than nBiSC and spectral clustering in terms of the error rate. We observe the following three phenomena. First, in both Scenario 1 and Scenario 2, the error rates of Bi-SCORE and nBiSC decrease notably as the sample size increases, while the performance of spectral clustering remains relatively unchanged. Specifically, Bi-SCORE achieves the lowest error rates and exhibits the fastest rate of improvement with increasing sample size. Second, in Scenario 3, increasing $\rho$ leads to greater degree heterogeneity in the network, which in turn improves the performance of both Bi-SCORE and nBiSC. In particular, when $\rho>0.2$, Bi-SCORE outperforms the other two methods, highlighting its adaptability to heterogeneous network structures. Third, in Scenario 4, the performance of all methods declines as the network sparsity increases. Although the advantage of Bi-SCORE over nBiSC becomes less pronounced under this setting, Bi-SCORE still maintains a slight advantage, especially in less complex cases.

Figure 4 displays the adjusted Rand index curves in four scenarios. In general, Bi-SCORE achieves higher adjusted Rand indices than both nBiSC and spectral clustering. We observe the following three phenomena. First, in Scenario 1 and Scenario 2, the adjusted Rand indices of Bi-SCORE and nBiSC increase substantially as the sample size grows, while the adjusted Rand index of spectral clustering remains relatively lower. In particular, Bi-SCORE attains the highest adjusted Rand index and demonstrates the fastest rate of improvement with increasing sample size. Second, in Scenario 3, both Bi-SCORE and nBiSC achieve higher adjusted Rand indices as the degree heterogeneity parameter $\rho$ increases. Specifically, Bi-SCORE achieves the highest adjusted Rand index among all three methods. Third, in Scenario 4, the adjusted Rand indices for all methods decrease as the network sparsity increases. Although the advantage of Bi-SCORE over nBiSC becomes less pronounced in this setting, Bi-SCORE still maintains a slight advantage, especially in less complex cases.

5 Knowledge Discovery of Statistical Journals

5.1 Knowledge Discovery of the Journal Citation Network

In this section, we apply the newly proposed Bi-SCORE to our journal citation network described in Section 2.2. Since our goal is to conduct knowledge source discovery in the eight core statistical journals, we only need to identify the community structure among the column nodes, i.e., cited journals. To this end, we first need to decide the number of communities $L$. We apply Bi-SCORE for $L=\{3,4,5,6,7,8\}$, and compare the clustering results for each $L$. The results show that the community structure is clearest when $L=6$. Therefore, we set $L=6$ in the subsequent analysis. The result of the community structure is further visualized in Figure 5. It reveals that the journal citation network is divided into 6 communities, which are summarized in detail in the following.

• •

Applied Statistics (Purple). This community contains 98 journals, which primarily focus on the practical application of statistical methods. Research in this community emphasizes real-world data analysis, interdisciplinary collaboration, and methodological adaptation. The most representative journal in this community is AOAS, as expected for a core journal. In addition, Bioinformatics, Nature, Neuroscience, and American Journal of Epidemiology (AJE) also have high weighted in-degree within this community.

• •

Statistical Methodology (Cyan). This community contains 73 journals, which primarily focus on methodological development and innovation in statistical science. Research topics include new statistical methods, model selection, and statistical inference. The main journals in this community are Biometrika, JASA, and JRSS-B (core journals). Furthermore, Biometrics, Econometrica, and Statistica Sinica also exhibit high weighted in-degree in this community.

• •

Computational Statistics (Blue). This community contains 50 journals, which emphasize computational methodology, statistical software, and data science techniques. These features highlight the central role of algorithmic and computational innovation in modern statistical research, bridging traditional statistical modeling with emerging data science applications. Specifically, JCGS, SC, CSDA, and Journal of Statistical Software have the highest weighted in-degree.

• •

Mathematical Statistics (Yellow). This community contains 45 journals, which focus primarily on mathematical analysis, probability theory, and statistical inference. Research in this community emphasizes theoretical innovation and foundational contributions. The leading journal in this community is AOS, which is the core journal defined in our journal citation network. In addition, Annals of Probability, Probability Theory and Related Fields, and Mathematical Methods of Statistics also have high weighted in-degree in this community.

• •

Econometrics and Business Statistics (Dark Green). This community contains 39 journals, with a primary focus on econometrics, finance, and quantitative economics. These citation patterns reflect a close integration of econometric methods with empirical research in economics and finance. Notably, JBES, Journal of Applied Econometrics (JAE), American Economic Review (AER), and Journal of Finance (JF) have the highest weighted in-degree within this community.

• •

Others (Red). This community contains 29 journals spanning a wide range of fields, including education, psychology, medicine, biology, economics, political science, engineering, and others. Research in this group emphasizes the application and advancement of statistical methods in emerging interdisciplinary fields. Representative journals in this community include Journal of Educational Psychology (JEP), Journal of Forecasting, Journal of Human Resources (JHR), and Journal of the Royal Statistical Society Series A (JRSS-A).

Figure 6 reveals the knowledge sources for the core journals in the 6 identified communities. The following three main phenomena can be observed. First, all core journals extensively refer to the “Statistical Methodology” community. This suggests the central role of methodological development, which connects different fields of statistics from theory to application. In contrast, all core journals rarely cite from the “Others” community, indicating a concentration of knowledge sources within statistics-related domains. Second, both AOAS and JASA cite heavily from the “Applied Statistics” community, reflecting their strong emphasis on empirical research and real-world applications. This may be due to the special section Journal of the American Statistical Association: Applications and Case Studies (JASA ACS), which is devoted to applied statistical analyses of real data. Specifically, JASA also cites extensively from the “Statistical Methodology” community, highlighting its integrated role in both methodological development and applied research. Third, SC and JCGS cite frequently from the “Computational Statistics” community. This reflects their emphasis on computational tools, graphical methods, and data-driven applications in diverse scientific domains. Moreover, AOS cites heavily from the “Mathematical Statistics” community, and JBES cites heavily from the “Econometrics and Business Statistics” community. This reveals the distinct orientations of the two journals: AOS emphasizes theoretical innovation in probability and statistical inference, while JBES focuses on the development and application of statistical methods in economics, finance, and business.

5.2 Knowledge Discovery of the Applied Statistics Community

In Section 5.1, we analyze the knowledge sources of our journal citation network and divide it into 6 communities. In addition, as the largest community in the journal citation network (98 journals), the “Applied Statistics” community contains the core journal AOAS and many other representative journals. To further investigate the knowledge sources of this community, we apply Bi-SCORE to its subnetwork. Similarly, we apply Bi-SCORE for $L=\{2,3,4,5\}$ and compare the clustering results. Among these alternatives, the partition obtained with $L=4$ provides the most coherent and interpretable community division. Consequently, $L=4$ is adopted in the subsequent investigation. The result of the community structure is visualized in Figure 7, and its characteristics are discussed in detail below.

• •

Interdisciplinary Research (Pink). This subcommunity contains 28 journals, where statistical methods are applied to complex systems in various academic fields. Research in this subcommunity focuses on the practical applications conducted in an interdisciplinary context. The main journals in this subcommunity are PLoS One, Nucleic Acids Research (NAR), and Journal of Agricultural, Biological and Environmental Statistics (JABES), which have high weighted in-degree in this subcommunity.

• •

Medical Science (Brown). This subcommunity contains 28 journals, with a main focus on medical imaging, neuroscience, clinical research, and healthcare management. Research in this subcommunity emphasizes the use of statistical modeling to advance medical technologies and improve patient outcomes. Representative journals in this subcommunity include NeuroImage, Journal of the American Medical Association (JAMA), and Journal of Clinical Oncology (JCO).

• •

Natural Science (Orange). This subcommunity contains 25 journals, which primarily focus on astronomy, ecology, environmental science, and geoscience, as shown in Table 5. These citation patterns highlight the central role of statistical methods in natural science applications. Specifically, Monthly Notices of the Royal Astronomical Society (MNRAS), Methods in Ecology and Evolution (MEE), and Atmospheric Environment have the highest weighted in-degree within this subcommunity.

• •

Biostatistics (Light Green). This subcommunity contains 17 journals, which emphasize biology-related fields such as biochemistry, biotechnology, and genetics. Research in this community focuses on biological discovery through the application of statistical methodologies in the life sciences. The most representative journal in this community is AOAS, as expected for a core journal. In addition, Bioinformatics, Biostatistics, and Nature Genetics also have high weighted in-degree.

Figure 8 reveals the knowledge sources for core journals in the 4 identified subcommunities. The following three main phenomena are observed. First, both AOAS and JASA cited heavily from the “Interdisciplinary Research” and “Medical Science” subcommunities. This reflects their emphasis on interdisciplinary applications of statistics, particularly in domains that combine methodological development with real-world problems. In addition, driven by the growing demand for data analysis in clinical research and health management, the application of statistical methods in medical science has grown substantially [48]. Second, except for AOAS, all other core journals rarely cite from the “Natural Science” subcommunity. This is because AOAS aims to provide a forum for all areas of applied statistics, including nature science. In contrast, other core journals focus primarily on computational applications, methodological development, or the theory of mathematical statistics. Third, almost all core journals cite extensively from the “Biostatistics” subcommunity. This reveals the prominence of biostatistics as one of the most active areas in contemporary applied statistical research [21]. Specifically, JBES primarily targets topics in economics and finance, which explains its limited citation of biostatistics-related work.

6 Discussion

In this study, we construct a weighted bipartite journal citation network using publication data from eight core journals over 2001–2023. To address the challenges posed by degree heterogeneity and network sparsity, we propose a novel community detection algorithm, Bi-SCORE. The method integrates spectral decomposition with a ratio-based normalization to extract latent community structures. Theoretical analysis establishes upper bounds for the numbers of misclustered row and column nodes, respectively. Extensive simulation studies under various network settings, including balanced and imbalanced sample sizes, varying degrees of heterogeneity, and different levels of sparsity, demonstrate that Bi-SCORE consistently outperforms existing methods such as nBiSC and spectral clustering in terms of clustering accuracy. In addition, we apply Bi-SCORE to the real-world citation network and identify 6 distinct communities that correspond to major academic fields. These communities include but are not limited to “Applied Statistics”, “Statistical Methodology”, “Computational Statistics”, and “Mathematical Statistics”. Specifically, we further apply Bi-SCORE within the “Applied Statistics” community and identify 4 subcommunities. The empirical results suggest that Bi-SCORE is well-suited for knowledge source discovery in large-scale citation networks with heterogeneous structure.

Several directions for future research are possible. First, the current analysis focuses on citation data from eight core statistical journals over the period 2001–2023. While these journals represent influential sources within the field of statistics, many relevant works are published in multidisciplinary outlets which also contribute to the diffusion of statistical knowledge. Incorporating a broader range of journals may yield a more comprehensive understanding of cross-disciplinary knowledge flows. Second, the Bi-SCORE method is developed for static bipartite networks. However, journal citation networks evolve over time, and it is of interest to extend Bi-SCORE to dynamic settings to better capture the temporal evolution of community structures. Third, our method assigns each journal to a single community, but in reality, some journals span multiple research areas. Future work could explore mixed-membership models that allow journals to simultaneously belong to multiple communities.

Appendix A Proofs

Appendix A provides supplementary material to support the main text. Sections A.1 - A.10 show the detailed proofs of the propositions, theorem, and lemmas discussed in the main paper.

A.1 Proof of Proposition 1

Proof.

We define the matrices $\mathcal{J}_{\theta}\in\mathbb{R}^{n\times K}$ and $\mathcal{J}_{\gamma}\in\mathbb{R}^{m\times L}$ as follows: for $1\leq i\leq n$, $1\leq j\leq m$, $1\leq k\leq K$ and $1\leq l\leq L$,

[TABLE]

Based on these definitions, along with the definitions of $\Psi_{\theta}$, $\Psi_{\gamma}$, the expectation matrix $\Omega$ is expressed as $\Omega=(\mathcal{J}_{\theta}\|\theta\|\Psi_{\theta})B(\mathcal{J}_{\gamma}\|\gamma\|\Psi_{\gamma})^{\top}$. By introducing the substitution $S\equiv\Psi_{\theta}B\Psi_{\gamma}^{\top}$, we can simplify the expression for $\Omega$ to

[TABLE]

Since the diagonal matrices $\Psi_{\theta}$ and $\Psi_{\gamma}$ are of full rank, the rank of $S$ is given by $\text{rank}(S)=\text{rank}(\Psi_{\theta}B\Psi_{\gamma}^{\top})=\kappa$. We then define the compact SVD of the $K\times L$ matrix $S$ as

[TABLE]

where $\Lambda_{S}$ is a $\kappa\times\kappa$ rectangular diagonal matrix containing the non-negative singular values of $S$ in descending order. The matrices $Y\in\mathbb{R}^{K\times\kappa}$ and $H\in\mathbb{R}^{L\times\kappa}$ are orthogonal matrices.

Substituting the SVD of $S$ from (16) into the expression for $\Omega$ from (15), we obtain:

[TABLE]

From the definitions of $\mathcal{J}_{\theta}$ and $\mathcal{J}_{\gamma}$, it follows that they possess orthonormal columns, specifically $\mathcal{J}_{\theta}^{\top}\mathcal{J}_{\theta}=I_{K\times K}$ and $\mathcal{J}_{\gamma}^{\top}\mathcal{J}_{\gamma}=I_{L\times L}$. Consequently,

[TABLE]

By (A.1), we observe that $\mathcal{J}_{\theta}Y$ and $\mathcal{J}_{\gamma}H$ have orthogonal columns. Thus, (17) constitutes the compact SVD of the matrix $\Omega$. Denoting the compact SVD of $\Omega$ as $\Omega=U\Lambda_{\Omega}V^{\top}$, we can identify

[TABLE]

From (19), we directly derive the relationships given in (10), (11), and (12). As $Y$ is an orthogonal matrix, its rows are orthonormal, meaning $\|Y_{c^{r}_{i}\cdot}\|=1$. Consequently, $\|U_{i\cdot}\|=\|(\theta_{i}/\|\theta^{(c_{i}^{r})}\|)Y_{c^{r}_{i}\cdot}\|=\theta_{i}/\|\theta^{(c_{i}^{r})}\|$. Furthermore, based on (8), we have $\|U_{i\cdot}\|\asymp\theta_{i}/\|\theta\|$. A similar derivation yields $\|V_{j\cdot}\|\asymp\gamma_{j}/\|\gamma\|$. This completes the proof of Proposition 1. ∎

A.2 Proof of Proposition 2

Proposition 2 aims to bound the distance between the singular vectors of $\Omega$ and those of $A$. To achieve this, we first establish several key lemmas. These include Lemma A.1, which demonstrates that the eigenvalues of $\Omega$ are on the order of $\|\theta\|^{2}\|\gamma\|^{2}$; Lemma A.2, which bounds the distance between the random adjacency matrix $A$ and its expected value $\Omega$; and Lemma A.3, which provides a positive lower bound for the singular value gap $\lambda_{1}(S^{\top}S)-\lambda_{2}(S^{\top}S)$. By combining these foundational lemmas, we then apply the Davis-Kahan Theorem (Lemma A.4) to establish Proposition 2.

Lemma A.1**.**

Under the weighted bipartite DCBM and in conjunction with Assumption 2, for $1\leq k\leq\kappa$, the eigenvalues of $\Omega^{\top}\Omega$ (and $\Omega\Omega^{\top}$) are of the order:

[TABLE]

The proof can be found in Appendix A.7.

Lemma A.2**.**

Under Assumption 1, for sufficiently large $n$ and $m$, with probability at least $1-1/n-1/m$, the spectral norm distance between the random adjacency matrix $A$ and its expectation $\Omega$ is bounded as:

[TABLE]

The proof can be found in Appendix A.8.

Lemma A.3**.**

For a matrix $S=Y\Lambda_{S}H^{\top}$, and for $k=1,\dots,K$ and $l=1,\dots,L$, under Assumptions 1 and 2, the following properties hold:

[TABLE]

The proof of Lemma A.3 can be found in Appendix A.9.

Lemma A.4** (Theorem 2 in [68]).**

Let $G,\widehat{G}\in\mathbb{R}^{n\times n}$ be symmetric matrices with eigenvalues $\lambda_{1}\geq\cdots\geq\lambda_{n}$ and $\widehat{\lambda}_{1}\geq\cdots\geq\widehat{\lambda}_{n}$, and corresponding eigenvectors $d_{1},\dots,d_{n}$ and $\widehat{d}_{1},\dots,\widehat{d}_{n}$, respectively. Fix indices $1\leq r\leq s\leq n$. Assume that the eigenvalue gap condition holds: $\min(\lambda_{r-1}-\lambda_{r},\lambda_{s}-\lambda_{s+1})\geq 0$, where we define $\lambda_{0}=\infty$ and $\lambda_{n+1}=-\infty$. Let $k=s-r+1$, $D=(d_{r},d_{r+1},\dots,d_{s})\in\mathbb{R}^{n\times k}$, and define the matrices $\widehat{D}=(\widehat{d}_{r},\widehat{d}_{r+1},\dots,\widehat{d}_{s})\in\mathbb{R}^{n\times k}$. Then, there exists an orthogonal matrix $O\in\mathbb{R}^{k\times k}$ such that

[TABLE]

Proof of Proposition 2.

We begin by bounding the spectral norm of the difference between $A^{\top}A$ and $\Omega^{\top}\Omega$. First, with a probability at least $1-1/n-1/m$, we have

[TABLE]

where the fifth step follows from Lemma A.2, and the sixth step follows from Lemma A.1, which implies $\|\Omega\|=\sqrt{\lambda_{1}(\Omega^{\top}\Omega)}\asymp\|\theta\|\|\gamma\|$.

Next, we apply Lemma A.4 to bound the distance between the top left singular vectors. We have

[TABLE]

where the second step follows from (20), the third step follows from (10) in Proposition 1, and Lemma A.3, which implies that

[TABLE]

and the last step follows from (9), and thus on the right-hand side of the inequality, the first term dominates the second term for large $n$ and $m$.

Similarly, we apply Lemma A.4 to bound the Frobenius norm distance between the subspaces spanned by the 2nd to $\kappa$-th left singular vectors.

[TABLE]

where the second step follows from Proposition 1, which implies that $\lambda_{(\kappa+1)}(\Omega\Omega^{\top})=0$, the third step follows from Lemma A.1, (20) and (21).

Following an analogous procedure for the right singular vectors, we can obtain similar bounds for $\widehat{V}$ and $V$:

[TABLE]

This completes the proof. ∎

A.3 Proof of Proposition 3

Proof.

For two row nodes $i_{1}$ and $i_{2}$, the distance is

[TABLE]

where the second step is derived from Proposition 2, and the fifth equality follows from Proposition 1. Consequently, (A.3) shows that $\|R^{r}_{i_{1}\cdot}-R^{r}_{i_{2}\cdot}\|^{2}=0$ if $c^{r}_{i_{1}}=c^{r}_{i_{2}}$. Otherwise, if $c^{r}_{i_{1}}\neq c^{r}_{i_{2}}$, we obtain

[TABLE]

where the third step holds because $H$ is an orthogonal matrix. The fourth step follows from Lemma A.3.

Similarly, it can be demonstrated that $\|R^{c}_{j_{1}\cdot}-R^{c}_{j_{2}\cdot}\|^{2}\geq 2$ when $c^{c}_{j_{1}}\neq c^{c}_{j_{2}}$, and is zero otherwise. This concludes the proof. ∎

A.4 Proof of Proposition 4

To prove Proposition 4, we first establish Lemma A.5 to bound the number of ill-behaviour nodes, and then present a useful technical inequality in Lemma A.6.

For a constant $0<C<1$, we define the following sets

[TABLE]

After this, Lemma A.5 helps us figure out how many nodes are not in these sets $\widehat{T}^{r}$ and $\widehat{T}^{c}$.

Lemma A.5**.**

For nodes in $\widehat{T}^{r}$ and $\widehat{T}^{c}$, and under Assumptions 1 and 2, the following equations hold

[TABLE]

Furthermore, with a probability at least $1-1/n-1/m$, the cardinality of $\mathcal{V}^{r}\setminus\widehat{T}^{r}$ and $\mathcal{V}^{c}\setminus\widehat{T}^{c}$ satisfy

[TABLE]

The proof of Lemma A.5 can be found in Appendix A.10.

Lemma A.6** (Lemma A.6 in [65]).**

For any two vectors $u,v\in\mathbb{R}^{n}$, and any two positive numbers $a,b\in\mathbb{R}$, the following inequality holds,

[TABLE]

Proof of Proposition 4.

We observe that

[TABLE]

where the last inequality comes from Proposition 1. Next, we want to show $\|R^{r}_{i\cdot}\|^{2}\leq C$. Based on the definition of $R^{r}$, we have

[TABLE]

where the third step follows from Proposition 1, and the fourth step uses Lemma A.3, which implies $Y_{k1}\geq C>0$.

Next, we want to bound the total squared Frobenius norm difference between our estimated representation $\widehat{R}^{r}$ and the true representation $R^{r}$. We split this total difference into two parts: the sum over “ill-behaved” nodes (those outside $\widehat{T}^{r}$) and the sum over “well-behaved” nodes (those inside $\widehat{T}^{r}$).

[TABLE]

For the first part, we have

[TABLE]

where the second step uses equations (3) and (26) to bound the individual terms. The third step states that as $n$ grows large enough, the term $\tau_{n}$ becomes dominant compared to the constant $C$. The last step applies Lemma A.5.

For the second part in (27), we have

[TABLE]

where the first step follows from the fact that $\|R_{i\cdot}^{r}\|\leq C$ (see (26)), and $\tau_{n}$ scales with $n$, which implies $\tau_{n}\geq C\geq\|R_{i\cdot}^{r}\|$ for $n$ large enough. Thus, although (3) shows that the elements in $R_{i\cdot}^{r}$ are truncated by $\tau_{n}$, we still have $\|\widehat{R}^{r}_{i\cdot}-R^{r}_{i\cdot}\|^{2}\leq\left\|\frac{\widehat{U}_{i2\sim\kappa}}{\widehat{U}_{i1}}-\frac{U_{i2\sim\kappa}O_{U}}{C_{U}U_{i1}}\right\|^{2}$ for large $n$, the second step uses Lemma A.6, the third step uses Lemma A.5 and (25), and the last step follows from Proposition 2.

By combining the bounds from (28) and (A.4), we get the total bound for the row representations $\|\widehat{R}^{r}-R^{r}\|_{F}^{2}\leq\frac{C\tau_{n}^{2}\left\{\sqrt{\log(nm)\mathcal{Z}}+\log(nm)\right\}^{2}}{\theta_{\min}^{2}\|\gamma\|^{2}}$. A similar line of reasoning applies to the column representations, leading to $\|\widehat{R}^{c}-R^{c}\|_{F}^{2}\leq\frac{C\tau_{m}^{2}\left\{\sqrt{\log(nm)\mathcal{Z}}+\log(nm)\right\}^{2}}{\gamma_{\min}^{2}\|\theta\|^{2}}$. This completes the proof. ∎

A.5 Proof of Proposition 5

Proof.

Recall that $\widehat{X}^{r}$ is defined as

[TABLE]

Note that $R^{r}$ and $R^{c}$ also belong to $\mathcal{X}_{n,K-1,K}$ and $\mathcal{X}_{m,K-1,L}$, respectively. Therefore,

[TABLE]

From this, we have

[TABLE]

where the last step follows from Proposition 4. Similarly, we can show that

[TABLE]

This completes the proof of Proposition 5. ∎

A.6 Proof of Theorem 1

Proof.

For row nodes $i_{1}$ and $i_{2}$ in set $\mathcal{W}^{r}$, if they belong to different communities, then

[TABLE]

Similarly, for column nodes $j_{1}$ and $j_{2}$ in set $\mathcal{W}^{c}$, if they are in different communities, we have

[TABLE]

By Proposition 3 and the definition of $\mathcal{W}^{r}$ and $\mathcal{W}^{c}$, for $i_{1},~i_{2}\in\mathcal{W}^{r}$ and $j_{1},~j_{2}\in\mathcal{W}^{c}$, we have

[TABLE]

Therefore, if nodes $i_{1}$ and $i_{2}$ belong to different communities, their corresponding rows in $\widehat{X}^{r}$ are sufficiently different. Since we assume that $\lvert\mathcal{V}^{r}\setminus\mathcal{W}^{r}\rvert<\min\{n_{1},\dots,n_{K}\}$, the set $\mathcal{W}^{r}$ contains at least one node from each community. Given that $\widehat{X}^{r}$ has exactly $K$ distinct rows, we conclude that if two nodes in $\mathcal{W}^{r}$ are in the same community, their rows in $\widehat{X}^{r}$ are identical; otherwise, their rows are sufficiently different. Thus, nodes in $\mathcal{W}^{r}$ are clustered correctly. By the same reasoning, nodes in $\mathcal{W}^{c}$ are also clustered correctly. From the definitions of $\mathcal{W}^{r}$, $\mathcal{W}^{c}$ and Proposition 5, we have

[TABLE]

Similarly, we obtain

[TABLE]

This completes the proof. ∎

A.7 Proof of Lemma A.1

Proof of Lemma A.1.

According to Proposition 1, we have $\sigma_{k}(\Omega)=\|\theta\|\|\gamma\|\sigma_{k}(S)$, for $1\leq k\leq\kappa$. Therefore, it is sufficient to show $0<C_{1}\leq\sigma_{\kappa}(S)\leq\sigma_{1}(S)\leq C_{2}$.

Recall $S=\Psi_{\theta}B\Psi_{\gamma}^{\top}$, where $\Psi_{\theta}$ and $\Psi_{\gamma}$ are diagonal matrices. The terms $\|\Psi_{\theta}\|$ and $\|\Psi_{\theta}\|_{\min}$ represent the largest and smallest absolute values of the diagonal entries of $\Psi_{\theta}$, respectively. From (8), we have

[TABLE]

Thus, there exist two constants $C_{m}>0$ and $C_{M}>0$ such that

[TABLE]

A similar argument shows that

[TABLE]

From the definition of $S$, we have

[TABLE]

because $B$ is a constant matrix (i.e., it does not change with $n$), and therefore there exists a constant $C$, such that $\|B\|\leq C$. On the other hand,

[TABLE]

where the last inequality follows from Assumption 1.

Combining (31) and (32), we get $\sigma_{k}\asymp C$ for $k=1,\dots,\kappa$. Therefore, we have

[TABLE]

This completes the proof of Lemma A.1. ∎

A.8 Proof of Lemma A.2

We first introduce the following lemma.

Lemma A.7** (Corollary 4 in [6]).**

Let $A$ be a $n\times m$ random matrix whose entries $A_{ij}$ follow a Poisson distribution with mean $\Omega_{ij}$. Then, for any $x>0$, we have

[TABLE]

where $\Omega=(\Omega_{ij})$, $\|\Omega\|_{1,\infty}=\max_{i}\|\Omega_{i\cdot}\|_{1}$, and $\|\Omega\|_{\infty,1}=\max_{j}\|\Omega_{\cdot j}\|_{1}$.

Proof of Lemma A.2.

From Assumption 1, we have $\Omega_{ij}=\theta_{i}\gamma_{j}B_{c_{i}^{r}c_{j}^{c}}\leq\theta_{i}\gamma_{j}$. Therefore,

[TABLE]

Let $\mathcal{Z}=\max(\theta_{\max},\gamma_{\max})\max(\|\theta\|_{1},\|\gamma\|_{1})$, where $\theta_{\max}=\max_{i}\theta_{i}$ and $\gamma_{\max}=\max_{j}\gamma_{j}$. Then, $\max(\|\Omega\|_{1,\infty},\|\Omega\|_{\infty,1})\leq\mathcal{Z}$. Applying Lemma A.7 with $x=\log(nm)$, we obtain that, with probability at least $1-1/n-1/m$,

[TABLE]

This completes the proof. ∎

A.9 Proof of Lemma A.3

We first present the following lemma.

Lemma A.8** (Theorem 8.4.4 in [32], Lemma A.7 in [65]).**

For every $K\times K$ irreducible nonnegative, and positive semidefinite matrix $M$, let $U_{\cdot 1}$ be the eigenvector corresponding to its largest eigenvalue. Then,

• (i)

$U_{\cdot 1}$* can be a positive vector.*

• (ii)

The largest eigenvalue is an algebraically simple eigenvalue.

Proof of Lemma A.3.

To prove Lemma A.3, it is sufficient to show that $S^{\top}S$ and $SS^{\top}$ are irreducible and nonnegative, and these properties do not depend on $n$ or $m$. Once this is established, part (i) of Lemma A.8 implies that $Y_{k1}\geq C>0$ and $H_{l1}\geq C>0$, where $Y_{\cdot 1}$ and $H_{\cdot 1}$ are the eigenvectors corresponding to $SS^{\top}$ and $S^{\top}S$, respectively. Moreover, $Y_{k1}\geq C>0$ leads to $U_{i1}>0$, since $U_{i\cdot}=(\theta_{i}/\|\theta^{(c_{i}^{r})}\|)Y_{c_{i}^{r}\cdot}$ for $i=1,\dots,n$. Similarly, we have $V_{j1}>0$ for $j=1,\dots,m$. Part (ii) of Lemma A.8 further gives $\lambda_{1}(S^{\top}S)-\lambda_{2}(S^{\top}S)=\lambda_{1}(SS^{\top})-\lambda_{2}(SS^{\top})\geq C$. This completes the proof of Lemma A.3.

Next, we verify that $S^{\top}S$ and $SS^{\top}$ are irreducible and nonnegative. Recall that $S=\Psi_{\theta}B\Psi_{\gamma}^{\top}$. From (29), (30), and the definition of the diagonal matrices $\Psi_{\theta}$ and $\Psi_{\gamma}$, we see that there exist constants $C_{1}$ and $C_{2}$ such that $0<C_{1}\leq\Psi_{\theta,kk}\leq C_{2}$ and $0<C_{1}\leq\Psi_{\gamma,ll}\leq C_{2}$. Therefore, for $k=1,\dots,K$ and $l=1,\dots,L$,

[TABLE]

For $1\leq l_{1},l_{2}\leq L$, we have

[TABLE]

Similarly, for $1\leq l_{1},l_{2}\leq L$, we obtain

[TABLE]

Since $B$ is a constant matrix with positive entries, $B^{\top}B$ is irreducible and nonnegative by Assumption 1. Hence, $S^{\top}S$ is also nonnegative and irreducible, and these properties do not depend on $n$ or $m$. By the same reasoning, $SS^{\top}$ is nonnegative and irreducible, and these properties do not depend on $n$ or $m$. This completes the proof.

∎

A.10 Proof of Lemma A.5

Proof of Lemma A.5.

We first prove $|C_{U}U_{i1}|\asymp|\theta_{i}|/\|\theta\|$ for $i=1,\dots,n$. By Proposition 1, we have $U_{i\cdot}=\frac{\theta_{i}}{\|\theta^{(c^{r}_{i})}\|}Y_{c^{r}_{i}\cdot}$, and thus $C_{U}U_{i1}=C_{U}\frac{\theta_{i}}{\|\theta^{(c^{r}_{i})}\|}Y_{c^{r}_{i}1}$, where $|C_{U}|=1$ by Proposition 2. From (8) and Lemma A.3, we obtain

[TABLE]

Next, by (24), nodes in the set $\widehat{T}^{r}$ satisfies $\left|\frac{\widehat{U}_{i1}}{C_{U}U_{i1}}-1\right|\leq C<1$, and therefore $|\widehat{U}_{i1}|\asymp|C_{U}U_{i1}|$. Thus, for $i\in\widehat{T}^{r}$ we have

[TABLE]

Similarly, $|\widehat{V}_{j1}|\asymp|C_{V}V_{j1}|\asymp\frac{|\gamma_{j}|}{\|\gamma\|}$ for $j\in\widehat{T}^{c}$.

Furthermore, by (33) and (5) in Assumption 2, we have $|C_{U}U_{i1}|\asymp\frac{|\theta_{i}|}{\|\theta\|}>0$, so using $C_{U}U_{i1}$ as the denominator for all $i$ is valid. We then derive

[TABLE]

where the second step follows from (33), and the last step follows from Proposition 2.

Since nodes in the set $\mathcal{V}^{r}\setminus\widehat{T}^{r}$ satisfy $\{\widehat{U}_{i1}/(C_{U}U_{i1})-1\}^{2}\geq C_{0}^{2}$, we have

[TABLE]

Similarly, we can show that

[TABLE]

This completes the proof. ∎

Appendix B Additional Results

Appendix B presents additional results from our analysis of the journal citation network.