Memory-Efficient Estimation of Large Correlation Matrices

DS3294 Data Science Practice, Team 9

Project Overview

This project studies how memory and runtime scale when computing Pearson correlation matrices as the number of data series (N) grows. We benchmark the naive in-memory approach (np.corrcoef), identify the practical memory wall where computation becomes infeasible, and implement four alternative strategies to work around it.

The Core Problem

Computing a Pearson correlation matrix for N time series produces an NxN output.

• Memory: O(N^2), an NxN float64 matrix needs 8 * N^2 bytes
• Runtime: O(N^2 * T), dominated by the matrix multiply

At N=10,000 the matrix alone is 763 MB. At N=50,000 it's 18.6 GB, which exceeds most machines' RAM. This is the memory wall.

Data Source

We use synthetically generated data, matrices of shape (N, T) filled with iid standard normal random numbers via np.random.default_rng(seed).standard_normal(). A fixed seed (42) ensures reproducibility.

Why synthetic: scaling behaviour depends on N and T, not on actual data values. Synthetic data lets us test N from 10 to 10,000+ without sourcing external datasets.

Pre-processing

Minimal pre-processing since data is synthetic:

1. Data is generated as an (N x T) float64 array, one row per time series
2. np.corrcoef handles centering and normalization internally
3. For manual validation, we z-score each row: subtract mean, divide by std (ddof=1)

Repository Structure

├── CoreProblem/                          # Baseline study
│   ├── main.py                           # Benchmarks, memory wall, correctness check
│   ├── utils.py                          # generate_data(), measure_correlation()
│   ├── plotting.py                       # Scaling and breakdown plots
│   ├── scaling_analysis.py               # Log-log slope fitting, O(N^2) verification
│   └── method_comparison.py              # Side-by-side comparison of all approaches
│
├── AdditionalMilestones/
│   ├── Blockwise-study/                  # Tiled NxN computation
│   │   ├── blockUtils.py                 # measure_correlation_blockwise()
│   │   ├── main.py                       # Benchmark + tile size sensitivity
│   │   └── plotting.py
│   │
│   ├── Streaming-study/                  # Chunked along sample axis
│   │   ├── streamUtils.py                # measure_correlation_streaming()
│   │   ├── main.py                       # Benchmark vs baseline at all N
│   │   └── plotting.py
│   │
│   ├── standardisation/                  # On-the-fly z-scoring
│   │   ├── standutils.py                 # correlation_on_the_fly()
│   │   ├── standardization.py            # Benchmark runner
│   │   └── plotting.py
│   │
│   ├── correlation_approximation/        # Random pair sampling
│   │   ├── corrutils.py                  # approximate_correlation_sampling()
│   │   ├── correlation_approximation.py  # Benchmark runner
│   │   └── plotting.py
│   │
│   ├── GarbageCollection-study/          # GC impact on measurement
│   │   ├── defUtils.py, main.py, plotting.py
│   │
│   ├── PSUTIL-study/                     # RSS-based memory measurement
│   │   ├── psutils.py, main.py, plotting.py
│   │
│   └── MemoryEstimationDashboard-UI/     # Early dashboard prototype
│
├── results/                              # Pre-computed benchmark plots
│   ├── *.png                             # Core problem plots
│   ├── Blockwise-study/                  # Blockwise plots
│   ├── Streaming-study/                  # Streaming plots
│   ├── Standardisation-study/            # On-the-fly plots
│   ├── Approximate-study/                # Sampling plots
│   ├── GarbageCollection-study/          # GC study plots
│   └── PSUTIL-study/                     # psutil study plots
│
├── Reports_Docs_etc/
│   ├── brief_report.txt                  # Summary of findings
│   ├── report/                           # LaTeX report + compiled PDF
│   └── DiscussionAndWorkflow.md
│
├── UI_testing.py                         # Interactive Streamlit dashboard
├── correlation_UI_README.md              # Dashboard documentation
├── Checklist.md                          # Task completion status (all done)
└── README.md                             # This file

Methods Studied

1. Baseline (np.corrcoef)

Standard NumPy. Computes the full NxN matrix in one call using optimised BLAS. Fastest option, but internally allocates 3 to 4 NxN temporary arrays (centered data, covariance, normalisation), so peak memory is roughly 3 to 4x the output matrix size.

2. Blockwise (Tiled)

Z-scores all rows once, then fills the NxN output tile by tile. Only one small tile is in memory at a time besides the output. Reduces peak memory from roughly 3 to 4x down to 1 to 2x the NxN matrix size. Slightly slower than baseline due to the Python tile loop, but more robust at scale.

3. Streaming (Two-Pass Chunked)

Processes data in chunks along the sample (column) axis using a two-pass approach:

• Pass 1: Computes the global row means by streaming through chunks.
• Pass 2: Centers each chunk relative to the global mean and accumulates the covariance matrix as M2 += centered @ centered.T.

After all chunks are processed, M2 is divided by (T-1) and normalised to a correlation matrix in place, avoiding extra NxN copies.

Important note on memory measurement: tracemalloc (our primary measurement tool) tracks Python-level memory allocations only. When we benchmark streaming vs baseline, tracemalloc reports higher peak memory for streaming. This is not because streaming uses more memory in practice; it is because np.corrcoef does its heavy lifting in C/BLAS (invisible to tracemalloc), while the streaming method's Python-level NxN accumulator is fully tracked. The psutil study (which measures total process RSS including C allocations) confirms that the real picture is more nuanced. This is itself an interesting finding about the limitations of Python-level memory profiling.

4. On-the-fly Standardisation

Z-scores each row only when needed for a dot product, then discards it. Eliminates the NxT z-score matrix entirely. Peak memory is NxN (output) plus two row buffers of length T. The tradeoff is a significant runtime penalty due to a Python-level double loop (O(N^2) iterations in Python, each doing O(T) numpy work).

5. Approximate (Random Sampling)

Instead of computing all N*(N-1)/2 pairwise correlations, randomly samples K=10,000 pairs. Memory is O(K), completely independent of N. Gives the distribution of correlations (mean, std, histogram) but not the full matrix.

Summary Table

Method	Speed	RAM Usage	Max N (16 GB)
Baseline (np.corrcoef)	Fastest	Extreme (~3-4x NxN)	~25,000
Blockwise (Tiled)	Medium	High (~1-2x NxN)	~30,000
Streaming (Chunked)	Medium	High (NxN accumulator)	~30,000
On-the-fly	Very slow	Moderate (NxN only)	~30,000
Approximate (K pairs)	Fast	Minimal (O(K))	Unlimited

Additionally, we study GC impact (with vs without gc.collect()) and psutil vs tracemalloc for memory measurement methodology.

Interactive Dashboard

A Streamlit dashboard (UI_testing.py) provides:

• Theoretical analysis: enter your RAM and see your memory wall
• Pre-computed results: browse benchmark plots from all 7 studies without waiting
• Live benchmarks: run methods on your own hardware at any N (with RAM safety guard)
• Correctness validation: verify each method against the NumPy baseline

pip install streamlit numpy matplotlib pandas psutil
streamlit run UI_testing.py

See correlation_UI_README.md for detailed documentation.

How to Run

Core study

cd CoreProblem
python main.py

Individual additional studies

cd AdditionalMilestones/Blockwise-study && python main.py
cd AdditionalMilestones/Streaming-study && python main.py
cd AdditionalMilestones/standardisation && python standardization.py
cd AdditionalMilestones/correlation_approximation && python correlation_approximation.py

Scaling analysis and method comparison

cd CoreProblem
python scaling_analysis.py
python method_comparison.py

Requirements: Python 3.8+, NumPy, Matplotlib, Pandas, psutil

Testing and Correctness

• Baseline validation: np.corrcoef output matches manual z-score computation to ~1e-16 (floating-point limit), confirmed via np.allclose
• All alternative methods validated against baseline at small N, all pass np.allclose
• Diagonal check: all correlation matrices have exactly 1.0 on the diagonal
• Memory measurement: tracemalloc with gc.collect() before each run gives the most reliable isolation of computation cost

Key Findings

1. Memory scales as O(N^2): log-log slope ~2.0 for N >= 500, matching theory
2. Runtime scales as O(N^2): exponent ~2.04 for large N
3. np.corrcoef overhead: measured peak is ~5% above theoretical minimum at N=10,000 due to internal temporaries
4. Blockwise reduces peak memory by eliminating intermediate NxN arrays, extending feasible N
5. Streaming processes data in chunks along the sample axis using a two-pass approach, reducing intermediate memory; the NxN accumulator is still required, so it does not eliminate the quadratic wall
6. Approximate sampling is the only method that completely bypasses O(N^2) memory, running in O(K) regardless of N
7. On-the-fly standardisation saves the NxT z-score matrix but is impractically slow due to Python loops
8. tracemalloc vs psutil: tracemalloc only sees Python-level allocations, missing C/BLAS internals. This makes np.corrcoef appear cheaper than it actually is, and makes Python-based alternatives appear relatively more expensive. Both tools have their place.

See Reports_Docs_etc/brief_report.txt for the full summary and Reports_Docs_etc/report/report.pdf for the complete LaTeX report.

Future Work

• Applying these techniques to real-world datasets (financial time series, LLM embedding spaces)
• Investigating lower-complexity algorithms (random projections) to reduce the O(N^2) bottleneck
• Non-stationarity handling for real-world time series (see Reports_Docs_etc/NonStationarityExplained.md)

Authors

1. Aditya Raj Chaudhary
2. Ananthakrishna S V
3. Sriram Kintada
4. Vaibhav Dhekawat