[PERF] Memory Profiling Baseline for O(N^3) EDoF Calculations

[hritikkumarpradhan]

Description

While preparing my GSoC 2026 proposal focused on optimizing pyGAM's matrix operations, I profiled the GAM._pirls optimization loop and isolated a severe $O(N^3)$ memory bottleneck.

When fitting overparameterized models (where the number of splines/features $p$ exceeds the number of samples $n$), the Effective Degrees of Freedom (EDoF) calculation creates a massive intermediate dense matrix that scales quadratically with $p$, causing massive RAM spikes and fatal MemoryError crashes on standard hardware.

The Bottleneck

The leak originates during the hat matrix trace extraction. The current codebase calculates this via cross-multiplication:
edof = np.diagonal(U1.dot(U1.T))

If $p = 10,000$, U1.dot(U1.T) forces the allocation of a $10,000 \times 10,000$ float64 matrix before immediately discarding everything except the diagonal.

Memory Profiling Proof

I wrote a synthetic benchmark ($n=500, p=10000$) to isolate this exact line using memory_profiler.

As shown below, the legacy operation triggers an immediate ~760 MiB memory spike. I have also profiled my proposed GSoC fix, which replaces the dense matrix creation with a vectorized $O(N)$ tensor-sum: (U1**2).sum(axis=1). This drops the overhead to near zero while preserving mathematical equivalence.

--- Running Legacy Calculation ---
Line #    Mem usage    Increment  Occurrences   Line Contents
=============================================================
    22     76.4 MiB      0.0 MiB           1       start_time = time.time()
    24                                             # The bottleneck:
    25    839.5 MiB    763.1 MiB           1       edof = np.diagonal(u1_matrix.dot(u1_matrix.T))
    27    839.5 MiB      0.0 MiB           1       end_time = time.time()

--- Running Proposed Vectorized Fix ---
Line #    Mem usage    Increment  Occurrences   Line Contents
=============================================================
    38     76.5 MiB      0.0 MiB           1       start_time = time.time()
    40                                             # The fix:
    41    114.6 MiB     38.1 MiB           1       edof = (u1_matrix**2).sum(axis=1)
    43    114.6 MiB      0.0 MiB           1       end_time = time.time()

Next Steps

I have already contributed this specific stress-test as an EDoFBenchmark class to the ASV automated benchmarking suite currently being built in PR #503.

I am opening this issue to formally track the algorithmic bottleneck itself. I plan to submit the vectorized (U1**2).sum(axis=1) fix as a core part of my GSoC 2026 project deliverables.

[hritikkumarpradhan]

Hi @dswah, following up from Discord. I have opened this issue to formally document and track the (N^3) memory bottleneck located in GAM._pirls (pygam.py:1033).

As discussed, my current PR against Karan's ASV suite (KaranSinghDev/pyGAM#1 -> #503) only benchmarks the legacy implementation to establish a formal baseline. I am now updating that PR to include time_optimized_edof() and peakmem_optimized_edof() so we can track the O(N) improvement side-by-side.

Once the benchmarking infrastructure is in place, I'll submit the actual vectorized fix to pygam.py to resolve this issue as part of my GSoC 2026 deliverables.

[KaranSinghDev]

@hritikkumarpradhan Thanks for the great contribution, I wnated to know whether you are going to add the optimized part in a seperate test file or in the same one which I merged ? I was also working on it, so Let me know how are you planning to proceed with it, looking forward to it . Also I suggest to keep in mind to use less size and mention the estimate spike the operation may generate (the current one in which the memory spike was about 800 Mib is fine) .

[hritikkumarpradhan]

Hey @KaranSinghDev! Thanks for the kind words, and great work refactoring the EDoF benchmark into #503 — locking the dimensions as constants and adding the random seed was exactly the right call.

I'm planning to add the optimized methods (time_optimized_edof / peakmem_optimized_edof) within the same bench_edof.py file, keeping both legacy and optimized benchmarks in the same class for a clean ASV side-by-side comparison. I'll keep the matrix at (10000, 500) to preserve the ~800 MiB legacy spike vs ~38 MiB optimized — that contrast is too compelling to change.

Since you were also working on this, let's sync before I open the follow-up PR to avoid any overlap. What were you thinking of adding?

[KaranSinghDev]

@hritikkumarpradhan Ok, I was working on putting the optimized solution as well so that we can show both in the test in the benchmark but since it was proposed by you first, so I left it out.