diff --git a/.jules/bolt.md b/.jules/bolt.md index 73e3fbaf9..8fc2021f8 100644 --- a/.jules/bolt.md +++ b/.jules/bolt.md @@ -33,3 +33,7 @@ ## 2025-05-19 - Dot product scalar gradients allocation **Learning:** During gradient calculation, `float((e * (-gamma * distance)).sum())` creates two full-size `(N, J)` arrays: one for the scaled distance and one for the element-wise multiplication before reduction. **Action:** Replace `(A * B).sum()` with `np.vdot(A, B)` when scalar reduction is needed over matrix multiplication (where `B` can incorporate scalars naturally like `-gamma * np.vdot(A, B)`). This entirely avoids the 2D array allocation overhead and yields order-of-magnitude improvements in scalar gradient components. + +## 2024-05-18 - Optimize grad_alpha computation via dot product +**Learning:** In NumPy, combining an element-wise multiplication with an axis sum (`(e * theta[:, factors]).sum(axis=0)`) creates a massive intermediate array allocation (N x J). For large datasets, this becomes a major memory and performance bottleneck. +**Action:** Replace element-wise multiplication and `.sum(axis=0)` reductions on 2D arrays with dense matrix multiplications and advanced indexing (`(e.T @ theta)[np.arange(J), factors]`). This computes the exact same result while drastically reducing memory allocation and leveraging highly optimized BLAS routines. diff --git a/python/fast_mlsirm/objective.py b/python/fast_mlsirm/objective.py index e43df9809..a3c429142 100644 --- a/python/fast_mlsirm/objective.py +++ b/python/fast_mlsirm/objective.py @@ -133,7 +133,7 @@ def neg_loglik_and_grad( grad_b = e.sum(axis=0) grad_alpha = np.zeros_like(params.alpha) if free_alpha: - grad_alpha = (e * params.theta[:, factors]).sum(axis=0) * a + grad_alpha = (e.T @ params.theta)[np.arange(e.shape[1]), factors] * a # Optimized gradient computation: replace loop over dimensions with matrix multiplication # We embed 'a' directly into the projection matrix to avoid a JxD intermediate array allocation during multiplication