NAICS-aware GraphSAGE for POI Prediction

• The paper introduces a GraphSAGE-based model that integrates learnable NAICS embeddings with multi-modal features for edge-level regression of POI co-visitation counts.
• Graph construction decomposes nationwide data into state-level, sparse spatio-temporal graphs while incorporating engineered spatial, temporal, and socioeconomic attributes for robust prediction.
• Empirical results show substantial performance gains (R² improving from 0.243 to 0.625) over baselines, highlighting the critical role of business taxonomy in modeling visitor behavior.

NAICS-aware GraphSAGE is a graph neural network for population-level prediction of point-of-interest (POI) co-visitation, introduced in “NAICS-Aware Graph Neural Networks for Large-Scale POI Co-visitation Prediction: A Multi-Modal Dataset and Methodology” (Alrubyli et al., 25 Jul 2025). It formulates co-visitation forecasting as edge-level regression on a spatio-temporal graph whose nodes are brand-level businesses and whose targets are future counts of distinct mobile devices visiting both brands within a one-hour window. Its defining characteristic is the integration of detailed business taxonomy through learnable embeddings of 6-digit NAICS codes, combined with spatial, temporal, and socioeconomic features in an end-to-end GraphSAGE-based framework. The method is designed for sparse, large-scale graphs and is evaluated on a multi-state dataset with tens of millions of co-visitation observations.

1. Problem formulation and predictive target

The task is population-level co-visitation prediction between business brands, cast as an edge-regression problem on a spatio-temporal graph (Alrubyli et al., 25 Jul 2025). For a brand pair $(i,j)$ and prediction horizon $t+\Delta$, the target is the count of distinct mobile devices that visited both brands within a one-hour time window during month $t+\Delta$:

$$y_{ij}^{t+\Delta} = \text{number of distinct devices that co-visited } i \text{ and } j \text{ in month } t+\Delta \text{ within a one-hour window.}$$

The forecasting objective is expressed as

$$\hat{y}_{ij}^{t+\Delta} = f\bigl(\mathcal{G}_{t-L:t}, \mathbf{X}_{t-L:t}, \mathbf{S}_t; \Theta \bigr),$$

where $$\mathcal{G}_{t-L:t} = (\mathcal{V}, \mathcal{E}_{t-L} \cup \cdots \cup \mathcal{E}_t)$$, $\mathbf{X}_{t-L:t}$ denotes node and edge feature sequences, $\mathbf{S}_t$ denotes socioeconomic context, and $\Theta$ are learnable parameters.

A notable modeling choice is that no log-transform or normalization is applied to the target $y$. The counts are described as skewed and sparse, yet mean squared error is minimized directly on raw counts. This distinguishes the method from approaches that normalize or transform highly skewed visitation volumes before training. A plausible implication is that the model is intended to preserve direct interpretability of predicted counts at inference time, even though this may accentuate head–tail imbalance.

The motivating premise is that geographic distance alone does not explain heterogeneity in co-visitation behavior. In the paper’s framing, business semantics matter: coffee shops, fine dining restaurants, fitness centers, banks, entertainment venues, and specialty food stores can induce distinct patterns of sequential or near-synchronous movement even when located close together. NAICS-aware GraphSAGE operationalizes that claim by treating industry taxonomy as a core component of node representation rather than as auxiliary metadata.

2. Graph construction, decomposition, and dataset scope

The national network is decomposed into 48 state-level graphs $t+\Delta$0 in order to preserve regional patterns and ensure scalability (Alrubyli et al., 25 Jul 2025). Each node represents a brand within a state, with multiple physical locations of the same brand aggregated into a single brand-level node. Edges are undirected and are included only when the monthly co-visitation count exceeds a minimum threshold of 5 device traces. The paper states that this threshold is intended to ensure statistical significance while controlling graph density.

The resulting graphs are extremely sparse, with densities below $t+\Delta$1. Scale varies substantially across states: large states such as Texas reach approximately 1.3 million edges, while small states such as Vermont have approximately 27 thousand edges; Vermont is described as having fewer than 2 thousand nodes. The observed co-visit counts are used as supervision targets on edges, but the convolutional layers themselves do not use explicit edge weights in message passing. Instead, mean aggregation is applied uniformly, and spatial, temporal, and socioeconomic information enters through engineered features used in the prediction head.

The underlying dataset, termed POI-Graph, is assembled from two proprietary datasets across the continental United States spanning January 2018 to June 2022, although modeling focuses on the pre-pandemic interval from January 2018 to March 2020. The POI monthly data contains approximately 286 million records. After cleaning, including removal of 0.03% missing NAICS annotations and assignment of the most frequent 6-digit NAICS code per brand, the initial coverage includes 236,814 brands and 1,020 NAICS categories. The modeling subset uses 92,486 brands across 48 states with 276 unique 6-digit NAICS codes represented as embeddings.

The co-visitation source consists of 251.5 million weekly records of $t+\Delta$2, where counts reflect devices co-visiting within one hour. These records are aggregated to the monthly level for stability and preservation of seasonality. Cleaning removes malformed timestamps, excludes post-March-2020 records, left-joins with POI data, and filters outliers such as anomalous monthly co-visitation counts above 40,000. The final dataset comprises 94.9 million monthly co-visit observations, 45.3 million unique edges, 92,486 brands, and 48 states.

Temporal splitting is chronological rather than random: training uses January 2018 to December 2019, validation uses January to February 2020, and testing uses March 2020. This design reflects the forecasting orientation of the task and avoids leakage across time.

3. Feature engineering and NAICS taxonomy integration

Node features are defined as

$t+\Delta$3

where $t+\Delta$4 is a learned embedding for the brand’s 6-digit NAICS code and $t+\Delta$5 is a state- and industry-stratified quantile-based popularity score (Alrubyli et al., 25 Jul 2025). The embedding table is $t+\Delta$6 with $t+\Delta$7 in the final dataset. The first vector, at index 0, is fixed to zeros; all other vectors use PyTorch defaults for initialization.

The model does not explicitly encode the hierarchical structure of NAICS at the 2-, 3-, 4-, 5-, and 6-digit levels. There are no hierarchical priors, regularizers, or parameter-sharing mechanisms across levels. Instead, the method relies on end-to-end training to learn semantic proximity among 6-digit categories from co-visitation supervision and graph context. This is an important clarification because the term “NAICS-aware” could be misconstrued as implying direct hierarchical modeling; in this method, taxonomy awareness is realized through learnable 6-digit embeddings only.

Edge features are 10-dimensional. Spatial distance contributes one dimension, standardized after log-transform:

$t+\Delta$8

Temporal seasonality contributes two dimensions via cyclical month encoding:

$t+\Delta$9

Popularity interactions contribute seven dimensions:

$t+\Delta$0

Socioeconomic context is encoded at the state level as $t+\Delta$1 using ACS/BEA indicators with lag-1 integration to reflect delayed impacts. Edge-level extended features are then formed as

$t+\Delta$2

The empirical role of NAICS information is especially prominent. In the ablation study, removing NAICS embeddings decreases $t+\Delta$3 from 0.625 to 0.258, an absolute drop of 0.367, which is the largest degradation among the tested components; replacing them with random embeddings yields $t+\Delta$4 (Alrubyli et al., 25 Jul 2025). The paper also reports that t-SNE visualizations of learned embeddings show meaningful industry clusters. This suggests that end-to-end supervision is sufficient to induce an economically interpretable embedding space even without explicit NAICS hierarchy constraints.

4. Model architecture and edge-level regression mechanism

The backbone is a 5-layer GraphSAGE with mean aggregation, hidden dimension $t+\Delta$5 at every layer, $t+\Delta$6 activation, and dropout 0.2 per layer (Alrubyli et al., 25 Jul 2025). Initialization is

$t+\Delta$7

At layer $t+\Delta$8, the neighborhood message and update are

$t+\Delta$9

$$y_{ij}^{t+\Delta} = \text{number of distinct devices that co-visited } i \text{ and } j \text{ in month } t+\Delta \text{ within a one-hour window.}$$0

with mean aggregation defined by

$$y_{ij}^{t+\Delta} = \text{number of distinct devices that co-visited } i \text{ and } j \text{ in month } t+\Delta \text{ within a one-hour window.}$$1

After $$y_{ij}^{t+\Delta} = \text{number of distinct devices that co-visited } i \text{ and } j \text{ in month } t+\Delta \text{ within a one-hour window.}$$2 layers, the node representation is $$y_{ij}^{t+\Delta} = \text{number of distinct devices that co-visited } i \text{ and } j \text{ in month } t+\Delta \text{ within a one-hour window.}$$3.

The architecture separates node representation learning from edge-context conditioning. Node-pair features are projected as

$$y_{ij}^{t+\Delta} = \text{number of distinct devices that co-visited } i \text{ and } j \text{ in month } t+\Delta \text{ within a one-hour window.}$$4

where $$y_{ij}^{t+\Delta} = \text{number of distinct devices that co-visited } i \text{ and } j \text{ in month } t+\Delta \text{ within a one-hour window.}$$5. Extended edge-context features are projected as

$$y_{ij}^{t+\Delta} = \text{number of distinct devices that co-visited } i \text{ and } j \text{ in month } t+\Delta \text{ within a one-hour window.}$$6

where $$y_{ij}^{t+\Delta} = \text{number of distinct devices that co-visited } i \text{ and } j \text{ in month } t+\Delta \text{ within a one-hour window.}$$7. Final prediction uses a linear layer over the concatenated 288-dimensional fused representation:

$$y_{ij}^{t+\Delta} = \text{number of distinct devices that co-visited } i \text{ and } j \text{ in month } t+\Delta \text{ within a one-hour window.}$$8

The stated rationale for this two-stage fusion is dimensionality reduction from 1072 to 288 dimensions, separation of node and edge-context transformations, and improved generalization and efficiency. The graph convolution itself does not consume edge features; those features influence predictions only through the head. This design choice matters conceptually: it keeps message passing inductive and scalable while allowing edge-level regression to remain sensitive to distance, seasonality, popularity interactions, and state socioeconomic conditions.

Neighborhood sampling is used for scalability and inductive training, with decreasing fanout $$y_{ij}^{t+\Delta} = \text{number of distinct devices that co-visited } i \text{ and } j \text{ in month } t+\Delta \text{ within a one-hour window.}$$9. The paper states that this stabilizes memory usage and reduces noise from distant weak neighbors while preserving broader context in earlier layers. Empirically, this achieves approximately 94% of full-neighborhood performance with approximately 75% lower memory. A plausible implication is that the architecture is tuned not only for accuracy but also for practical deployability on million-edge state graphs.

5. Optimization, sparsity handling, and computational scaling

Training is performed as edge-level regression on counts with mini-batch mean squared error loss:

$$\hat{y}_{ij}^{t+\Delta} = f\bigl(\mathcal{G}_{t-L:t}, \mathbf{X}_{t-L:t}, \mathbf{S}_t; \Theta \bigr),$$0

Because graph density is below $$\hat{y}_{ij}^{t+\Delta} = f\bigl(\mathcal{G}_{t-L:t}, \mathbf{X}_{t-L:t}, \mathbf{S}_t; \Theta \bigr),$$1, the training procedure uses balanced mini-batches containing equal numbers of positives, meaning observed co-visits, and negatives, meaning randomly sampled zero co-visit pairs (Alrubyli et al., 25 Jul 2025). This is the principal mechanism for handling extreme sparsity in supervision.

Optimization uses AdamW with learning rate $$\hat{y}_{ij}^{t+\Delta} = f\bigl(\mathcal{G}_{t-L:t}, \mathbf{X}_{t-L:t}, \mathbf{S}_t; \Theta \bigr),$$2, weight decay $$\hat{y}_{ij}^{t+\Delta} = f\bigl(\mathcal{G}_{t-L:t}, \mathbf{X}_{t-L:t}, \mathbf{S}_t; \Theta \bigr),$$3, and $$\hat{y}_{ij}^{t+\Delta} = f\bigl(\mathcal{G}_{t-L:t}, \mathbf{X}_{t-L:t}, \mathbf{S}_t; \Theta \bigr),$$4. Additional regularization and scheduling choices are dropout 0.2 per layer, early stopping on validation MAE with patience 20, and learning-rate halving after 10 epochs of plateau. Mixed precision (AMP) and gradient accumulation are supported and used when needed.

The implementation hyperparameters reported for the selected model are a 5-layer GNN, hidden size 512, $$\hat{y}_{ij}^{t+\Delta} = f\bigl(\mathcal{G}_{t-L:t}, \mathbf{X}_{t-L:t}, \mathbf{S}_t; \Theta \bigr),$$5 activations, dropout 0.2, mean aggregation, and neighbor sampling fanout $$\hat{y}_{ij}^{t+\Delta} = f\bigl(\mathcal{G}_{t-L:t}, \mathbf{X}_{t-L:t}, \mathbf{S}_t; \Theta \bigr),$$6. The fusion head uses $$\hat{y}_{ij}^{t+\Delta} = f\bigl(\mathcal{G}_{t-L:t}, \mathbf{X}_{t-L:t}, \mathbf{S}_t; \Theta \bigr),$$7, $$\hat{y}_{ij}^{t+\Delta} = f\bigl(\mathcal{G}_{t-L:t}, \mathbf{X}_{t-L:t}, \mathbf{S}_t; \Theta \bigr),$$8, ReLU activations, and a final linear layer over 288 dimensions. Model selection is based on grid search over hidden sizes $$\hat{y}_{ij}^{t+\Delta} = f\bigl(\mathcal{G}_{t-L:t}, \mathbf{X}_{t-L:t}, \mathbf{S}_t; \Theta \bigr),$$9, layer counts $$\mathcal{G}_{t-L:t} = (\mathcal{V}, \mathcal{E}_{t-L} \cup \cdots \cup \mathcal{E}_t)$$0, learning rates $$\mathcal{G}_{t-L:t} = (\mathcal{V}, \mathcal{E}_{t-L} \cup \cdots \cup \mathcal{E}_t)$$1, dropout rates $$\mathcal{G}_{t-L:t} = (\mathcal{V}, \mathcal{E}_{t-L} \cup \cdots \cup \mathcal{E}_t)$$2, and NAICS embedding dimensions $$\mathcal{G}_{t-L:t} = (\mathcal{V}, \mathcal{E}_{t-L} \cup \cdots \cup \mathcal{E}_t)$$3, with 16 selected as best. Baselines are stated to be tuned comparably.

The computational analysis gives per-layer complexity $$\mathcal{G}_{t-L:t} = (\mathcal{V}, \mathcal{E}_{t-L} \cup \cdots \cup \mathcal{E}_t)$$4 with $$\mathcal{G}_{t-L:t} = (\mathcal{V}, \mathcal{E}_{t-L} \cup \cdots \cup \mathcal{E}_t)$$5, and total per-epoch complexity $$\mathcal{G}_{t-L:t} = (\mathcal{V}, \mathcal{E}_{t-L} \cup \cdots \cup \mathcal{E}_t)$$6 with $$\mathcal{G}_{t-L:t} = (\mathcal{V}, \mathcal{E}_{t-L} \cup \cdots \cup \mathcal{E}_t)$$7 and $$\mathcal{G}_{t-L:t} = (\mathcal{V}, \mathcal{E}_{t-L} \cup \cdots \cup \mathcal{E}_t)$$8 the average sampled neighborhood factor. Memory complexity is $$\mathcal{G}_{t-L:t} = (\mathcal{V}, \mathcal{E}_{t-L} \cup \cdots \cup \mathcal{E}_t)$$9. State-wise decomposition enables training per-state models or batched state subgraphs. The fanout schedule controls degree explosion and enables graphs with more than one million edges to run on a single GPU.

The reported throughput is approximately 25,000 edge predictions per second after embedding computation. Training time per epoch ranges from 12 minutes for Vermont-sized graphs with 27 thousand edges to 3.2 hours for Texas-sized graphs with 1.3 million edges. Peak memory reaches up to 28 GB on the largest state. Hardware descriptions vary by section: the experiments report NVIDIA RTX 3090 (24GB), PyTorch 2.0, and PyG 2.3; implementation details report NVIDIA A100 (40GB); the appendix reports NVIDIA L40 (40GB). The paper’s description indicates that the pipeline runs across these GPU classes, with A100 and L40 used for large-scale runs and RTX 3090 for reported benchmark experiments.

6. Empirical performance, interpretation, and limitations

Evaluation uses MAE, RMSE, MSE, and $\mathbf{X}_{t-L:t}$0 for regression accuracy, together with NDCG@10 and MRR for ranking quality. Results are averaged over 5 seeds with 95% confidence intervals, and paired $\mathbf{X}_{t-L:t}$1-tests with Bonferroni correction and Cohen’s $\mathbf{X}_{t-L:t}$2 are reported (Alrubyli et al., 25 Jul 2025).

The held-out test results identify STHGCN as the strongest baseline among those reported, but NAICS-aware GraphSAGE substantially outperforms it.

Model	Key regression result	Key ranking result
STHGCN	$\mathbf{X}_{t-L:t}$3	NDCG@10 $\mathbf{X}_{t-L:t}$4, MRR $\mathbf{X}_{t-L:t}$5
NAICS-aware GraphSAGE	MAE $\mathbf{X}_{t-L:t}$6, RMSE $\mathbf{X}_{t-L:t}$7, MSE $\mathbf{X}_{t-L:t}$8, $\mathbf{X}_{t-L:t}$9	NDCG@10 $\mathbf{S}_t$0, MRR $\mathbf{S}_t$1

The paper summarizes the relative improvements over STHGCN as $\mathbf{S}_t$2 for a 157% relative gain, NDCG@10: $\mathbf{S}_t$3 for a 31–32% gain, and MRR: $\mathbf{S}_t$4 for a 24–25% gain. Statistical significance is reported as $\mathbf{S}_t$5, with very large effect sizes; for $\mathbf{S}_t$6, Cohen’s $\mathbf{S}_t$7 is approximately 25.4. Simpler baselines perform poorly by comparison: Gravity has $\mathbf{S}_t$8, GeoMF $\mathbf{S}_t$9, LightGBM $\Theta$0, GAT $\Theta$1, and GCN $\Theta$2.

Geographic breakdown indicates strong performance across several large states, with California at MAE 4.87 and $\Theta$3 0.698, Texas at MAE 5.01 and $\Theta$4 0.672, New York at MAE 5.34 and $\Theta$5 0.645, and Florida at MAE 5.28 and $\Theta$6 0.651. Lower but still positive performance is reported for Wyoming at MAE 6.12 and $\Theta$7 0.523 and Vermont at MAE 6.08 and $\Theta$8 0.534. This aligns with the stated limitation that rural states show lower $\Theta$9 due to data sparsity.

The ablation study is central to the model’s interpretation. Removing NAICS embeddings reduces $y$0 to 0.258; removing socioeconomic features reduces it to 0.511; using 3 layers instead of 5 reduces it to 0.460; reducing hidden dimension from 512 to 256 reduces it to 0.402; removing temporal encoding reduces it to 0.448; and removing popularity features reduces it to 0.594. These results support the paper’s argument that business semantics matter beyond proximity and that multi-modal fusion contributes materially to predictive performance.

Error analysis indicates that the model is well-calibrated over five orders of magnitude in a predicted-versus-actual log-log view, though it slightly underestimates head edges. The appendix reports best accuracy in mid-frequency ranges. The paper suggests that specialized losses could improve tail behavior. Failure modes include head–tail imbalance, holiday anomalies, higher variance in professional services under NAICS 54, and lower performance in rural states. The scope is explicitly limited to pre-pandemic, US-only modeling and to 6-digit NAICS semantics without explicit hierarchy.

Interpretability is illustrated qualitatively through case studies. In Austin, a coffee chain is reported to have higher co-visitation with fitness centers (NAICS 713940) than with banks (NAICS 522110), beyond equal proximity and foot traffic. In Denver mixed-use planning, strong complementarities are reported between restaurants (NAICS 722511) and entertainment (NAICS 711190), with weaker association to automotive services (NAICS 811111). Specialty food stores (NAICS 445299) are described as exhibiting dependency on restaurants, with potential relevance for policy targeting during disruptions. These examples are consistent with the model’s central premise that semantic business relationships complement spatial structure in co-visitation prediction.

Future directions listed in the paper include temporal graph networks for finer seasonality and drift, focal or quantile losses for skewed targets, transfer learning for data-sparse regions, international validation, federated learning for privacy, and faster real-time inference strategies. This suggests that NAICS-aware GraphSAGE is best understood not as a complete solution to POI interaction modeling, but as a scalable framework that demonstrates the predictive value of integrating business taxonomy into inductive graph learning.