Making Corr-Solver Fast

---

# 🚀 Making Corr-Solver Fast

## A JIT oracle, a silent termination bug, and the same bug in two languages

### 🐍 → 🐍 → ⚙️

@luk036 👨‍💻 · 2026 📅

---

### 📋 Agenda

**Part 1 — The Setup** 🧭
- What corr-solver solves
- The cutting-plane architecture
- How we measure (and why wall time lies) 📏

**Part 2 — Python: Where the Time Goes** 🐍
- The profile 🔍
- 72% is one factorization
- Fix 1: JIT the oracle ⚡

**Part 3 — The Bug Behind the Iterations** 🐛
- 2000 iterations for 4 variables 😳
- An unreachable tolerance
- 2000 → 61, bit-identical ✅

]

**Part 4 — Three Traps** 🪤
- The vectorization that passed 300 spot-checks 🙃
- The ellipsoid that returns *wrong* answers 💥
- Parallel cuts, and a claim I had to retract 📉

**Part 5 — C++: Same Bug, Different Physics** ⚙️
- The same bug in `ellalgo-cpp` 🐛
- Allocation churn is the new interpreter 🗑️
- $O(n^4) \to O(n^3)$ on the MLE gradient 📉

**Part 6 — Results & Lessons** 🎓

]

---

## 🧭 Part 1 — The Setup

---

### 🎯 What Corr-Solver Solves

Fits a **correlation polynomial** $\Omega(p)$ to a **biased sample covariance** $Y$ — recovering the underlying spatial correlation from noisy data. 📡

$$ \min_{p} \; \log\det {\color{#5e81ac}\Omega(p)} \;+\; \operatorname{Tr}\!\bigl({\color{#5e81ac}\Omega(p)}^{-1} {\color{#bf616a}Y}\bigr)
\qquad \text{s.t.} \quad 2{\color{#bf616a}Y} \succeq {\color{#5e81ac}\Omega(p)} \succeq 0 $$

with $\Omega(p) = \sum_{k} p_k \Sigma_k$ built from basis matrices $\Sigma_k = D_1^{\circ k}$ (element-wise powers of the distance matrix). 🧮

- $n = 20$ sites, $m = 4$ coefficients, $K = 4$ basis matrices
- Constraint is a **linear matrix inequality** — solved by an **ellipsoid method** 🥚
- The objective lives **on the boundary** of the PSD cone 😬

> Small problem. The interesting part is not the size — it is *where the time goes*. 🔍

---

### 🔄 The Cutting-Plane Architecture

<pre>
graph LR
    CP["cutting_plane_optim\n(ellalgo)"] --> SP["Ell\nsearch space"]
    CP --> OR["Oracle\nassess_optim"]
    OR --> L0["LMI0Oracle\nSigma x_k >= 0"]
    OR --> QM["QMIOracle\nt*I - F F^T >= 0"]
    L0 --> LD["LDLTMgr\nfactor / witness"]
    QM --> LD
    LD --> CUT["cut (g, beta)"]
    CUT --> CP
    style CP fill:#e3f2fd,stroke:#1565c0,stroke-width:3px
    style OR fill:#fff9c4,stroke:#f57f17,stroke-width:3px
    style LD fill:#bf616a,stroke:#8b2c3a,stroke-width:4px,color:#ffffff
    style SP fill:#d1c4e9,stroke:#6a1b9a,stroke-width:3px
</pre>

]

The oracle answers *"is $x$ feasible?"* and, if not, returns a **separating hyperplane**. 🔪

> Every iteration calls the oracle at least once. Every oracle call runs at least one **LDL^T factorization**. That is the whole story. 🎯

---

### 📏 How We Measure — and Why Wall Time Lies

First attempt: run the benchmark, read the number. Then run it again. 😵

<text x="200" y="134" text-anchor="end">run B (same code!)</text>
    <rect x="210" y="118" width="453" height="22" fill="#bf616a" rx="3"/>
    <text x="594" y="134" fill="#bf616a" font-weight="bold" text-anchor="end">54.9 s · 2.4× slower</text>
  </g>
  <text x="16" y="186" font-family="sans-serif" font-size="14" fill="#4c566a">Nothing changed. Only the machine state.</text>
  <text x="16" y="210" font-family="sans-serif" font-size="14" fill="#4c566a">Absolute timings are <tspan font-weight="bold" fill="#bf616a">not</tspan> comparable across runs.</text>
</svg>

</p>

**The protocol that does work:** 🧪
- Alternate baseline and candidate **in one process**, take the **min of 3** 🔁
- Report **CPU time** as a cross-check — if it disagrees with wall time, distrust the result 🖥️
- Trust **iteration counts** most of all: they are deterministic and machine-independent ✅

---

## 🐍 Part 2 — Python: Where the Time Goes

---

### 🔍 The Profile

`python -m cProfile -m pytest tests/test_corr_fn.py` — **45.3 s**, **12.6 M** function calls:

<text x="230" y="110" text-anchor="end">ellipsoid update + drivers</text>
    <rect x="240" y="94" width="131" height="22" fill="#5e81ac" rx="3"/>
    <text x="379" y="110" fill="#5e81ac">11.9 s · 26%</text>

<text x="230" y="156" text-anchor="end">copy.deepcopy (Ell clone)</text>
    <rect x="240" y="140" width="5" height="22" fill="#ebcb8b" rx="2"/>
    <text x="253" y="156" fill="#b58900">0.43 s · 1%</text>

<text x="230" y="202" text-anchor="end">data setup</text>
    <rect x="240" y="186" width="2" height="22" fill="#a3be8c" rx="1"/>
    <text x="250" y="202" fill="#2e7d32">0.05 s · 0.1%</text>
  </g>
  <text x="16" y="240" font-family="monospace" font-size="13" fill="#bf616a">of the 32.7 s</text>
  <text x="140" y="240" font-family="sans-serif" font-size="13" fill="#4c566a">— 18.3 s is Python-level element access (56%), 8.5 s loop bookkeeping.</text>
  <text x="16" y="262" font-family="sans-serif" font-size="13" fill="#4c566a">18,907 factorizations · 2.27 M element accesses · <tspan font-weight="bold" fill="#bf616a">~120 interpreter calls each</tspan></text>
</svg>

</p>

> The bottleneck is **not** a numpy operation. It is the **interpreter**, one matrix element at a time. 🐢

---

### 🐌 The Root Cause: A Callback Per Matrix Element

The LDL^T factorization is **lazy** — it asks for one element at a time through a Python closure:

```python
def get_elem(i, j):
    return sum(self.mat_f[k][i, j] * x[k] for k in range(n))
```

$$ \text{cost per factorization} \;=\; \underbrace{\frac{n(n+1)}{2}}_{\text{elements}} \times \underbrace{K}_{\text{Python } \texttt{sum}} \;=\; 210 \times 4 $$

- $n = 20$, $K = 4$ → each factorization makes **~120 interpreter calls** 🐍
- The matrices are **tiny**. There are no FLOPs to save — only **overhead to remove** 📦
- Measured: `get_elem` + `eval` is ${\color{#bf616a}56\%}$ of the factorization time

> The work is $O(n^3)$ flops. The *cost* is $O(n^2 K)$ **bytecode dispatches**. 🎯

---

### ⚡ Fix 1: JIT the Oracle

Reimplement the oracles as Numba `@njit` kernels — the LDL^T, the basis assembly, the gradient. Same arithmetic order, compiled. 🔥

<text x="230" y="176" text-anchor="end">mle_corr_bspline 🐍</text>
    <rect x="240" y="160" width="92" height="22" fill="#5e81ac" rx="3"/>
    <text x="340" y="176" fill="#5e81ac" font-weight="bold">13× · 1.22 s → 0.09 s</text>
  </g>
  <text x="16" y="222" font-family="sans-serif" font-size="13" fill="#4c566a">Bar width ∝ speedup. The pure-Python path stays available — the kernels are an</text>
  <text x="16" y="244" font-family="sans-serif" font-size="13" fill="#4c566a">opt-in module behind a <tspan font-family="monospace" fill="#bf616a">numba</tspan> extra, not a replacement.</text>
</svg>

</p>

---

### 🧪 Verifying a Rewrite That Must Not Change the Answer

A faster oracle is worthless if it takes a different path. 😬

- **Kernel-level:** the compiled LDL^T matches `LDLTMgr` on `pos`/`D`/`L`/witness/`ep` over **6,000 random symmetric matrices** ✅
- **Oracle-level:** outputs agree to `1e-8` on random points ✅
- **Trajectory-level:** **0 feasibility-verdict mismatches** across **1,437** real solver-trajectory points ✅

Why the care? The QMI is evaluated **at the PSD boundary**, where the smallest pivot is $\approx 0$:

$$ d_i \;=\; A_{ii} - \sum_{k<i} L_{ik}^2 D_k \;\approx\; {\color{#bf616a}0} \quad\Longrightarrow\quad \operatorname{sign}(d_i) \text{ is decided by rounding} $$

Any reordering of the arithmetic can flip **which pivot** goes negative first → a different cut → a different trajectory. 🎲

> Correct answer, different iteration count, is *expected*. Different **answer** is not. ✅

---

## 🐛 Part 3 — The Bug Behind the Iterations

---

### 😳 2000 Iterations for a 4-Variable Problem

A 4-variable problem should not need 2000 cutting-plane steps.

| case | iterations | stop reason |
|:-----|-----------:|:------------|
| `lsq_corr_poly` | **2000** | `max_iters` 😬 |
| `lsq_corr_poly2` | 1089 | tolerance |
| `mle_corr_bspline` | 371 | `NoEffect` ✅ |

And the test assertion was **vacuous**:

```python
assert num_iters <= 2000   # the old code ALWAYS returned exactly 2000
```

> An assertion that cannot fail is not a test. 🧪

---

### 📉 The Trace: `tau` Pinned at 1.11e-16

Instrumented bisection: `tau` reaches float resolution after ~62 steps — then **stays there for 1938 more**.

<line x1="70" y1="217" x2="610" y2="217" stroke="#bf616a" stroke-width="2" stroke-dasharray="6,4"/>
  <text x="612" y="221" font-family="sans-serif" font-size="11" fill="#bf616a" text-anchor="end" dx="-4">tolerance = 1e-20 — never reached ⛔</text>

<circle cx="370" cy="187" r="5" fill="#a3be8c" stroke="#2e7d32" stroke-width="2"/>
  <text x="366" y="176" font-family="sans-serif" font-size="12" fill="#2e7d32" text-anchor="end" font-weight="bold">62 steps: converged</text>

<line x1="372" y1="45" x2="372" y2="235" stroke="#2e3440" stroke-width="1" stroke-dasharray="3,3"/>
  <line x1="378" y1="45" x2="378" y2="235" stroke="#2e3440" stroke-width="1" stroke-dasharray="3,3"/>
  <text x="375" y="42" font-family="sans-serif" font-size="11" fill="#2e3440" text-anchor="middle">axis break</text>

<text x="490" y="178" font-family="sans-serif" font-size="12" fill="#bf616a" font-weight="bold">1938 steps doing nothing 🗑️</text>
  <text x="16" y="262" font-family="sans-serif" font-size="13" fill="#4c566a">tau bottoms out at the <tspan font-weight="bold">ulp of the threshold value</tspan> (1.11e-16) and cannot shrink further.</text>
  <text x="16" y="284" font-family="sans-serif" font-size="13" fill="#4c566a">The loop is not converging — it is <tspan font-weight="bold" fill="#bf616a">spinning</tspan>.</text>
</svg>

</p>

---

### 🧮 Why: An Absolute Tolerance on a Scale-Dependent Width

`bsearch` compared a **width** against an **absolute** constant:

```python
tau = (upper - lower) / 2
if tau < options.tolerance:      # tolerance = 1e-20
    return upper, niter
```

$$ \text{stop} \iff \tau < 10^{-20}, \qquad \tau \;\ge\; \operatorname{ulp}\!\bigl(\max(|{\color{#5e81ac}\text{lower}}|, |{\color{#bf616a}\text{upper}}|)\bigr) \approx {\color{#bf616a}10^{-16}} $$

The bracket converges to $t^\* = 1.4475$, not to $0$ — so $\tau$ bottoms out at $\operatorname{ulp}(1.4475) \approx 1.11 \times 10^{-16}$. 📌

$$ {\color{#bf616a}1.11 \times 10^{-16}} \;\not<\; {\color{#5e81ac}10^{-20}} \quad\Longrightarrow\quad \text{unreachable, for any realistic scale} $$

> An absolute tolerance on a scale-dependent quantity is a bug waiting for a scale. 💣

---

### ✅ The Fix: A Stall Guard

Three lines, no new options, **bit-identical result**:

```python
gamma = lower + tau
if not (lower < gamma < upper):
    return upper, niter        # the midpoint can no longer move
```

| | before | after |
|:--|-------:|------:|
| iterations | 2000 | **61** |
| objective $t$ | 1.4475055597437385 | **1.4475055597437385** ✅ |
| stop reason | `max_iters` | stall |

The guard **cannot** change the answer: it fires only when no further refinement is *representable*. 🎯

> 61 is the information-theoretic floor: $\log_2\!\bigl(939 / \operatorname{ulp}(1.4475)\bigr) \approx 62$. 📐

---

### 💥 The Inner Work Collapsed Too

The wasted outer steps were each driving **real subproblem solves**. Removing them removes those too:

<text x="150" y="194" text-anchor="end">site=10x8 m=4</text>
    <rect x="160" y="180" width="360" height="18" fill="#bf616a" rx="3"/>
    <text x="526" y="194" fill="#bf616a" font-weight="bold" text-anchor="end">901,861</text>
    <rect x="160" y="200" width="222" height="18" fill="#a3be8c" rx="3"/>
    <text x="388" y="214" fill="#2e7d32" font-weight="bold">16,195 · 56×</text>
  </g>
  <text x="16" y="244" font-family="sans-serif" font-size="13" fill="#4c566a">Outer iterations went 2000 → ~62 on every instance in a 10-case battery.</text>
</svg>

</p>

---

### 📦 Root Cause Lives in `ellalgo` → v0.8

The bug is in the **shared driver**, not in corr-solver — so it was fixed at the source. 🔧

- `ellalgo` is a separate repo; corr-solver installs it from git 🔗
- Fix landed with a **regression test that fails without the guard** (2000 vs <200) 🧪
- Released as **`ellalgo` v0.8**, then corr-solver's bound tightened `2000 → 100` 🔒

```python
# tests/test_corr_fn.py
assert num_iters <= 100    # fails against an unpatched ellalgo
```

> The vacuous assertion is now a real one — and it catches a **cross-repo** regression. 🌉

---

## 🪤 Part 4 — Three Traps

---

### 🙃 Trap 1: The Vectorization That Passed 300 Spot-Checks

I replaced the row-wise LDL^T with a vectorized one. It agreed on **300/300** random points. ✅

It was **wrong**: it read `L[i, :i]` *before* those entries were computed. 🐛

| | correct | "optimized" |
|:--|--------:|-----------:|
| `lsq_corr_poly` iterations | 2000 | **63** ❌ |
| `ep` (violation) | 1428922.09 | 198.25 ❌ |
| verdicts on random points | — | 300/300 agree 😬 |

The fix: the vectorized LDL^T must be **column-wise (left-looking)**, not row-wise. 🧭

> **Spot-checking random points is not validation.** Assert the *iteration count*. 🎯

---

### 💥 Trap 2: The Ellipsoid That Returns *Wrong* Answers

"Tighten the initial ellipsoid" looked like free speed. It is not. 🪤

| initial box | iterations | objective $\gamma$ |
|:--|-----------:|------------------:|
| `val[-1] = normY2**2` (current) | 1059 | **1.4475** ✅ |
| `val[-1] = normY2**2 * 1e-6` | **181** ⬅ faster! | **14961** ❌ |

The true optimum falls **outside** the shrunken box, so the solver returns the best point *inside* it — a valid-looking, **silently wrong** answer. 😱

For the MLE path, $\alpha = 0.1$ fails outright: $\gamma = \texttt{inf}$. 💀

> A speedup that changes the answer is not a speedup. It is a bug with a nice benchmark. 🚨

---

### 📉 Trap 3: Parallel Cuts — and a Claim I Had to Retract

**Parallel cuts:** `EllCalc.use_parallel_cut` is already `True`, but it only activates when the oracle returns a **pair** of betas. These oracles return **scalars**, so it never fires. Enabling it means changing oracle semantics for ~2×. **Rejected.** ❌

**A claim I published and then had to retract:** 📝

$$ \text{I claimed } \Sigma_k = D_1^k \text{ are matrix powers} \implies \text{they commute} \implies \text{shared eigenbasis} \implies O(n^2) \text{ per iteration} $$

But the code says:

```python
D = np.multiply(D, D1)     # element-wise (Hadamard) powers!
```

Hadamard powers of a symmetric matrix are symmetric — and **generally do not commute**. The shortcut does not apply. 🚫

> I filed the wrong claim in a public issue, then corrected it. Retracting beats rationalising. ✅

---

## ⚙️ Part 5 — C++: Same Bug, Different Physics

---

### 🐛 The Same Bug in `ellalgo-cpp` → v1.6.14

The C++ port had the **identical** defect — an absolute `tolerance` against a scale-dependent `tau`:

```cpp
auto tau = algo::half_nonnegative(upper - lower);
if (tau < options.tolerance) return {upper, niter};
// ... spins to max_iters for the same reason
```

Same guard, same proof:

| | before | after |
|:--|-------:|------:|
| iterations | 2000 | **< 200** |
| threshold found | 500.0 | 500.0 ✅ |

- Regression test **fails without the guard**: `CHECK( 2000 < 200 )` 🧪
- Released as **`ellalgo-cpp` v1.6.14**; 138 cases / 502 assertions ✅

> Interesting: corr-solver-cpp's LSQ driver uses `cutting_plane_optim`, **not** `bsearch` — so it never hit this. The bug was live for every *other* consumer. 🌍

---

### 🗑️ In C++, the Interpreter Is the Allocator

There is no bytecode to remove. The analogous overhead is **allocation churn in the hot path**:

- `MleOracle::assess_optim` allocated **five $n \times n$ buffers** per call 🗑️
- `QmiOracle::assess_feas` allocated four more, including an $O(n^2)$ scratch 📦
- `trace()` built a **vector** just to sum its own diagonal 🤦
- `norm()` built a full `a * a` **copy** 📋

And the matrix multiply walked `B` **column by column** (stride $n$) — cache-hostile and unvectorizable:

```cpp
// i-k-j order: B(t,:) and out(i,:) are contiguous → the inner loop vectorizes
for (i) for (t) { const auto a = A(i,t); for (j) out(i,j) += a * B(t,j); }
```

> Same lesson, different physics: find the per-iteration overhead, then remove it. 🎯

---

### 📉 MLE Gradient: $O(n^4) \to O(n^3)$

The gradient computed $n$ separate $n \times n$ products **per oracle call** — $O(n^4)$ overall. 😱

But the trace has a cyclic-permutation identity:

$$ g_i \;=\; \operatorname{Tr}\!\bigl(S{\color{#5e81ac}\Sigma_i}\bigr) - \operatorname{Tr}\!\bigl({\color{#5e81ac}\Sigma_i} S Y S\bigr)
\;=\; \bigl\langle\, {\color{#bf616a}S - SY S},\; {\color{#5e81ac}\Sigma_i} \,\bigr\rangle_F $$

Forming ${\color{#bf616a}V = S - SY S}$ **once** turns a per-$i$ $O(n^3)$ product into an $O(n^2)$ inner product. 🪄

<text x="170" y="134" text-anchor="end">after ⚙️</text>
    <rect x="180" y="118" width="273" height="22" fill="#a3be8c" rx="3"/>
    <text x="461" y="134" fill="#2e7d32" font-weight="bold">110.8 ms · −24% ✅</text>
  </g>
  <text x="16" y="180" font-family="sans-serif" font-size="13" fill="#4c566a">Iteration count unchanged at 200 — the algebra is identical, only the evaluation order changed.</text>
</svg>

</p>

---

### 🧹 De-duplicating the Benchmark (−230 Lines)

The benchmark contained a **verbatim copy** of the implementation. 📄📄

- Any fix in `source/lsq_corr_ell.cpp` would **not** reach the benchmark 🚫
- I was about to fix the same bug **twice** — the signal that something is wrong 🚩
- Now the benchmark declares `extern` and links the real library ✅

```lua
-- xmake.lua
add_files("benchmark/source/*.cpp")
add_files("source/*.cpp")          -- the real implementation, once
```

> If a benchmark copies the code it measures, it measures a **different program**. 🎯

---

## 🎓 Part 6 — Results & Lessons

---

### 🏁 Cross-Language Scoreboard

Python, same instance, controlled alternating A/B (min of 3, CPU time agrees):

<text x="240" y="180" text-anchor="end">mle_corr_bspline 🐍</text>
    <rect x="250" y="164" width="77" height="22" fill="#5e81ac" rx="3"/>
    <text x="335" y="180" fill="#5e81ac" font-weight="bold">13× · 1.2 s → 0.09 s</text>

<text x="240" y="236" text-anchor="end">MLE_corr ⚙️</text>
    <rect x="250" y="220" width="139" height="22" fill="#d1c4e9" rx="3"/>
    <text x="397" y="236" fill="#6a1b9a" font-weight="bold">1.24× · 146.0 → 110.8 ms</text>
  </g>
  <text x="16" y="274" font-family="sans-serif" font-size="13" fill="#4c566a">Bar width ∝ speedup. Iteration counts are the machine-independent evidence. ✅</text>
</svg>

</p>

---

### 📜 Artifacts From This Session

| Repo | Change |
|:-----|:-------|
| `py/corr-solver` | `numba_oracles.py`: JIT oracles (LDL^T, QMI, LMI0, LMI, lsq, mle) ⚡ |
| `py/corr-solver` | 10 new tests: kernel equivalence, oracle equivalence, trajectory replay 🧪 |
| `py/corr-solver` | bound tightened `2000 → 100`; **15 tests, 100% coverage** ✅ |
| `py/ellalgo` | `bsearch`: floating-point stall guard + regression test 🐛 |
| `py/ellalgo` | **released v0.8** — 164 tests, `main` + tag + GitHub Release 📦 |
| `cpp/ellalgo-cpp` | same guard + regression test; **released v1.6.14** ⚙️ |
| `cpp/corr-solver-cpp` | MLE gradient $O(n^4) \to O(n^3)$; `matmul` i-k-j; no-alloc `trace`/`norm` 🚀 |
| `cpp/corr-solver-cpp` | benchmark de-duplicated (−230 lines), now links the real library 🧹 |
| `docs` | issue #1 filed — and corrected after the Hadamard-powers retraction 📝 |

]

> Two releases, one bug fixed in two languages, and a wrong claim retracted. 🎉

---

### 🧠 Key Takeaways

**Performance** ⚡
- **Profile first** — 72% was one factorization, 56% of *that* was element access 🔍
- Match the fix to the physics: Python → **JIT**; C++ → **stop allocating** 🐍⚙️
- An absolute tolerance on a scale-dependent width is a **latent bug** 💣
- `2000 → 61` iterations came from **termination**, not from a better algorithm 📉

**Correctness** 🐛
- A speedup that changes the answer is a **bug with a nice benchmark** 🚨
- **Spot-checking is not validation** — assert the iteration count 🧪
- Bit-identical is a *feature*: it proves you changed the cost, not the math ✅
- Validate the invariant that matters — here, **feasibility verdicts** (0/1437 mismatches) 🎯

**Process** 🔁
- Machine state drifts **2.4×**; alternate A/B and use CPU time 📏
- If a fix must be applied **twice**, that is the bug 🚩
- Retracting a published claim beats rationalising it 📝

---

# 🙋 Q&A

**Making Corr-Solver Fast — and Correct**

### Questions? Discussion? 💬

---

# 👏 Thank You

## Profile. Measure. Verify the *answer*, not just the clock. 🎯

Slides built with Remark.js 🎉 | KaTeX 📐 | Mermaid 🧩 | Nord Theme 🌙