perf(objective): compute linear_predictor distances with in-place NumPy operations

.jules/bolt.md

@@ -33,3 +33,6 @@
 ## 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-07-18 - NumPy In-Place Operations for Distance Calculation
+**Learning:** In the `fast_mlsirm` distance calculation, creating multiple intermediate arrays of size `(N, J)` during arithmetic operations (addition, `np.maximum`, `np.sqrt`) is a significant bottleneck. Using `dist_sq += ...` and `out=` kwargs (e.g. `np.sqrt(dist_sq, out=dist_sq)`) reduces memory overhead and improves performance drastically for large matrices.
+**Action:** Always prefer in-place NumPy operations (like `+=`, `-=`, and `out=`) when calculating large pairwise metrics if it is safe to overwrite the array, ensuring memory efficiency and faster execution times without introducing regressions.

python/fast_mlsirm/objective.py

@@ -97,11 +97,12 @@ def linear_predictor(
         interaction = np.dot(params.xi, params.zeta.T)
     elif uses_space:
         # Optimized distance computation: replace O(N*J*D) 3D broadcast with O(N*J) 2D dot product
-        xi_sq = np.einsum('ij,ij->i', params.xi, params.xi)
-        zeta_sq = np.einsum('ij,ij->i', params.zeta, params.zeta)
-        dist_sq = xi_sq[:, None] + zeta_sq[None, :] - 2 * np.dot(params.xi, params.zeta.T)
-        dist_sq = np.maximum(dist_sq, 0.0)
-        distance = np.sqrt(dist_sq + eps_distance)
+        dist_sq = -2.0 * np.dot(params.xi, params.zeta.T)
+        dist_sq += np.einsum('ij,ij->i', params.xi, params.xi)[:, None]
+        dist_sq += np.einsum('ij,ij->i', params.zeta, params.zeta)[None, :]
+        np.maximum(dist_sq, 0.0, out=dist_sq)
+        dist_sq += eps_distance
+        distance = np.sqrt(dist_sq, out=dist_sq)
         interaction = -params.gamma * distance
     else:
         distance = np.zeros((params.theta.shape[0], len(factor_id)), dtype=np.float64)