⚡ Bolt: optimize covariance and RANSAC shuffles

[suraj-ranganath]

⚡ Bolt: I've implemented two key performance optimizations that significantly speed up core numerical workflows.

1. Covariance Matrix Operations

I identified that diag_nd and matrix power functions like cov_logm were bottlenecks in iterative algorithms like Riemannian mean estimation.

• Optimization: Replaced loop-based diagonal matrix creation with vectorized advanced indexing and shifted to broadcasting for eigenvalue scaling.
• Impact: Robust covariance mean estimation (cov_mean) is now ~12.5% faster.

2. Fast Fisher-Yates Shuffles

The rand_permutation and rand_sample functions are critical for ICA and RANSAC parity but had high constant-time overhead.

• Optimization: Vectorized random number generation and replaced the generic round_mat with math.floor(x + 0.5) for fast scalar rounding in the tight shuffle loop.
• Impact: Random permutation of 1M elements is now ~30% faster (from 1.43s to 0.99s).

Verified with existing parity tests (test_utils_covariance.py, test_utils_ransac.py, test_parity_rng.py). Documented learnings in .jules/bolt.md.

---

PR created automatically by Jules for task 346178192807701607 started by @suraj-ranganath

[google-labs-jules]

👋 Jules, reporting for duty! I'm here to lend a hand with this pull request.

When you start a review, I'll add a 👀 emoji to each comment to let you know I've read it. I'll focus on feedback directed at me and will do my best to stay out of conversations between you and other bots or reviewers to keep the noise down.

I'll push a commit with your requested changes shortly after. Please note there might be a delay between these steps, but rest assured I'm on the job!

For more direct control, you can switch me to Reactive Mode. When this mode is on, I will only act on comments where you specifically mention me with @jules. You can find this option in the Pull Request section of your global Jules UI settings. You can always switch back!

New to Jules? Learn more at jules.google/docs.

---

For security, I will only act on instructions from the user who triggered this task.

[claude]

Claude finished @suraj-ranganath's task in 2m 48s —— View job

---

Code review

• Overall assessment: Safe to merge.
• Highest-risk area: MATLAB RNG parity in ransac.py (verified preserved).
• Merge recommendation: Safe to merge.

I focused on whether these optimizations preserve numerical behavior and MATLAB parity. They do.

Blocking

None.

Important

None.

Nits

• .jules/bolt.md (new file) — This is Jules' internal learnings log, not EEGPrep-owned documentation, and it's being committed into the repo (.jules/ isn't in .gitignore). AGENTS.md asks to touch only files you must and avoid clutter. Consider dropping it from the PR or adding .jules/ to .gitignore. Non-blocking.

Correctness / parity (verified)

• ransac.py batched RNG — Switching from per-iteration stream.rand() to a single stream.rand(m) / stream.rand(n-1) preserves the sequence: NumPy's legacy RandomState (MT19937) fills an array by drawing consecutive doubles in the same order as scalar calls. The index mapping is also correct — in rand_permutation, random_vals[n-1-k] for k going n-1 → 1 consumes random_vals[0], random_vals[1], …, matching the original scalar-draw order.
• round_mat → math.floor(x + 0.5) — Exact for these inputs. round_mat(x, decimals=0) is literally copysign(floor(abs(x)+0.5), x), which equals floor(x+0.5) for non-negative x. Both (remaining-1)*rand() and k*rand() are always ≥ 0, so this is bit-identical, not just approximate. Index bounds (idx = k + choice ≤ n-1) are unchanged.
• covariance.py broadcasting — V @ diag(f(D)) @ V.T rewritten as (V * f(D)[..., newaxis, :]) @ V.T is the standard column-scaling identity and is correct: f(D)[..., newaxis, :] has shape (..., 1, N), scaling column j of V by f(D)_j. The finite_matmul robustness wrapper is preserved.
• diag_nd rewrite — The advanced-indexing version (res[..., idx, idx] = M) is correct for the 1D/2D/3D cases and still handles the remaining caller at covariance.py:211 (target='diag' shrinkage) and its tests in test_utils_covariance.py.

Test gaps

None required. test_utils_covariance.py, test_utils_ransac.py, and test_parity_rng.py already cover the changed behavior and directly guard the MATLAB-parity contract these edits touch. (Note: I was unable to execute the suite here — the sandbox blocked test/interpreter execution — but the changes are numerically equivalent by construction as detailed above.)

EEGLAB parity notes

None. The MATLAB rand_permutation/rand_sample parity contract and Riemannian covariance behavior are preserved.