Block-Preconditioned PageRank for Centrality and Intervention Analysis In Commercial Networks

Summary

This protocol computes PageRank-based centrality and evaluates targeted interventions in commercial networks by integrating principal component analysis–based directional gain, multi-source edge weighting, and a block-preconditioned Krylov-subspace linear solver to improve numerical stability and reproducibility.

Abstract

Commercial networks are shaped by heterogeneous drivers (e.g., geography, category structure, and operating performance), so a single financial indicator may not reflect a node’s structural role in the overall system. This protocol extracts low-dimensional feature directions using principal component analysis and applies a paired-spectrum consistency check based on a structured symplectic block matrix to improve the reliability of the selected directions. Using these directions, the workflow constructs a weighted directed graph that fuses spatial similarity, business-category synergy, and a directional-gain term derived from feature gradients. Centrality is then computed by solving the PageRank linear system with a block-preconditioned generalized minimal residual method, with explicit convergence and diagnostic checkpoints reported for reproducibility. The workflow is demonstrated on state-level retail aggregates (49 nodes) from a public dataset, where a small number of principal components captures most of the feature variance and supports stable directional weighting. Finally, the protocol evaluates a targeted strong-to-weak intervention by reallocating a quantified fraction of edge weight from high-centrality nodes to low-centrality nodes and recomputing PageRank under the same personalization setting. Overall, this protocol enables users to build an interpretable multi-source commercial network, compute numerically verified PageRank centrality, and test intervention policies with clearly defined diagnostics.

Introduction

Commercial networks (e.g., state-level retail systems and regional business districts) consist of heterogeneous actors that both cooperate and compete. Their interactions are shaped by spatial location, category structure, and directional movement of customers and goods. As a result, single-point indicators (e.g., sales alone) may not objectively reflect a node’s structural role in the network^1,2,3. Firms within the same district can differ markedly in scale and operating conditions, and their interactions are further conditioned by spatial context and temporal footfall. Therefore, heuristic indicators and simple sales rankings can be inadequate for characterizing network position and cross-node reinforcement or inhibition³.

Complex-network theory provides a holistic perspective for quantifying node importance. PageRank (PR) estimates importance by simulating a random walk and has been widely used in ranking and diffusion models. Prior work shows that PR can be generalized to weighted settings by replacing the adjacency matrix with a weight matrix and by using node strength rather than degree⁴. Weighted PR formulations further balance degree- and strength-based contributions through tunable parameters and have been applied to large-scale economic dependency networks⁵. Mobility- and spending-derived place-to-place dependencies have also been used to construct behavior-based networks that link centrality patterns to urban economic resilience⁶. Recent studies integrating multi-source urban data highlight that single-indicator rankings often fail to represent multi-factor mechanisms, motivating network constructions that jointly model structure and attributes⁷. Evidence at business-district granularity similarly suggests that quantity, category diversity, and category structure jointly shape commercial vitality and should be modeled together². Related analyses connecting centrality to facility patterns and reachability further reinforce the need to capture both relational structure and contextual features in commercial settings^8˒9.

Methodologically, the literature points to two practical needs that are often under-addressed in applied commercial-network studies. First, when multivariate attributes are used to construct directional or feature-informed weights, the stability of the extracted directions should be explicitly checked rather than assumed⁷. Second, centrality computation should report reproducible diagnostic checkpoints (e.g., sparsity/connectivity checks and solver convergence) so results are verifiable beyond a single ranking output^4˒5.

Despite this progress, many applied studies still report rankings without (i) an explicit stability check for directional information extracted from multivariate features, (ii) reproducible checkpoints that verify graph sparsity/connectivity and numerical convergence, and (iii) an intervention evaluation that is parameterized and comparable across settings^1,2,3. This protocol addresses these gaps by making the full workflow auditable: it extracts low-dimensional feature directions, verifies direction consistency, constructs a multi-source weighted directed graph, and computes PR centrality using a numerically stable block-preconditioned iterative solver with clearly reported diagnostics. It further evaluates a targeted strong-to-weak intervention under a fixed personalization setting to quantify how interventions comparably redistribute centrality.

Compared with PR variants that only reweight edges or add layers, this workflow provides measurable outputs that users can directly verify and compare: (i) an interpretable directional contribution at the edge level, (ii) explicit intermediate diagnostics (variance-explained targets, sparsity/connectivity checks, and solver convergence flags), and (iii) an intervention-effect metric that quantifies changes in PR for specified source/target sets under a fixed personalization setting. This protocol is most appropriate when nodes represent spatially located commercial units (e.g., states, cities, districts, platforms) with coordinates and multivariate operating features, and when the feature space exhibits a stable low-dimensional structure rather than extreme sparsity or low-signal noise.

Protocol

1. Data preparation and normalization

1. Import and filter the dataset. Import a public retail orders dataset and retain records from 2014–2017¹⁸.
2. Define nodes. Define the node set V as U.S. states and retain states with sufficient observations so that |V|=49. Let n=|V|.
3. Construct node-level feature vectors. For each node i\in V, aggregate transaction records and compute a 16-dimensional feature vector x_i∈R¹⁶. The feature vector includes sales statistics (total, mean, standard deviation, and range), profit statistics (total, mean, standard deviation, and range), quantity statistics (total, mean, standard deviation), discount statistics (mean and standard deviation), and three additional indicators consisting of profit margin (ProfitMargin), order count (OrderCount), and average profit per order (AvgOrderProfit).
4. Deterministic handling of missing values and numerical stabilizer. Replace any missing derived feature entry with 0 and record the number of replacements. Define a fixed stabilizer δ=10^-12 for denominators that may approach zero in later computations (used in Equation 4).
5. Standardize features. Apply Z-score normalization across nodes to obtain the standardized feature matrix X _std∈R^49×16.
   CHECKPOINT 1: Verify that X_{\text{std}} is 49\times 16 with no missing entries after standardization and confirm that the subsequent class counts sum to 49.
6. Define sales classes for stratification. Partition nodes into three classes (Low/Medium/High) by tertiles of total sales and record class sizes (Low: 16; Medium: 17; High: 16)¹⁸. Denote class membership by c(i) ∈ {L,M,H}.

2. JRS-validated PCA dimensionality reduction (JRS-PCA)

1. Covariance estimation. Compute the covariance matrix C from X_{\text{std}} as follows.
   C = cov(X_std)  (1) 
2. Principal component extraction. Compute the first k eigenpairs of C using sparse eigendecomposition and set k = 4 by default.
3. Variance explained. Record the variance explained by each principal component and the cumulative variance explained by PCs 1–4².
4. PC score representation. For each node i, compute the PC score vector z_i\in\mathbb{R}^k and record (z_i,1,z_i,2) for visualization in a PC1–PC2 scatter colored by sales class.
5. JRS block construction. Construct the JRS-structured symplectic block matrix as follows.
   S = diag(C,-C)  (2) 
6. Paired-spectrum consistency with explicit matching. Compute eigenvalues of S, partition them into positive and negative sets, sort the positive eigenvalues in descending order and sort the absolute values of the negative eigenvalues in descending order, pair them by index, and compute the maximum pairing deviation \max_j|\lambda^+_j+\lambda^-_j|. This matching rule fixes the implementation detail for reproducibility².
   CHECKPOINT 2: Report the cumulative variance explained by PCs 1–4 and report  to confirm a stable low-dimensional structure prior to directional gain construction².

3. Construction of a multi-source weighted directed graph

1. Spatial similarity weights. For each ordered pair (i,j), compute the spatial Gaussian kernel weight as follows:
      (3) 
   Here p_i = (x_i,y_i)  denotes the centroid coordinate of node i and σ is the bandwidth, set to 10.0 by default.
2. Class-based synergy weights. Define the synergy baseline w_syn(i,j) using the sales-class relation from Step 1.6.
   Set w_syn(i,j)=0.6 when c(i) = c(j),setw_syn(i,j)=0.4 when classes are adjacent (Low–Medium or Medium–High), and set w_syn(i,j) = 0.2  when classes are non-adjacent (Low–High).
3. PCA directional gain. Compute the stabilized unit direction in PC score space by
      (4) 
   with δ=10^-12 . Define the nonnegative alignment score along PC1 using the reference direction r = (1,0,…,0)^⊤ as
      (5) 
   Define the same-class penalty as
      (6) 
   and compute the directional-gain weight by
   w_ori(i,j) = β⋅w_syn(i,j) ⋅ g(i,j) ⋅ π(i,j)   (7) 
   with β=1.0. The penalty π(i,j)  implements a mild same-class down-weighting to encourage cross-class links while preserving within-class connectivity⁵.
4. Weight fusion and self-loop removal. Fuse multi-source weights to form the directed weight matrix
      (8) 
   and set diag(W_total) = 0 to remove self-loops.
5. Symmetric inspection matrix (optional). Compute the symmetric proxy matrix
      (9) 
   for inspection only.
6. Sparsification with quantitative validation. Retain the top q% of nonzero entries in W_total(defaultq=35%) to obtain the sparse directed adjacency matrix A_dir. Record nnz(A_dir) and density², and optionally apply the same sparsification rule toW_sym to obtain A_sym for inspection².
   CHECKPOINT 3: After sparsification, report nnz(A_dir) and density, report the number of weakly connected components (target: 1), and check whether any column sums are zero. If zero-sum columns exist, apply dangling-node handling in Step 4.1². If sparsification disconnects the graph, relax the sparsification stringency sequentially from 35% to 30% and then to 25%, repeating the checkpoint after each adjustment.

4. Category-block preconditioned LGMRES PageRank solution

1. Column-stochastic transition matrix with dangling-node handling. Compute column sums of A_dir and form the diagonal matrix D. For any zero-sum column, replace that column with the personalization vector u before normalization, and form the column-stochastic transition matrix by
   W_c = A_dirD^-1  (10) 
2. PageRank as a linear system. Define PageRank with damping factor α(defaultα=0.85) and personalization vector u by
   p = (1-α)u + αW_c^⊤p   (11) 
   and solve the equivalent linear system
   (I-( I -αW_c^⊤)p = (1 - α)u  (12) 
   (II-withM = I - αW_c^⊤ andb = (1-α)u. 
3. Category-block right preconditioner. Extract diagonal subblocks of M corresponding to Low/Medium/High classes and assemble a block-diagonal approximation P^{-1} using pseudoinverses for each diagonal subblock. For indices not covered by a block, use the reciprocal of the diagonal entry of M as a diagonal fallback¹⁹.
4. Iterative solution and diagnostic reporting. Solve Mp = b using a limited-memory generalized minimal residual solver with right preconditioner P^-1 , configured as rtol = 1e−6, atol = 0, maxiter = 500, inner_m = 30, outer_k = 3¹⁹. Record the solver info flag (target: 0), the final relative residual ‖Mp - b‖₂/‖b‖₂ , and the outer iteration counts¹⁹.
5. Normalization and validity checks. If small negative values occur, clip them to 0 and record the maximum clipped magnitude, then renormalize so that ∑_ip_i = 1 .
   CHECKPOINT 4: Report info, report ‖Mp - b‖₂/‖b‖₂ , report , and report ∑_ip_i as reproducibility diagnostics rather than qualitative statements¹⁹.

5. Intervention experiment and result comparison

1. Source and target selection. Using a fixed personalization vector u, define the top three nodes by baseline p as the source set S and the bottom three nodes as the target set T, and record the selected nodes².
2. Intervention strength with a cap constraint. Choose ε  in[0.05,0.15] such that the injected weight from each source node does not exceed 10% of that node’s original outgoing weight sum after sparsification, and record ε along with the injected-to-original ratio per source node².
3. Deterministic edge-weight update. Compute the mean positive edge weight  in A_dir
   , then update each ordered pair (s,t) with s ∈ Sandt ∈ T  by
   A'_dir(s,t) = A_dir(s,t) + εw  (13) 
   and recompute W_c from A'_{\text{dir}} as in Step 4.1 while keeping α  and u unchanged.
4. Recompute PageRank and quantify change. Repeat Steps 4.2–4.5 to obtain p', compute Δp=p'-p , and report absolute and relative changes for nodes in S and T². Verify that solver diagnostics remain within targets (info = 0; residual ≤ rtol) after intervention².
   CHECKPOINT 5: Report \varepsilon, report injected-to-original ratio per source node, report solver info and residual before and after intervention, and report a Δ  p summary for S and T².

Results

Experimental conditions and data

All analyses were run on a desktop computer with a Windows operating system. The workflow was implemented in Python 3.12 using NumPy, Pandas, SciPy, and Matplotlib, and the random seed was fixed at 2025 to support reproducibility¹⁸. A public retail orders dataset was filtered to retain records from 2014–2017²². Orders were aggregated to the state level to define 49 nodes (U.S. states with sufficie...

Discussion

This protocol provides a reproducible workflow for assessing centrality and evaluating targeted interventions in state-level commercial networks by integrating multi-source weighted graph construction, PCA-derived directional gain, and a category-block preconditioned LGMRES PageRank solver^18˒23. Rather than presenting centrality as a single ranking output, the workflow foregrounds intermediate diagnostics—variance-explained checks, sparsity/connectivit...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Matplotlib	Matplotlib Development Team		v3.8+ Figure generation and visualization of centrality/intervention results
NumPy	NumPy Developers		v1.26+ Numerical array operations and matrix computation
Pandas	Pandas Development Team		v2.2+ Data cleaning, aggregation, and tabular data management
Python	Python Software Foundation		v3.12 Core programming environment for data processing and model execution
SciPy	SciPy Community		v1.13+ Sparse linear algebra and iterative solver support (e.g., LGMRES)