⚡ Bolt: [performance improvement] 배열 할당 병목 제거 및 행렬 곱 최적화

.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.
+
+## 2026-07-30 - NumPy 2.0 boolean array astype copy error
+**Learning:** Changing a `bool` dtype to `float64` using `astype(..., copy=False)` works in NumPy 1.x by silently falling back to a copy. However, NumPy 2.0+ strictly enforces the `copy=False` argument and throws a `ValueError` because changing boolean bits to a 64-bit float physically requires a new memory allocation.
+**Action:** When casting boolean arrays to float arrays for matrix multiplication, simply use `.astype(dtype)` without specifying `copy=False`, or use `copy=None` to allow copying when necessary and prevent hard crashes in NumPy 2.0+.

python/fast_mlsirm/estimators/marginal.py

@@ -1385,7 +1385,10 @@ def estep(current_params):
     stopping_tolerance = float(tol * (1.0 + abs(ll)))
     for it in range(1, max_iter + 1):
         for i in range(n_items):
-            r = np.stack([post[y[:, i] == k].sum(axis=0) for k in range(k_cat)], axis=1)
+            # Optimized M-step aggregation: cast boolean mask and use matrix multiplication
+            # to entirely skip intermediate allocations across persons and categories.
+            mask = y[:, i, None] == np.arange(k_cat)
+            r = post.T @ mask.astype(post.dtype)
             params[i] = _gpcm_m_step_item(params[i], nodes, r)
         next_ll, post = estep(params)
         if not np.isfinite(next_ll):

python/fast_mlsirm/objective.py

@@ -133,7 +133,9 @@ 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
+        # Optimized alpha gradient: use dense matrix multiplication and advanced indexing
+        # to avoid O(N*J) intermediate array allocation when multiplying errors by theta.
+        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