Scalable Exact Densest P-Partite Subgraph Search in Heterogeneous Information Networks

Scalable Exact Densest 𝑃-Partite Subgraph Search in Heterogeneous Information Networks Jiaming Yang∗

Jiadong Xie∗

The Chinese University of Hong Kong jdxie@se.cuhk.edu.hk

Beijing Institute of Technology yumoym456@gmail.com

Kangfei Zhao

Jeffrey Xu Yu

Beijing Institute of Technology zkf1105@gmail.com

HKUST (Guangzhou) jeffreyxuyu@hkust-gz.edu.cn same predecessors

arXiv:2608.03061v1 [cs.DB] 4 Aug 2026

ABSTRACT Heterogeneous information networks (HINs) model typed entities and typed relations, where dense cross-type structures can reveal cohesive semantic patterns such as prolific author-paper-venue groups. Given a query meta-path, the densest 𝑃-partite subgraph search (DPpS) problem jointly selects a nonempty vertex set at each typed position and maximizes the number of induced metapath instances normalized by the geometric mean of the selected set sizes. Existing exact methods solve DPpS by searching over iRM-sets and reducing each fixed-M problem to minimum-cut computations. However, their scalability is limited by the large number of candidate iRM-sets and the high cost of repeatedly solving large auxiliary networks. In this paper, we propose BoxDPpS, an efficient exact approach that reduces both sources of cost. It performs box-level search with safe region pruning, eliminates redundant representations of the same iRM-set, improves early pruning through bounded warm-up, and compresses each fixed-M auxiliary network for exact parametric pseudoflow solving. Experiments on seven real-world datasets show that BoxDPpS preserves the exact DPpS optimum while achieving an average speedup of 27.04× over the state-of-the-art method.

1

INTRODUCTION

Heterogeneous information networks (HINs) model data whose entities and relations have distinct types. For example, a bibliographic network connects authors, papers, and venues, whereas a recommendation network relates users, movies, genres, and tags. In such networks, meaningful structures often arise from the joint participation of several entity types rather than from the connectivity of a single type. A group of authors, for instance, is better characterized together with the papers they write and the venues in which those papers are published. Meta-paths provide a standard mechanism for specifying such cross-type semantics [13, 27, 29]. A meta-path describes a sequence of vertex and relation types that a valid interaction must follow. For example, the meta-path Author → Paper → Venue in Fig. 1(a) represents the semantic relation that an author writes a paper published in a venue, and each concrete author–paper–venue sequence satisfying these relations is an instance of the meta-path. A large body of existing HIN research studies typed structures from either an object-level or a query-driven perspective. Metapath similarity measures quantify the semantic relevance between ∗ Jiadong Xie and Jiaming Yang are the joint first authors.

𝑣2

𝑣3

𝑝3

𝑝4

𝑝5

𝑎1

𝑎2

𝑣1 A

write

P

pubIn

V

𝑝1

𝑝2

P :𝐴 →𝑃 →𝑉 Author

Paper

Venue

𝑎3

𝑎4

14 selected instances

(a) Schema and meta-path

(b) Running HIN instance

Figure 1: An Example of Author–Paper–Venue HIN individual objects [27, 29], whereas HIN community-search models retrieve subgraphs satisfying query, degree, core, truss, motif, relation, or attribute constraints [7, 14, 16, 18, 35]. These formulations are effective when an analyst has designated query vertices, a target vertex type, or a desired cohesiveness threshold. However, they are not designed to globally identify a collection of typed vertex sets that are jointly dense with respect to complete cross-type interactions. The densest 𝑃-partite subgraph search (DPpS) problem [4] offers a global-density perspective. Given a query meta-path P, DPpS jointly selects a nonempty vertex set for every typed position along P and maximizes the number of complete meta-path instances induced by these sets, normalized by the geometric mean of their cardinalities. Fig. 1(b) illustrates the densest vertex sets selected for an Author–Paper–Venue query. By counting complete meta-path instances rather than individual edges or pairwise similarities, DPpS directly measures how intensively all participating types interact. Its geometric-mean normalization treats the typed positions symmetrically and discourages enlarging a selected set unless the additional vertices contribute sufficiently many complete instances. Thus, DPpS requires neither a query vertex nor a designated target type or cohesiveness threshold, making it suitable for exploratory discovery of compact cross-type groups. Such groups may correspond, for example, to authors publishing intensively through a collection of papers and venues, as shown in Fig. 4. It is a huge challenge to solve the DPpS problem exactly. The DPpS objective is a ratio: its numerator depends jointly on all selected vertex sets, while its denominator is a nonlinear geometric mean whose optimal type-wise cardinalities are unknown in advance. A principled route to exactness is to combine the iRM-set formulation with a minimum-cut reduction [4]. An iRM-set M represents the relative cardinalities of the vertex sets selected at the different positions of P. Once M is fixed, the geometric-mean

Jiadong Xie, Jiaming Yang, Kangfei Zhao, and Jeffrey Xu Yu

denominator can be expressed as a linear weighted sum of the selected type-wise cardinalities. For a density threshold 𝛾, determining whether a denser solution exists then becomes a problem that can be solved exactly through a minimum 𝑠-𝑡 cut. The state-ofthe-art exact algorithm, AdvExactGVIt [4], builds on this decomposition by searching over admissible iRM-sets and invoking iterative minimum-cut computations for the corresponding fixed-M subproblems. However, such methods still have two issues. First, since the optimal cardinalities are unknown, an exact algorithm must reason over a potentially large collection of admissible iRM-sets. The underlying count domain grows as the product of the number of available vertices at the different typed positions of the query meta-path, and even an iRM-set that is eventually pruned can incur nontrivial generation and condition-evaluation costs. Second, every surviving iRM-set requires conducting minimum-cuts, which creates one auxiliary node for every meta-path instance and connects it to its participating vertices. Because the number of meta-path instances can be much larger than the number of distinct vertices, repeatedly constructing and solving these networks is expensive. To address the efficiency bottlenecks of AdvExactGVIt, we propose BoxDPpS, an efficient exact approach that reduces both the number of iRM-sets that need to be processed and the cost of solving each surviving fixed-M problem. The main contributions are summarized as follows. ➊ We develop an exact box-level search to avoid exhaustively enumerating and processing the potentially large number of iRM-sets. Even if an iRM-set is eventually pruned and never reaches the fixed-M solver, generating it and evaluating the pruning conditions still incur substantial overhead due to point-by-point enumeration. Thus, BoxDPpS organizes count vectors into integer boxes and prunes an entire box whenever a safe upper bound proves that none of the solutions represented by the box can improve the current best density. ➋ To reduce the number of iRM-sets processed, we eliminate redundant representations of the same iRM-set. We prove a one-to-one correspondence between the resulting canonical representations and admissible iRM-sets, ensuring that each fixed-M problem is processed at most once. ➌ To further increase the number of iRM-sets that can be pruned, we design a bounded warm-up strategy. It evaluates a small set of representative ratios to obtain feasible solutions and schedules the most promising ratios for early exact processing. ➍ To accelerate each fixed-M problem, we compress its auxiliary network and prove that the compression preserves the original objective. We then solve the resulting compact network exactly using parametric Hochbaum pseudoflow [11], reducing the per-iRM-set complexity 2 log 𝑛 ) in BoxDPpS, from 𝑂 (|𝑉 (𝐷)| 3 ) in AdvExactGVIt to 𝑂 (𝑛𝑀 𝑀 where 𝑛𝑀 ≤ |𝑉 (𝐷)|. ➎ We conduct extensive experiments on seven real-world HIN datasets. On the query meta-paths completed by both exact methods, BoxDPpS achieves a speedup of 27.04× over AdvExactGVIt on average. It also completes several large workloads on which AdvExactGVIt times out or exhausts memory, while preserving the exact DPpS optimum. Roadmap. The rest of this paper is organized as follows. We introduce the preliminaries in § 2. § 3 introduces the existing algorithm and its issues. We present an overview of our algorithm BoxDPpS in § 4. § 5, § 6, § 7, and § 8 elaborate on the deduplication of iRM-set, box-level search, bounded warm-up, and acceleration of fixed-M

Table 1: Summary of Notations Notation

G 𝑉,𝐸 A, R 𝜙,𝜓 𝑇G P = (𝐴1 , . . . , 𝐴𝑘 ) 𝑈𝑖 , 𝑁𝑖 FP V V 𝑑𝑖 (𝑣; I ) c( V ) c̄ D𝑐 𝐺 (c) M(c) M 𝜌 ( V ), 𝜌∗ ∗ 𝜌M 𝜌 best 𝛾

Definition input heterogeneous information network vertex set and edge set of G vertex-type set and edge-type set vertex-type and edge-type mappings network schema of G query meta-path with 𝑘 typed positions vertices admissible at position 𝑖 and 𝑁𝑖 = |𝑈𝑖 | all instances of P in G 𝑃 -family (𝑉1 , . . . , 𝑉𝑘 ) set of all feasible 𝑃 -families support of a position-𝑖 vertex in an instance collection I count vector ( |𝑉1 |, . . . , |𝑉𝑘 | ) primitive count vector obtained by gcd normalization count domain [𝑁 1 ] × · · · × [𝑁𝑘 ] Î geometric mean ( 𝑖 𝑐𝑖 ) 1/𝑘 iRM-set induced by count vector c set of all admissible iRM-sets density of V and the global optimum density optimum density restricted to families conforming to M best feasible density found so far by an algorithm density threshold used in a fixed-M auxiliary problem

solving, respectively. We report our experimental studies in § 9 and review the related works in § 10. § 11 concludes the paper.

2

PRELIMINARIES

This section introduces the notation used throughout the paper and the densest 𝑃-partite subgraph search problem. Table 1 summarizes the frequently used notations. A heterogeneous information network (HIN) is a graph G = (𝑉 , 𝐸, A, R, 𝜙,𝜓 ), where 𝑉 and 𝐸 are the vertex and edge sets, A and R are the vertex-type and edge-type sets, 𝜙 : 𝑉 → A maps each vertex to a type, and 𝜓 : 𝐸 → R maps each edge to a relation type. The network is heterogeneous if |A| > 1 or |R| > 1. The network schema 𝑇 G is the type-level graph induced by A and R. Each schema vertex represents a vertex type, and each schema edge represents an admissible relation type between two vertex types. A query over an HIN is usually specified on this schema rather than on individual vertices. Meta-path and Instance: A query meta-path P is a typed path 𝑅1

𝑅2

on 𝑇 G : P = 𝐴1 −−→ 𝐴2 −−→ · · · −−−−→ 𝐴𝑘 , where 𝐴𝑖 ∈ A and 𝑅𝑖 ∈ R. For position 𝑖 ∈ [𝑘] = {1, . . . , 𝑘 }, let 𝑈𝑖 = {𝑣 ∈ 𝑉 | 𝜙 (𝑣) = 𝐴𝑖 } and 𝑁𝑖 = |𝑈𝑖 |. When relation types are clear from context, we write P = (𝐴1, . . . , 𝐴𝑘 ) and call 𝑘 the number of vertex-type positions in P. A vertex sequence 𝑝 = (𝑣 1, . . . , 𝑣𝑘 ) ∈ 𝑈 1 × · · · ×𝑈𝑘 is an instance of P if every consecutive pair follows the required relation, i.e., (𝑣𝑖 , 𝑣𝑖+1 ) ∈ 𝐸 and 𝜓 ((𝑣𝑖 , 𝑣𝑖+1 )) = 𝑅𝑖 for all 𝑖 < 𝑘. We write 𝑝𝑖 = 𝑣𝑖 and denote by FP the set of all such instances. For example, consider the meta-path P = (Author, Paper, Venue) in Fig. 1. The vertex sequences (𝑎 1, 𝑝 1, 𝑣 1 ) and (𝑎 2, 𝑝 4, 𝑣 2 ) are instances of P. Following the setting in existing work [4], we assume that the query positions have pairwise distinct vertex types. Hence, the sets 𝑈 1, . . . , 𝑈𝑘 are disjoint, and every query vertex has an unambiguous position. For an instance collection I ⊆ FP , the support of a position-𝑖 vertex 𝑣 ∈ 𝑈𝑖 is 𝑑𝑖 (𝑣; I) = {𝑝 ∈ I | 𝑝𝑖 = 𝑣 } . 𝑅𝑘 −1

𝑃-family and Induced Instances: Given P, a 𝑃-family is a tuple V = (𝑉1, . . . , 𝑉𝑘 ) such that ∅ ≠ 𝑉𝑖 ⊆ 𝑈𝑖 for every 𝑖 ∈ [𝑘]. Let V denote the set of all 𝑃-families. The instances induced by V are F (V) = {𝑝 ∈ FP | 𝑝𝑖 ∈ 𝑉𝑖 for all 𝑖 ∈ [𝑘]}. The induced 𝑃-partite subgraph consists of the selected vertices and the query-compatible

Scalable Exact Densest 𝑃 -Partite Subgraph Search in Heterogeneous Information Networks

edges between consecutive positions. Since every objective, pruning rule, and flow construction in this paper depends on the complete instances in F (V), we use the family and its induced instance set as the primary representation. For V ∈ V, its count vector is c(V) = (𝑐 1, . . . , 𝑐𝑘 ), where 𝑐𝑖 = Î  1/𝑘 𝑘 |𝑉𝑖 |. Let 𝐺 (c) = be the geometric mean of the position𝑖=1 𝑐 𝑖 wise counts. The admissible count domain is D𝑐 = [𝑁 1 ] ×· · ·× [𝑁𝑘 ], where [𝑁𝑖 ] = {1, . . . , 𝑁𝑖 }. Density: The density of a 𝑃-family V is defined as 𝜌 (V) =

|F (V)| |F (V)| = Î𝑘 . 𝐺 (c(V)) ( 𝑖=1 |𝑉𝑖 |) 1/𝑘

The numerator counts complete query instances rather than individual edges, while the geometric-mean denominator normalizes all query positions symmetrically. For example, the highlighted 𝑃family in Fig. 1 has count vector (4, 2, 2) and induces 14 meta-path 14 14 instances. Its density is (4×2×2) . 3 1/3 = √ 16 Based on the above definitions, we can formulate the Densest 𝑃-Partite Subgraph Search (DPpS) problem as follows. Densest 𝑃-Partite Subgraph Search Problem: Given an HIN G and a query meta-path P = (𝐴1, . . . , 𝐴𝑘 ), the DPpS problem aims to find a 𝑃-family V ∗ ∈ arg max V ∈V 𝜌 (V). The subgraph induced by V ∗ is called a densest 𝑃-partite subgraph, and we denote its density as 𝜌 ∗ = max V ∈V 𝜌 (V). To explain how exact search checks whether a target density is attainable, we next introduce the iRM-set formulation from previous work [4]. An iRM-set encodes the relative position-wise cardinalities of a 𝑃-family. Once an iRM-set is fixed, the geometric-mean denominator can be rewritten as a linear weighted sum. iRM-set: For a positive count vector c = (𝑐 1, . . . , 𝑐𝑘 ) ∈ D𝑐 , its induced iRM-set is M(c) = (𝑚 1, . . . , 𝑚𝑘 ), where 𝑚𝑖 = 𝐺𝑐(c) . The 𝑖 entries of iRM-set are indexed by the positions of P, we therefore write it as an ordered vector. A 𝑃-family V conforms to M if M(c(V)) = M. Let M = {M(c) | c ∈ D𝑐 } be the set of all admissible iRM-sets. For example, the highlighted√family in Fig. 1 has c = (4, 2, 2). Hence, 𝐺 (c)√= (4√× 2 × 2) 1/3 = 3 16, and its induced √ 3 16 3 16 3 16 iRM-set is M(c) = 4 , 2 , 2 .

𝑚𝑖 For every 𝑖, 𝑗 ∈ [𝑘], the definition gives 𝑚𝑖 𝑐𝑖 = 𝐺 (c), 𝑚 = 𝑐𝑖𝑗 , 𝑗 Î𝑘 Í 𝑘 1 𝑖=1 𝑚𝑖 = 1, and 𝐺 (c) = 𝑘 𝑖=1 𝑚𝑖 𝑐 𝑖 . The last identity converts the geometric-mean denominator into a linear weighted sum once ∗ = M is fixed. For M ∈ M, define the ratio-restricted optimum 𝜌 M ∗ ∗ max{𝜌 (V) | V ∈ V, M(c(V)) = M}. Thus, 𝜌 = maxM∈M 𝜌 M . 𝑐

3

EXISTING APPROACH AND ISSUES

The DPpS objective is difficult to optimize exactly because its denominator is a nonlinear geometric mean and the optimal cardinalities of the typed vertex sets are unknown. Existing exact search [4], named AdvExactGVIt, addresses this difficulty by combining iRMsets with minimum-cut computations. At a high level, it first fixes a relative-cardinality pattern, checks whether a given density can be achieved under this pattern, and repeatedly increases the density until it can no longer be improved. Since the optimal relativecardinality pattern is unknown, this process must be performed over multiple iRM-sets.

From an iRM-Set to a Density Check: An iRM-set represents the relative cardinalities of the selected vertex sets. Consider the highlighted Author–Paper–Venue family in Fig. 1, whose count vector is (4, 2, 2). This vector means that the selected numbers of authors, papers, and venues follow the ratio 2 : 1 : 1. Its iRM-set assigns weights inversely proportional to these cardinalities. Therefore, papers and venues receive twice the cardinality weight of authors. The benefit of fixing an iRM-set is that it converts the nonlinear denominator of the DPpS density into a linear weighted sum. Specifically, for every family conforming to M(c) = (𝑚 1, . . . , 𝑚𝑘 ), Í its geometric-mean denominator satisfies 𝐺 (c) = 𝑘1 𝑘𝑖=1 𝑚𝑖 |𝑉𝑖 |. Therefore, checking whether a density threshold 𝛾 is attainable becomes a linear gain–cost problem: the selected meta-path instances provide the gain, whereas the selected vertices contribute a weighted cost proportional to 𝛾. The resulting problem can be solved exactly by a minimum 𝑠-𝑡 cut. For example, the highlighted family in Fig. 1 has a count vector (4, 2, 2), so the paper and venue positions have twice the weight√of the author position. It induces 14 instances and has density 14/ 3 16 ≈ 5.56. Hence, it demonstrates that the threshold 𝛾 = 5 is attainable. The minimum-cut computation performs the same test over all possible vertex selections under the fixed iRM-set. Iterative Density Checking: For one fixed iRM-set, the exact solver starts from a feasible density threshold 𝛾 and asks whether a denser family exists. If the minimum cut returns a family with density larger than 𝛾, the threshold is raised to the density of that family and checked again. Otherwise, no family associated with this iRM-set can improve the current threshold, and the fixed iRM-set solve terminates. Search over iRM-Sets: A fixed iRM-set covers one relative-cardinality pattern. Since the optimal pattern is unknown, AdvExactGVIt generates the iRM-sets induced by feasible count vectors and processes the candidates that cannot be pruned. For each such candidate, it may first obtain a feasible lower bound and reduce the graph, and then run the iterative density checks described above. A completed exact solve for an iRM-set M𝑎 = (𝑚𝑎1 , . . . , 𝑚𝑘𝑎 ) may improve the global best density and produce a point-wise certificate for pruning other iRM-sets. In particular, a remaining iRM-set M𝑐 = Í 𝑚𝑎 (𝑚𝑐1, . . . , 𝑚𝑐𝑘 ) can be removed when 𝑘 ≤ 𝑘𝑖=1 𝑚𝑖𝑐 ≤ 𝐵𝑎 . Here, 𝐵𝑎 is 𝑖 computed from the output of the completed fixed-M𝑎 solve rather 𝑎 than from M alone. If the solve returns a family whose induced Í iRM-set is M′ = (𝑚 1′ , . . . , 𝑚𝑘′ ), the certificate uses 𝐵𝑎 = 𝑘𝑖=1 𝑚𝑖𝑎 /𝑚𝑖′ . Alternatively, if the solve derives a certified upper bound 𝜌b𝑎 on ∗ and obtains a feasible family V ′ satisfying 𝜌 (V ′ ) > 𝜌 𝜌M b𝑎 , the 𝑎 certificate uses 𝐵𝑎 = 𝑘𝜌 (V ′ )/b 𝜌𝑎 [4]. Overall, the existing exact algorithm consists of an outer search over iRM-sets and an inner sequence of density checks for every surviving iRM-set. The overall cost contains three parts: 𝑇AdvExactGVIt = 𝑇generation + 𝑇certificate + 𝑁 solve𝑇flow, where 𝑇generation is the cost of enumerating count vectors and constructing their iRM-set representations, 𝑇certificate is the cost of checking the available crossiRM-set pruning conditions, and 𝑁 solve is the number of candidates that reach exact fixed-M solving. The dominant exact-solving term is commonly summarized as Θ(|M|𝑇flow ) in the existing analysis, where |M| is the number of iRM-sets processed and 𝑇flow is the

Jiadong Xie, Jiaming Yang, Kangfei Zhao, and Jeffrey Xu Yu

cost of solving one fixed-M problem. The existing analysis gives 𝑇flow = 𝑂 (|𝑉 (𝐷)| 3 ) for the flow network 𝐷, resulting in a worst-case total complexity of 𝑂 ((𝑛/𝑘) 4𝑘 ) [4].

Issues: The time complexity Θ(|M|𝑇flow ) exposes two multiplicative sources of cost: the number of candidate iRM-sets |M| that the search must process, and the cost 𝑇flow of solving a single fixed-M subproblem. Both factors grow rapidly with the graph size and the meta-path length 𝑘. We thus analyze the two dimensions separately. (1) iRM-set Search Issue. The number of candidate iRM-sets can be very large, leading to the following three issues. First, AdvExactGVIt enumerates the entire count domain [1, 𝑁 1 ]× · · · × [1, 𝑁𝑘 ] point by point, maps each vector to its iRM-set, and only then evaluates the pruning conditions. As a consequence, an iRM-set must be generated before it can be discarded: even candidates that are ultimately pruned are individually generated and Î tested, so the search pays at least Ω( 𝑖 𝑁𝑖 ) work irrespective of how few candidates survive. The pruning is itself point-wise; a certificate obtained from one completed solve is re-tested against every remaining candidate, so in the worst case, certificate handling grows to Θ(|M| 2 ). Second, since an iRM-set encodes only relative sizes, a count vector c and every integer multiple 𝑔c induce the same M. For example, the count domain of Fig. 1 contains c = (2, 1, 1) and c′ = (4, 2, 2) = 2c. Doubling every coordinate also doubles the geometric mean, i.e., 𝐺 (c′ ) = 2𝐺 (c). Hence, the common factor cancels in every ratio 𝑚𝑖 = 𝐺 (c)/𝑐𝑖 , and M(c′ ) = M(c), and pointwise countvector enumeration in AdvExactGVIt can repeatedly encounter the same iRM ratio. Even if duplicate iRM objects are removed by a set data structure, the enumeration and lookup overhead remain high due to the large number of candidate iRM-sets. Third, effective pruning may become available only after sufficiently good feasible solutions or useful exact certificates have been obtained. However, AdvExactGVIt initializes 𝜌 best = 0 and improves it only as iRM-sets are processed in the order induced by Cartesian enumeration. If promising ratios appear late in this order, the incumbent remains weak and few useful certificates are available during the early stages of the search, when the candidate space is still large. Consequently, a substantial portion of M may be generated and examined before effective pruning becomes possible. (2) Fixed-M Solving Issue. For each surviving iRM-set, the fixed-M solver builds a flow network with one auxiliary node per metapath instance and 𝑘 infinite-capacity arcs from that node to its endpoint vertices. In dense regions, the number of instances F0 is combinatorially larger than the number of distinct vertices—a block fully connected across positions induces up to (𝑛/𝑘)𝑘 instances over only 𝑛 vertices—so the network is dominated by |F0 |, and its per-solve cost 𝑇flow = 𝑂 (|𝑉 (𝐷)| 3 ) is cubic in this inflated size. These issues motivate the design of our approach BoxDPpS in the next section. BoxDPpS keeps the same DPpS objective and the same fixed-M auxiliary formulation, but changes how count ratios are generated and how to represent and solve fixed-M closure.

4

OVERVIEW OF BoxDPpS

§ 3 shows that the cost of exact DPpS search is determined by two multiplicative factors: the number of iRM-sets processed and the cost of the repeated density checks for each surviving iRM-set. The

iRM-set search side further suffers from point-wise enumeration, duplicate ratio representations, and weak pruning at the beginning of the search. BoxDPpS addresses these issues in the same order while preserving the original iRM-set formulation and exactness guarantee. Avoiding Point-Wise Count-Vector Enumeration: To avoid the cost incurred by AdvExactGVIt in enumerating every count vector individually, BoxDPpS organizes the count domain into integer Î boxes and searches these boxes directly. A box X = 𝑘𝑖=1 [ℓ𝑖 , ℎ𝑖 ] represents all count vectors whose 𝑖-th coordinate lies between ℓ𝑖 and ℎ𝑖 . For example, the box X = [1, 4] × [1, 2] × [1, 2] represents 16 Author–Paper–Venue count vectors. A point-wise search must generate these vectors before deciding whether their iRM-sets can be discarded. In contrast, for each box, BoxDPpS computes a safe upper bound 𝑈 box (X) on the maximum density attainable by any 𝑃-family whose count vector belongs to X. If 𝑈 box (X) ≤ 𝜌 best , the entire box is pruned without materializing its individual count vectors or their induced iRM-sets. Thus, a large, unpromising region can be eliminated with a single box-level test rather than by pointwise generation and evaluation. Completed exact solves may additionally produce point-wise certificates for pruning other iRM-sets. BoxDPpS lifts such a certificate to the box level whenever it covers every count vector represented by the box, allowing all corresponding candidates to be removed together. A box that cannot be pruned is divided into smaller disjoint boxes, and only a surviving singleton is converted into its induced iRM-set and sent to fixed-M processing. Section 6 presents the safe upper bounds, certificate lifting, box-splitting procedure, and correctness proof. Eliminating Scale-Equivalent iRM-Sets: The second issue is that proportional count vectors induce the same iRM-set. For example, (2, 1, 1) and (4, 2, 2) describe the same Author–Paper–Venue cardinality ratio and therefore require the same fixed-M solve. To avoid processing such duplicates, every surviving count vector is mapped to its primitive representative: Primitive(c) = c/gcd(𝑐 1, . . . , 𝑐𝑘 ). All vectors on the same scaling ray therefore share one canonical key. For example, (4, 2, 2) is mapped to (2, 1, 1), and the corresponding iRM-set is solved only once. Section 5 proves that primitive count vectors are in one-to-one correspondence with admissible iRM-sets, so canonicalization removes only redundant representations and preserves every feasible relative-cardinality pattern. Strengthening Early iRM-Set Pruning: The third issue identified in Section 3 is that effective pruning may become available late because both the initial incumbent and the set of exact certificates are weak. This issue is particularly important for the box search of BoxDPpS, because a safe box upper bound can prune a box only when it does not exceed the current best feasible density 𝜌 best . A stronger incumbent therefore allows more unpromising boxes to be removed early. To obtain a useful lower bound without exhaustively solving the search space, BoxDPpS performs a bounded warmup before the main box traversal. It evaluates a small collection of representative cardinality ratios and uses the returned feasible families to improve 𝜌 best . The ratios producing the best feasible solutions are then processed exactly as early seeds. Besides further improving 𝜌 best , these completed exact solves may generate iRMset certificates that become available to the subsequent box search. Approximate results affect only the feasible lower bound and the

Scalable Exact Densest 𝑃 -Partite Subgraph Search in Heterogeneous Information Networks

Algorithm 1: BoxDPpS (G, P)

𝜔( 𝑣ˆ𝑇 ) = 2

Figure 2: Reducing the Network for A Fixed-M Density Check processing order; they are never used as pruning certificates. § 7 presents the candidate-selection procedure and its safety argument. Accelerating Fixed-M Density Checks: The final issue is that every surviving iRM-set requires repeated density checks on a flow network containing one auxiliary node per meta-path instance. BoxDPpS retains the same fixed-M gain–cost objective but removes repeated structures from its network representation. Fig. 2 illustrates the idea using the subgraph induced by authors {𝑎 1, 𝑎 2 }, papers {𝑝 4, 𝑝 5 }, and venues {𝑣 2, 𝑣 3 } in Fig. 1. This subgraph contains six Author–Paper–Venue instances, each represented separately in the existing network. However, 𝑣 2 and 𝑣 3 have identical predecessor multisets at the terminal position in the current instance universe and can be represented by one weighted vertex. Moreover, several instances share the same projections on the selected groups of meta-path positions. They can therefore be aggregated into weighted interactions. In this example, six individual instance representations are reduced to three grouped interactions. We prove that these transformations preserve the fixed-M gain–cost objective for every density threshold. The compact network therefore returns exactly the same answer as the original network, while its size depends on the number of distinct structural signatures rather than directly on the total number of meta-path instances. Recall that a fixed-M solve repeatedly increases 𝛾 and checks whether the new density can be exceeded. Across these checks, the structure of the compact network remains unchanged; only the vertex costs increase linearly with 𝛾. We exploit this monotone structure using exact parametric Hochbaum pseudoflow [11]. Under the sparse compact-network conditions analyzed in Section 8, the per-iRM-set complexity is reduced from 𝑂 (|𝑉 (𝐷)| 3 ) to 2 log 𝑛 ), where 𝑛 ≤ |𝑉 (𝐷)|. 𝑂 (𝑛𝑀 𝑀 𝑀 As a summary, Algorithm 1 is the pseudocode for the overall framework of BoxDPpS. It first initializes the best feasible solution and the two sets used to record processed primitive keys and exact pruning certificates (lines 1–2). The bounded warm-up then improves the initial lower bound, after which several promising ratios are processed exactly to obtain early solutions and certificates (lines 3–4). The complete count domain is initially represented by a single box (line 5). During the main traversal, BoxDPpS repeatedly removes a box from the frontier (lines 6–7). If a safe upper bound or an exact certificate proves that the box cannot contain a better solution, the entire box is discarded without generating its individual count vectors (lines 8–9). For a surviving singleton box, BoxDPpS maps its count vector to the corresponding primitive key

Input : HIN G and query meta-path P Output : The densest 𝑃 -family Vbest 1 Initialize the best feasible family Vbest and its density 𝜌 best ; 2 Initialize the processed primitive keys K and exact certificates C; 3 Run the bounded warm-up and update 𝜌 best using its feasible families; // Section 7 4 Process the selected exact seeds and add any resulting exact certificates to C; // Sections 6, 7, and 8 5 Initialize the search frontier Q with the box X0 ; 6 while Q is not empty do 7 Remove a box X from Q; 8 if X is covered by a safe upper bound or an exact certificate then 9 Continue; // Section 6 10 11 12 13

14 15 16 17

else if X is a singleton {c} then c̄ ← Primitive(c); // Section 5 if c̄ ∉ K then Solve the fixed-M(c̄) problem exactly using the compact parametric solver; // Section 8 Update Vbest , 𝜌 best , K, and C; // Sections 6 and 8 else

Divide X into disjoint child boxes and insert them into Q;

return Vbest ;

(lines 10–11). If the key has not been processed, its fixed-M problem is solved using the compact parametric solver, and the resulting solution, key, and certificate are used to update the global state (lines 12–14). A surviving non-singleton box is divided into disjoint child boxes for further processing (line 16). When the frontier is exhausted, every admissible iRM-set has either been safely pruned or exactly processed through its unique primitive key. The algorithm therefore returns the globally optimal 𝑃-family (line 17).

5

PRIMITIVE COUNT VECTOR

§ 4 uses gcd normalization as the key of an iRM ratio. This section formally proves that the normalization is exact and quantifies how many duplicate count-vector representations it removes. For a positive integer count vector c = (𝑐 1, . . . , 𝑐𝑘 ), let 𝑔(c) = gcd(𝑐 1, . . . , 𝑐𝑘 ) and c̄ = c/𝑔(c). We call c̄ the primitive representative of c and let Dprim = {c ∈ D𝑐 | 𝑔(c) = 1}. This section keeps only the two properties needed by the algorithm: canonicalization is exact, and the number of canonical keys is smaller than the absolute count domain. We first prove that processing each primitive key once preserves the global DPpS optimum. Theorem 5.1: For positive vectors x and y, M(x) = M(y) if and only if y = 𝛼x for some 𝛼 > 0. Consequently, the mapping c ↦→ M(c) is a bijection from Dprim to the admissible iRM-set space M, and ∗ 𝜌 ∗ = maxc̄∈ Dprim 𝜌 M(c̄) .

Proof Sketch: If y = 𝛼x, then for every 𝑖, 𝑚𝑖 (y) = Î ( 𝑗 𝛼𝑥 𝑗 ) 1/𝑘 /(𝛼𝑥𝑖 ) = 𝑚𝑖 (x). Conversely, equality of the 𝑖-th iRM coordinates gives 𝑦𝑖 /𝑥𝑖 = 𝐺 (y)/𝐺 (x) for every 𝑖, so one positive scalar relates all coordinates. Every admissible iRM-set therefore has a primitive representative obtained by gcd normalization. If

Jiadong Xie, Jiaming Yang, Kangfei Zhao, and Jeffrey Xu Yu

two primitive integer vectors induce the same iRM-set, they are positive scalar multiples. Writing the scalar in lowest rational terms, integrality forces its denominator to divide every coordinate of the first vector, and primitivity forces that denominator to be one; primitivity of the second vector then forces the integer numerator to be one. Hence, the primitive representative is unique. The displayed optimization identity follows because the original exhaustive search is a maximum over M. □ We next consider the time complexity when utilizing the primitive count-vector canonicalization. Ímin𝑖 𝑁𝑖 Î j k Theorem 5.2: For 𝑘 ≥ 2, |M| = |Dprim | = 𝑑=1 𝜇 (𝑑) 𝑘𝑖=1 𝑁𝑑𝑖 , where 𝜇 is the Möbius function. If the 𝑁𝑖 grow proportionally, then Î Î|M| −→ 1 . Hence, canonicalization removes exactly 𝑖 𝑁𝑖 − |M| 𝜁 (𝑘 ) 𝑖 𝑁𝑖 duplicate absolute count vectors; asymptotically, the removed fraction is 1 − 1/𝜁 (𝑘).

Proof Sketch: Use the identity [gcd(𝑐 1, . . . , 𝑐𝑘 ) = 1] = Í the order 𝑑 |gcd(𝑐 1 ,...,𝑐𝑘 ) 𝜇 (𝑑) and sum it over D𝑐 . After exchanging Î of summation, a fixed 𝑑 contributes to exactly 𝑖 ⌊𝑁𝑖 /𝑑⌋ vectors. The first formula follows, and Theorem 5.1 identifies this count Î with |M|. Dividing by 𝑖 𝑁𝑖 and taking proportional limits gives Í 𝑘 □ 𝑑 ≥1 𝜇 (𝑑)/𝑑 = 1/𝜁 (𝑘).

A key is normalized in 𝑂 (𝑘 log 𝑁 max ) time by repeated gcd operations (or equivalently by a precomputed prime-divisor table), where 𝑁 max = max𝑖 ∈ [𝑘 ] 𝑁𝑖 . Canonicalization is performed only when a singleton or an exact seed is processed; it is not used to delete a non-singleton box, whose safe bounds remain count-specific. Primitive canonicalization removes duplicate fixed-M solves, but it does not by itself justify pruning a non-singleton count box. A box may contain count vectors from many different scaling rays, and its density bounds depend on the absolute counts represented by the box. Therefore, BoxDPpS applies canonicalization only when processing an exact seed or a surviving singleton, while all nonsingleton boxes remain in the original count coordinates and are pruned only by the safe conditions developed in § 6.

6

BOX-LEVEL IRM-SET SEARCH

§ 4 explains how BoxDPpS replaces point-wise iRM-set enumeration with a search over count boxes. This section formalizes the box frontier, the two box-level pruning mechanisms, and the correctness of the resulting search.

6.1

Box Frontier and Safe Deletion

For a box X as defined in § 4, a box is a singleton when ℓ𝑖 = ℎ𝑖 for every 𝑖. Let 𝜌 max (X) be the largest DPpS density among feasible families whose count vector lies in X. A function 𝑈 (X) is safe when 𝑈 (X) ≥ 𝜌 max (X).

Lemma 6.1: If 𝑈 (X) is a safe upper bound and 𝑈 (X) ≤ 𝜌 best , removing the entire box cannot remove a solution that improves the current incumbent. Proof Sketch: Every feasible family represented by the box has density at most 𝑈 (X), which is at most the density of an already

feasible family. Hence no family in the box can strictly improve the incumbent. □ For a non-singleton box, define its multiplicative span at dimension 𝑖 as 𝜎𝑖 (X) = log(ℎ𝑖 + 1) − log ℓ𝑖 . The implementation chooses 𝑟 ∈ arg max𝑖:ℓ𝑖 <ℎ𝑖 𝜎𝑖 (X), breaking ties by the raw width ℎ𝑖 − ℓ𝑖 . With 𝑞 = ⌊(ℓ𝑟 + ℎ𝑟 )/2⌋, it creates X − = {c ∈ X : 𝑐𝑟 ≤ 𝑞}, and X + = {c ∈ X : 𝑐𝑟 ≥ 𝑞 + 1}. All other coordinate intervals are unchanged. Lemma 6.2: For every non-singleton box X, the two children X − and X + constructed above satisfy (i) X − ∩ X + = ∅, (ii) X − ∪ X + = X, and (iii) repeated splitting terminates after finitely many operations at singleton boxes. Proof Sketch: Because ℓ𝑟 < ℎ𝑟 , the midpoint satisfies ℓ𝑟 ≤ 𝑞 < ℎ𝑟 . Every integer 𝑐𝑟 ∈ [ℓ𝑟 , ℎ𝑟 ] satisfies exactly one of 𝑐𝑟 ≤ 𝑞 and 𝑐𝑟 ≥ 𝑞+ 1, proving disjointness and union. Each child has a strictly smaller width in dimension 𝑟 , while no width increases. The nonnegative Í integer potential 𝑖 (ℎ𝑖 − ℓ𝑖 ) therefore strictly decreases along every root-to-leaf path and reaches zero exactly at a singleton. □ The active boxes form a disjoint frontier because the root represents the whole domain and every replacement obeys Lemma 6.2. We use a stack and insert the upper child before the lower child, so the lower-count child is examined next. This order can improve the incumbent earlier but has no role in correctness. Every box is tested with the current 𝜌 best and all certificates available when it is popped, so later improvements automatically strengthen the tests for the remaining frontier.

6.2

Safe Box Upper Bounds

All box-level statistics are computed from the complete queryinduced instance multiset F0 = FP . An exact weighted representation may be used internally, but each stored weight equals the number of original instances it represents. These summaries are independent of the processed box and are reused throughout the traversal. We combine several complementary safe upper bounds. By Lemma 6.1, each bound can safely prune X when 𝑈 (X) ≤ 𝜌 best . Total-instance Bound: Let F0 be the set of all meta-path instances in the current search universe, with multiplicity included when compression is used. No family can contain more than |F0 | total instance weight, and every c ∈ X satisfies 𝑐𝑖 ≥ ℓ𝑖 . Thus 𝑈 all (X) = | F0 | is safe. Î𝑘 1/𝑘 ( 𝑖=1 ℓ𝑖 )

Degree Bound: Let 𝑑𝑖0 (𝑣) = 𝑑𝑖 (𝑣; F0 ). Sort vertices at position 𝑖 by non-increasing 𝑑𝑖0 (𝑣) and let 𝐷𝑖 (𝑡) be the prefix sum of the largest 𝑡 degrees. For any family with |𝑉𝑖 | = 𝑡, its induced instance count is at most 𝐷𝑖 (𝑡). Therefore, 𝑈 deg,𝑖 (X) = maxℓ𝑖 ≤𝑡 ≤ℎ𝑖 (𝑡 Î𝐷𝑖 (𝑡ℓ ) ) 1/𝑘 is safe for position 𝑖, and 𝑈 deg (X) = min𝑖 𝑈 deg,𝑖 (X) is safe.

𝑗 ≠𝑖 𝑗

Projection-aware Bound: For a position subset 𝑆 ⊂ [𝑘] with 2 ≤ |𝑆 | < 𝑘, group instances by the projection 𝜋𝑆 (𝑝) = (𝑝𝑖 )𝑖 ∈𝑆 . Let 𝑤 𝑆 (𝑎) be the total weight of group 𝑎, and let 𝑇𝑆 (𝑞) be the sum of the largest 𝑞 group weights. A family in X can activate at most Î keys 𝑖 ∈𝑆 ℎ𝑖 such projection keys. Counting all paths under those Î 𝑇 ( ℎ ) only overestimates the numerator, so 𝑈 proj,𝑆 (X) = 𝑆( Î 𝑖ℓ∈𝑆) 1/𝑘𝑖 is 𝑖 𝑖 safe.

Scalable Exact Densest 𝑃 -Partite Subgraph Search in Heterogeneous Information Networks

Pair-projection bound. For positions (𝑟, 𝑐), form a weighted ma𝑟,𝑐 trix 𝐴𝑟,𝑐 where 𝐴𝑢𝑣 is the number of meta-path instances with Í Í 𝑟,𝑐 𝑝𝑟 = 𝑢 and 𝑝𝑐 = 𝑣. Let 𝐻𝑟,𝑐 (𝑎, 𝑏) = max |𝑅 | ≤𝑎,|𝐶 | ≤𝑏 𝑢 ∈𝑅 𝑣 ∈𝐶 𝐴𝑢𝑣 . 𝐻𝑟,𝑐 (ℎ𝑟 ,ℎ𝑐 ) Then 𝑈 pair,𝑟,𝑐 (X) = ( Î ℓ ) 1/𝑘 is safe because every family in the 𝑖 𝑖 box selects at most ℎ𝑟 rows and ℎ𝑐 columns at the two positions. Theorem 6.1: Let 𝑈 (X) be the minimum of any subset of the totalinstance, degree, projection-aware, and pair-projection bounds above. Then 𝑈 (X) is a safe upper bound for X.

Proof Sketch: The total-instance bound uses the global numerator upper bound and the minimum possible denominator in the box. The degree bound follows because each induced path contributes to exactly one selected vertex at any fixed position, so the numerator is at most the sum of the largest selected-position degrees. The projection-aware bound relaxes the family to choose the heaviest possible projection keys, and the pair-projection bound further relaxes it to choose the heaviest row/column submatrix. Each relaxation can only increase the numerator and each denominator is Î lower-bounded by ( 𝑖 ℓ𝑖 ) 1/𝑘 . Therefore, every listed bound is safe, and its minimum remains safe. □

6.3

Lifting Exact iRM Certificates to Boxes

A completed exact fixed-M𝑎 solve may generate a pointwise pruning Í 𝑚𝑎 condition of the form 𝑘 ≤ 𝑘𝑖=1 𝑚𝑖𝑖 ≤ 𝐵𝑎 , under the premises of the exact iRM pruning lemmas. We use such a condition only when it is backed by a strict full exact For a count vector c, Í certificate. Í 𝑚𝑖𝑎 𝑚𝑎 𝑐 𝑖 define 𝐶𝑎 (c) = 𝑘𝑖=1 𝑚𝑖 (c) = ( Î𝑖 𝑐 𝑖) 1/𝑘 . The lower bound 𝐶𝑎 (c) ≥ 𝑘 𝑗

𝑗

follows from AM–GM because both iRM-sets have product one. To certify an entire box, we maximize 𝐶𝑎 over a log-box relaxation. Let 𝑐𝑖 = 𝑒 𝑦𝑖 with 𝑦𝑖 ∈ [log ℓ𝑖 , log ℎ𝑖 ]. Then 𝐶𝑎 (𝑒 y ) = Í𝑘 1 Í𝑘 𝑎 𝑖=1 𝑚𝑖 exp 𝑦𝑖 − 𝑘 𝑗=1 𝑦 𝑗 , which is a convex function of y.

Theorem 6.2: Given a valid point-wise certificate (M𝑎 , 𝐵𝑎 ), if every corner of the log-box satisfies 𝐶𝑎 (𝑒 y ) ≤ 𝐵𝑎 , then every integer count vector in X satisfies the pointwise iRM pruning condition, and the whole box can be safely pruned. Proof Sketch: The log-box is the convex hull of its 2𝑘 corners. Since 𝐶𝑎 (𝑒 y ) is convex in y, Jensen’s inequality implies that its maximum over the log-box is attained at a corner. Thus all points in the box satisfy 𝐶𝑎 ≤ 𝐵𝑎 . The lower inequality 𝐶𝑎 ≥ 𝑘 holds by AM–GM. The strict exact certificate then proves that no iRM-set induced by the box can improve 𝜌 best . □

6.4

Processing a Surviving Canonical Key

A singleton reached by the frontier and an exact seed produced by warm-up are handled identically. Algorithm 2 first canonicalizes the count vector, suppresses a duplicate key, invokes the exact compact fixed-M solver, and updates the two pieces of global state that can change: the feasible incumbent and the strict-certificate set. The gcd normalization in line 1 maps every point on a scaling ray to the unique key proved in Theorem 5.1; the membership test then prevents a second solve for that iRM-set (line 2). The iRM weights are materialized only for a new key (line 3), which is passed to the exact compact solver in line 4. The call returns the best feasible

Algorithm 2: ProcessKey (c, G, P, K, C, Vbest, 𝜌 best )

: Count vector c, query ( G, P ), and mutable state ( K, C, Vbest , 𝜌 best ) Output : Updated mutable state ( K, C, Vbest , 𝜌 best ) 1 c̄ ← c/gcd(𝑐 1 , . . . , 𝑐 𝑘 ); 2 if c̄ ∈ K then return; 3 M ← M(c̄); c 𝛾 final , cert) ← FixMSolve( G, P, M, 𝜌 best ); 4 ( V, // Algorithm 4 c ≠ ∅ and 𝜌 ( V c) > 𝜌 best then Vbest ← V; c 𝜌 best ← 𝜌 ( V c); 5 if V 6 insert c̄ into K; 7 if cert is a strict full exact certificate then add cert to C; Input

family it discovers, the final threshold of the completed exact fixedM iteration, and an optional cross-ratio certificate. Only a feasible family can raise the global incumbent (line 5). After the exact call, the key is memoized (line 6). A certificate enters C only when the full exact premises of the original iRM pruning result hold (line 7); otherwise cert = ⊥ and no cross-ratio deletion is enabled.

6.5

Correctness of the Box Search

Theorem 6.3: Assume that every box deletion uses either a safe upper bound or a strict full exact iRM certificate lifted by Theorem 6.2, and that every surviving singleton is normalized according to Theorem 5.1 and processed by Algorithm 2. Then the box search preserves the global DPpS optimum. Proof Sketch: Let V ∗ be a globally optimal 𝑃-family, with count vector c∗ = c(V ∗ ). Consider the unique chain of frontier boxes containing c∗ . While the incumbent is strictly smaller than 𝜌 ∗ , no box on this chain can be removed by a safe upper bound, because that bound must be at least 𝜌 (V ∗ ) = 𝜌 ∗ . Nor can the chain be removed by a lifted strict full exact certificate, because such a certificate proves that every represented iRM-set cannot improve the incumbent, whereas M(c∗ ) can improve it. Hence, unless the incumbent is already optimal, the chain is divided down to the singleton {c∗ }. At that singleton, Algorithm 2 computes the canonical key c̄∗ . If the key is new, Algorithm 4 solves its fixed-M problem exactly and therefore considers V ∗ among its feasible families. If the key is already in K, the same iRM-set was exactly processed earlier from another point on the scaling ray. In either case, the incumbent has already reached 𝜌 ∗ or is raised to 𝜌 ∗ . Therefore, the box search cannot lose the global optimum. □

7

REPRESENTATIVE-RATIO WARM-UP

§ 4 uses warm-up only to obtain an early feasible lower bound and schedule a few exact keys. This section specifies that bounded procedure. Let 1 = (1, . . . , 1) and N = (𝑁 1, . . . , 𝑁𝑘 ). For position 𝑖, let a𝑖 (𝑞) have coordinate 𝑞 at 𝑖 and one elsewhere, and √ define 𝑄𝑖 = 𝑁𝑖 , max{1, ⌊𝑁𝑖 /2⌋}, max{1, ⌊𝑁𝑖 /4⌋}, max{1, round( 𝑁𝑖 )} . The budget 𝐵 𝑤 bounds approximate evaluations, and 𝐵𝑒 bounds exact seeds. Algorithm 3 shows the details. The deterministic representatives in line 1 cover balanced, full-size, and axis-dominant regimes. Lines 2–6 initialize the queue and remove proportional duplicates before they consume budget. If these candidates and their adaptive descendants are exhausted, the lower-child-first singleton iterator

Jiadong Xie, Jiaming Yang, Kangfei Zhao, and Jeffrey Xu Yu

Algorithm 3: Warm − Up (G, P, X0, 𝐵 𝑤 , 𝐵𝑒 )

Input : HIN G, meta-path P, root box X0 , and budgets 𝐵 𝑤 , 𝐵𝑒 Output : Warm-up records R and ordered exact-seed keys S 1 B ← {1, N} ∪ {a𝑖 (𝑞) : 𝑞 ∈ 𝑄 𝑖 , 𝑖 ∈ [𝑘 ] }; 2 Q𝑤 ← an empty FIFO queue; H ← ∅; R ← [ ]; 3 Ileaf ← a lazy lower-child-first singleton iterator over X0 ; 4 for each c ∈ B in the displayed order do 5 c̄ ← Primitive(c); 6 if c̄ ∉ H then insert c̄ into H and enqueue it into Q𝑤 ;

𝑗 ← 0; while 𝑗 < 𝐵 𝑤 do 9 while Q𝑤 = ∅ and Ileaf has a next singleton do 10 c′ ← next( Ileaf ); c̄′ ← Primitive(c′ ); 11 if c̄′ ∉ H then insert c̄′ into H and enqueue it into Q𝑤 ; 7

8

12 13 14 15 16 17 18 19 20 21 22 23 24

if Q𝑤 = ∅ then break; c̄ ← dequeue( Q𝑤 ); 𝑠 ← 0; f Run the approximation algorithm in [4] on M(c̄) and obtain V; f ≠ ∅ then if V f); b f); A ← {b 𝑠 ← 𝜌 (V c ← c( V c}; for 𝑖 ← 1 to 𝑘 do c− ← b c; 𝑐𝑖− ← max{1, ⌊0.8b 𝑐𝑖 ⌋ }; c+ ← b c; 𝑐𝑖+ ← min{𝑁𝑖 , ⌈1.25b 𝑐𝑖 ⌉ }; insert c − and c+ into A;

for each c′ ∈ A do c̄′ ← Primitive(c′ ); if c̄′ ∉ H then insert c̄′ into H and enqueue it into Q𝑤 ;

Algorithm 4: FixMSolve (G, P, M, 𝛾 0 )

: HIN G, meta-path P, iRM-set M, and certified initial threshold 𝛾 0 c final threshold 𝛾 final , and optional Output : Best feasible family V, strict certificate cert 1 F𝜔 ← replace every verified terminal-twin class by one weighted representative; 2 (𝐿, 𝑅) ← arg min over admissible nonempty bipartitions of the projection-size estimate; 3 {𝑤𝑎𝑏 }, {𝑊𝑎 },𝑊all ← aggregate F𝜔 by (𝑎, 𝑏 ) = (𝜋𝐿 (𝑝 ), 𝜋𝑅 (𝑝 ) ); 4 Construct the grouped closure 𝐷 M (𝛾 ) with capacities 𝑊𝑎 , 𝑤𝑎𝑏 , implication arcs, and 𝑐 𝑣 (𝛾 ) = 𝜔 (𝑣)𝑚 pos(𝑣) 𝛾 /𝑘; 5 𝑈 𝑀 ← 𝑊all ; 6 Invoke the exact parametric HPF algorithm [11] on 𝐷 M (𝛾 ) for 𝛾 ∈ [0, 𝑈𝑀 ], using 𝛾 0 only as a pivot, and store the results in BM ; c ← ∅; 𝛾 ← 𝛾 0 ; 7 V 8 while 𝛾 ≤ 𝑈 𝑀 do 9 𝑆 ← the canonical source set in BM that is optimal at 𝛾 ; 10 V𝑆 ← recover the original weighted 𝑃 -family represented by 𝑆; Í 11 𝑧 ← | F ( V𝑆 ) | − (𝛾 /𝑘 ) 𝑖 𝑚𝑖 |𝑉𝑖𝑆 |; 12 if V𝑆 = ∅ or 𝑧 ≤ 0 then break; c = ∅ or 𝜌 ( V𝑆 ) > 𝜌 ( V c) then V c ← V𝑆 ; 13 if V 14 𝛾 ← 𝜌 ( V𝑆 ); Input

15 𝛾 final ← 𝛾 ;

cert ← the strict full exact iRM certificate implied by the completed iteration, or ⊥ when its premises do not hold; c 𝛾 final , cert); 17 return ( V, 16

f 𝑠, 𝑗 ) to R; 𝑗 ← 𝑗 + 1; append (c̄, V,

sort R by non-increasing 𝑠; S ← the keys in the first min{𝐵𝑒 , | R | } records; 27 return ( R, S);

25

8

26

§ 4 outlined the multiplicity-aware solver. This section gives its algorithm and proves only the two representation equivalences and the HPF bound needed by the global correctness argument. The call receives the feasible threshold 𝛾 0 = 𝜌 best , used only as a pivot or warm start. Let 𝑊all be the total aggregated path weight; 𝑈𝑀 = 𝑊all is a safe upper endpoint because every nonempty family has a geometric-mean denominator at least one.

supplies unseen primitive keys without invoking the certificatebased cross-iRM-set pruning rule (lines 9–11); line 12 stops when no key remains. Each dequeued key is evaluated by the existing peeling approximation (lines 13–14). A returned witness is scored by its actual feasible density and exposes an observed count vector (lines 13–16). The observed vector and clipped one-coordinate 0.8 contractions and 1.25 expansions form the adaptive neighborhood (lines 16–20); they are again canonicalized and deduplicated before enqueueing (lines 21–23). Finally, records are ranked by feasible density and the first 𝐵𝑒 keys become exact seeds (lines 25–27). Theorem 7.1: Algorithm 3 does not affect exactness provided that an approximate witness updates 𝜌 best only through its feasible density, the approximate output creates no certificate, and every selected seed is passed to Algorithm 2. Proof Sketch: A feasible witness is a lower bound on 𝜌 ∗ , so replacing the incumbent by its density cannot make any safe upper-bound comparison unsound. Candidate ranking changes only processing order. Since approximation emits no certificate, every deletion is still justified by the same safe bound or strict exact certificate as without warm-up. Exact seeds are ordinary canonical keys processed earlier, and memoization merely prevents a duplicate solve later. □

COMPACT FIXED-M SOLVING

Overall Procedure: Algorithm 4 shows the overall procedure. Line 1 replaces only terminal-twin classes that satisfy the exact condition of Theorem 8.1. It then selects a projection bipartition by an estimated compact-network size (line 2) and merges all instances with the same endpoint projections into one weighted pair (line 3). These steps change representation and multiplicity, not the objective. The grouped network is built once with constant path-profit and implication capacities and linear vertex-to-sink capacities (line 4). The total path weight gives a safe parameter endpoint (line 5). The exact HPF call materializes the nested source-set representation over the certified range and uses 𝛾 0 only as a pivot (line 6); hence, the incumbent can guide execution without removing a parameter regime needed by the completed exact solve. The routine starts from the certified feasible threshold (line 7). At every iteration, it retrieves the source set that minimizes the grouped cut at the current threshold, recovers the represented family, and evaluates the corresponding auxiliary objective (lines 8–11). A nonpositive value proves that no conforming family can improve the current threshold, so the iteration stops (line 12). Otherwise, the feasible witness may replace the local best family and strictly raises

Scalable Exact Densest 𝑃 -Partite Subgraph Search in Heterogeneous Information Networks

𝛾 to its actual density (lines 13–14). The final threshold is recorded after termination (line 15). A cross-ratio certificate is emitted only if the full exact premises of the existing iRM pruning result hold (line 16); otherwise, it is ⊥. The completed result is returned to Algorithm 2 (line 17).

8.1

Terminal-Twin Compression

At terminal position 𝑘, vertices 𝑢 and 𝑣 are terminal twins when they have the same predecessor multiset at position 𝑘 − 1 and occur only at the terminal position. A class 𝑇 is represented by one node with multiplicity 𝜔 (𝑣ˆ𝑇 ) = |𝑇 |.

Theorem 8.1: Replacing every terminal-twin class by one representative whose vertex cost and incident path multiplicities are multiplied by its class size preserves both the global DPpS optimum and every fixed-M auxiliary optimum.

Proof Sketch: Fix all selections outside one class 𝑇 and let 𝑞 be the number of selected twins. Since all twins complete the same predecessor multiset, the induced path contribution is affine in 𝑞, say 𝐴+𝐵𝑞. The fixed-M vertex penalty is also affine, so the auxiliary objective is affine and has an optimum at 𝑞 = 0 or 𝑞 = |𝑇 |. For the density objective, let 𝑏 be the number of already selected non-twin terminal vertices and let 𝐶 be the product of the fixed counts at the other positions. If 𝑏 > 0, the density as a continuous function of 𝑞 𝐴+𝐵𝑞 is 𝑓 (𝑞) = (𝐶 (𝑏+𝑞) ) 1/𝑘 . After removing the positive common factor, the sign of 𝑓 ′ (𝑞) is the sign of 𝐵(𝑏 + 𝑞) − (𝐴 + 𝐵𝑞)/𝑘, an affine nondecreasing function of 𝑞. Thus 𝑓 can be monotone or decrease and then increase, but cannot have a strict interior maximum; an endpoint is optimal. If 𝑏 = 0, feasibility implies 𝐴 = 0, and 𝑓 (𝑞) is proportional to 𝑞 1−1/𝑘 for 𝑞 ≥ 1, so selecting all twins is optimal. Hence, some global optimum is also all-or-none on every class. The weighted representative has exactly these two states. Its multiplied path gains and vertex cost equal those of selecting all original twins, while its unselected state equals selecting none. Applying the argument class by class preserves both optimum values. □

8.2

Projection-Grouped Closure Network

After terminal compression, choose a nonempty bipartition 𝐿 ∪ 𝑅 = [𝑘]. For an instance 𝑝, let 𝑎 = 𝜋𝐿 (𝑝) and 𝑏 = 𝜋𝑅 (𝑝); all instances with the same pair are aggregated into weight 𝑤 𝑎𝑏 , and 𝑊𝑎 = Í 𝑏 𝑤 𝑎𝑏 . The implementation chooses the admissible bipartition with the smallest estimated grouped-network size; this changes representation only. The grouped closure contains source 𝑠, sink 𝑡, compressed vertex nodes, and one node for each distinct left and right projection. It has arcs 𝑠 → 𝑎 of capacity 𝑊𝑎 , arcs 𝑎 → 𝑏 of capacity 𝑤 𝑎𝑏 , infinitecapacity implications from each projection node to its member vertices, and vertex-to-sink arcs of capacity 𝜔 (𝑣)𝛾𝑚𝑖 /𝑘. A singletonside projection node may be reused as its vertex node. Theorem 8.2: For any fixed M and threshold 𝛾, the minimum cut in the projection-grouped closure maximizes exactly the same auxiliary objective 𝜁 (V, 𝛾, M) as the standard one-path-node network.

Proof Sketch: Fix a set 𝑆 of selected original vertices. Infinitecapacity implications allow a projection node on the source side only when all of its member vertices are in 𝑆. For a compatible left

projection 𝑎, placing 𝑎 on the source side cuts only weights leading to incompatible right projections, which costs no more than cutting 𝑠 → 𝑎 with capacity 𝑊𝑎 ; therefore an optimal finite cut places every compatible left projection on the source side. The same argument applies to right projections. Consequently, the only path weights not cut are those whose two projections are compatible with 𝑆, and the Í path-related cut cost is 𝑊all − 𝑎⊆𝑆, 𝑏 ⊆𝑆 𝑤 𝑎𝑏 . The vertex-to-sink arcs Í add (𝛾/𝑘) 𝑖 𝑚𝑖 |𝑉𝑖 |, with terminal multiplicities already preserved by Theorem 8.1. Hence the complete cut value is 𝑊all − 𝜁 (V, 𝛾, M), exactly as in the standard construction. Minimizing the cut therefore maximizes the same auxiliary objective. □ Network Size: If 𝑛𝑐 is the number of compressed vertices, 𝑔𝐿 and 𝑔𝑅 are the numbers of distinct left and right projections, and 𝑚𝐿𝑅 is the number of nonzero projection pairs, then the grouped network has 𝑂 (𝑛𝑐 + 𝑔𝐿 + 𝑔𝑅 ) nodes and 𝑂 (𝑛𝑐 + 𝑘𝑔𝐿 + 𝑘𝑔𝑅 + 𝑚𝐿𝑅 ) arcs. In contrast, the standard network has |F0 | path nodes and Θ(𝑘 |F0 |) path-related arcs. The grouped construction therefore depends on distinct signatures rather than total path multiplicity.

8.3

HPF-Compatible Parametric Solving

For fixed M, all path-profit and implication capacities in the grouped closure are constant. Only vertex-to-sink capacities vary with the 𝑚 𝛾 . These capacities are linear density threshold: 𝑐 𝑣 (𝛾) = 𝜔 (𝑣) pos(𝑣) 𝑘 and monotone in 𝛾. Therefore, the compact closure satisfies the monotone parametric cut condition required by HPF [11]. Theorem 8.3: Let 𝑛𝑀 and 𝑚𝑀 be the number of nodes and arcs in the grouped closure network for a fixed M, and let 𝑏 𝑀 be the number of materialized breakpoints. A fully-parametric solves   HPF backend   𝑛2

the monotone fixed-M closure in 𝑂 𝑚𝑀 𝑛𝑀 log 2 + 𝑚𝑀𝑀 + 𝑏 𝑀 𝑛𝑀 .

If the grouped closure is sparse, 𝑚𝑀 = 𝑂 (𝑛𝑀 ), and 𝑏 𝑀 ≤ 𝑛𝑀 , this 2 log 𝑛 ). becomes 𝑂 (𝑛𝑀 𝑀

Proof Sketch: The monotonicity condition follows from the linear nondecreasing vertex-to-sink capacities above. Hochbaum’s pseudoflow analysis gives a parametric cut bound of the same asymptotic core order as one pseudoflow run; with the dynamic-tree
2 /𝑚 refinement this is 𝑂 (𝑚𝑀 𝑛𝑀 log 2 + 𝑛𝑀 𝑀 ). If the implementation explicitly materializes 𝑏 𝑀 breakpoint cuts and scans 𝑛𝑀 nodes for each, an output cost 𝑂 (𝑏 𝑀 𝑛𝑀 ) is added. In a monotone parametric cut, canonical source sets are nested, and each strict breakpoint changes at least one node side, so 𝑏 𝑀 ≤ 𝑛𝑀 . Substituting 2 log 𝑛 ). 𝑚𝑀 = 𝑂 (𝑛𝑀 ) yields 𝑂 (𝑛𝑀 □ 𝑀 By comparison, the original DPpS analysis uses the cubic per-M bound 𝑂 (|𝑉 (𝐷)| 3 ) for its explicit flow network 𝐷. Thus, under the sparse grouped-closure conditions of Theorem 8.3, the theoretical dependence on the compact network becomes near quadratic with one logarithmic factor, in addition to the reduction from |𝑉 (𝐷)| to 𝑛𝑀 . This statement does not imply that every HPF binary is faster in practice; library binding, memory copying, and breakpoint materialization can dominate on small or moderately sized networks.

Jiadong Xie, Jiaming Yang, Kangfei Zhao, and Jeffrey Xu Yu

Table 2: Dataset Statistics |A| 5 5 6 417 21,218 4 4

|R| 4 4 6 634 1,152 3 3

|𝑉 | 2,672 37,795 37,595 4,767,652 36,791,318 2,405 47,490

AdvExactGVIt

|𝐸 | 104,747 174,851 1,714,941 12,160,769 94,846,085 7,704 110,730

104 103 102 101 100 10−1

ML DBLP DBN DBP

FB

(a) 𝑘 = 3

9

G21 G22

ML DBLP DBN DBP

FB

(b) 𝑘 = 4

G21 G22

EXPERIMENTAL STUDIES

We evaluate BoxDPpS from four perspectives: ➊ Compare BoxDPpS end-to-end efficiency with AdvExactGVIt; ➋ Evaluate the returned subgraphs under the DPpS density objective and in a cybersecurity user-grouping application; ➌ Quantify the contribution of each proposed technique through an ablation study; ➍ Examine scalability with respect to both data size and meta-path length. We first introduce the datasets, query workloads, compared methods, and experimental environment.

9.1

BoxDPpS

OOM/OOT 105

Time (s)

Dataset (Abbr. Name) MovieLens (ML) DBLP Douban (DBN) DBpedia (DBP) Freebase (FB) Cisco g21 (G21) Cisco g22 (G22)

Experimental Setup

Datasets. We use seven real-world HIN datasets from different application domains: MovieLens, DBLP, and Douban [1], DBpedia and Freebase [35], and Cisco g21 and Cisco g22 [24]. The first five datasets are used for the general performance evaluation, while Cisco g21 and g22 contain ground-truth user groups for the cybersecurity evaluation. Table 2 summarizes the numbers of vertex types, relation types, vertices, and edges in each dataset. Query meta-paths. For each dataset and each meta-path length considered in an experiment, we generate candidate meta-paths with distinct vertex-type positions. We use 20 query meta-paths when at least 20 candidates are available and retain all candidates otherwise. Compared Algorithms. We compare our exact framework, BoxDPpS, with AdvExactGVIt [4], whose source code implementation is provided by the authors. For the solution-quality evaluation, we reimplement MAvgP and iBF following the multipartite baseline formulations [4]. MAvgP maximizes the average number of metapath instances, whereas iBF adapts the butterfly-core model [6] to a query meta-path. Each algorithm is run five times on every selected query, and we report the average running time. All compared methods are evaluated using the same queries. Environment. All algorithms are implemented in C++17 and compiled using g++ 11.4.0 with the -O3 optimization flag. Experiments are conducted in single-thread mode on an Ubuntu 22.04.5 server equipped with two AMD EPYC 9554 processors and 512 GB of main memory. We exclude only graph loading from reported running time; all algorithm-specific preprocessing and solving costs are included. A run that exceeds 48-hour limit is marked as OOT, whereas a run terminated by memory exhaustion is marked as OOM. Parameter Settings. For the iRM-set pruning rules in Section 3, we use 10−4 as the certificate margin and retain a certificate only when 𝐵𝑎 − 𝑘 > 10−4 . For the representative-ratio warm-up in Section 7, we set 𝐵 𝑤 = 64 and use 𝐵𝑒 = 4 by default, increasing it to 8 when Î𝑘 8 𝑖=1 𝑁𝑖 ≥ 10 or max𝑖 𝑁𝑖 ≥ 5000.

Figure 3: Running Time of BoxDPpS and AdvExactGVIt Table 3: Normalized Density Gain of the Compared Methods Method

MAvgP iBF AdvExactGVIt BoxDPpS

9.2

MovieLens 𝑘 =3 𝑘 =4 1.27 1.05 1.03 1.00 1.46 1.31 1.46 1.31

DBLP 𝑘 =3 𝑘 =4 0.77 0.37 1.13 1.00 1.94 1.37 1.96 1.38

DBpedia 𝑘 =3 𝑘 =4 2.12 2.00 2.71 1.61 5.80 5.08 5.82 5.15

Freebase 𝑘 =3 𝑘 =4 7.90 2.23 4.98 1.37 32.90 13.93 33.21 14.01

Experimental Results

Exp-1: Overall Efficiency. We compare the end-to-end running time of BoxDPpS and AdvExactGVIt on all seven datasets with 𝑘 = 3 and 𝑘 = 4. Fig. 3 reports the running-time results. BoxDPpS substantially outperforms AdvExactGVIt on larger graphs. Among the query meta-paths completed by both algorithms, BoxDPpS achieves an average speedup of 27.04×. On DBLP, DBpedia, and Freebase, the speedups range from 3.59× to 40.74× for 𝑘 = 3 and from 22.61× to 131.77× for 𝑘 = 4. Moreover, AdvExactGVIt times out on both Douban settings and runs out of memory on Cisco g22 when 𝑘 = 4, whereas BoxDPpS completes all these workloads. On the small Cisco g21 graph, AdvExactGVIt is faster because the additional search and preprocessing overhead of BoxDPpS is not fully amortized. Overall, the advantage becomes more remarkable as the query workload grows, demonstrating the effectiveness of box-based search and compact fixed-M solving. Exp-2: Effectiveness of Normalized Density Gain. We next evaluate the densities of the subgraphs returned by the compared methods. Since absolute densities are not directly comparable across query meta-paths, we measure the density gain over the complete query-induced 𝑃-partite graph. For an algorithm 𝑎 and a query meta-path 𝑝, let 𝜌𝑎,𝑝 denote the density of the returned subgraph, and let 𝜌 whole,𝑝 denote the density before any vertices are removed. We first compute the relative density 𝑟 𝑎,𝑝 = 𝜌𝑎,𝑝 /𝜌 whole,𝑝 . For a group 𝐺 consisting of one dataset and one value of 𝑘, let P𝐺 be its set of query meta-paths. We report their geometric mean as
the Í normalized density gain NGD𝑎,𝐺 = exp | P1𝐺 | 𝑝 ∈ P𝐺 log 𝑟 𝑎,𝑝 . Table 3 shows that BoxDPpS matches AdvExactGVIt on MovieLens after rounding and yields slightly larger density gains in all six datasets. For example, on DBpedia, the gains increase from 5.80 to 5.82 for 𝑘 = 3 and from 5.08 to 5.15 for 𝑘 = 4. Both methods optimize the same DPpS density and should therefore return the same best density. The small difference in Table 3 comes from floating-point threshold updates in the implementation of AdvExactGVIt. The resulting precision loss can stop its fixed-M iteration before the best

Scalable Exact Densest 𝑃 -Partite Subgraph Search in Heterogeneous Information Networks Papers (18)

Authors (5) Venues (3) Jeffrey Xu Yu ICDE Lu Qin

Papers (20)

Authors (5) Venues (4) Jeffrey Xu Yu ICDE Lu Qin

PVLDB Lijun Chang Xuemin Lin

PVLDB Lijun Chang

SIGMOD

Xuemin Lin

Table 4: F1 Scores on G21 under Different Packet Thresholds Method MAvgP iBF AdvExactGVIt BoxDPpS

1 0.435 0.435 0.615 0.615

10 0.435 0.435 0.615 0.615

100 0.435 0.435 0.615 0.615

1,000 0.375 0.402 0.519 0.519

5,000 0.241 0.402 0.505 0.633

10,000 0.248 0.402 0.514 0.514

SIGMOD

Table 5: Arithmetic Mean of Per-query Running Times (s) Hong Cheng

Hong Cheng KDD

(a) iBF, 𝜌 = 25.529

(b) BoxDPpS, 𝜌 = 26.873

Figure 4: Case study on a DBLP Paper–Author–Venue query

density is reached. We verify the solutions returned by BoxDPpS using exact pruning and fixed-M certificates, indicating that BoxDPpS recovers solutions missed by the reference implementation. Compared with MAvgP and iBF, BoxDPpS achieves higher NGD because it directly optimizes the DPpS density. MAvgP maximizes the average number of meta-path instances without normalizing by the geometric mean of the selected vertex-set sizes, whereas iBF focuses on butterfly-core cohesion and local support. Their returned subgraphs therefore need not maximize the density of complete meta-path instances over all positions of the query meta-path. This difference is more evident on the larger and more heterogeneous DBpedia and Freebase graphs. Exp-3: Effectiveness of Qualitative Case Study. As a qualitative case study, Fig. 4 visualizes BoxDPpS and iBF on an ego-centric Paper–Author–Venue subgraph around author ‘Jeffrey Xu Yu’, extracted from 2010–2024 DBLP. The subgraphs are searched from DBLP records from 2010 to 2024 for the four venues: SIGMOD, VLDB, ICDE, and KDD. We report the iBF result with the highest DPpS density over all butterfly thresholds. Both methods deliver the same five-author core. However, BoxDPpS additionally retains the venue KDD and two papers, adding 33 meta-path instances and increasing the density from 25.529 to 26.873. This case study illustrates that vertices with limited local support may still jointly improve the global multipartite density. Exp-4: Cybersecurity User Grouping. We further evaluate application-level group quality on Cisco g21, which contains nonoverlapping ground-truth user groups in a User–Port–Protocol– Server HIN. We repeatedly extract the users in the densest 𝑃-partite subgraph as a discovered group and report the average F1 score against the ground-truth groups. We vary the packet threshold 𝜏 and retain communication relations involving at least 𝜏 packets. Table 4 shows that BoxDPpS matches AdvExactGVIt at five of the six packet thresholds and improves the F1 from 0.505 to 0.633 at 𝜏 = 5,000. These results indicate that the efficiency-oriented transformations of BoxDPpS do not sacrifice application-level quality. Exp-5: Ablation Studies. To isolate the effect of each proposed component, we conduct the ablation study on a controlled variant of BoxDPpS that uses the same network-flow algorithm as AdvExactGVIt for all fixed-M subproblems. The configuration with all evaluated components enabled is denoted by Full. Each ablated variant disables one component from Full while keeping all other components unchanged. Specifically, No-Prim disables primitive

Method

Full No-Prim No-BoxUB No-W/S No-Group

MovieLens 𝑘 =3 𝑘 =4 11.70 15.36 12.08 17.23 11.73 16.10 12.07 17.97 20.15 96.32

DBLP 𝑘 =3 𝑘 =4 23.31 19.54 55.50 31.61 30.15 24.48 23.62 441.80 101.42 51.15

DBpedia 𝑘 =3 𝑘 =4 22.48 68.32 27.84 93.89 26.10 104.66 24.54 221.29 46.21 79.05

Freebase 𝑘 =3 𝑘 =4 812.20 188.77 923.44 235.87 1152.36 299.52 1077.63 3045.11 1408.02 299.78

count-vector canonicalization; No-BoxUB disables all box-level upper bounds; No-W/S disables both warm-up and exact seed scheduling; and No-Group disables projection grouping in fixed-M solving. Table 5 reports the average running time over the meta-path queries for each dataset and each value 𝑘 ∈ {3, 4}. To summarize the overall effect of a component, we additionally compare the accumulated running time of its variant with that of Full. The three search-side techniques prune the search space at different stages. Warm-up and exact seed scheduling make useful lower bounds and strict certificates available early in the search. Their effect is limited on most 𝑘 = 3 queries but becomes remarkable on the longer 𝑘 = 4 queries: No-W/S is the slowest variant on DBLP, DBpedia, and Freebase for 𝑘 = 4. Aggregated over all evaluated queries, disabling the two techniques increases the running time to 3.76× that of Full. Box upper bounds subsequently use the current lower bound to discard count-vector regions without examining their members individually, whereas primitive canonicalization avoids reprocessing count vectors that represent the same iRM-set. The more moderate slowdowns of No-BoxUB and No-Prim in Table 5 reflect their complementary roles: region-level pruning removes unpromising portions of the search space, while canonicalization removes duplicate representations. A separate source of cost arises after a count vector survives the search-side pruning and requires an exact fixed-M solve. We observe the benefit of projection grouping across all four datasets for 𝑘 = 3, which remains substantial on several 𝑘 = 4 queries. Disabling projection grouping increases the accumulated running time to 1.70× that of Full. This result confirms that reducing the networkflow instance for each remaining exact solve is important, even after the number of examined count vectors has been reduced. Although some individual ablations cause only moderate slowdowns, we find that removing all evaluated components together causes most query workloads to exceed the time limit. This result indicates that these components are complementary and have a substantially stronger combined effect. Overall, the ablated results show that our proposed components both avoid unnecessary fixed-M subproblems and reduce the cost of fixed-M subproblems that remain. Exp-6: Scalability w.r.t. data size. We study the scalability of BoxDPpS, by varying the graph size. We independently sample percentages of {20%, 40%, 60%, 80%, 100%} vertices from each vertex type and construct the corresponding induced HIN. For each percentage below 100%, we repeat the sampling using ten fixed random seeds; the 100% setting uses the complete dataset. We evaluate the

Jiadong Xie, Jiaming Yang, Kangfei Zhao, and Jeffrey Xu Yu

DBpedia

Freebase

104

Time (s)

103 102 101 100 10−1

20

40

60

80

(a) Retained vertices (%)

100

5

6

7

8

9

(b) Meta-path length 𝑘

Figure 5: Scalability of BoxDPpS w.r.t. data size and 𝑘 𝑘 = 3 query set at every graph scale. Each point in Fig. 5(a) reports the average per-query running time. Fig. 5(a) shows that the running time of BoxDPpS increases steadily with the retained data size. Freebase requires more time because its sampled query graphs contain substantially more vertices and edges than those of DBpedia. Nevertheless, BoxDPpS completes all queries at every tested scale, including the complete datasets. This result demonstrates that BoxDPpS scales to the full DBpedia and Freebase workloads considered in this experiment. Exp-7: Scalability w.r.t. meta-path length. We further evaluate BoxDPpS on DBpedia and Freebase using meta-paths with 𝑘 ∈ {5, 6, 7, 8, 9}. For each dataset and each value of 𝑘, we select structurally diverse queries using fixed criteria based on the induced graph size, the number of meta-path instances, and the estimated count-vector search space. Each point in Fig. 5(b) reports the average per-query running time. The running time is non-monotonic in 𝑘. Recall that for DPpS, longer meta-paths tend to induce smaller connected 𝑃-partite graphs because distant vertex types are more weakly related [4]; the decreasing vertex and edge counts in our workloads exhibit the same size effect. For BoxDPpS, however, each additional type also raises the dimension of the count-vector search. The frontier may consequently visit more boxes, and testing a lifted exact certificate on one box may inspect up to 2𝑘 log-box corners. The reduced graph size lowers instance-processing and fixed-M solving costs, whereas the higher-dimensional box search increases frontier and pruning costs. The interaction between these opposing factors explains the non-monotonic curves and shows that the induced query structure, rather than 𝑘 alone, determines the practical workload. Overall, BoxDPpS processes all tested meta-paths with up to nine vertex-type positions.

10

RELATED WORK

Densest Subgraph Search. Densest subgraph search is a fundamental graph mining problem that seeks a vertex-induced subgraph with maximum density, which has attracted extensive attention from the community [5, 9, 15, 17, 19–21, 25, 26, 28, 31, 33, 36]. The classical edge-density formulation admits exact max-flow-based solutions [10] and greedy peeling approximations [2]. Subsequent studies have improved scalability and broadened the objective family. For example, core-based reduction and localization techniques accelerate exact and approximate densest subgraph discovery on large graphs [8], while supermodularity, peeling, and flow have been studied in a unified framework [3]. Directed densest subgraph

search is closely related to our setting because it optimizes density over two vertex sets instead of one. Recent directed methods improve efficiency using core-style pruning, exact/approximate algorithms, and convex-programming formulations [22, 23]. Beyond edge density, higher-order density objectives use motifs or cliques to reveal stronger cohesion. The 𝑘-clique densest subgraph problem generalizes edge density to clique counts [32], and recent counting-based approaches reduce expensive clique enumeration when optimizing such objectives [38]. However, these studies are designed mainly for homogeneous, directed, or higher-order homogeneous graphs. The most relevant work is the densest 𝑃-partite subgraph search problem in HINs [4], which provides exact and approximate algorithms based on iRM-set enumeration and fixed-M min-cut subproblems. Our work keeps the same DPpS objective and exactness guarantee, but provides an efficient solution by pruning count space and reducing the cost of surviving fixed-M solves. Cohesive Subgraph Search in HINs. HINs model multi-typed objects and relations, where meta-paths are widely used to express typed semantics. Early HIN studies [27, 29] define meta-path-based similarity and relevance measures as the semantic foundations of cohesive subgraphs over HINs. Cohesive subgraph and community search in HINs adapt classic notions such as core, truss, clique, and relational constraints to typed networks. Meta-path-based community search identifies cohesive communities containing query vertices [7], and star-schema HINs further support hub-centered community models [18]. Truss-based HIN models strengthen cohesion using triangle-like typed structures [35], while relational community search handles dynamic HINs with relation constraints [14]. Motif-clique models capture richer typed patterns beyond simple paths [12]. Recent extensions incorporate additional semantics, including influence [37], structural similarity and vertex roles [34], and temporal proximity [30]. These models are effective for querydriven or constraint-based HIN mining, but they usually require query vertices, thresholds, influence or attribute constraints, temporal windows, or predefined structural roles. They also often focus on one target type or on satisfying a cohesiveness constraint. In contrast, DPpS is a global density maximization problem over all vertex types along a meta-path.

11

CONCLUSION

We study exact densest 𝑃-partite subgraph mining in HINs under the DPpS density objective. The main bottleneck is the combination of high-dimensional count-vector search and repeated fixed-M solves over many meta-path instances. We propose a box-level exact search that prunes count regions by safe upper bounds and lifted exact pruning conditions, deduplicates repeated iRM-sets by primitive count-vector keys, and uses bounded warm-up and seed scheduling to improve 𝜌 best early. We further introduce a multiplicity-aware exact solution via terminal-twin compression and projection grouping, thereby reducing the flow network while preserving the original auxiliary objective. Together, these techniques provide an exact solution that reduces redundant search without weakening the optimality guarantee. For the query metapaths completed, BoxDPpS achieves an arithmetic mean speedup of 27.04× over the state-of-the-art approach AdvExactGVIt.

Scalable Exact Densest 𝑃 -Partite Subgraph Search in Heterogeneous Information Networks

REFERENCES [1] [n. d.]. https://github.com/librahu/HIN-Datasets-for-Recommendation-andNetwork-Embedding [2] Moses Charikar. 2000. Greedy approximation algorithms for finding dense components in a graph. In Approximation Algorithms for Combinatorial Optimization, Third International Workshop, APPROX 2000, Saarbrücken, Germany, September 5-8, 2000, Proceedings (Lecture Notes in Computer Science, Vol. 1913), Klaus Jansen and Samir Khuller (Eds.). Springer, 84–95. [3] Chandra Chekuri, Kent Quanrud, and Manuel R. Torres. 2022. Densest Subgraph: Supermodularity, Iterative Peeling, and Flow. In Proceedings of the 2022 ACM-SIAM Symposium on Discrete Algorithms, SODA 2022, Virtual Conference / Alexandria, VA, USA, January 9 - 12, 2022, Joseph (Seffi) Naor and Niv Buchbinder (Eds.). SIAM, 1531–1555. [4] Lu Chen, Chengfei Liu, Rui Zhou, Kewen Liao, Jiajie Xu, and Jianxin Li. 2023. Densest Multipartite Subgraph Search in Heterogeneous Information Networks. Proc. VLDB Endow. 17, 4 (2023), 699–711. [5] Xin Chen, Wenqing Lin, Haoxuan Xie, Sibo Wang, and Siqiang Luo. 2025. Finding Near-Optimal Maximum Set of Disjoint $k$-Cliques in Real-World Social Networks. In 41st IEEE International Conference on Data Engineering, ICDE 2025, Hong Kong, May 19-23, 2025. IEEE, 3316–3328. [6] Zheng Dong, Xin Huang, Guorui Yuan, Hengshu Zhu, and Hui Xiong. 2021. Butterfly-Core Community Search over Labeled Graphs. Proc. VLDB Endow. 14, 11 (2021), 2006–2018. [7] Yixiang Fang, Yixing Yang, Wenjie Zhang, Xuemin Lin, and Xin Cao. 2020. Effective and Efficient Community Search over Large Heterogeneous Information Networks. Proc. VLDB Endow. 13, 6 (2020), 854–867. [8] Yixiang Fang, Kaiqiang Yu, Reynold Cheng, Laks V. S. Lakshmanan, and Xuemin Lin. 2019. Efficient Algorithms for Densest Subgraph Discovery. Proc. VLDB Endow. 12, 11 (2019), 1719–1732. [9] Edoardo Galimberti, Francesco Bonchi, and Francesco Gullo. 2017. Core Decomposition and Densest Subgraph in Multilayer Networks. In Proceedings of the 2017 ACM on Conference on Information and Knowledge Management, CIKM 2017, Singapore, November 06 - 10, 2017. ACM, 1807–1816. [10] A. V. Goldberg. 1984. Finding a Maximum Density Subgraph. Technical Report UCB/CSD-84-171. EECS Department, University of California, Berkeley. https: //www2.eecs.berkeley.edu/Pubs/TechRpts/1984/5956.html [11] Dorit S. Hochbaum. 2008. The Pseudoflow Algorithm: A New Algorithm for the Maximum-Flow Problem. Oper. Res. 56, 4 (2008), 992–1009. [12] Jiafeng Hu, Reynold Cheng, Kevin Chen-Chuan Chang, Aravind Sankar, Yixiang Fang, and Brian Y. H. Lam. 2019. Discovering Maximal Motif Cliques in Large Heterogeneous Information Networks. In 35th IEEE International Conference on Data Engineering, ICDE 2019, Macao, China, April 8-11, 2019. IEEE, 746–757. [13] Zhipeng Huang, Yudian Zheng, Reynold Cheng, Yizhou Sun, Nikos Mamoulis, and Xiang Li. 2016. Meta Structure: Computing Relevance in Large Heterogeneous Information Networks. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, San Francisco, CA, USA, August 13-17, 2016. ACM, 1595–1604. [14] Xun Jian, Yue Wang, and Lei Chen. 2020. Effective and Efficient Relational Community Detection and Search in Large Dynamic Heterogeneous Information Networks. Proc. VLDB Endow. 13, 10 (2020), 1723–1736. [15] Jiawei Jiang, Hao Huang, Zhigao Zheng, Yi Wei, Fangcheng Fu, Xiaosen Li, and Bin Cui. 2025. Detecting and Analyzing Motifs in Large-Scale Online Transaction Networks. IEEE Trans. Knowl. Data Eng. 37, 2 (2025), 584–596. [16] Jiaxin Jiang, Siyuan Yao, Yuhang Chen, Bingsheng He, Yudong Niu, Yuchen Li, Shixuan Sun, and Yongchao Liu. 2025. Community Detection in Heterogeneous Information Networks Without Materialization. Proc. ACM Manag. Data 3, 3 (2025), 139:1–139:27. [17] Jiaxin Jiang, Siyuan Yao, Yuchen Li, Qiange Wang, Bingsheng He, and Min Chen. 2025. Dupin: A Parallel Framework for Densest Subgraph Discovery in Fraud Detection on Massive Graphs. Proc. ACM Manag. Data 3, 3 (2025), 150:1–150:26. [18] Yangqin Jiang, Yixiang Fang, Chenhao Ma, Xin Cao, and Chunshan Li. 2022. Effective Community Search over Large Star-Schema Heterogeneous Information Networks. Proc. VLDB Endow. 15, 11 (2022), 2307–2320. [19] Song Kim, Dahee Kim, Taejoon Han, Junghoon Kim, Hyun Ji Jeong, and Jungeun Kim. 2026. Efficient Locality-based Indexing for Cohesive Subgraphs Discovery in Hypergraphs. In Proceedings 29th International Conference on Extending

Database Technology, EDBT 2026, Tampere, Finland, March 24-27, 2026. OpenProceedings.org, 170–182. [20] Xuanming Liu, Tingjian Ge, and Yinghui Wu. 2022. A Stochastic Approach to Finding Densest Temporal Subgraphs in Dynamic Graphs. IEEE Trans. Knowl. Data Eng. 34, 7 (2022), 3082–3094. [21] Chengyang Luo, Qing Liu, Yunjun Gao, and Jianliang Xu. 2025. Synergetic Community Search over Large Multilayer Graphs. Proc. VLDB Endow. 18, 5 (2025), 1412–1424. [22] Chenhao Ma, Yixiang Fang, Reynold Cheng, Laks V. S. Lakshmanan, and Xiaolin Han. 2022. A Convex-Programming Approach for Efficient Directed Densest Subgraph Discovery. In SIGMOD ’22: International Conference on Management of Data, Philadelphia, PA, USA, June 12 - 17, 2022, Zachary G. Ives, Angela Bonifati, and Amr El Abbadi (Eds.). ACM, 845–859. [23] Chenhao Ma, Yixiang Fang, Reynold Cheng, Laks V. S. Lakshmanan, Wenjie Zhang, and Xuemin Lin. 2020. Efficient Algorithms for Densest Subgraph Discovery on Large Directed Graphs. In Proceedings of the 2020 International Conference on Management of Data, SIGMOD Conference 2020, online conference [Portland, OR, USA], June 14-19, 2020, David Maier, Rachel Pottinger, AnHai Doan, Wang-Chiew Tan, Abdussalam Alawini, and Hung Q. Ngo (Eds.). ACM, 1051–1066. [24] Omid Madani, Sai Ankith Averineni, and Shashidhar Gandham. 2022. A Dataset of Networks of Computing Hosts. In IWSPA@CODASPY 2022: Proceedings of the 2022 ACM on International Workshop on Security and Privacy Analytics, Baltimore, MD, USA, April 27, 2022, Andrew H. Sung, Rakesh M. Verma, and Roland H. C. Yap (Eds.). ACM, 100–104. [25] Lutz Oettershagen, Honglian Wang, and Aristides Gionis. 2024. Finding Densest Subgraphs with Edge-Color Constraints. In Proceedings of the ACM on Web Conference 2024, WWW 2024, Singapore, May 13-17, 2024. ACM, 936–947. [26] Arkaprava Saha, Xiangyu Ke, Arijit Khan, and Cheng Long. 2023. Most Probable Densest Subgraphs. In 39th IEEE International Conference on Data Engineering, ICDE 2023, Anaheim, CA, USA, April 3-7, 2023. IEEE, 1447–1460. [27] Chuan Shi, Xiangnan Kong, Yue Huang, Philip S. Yu, and Bin Wu. 2014. HeteSim: A General Framework for Relevance Measure in Heterogeneous Networks. IEEE Trans. Knowl. Data Eng. 26, 10 (2014), 2479–2492. [28] Longxu Sun, Xin Huang, Rong-Hua Li, Byron Choi, and Jianliang Xu. 2022. Index-Based Intimate-Core Community Search in Large Weighted Graphs. IEEE Trans. Knowl. Data Eng. 34, 9 (2022), 4313–4327. [29] Yizhou Sun, Jiawei Han, Xifeng Yan, Philip S. Yu, and Tianyi Wu. 2011. PathSim: Meta Path-Based Top-K Similarity Search in Heterogeneous Information Networks. Proc. VLDB Endow. 4, 11 (2011), 992–1003. [30] Yifu Tang, Chengfei Liu, Lu Chen, Rui Zhou, and Jianxin Li. 2025. Finding TimeProximity Communities in Temporal Heterogeneous Information Networks. Proc. VLDB Endow. 18, 13 (2025), 5740–5752. [31] Siyi Teng, Jiadong Xie, Fan Zhang, Can Lu, Juntao Fang, and Kai Wang. 2024. Optimizing Network Resilience via Vertex Anchoring. In Proceedings of the ACM on Web Conference 2024, WWW 2024, Singapore, May 13-17, 2024. ACM, 606–617. [32] Charalampos E. Tsourakakis. 2015. The K-clique Densest Subgraph Problem. In Proceedings of the 24th International Conference on World Wide Web, WWW 2015, Florence, Italy, May 18-22, 2015, Aldo Gangemi, Stefano Leonardi, and Alessandro Panconesi (Eds.). ACM, 1122–1132. [33] Kaixin Wang, Kaiqiang Yu, and Cheng Long. 2024. Efficient k-Clique Listing: An Edge-Oriented Branching Strategy. Proc. ACM Manag. Data 2, 1 (2024), 7:1–7:26. [34] Shu Wang, Yixiang Fang, and Wensheng Luo. 2025. Searching and Detecting Structurally Similar Communities in Large Heterogeneous Information Networks. Proc. VLDB Endow. 18, 5 (2025), 1425–1438. [35] Yixing Yang, Yixiang Fang, Xuemin Lin, and Wenjie Zhang. 2020. Effective and Efficient Truss Computation over Large Heterogeneous Information Networks. In 36th IEEE International Conference on Data Engineering, ICDE 2020, Dallas, TX, USA, April 20-24, 2020. IEEE, 901–912. [36] Fan Zhang, Jiadong Xie, Kai Wang, Shiyu Yang, and Yu Jiang. 2022. Discovering key users for defending network structural stability. World Wide Web 25, 2 (2022), 679–701. [37] Yingli Zhou, Yixiang Fang, Wensheng Luo, and Yunming Ye. 2023. Influential Community Search over Large Heterogeneous Information Networks. Proc. VLDB Endow. 16, 8 (2023), 2047–2060. [38] Yingli Zhou, Qingshuo Guo, Yixiang Fang, and Chenhao Ma. 2024. A Countingbased Approach for Efficient k-Clique Densest Subgraph Discovery. Proc. ACM Manag. Data 2, 3 (2024), 119.

Related documents

• Algorithmic Reduction of Biological Networks With Multiple Time Scales#617719
• Algorithmic Reduction of Biological Networks With Multiple Time Scales#617219
• Graphs as focal interaction systems in higher-order dynamics#1076043
• From coincidence to purposeful flow? properties of transcendental information cascades#749542
• ConRAG: Lightweight inference of multi-hop relations#1129790