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 ~20n2 to ~4n2 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(n2), 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.

Factorization at scale[](#factorization_at_scale)

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 20n2 bytes to 4n2 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(n2) 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 temporal setup[](#the_temporal_setup)

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.

Choosing the rank k[](#choosing_the_rank_k)

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