GPU-Accelerated Clustering for Financial Instruments at Scale

Quant strategies routinely group instruments for portfolio construction, risk aggregation, statistical arbitrage, and trade surveillance. Incorrect groupings can make concentrated positions appear diversified, obscure risk shared across nominal boundaries, and select statistical-arbitrage pairs whose relationships fail under stress.

The practical difficulty is that the right groupings are neither directly observable nor stable. Factor exposures drift, instruments change classifications, and dependencies can change sharply during market stress. A clustering pipeline must therefore separate routine variation from structural change and be inexpensive enough to rerun as new returns arrive.

There are two common ways to group instruments from a dependence matrix. Hard clustering methods are computationally cheap but assign every instrument to exactly one group, which breaks down at sector boundaries and masks the graded exposures that matter for risk budgeting. Soft factorization methods like SymNMF handle boundary instruments and produce usable factor loadings, but their dense matrix objectives have historically limited practical use to moderate instrument counts rather than the scale at which this problem actually lives.

This post covers a workflow that addresses both limitations. The workflow starts with rolling return windows and constructs two complementary inputs: absolute Pearson correlation for broad co-movement and the tail pairwise dependence matrix (TPDM) for joint behavior during extreme observations. SymNMF represents each instrument through a row of nonnegative factor loadings. Retaining the row provides a soft representation; taking its argmax produces a hard label.

A memory-efficient SymNMF formulation reduces peak storage from ~20n² to ~4n² bytes, which is what makes ~100,000 instruments fit on a single NVIDIA GB200. For larger problems, a distributed implementation row-shards the dependence matrix and reduces communication to O(nk) rather than O(n²), enabling factorization of one million instruments across 16 nodes. A single adaptive solver, AdaptGrow, handles both correlation and tail-dependence inputs by reading the eigenspectrum to choose between full-batch and block-stochastic gradients, eliminating the need to select or tune separate solvers for different input structures.

The result is a clustering pipeline that produces hard labels, soft factor loadings, and structural-break signals, reruns cheaply as new returns arrive, and scales from a single GPU to multi-node infrastructure without changing the solver interface.

A companion notebook, linked below, implements the full pipeline and reproduces all results in this post.

Scale is constrained first by memory. A dense FP32 dependence matrix requires ~40 GB for 100,000 instruments and ~4 TB for 1 million instruments. A naive SymNMF implementation also materializes several additional n x n intermediates. The trace-based formulation used here eliminates those intermediates, reducing estimated peak storage from approximately 20n² bytes to 4n² bytes plus smaller factor buffers. This change is what makes approximately 100,000 instruments fit on one high-memory GPU.

NVIDIA acceleration enters at each stage of the pipeline. PyTorch dispatches the dominant SH matrix multiplications to cuBLAS. cuSOLVER performs the spectral probe used for rank and solver selection. cuDF keeps optional Parquet ingestion and preprocessing on the GPU. For scale-out, PyTorch Distributed row-shards S while keeping a replica of H on each worker. NCCL all-gathers the row-sharded S H products and all-reduces the gradients, so communication operates on O(nk) data rather than the full O(n²) matrix. The environment is packaged with an NVIDIA NGC PyTorch container and cudf-cu13.

In the companion paper, the 100,000-instrument matrix was distributed across four NVIDIA GB200 GPUs for faster execution, although its 40 GB input fits on one GB200. Across three seeds in FP32, AdaptGrow converged in 13.0 seconds on correlation and 12.4 seconds on TPDM. At 1 million instruments, the 4 TB matrix was row-sharded across 64 GB200 GPUs on 16 nodes; full-batch AdaGrad completed the correlation factorization in approximately 2 minutes, while AdaptGrow completed the TPDM factorization in approximately 4 minutes. These are individual factorization measurements, not end-to-end timings for all 250 temporal windows.

The workflow evaluates 250 rolling windows, approximating daily re-clustering over one trading year. The synthetic return stream contains two controlled events: instruments changing their planted group membership and several groups experiencing a joint tail-stress episode without changing membership.

This controlled setup verifies two different behaviors. Adjusted Rand index (ARI) should identify the membership change, while TPDM should expose the co-crash that ordinary correlation largely misses. For production use, replace the synthetic generator with a returns table while preserving the same windowing, dependence-estimation, factorization, and monitoring stages.

The million-instrument results are separate distributed scale tests and require infrastructure comparable to the published 16-node configuration.

Start by inspecting the leading eigenvalues from an initial representative window. Choose k at the clearest separation between signal eigenvalues and the noise floor, then keep k fixed across subsequent windows so that stability scores remain comparable. The synthetic data used here has a planted rank of 24. Production data may not contain a sharp gap, so rank selection should also be checked against cluster interpretability and stability.

Figure 1. The correlation and TPDM spectra both recover the planted rank k=24 from the eigenvalue gap

For each rolling window, the workflow passes the dependence matrix S and selected rank k to AdaptGrow. The solver returns H, where each row of H contains an instrument's soft factor loadings, and taking the row-wise argmax produces a hard cluster label.

Each window uses the same fixed-seed initialization and is fitted independently rather than warm-started, preventing previous labels from masking a genuine reclassification. AdaptGrow is a single adaptive solver that auto-configures from the matrix's eigenvalue spectrum. Both regimes use the same per-coordinate AdaGrad preconditioner and differ only in how the gradient is computed. AdaptGrow seeds its batch fraction from the post-rank eigenvalue gap described in more detail below.

A clean gap selects the full gradient. A flatter post-rank spectrum starts with a cheaper block-sampled gradient corrected using Stochastic Variance Reduced Gradient (SVRG), then expands the sample toward the full gradient if progress stalls. The same solver runs unchanged from one GPU to many.

where g is either the full gradient ∇f(H) (full-batch branch) or a block-sampled estimate of it (stochastic branch), with all operations above being element-wise.

After selecting the factorization rank k, AdaptGrow examines the post-rank gap ratio \(\gamma_{k+1} = \lambda_{k+1} / |\lambda_{k+2}|\). The extra eigenvalue allows for a common factor, such as a market factor alongside the group-level factors. In the companion experiments, a large post-rank gap corresponded to shorter correlation-matrix runs, for which full-batch AdaGrad was most efficient. A flatter spectrum corresponded to longer TPDM runs, where lower-cost block-sampled gradients were more useful. AdaptGrow therefore starts with sampled SVRG updates when the gap is small and increases the sampled fraction toward the full matrix when progress stalls. The threshold of five used here is an empirical setting from these experiments, not a universal statistical cutoff. AdaptGrow algorithm sketch:

When only one label per instrument is required, spherical k-means provides a lower-cost baseline. It clusters the L2-normalized rows of S using cosine similarity, making it the appropriate comparison for