Making the Multiplierless FIR Designer Fast

---

# ⚡ Making the Multiplierless FIR Designer Fast

## A Python 94× win, a C++ 2× regression, and a Rust 1.3× kernel

### 🐍 → ⚙️ → 🦀

@luk036 👨‍💻 · 2026 📅

---

### 📋 Agenda

**Part 1 — The Setup** 🧭
- What "multiplierless" actually does
- The ellipsoid optimization loop
- How we measure 📏

**Part 2 — Python: The 94× Find** 🐍
- The profile surprise 🔍
- Vectorizing the constraint scan
- Results & equivalence ✅

]
.pull-right[

**Part 3 — C++: When the Obvious Port Fails** ⚙️
- The naive port regressed 2× 📉
- Early exit beats a full matvec
- Raw pointers, not SIMD 🎯

**Part 4 — Rust: The Same Lesson** 🦀
- The scalar `.sum()` trap 🕳️
- A 4-accumulator kernel ⚡
- 1.3× on the oracle, flat end-to-end

]

---

## 🧭 Part 1 — The Setup

---

### 🎛️ What Is "Multiplierless"?

A **multiplierless FIR filter** replaces every multiply by **shift + add**, using the
**Canonical Signed Digit (CSD)** representation. 🪄

<pre>
graph LR
    SPEC["Filter spec\nN, wpass, wstop"]
    OPT["Ellipsoid optimizer\ncutting-plane"]
    SF["Spectral factor\nR(w) → h"]
    CSD["CSD quantize\nh → ±2^k"]
    VER["Verilog\nshift + add"]
    SPEC --> OPT
    OPT -->|"feasible r"| SF
    SF --> CSD
    CSD -->|"residual cut"| OPT
    CSD --> VER
    style SPEC fill:#e3f2fd,stroke:#1565c0,stroke-width:3px
    style OPT fill:#fff9c4,stroke:#f57f17,stroke-width:3px
    style SF fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px
    style CSD fill:#f8bbd0,stroke:#c2185b,stroke-width:3px
    style VER fill:#d1c4e9,stroke:#6a1b9a,stroke-width:3px
</pre>

]

> The output is an **intermediate step** in an iterative design loop. 🔁

---

### 🧮 The Optimization Problem

Minimize the worst stopband response, subject to passband bounds:

$$ \min_{\gamma} \; \gamma \quad \text{s.t.} \quad
{\color{green} L^2 \le R(\omega) \le U^2},\; \omega \in \Omega_p,
\quad {\color{red} R(\omega) \le \gamma},\; \omega \in \Omega_s $$

where the squared-magnitude response is a cosine series in the autocorrelation $r$:

$$ R(\omega) \;=\; r_0 \;+\; 2\sum_{k=1}^{N-1} r_k \cos(k\omega)
   \;=\; {\color{#1565c0} A(\omega)\, r} $$

> Every ellipsoid iteration asks the oracle: **"is this $r$ feasible? if not, give me a cut."** ✂️

---

### 🔬 How We Measure

- `python -m cProfile -m multiplierless.fir_design sample_filter.json`
- `pytest --benchmark-only` (pytest-benchmark) 📊
- Canonical case: $N = 32$, `csd_nnz = 4`

| Signal | Value |
|:-------|:------|
| Iterations | 1 482 |
| Wall clock | ≈ 3.0 s |
| Import cost | ≈ 1.4 s (numpy + ellalgo) |

> Profiling before optimizing is not optional. 🔎

---

## 🐍 Part 2 — Python: The 94× Find

---

### 🔍 The Profile Surprise

The FFT — the "obviously expensive" part — was **only 6%**. The scan was **83%**.

<text x="180" y="104" text-anchor="end">spectral_fact_fft</text>
    <rect x="190" y="90" width="29" height="20" fill="#5e81ac" rx="3"/>
    <text x="227" y="105" fill="#5e81ac">0.22 s · 6%</text>

<text x="180" y="144" text-anchor="end">csd_quantize</text>
    <rect x="190" y="130" width="20" height="20" fill="#ebcb8b" rx="3"/>
    <text x="218" y="145" fill="#b58900">0.15 s · 4%</text>

<text x="180" y="184" text-anchor="end">ellipsoid update</text>
    <rect x="190" y="170" width="3" height="20" fill="#a3be8c" rx="1"/>
    <text x="201" y="185" fill="#2e7d32">0.01 s</text>
  </g>
  <text x="16" y="238" font-family="monospace" font-size="13" fill="#bf616a">553 912</text>
  <text x="86" y="238" font-family="sans-serif" font-size="13" fill="#4c566a">tiny per-row ndarray.dot calls — ~373 per iteration</text>
</svg>

</p>

---

### 🐌 The Offending Loop

```python
def _scan_constraints(rr, count, x, spectrum, check):
    for _ in range(count):              # up to 480 rows
        k = rr.next()
        col_k = spectrum[k]
        v = col_k.dot(x)                # 553 912 calls!
        if cut := check(k, col_k, v):
            return cut
    return None
```

- `assess_feas` calls it **3×** per iteration (passband, stopband, transition)
- Each dot is length $N = 32$ — **pure Python dispatch overhead** dominates
- 4 217 scan calls, **553 912** dot calls total 🐢

---

### 💡 The Insight: Compute Once, Mask Once

Instead of ~400 tiny dots, do **one** matrix–vector product and test bounds with masks:

$$ v = {\color{#1565c0} A}\, x, \qquad
   v_i = \sum_{j} A_{ij} x_j $$

<pre>
graph TD
    subgraph BEFORE["Before · per-row"]
        B1["row k"] --> B2["dot(row, x)"] --> B3["check"] --> B4["row k+1"]
    end
    subgraph AFTER["After · vectorized"]
        A1["A @ x\none BLAS call"] --> A2["boolean masks\n+ argmax"]
    end
    style B1 fill:#f8bbd0,stroke:#c2185b,stroke-width:3px
    style B2 fill:#f8bbd0,stroke:#c2185b,stroke-width:3px
    style B3 fill:#f8bbd0,stroke:#c2185b,stroke-width:3px
    style B4 fill:#f8bbd0,stroke:#c2185b,stroke-width:3px
    style A1 fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px
    style A2 fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px
</pre>

]

---

### ⚡ The Micro-Benchmark

<text x="150" y="126" text-anchor="end">S @ x matvec</text>
    <rect x="160" y="110" width="5" height="24" fill="#a3be8c" rx="2"/>
    <text x="175" y="127" fill="#2e7d32" font-weight="bold">8.51 µs</text>
  </g>
  <text x="160" y="178" font-family="sans-serif" font-size="18" font-weight="bold" fill="#5e81ac">≈ 94× faster on the primitive</text>
</svg>

</p>

---

### 🧩 The Vectorized `assess_feas`

```python
def assess_feas(self, x):
    v = self.spectrum @ x                       # one BLAS matvec

seg = v[:nwpass]                            # passband
    mask = (seg > self.up_sq) | (seg < self.lp_sq)
    k = _first_true_cyclic(mask, _peek(self.idx1, 0, nwpass), 0, nwpass)
    if k >= 0:
        self.idx1._cur = k                      # preserve round-robin cursor
        ...
```

- **Same cut selection**, same cursor state
- `_first_true_cyclic` picks the first violation in **cyclic order** from the cursor 🔄

---

### 🔄 Preserving Round-Robin Semantics

The round-robin cursor is what stops the ellipsoid from re-adding the same cut. It must survive.

<pre>
graph LR
    PEEK["_peek(rr)\nread cursor\n(no advance)"]
    MASK["mask.any()\nfirst True\ncyclic from start"]
    SET["rr._cur = k\nland on violation"]
    NONE["no violation\ncursor unchanged"]
    PEEK --> MASK
    MASK -->|"found"| SET
    MASK -->|"none"| NONE
    style PEEK fill:#e3f2fd,stroke:#1565c0,stroke-width:3px
    style SET fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px
    style NONE fill:#fff9c4,stroke:#f57f17,stroke-width:3px
</pre>

]

> `RoundRobin` exposes no getter — the cursor slot is read directly. 🔧

---

### 📊 Python Results

<text x="120" y="202" text-anchor="end">end-to-end wall</text>
    <rect x="130" y="184" width="329" height="16" fill="#bf616a" rx="2"/>
    <text x="466" y="197" fill="#bf616a">3.01 s</text>
    <rect x="130" y="202" width="204" height="16" fill="#a3be8c" rx="2"/>
    <text x="341" y="215" fill="#2e7d32" font-weight="bold">1.87 s · 1.6×</text>
  </g>
  <text x="16" y="262" font-family="sans-serif" font-size="13" fill="#4c566a">After the scan fix the run is <tspan font-weight="bold" fill="#5e81ac">import-bound</tspan> (~1.4 s of numpy/ellalgo import).</text>
</svg>

</p>

---

### ✅ Behavior Equivalence

The vectorized version is **logically identical** — only ~$10^{-15}$ float noise differs.

| Property | Old | New |
|:---------|:----|:----|
| Gradient returned | same | **same** |
| Round-robin cursor | same | **same** |
| `kmax` tie-break | first in scan order | **same** |
| Iterations | 1 483 | 978 |
| Passband $|H|$ | [0.8827, 1.1249] | [0.8838, 1.1234] |
| Stopband max $|H|$ | 0.0071 | 0.0072 |

]

> `matvec` vs `dot` differs in the last bits; the chaotic iteration amplifies it into a **different, equally valid** trajectory. 🎲

---

### 📉 What Moved After the Fix

<pre>
graph LR
    SCAN["constraint scan\n0.082 s ✅"]
    FFT["spectral_fact\n0.15 s ⏭️ next"]
    CSD["csd_quantize\n0.14 s ⏭️ next"]
    IMP["imports\n1.4 s 🐘"]
    IMP --> SCAN
    SCAN --> FFT
    FFT --> CSD
    style SCAN fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px
    style FFT fill:#fff9c4,stroke:#f57f17,stroke-width:3px
    style CSD fill:#fff9c4,stroke:#f57f17,stroke-width:3px
    style IMP fill:#f8bbd0,stroke:#c2185b,stroke-width:3px
</pre>

]

> The bottleneck always moves. Next targets: FFT buffers and CSD quantization. 🎯

---

## ⚙️ Part 3 — C++: When the Obvious Port Fails

---

### 👥 The Sibling Project

`multiplierless-cpp` mirrors the Python package — same oracle, same math, **compiled**.

<pre>
graph LR
    PY["multiplierless.py\noracle + FFT + CSD"]
    CPP["multiplierless-cpp\nArr · RoundRobin · FFTW"]
    PY -.->|"same algorithm"| CPP
    CM["CMake + CPM\npinned deps"]
    XM["xmake\nlocal siblings"]
    CM --> CPP
    XM --> CPP
    style PY fill:#e3f2fd,stroke:#1565c0,stroke-width:3px
    style CPP fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px
    style CM fill:#d1c4e9,stroke:#6a1b9a,stroke-width:3px
    style XM fill:#ffe0b2,stroke:#e65100,stroke-width:3px
</pre>

]

> "Enhance the sibling accordingly" — sounds simple. It wasn't. 😅

---

### 📉 The Naive Port Regressed 2×

I ported the Python idea literally: one `dot(mat, x)` per segment, then scan.

<text x="150" y="136" text-anchor="end">full-matvec port</text>
    <rect x="160" y="118" width="400" height="24" fill="#bf616a" rx="3"/>
    <text x="568" y="135" fill="#bf616a" font-weight="bold">0.327 s</text>
  </g>
  <text x="160" y="186" font-family="sans-serif" font-size="16" font-weight="bold" fill="#bf616a">❌ 1.9× slower — a regression, not an optimization</text>
</svg>

</p>

---

### 🤔 Why? Early Exit Beats a Full Matvec

In C++ the per-row dot is **already compiled** — and the scan **stops at the first violation**.

$$ \underbrace{\textstyle\sum_{\text{scanned}} \|\text{row}\|}_{\text{early exit:} \approx \tfrac{1}{2}m}
   \;<\; \underbrace{\textstyle\sum_{i=1}^{m} \|\text{row}\|}_{\text{full matvec: } m} $$

- Python's win came from killing **interpreter overhead**, not from FLOPs 🐍
- C++ has no interpreter overhead to kill
- Computing every row when half suffice = **2× the work** 🔁

> The lesson: port the *measurement*, not the *conclusion*. 🔬

---

### 🎯 The Real Win: Raw Pointer Access

`mat(row, j)` hides `_data[row*cols + j]`. Exposing pointers lets the compiler vectorize.

```cpp
auto dot_row(const Arr& mat, std::size_t row, const Arr& x) -> double {
    const double* row_data = mat.data() + row * mat.cols();
    const double* x_data = x.data();
    const auto n = x.size();
    double sum = 0.0;
    for (std::size_t j = 0; j < n; ++j) sum += row_data[j] * x_data[j];
    return sum;
}
```

</p>

---

### 🧪 SIMD Was Not the Bottleneck

I tried the usual suspects. None helped.

| Build change | Result |
|:-------------|:-------|
| `/arch:AVX2` only | 0.170 s — no change |
| `/fp:fast` only | 0.157 s — marginal |
| `/fp:fast` + `/arch:AVX2` | 0.155 s — marginal |
| SSE2 vs AVX2 | **identical** → not FP-throughput bound |

]

> The limit was **loop/pointer overhead**, not arithmetic width. 🚫🏎️

---

### 🖥️ End-to-End: The CLI Didn't Change

The CLI uses the **quantized** path, dominated by FFT + CSD — not the scan.

<text x="150" y="126" text-anchor="end">C++ ⚙️</text>
    <rect x="160" y="108" width="42" height="24" fill="#a3be8c" rx="3"/>
    <text x="210" y="125" fill="#2e7d32" font-weight="bold">0.197 s · ~9× faster</text>
  </g>
  <text x="16" y="182" font-family="sans-serif" font-size="13" fill="#4c566a">The C++ sibling was already fast — so the Python optimization had little to transfer.</text>
</svg>

</p>

---

## 🦀 Part 4 — Rust: The Same Lesson

---

### 🦀 The Third Sibling

`multiplierless-rs` mirrors the same oracle — same cosine matrix, same
round-robin, same early-exit scan.

<pre>
graph LR
    PY["multiplierless.py\noracle + FFT + CSD"]
    CPP["multiplierless-cpp\nArr · RoundRobin · FFTW"]
    RS["multiplierless-rs\nArr · RoundRobin · rustfft"]
    PY -.->|"same algorithm"| CPP
    PY -.->|"same algorithm"| RS
    style PY fill:#e3f2fd,stroke:#1565c0,stroke-width:3px
    style CPP fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px
    style RS fill:#ffe0b2,stroke:#e65100,stroke-width:3px
</pre>

]

> Third language, same scan — so we already knew what to expect. 🤔

---

### 🕳️ The Scalar `.sum()` Trap

The row dot product was a **single dependent f64 chain**:

```rust
let v: f64 = row.iter().zip(x).map(|(a, b)| a * b).sum();
```

LLVM **cannot reassociate** floating-point without fast-math, so it cannot
vectorize the reduction. The chain latency grows with $N$:

$$ T_{\text{chain}} \;\approx\; N \times \text{latency}(\text{add}) $$

> A 32-wide dot becomes ~32 dependent adds — the arithmetic units sit idle. 😴

---

### ⚡ The 4-Accumulator Kernel

```rust
#[inline]
fn dot_row(row: &[f64], x: &[f64]) -> f64 {
    let (mut a0, mut a1, mut a2, mut a3) = (0.0, 0.0, 0.0, 0.0);
    let mut r = row.chunks_exact(4);
    let mut xc = x.chunks_exact(4);
    for (r4, x4) in r.by_ref().zip(xc.by_ref()) {
        a0 += r4[0] * x4[0];
        a1 += r4[1] * x4[1];
        a2 += r4[2] * x4[2];
        a3 += r4[3] * x4[3];
    }
    (a0 + a1) + (a2 + a3)
}
```

`chunks_exact(4)` removes bounds checks; four independent partial sums expose
the parallelism **without** asking the compiler to reassociate anything.

<pre>
graph LR
    subgraph ONE["1 accumulator · 32 deep"]
        I1["x0"] --> I2["x1"] --> I3["x2"] --> I4["... x31"]
    end
    subgraph FOUR["4 accumulators · 8 deep"]
        F0["a0: x0, x4, ..."]
        F1["a1: x1, x5, ..."]
        F2["a2: x2, x6, ..."]
        F3["a3: x3, x7, ..."]
    end
    style I1 fill:#f8bbd0,stroke:#c2185b,stroke-width:3px
    style I2 fill:#f8bbd0,stroke:#c2185b,stroke-width:3px
    style I3 fill:#f8bbd0,stroke:#c2185b,stroke-width:3px
    style I4 fill:#f8bbd0,stroke:#c2185b,stroke-width:3px
    style F0 fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px
    style F1 fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px
    style F2 fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px
    style F3 fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px
</pre>

]

---

### 📊 Rust Results

<text x="190" y="131" text-anchor="end">assess_optim (full scan)</text>
    <rect x="200" y="115" width="360" height="20" fill="#bf616a" rx="3"/>
    <text x="568" y="130" fill="#bf616a">13.9 µs</text>
    <rect x="200" y="137" width="278" height="20" fill="#a3be8c" rx="3"/>
    <text x="486" y="152" fill="#2e7d32" font-weight="bold">10.7 µs · 1.3×</text>

</p>

> The kernel wins big on the **oracle hot path**, but end-to-end is **unchanged** — 8089 iterations both ways. 🎯

---

### 🎓 The Same Lesson, Third Time

<pre>
graph LR
    Q["Where is the time?"]
    Q -->|"interpreter overhead"| PY["Python\nhuge win ✅"]
    Q -->|"compiled, early exit"| CC["C++ / Rust\nscan not the bottleneck"]
    CC --> K["kernel still faster\n(oracle API) ✅"]
    CC --> E["end-to-end flat ⚠️"]
    style Q fill:#e3f2fd,stroke:#1565c0,stroke-width:3px
    style PY fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px
    style CC fill:#fff9c4,stroke:#f57f17,stroke-width:3px
    style K fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px
    style E fill:#f8bbd0,stroke:#c2185b,stroke-width:3px
</pre>

]

- The real trajectory is dominated by **short early-exit scans**, not the dot 🐢
- So the compiled ports stay flat even when the kernel improves 📏
- **Port the measurement, not the conclusion** — again 🔬

---

### 🧠 Key Takeaways

**Performance** ⚡
- **Profile first** — the bottleneck is rarely where you guess
- Python's win was *interpreter overhead*, not FLOPs
- Port the **measurement**, not the conclusion
- Early exit can beat vectorization in compiled code
- Raw pointers + alias info enable the optimizer 🎯
- Rust's `.sum()` is a scalar chain — 4 accumulators fix it ⚡

---

### 📜 Artifacts From This Session

| Artifact | Link |
|:---------|:-----|
| Python issue | `luk036/multiplierless#3` |
| Python change | `lowpass_oracle.py` vectorized scan |
| C++ change | `lowpass_oracle.cpp` raw-pointer dot |
| Rust change | `lowpass_oracle.rs` 4-accumulator dot 🦀 |
| Benchmarks | `test_benchmarks.py` · `benches/lowpass.rs` |

> Python: 94× primitive, 1.6× end-to-end · C++: 7% · Rust: 1.3× oracle 🎉

---

# 🙋 Q&A

**Making the Multiplierless FIR Designer Fast**

### Questions? Discussion? 💬

---

# 👏 Thank You

## Profile. Measure. Port the *idea*, not the code. 🎯

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