Algebraic Subgraph Counting Qiuyu Guo
Jianye Yang∗
Wenjie Zhang
University of New South Wales Sydney, Australia [email protected]
Guangzhou University PengCheng Laboratory Guangzhou, China [email protected]
University of New South Wales Sydney, Australia [email protected]
Hanchen Wang
Ying Zhang
Xuemin Lin
University of Technology Sydney Sydney, Australia [email protected]
Zhejiang Gongshang University Hangzhou, China [email protected]
The Chinese University of Hong Kong, Shenzhen Shenzhen, China [email protected]
Abstract Subgraph isomorphism counting is a fundamental problem in graph analytics, which aims to find the number of subgraph isomorphisms of a query graph in a data graph. The candidate tree-based framework provides a promising foundation for subgraph counting tasks, offering a unified counting paradigm that can be extended beyond tree patterns. However, supporting subgraph isomorphism within this framework remains challenging, as it requires handling both the non-tree edge constraint and the injective mapping constraint. Although existing solutions employ sampling or learning techniques to address these constraints in this framework, they still either suffer from inherent sampling failures or rely heavily on supervision. In this paper, we propose ASC, an algebraic subgraph counting approach built on the candidate tree-based counting framework. In our method, the non-tree edge constraint is directly incorporated into the candidate tree-based counting process through a matrix-based computation method, enabling subgraph homomorphism counting with high accuracy in polynomial time. Based on the resulting subgraph homomorphism count, we further apply a local sampling method to address the injective mapping constraint, thereby obtaining the final subgraph isomorphism count. Extensive experiments show that ASC can achieve substantially better and more stable performance over the baselines across various datasets, while scaling to billion-edge graphs. Most impressively, as a non-learning method, ASC can even achieve more than an order of magnitude higher average accuracy than the state-of-the-art learning-based method FlowSC with similar efficiency. This paper is the full version of the work accepted at SIGMOD 2027. The code is available at https://github.com/EricaGuoQiuyu/AlgebraicSubgraphCounting.
1
Introduction
Given a query graph 𝑄 and a data graph 𝐺, the problem of subgraph counting is to find the number of subgraph isomorphisms of 𝑄 in 𝐺. It is a fundamental problem with a wide range of applications, such as query optimization [4, 8, 17, 21], biological network analysis [12, 41, 42, 54], and fraud detection [44, 45]. Despite its importance, subgraph counting is computationally challenging due to its #phardness [18, 47]. Exact counting via subgraph matching [5, 14, 27, 28, 35, 36, 51, 58] is often intractable due to the inherent complexity of enumeration. Even counting a 32-clique query in a 32-clique ∗ Corresponding author
graph can exceed a two-hour time limit in our experiments. To broaden applicability, recent efforts have focused on approximate counting [18, 30, 39, 47]. In this paper, we focus on approximate subgraph counting. Existing Solutions and Limitations. Current solutions for approximate subgraph counting can generally be classified into three categories, including summarization-based, sampling-based, and learning-based methods. Summarization-based methods [11, 16, 39, 48] adopt a decomposition-aggregation paradigm, decomposing the query graph into substructures and aggregating their individual counts to estimate the total count. However, their reliance on the independence assumption between substructures often leads to unsatisfactory accuracy in real-world graphs. Sampling-based methods [10, 13, 30, 47] estimate counts via match frequencies in sampled subgraphs. While sampling can achieve high accuracy for simple queries on small graphs, they suffer from inherent sampling failures as the sample space expands. To address this, Fastest [47] employs advanced filtering and candidate tree-based sampling to prune the search space and improve sampling success rates. Nonetheless, eliminating sampling failures remains challenging for complex queries on large graphs. Learning-based methods [18, 24, 55, 60] predict counts by regressing from learned graph features. Earlier works [55, 60] struggle to capture the relationship between the graph features and counts or face training instabilities. The state-of-the-art FlowSC [18] mitigates these issues by simulating data flow within a candidate-tree framework and using a pretraining scheme. However, the performance of these techniques is heavily reliant on the labeled data, which is often impractical for hard query workloads where obtaining labels itself is computationally prohibitive. Matrix-based Pattern Counting. The connection between matrix computation and pattern counting has deep roots. Early seminal works [3, 25, 33] exploit fast matrix multiplication to exactly count specific structured patterns like small cycles and cliques. The Functional Aggregate Query (FAQ) framework [2] provides a unified theoretical foundation for these concepts, formulating general subgraph counting as a sum-product query:
∑︂
∏︂
𝑥 1 ,...,𝑥𝑘 ∈𝑉𝐺 (𝑖,𝑗 ) ∈𝐸𝑄
𝐸 (𝑥𝑖 , 𝑥 𝑗 ),
(1)
Guo et al.
where each variable 𝑥𝑖 denotes the mapping of a query vertex to a data vertex, and each factor 𝐸 (𝑥𝑖 , 𝑥 𝑗 ) enforces edge consistency. Under this framework, matrix multiplication can be seen as a special case of FAQ evaluation, which establishes a natural connection between subgraph counting and matrix-based computation. To solve such queries generally, FAQ proposes InsideOut, an exact algorithm based on variable elimination. Crucially, the FAQ perspective also reveals the fundamental bottleneck shared by all these exact methods: their efficiency is strictly bounded by 𝑂 (𝑁 faqw ), where 𝑁 is the data size and faqw is the fractional FAQ-width. While early matrix approaches and InsideOut are effective for low-width topologies, faqw is large for arbitrary or randomly generated subgraphs. Without an exploitable structure, exact evaluation becomes computationally impractical in general, motivating scalable alternatives beyond exact evaluation. Candidate Tree-based Counting. The candidate tree-based framework supports polynomial-time subgraph homomorphism counting against tree queries [18, 47]. However, there are two major issues in extending it for subgraph isomorphism counting against generic query graphs. First, it ignores the non-tree edges, i.e., edges exist in the query graph but not in its spanning tree. Ignoring the matching constraints defined by such edges would lead to significant overestimation. Second, it cannot handle the injective mapping constraint required by isomorphism. Because the candidate tree-based counting only maintains the weight of matching rather than matching states, it cannot enforce injectivity and thus overestimates the counts, as a data vertex might be used multiple times in a matching. To apply the candidate tree-based counting framework to subgraph isomorphism counting, existing methods develop effective techniques to handle the non-tree edge and injective constraints. Fastest treats the tree homomorphism count as an upper bound and adopts sampling to approximate both the non-tree and injectivity constraints. FlowSC simulates the data flow pattern in the candidate tree counting process and learns the relationship between graph features and subgraph counts, leveraging query features and supervision to implicitly enforce the constraints. However, these methods often suffer from sampling failures in excessive sampling spaces or from insufficient supervision. Our Approach and Contributions. In this paper, we propose a candidate tree-based counting method, called Algebraic Subgraph Counting (ASC). Since the candidate tree framework inherently satisfies tree-edge constraints, we decouple the remaining requirements for subgraph isomorphism into non-tree edge constraints and isomorphism (injectivity mapping) constraints. To handle the non-tree edges, we divide them into two categories: those connecting vertices with the same parent (Type I) and those connecting vertices with different parents (Type II). We first show that it takes exponential time to compute the exact homomorphism count if the Type-I non-tree edges are preserved in the candidate tree framework. We propose a matrix-based method to incorporate Type-I non-tree edge constraints into candidate tree-based counting, delivering high-quality approximations in polynomial time. For the Type-II non-tree edge constraint, we have a key observation that the compatibility of the candidates of the two endpoints of a Type-II non-tree edge with a valid matching depends on their lowest common ancestor (LCA) in the tree. We incorporate the Type-II
non-tree edge constraints into the candidate tree-based counting process by computing the LCA ratios, which are used to discard the invalid weights accumulated at LCAs in the counting phase. Together with the Type-I non-tree edge constraint computation, this leads to an algebraic method for subgraph homomorphism counting. Based on the resulting homomorphism count, we further apply local sampling to enforce injective mapping constraints, producing the final subgraph isomorphism count. We conduct extensive experiments on 10 real-world datasets, including two billion-edge graphs. The results demonstrate that ASC can achieve consistently better and more stable performance over the existing solutions under the majority of query settings.
2 Preliminary 2.1 Problem Definition In this paper, we focus on a simple vertex-labeled, undirected, and connected graph 𝑔 = (𝑉 , 𝐸, 𝐿, Σ). Here, 𝑉 is a set of vertices of 𝑔, 𝐸 ⊆ 𝑉 × 𝑉 is a set of edges of 𝑔, 𝐿 is a labeling function that assigns a label 𝐿(𝑣) ∈ Σ to each vertex 𝑣 in 𝑔, and Σ is a set of labels. We use 𝑁 (𝑣) = {𝑣 ′ ∈ 𝑉 | (𝑣, 𝑣 ′ ) ∈ 𝐸} to denote the neighbors of a vertex 𝑣, and 𝑑 (𝑣) = |𝑁 (𝑣)| to denote the degree of 𝑣. In this paper, we use 𝑄 = (𝑉𝑄 , 𝐸𝑄 , 𝐿, Σ) and 𝐺 = (𝑉𝐺 , 𝐸𝐺 , 𝐿, Σ) to denote the query graph and the data graph, respectively. Besides, we use 𝑢 to represent a vertex in 𝑄 and 𝑣 to represent a vertex in 𝐺. Definition 2.1 (Subgraph Isomorphism and Homomorphism). Given a query graph 𝑄 and a data graph 𝐺, subgraph isomorphism is an injective mapping 𝑓 : 𝑉𝑄 ↦→ 𝑉𝐺 , such that ∀𝑢 ∈ 𝑉𝑄 , 𝐿(𝑢) = 𝐿(𝑓 (𝑢)) and ∀(𝑢, 𝑢 ′ ) ∈ 𝐸𝑄 , (𝑓 (𝑢), 𝑓 (𝑢 ′ )) ∈ 𝐸𝐺 . When injectivity is not required, the mapping 𝑓 defines a subgraph homomorphism. In this paper, we focus on the subgraph isomorphism. An injective mapping from 𝑄 to 𝐺 is called a subgraph isomorphic embedding, also referred to as an embedding or a match when context is clear. Definition 2.2 (Subgraph Matching and Counting). Given a query graph 𝑄 and a data graph 𝐺, subgraph matching is to find all embeddings of 𝑄 in 𝐺, while subgraph counting is to compute the number of embeddings. Given the #P-hardness nature of exact subgraph counting, we focus on approximate solutions for subgraph counting in this paper. Problem Statement. Given a query graph 𝑄 and a data graph 𝐺, we aim to approximate the number of all embeddings of 𝑄 in 𝐺.
2.2
Candidate Tree-based Counting
Definition 2.3 (Candidate Space). Given a query graph 𝑄 and a data graph 𝐺, a candidate space 𝐶𝑆 maintains a candidate vertex (Resp. edge) set for each vertex (Resp. edge) in 𝑄, such that all matches of 𝑄 in 𝐺 are preserved in the solution space defined by 𝐶𝑆. • Given a vertex 𝑢 ∈ 𝑉𝑄 , its candidate vertex set 𝐶 (𝑢) ⊆ 𝑉𝐺 satisfies that, for any embedding 𝑓 , 𝑓 (𝑢) ∈ 𝐶 (𝑢) must hold, but a vertex in 𝐶 (𝑢) is not necessarily mapped to 𝑢 in any valid embedding. • Given an edge (𝑢, 𝑢 ′ ) ∈ 𝐸𝑄 , its candidate edge set 𝐶 (𝑢, 𝑢 ′ ) ⊆ 𝐸𝐺 satisfies that, for any embedding 𝑓 , (𝑓 (𝑢), 𝑓 (𝑢 ′ )) ∈ 𝐶 (𝑢, 𝑢 ′ ) must hold, but an edge in 𝐶 (𝑢, 𝑢 ′ ) is not necessarily mapped to (𝑢, 𝑢 ′ ) in any valid embedding.
Algebraic Subgraph Counting
u0 A u1 B
v0 A C u2 D u3
u4 A
v1 B
B v2 v3 C
v0 A C v4
D v5
u0 A
v0 A A v6
(b) Data graph G
(c) Candidate space CS
v0 A
W(u4, v0) = 1
u1 B
C u2 D u3 u4 A
(d) Query tree TQ
v1 B B v2 v3 C C v4
D v5
v3 C C v4
A v6
(a) Query graph Q
v1 B B v2
D v5
v0 A A v6 (e) Candidate tree TC
W(u4, v6) = 1
W(u3, v5) = 1 W(u2, v3) = 1
W(u2, v4) = 1
W(u1, v1) = 2
W(u1, v2) = 2
W(u0, v0) = 8
(f) Candidate tree counting
weights of all root matching pairs are collected to obtain the homomorphism count for the 𝑇𝑄 . More specifically, for a vertex 𝑢 ∈ 𝑉𝑄 , we use 𝑁𝑐 (𝑢) to denote its children in 𝑇𝑄 . For each 𝑣 ∈ 𝐶 (𝑢), we use a weight 𝑊 (𝑢, 𝑣) to represent the number of homomorphisms of the subtree 𝑇𝑄 (𝑢) (rooted at 𝑢) contained in the candidate subtree 𝑇𝐶 (𝑣) (rooted at 𝑣). Then, 𝑊 (𝑢, 𝑣) is updated as follows: ∏︂ ∑︂ 𝑊 (𝑢, 𝑣) = 𝑊 (𝑢𝑐 , 𝑣𝑐 ). (2) 𝑢𝑐 ∈𝑁𝑐 (𝑢 ) 𝑣𝑐 ∈𝐶 (𝑢𝑐 |𝑢,𝑣)
Each matching pair updates its weight with the weights of its child matching pairs, and 𝑊 (𝑢, 𝑣) = 1 if 𝑢 is a leaf vertex. The counting process proceeds from the leaves to the root, with an overall 𝑂 (|𝐸𝑄 ||𝐸𝐶 |) time [47]. v0 A
• Given an edge (𝑢, 𝑢 ′ ) ∈ 𝐸𝑄 and a candidate vertex 𝑣 ∈ 𝐶 (𝑢), the conditional candidate vertex set 𝐶 (𝑢 ′ | 𝑢, 𝑣) ⊆ 𝐶 (𝑢 ′ ) satisfies that, for any embedding 𝑓 with 𝑓 (𝑢) = 𝑣, 𝑓 (𝑢 ′ ) ∈ 𝐶 (𝑢 ′ | 𝑢, 𝑣) must hold, but a vertex in 𝐶 (𝑢 ′ | 𝑢, 𝑣) is not necessarily mapped to 𝑢 ′ in any valid embedding.
Example 2.7. Continuing the Example 2.6, we compute the candidate tree count inv1 Figure andC the weights for each matching v4 updated B B v1(e) D v5 2 v3 C pair are shown in Figure 1(f). Candidate tree counting proceeds in a bottom-up manner, from the leaves. All leaf matching pairs v0 A starting A v6 have a weight of 1, e.g., 𝑊 (𝑢 4, 𝑣 0 ) = 1, and 𝑊 (𝑢 4, 𝑣 6 ) = 1. Following the counting rule in Equation 2, 𝑊 (𝑢 1, 𝑣 1 ) = (𝑊 (𝑢 4, 𝑣 0 ) +𝑊 (𝑢 4, 𝑣 6 )), so we have 𝑊 (𝑢 1, 𝑣 1 ) = 2. Similarly, 𝑊 (𝑢 1, 𝑣 2 ) = 2. At the root level, 𝑊 (𝑢 0, 𝑣 0 ) = (𝑊 (𝑢 1, 𝑣 1 ) + 𝑊 (𝑢 1, 𝑣 2 )) × (𝑊 (𝑢 2, 𝑣 3 ) + 𝑊 (𝑢 2, 𝑣 4 )) × 𝑊 (𝑢 3, 𝑣 5 ) = 8, so the tree homomorphism count is 8.
Example 2.4. Consider the query graph 𝑄 and data graph 𝐺 shown in Figure 1(a) and (b). We construct the corresponding candidate space 𝐶𝑆 shown in Figure 1(c). In this candidate space, 𝑣 1 and 𝑣 2 belong to 𝐶 (𝑢 1 ), (𝑣 0, 𝑣 1 ) ∈ 𝐶 (𝑢 0, 𝑢 1 ), and 𝑣 1 ∈ 𝐶 (𝑢 1 | 𝑢 2, 𝑣 3 ) as 𝑢 1 can match 𝑣 1 under the condition 𝑢 2 ↦→ 𝑣 3 .
Issues for Subgraph Isomorphism. The candidate tree-based counting method offers an efficient paradigm for subgraph homomorphism counting against tree queries. However, there are two major issues in extending it for subgraph isomorphism counting against generic query graphs.
To improve the efficiency and accuracy of subgraph counting, existing methods often build a compact candidate space by removing the invalid vertices and edges through candidate filtering algorithms [7, 18, 22, 47]. Then, the subsequent counting process is performed against a more compact candidate space instead of the original data graph, without sacrificing completeness. The candidate tree structure organizes the candidate space hierarchically to assist subgraph counting.
• Issue 1: Non-tree Edge Constraint. As shown in Equation 2, for each matching pair (𝑢, 𝑣), we compute, for each child vertex 𝑢𝑐 ∈ 𝑁𝑐 (𝑢), the sum of 𝑊 (𝑢𝑐 , 𝑣𝑐 ) over all 𝑣𝑐 ∈ 𝐶 (𝑢𝑐 | 𝑢, 𝑣), and then take the product of these per-child sums to obtain 𝑊 (𝑢, 𝑣). However, in the processing, the non-tree edges are ignored, which may significantly affect the accumulated weight for the parent vertex, since every edge in the query graph must be matched. Here, a non-tree edge refers to an edge existing in the query graph but not in the query tree.
Figure 1: Running example of candidate tree-based counting.
Definition 2.5 (Candidate Tree). Let 𝑇𝑄 be a spanning tree of the query graph 𝑄. The candidate tree 𝑇𝐶 is obtained by organizing the candidate space according to the structure of 𝑇𝑄 : each query vertex 𝑢 in 𝑇𝑄 is associated with its candidate set 𝐶 (𝑢), and a tree edge (𝑢, 𝑢 ′ ) in 𝑇𝑄 is replaced by the candidate edges between 𝐶 (𝑢) and 𝐶 (𝑢 ′ ). Consequently, 𝑇𝐶 encodes all tree homomorphisms of 𝑇𝑄 in the candidate space. Example 2.6. For the query graph 𝑄 in Figure 1(a), we construct its BFS spanning tree 𝑇𝑄 as shown in Figure 1(d). Its corresponding candidate tree𝑇𝐶 is illustrated in Figure 1(e). In particular, (⟨𝑢 0, 𝑣 0 ⟩, ⟨𝑢 1, 𝑣 1 ⟩, ⟨𝑢 2, 𝑣 3 ⟩, ⟨𝑢 3, 𝑣 5 ⟩, ⟨𝑢 4, 𝑣 0 ⟩) forms a homomorphism embedding of 𝑇𝑄 . General Idea of Candidate Tree-based Counting. The candidate tree-based counting framework utilizes dynamic programming to compute the number of homomorphisms of 𝑇𝑄 in polynomial time, providing a tight upper bound for subgraph isomorphism counting of 𝑄. In this framework, homomorphism counts are accumulated in a bottom-up manner along the tree edges of 𝑇𝐶 . During the processing, each matching pair maintains a weight as the number of homomorphism matches of its subtrees. At the root level, the
Example 2.8. Comparing the 𝑄 in Figure 1(a) with its BFS tree 𝑇𝑄 in Figure 1(d), there are 2 non-tree edges (𝑢 1, 𝑢 2 ) and (𝑢 2, 𝑢 4 ), that are not included in 𝑇𝑄 . As in Example 2.7, the resulting candidate-tree count is 8. However, when graph homomorphism constraints are considered, the true count is 4. The discrepancy is caused by non-tree constraints. Considering the constraints imposed by the non-tree edge (𝑢 1, 𝑢 2 ), the mappings 𝑢 1 ↦→ 𝑣 2 and 𝑢 2 ↦→ 𝑣 3 are conflicting. Nevertheless, the weight update of 𝑊 (𝑢 0, 𝑣 0 ) implicitly assumes that every candidate in 𝐶 (𝑢 1 | 𝑢 0, 𝑣 0 ) connects to every candidate in 𝐶 (𝑢 2 | 𝑢 0, 𝑣 0 ), which does not hold in practice. Similarly, the constraint imposed by (𝑢 2, 𝑢 4 ) further restricts valid combinations and should also limit the weights propagation in the candidate tree counting. • Issue 2: Isomorphism Constraint. In the processing of candidate tree-based counting, only the weights of matching pairs are maintained, rather than the explicit matching states. This implies that a data vertex may correspond to multiple query vertices, resulting in a non-injective mapping, which overestimates the number of subgraph isomorphism embeddings. We use the term isomorphism (iso) constraint to refer to the injectivity mapping constraint.
Guo et al.
Example 2.9. Continuing the Example 2.8, when vertex labels and vertex degrees are considered as filtering criteria [7, 53], 𝑣 0 is contained in both 𝐶 (𝑢 0 ) and 𝐶 (𝑢 4 ) as shown in Figure 1(c). Consequently, in the current candidate tree, the weight of leaf 𝑣 0 in 𝑊 (𝑢 4, 𝑣 0 ) is propagated upward to the root 𝑣 0 in 𝑊 (𝑢 0, 𝑣 0 ), causing 𝑣 0 to be counted twice and resulting in an overestimation in subgraph isomorphism count. Remedy of Existing Solutions. To support generic subgraph counting, prior works extend candidate tree counting by approximating non-tree and injective constraints. Fastest [47] estimates the final count by sampling candidate trees and evaluating whether they form valid query matches. However, its performance deteriorates rapidly when the sampling space becomes prohibitively large, which often happens on larger datasets and leads to sampling failures. Alternatively, FlowSC [18] employs supervised learning with cross-graph attention to implicitly emulate matching checks. However, its effectiveness relies on sufficient ground truth, which is costly to obtain for hard query workloads.
Algorithm 1: Algebraic Subgraph Counting Input: A query graph 𝑄, a data graph 𝐺 Output: Estimated subgraph count of 𝑄 in 𝐺 1 CS ← BipartitePlus (𝑄, 𝐺); 2 𝑇𝑄 ← BuildTree(𝑄); 3 LCARatio(𝐶𝑆, 𝑇𝑄 ); 4 foreach 𝑢 in reversed BFS Traversal of 𝑇𝑄 do 5 𝑠 ← BuildSequence(𝑄 [𝑁𝑐 (𝑢)]); 6 foreach 𝑣 ∈ 𝐶 (𝑢) do 7 𝑊 (𝑢, 𝑣) ← SequenceCounting(𝑢, 𝑣,𝐶𝑆,𝑠); 8 if 𝑢 is a LCA vertex then 9 𝑊 (𝑢, 𝑣) ← 𝑊 (𝑢, 𝑣) × 𝐿𝐶𝐴(𝑢, 𝑣); HomoCounts ← 0; foreach 𝑣𝑟 ∈ 𝐶 (𝑢𝑟 ) where 𝑢𝑟 is the root of 𝑇𝑄 do 12 HomoCounts + = 𝑊 (𝑢𝑟 , 𝑣𝑟 );
10 11
IsoRatio ← LocalSampling(𝐶𝑆, 𝑄); IsoCounts ← HomoCounts × IsoRatio; 15 return IsoCounts 13 14
3
Our Approach
Overview. Built upon the candidate tree framework, ASC first incorporates Type-I constraints via a matrix-based scheme (section 4), and utilizes an LCA ratio mechanism to filter out invalid combinatorial weights induced by Type-II constraints (section 5). Based on the resulting algebraic homomorphism count, we finally enforce the injectivity constraints via local sampling on subgraphs induced by same-label vertices (section 6). Algorithm Details. Algorithm 1 summarizes the counting process of ASC. Given a query graph 𝑄 and a data graph 𝐺, we first construct and refine the candidate space 𝐶𝑆 using the filtering method BipartitePlus [18] (line 1), and build a BFS spanning tree 𝑇𝑄 of 𝑄. We choose BFS rather than DFS spanning trees because the former generates more Type-I non-tree edges, whereas the latter generates only Type-II non-tree edges. Type-I non-tree edges are easier to handle than Type-II ones, as discussed in the following two sections. Then we compute the LCA ratio for each LCA candidate (line 3). In the counting phase, the parent weight 𝑊 (𝑢, 𝑣) is updated from its child vertices 𝑁𝑐 (𝑢) using a matrix computation-based counting method SequenceCounting that incorporates Type-I non-tree constraints in 𝑄 [𝑁𝑐 (𝑢)] (line 7), where 𝑄 [𝑁𝑐 (𝑢)] denotes the subgraph of 𝑄 induced by 𝑁𝑐 (𝑢). The Type-II non-tree edge constraints are handled by applying the LCA ratios at LCA vertices (line 9). We update weights in reversed BFS order and aggregate them at the root (line 12) to obtain the subgraph homomorphism count1 . An isomorphism ratio is then estimated using a local sampling method (line 13), and the final subgraph isomorphism count is obtained by multiplying the homomorphism count by this ratio (line 14).
4
Handling Type-I Non-tree Edge Constraints
In this section, we introduce the key component of ASC, which is a matrix-based counting model that can approximate the Type-I non-tree edge constraint in candidate tree-based counting. The tree counting formulation in Equation 2 implicitly assumes child vertices are mutually independent, leading to error accumulation when inter-dependencies (i.e., non-tree edge constraints) exist. 1 As BipartitePlus applies injectivity filtering, this is not a strict homomorphism count.
While exact counting for these constraints among child vertices typically incurs exponential complexity, we develop an algebraic method to simulate this process in polynomial time. By transforming candidate connectivity into transition matrices and candidate weights into weight vectors, we propagate weights via matrix multiplications to simulate weight accumulation in the Cartesian product space (subsection 4.1). For complex structures like cycles among the children, we utilize weight matrices to encode both single and multiple indirect constraints (subsection 4.2). By traversing the children through carefully designed matrix operations, our method incorporates Type-I non-tree constraints into the tree counting, with the resulting vector sum yielding the weight for a parent vertex.
4.1
A Simplified Counting Case
Fixing Candidate Tree-based Counting. Equation 2 formulates a tree-based counting paradigm where weight updates are computed by decoupling each child branch, implicitly assuming that child vertices are mutually isolated. To accommodate it to the case where the child vertices are connected (i.e., non-tree edges), we need to inspect each matching combination for the child vertices. Specifically, given a query vertex 𝑢 and its matching vertex 𝑣, let Ω = ×𝑢𝑐 ∈𝑁𝑐 (𝑢 ) 𝐶 (𝑢𝑐 | 𝑢, 𝑣) be the Cartesian product of the conditional candidates of 𝑢’s child vertices. Clearly, Ω defines a complete matching space for the child vertices of 𝑢 when 𝑢 ↦→ 𝑣. Based on Ω, we can compute the weight as follows: 𝑊 (𝑢, 𝑣) =
∑︂ ⎛ ∏︂ ⎞ 𝑊 (𝑢𝑐 , v[𝑢𝑐 ]) ⎟ . ⎜𝜎 (v) · v∈Ω ⎝ 𝑢𝑐 ∈𝑁𝑐 (𝑢 ) ⎠
(3)
Here, each combination v ∈ Ω represents a possible matching for vertices in 𝑁𝑐 (𝑢), and v[𝑢𝑐 ] denotes the candidate vertex mapped to 𝑢𝑐 . The term 𝜎 is an indicator function with value 1 if the subgraph induced by v is a homomorphism of 𝑄 [𝑁𝑐 (𝑢)], and 0 otherwise.
Algebraic Subgraph Counting
Equation 3 follows the same sum-product structure as the FAQ based subgraph counting (Equation 1), where the summation enumerates possible vertex mappings and the product enforces structural consistency. Compared with the standard FAQ expression, Equation 3 further incorporates candidate-tree-based weights, transforming the generic sum-product query into a recursive weighted aggregation process. Direct evaluation of Equation 3 requires enumerating all combinations in the matching space Ω, resulting in 𝑂 (𝑛 |𝑁𝑐 (𝑢 ) | ) time, which is exponential in |𝑁𝑐 (𝑢)|, assuming that the maximum size of 𝐶 (𝑢𝑐 | 𝑢, 𝑣) is 𝑛. In the FAQ framework, the InsideOut algorithm reduces the cost through optimized variable elimination, but its complexity still depends exponentially on the fractional FAQwidth of the subquery. To address this challenge, we propose a matrix-based approximation that avoids exhaustive enumeration and produces high-quality counting results in polynomial time. Intuition of Matrix Computation. The main innovation of our method is transforming the subgraph homomorphism counting into a sequence of matrix multiplication operations. Intuitively, adjacency matrix multiplication provides an efficient mechanism for evaluating connectivity between vertex sets. Moreover, its computational rule naturally aligns with Cartesian product operations across vertex sets, where each resulting entry counts the number of connecting paths between the corresponding pair of vertices in the candidate space. With a carefully designed computation strategy, we can incorporate non-tree constraints via algebraic methods. Constrained Sequence. Given a query vertex 𝑢, its child vertices can induce a set of connected components. For example, in the 𝑄 of Figure 1(a), {𝑢 1, 𝑢 2 } and {𝑢 3 } are two connected components. According to the above discussion, each connected component can be processed independently. In the following, we focus on the nontrivial connected components (i.e., containing at least 2 vertices), since the non-tree edges only exist in such components. Definition 4.1 (Constrained Sequence). Given a query vertex 𝑢 in the query tree 𝑇𝑄 , let 𝑆 be a non-trivial connected component induced by a set of vertices in 𝑁𝑐 (𝑢) from 𝑄. A constrained sequence 𝑠 is an ordered arrangement of vertices in 𝑆 with the following properties. • Given a vertex 𝑢𝑐 ∈ 𝑆, a vertex appearing before (Resp., after) 𝑢𝑐 in 𝑠 is called the predecessor (Resp., successor) of 𝑢𝑐 . • Given two vertices 𝑢𝑐′ , 𝑢𝑐 ∈ 𝑆, 𝑢𝑐′ is called the direct predecessor (Resp., direct successor) of 𝑢𝑐 if it appears immediately before (Resp., after) 𝑢𝑐 in 𝑠; Otherwise, 𝑢𝑐′ is an indirect predecessor (Resp., indirect successor) of 𝑢𝑐 if 𝑢𝑐′ is a neighbor of 𝑢𝑐 . Based on the second property, we note that there must be a non-tree edge between the indirect predecessor-successor vertices, while it is not necessary between the direct ones.
Equation 3, we maintain a weight vector for each child vertex 𝑢𝑐 of 𝑢 and a transition matrix for each direct predecessor-successor vertex pair of 𝑢, which are defined as follows. Definition 4.3 (Weight Vector). Given a vertex 𝑢𝑐 in the constrained sequence 𝑠, the weight vector of 𝑢𝑐 , denoted as V, is a column vector with |𝐶 (𝑢𝑐 | 𝑢, 𝑣)| entries, which are initialized as 𝑊 (𝑢𝑐 , 𝑣𝑐 ) for each 𝑣𝑐 ∈ 𝐶 (𝑢𝑐 | 𝑢, 𝑣). Definition 4.4 (Transition Matrix). Given a direct predecessorsuccessor vertex pair 𝑢𝑐′ and 𝑢𝑐 in the constrained sequence 𝑠, the transition matrix between 𝑢𝑐′ and 𝑢𝑐 , denoted as A, is constructed with rows and columns corresponding to 𝐶 (𝑢𝑐′ | 𝑢, 𝑣) and 𝐶 (𝑢𝑐 | 𝑢, 𝑣), respectively. This matrix A encodes the non-tree constraint between 𝑢𝑐′ and 𝑢𝑐 within the candidate space. If a non-tree edge exists between 𝑢𝑐′ and 𝑢𝑐 , then A𝑖 𝑗 = 1 if there exists a candidate edge between the 𝑖-th candidate of 𝑢𝑐′ and the 𝑗-th candidate of 𝑢𝑐 , and A𝑖 𝑗 = 0 otherwise. If no non-tree edge exists, A contains all ones, representing that all candidate combinations between 𝑢𝑐′ and 𝑢𝑐 are considered valid. Example 4.5. Consider the query tree 𝑇𝑄 in Figure 2(a) and the candidate trees 𝑇𝐶 in Figure 2(b), where ⟨𝑢 1, 𝑢 2, 𝑢 3, 𝑢 4 ⟩ is a constrained sequence. Note here that we ignore the red edges at this moment. The weight vectors of vertices are shown in Figure 2(c). Taking 𝑢 2 for example, since it is a leaf vertex and has two candidates, i.e., 𝑣 2 and 𝑣 3 , its weight vector is initialized as [[1], [1]]. The transition matrices are shown in Figure 2(d). With the weight vector and transition matrix, we can compute 𝑊 (𝑢, 𝑣) as follows. We start from the last vertex in the constrained sequence and propagate backwards. At each step, the weight vector of the successor vertex is multiplied by the transition matrix and then combined with the predecessor’s own weight vector via a Hadamard product. Therefore, the counts are propagated along the non-tree candidate edges and accumulated into the predecessor’s weights. We repeat this processing until the head is reached, where the sum of the final vector’s entries is used to update 𝑊 (𝑢, 𝑣). Example 4.6. The constrained sequence counting starts from the end of sequence, proceeds backward along the sequence, as shown in 1 the weight vector V4 of 𝐶 (𝑢 4 | 𝑢 0, 𝑣 0 ) is multiFigure 2(e): in Step ○, plied by A34 to produce the input vector V𝑖𝑛 3 for 𝐶 (𝑢 3 | 𝑢 0 , 𝑣 0 ), which are then multiplied element-wise with the original weight vector V3 2 to obtain the updated vector V𝑛𝑒𝑤 in Step ○ 3 . This process continues 3 ○. 6 Then, at backwards along the sequence, following the Steps ○the head 𝑢 1 , the sum of the entries in V𝑛𝑒𝑤 yields the count of the con1 strained sequence, which is 4 and exactly equals the homomorphism count of the constrained sequence.
Example 4.2. Consider the 𝑇𝑄 1 in Figure 2(a), where ⟨𝑢 1, 𝑢 2, 𝑢 3, 𝑢 4 ⟩ is a constrained sequence generated by vertices in 𝑁𝑐 (𝑢 0 ). Clearly, 𝑢 1 is a direct predecessor of 𝑢 2 and an indirect predecessor of 𝑢 4 .
Lemma 4.7. Given a constrained sequence 𝑠 with parent vertex 𝑢 matched to a candidate 𝑣, the matrix-based method returns the exact homomorphism count of the induced subgraph 𝑄 [𝑠] in the candidate space conditioned on 𝑢 ↦→ 𝑣, assuming there are no indirect predecessor-successor connections and no false positive candidates.
Counting via Matrix Multiplication Propagation. For ease of exposition, we start by considering a simplified case where 𝑠 has no indirect predecessor-successor vertices. Given a query tree 𝑇𝑄 and its corresponding candidate trees 𝑇𝐶 , let 𝑢 and 𝑣 be a pair of matching vertices in the trees. To compute the weight 𝑊 (𝑢, 𝑣) in
Proof. For any predecessor-successor pair (𝑢𝑐′ , 𝑢𝑐 ) in the sequence 𝑠, the matrix-vector product V𝑢𝑖𝑛′ = A𝑢𝑐′ 𝑢𝑐 × V𝑢𝑐 exclusively 𝑐 collects weights from candidates in 𝐶 (𝑢𝑐 | 𝑢, 𝑣) that form valid ′ matching with candidates in 𝐶 (𝑢𝑐 | 𝑢, 𝑣) under the non-tree constraint of (𝑢𝑐′ , 𝑢𝑐 ). The Hadamard product V𝑢𝑛𝑒𝑤 = V𝑢𝑐′ ⊙V𝑢𝑖𝑛′ updates ′ 𝑐
𝑐
Guo et al.
0# :
u4
0% :
(
&$ 1 u u u B 1 C&% 21 D 3 0& : && 1
&" 1
1 × 10 : = &2! 1 1# u4 E 5%& 00& : &0" '( %1 1
%
&#
E
v6
&%
5%& : &" 1
B &&
C
C&$ 1 D
v5
E
E
2 v62 0'( $
1
(d) Transition(b) matrices (!"tree ) !" (!" ) Candidate #
1
%&'((1 1 1) ×× 1! 1 0 1 == 11 11 11 " ! 0 1 " JKLM(05 %) %&
GG'( %&
GG'( P()* % %
②①
()* ()* '('( JKLM(0 0 0P()* !!"#$ P $ $ ) 0G $G$ $$
'( &0 % %1 && 1
0& :
⨀
&# &$ 1 = 2 5 1#$ : &! 21 1
=
21 1 1 × 1 0()* %5%&
0%
1 × 2 5$% : &# 1 2 &$ 1
21 1
=
×
50$%'( %
0&
②
&"
1
4 2 2
⨀
1 1
2 2 2
0()* %
0'( $0'( %
⨀
24 2
=
⨀
1
&!
=
②
1 = 14 1 × 2 2
&#
2 &% 1 && 1 0()* %
0%
③
① =
=
2
1 0
&$
1 1
0 1
1 × 11 ) 1 1 1 1= 11 1 1 1 1 1 %&'(( = × 1 1 1 1 5$% % ) G()* JKLM(0 GP'( % %
G'( $
()* ()* 5JKLM(0 G "#$%(! #$ $%$) )
⑦
()* GP %
G'( $5$%
'( GP'( $G#
'(()* 0'( GG $ $P %
1 1
2 1⨀ 11 ×= 1 2 = 1
G'( P $
(f) Initial weight matrix 4 4 of F& (!"# )
()* GP $
1 = 4
⨀
&! &% ! " ! && "
0#'( 5#$0# G()* 0#()* PG#'( (!!"#$ ) $P
⑥
⑧⑦
(i) Constrained sequence counting (!!# → !"# ) (i) Constrained sequence counting (!!# → !"# )
&#
1
1 1
0
+
1 1
'( GP'( !0!"#$ $ $
④
&$
0 2 1
0#()*
0# ⑥
between F& and F$ (! ) '( "# GP $
2
0$'(
&! &% 1 && 1
⨀
1
0#'( ) 0# G& P (!!"#$
(h) Regularized weight matrix of F& (!"# )
⑧
=
2 ()* 0P #
&#
&$
1 0 0 1 G'('+ &
③
=
④③
2 2
=
0()* %
1 1 1 1 11 11 1 1 1 1 1 1 + = × 1 1 1 1 1 11 1 1 1 (g) Reachability1matrix
1 1
②
1 %&'(( 1 × 1 1) ×= 1 2 = 11 1
2
×
5$%
G'('+ &
'( && 0$ 0&()* 5#$ 0()* 0#'( 0'( 0# 05#()* 0$ 00()* 0()* $% $ 0#'( &# &$ 0#'( #$ $# $ $ 5 : &" 1 1 &% 1 1 ④%& ⑤ ⑥ ④ ⑤ && 1 1 (e) Constrained (d) Transition matrices (!" )sequence counting (!! → !" ) (e) Constrained sequence counting (!! → !" ) P.$ [F$ ]
③
1 1 11 %&'(( ) ⨀x× 1 = = 1 11 11
⑥⑤
&$ 1
"
v0 A &"
&$
5#$ : &! 1 v1 1 v25$% : &#v3 1v4
2 1 0$ : &# 1 ⨀
① (c) Weight vectors (! )
(c) Weight vectors (!"tree ) ! (! ) (a) Query ! !#
) #
1
u
0 A 0$ : &# 1
&! 1
(f) Initial weight matrix of F& (!"# ) &# &% 1 && 1
&$
1 1
P.$ [F$ ] (g) Reachability matrix between F& and F$ (!"# ) &! &% ! " ! && "
&#
&$
1 0 0 1 G&
(h) Regularized weight matrix of F& (!"# )
Figure 2: Example of matrix-based counting: dashed edges denote non-tree edges and red edges present only in 𝑇𝑄 1 and 𝑇𝐶1 . the weight vector of 𝑢𝑐′ by multiplying each candidate’s weight in V𝑢𝑐′ with the sum of weights propagated from the candidates in 𝐶 (𝑢𝑐 | 𝑢, 𝑣) via A𝑢𝑐′ 𝑢𝑐 . Since 𝐶 (𝑢𝑐′ | 𝑢, 𝑣) and 𝐶 (𝑢𝑐 | 𝑢, 𝑣) are restricted to candidates under the condition that 𝑢 matches 𝑣, this weight accumulation simultaneously satisfies both constraints imposed by the tree-edges from 𝑢 and the non-tree edge of (𝑢𝑐′ , 𝑢𝑐 ). Induction: By repeating the propagation step backwards through the sequence, each step incorporates exactly one non-tree edge constraint. The final value is a sum of product terms where each term represents a combination of candidates that satisfies all tree and non-tree constraints. Since each valid combination is counted exactly once through these algebraic operations, the final summation of entries in the head vector yields the exact total count. □
4.2
Extending to Generic Case
Motivation. When a query vertex has an indirect predecessor, the indirect and direct edges along the sequence jointly impose constraints on the candidate matching. For example, in the query tree 𝑇𝑄 1 shown in Figure 2(a), where the red edges are included in the tree, 𝑢 4 has 2 indirect predecessors 𝑢 1 and 𝑢 2 . The indirect edges (𝑢 1, 𝑢 4 ), (𝑢 2, 𝑢 4 ) and the sequence ⟨𝑢 1, 𝑢 2, 𝑢 3, 𝑢 4 ⟩ jointly impose matching constraints on this sequence. We propose to use a weight matrix to encode the indirect constraints on a vertex, and further approximate the interactions among these constraints through an enhanced regularized weight matrix. With carefully designed propagation rules, the weight matrix accumulates valid weights during the traversal of the constrained sequence to support approximate counting with indirect constraints. Incorporating Indirect Predecessor Constraints. In the following, we proceed to deal with the indirect predecessor-successor relationships between vertices presented in the constrained sequence. For a vertex with indirect predecessors, we need to properly take
the connection constraints brought by such predecessors into consideration during the matrix propagation. To this end, we introduce the concept of weight matrix, which is essentially an advanced version of the weight vector by integrating both the candidate weights and the connection status with the indirect predecessors. Definition 4.8 (Weight Matrix). Given a vertex 𝑢𝑐 in the constrained sequence 𝑠 with 𝑚 indirect predecessors ⟨𝑢𝑝 1 , ..., 𝑢𝑝𝑚 ⟩ ordered as in 𝑠, the weight matrix of 𝑢𝑐 , F, consists of 𝑚 column-wise concatenated blocks, where the 𝑘-th block F[𝑢𝑝𝑘 ] corresponds to 𝑢𝑝𝑘 : • For an indirect predecessor 𝑢𝑝 , the corresponding block F[𝑢𝑝 ] is a matrix with |𝐶 (𝑢𝑐 | 𝑢, 𝑣)| rows and |𝐶 (𝑢𝑝 | 𝑢, 𝑣)| columns, where an entry F[𝑢𝑝 ] 𝑖 𝑗 = 1 if there is a candidate edge between the 𝑖-th candidate of 𝑢𝑐 and the 𝑗-th candidate of 𝑢𝑝 , and 0 otherwise. • For the block corresponding to the first indirect predecessor 𝑢𝑝 1 , we further set its nonzero entries to the corresponding candidate weights 𝑊 (𝑢𝑐 , 𝑣𝑐 ) of 𝑢𝑐 . The weight accumulation in the propagation process is performed through this primary block, while the remaining blocks are used to verify connectivity for non-tree constraints. This ensures that the candidate weights of 𝑢𝑐 are accumulated exactly once in the Cartesian-product space. Example 4.9. For the constrained sequence ⟨𝑢 1, 𝑢 2, 𝑢 3, 𝑢 4 ⟩ in 𝑇𝑄 1 shown in Figure 2(a) and its corresponding candidate trees 𝑇𝐶1 in Figure 2(b), the initial weight matrix F𝑖𝑛𝑖𝑡 of 𝑢 4 is presented in Figure 2(f), 4 where the red dashed line separates the blocks corresponding to different indirect predecessors. Since 𝑢 4 has two indirect predecessors 𝑢 1 and 𝑢 2 , its weight matrix has two blocks. Considering 𝑢 1 , since there is only one vertex 𝑣 1 in 𝐶 (𝑢 1 | 𝑢 0, 𝑣 0 ) which connects to both candidates 𝑣 5 and 𝑣 6 in 𝐶 (𝑢 4 | 𝑢 0, 𝑣 0 ), the corresponding block is a 2 × 1 matrix of ones, representing both connectivity and candidate weights. The block for 𝑢 2 is constructed similarly based on connectivity only. Regularizing Weight Matrix. Recall that the overall counting process is formulated as a matrix-based propagation backwards
Algebraic Subgraph Counting
along the constrained sequence. Based on Definition 4.8, the candidate weights of the successor and the non-tree constraint from its first indirect predecessor are embedded in the first block of the successor’s weight matrix F, where weight accumulation in propagation takes place. The constraints from other indirect predecessors are encoded in separate blocks. These constraints from different indirect predecessors cannot be applied to the successor’s candidates simultaneously for valid weight accumulation during the propagation process. Example 4.10. In 𝑇𝑄 1 shown in Figure 2(a), both 𝑢 1 and 𝑢 2 are indirect predecessors of 𝑢 4 . In the candidate tree 𝑇𝐶1 in Figure 2(b), the weights of 𝑢 4 ’s candidates 𝑣 5 and 𝑣 6 can propagate backwards along the sequence to 𝑣 2 and 𝑣 3 . This propagation satisfies the indirect constraint between 𝑢 1 and 𝑢 4 , since 𝑣 1 connects to both 𝑣 5 and 𝑣 6 , but it violates the indirect constraint between 𝑢 2 and 𝑢 4 , because 𝑣 2 is not adjacent to 𝑣 6 , and 𝑣 3 is not adjacent to 𝑣 5 . The indirect constraints of (𝑢 1, 𝑢 4 ) and (𝑢 2, 𝑢 4 ) are encoded in different blocks of the weight matrix, which may not be satisfied simultaneously on 𝑣 5 and 𝑣 6 . To jointly enforce the constraints from multiple indirect predecessors on a successor without introducing complex multi-constraint checks during propagation, we develop regularization techniques for the weight matrix. Specifically, we apply a scaling ratio matrix to the block of the first indirect predecessor in the weight matrix to approximately offset invalid weights before the propagation begins. Essentially, each ratio estimates how much of the weight propagated from a successor candidate to the corresponding candidate of the first indirect predecessor remains valid after jointly enforcing all indirect-predecessor constraints along the sequence. To facilitate the computation, we make use of the following definition reachability matrix, which records the reachability situation between the candidates of two query vertices. Definition 4.11 (Reachability Matrix). Given two vertices 𝑢𝑐 and 𝑢𝑐′ in the constrained sequence, the reachability matrix between 𝑢𝑐 and 𝑢𝑐′ , denoted as R𝑢𝑐 [𝑢𝑐′ ], maintains the reachability from vertices in 𝐶 (𝑢𝑐 | 𝑢, 𝑣) to 𝐶 (𝑢𝑐′ | 𝑢, 𝑣) along the constrained sequence 𝑠, which can be easily obtained via successive multiplying the transition matrices between consecutive vertices in the constrained sequence. Example 4.12. Figure 2(g) shows the reachability matrix R𝑢4 [𝑢 2 ] from candidate of 𝑢 4 to the candidates of 𝑢 2 , which is obtained by binarizing (A23 × A34 ) ⊤ , where all non-zero entries are set to 1. It indicates that both candidates 𝑣 5 and 𝑣 6 can reach 𝑣 2 and 𝑣 3 in the candidate space along the constrained sequence. Now, given a vertex 𝑢𝑐 in the constrained sequence with indirect predecessors ⟨𝑢𝑝 1 , ..., 𝑢𝑝𝑚 ⟩ and its initial weight matrix F, we use F[𝑢𝑝𝑖 ] to denote the matrix block corresponding to the indirect predecessor 𝑢𝑝𝑖 . With the reachability matrix, we update the weight matrix as follows. F[𝑢𝑝𝑖 ] ←F[𝑢𝑝𝑖 ] ⊙ (diag(R𝑢𝑐 [𝑢𝑝𝑖+1 ] · 1) −1 · (4) ((F[𝑢𝑝𝑖+1 ] ⊙ R𝑢𝑐 [𝑢𝑝𝑖+1 ]) · R𝑢𝑝𝑖+1 [𝑢𝑝𝑖 ])). In Equation 4, the state at 𝑢𝑝𝑖 depends on 𝑢𝑝𝑖+1 . We examine the reachability of each candidate of 𝑢𝑐 to candidates at 𝑢𝑝𝑖+1 with R𝑢𝑐 [𝑢𝑝𝑖+1 ]. Reachable candidates at 𝑢𝑝𝑖+1 receive the accumulated weights propagated from 𝑢𝑐 ’s candidates. Accordingly, multiplying
by an all-ones vector 1 yields the number of reachable candidates at 𝑢𝑝𝑖+1 from each candidate of 𝑢𝑐 . The resulting vector is broadcast via the 𝑑𝑖𝑎𝑔(·) operation to form a diagonal matrix, whose diagonal entries serve as the denominator for the corresponding row in the ratio computation. The numerator accounts for how much of the accumulated weight propagated from 𝑢𝑐 ’s candidates satisfies the indirect constraints at 𝑢𝑝𝑖+1 and can further propagate to 𝑢𝑝𝑖 . This is computed in the second line of Equation 4: F[𝑢𝑝𝑖+1 ] encodes the indirect constraints of 𝑢𝑝𝑖+1 , while R𝑢𝑐 [𝑢𝑝𝑖+1 ] captures reachability from 𝑢𝑐 to 𝑢𝑝𝑖+1 . Candidates that are both reachable and indirect constraint-satisfying are propagated to 𝑢𝑝𝑖 through R𝑢𝑝𝑖+1 [𝑢𝑝𝑖 ], forming the numerator of the ratio computation. The resulting ratio matrix is combined with F[𝑢𝑝𝑖 ] via the Hadamard product, thereby incorporating the indirect constraints at 𝑢𝑝𝑖 and producing an updated matrix block, which is then used to update 𝑢𝑝𝑖 −1 . Since only the first block of F contains candidate weights of 𝑢𝑐 , we retain only the regularized update of this block in Equation 4. The remaining blocks encode connectivity, where only the zero/non-zero distinction matters, so their original values are kept. Example 4.13. Based on the initial feature matrix in Example 4.9, we apply a scaling ratio matrix on the block F𝑖𝑛𝑖𝑡 4 [𝑢 1 ] of 𝑢 1 . For the reachability matrix R𝑢4 [𝑢 2 ] shown in Figure 2(g), we compute the row sums, resulting in the vector [[2], [2]], which serves as the perrow denominators of the ratio matrix. Next, we take the Hadamard product of R𝑢4 [𝑢 2 ] and F𝑖𝑛𝑖𝑡 4 [𝑢 2 ], obtaining [[1, 0], [0, 1]], which incorporates the indirect constraint on (𝑢 2, 𝑢 4 ). To evaluate how much of the accumulated weight can further propagate to 𝐶 (𝑢 1 | 𝑢 0, 𝑣 0 ), we multiply the obtained matrix by R𝑢2 [𝑢 1 ] = [[1], [1]], yielding [[1], [1]]. Dividing each entry by its corresponding row denominator produces the ratio matrix [[1/2], [1/2]]. Finally, the ratio matrix is applied to F𝑖𝑛𝑖𝑡 4 [𝑢 1 ] via Hadamard product, which further incorporates the indirect constraint between (𝑢 1, 𝑢 4 ). The resulting matrix is the final feature matrix F4 , as shown in Figure 2(h). Propagating Weight Matrix. After setting the weight matrix properly, we propagate it via matrix multiplication along the reverse direction of the constrained sequence until reaching the first vertex, where the accumulated weights are extracted as the homomorphism count. Because each vertex in the constrained sequence might have its own weight matrix, to avoid ambiguity, we call the propagated weight matrix as passing matrix P. During the propagation, we need to execute the following two types of operations accordingly, including weight aggregation and weight extraction. • Weight Aggregation. Let 𝑢𝑐 be the current vertex in the constrained sequence, and 𝑢𝑐′ its direct predecessor during the propagation. We incorporate the weight matrix F of 𝑢𝑐′ into the current passing matrix P to update the accumulated weights, as detailed in algorithm 2. Initially, we handle the boundary cases where either matrix acts as a vector because it lacks indirect predecessors. Case 1: P and F have no associated indirect predecessors. In this scenario, since P has not accumulated any indirect predecessors along the propagation path (line 1), both P and F act purely as weight vectors (Definition 4.3). Consequently, we simply compute their 2 of Example 4.6. Hadamard product (line 3), identical to step ○ Case 2: F has associated indirect predecessors but P does not. F functions as a weight matrix while P remains a weight vector (line 5)
Guo et al.
under this condition. We row-scale the first block of F by P, and then F takes over as the updated passing matrix (line 7–line 8). Case 3: P has associated indirect predecessors, but F does not. P has already become a matrix but F remains a vector (line 11) in this scenario, we directly apply F to the first block of P via row-wise scaling to accumulate the candidate weights (line 12). Case 4: both P and F have associated indirect predecessors. P and F both function as weight matrices (Definition 4.8) in this scenario. We identify their predecessor sets U𝑃 and U𝐹 , along with their earliest predecessor 𝑢𝑃,1 and 𝑢 𝐹,1 . We first process any overlapping predecessors (𝑢 ∈ U𝑃 ∩ U𝐹 ) by performing an element-wise Hadamard product on their corresponding blocks. This step merges their connectivity and weights simultaneously ( line 16–line 18). A critical requirement is that only the block corresponding to the earliest indirect predecessor in the constrained sequence acts as the primary block carrying the accumulated weights. To ensure this, we compare the sequence positions of 𝑢 𝐹,1 and 𝑢𝑃,1 . If 𝑢 𝐹,1 appears earlier (𝑢 𝐹,1 ≺ 𝑢𝑃,1 ) in the constrained sequence, F[𝑢 𝐹,1 ] becomes the new primary block. We extract the row weights vrow from the former primary block P[𝑢𝑃,1 ], row-broadcast them into F[𝑢 𝐹,1 ], and degrade P[𝑢𝑃,1 ] into a pure 0/1 connectivity block via binarization (line 19–line 23). Conversely, if 𝑢𝑃,1 remains the primary block, we extract the weights from F[𝑢 𝐹,1 ], apply them to P[𝑢𝑃,1 ], and binarize F[𝑢 𝐹,1 ] (line 24–line 29). Finally, any remaining non-overlapping blocks belonging solely to the weight matrix (𝑢 ∈ U𝐹 \ U𝑃 ) are appended to the passing matrix P column-wise, expanding its predecessor set accordingly (line 30–line 32). • Weight Extraction. During backward propagation along the constrained sequence, if the predecessor 𝑢𝑐′ recorded in the current passing matrix P is encountered, it signifies a cycle closure. We perform a weight extraction operation as detailed in algorithm 3. At this stage, the block P[𝑢𝑐′ ] forms a square matrix. Its diagonal elements indicate whether the corresponding candidates of 𝑢𝑐′ can successfully complete a cycle traversal under the indirect constraints, as in the cycle (𝑢 2 → 𝑢 4 → 𝑢 3 → 𝑢 2 ) of 𝑇𝑄 1 in Figure 2(a). We extract these diagonal elements into a vector vdiag (line 3) to capture the valid states of 𝑢𝑐′ ’s candidates. Since the information at 𝑢𝑐′ has now been fully evaluated and aggregated, we simplify the passing matrix by completely removing the block corresponding to 𝑢𝑐′ (line 5). The subsequent update of P depends on whether other indirect predecessors remain. Case 1: Other blocks remain in P. In this scenario (line 6), the presence of remaining blocks implies that 𝑢𝑐′ ’s block was not the primary carrier of the accumulated weights. Nevertheless, its constraint validation is crucial. Therefore, we binarize the extracted vector to retain only connectivity information (line 8) and apply it via rowwise scaling to the first block of P (line 9). This operation effectively filters out the accumulated weights of candidates that fail to satisfy the cycle constraints. Case 2: 𝑢𝑐′ ’s block was the only block in P. In this scenario (line 10), the extracted vector vdiag fully represents the newly accumulated weights. In this case, vdiag entirely becomes the updated passing matrix P (line 11) and will continue propagating backward along the constrained sequence.
Algorithm 2: WeightAggregation Input: Passing matrix P at 𝑢𝑐 , Weight matrix F of current direct predecessor 𝑢𝑐′ Output: Updated passing matrix P 1 if P has no associated indirect predecessors then 2 if 𝑢𝑐′ has no indirect predecessor then // Case 1 3 P ← P ⊙ F; 4 return P 5 6 7 8 9
U𝐹 ← indirect predecessors of blocks in F ; // Case 2 𝑢 𝐹,1 ← the first indirect predecessor in U𝐹 ; F[𝑢 𝐹,1 ] ← F[𝑢 𝐹,1 ] ⊗row P; return F U𝑃 ← indirect predecessors of blocks in P;
10 𝑢 𝑃,1 ← the first indirect predecessor in U𝑃 ;
if 𝑢𝑐′ has no indirect predecessor then 12 P[𝑢𝑃,1 ] ← P[𝑢𝑃,1 ] ⊗row F; 13 return P
11
14
U𝐹 ← indirect predecessors of blocks in F ;
// Case 3
// Case 4
15 𝑢 𝐹,1 ← the first indirect predecessor in U𝐹 ;
foreach 𝑢 ∈ U𝑃 ∩ U𝐹 do F[𝑢] ← P[𝑢] ⊙ F[𝑢]; 18 P[𝑢] ← F[𝑢];
16
17
if 𝑢 𝐹,1 ≺ 𝑢𝑃,1 then vrow ← ExtractRowWeights(P[𝑢𝑃,1 ]); 21 F[𝑢 𝐹,1 ] ← F[𝑢 𝐹,1 ] ⊗row vrow ; 22 P[𝑢𝑃,1 ] ← Binarize(P[𝑢𝑃,1 ]); 23 P[𝑢 𝐹,1 ] ← F[𝑢 𝐹,1 ];
19
20
else vrow ← ExtractRowWeights(F[𝑢 𝐹,1 ]); 26 P[𝑢𝑃,1 ] ← P[𝑢𝑃,1 ] ⊗row vrow ; 27 F[𝑢 𝐹,1 ] ← Binarize(F[𝑢 𝐹,1 ]); 28 if 𝑢 𝐹,1 ∈ U𝑃 then 29 P[𝑢 𝐹,1 ] ← F[𝑢 𝐹,1 ];
24
25
foreach 𝑢 ∈ U𝐹 \ U𝑃 do Append(P, F[𝑢]); 32 U𝑃 ← U𝑃 ∪ {𝑢};
30
31
33
return P
Algorithm 4 outlines the overall procedure for computing the candidate weight𝑊 (𝑢, 𝑣) under both sequence constraints and standard tree constraints. The algorithm performs a backward propagation along the constrained sequence 𝑠. We initialize the process at the last vertex 𝑢𝑐 of sequence 𝑠, setting its weight matrix as the initial state of the passing matrix P (line 1–line 2). During the reversed traversal of 𝑠, for each adjacent pair (𝑢𝑐′ , 𝑢𝑐 ), we first propagate the passing matrix P from 𝑢𝑐 to its direct predecessor 𝑢𝑐′ by multiplying it with the transition matrix A𝑢𝑐′ ,𝑢𝑐 (line 4). Then, we generate the weight matrix F for 𝑢𝑐′ (line 5). If 𝑢𝑐′ closes a cycle—meaning it is an indirect predecessor previously visited—we refine the passing matrix via the WeightExtraction operation to capture cycle constraints (line 7). Subsequently, the
Algebraic Subgraph Counting
Algorithm 3: WeightExtraction Input: Passing matrix P at 𝑢𝑐 , Current direct predecessor 𝑢𝑐′ Output: Updated passing matrix P 1 U𝑃 ← indirect predecessors of blocks in P; ′ 2 if 𝑢𝑐 ∈ U𝑃 then 3 vdiag ← ExtractDiag(P[𝑢𝑐′ ]); 4 U𝑃 ← U𝑃 \ {𝑢𝑐′ }; 5 RemoveBlock(P, 𝑢𝑐′ ); 6 if U𝑃 ≠ ∅ then // Case 1 7 𝑢𝑃,1 ← the first indirect predecessor in U𝑃 ; 8 vdiag ← Binarize(vdiag ); 9 P[𝑢𝑃,1 ] ← P[𝑢𝑃,1 ] ⊗row vdiag ; 10 11
12
else P ← vdiag ;
// Case 2
return P
Algorithm 4: SeqenceCounting Input: 𝑢, 𝑣, 𝐶𝑆, 𝑠 Output: 𝑊 (𝑢, 𝑣) 1 𝑢𝑐 ← the last vertex in 𝑠; 2 P ← WeightMatrix(𝑢𝑐 ); ′ 3 foreach (𝑢𝑐 , 𝑢𝑐 ) in the reversed traversal of 𝑠 do 4 P ← A𝑢𝑐′ 𝑢𝑐 P; 5 F ← WeightMatrix(𝑢𝑐′ ); 6 if 𝑢𝑐′ closes a cycle then 7 P ← WeightExtraction(P, 𝑢𝑐′ ); 8
P ← WeightAggregation(P, F, 𝑢𝑐′ );
9 𝑤 con ← Sum(P); 10 𝑤 norm ← IsolatedChildrenWeight(𝑢, 𝑣);
3 Since 𝑢 2 is the indirect predecessor of 𝑢 4 , whose A23 , obtaining P in ○. constraint on 𝑢 4 is enforced in the passing matrix, we extract the block 4 the diagonal (circled in red) of the block corresponding of 𝑢 2 . In Step ○, to 𝑢 2 is extracted, obtaining V𝑑𝑖𝑎𝑔 (line 3 of algorithm 3), and then this block is removed from P (line 5 of algorithm 3). V𝑑𝑖𝑎𝑔 is first binarized (line 8 of algorithm 3) and applied to the first block of P (line 9 of algorithm 3), yielding the updated P after weight extraction as shown 5 In Step ○, 6 the initial weight vector V2 of 𝑢 2 is row-scaled in Step ○. into P, which completes the weight aggregation at 𝑢 2 and produces the updated P. After that, the obtained P continues to propagate backward to 𝑢 1 through A12 . At this point, 𝑢 1 is the only remaining indirect predecessor encoded in the passing matrix, the diagonal (circled in 7 red) of P is extracted to form the vector V𝑑𝑖𝑎𝑔 as shown in Step ○, which is then combined with V1 via a Hadamard product. Finally, 8 yields the constrained sequence summing the entries of P in Step ○ count, which equals to 2. This value is exactly the true homomorphism count of the constrained sequence. Remark. It should be remarked that the homomorphism counts obtained from the constrained sequence counting are not exact. When the constrained sequence contains indirect constraints, multiple indirect constraints can simultaneously affect a single vertex. Although the regularized weight matrix is introduced to precompute such joint constraint effects before propagation, it can only approximate the proportion of valid weights propagated based on the connected paths. In fact, the valid weight contributed by each connected path is different. The counting rules for such a complicated structure might be further optimized in future work.
4.3
Complexity Analysis
Theorem 4.15. Given a connected component 𝑆 under a query vertex in the query tree, the time complexity of the overall counting process for 𝑆 is polynomial to the size of 𝑆 and the candidate space.
11 𝑊 (𝑢, 𝑣) ← 𝑤 con × 𝑤 norm ; 12
return 𝑊 (𝑢, 𝑣)
state is updated by fusing P and F via the WeightAggregation operation (line 8). Once the propagation reaches the head vertex of 𝑠, the finalized passing matrix P degrades into a vector containing the accumulated valid weights. We sum all elements of this vector to obtain the total sequence-constrained weight 𝑤 con (line 9). Furthermore, the query vertex 𝑢 may also possess isolated children that are completely independent of any Type-I non-tree constraints. We compute their normal candidate weight 𝑤 norm using Equation 2 (line 10). Finally, the product of 𝑤 con and 𝑤 norm yields the final matching weight 𝑊 (𝑢, 𝑣) (line 11). Example 4.14. Continuing Example 4.13, where the weight matrix F4 is constructed in Figure 2(h). The propagation process is illustrated in Figure 2(i). Starting from 𝑢 4 , the weight matrix F4 is propagated to 1 The weight vector 𝑢 3 by multiplication with A34 , illustrated in Step ○. of 𝑢 3 , which encodes the weights of its candidates, is then incorporated into the first block of P through row-broadcast as specified in line 12 2 The of algorithm 2, producing the updated passing matrix in Step ○. passing matrix is subsequently propagated backward to 𝑢 2 through
Proof. In the counting process, we only conduct two rounds of matrix propagation. At each step of the propagation, we conduct a constant number of matrix operations (i.e., matrix multiplication), which can be done in cubic time of the candidate size of query vertices. Since the number of steps in each propagation is bounded by the query graph size, the theorem is immediate. □ To assess the efficiency of the constrained sequence counting, we conduct a preliminary experiment on a synthetic dataset. Specifically, we use a 32-clique query graph and a 32-clique data graph, where vertices share the same label. By definition, the exact count is 32!. We evaluate this instance using two representative subgraph matching algorithms recommended in recent survey papers [50, 57]. None of them can finish within two hours. The state-of-the-art sampling-based subgraph counting method Fastest [47] also fails to return a valid result within 100 seconds. In contrast, our method computes the count in 7.11 seconds with a q-error of 1.98, where the majority of the runtime is spent on candidate space construction and refinement, and the counting phase only spends 0.52 seconds.
5
Handling Type-II Non-tree Edge Constraints
In this section, we introduce the lowest common ancestor (LCA for short) to handle the Type-II non-tree edges constraints.
Guo et al.
Recall that a Type-II non-tree edge connects two vertices with different parents in the query tree. However, the two vertices must have an LCA vertex in the tree. Clearly, the paths from the LCA vertex to the two endpoints form a triangle together with this nontree edge. A necessary condition for a valid match of the query graph is that such a triangle must also exist in the candidate space. However, in the candidate tree-based counting method defined by Equation 2, the weights of candidate vertices are propagated upward along the tree paths without considering such non-tree edge constraints. To remove the false positive counting, we develop an LCA ratio computation method. Given a Type-II non-tree edge, we check if the two endpoints of its matching edge can reach the same LCA candidate. By calculating the positive ratio, we can therefore estimate the valid subgraph homomorphism count. To compute the LCA ratio, we first find the descendants for each candidate of an LCA vertex using a dynamic programming method. Then, based on the LCA descendants, we compute the LCA ratio for each candidate by using set intersections to count descendant pairs that form valid triangles across endpoints.
Algorithm 5: Find LCA Descendants Input: The query tree 𝑇𝑄 , corresponding 𝐶𝑆 Output: LCA descendants for each LCA path 1 foreach 𝑢 in BFS Traversal of 𝑇𝑄 do 2 if 𝑢 is involved in any LCA path then 3 𝑢𝑝 ← Parent(𝑢); 4 if 𝑢𝑝 exists then 5 All_current_descendants ← {𝑣 𝑝 | 𝑣 𝑝 ∈ 𝐶𝑢𝑟𝑟𝑒𝑛𝑡_𝑑𝑒𝑠𝑐𝑒𝑛𝑑𝑎𝑛𝑡𝑠 [𝑃], ∀𝑃 covering 𝑢𝑝 }; 6 foreach 𝑣 𝑝 ∈ All_current_descendants do 7 Cache[𝑣 𝑝 ] = 𝐶 (𝑢 | 𝑢𝑝 , 𝑣 𝑝 ); 8 9
10 11 12 13
5.1
foreach LCA path 𝑃 passing through 𝑢 do Current_descendants[𝑃] ← {𝐶𝑎𝑐ℎ𝑒 [𝑣 𝑝 ] | 𝑣 𝑝 ∈ 𝐶𝑢𝑟𝑟𝑒𝑛𝑡_𝑑𝑒𝑠𝑐𝑒𝑛𝑑𝑎𝑛𝑡𝑠 [𝑃]}; if 𝑢 is the LCA of path 𝑃 then Current_descendants[𝑃] ← 𝐶 (𝑢); if 𝑢 is the Endpoint of path 𝑃 then Save Current_descendants[𝑃];
Finding LCA Descendants
For each LCA candidate, we compute its descendant candidates located at the endpoints of a non-tree edge. This is done by propagating reachable candidates downward along the tree path from the LCA to the endpoint. To speed up the computation, we adopt a level-by-level manner where the results computed in the previous level can be reused in the current level. This can be implemented as a dynamic programming on the candidate tree. Algorithm Details. The computation details are summarized in algorithm 5. For each Type-II non-tree edge, we form two LCA paths, each connecting the LCA to one of its endpoints. We perform a one-pass top-down traversal of the query tree 𝑇𝑄 in BFS order, and compute the descendants of all LCA paths simultaneously. For each LCA path 𝑃, we maintain a path state 𝐶𝑢𝑟𝑟𝑒𝑛𝑡_𝑑𝑒𝑠𝑐𝑒𝑛𝑑𝑎𝑛𝑡𝑠 [𝑃], which functions as a mapping from the LCA candidates to their descendants at the current query vertex. As we traverse downwards, the path states at the parent 𝑢𝑝 are propagated to the current vertex 𝑢. Specifically, when visiting 𝑢, we first identify all 𝑣 𝑝 ∈ 𝐶 (𝑢𝑝 ) involved in any LCA path passing through its parent 𝑢𝑝 (line 5) and cache their child candidates 𝐶 (𝑢 | 𝑢𝑝 , 𝑣 𝑝 ) (line 7) to ensure each candidate edge is processed only once. Subsequently, for each path 𝑃 passing through 𝑢, we update the path state 𝐶𝑢𝑟𝑟𝑒𝑛𝑡_𝑑𝑒𝑠𝑐𝑒𝑛𝑑𝑎𝑛𝑡𝑠 [𝑃] by replacing the parent candidates 𝑣 𝑝 with their cached child candidates 𝐶𝑎𝑐ℎ𝑒 [𝑣 𝑝 ] (line 9). If 𝑢 is the LCA vertex of a path 𝑃, its state is initialized as an identity mapping of 𝐶 (𝑢) (line 11). Once 𝑢 reaches the endpoint of a path 𝑃, the current state is saved (line 13) as the mapping from the LCA candidates to the endpoint descendants. Lemma 5.1. The complexity of Algorithm 5 is 𝑂 (|𝐸𝐶 |+|𝑉𝑄 ||𝐸𝑄 ||𝑉𝐶 |), where 𝑉𝑄 and 𝐸𝑄 are the query graph vertex and edge sets, and 𝑉𝐶 and 𝐸𝐶 are the vertex set and edge set in the candidate space. Proof. The process of finding LCA descendants traverses the candidate tree in a top-down manner. By employing the cache mechanism, each candidate edge is scanned only once, resulting in a complexity of 𝑂 (|𝐸𝐶 |). During the traversal, the algorithm propagates path states at each 𝑢. The number of LCA paths passing
through a single vertex 𝑢 is proportional to the tree depth 𝑂 (|𝐸𝑄 |). For each LCA path, the path state maintains a mapping from LCA candidates to the candidates of the current vertex 𝑢, where the total number of involved candidates is 𝑂 (|𝑉𝐶 |). Therefore, the complexity of path states propagation is (|𝑉𝑄 ||𝐸𝑄 ||𝑉𝐶 |). The overall complexity of the algorithm is 𝑂 (|𝐸𝐶 | + |𝑉𝑄 ||𝐸𝑄 ||𝑉𝐶 |). □ Discussion. The high efficiency of Algorithm 5 is achieved by reusing the intermediate results as follows. • For two non-tree edges sharing the same LCA-endpoint path, the descendant results along a shared path only need to be computed once because the endpoints of a non-tree edge are decomposed. • For two LCA-endpoint paths sharing the LCA, if one endpoint is the ancestor of the other, the descendant candidates computed for the ancestor endpoint can be reused for the other. • In one LCA-endpoint path, multiple LCA candidates may map to the same descendant candidate at an intermediate vertex. The computation of that descendant candidate’s child candidates is performed only once and shared across all such LCA candidates.
5.2
Computing LCA Ratio
We view the tree paths from the LCA to the two non-tree endpoints as virtual edges. Together with the non-tree edge, these two virtual edges form a triangle. Our goal is to count, for each LCA candidate, how many combinations of endpoint candidates form valid triangles. Specifically, for a given LCA candidate, we first extract the candidate sets of the two endpoints of a Type-II non-tree edge that connected to it through the candidate tree paths. The total number of possible triangles is given by the product of the sizes of these two candidate sets, which aligns with the counting logic of Equation 2. We then explicitly examine how many candidate edges exist between the two endpoint candidate sets. Each such candidate edge, together with the connections from its endpoints to the LCA candidate, forms a valid triangle.
Algebraic Subgraph Counting
A single LCA vertex may be associated with multiple non-tree edges, some of which share a common endpoint. We group such nontree edges by their shared endpoint and treat this endpoint as an anchor. Within each group, we fix one candidate path from the LCA to the anchor, and compute the fraction of candidate combinations of the other endpoints that satisfy all triangle constraints induced by the group. The final LCA ratio for an LCA candidate is obtained by multiplying the ratios over all anchor groups. Lemma 5.2. Given an LCA-non-tree triangle, the LCA ratio can be computed in 𝑂 (|𝐸𝐶 ||𝑉𝐶 |) time. Proof. For each candidate edge between the LCA and the anchor endpoint, we compute the intersection of the candidate neighbors of the LCA candidate and the anchor candidate. We conduct such intersections for 𝑂 (|𝐸𝐶 |) anchored candidate edges, each consuming 𝑂 (|𝑉𝐶 |) time. Thus, the lemma is complete. □ Remark. By integrating Type-I and Type-II non-tree constraints with the existing tree-edge constraints in the candidate tree framework, our structural module bridges the gap between tree homomorphism and graph homomorphism, extending the candidate tree framework to full subgraph homomorphism counting. This extension is non-trivial, since tree homomorphism counting is polynomial-time solvable, whereas subgraph homomorphism counting itself is a #P-hard problem. Compared with Fastest, which relies on sampling over simple tree homomorphism counts, our approach moves the more complex structural-constraint computation into the counting stage and leaves only the injectivity constraints to the sampling module. For queries that do not involve injectivity violations, the sampling stage can even be skipped entirely.
6
Local Sampling
While the non-tree edge constraints can be processed by the matrix computation-based method, the isomorphism constraint remains unsolved. This is because the candidate tree-based counting only propagates and aggregates vertex weights, rather than the specific matching states. In this section, we introduce a local sampling approach to estimate the proportion of isomorphism count within the homomorphism count. Unlike previous sampling-based methods, which usually conduct sampling over the entire candidate space, we conduct sampling in a local candidate space, where isomorphism constraints arise. In particular, isomorphism violations can only occur among candidates of the same-label query vertices; even within this set, clique-structured matches cannot violate injectivity, since no self-loops exist. Query vertices with the same label may appear at different positions in the query graph and are not necessarily adjacent. Extracting them in isolation would break the structural connectivity constraints among them. To preserve these inherent structural constraints, we directly identify a connected subgraph within the original query graph that encompasses all these same-label vertices. This connected subgraph naturally consists of both tree edges and non-tree edges from the query graph, seamlessly maintaining the exact topological relationships. As a result, the target same-label vertices form a locally connected structure where each vertex is properly constrained by its neighbors. We then perform sampling
Table 1: Statistics of Datasets. Dataset Yeast HPRD Human WordNet DBLP Youtube Eu2005 Patents Twitter Friendster
|𝑉𝐺 |
|𝐸𝐺 |
|Σ|
𝑑
|𝑉𝑄 |
3,112 9,460 4,674 76,853 317,080 1,134,890 862,664 3,774,768 41,652,230 65,608,366
12,519 34,998 86,282 120,399 1,049,866 2,987,624 16,138,468 16,518,947 1,202,513,344 1,806,067,135
71 307 44 5 15 25 40 20 100 100
9.0 7.4 36.9 3.1 6.6 5.3 37.4 8.8 57.7 55.1
4 to 32 4 to 32 4 to 20 4 to 20 4 to 16 4 to 32 4 to 8 4 to 32 4 4 to 8
exclusively within the candidate space induced by this locally connected subgraph to estimate the isomorphism ratio. The sampling procedure proceeds by following the traversal order of the connected subgraph. When sampling a query vertex, its feasible sampling region is bounded by its already visited neighbors, meaning that earlier sampled vertices directly constrain the subsequently sampled ones. To optimize this process, we apply a greedy selection rule: among its already visited neighbors, we select the one with the strongest constraint as its tree parent to induce the feasible sampling region, from which a candidate vertex is then uniformly sampled, with the inverse sampling probability accumulated throughout the sampling process. The pre-computed homomorphism weights 𝑊 (𝑢, 𝑣) from the prior counting phase, used in a similar way to [47], are utilized to accurately reflect the matching distribution and continuously monitor the statistical variance of the random walks. The isomorphism ratio is then estimated as the weighted fraction of isomorphic samples among the homomorphic ones. For queries where the empirical samples fail to reach statistical confidence, we approximate this ratio by extrapolating the joint constraint-satisfaction probabilities for homomorphic and isomorphic samples separately from their average marginal survival rates. Additionally, a weight cutoff strategy is applied to truncate overly large sample weights, mitigating extreme variance caused by overly dense homomorphic regions.
7 Experiment 7.1 Experiment Setup Datasets. We conduct experiments on 10 real-world datasets. The first 8 datasets have publicly available query sets and are widely used in previous works like [18, 47, 55], while the remaining 2 are large-scale graphs where we generate queries via random walks on data graphs following [18, 50]. To obtain the exact subgraph counts for each query, we use the recommended method in [50]. Given the extreme difficulty of subgraph counting, we set a timeout of 30 minutes [18, 55] for computing the ground-truth counts on the first 8 datasets, and a timeout of 2 hours [18] for the two large-scale graphs. Our experiments are evaluated on queries with obtainable ground-truth counts, which are further grouped into subsets by size for evaluation. Dataset statistics are reported in Table 1.
Guo et al.
Compared Algorithms. We compare ASC with the following approximate algorithms: (1) FlowSC [18], the state-of-the-art learningbased method; (2) Fastest [47], the state-of-the-art sampling-based method; (3) LearnSC [24]; (4) NeurSC [55]; and (5) LSS [60]. We also evaluate subgraph homomorphism counting against three baselines: Fastest [47], Alley [30], and its index-optimized variant AlleyTPI [30]. The source code of FlowSC, Fastest, NeurSC, Alley, AlleyTPI, and LSS are publicly available, and LearnSC is obtained from the authors. All parameter settings follow their default values. Query Setting. Following the common setting [18, 50, 55], we load a data graph and a query graph into memory each time, and report the elapsed time per data-query pair and the q-error max(1,𝑐 ) max(1,𝑐ˆ ) = max( max(1, 𝑐ˆ ) , max(1,𝑐 ) ), where 𝑐ˆ and 𝑐 are the estimated and ground truth counts, respectively. Experiments are conducted on a Ubuntu 22.04.1 system, equipped with an Intel Xeon Silver 4314 CPU @ 2.40GHz with 64 cores.
7.2
Performance Evaluation
Accuracy Comparison. Figure 3 reports the accuracy performance over all datasets. The boxplot presents the distribution of log q-error of a query graph set. The upper bound (Resp. lower bound) of a box is 75% (Resp. 25%) quantile of the q-error, the whiskers contain all q-errors from minimum to maximum, and the line in the box is the median. Values exceeding 105 are clipped, with bold red caps indicating truncated whiskers. Overall, ASC demonstrates consistently higher and more stable accuracy. On the small dataset Yeast, ASC achieves stable performance across query sets of different sizes, attaining the lowest average q-error among all methods and yielding a 3.24× improvement over the state-of-the-art method FlowSC. Although the samplingbased method Fastest performs nearly perfectly on small queries in Yeast, it begins to suffer from sampling failures when the query size reaches 16, resulting in a sharp degradation in accuracy. On the other two small datasets, HPRD and Human, ASC outperforms all learning-based methods, but is slightly worse than Fastest, as Fastest is particularly effective when the sampling space is small. On large datasets such as Youtube and Patents, Fastest performs well only for small query sizes. In contrast, ASC retains high accuracy on larger datasets and for larger queries. For example, on Youtube_32, ASC achieves the best performance with a low q-error, while FlowSC degrades to a substantially higher error and Fastest fails entirely, yielding a q-error about 7 orders of magnitude larger than that of ASC. In terms of learning-based methods, NeurSC and LearnSC are difficult to train on large datasets. Although LSS can handle Twitter, it is significantly outperformed by ASC. FlowSC can process all datasets. While it achieves higher accuracy than other learningbased methods, but returns lower accuracy than ASC on most datasets, except for WordNet. The exception is largely due to the low variance of the true counts in WordNet, where many queries share identical counts. Such distributions favor regression-based learning methods because of their tendency to regress toward the mean [18]. Consequently, while LearnSC and NeurSC do not benefit from this property due to their difficulty in training, FlowSC and LSS perform particularly well on WordNet.
In addition, these learning-based methods rely heavily on supervision, where the performance is greatly influenced by the number of training samples. As a result, even FlowSC, which performs best among learning-based methods, obtains poor accuracy on the Human dataset. Moreover, on Human_20 and DBLP_16, the scarcity of effective training samples prevents the models from learning. In contrast, the performance of ASC and Fastest is not constrained by the number of training samples, and both methods perform well even on datasets with limited samples. Moreover, only ASC and FlowSC can scale to the billion-scale dataset, where ASC still achieves higher accuracy than FlowSC. Counting Range Analysis. We conduct an experiment on a selected dataset to evaluate how the ground truth count of query graphs influences the distribution of q-error. As shown in Figure 4, our algorithm performs consistently across different count ranges. In the range with larger counts, the accuracy becomes slightly worse, mainly due to the scale effect: queries with larger candidate spaces naturally accumulate more approximation error. The learning-based methods FlowSC, LearnSC, and NeurSC tend to overestimate small-count queries, and underestimate large-count queries, which is consistent with the regression-to-the-mean effect commonly observed in regression-based approaches [18]. The active learning strategy of LSS alleviates this issue. The accuracy of Fastest does not show a clear pattern across different counting ranges. Its performance depends on sampling difficulty, and sampling failures are observed in every interval from 102 to 1010 . Query Processing Time. Figure 5 reports the average processing time. A cross indicates that a method encounters training failure or out-of-memory. In general, ASC is highly competitive on both small and large datasets. On small datasets such as Yeast, HPRD, and Human, ASC is faster than Fastest and achieves runtime comparable to the fastest learning-based method, FlowSC. ASC first incorporates structural constraints into a polynomial-time counting stage, which is efficient in practice. Its sampling stage is then performed only on local subspaces where injectivity violations may occur. If a query does not induce any injectivity violations, the sampling stage can be skipped entirely. Fastest performs sampling based on tree homomorphisms, its sampling is conducted over the full candidate space, which is often much larger than the solution space, making isomorphic samples difficult to obtain. For example, in our preliminary experiment of counting a same-label 32-clique in a same-label 32-clique data graph, Fastest did not return a result within 100 seconds, while ASC skipped sampling entirely because this query does not produce injectivity violations. On the WordNet dataset, NeurSC and LearnSC cannot be applied due to excessive training time. Both ASC and FlowSC are slower than Fastest because their filtering technique BipartitePlus is less efficient on graphs with highly concentrated label distributions. On larger datasets such as Youtube, Eu2005, and Patents, ASC and FlowSC exhibit similar runtime and show a clear advantage over Fastest. As the graph size increases, the full candidate space becomes much larger, causing Fastest to spend more time obtaining isomorphic samples. Due to their complex network architectures, NeurSC and LearnSC cannot scale to large datasets. Although LSS uses a lightweight model and is therefore more scalable, it still cannot process the largest dataset, Friendster. Among learning-based approaches, FlowSC offers the
Algebraic Subgraph Counting
ASC
FlowSC
Fastest
LearnSC
NeurSC
LSS
Yeast
HPRD
Human
WordNet
DBLP
Youtube
Eu2005
Patents
Friendster
Figure 3: Evaluating accuracy with the x-axis representing the number of query graph vertices. ASC
FlowSC
Fastest
LearnSC
NeurSC
LSS
ASC
FaSTest
FlowSC
LSS
NeurSC
LearnSC
Table 2: Homomorphism evaluation. Alg.
ASC
Fastest
Alley
AlleyTPI
Data.
xx
xx
xx
xxx
Youtube
Figure 4: Effect of counting range.
Figure 5: Efficiency Performance.
best efficiency and scalability, mainly due to its highly efficient one-pass bottom-up flow-learning model. Overall, ASC achieves efficiency comparable to FlowSC on large and billion-scale datasets such as Twitter and Friendster.
7.3
Testing on Homomorphism
We conduct experiments for subgraph homomorphism counting to evaluate the effectiveness of the main counting module, i.e., matrix computation-based counting (section 4) together with the LCA ratio (section 5). Since computing the exact subgraph homomorphism counts is also computationally expensive, we select those queries with existing subgraph isomorphism ground truth that do not induce isomorphism constraints. Specifically, a query graph will not violate isomorphism constraints if its vertex labels are unique, or if vertices with the same label form a clique. In such cases, the subgraph isomorphism count coincides with the subgraph homomorphism count. We select all queries satisfying these conditions from each dataset for homomorphism testing. We compare ASC with Fastest and the homomorphism counting method Alley and its variant AlleyTPI. As shown in Table 2, ASC provides a better accuracy-efficiency trade-off, it achieves accuracy comparable to Fastest with substantially faster runtime, and clearly outperforms Alley and AlleyTPI in both accuracy and efficiency. Alley and AlleyTPI frequently suffer
xxxx
Yeast HPRD Human WordNet DBLP Youtube Eu2005 Patents Twitter Friendster
q-error
Time
q-error
Time
q-error
Time
q-error
Time
1.05 1.01 1.18 1.00 1.09 2.74 6.48 1.22 1.04 1.15
0.02 0.03 0.08 2.23 0.46 1.32 3.92 8.83 505.05 1×103
1.02 1.00 1.01 1.08 1.01 1.32 1.14 1.05 -
0.04 0.11 0.09 0.25 1.64 6.10 16.69 38.52 -
1×104 22.45 7×108 11.85 -
0.17 0.01 6×104 7×103 -
1×103 13.15 7×108 5.90 -
0.17 0.01 6×104 4×104 -
from sampling failures, which leads to unsatisfactory performance even on small datasets. In addition, AlleyTPI incurs extra overhead for index construction, making it unsuitable for random data-query pairs. We omit their results on larger datasets. Fastest achieves high accuracy but at the cost of substantial runtime. In contrast, ASC attains accuracy comparable to Fastest on most datasets while requiring significantly less time. Although ASC is less accurate than Fastest on Youtube and Eu2005, where the candidate space is large, its accuracy is acceptable given the substantially lower runtime. Although WordNet also has a larger candidate space, it contains only four 4−size queries without isomorphism constraints, and ASC achieves near-zero error on these queries.
7.4
Evaluation of Individual Techniques
Error Characterization of Type-I Module. We conduct error characterization experiments on small datasets Yeast and Human to analyze the factors that affect algebraic counting. However, error characterization requires exact counting of graph homomorphisms, which is extremely time-consuming even on small datasets. To tackle the computational difficulty, we reduce the exact counting to the child neighborhood of a single candidate 𝑣, as defined in Equation 3, and compare 𝑊 (𝑢, 𝑣) with the approximate count produced by the algebraic method. To better investigate the underlying characteristics of our method, we consider 4 scale-insensitive high-order
Type I Constraint
Type II Constraint
Iso Constraint
Guo et al.
Table 3: The impact of different orders. Order
Random
Min-Candidates
Max-Degree
Max-Candidates
Min-Degree
Heuristic
Yeast Human
1.23 2.23
1.24 4.64
1.44 9.28
1.47 9.61
1.27 2.88
1.21 2.13
Table 4: Effect of different modules.
(a) Yeast
(b) Human
Module
(-) Type-I
(-) Type-II
(-) Iso
Iso-only
Fastest
ASC
q-error Time (s)
6.8 ∗ 109
7.3 ∗ 1015
1.2 ∗ 106
2.1 ∗ 1016
7.5 ∗ 106
1.36
1.11
1.30
1.09
3.37
97.69 1.37
Figure 6: Correlation heatmap.
features of the constrained sequences: edge density (the number of internal edges/sequence length), indirect ratio (the number of indirect edges/internal edges), conflict intensity (the number of nested or intersecting cycles/indirect edge), and norm conflict (the number of nested or intersecting cycles/sequence length). The correlations are quantified using Spearman’s Rho as shown in Figure 6. As observed, the norm conflict and conflict intensity are the most significant factors of q-error on Yeast and Human, respectively. Since the two factors are highly correlated, the results suggest that the conflict pairs (i.e., nested or intersecting indirect edges within the sequence) are the main error source in algebraic counting. We also observe that the strongest correlation with q-error reaches around 0.5 on Yeast and 0.8 on Human. This is because our analysis focuses on the structural features of constrained sequences in the query graph, while the corresponding candidate space also substantially affects accuracy. Interestingly, it is observed that, for a sequence 𝑠 under 𝑢, even if we cannot guarantee the exact solution due to nested or indirect edges in 𝑠, the weights 𝑊 (𝑢, 𝑣) for some candidates 𝑣 may still achieve perfect accuracy with q-error = 1. This is attributed to the effort that we introduce a simple penalty based on precomputed reachable weight ratios along the reverse sequence, which can suppress the accumulation of invalid weights during propagation. Although this penalty cannot guarantee counting that exactly conforms to the matching states in the candidate space, it performs very well in practice, which is frequently observed in our experiments. Impact of Sequence Order. From the perspective of the FAQ framework [2], our constrained sequence defines a variable elimination order for the underlying subquery. Consequently, the resulting matrix chain multiplication can be viewed as a specialized execution plan for the associated sum-product query. We evaluate 6 constrained sequence ordering strategies, namely (1) Random ordering; (2) Min-Candidates, prioritizing query vertices with the smallest candidate set; (3) Max-Degree, prioritizing vertices with the highest degree; (4) Max-Candidates, prioritizing vertices with the largest candidate set; (5) Min-Degree, prioritizing vertices with the lowest degree; and (6) Heuristic, a greedy strategy that minimizes nested or intersecting indirect edges within the sequence. The average q-error of each 𝑊 (𝑢, 𝑣) over queries from Yeast and Human for each strategy is presented in Table 3. The results show that Heuristic performs the best, which aligns with our previous
analysis, i.e., fewer conflicting indirect edges bring less accumulated error. In addition, we observe that different ordering strategies consume comparable query time as expected. This is because the number of vertices and edges in the sequence remains unchanged across different orders. Analysis on Type-II module. For the Type-II module, when computing the LCA Ratio, we only account for the number of valid combinations in the LCA-endpoints triangle. For the two virtual edges between the LCA and the endpoints, we consider only reachability and ignore the actual number of paths in the candidate space for efficiency. In other words, regardless of how many tree paths connect an endpoint candidate to an LCA candidate, we treat them as a single path. Due to the infeasibility of exact algorithms, we evaluate the accuracy of the LCA Ratio for each LCA candidate. In particular, for each LCA candidate, we compute the variance in the number of paths from endpoint candidates for each non-tree edge to that LCA candidate, and analyze its correlation with estimation accuracy. The experiment results report that Spearman’s rho reaches 0.65 and 0.42 on Yeast and Human, respectively, both showing significant p-values, which confirms our hypothesis. Ablative Analysis. We conduct ablation studies on 3 modules for handling Type-I non-tree edge, Type-II non-tree edge, and Iso constraints. All variants follow the framework of candidate tree counting. (-) Type-I denotes the model without Type-I non-tree constraint computation, while (-) Type-II and (-) Iso have similar meanings. Iso-only removes both Type-I and Type-II non-tree constraint computations. We conduct experiments on 4 representative datasets, namely Yeast, Human, WordNet, and Eu2005, with different characteristics on average degree and label distributions. Accuracy. It is reported on Table 4 that ASC achieves the best accuracy performance with 97.69 q-error, whereas others perform much worse. For example, Iso-only generates an extremely high qerror of 2.1∗1016 . Surprisingly, the model without the Iso constraint module (i.e., (-) Iso) performs fairly well and even outperforms Fastest. This implies that our Type-I and Type-II modules can effectively capture structural constraints of the graph. In contrast, without Type-I and Type-II modules, the tree-homomorphism-based Iso module alone cannot achieve satisfactory accuracy. Query Time. As reported in Table 4, all variants consume almost the same amount of query time. This indicates that removing a specific module does not bring a clear efficiency improvement. The minor time differences are mainly attributed to the distribution of constraints in the query graph.
Algebraic Subgraph Counting
Constraint Analysis. To study how the constraint intensity of each module affects the counting performance, we vary the constraint intensity within the candidate tree counting framework. Specifically, for each module, we control the proportion of its constraints involved in the counting across 5 incremental levels: 0%, 25%, 50%, 75%, and 100%. For Type-I non-tree constraints, which arise in the non-tree edges of each constrained sequence in the candidate tree, an 𝑥% coverage rate signifies that only 𝑥% of these non-tree edges are considered during computation, while the left 2 modules are kept intact. Type-II constraints involve non-tree edges with different parent vertices, and we group these edges by the LCA vertex of their endpoints. Similarly, we restrict the number of nontree edges used in LCA Ratio computation within each LCA group by utilizing only 𝑥% of the edges, and the remaining 2 modules are kept intact. For iso constraints, which target groups of same-label vertices that may violate injectivity, we randomly select 𝑥% of the vertices in the group for constraint computation. We conduct experiments on Eu2005 and Human, which have high average degrees, and Yeast, which has a moderate average degree. Although we preset 5 coverage ratio levels, the actual constraint coverage often deviates significantly. Within the candidate tree counting framework, not every parent vertex in the query tree poses Type-I non-tree constraints in its child vertices. For the parent vertices without Type-I constraints, the weight updates degrade to Equation 2, resulting in exact computation. If only one parent vertex involves Type-I non-tree constraints while others remain exact, the overall accuracy impact remains limited, even if that specific parent ignores its constraints entirely. Consequently, evaluating performance based on preset coverage yields erratic, irregular fluctuations. To ensure a rigorous analysis, we utilize actual constraint coverage instead. Each vertex in the query tree is assigned a local constraint coverage, representing the actual coverage ratio among its child vertices. For vertices that have no children or lack non-tree edges between their children, the local coverage is set to 1.0. For all other vertices, we randomly select non-tree edges based on the preset coverage and calculate the resulting actual coverage. The global actual coverage for each query is then defined as the average of the local coverage across all its vertices. Similarly, we calculate the actual Type-II constraint coverage for each vertex. For Iso constraints, vertices are grouped by their labels. Groups that do not violate injectivity are assigned a local coverage of 1.0. For the remaining groups, the local coverage is derived from the vertices retention rate within the connected subgraph used for sampling. The overall actual coverage is then determined by averaging the coverage across all groups. Figure 7 illustrates the accuracy performance of all queries across the three datasets relative to the varying coverage of the three constraint types. The actual constraint coverage is determined by both the query graph structure and the preset ratios. Overall, for all three constraint types, the accuracy for most queries improves as the coverage increases. However, a small number of queries do not achieve optimal accuracy at 100% coverage. For the Type-I constraint, the computation is not perfectly exact. For constrained sequences containing nested cycles, we apply penalties to potential invalid weights based on pre-calculated reachable weight ratios. Depending on the actual matching states within the Cartesian space, the penalty can be either excessive or insufficient. In such cases,
removing certain non-tree edges may inadvertently mitigate these over- or under-penalizations, thereby leading to improved accuracy. Type-II constraint computation is performed at the LCA-group level. Within a single LCA group, the LCA ratio only considers whether a candidate of a non-tree edge endpoint is connected to the LCA candidate, while actually multiple valid tree paths may exist between them. Since our estimation does not capture the variation in the number of such paths, removing certain non-tree edges can occasionally compensate for part of the estimation error. In the Iso constraint sampling process, certain regions within the query subgraphs and their corresponding candidate space can be particularly difficult to sample. For instance, the local area of the query subgraph may exhibit high symmetry, resulting in a high density of homomorphisms but sparse isomorphisms within the candidate space. By limiting the constraint coverage, this symmetry may be broken, which can lead to a slight improvement in accuracy. In summary, while limiting constraint coverage may occasionally offset inherent model errors in some queries, a 100% coverage rate remains the optimal choice for most queries. Furthermore, as shown in Table 5, increasing the coverage from 0% to 100% results in negligible changes to the average query processing time.
7.5
Discussion
ASC achieves high accuracy in homomorphism counting as shown in Table 2, benefiting from its precise handling of structural constraints, especially the non-tree ones. While it also performs well for isomorphism counting, the iso ratio module remains the primary error source, as estimating the isomorphism ratio via sampling is inherently less accurate than the algebraic counting. We see two directions to improve it. First, the current sampling could be replaced by more robust alternatives, since sampling still faces inherent limitations: global sampling (e.g., in Fastest and Alley) often becomes inaccurate or fails in vast search spaces, and although our local sampling alleviates this by skipping query regions that cannot violate injectivity and sampling only over the conflict-prone regions, it does not fundamentally escape the accuracy limitations of sampling. Second, although the structural constraints, including non-tree edges, are already incorporated into the counting, their effect is currently confined to the weight accumulation in the counting stage. Another promising direction is thus to let such constraints act directly on the Cartesian-product space underlying the counting, pruning its enumerable states rather than only correcting the accumulated weights as in the current design.
8
Related Work
Subgraph Matching. Subgraph matching research has been dominated by the filtering-ordering-enumeration paradigm. Recent studies have enhanced this framework through efficient filtering [6, 7, 18–20, 22, 47], optimized ordering strategies [7, 9, 15, 46], and advanced enumeration techniques [50, 57]. Existing enumeration techniques primarily fall into two categories: exploration-based and decomposition-based methods. Within exploration-based approaches, performance is often improved via search space pruning [5, 19, 58], intermediate result reuse [26, 28, 29, 36, 37], and recursion reduction [27, 36, 49, 56]. Alternatively, decomposition-based
Guo et al.
Table 5: Time growth. Type Type-I Type-II Iso Type I Constraint
Type II Constraint
Cov. = 0.0
Cov. = 1.0
0.82 0.83 0.81
0.83 0.83 0.83
Iso Constraint
Figure 7: Accuracy against constraint coverage.
methods [1, 31, 32, 35, 38, 43, 51] partition the query graph and employ join techniques, such as binary joins [31, 32, 52] and worst-case optimal joins (WCOJ) [1, 38, 40, 51]. Subgraph Counting. Subgraph counting primarily employs three approximate paradigms. Summarization-based methods [16, 23, 39, 48] decompose queries and aggregate substructure counts, yet often struggle with accuracy due to independence assumptions. Sampling-based approaches [10, 13, 34, 47] are widely adopted but face inherent sampling failures, despite recent efforts in searchspace compression and strategy optimization. Finally, learningbased methods [18, 55, 59, 60] leverage GNNs to capture graph interactions, but their effectiveness is limited by the complex modeling between graph features and counts, high training overhead, and heavy dependence on ground-truth supervision.
9
Conclusion
In this paper, we study the problem of subgraph counting. We propose an algebraic counting method based on the candidate treebased counting framework. Using fast matrix computation operations, we can obtain accurate subgraph homomorphism counting with high efficiency. Based on this, we obtain the final subgraph isomorphism count via a local sampling strategy. Extensive experiments demonstrate the superior performance of our proposals.
Acknowledgements This work was supported by China Scholarship Council (20230844 0228), National Natural Science Foundation of China (62572137), Guangdong Basic and Applied Basic Research Foundation (2025A15 15011716), Major Key Project of PCL (PCL2024A05 and PCL2025A16), Australian Research Council Centre of Excellence for Mathematical Modelling of Cellular Systems (CE230100001), Australian Research Council Discovery Project (DP260100689), Australian Research Council Discovery Early Career Researcher Awards (DE250100226), and National Natural Science Foundation of China (U2241211).
References [1] Christopher R Aberger, Andrew Lamb, Susan Tu, Andres Nötzli, Kunle Olukotun, and Christopher Ré. 2017. Emptyheaded: A relational engine for graph processing. ACM Transactions on Database Systems (TODS) 42, 4 (2017), 1–44. [2] Mahmoud Abo Khamis, Hung Q Ngo, and Atri Rudra. 2016. FAQ: questions asked frequently. In Proceedings of the 35th ACM SIGMOD-SIGACT-SIGAI Symposium on Principles of Database Systems. 13–28. [3] Noga Alon, Raphael Yuster, and Uri Zwick. 1997. Finding and counting given length cycles. Algorithmica 17, 3 (1997), 209–223. [4] Renzo Angles, Marcelo Arenas, Pablo Barceló, Peter Boncz, George Fletcher, Claudio Gutierrez, Tobias Lindaaker, Marcus Paradies, Stefan Plantikow, Juan Sequeda, et al. 2018. G-CORE: A core for future graph query languages. In Proceedings of SIGMOD. 1421–1432. [5] Junya Arai, Yasuhiro Fujiwara, and Makoto Onizuka. 2023. Gup: Fast subgraph matching by guard-based pruning. Proceedings of SIGMOD 1, 2 (2023), 1–26. [6] Bibek Bhattarai, Hang Liu, and H. Howie Huang. 2019. CECI: Compact Embedding Cluster Index for Scalable Subgraph Matching. In Proceedings of SIGMOD. 1447–1462. [7] Fei Bi, Lijun Chang, Xuemin Lin, Lu Qin, and Wenjie Zhang. 2016. Efficient Subgraph Matching by Postponing Cartesian Products. In Proceedings of SIGMOD. 1199–1214. [8] Angela Bonifati, Wim Martens, and Thomas Timm. 2020. An analytical study of large SPARQL query logs. The VLDB Journal 29, 2 (2020), 655–679. [9] Vincenzo Bonnici, Rosalba Giugno, Alfredo Pulvirenti, Dennis Shasha, and Alfredo Ferro. 2013. A subgraph isomorphism algorithm and its application to biochemical data. BMC Bioinformatics (2013). [10] Marco Bressan, Stefano Leucci, and Alessandro Panconesi. 2019. Motivo: Fast Motif Counting via Succinct Color Coding and Adaptive Sampling. In Proc. VLDB Endow.. 1651–1663. [11] Walter Cai, Magdalena Balazinska, and Dan Suciu. 2019. Pessimistic Cardinality Estimation: Tighter Upper Bounds for Intermediate Join Cardinalities. In Proceedings of SIGMOD. 18–35. [12] Mario Cannataro and Pietro H Guzzi. 2012. Data management of protein interaction networks. Vol. 17. John Wiley & Sons. [13] Xiaowei Chen and John C. S. Lui. 2016. Mining Graphlet Counts in Online Social Networks. In Proceedings of ICDM. 71–80. [14] Yunyoung Choi, Kunsoo Park, and Hyunjoon Kim. 2023. Bice: Exploring compact search space by using bipartite matching and cell-wide verification. Proc. VLDB Endow. 16, 9 (2023), 2186–2198. [15] L.P. Cordella, P. Foggia, C. Sansone, and M. Vento. 2004. A (sub)graph isomorphism algorithm for matching large graphs. IEEE Transactions on Pattern Analysis and Machine Intelligence 26, 10 (2004), 1367–1372. [16] Kyle Deeds, Diandre Sabale, Moe Kayali, and Dan Suciu. 2025. Color: A framework for applying graph coloring to subgraph cardinality estimation. Proc. VLDB Endow. 2 (2025), 130–143. doi:10.14778/3705829.3705834 [17] Kyle B Deeds, Dan Suciu, and Magdalena Balazinska. 2023. Safebound: A practical system for generating cardinality bounds. Proceedings of SIGMOD 1, 1 (2023), 1–26. [18] Qiuyu Guo, Jianye Yang, Wenjie Zhang, Hanchen Wang, Ying Zhang, and Xuemin Lin. 2025. Efficient and Accurate Subgraph Counting: A Bottom-up Flow-learning Based Approach. Proc. VLDB Endow. 18, 8 (2025), 2695–2708. doi:10.14778/3742728. 3742758 [19] Myoungji Han, Hyunjoon Kim, Geonmo Gu, Kunsoo Park, and Wook-Shin Han. 2019. Efficient Subgraph Matching: Harmonizing Dynamic Programming, Adaptive Matching Order, and Failing Set Together. In Proceedings of SIGMOD. 1429–1446. [20] Wook-Shin Han, Jinsoo Lee, and Jeong-Hoon Lee. 2013. Turboiso: towards ultrafast and robust subgraph isomorphism search in large graph databases. In Proceedings of SIGMOD. 337–348. [21] Stephen Harris and Nigel Shadbolt. 2005. SPARQL query processing with conventional relational database systems. In International Conference on Web Information Systems Engineering. Springer, 235–244.
Algebraic Subgraph Counting
[22] Huahai He and Ambuj K. Singh. 2008. Graphs-at-a-Time: Query Language and Access Methods for Graph Databases. In Proceedings of SIGMOD. 405–418. [23] Tomaž Hočevar and Janez Demšar. 2014. A combinatorial approach to graphlet counting. Bioinformatics 30, 4 (2014), 559–565. [24] Wenzhe Hou, Xiang Zhao, and Bo Tang. 2024. LearnSC: An Efficient and Unified Learning-Based Framework for Subgraph Counting Problem. In Proceedings of ICDE. IEEE, 2625–2638. [25] Alon Itai and Michael Rodeh. 1977. Finding a minimum circuit in a graph. In Proceedings of the ninth annual ACM symposium on Theory of computing. 1–10. [26] Xun Jian, Zhiyuan Li, and Lei Chen. 2023. Suff: Accelerating subgraph matching with historical data. Proc. VLDB Endow. 16, 7 (2023), 1699–1711. [27] Tatiana Jin, Boyang Li, Yichao Li, Qihui Zhou, Qianli Ma, Yunjian Zhao, Hongzhi Chen, and James Cheng. 2023. Circinus: Fast redundancy-reduced subgraph matching. Proceedings of SIGMOD 1, 1 (2023), 1–26. [28] Hyunjoon Kim, Yunyoung Choi, Kunsoo Park, Xuemin Lin, Seok Hee Hong, and Wook Shin Han. 2021. Versatile Equivalences: Speeding up Subgraph Query Processing and Subgraph Matching. In Proceedings of SIGMOD. 925–937. [29] Hyunjoon Kim, Yunyoung Choi, Kunsoo Park, Xuemin Lin, Seok-Hee Hong, and Wook-Shin Han. 2023. Fast subgraph query processing and subgraph matching via static and dynamic equivalences. The VLDB journal 32, 2 (2023), 343–368. [30] Kyoungmin Kim, Hyeonji Kim, George Fletcher, and Wook Shin Han. 2021. Combining Sampling and Synopses with Worst-Case Optimal Runtime and Quality Guarantees for Graph Pattern Cardinality Estimation. In Proceedings of SIGMOD. 964–976. [31] Longbin Lai, Lu Qin, Xuemin Lin, and Lijun Chang. 2017. Scalable subgraph enumeration in mapreduce: a cost-oriented approach. The VLDB Journal 26, 3 (2017), 421–446. [32] Longbin Lai, Lu Qin, Xuemin Lin, Ying Zhang, Lijun Chang, and Shiyu Yang. 2016. Scalable distributed subgraph enumeration. Proc. VLDB Endow. 10, 3 (2016), 217–228. [33] John M Lewis. 1978. On the complexity of the maximum subgraph problem. In Proceedings of the tenth annual ACM symposium on Theory of computing. 265–274. [34] Feifei Li, Bin Wu, Ke Yi, and Zhuoyue Zhao. 2016. Wander Join: Online Aggregation via Random Walks. In Proceedings of SIGMOD. 615–629. [35] Qiyan Li, Jeffrey Xu Yu, and Zongyan He. 2025. Subgraph Matching: A New Decomposition Based Approach. Proc. VLDB Endow. 18, 11 (2025), 4282–4294. S [36] Yujie Lu, Zhijie Zhang, and Weiguo Zheng. 2025. B○X: Subgraph Matching with Batch Backtracking Search. Proceedings of SIGMOD 3, 1 (2025), 1–27. [37] Yujie Lu, Zhijie Zhang, Weiguo Zheng, and Lei Zou. 2025. Accelerating Subgraph Matching through Fine-grained and Powerful Equivalences. Proc. VLDB Endow. 18, 11 (2025), 3896–3909. [38] Amine Mhedhbi and Semih Salihoglu. 2019. Optimizing Subgraph Queries by Combining Binary and Worst-Case Optimal Joins. Proc. VLDB Endow. 12, 11 (2019). [39] Thomas Neumann and Guido Moerkotte. 2011. Characteristic sets: Accurate cardinality estimation for RDF queries with multiple joins. In Proceedings of ICDE. 984–994. [40] Hung Q. Ngo, Ely Porat, Christopher Ré, and Atri Rudra. 2012. Worst-Case Optimal Join Algorithms: [Extended Abstract]. In Proceedings of the ACM PODS. 37–48. [41] Himchan Park and Min-Soo Kim. 2018. EvoGraph: An effective and efficient graph upscaling method for preserving graph properties. In Proceedings of SIGKDD. 2051–2059. [42] N Pržulj, Derek G Corneil, and Igor Jurisica. 2006. Efficient estimation of graphlet frequency distributions in protein–protein interaction networks. Bioinformatics 22, 8 (2006), 974–980. [43] Miao Qiao, Hao Zhang, and Hong Cheng. 2017. Subgraph matching: on compression and computation. Proc. VLDB Endow. 11, 2 (2017), 176–188. [44] Xiafei Qiu, Wubin Cen, Zhengping Qian, You Peng, Ying Zhang, Xuemin Lin, and Jingren Zhou. 2018. Real-time constrained cycle detection in large dynamic graphs. Proceedings of the VLDB Endowment 11, 12 (2018), 1876–1888. [45] Siddhartha Sahu, Amine Mhedhbi, Semih Salihoglu, Jimmy Lin, and M. Tamer Özsu. 2017. The ubiquity of large graphs and surprising challenges of graph processing. In Proc. VLDB Endow.. 420–431. [46] Haichuan Shang, Ying Zhang, Xuemin Lin, and Jeffrey Xu Yu. 2008. Taming verification hardness: an efficient algorithm for testing subgraph isomorphism. In Proc. VLDB Endow., Vol. 1. VLDB Endowment, 364–375. [47] Wonseok Shin, Siwoo Song, Kunsoo Park, and Wook-Shin Han. 2024. Cardinality Estimation of Subgraph Matching: A Filtering-Sampling Approach. In Proc. VLDB Endow.. 1697–1709. [48] Giorgio Stefanoni, Boris Motik, and Egor V. Kostylev. 2018. Estimating the Cardinality of Conjunctive Queries over RDF Data Using Graph Summarisation. In Proceedings of WWW. 1043–1052. [49] Shixuan Sun, Yulin Che, Lipeng Wang, and Qiong Luo. 2019. Efficient parallel subgraph enumeration on a single machine. Proceedings of ICDE (2019), 232–243. [50] Shixuan Sun and Qiong Luo. 2020. In-Memory Subgraph Matching: An In-depth Study. In Proceedings of SIGMOD. 1083–1098.
[51] Shixuan Sun, Xibo Sun, Yulin Che, Qiong Luo, and Bingsheng He. 2020. Rapidmatch: A holistic approach to subgraph query processing. In Proc. VLDB Endow., Vol. 14. 176–188. [52] Zhao Sun, Hongzhi Wang, Haixun Wang, Bin Shao, and Jianzhong Li. 2012. Efficient Subgraph Matching on Billion Node Graphs. Proc. VLDB Endow. 5, 9 (2012). [53] J. R. Ullmann. 1976. An Algorithm for Subgraph Isomorphism. J. ACM (1976), 31–42. [54] Vladimir Vacic, Lilia M Iakoucheva, Stefano Lonardi, and Predrag Radivojac. 2010. Graphlet kernels for prediction of functional residues in protein structures. Journal of Computational Biology 17, 1 (2010), 55–72. [55] Hanchen Wang, Rong Hu, Ying Zhang, Lu Qin, Wei Wang, and Wenjie Zhang. 2022. Neural Subgraph Counting with Wasserstein Estimator. In Proceedings of SIGMOD. 160–175. [56] Rongjian Yang, Zhijie Zhang, Weiguo Zheng, and Jeffrey Xu Yu. 2023. Fast continuous subgraph matching over streaming graphs via backtracking reduction. Proceedings of SIGMOD 1, 1 (2023), 1–26. [57] Zhijie Zhang, Yujie Lu, Weiguo Zheng, and Xuemin Lin. 2024. A Comprehensive Survey and Experimental Study of Subgraph Matching: Trends, Unbiasedness, and Interaction. In Proceedings of SIGMOD, Vol. 2. ACM New York, NY, USA, 1–29. [58] Zhijie Zhang and Weiguo Zheng. 2025. BEE: Towards Redundancy Reduction via Block-Separator Decomposition for Subgraph Matching. Proceedings of SIGMOD 3, 4 (2025), 1–27. [59] Kangfei Zhao, Jeffrey Xu Yu, Qiyan Li, Hao Zhang, and Yu Rong. 2023. Learned sketch for subgraph counting: a holistic approach. The VLDB Journal 32, 5 (2023), 937–962. [60] Kangfei Zhao, Jeffrey Xu Yu, Hao Zhang, Qiyan Li, and Yu Rong. 2021. A Learned Sketch for Subgraph Counting. In Proceedings of SIGMOD. 2142–2155.