From 53fac6c5ba06f96f9058126e519b41df3861a51a Mon Sep 17 00:00:00 2001
From: "google-labs-jules[bot]"
 <161369871+google-labs-jules[bot]@users.noreply.github.com>
Date: Fri, 12 Jun 2026 21:09:28 +0000
Subject: [PATCH] perf: vectorize lattice fourier transform using quadratic
 form

Co-authored-by: makskliczkowski <48489493+makskliczkowski@users.noreply.github.com>
---
 .jules/bolt.md                      |  4 ++++
 physics/spectral/greens_function.py | 15 ++++++---------
 2 files changed, 10 insertions(+), 9 deletions(-)

diff --git a/.jules/bolt.md b/.jules/bolt.md
index 96776c7..c73b3cc 100644
--- a/.jules/bolt.md
+++ b/.jules/bolt.md
@@ -5,3 +5,7 @@
 ## 2025-11-01 - [Avoid Temporary Array Allocations in Reductions]
 **Learning:** In loops over parameters like temperature (`thermal_scan`), doing `np.sum(arr1 * arr2)` where array sizes match the system Hilbert space creates huge temporary arrays per loop iteration. Memory allocation dominates execution time.
 **Action:** Use `np.dot(arr1, arr2)` instead of `np.sum(arr1 * arr2)` for 1D arrays to evaluate reductions in C-level BLAS/NumPy functions, bypassing Python/NumPy array allocations and boosting performance significantly with less peak memory.
+
+## 2025-11-01 - [Vectorize Lattice Fourier Transforms]
+**Learning:** O(N^2) Python loops computing phase factors for lattice transformations \sum_{i,j} G_{ij} \exp(i k \cdot (r_i - r_j)) are extremely slow (taking ~100+ ms for N=200).
+**Action:** Vectorize this by replacing the loop with matrix algebra. Since \exp(i k \cdot (r_i - r_j)) = \exp(i k \cdot r_i) \exp(-i k \cdot r_j), we can construct a 1D phase array `phases = np.exp(1j * np.dot(r_vectors, k))` and then compute the final result as the quadratic form `phases @ G @ phases.conj()`. This reduces execution time to sub-milliseconds, bypassing the Python overhead entirely.
diff --git a/physics/spectral/greens_function.py b/physics/spectral/greens_function.py
index 2309e1c..4c264e9 100644
--- a/physics/spectral/greens_function.py
+++ b/physics/spectral/greens_function.py
@@ -127,15 +127,12 @@ def fourier_transform_lattice(greens_function: Array, lattice_k_vectors: Array,
     greens_function = np.asarray(greens_function, dtype=complex)
     k               = np.asarray(lattice_k_vectors)
     r_vectors       = np.asarray(lattice_r_vectors)
-    N               = len(r_vectors)
-    result          = 0.0 + 0.0j
-    for i in range(N):
-        r_i = r_vectors[i]
-        for j in range(N):
-            r_j = r_vectors[j]
-            phase = np.exp(1j * np.dot(k, r_i - r_j))
-            result += greens_function[i, j] * phase
-    return result
+
+    # Vectorized: \sum_{i,j} G_{ij} exp(i k\cdot (r_i - r_j))
+    # Let P_i = exp(i k \cdot r_i). Then phase_{i,j} = P_i * P_j^*
+    # G(k) = P^T G P^*
+    phases = np.exp(1j * np.dot(r_vectors, k))
+    return complex(phases @ greens_function @ phases.conj())
 
 def fourier_transform_lattice_translational(greens_function: Array, lattice_k_vectors: Array, lattice_r_vectors: Array) -> Array:
     r"""