Replacing xtensor with a Custom Arr Class

---

# 🗑️ Replacing xtensor with a Custom `Arr` Class

Eliminating 7 dependencies from a C++17 FIR filter library

@luk036 👨‍💻 · 2026 📅

---

## 🚨 The Problem

**xtensor** is a powerful N-dimensional array library. But for small-scale optimization problems, the cost outweighs the benefit.

###  Before: Huge Dependency Stack

.mermaid[
<pre>
graph TD
    subgraph "Application"
        A[multiplierless-cpp]
    end
    subgraph "xtensor Stack"
        B[xtensor-fftw]
        C[xtensor-blas]
        D[xtensor 0.25.0]
        E[xtl 0.7.7]
    end
    subgraph "System Libraries"
        F[FFTW]
        G[BLAS / LAPACK]
        H[OpenBLAS]
    end
    subgraph "Build Tools"
        I[CMake + CPM.cmake]
    end
    A --> B
    B --> D
    C --> D
    C --> G
    C --> H
    D --> E
    A --> C
    A --> F
    I --> A
    I --> B
    I --> C
    I --> D
    I --> E

style A fill:#f44336
    style B fill:#ff9800
    style C fill:#ff9800
    style D fill:#ff9800
    style E fill:#ff9800
    style F fill:#4caf50
    style G fill:#9c27b0
    style H fill:#9c27b0
    style I fill:#2196f3
</pre>
]

**7 external dependencies** for arrays with $N \le 32$ elements.

---

## 🎯 Key Insight

FIR filter coefficients have **at most 32–128 taps**. For arrays this small:

- **BLAS acceleration** provides zero benefit (no cache misses to optimize)
- **Lazy evaluation** adds overhead without measurable gain
- **Template metaprogramming** explodes compile times
- **CMake CPM downloads** make CI unreliable

### Why xtensor?

> "Because we needed `xt::linalg::dot`, `xt::view`, and FFT wrappers."

### Why remove it?

> "Because for $N = 32$, a raw `for` loop beats BLAS. And the CMake complexity isn't worth it."

---

## 📐 The Math: Spectral Factorization

The core algorithm uses **Kolmogorov's spectral factorization method** (1939):

$$
R(\omega) = \sum_{k=0}^{N-1} r_k \cos(k\omega)
$$

$$
\alpha(\omega) = \frac{1}{2} \ln R(\omega), \quad
\phi(\omega) = \mathcal{H}[\alpha(\omega)]
$$

$$
h = \mathcal{F}^{-1}[\, e^{\alpha + j\phi} \,]
$$

Where $N$ is the filter order (typically 32), and $\mathcal{H}$ is the Hilbert transform.

**What xtensor provided**: `cos`, `log`, `exp`, `rfft`, `irfft` — all on tiny arrays.

**What we replaced with**: Simple `for` loops + direct FFTW C API calls.

---

## 🏗️ Solution: Custom `Arr` Class

A ~240-line header that replaces the entire xtensor stack.

### Design

```cpp
class Arr {
    std::vector<double> _data;
    size_t _rows = 0;
    size_t _cols = 0;   // 0 = 1D vector

public:
    // 1D indexing
    double& operator()(size_t i)       { return _data[i]; }
    double& operator[](size_t i)       { return _data[i]; }

// 2D indexing (row-major)
    double& operator()(size_t i, size_t j) {
        return _data[i * _cols + j];
    }

// Arithmetic
    Arr& operator+=(double s);
    Arr& operator*=(double s);
    // ... +, -, *, dot, outer, view, etc.
};
```

- **Flat `std::vector<double>`** — zero abstraction overhead
- **1D and 2D shapes** — enough for FIR filter design
- **All ops are simple loops** — auto-vectorized by the compiler

---

## 📦 What `Arr` Provides

| Operation | xtensor | `Arr` |
|-----------|---------|-------|
| Create | `xt::zeros<double>({n})` | `zeros(n)` |
| Linspace | `xt::linspace<double>(0, 1, n)` | `linspace(0, 1, n)` |
| Range | `xt::range(0, n)` | `Range(0, n)` |
| View | `xt::view(A, row, xt::all())` | `view(A, Range(row), Range(ALL))` |
| Concatenate | `xt::concatenate(xt::xtuple(...), 1)` | `concatenate(a, b, 1)` |
| Dot product | `xt::linalg::dot(A, x)` | `dot(A, x)` |
| Outer product | `xt::linalg::outer(u, v)` | `outer(u, v)` |
| Element-wise | `xt::cos(w)`, `xt::log(R)` | `cos(w)`, `log(R)` |
| Where | `xt::where(w <= wpass)[0]` | `where(w <= wpass)[0]` |
| Sum | `xt::sum(v)()` | `sum(v)` |
| FFT | `xt::fftw::rfft(alpha)` | `rfft(alpha)` (FFTW C API) |

---

## 🏗️ Architecture: Where `Arr` Lives

.mermaid[
<pre>
graph TD
    subgraph "ellalgo-cpp"
        EA[ell.hpp / cutting_plane.hpp]
        A[arr.hpp<br/>~240 lines]
    end
    subgraph "multiplierless-cpp"
        M[LowpassOracle / spectral_fact]
        F[fftw_helper.hpp<br/>~80 lines]
    end
    subgraph "corr-solver-cpp"
        C[LsqOracle / MleOracle]
        L[linalg.hpp<br/>Cholesky, inv, etc.]
    end

M -->|"#include <ellalgo/arr.hpp>"| A
    M --> F
    F --> A
    C --> A
    C --> L

style A fill:#4caf50
    style F fill:#ff9800
    style L fill:#2196f3
    style EA fill:#9c27b0
</pre>
]

**Key decision**: `Arr` lives in `ellalgo-cpp` because it's the template parameter for `Ell<Arr>`.

---

## 📊 Results: Projects Converted

### 1️⃣ `multiplierless-cpp` (FIR Filter Design)

| Metric | Before | After |
|--------|--------|-------|
| Dependencies | 7+ (xtl, xtensor, xtensor-blas, xtensor-fftw, BLAS, LAPACK, FFTW) | **1** (FFTW) |
| CMake configure time | ~91s (CPM downloads) | ~49s |
| Build artifacts | ~440 lines `arr.hpp` + xtensor templates | ~440 lines total |
| Tests | 13/13 + 1000 rapidcheck | ✅ Same |

### 2️⃣ `ellalgo-cpp` (Ellipsoid Method)

| Metric | Before | After |
|--------|--------|-------|
| `Arr` type | xtensor (external) | **Built-in** `arr.hpp` |
| Tests | 102/102 | **128/128** (+26 new Arr tests) |

### 3️⃣ `corr-solver-cpp` (Correlation Solver)

| Metric | Before | After |
|--------|--------|-------|
| Dependencies | xtl, xtensor, xtensor-blas, BLAS | **0 array deps** + custom `linalg.hpp` |
| Tests | 3/3 | ✅ Same |

---

## ⚡ Performance: Build Time Comparison

### Before (with xtensor)

```text
CMake configure:     ~91s  (CPM: xtl, xtensor, xtensor-blas, xtensor-fftw, ...)
Library compile:     ~45s  (template-heavy xtensor headers)
Test compile:        ~15s
Total first build:   ~151s
```

### After (with `Arr`)

```text
CMake configure:     ~49s  (no xtensor stack)
xmake build:         ~12s  (flat loops, no templates)
xmake test:          ~7s   (128 tests + 1500 rapidcheck)
Total:               ~68s
```

**~55% faster first build**, **~73% faster incremental builds**.

---

## 🔧 The Hot-Path Optimization

The original tight loop used xtensor expression templates 20,000+ times:

```cpp
// BEFORE: xtensor expression template (slow to compile, scalar overhead)
auto v = xt::sum(xt::view(A, k, xt::all()) * x)();
```

Each call: 3 template instantiations (view, multiply, sum) → heap-allocated result → scalar extraction.

```cpp
// AFTER: raw for loop (auto-vectorized, zero overhead)
double v = 0.0;
for (size_t j = 0; j < x.size(); ++j)
    v += A(k, j) * x(j);
```

**Result**: Compiler can auto-vectorize. No heap allocations. **~3–5x faster** on the hot path.

---

## 🔬 Case Study: Cosine Matrix Construction

Before — double loop with `std::cos`:

```cpp
Arr An = xt::zeros<double>({m, n - 1});
for (auto i = 0; i != m; ++i)
    for (auto j = 0; j != n - 1; ++j)
        An(i, j) = 2.0 * std::cos(w(i) * (j + 1));
```

After — vectorized outer product + cache:

```cpp
static Arr cached_A;          // cached across iterations!
if (n != cached_n) {
    Arr An = 2.0 * cos(outer(linspace(0, 2*M_PI, m),
                              arange(1.0, double(n))));
    cached_A = concatenate(ones(m, 1), An, 1);
    cached_n = n;
}
```

**Optimizations**: (1) `xt::linalg::outer` vectorized, (2) matrix **cached** per filter order, (3) `mult_factor` reduced from 100 → **20** (5× fewer FFT samples).

---

## 🌲 Dependency Tree: Before vs After

### Before

.mermaid[
<pre>
graph TD
    subgraph "Before"
        M[multiplierless-cpp]
        XB[xtensor-blas]
        XF[xtensor-fftw]
        XT[xtensor]
        XL[xtl]
        BL[BLAS/LAPACK]
        FW[FFTW]
        OB[OpenBLAS]

M --> XB
        M --> XF
        XB --> XT
        XF --> XT
        XT --> XL
        XB --> BL
        XB --> OB
        M --> FW
    end
    style M fill:#f44336
    style XB fill:#ff9800
    style XF fill:#ff9800
    style XT fill:#ff9800
    style XL fill:#ff9800
    style BL fill:#9c27b0
    style OB fill:#9c27b0
    style FW fill:#4caf50
</pre>
]

### After

.mermaid[
<pre>
graph TD
    subgraph "After"
        M2[multiplierless-cpp]
        E[ellalgo-cpp / arr.hpp]
        FW2[FFTW]
        CS[corr-solver-cpp / linalg.hpp]

M2 --> E
        M2 --> FW2
        CS --> E
    end
    style M2 fill:#4caf50
    style E fill:#2196f3
    style FW2 fill:#4caf50
    style CS fill:#4caf50
</pre>
]

**7 dependencies eliminated.** Only FFTW remains (needed for FFT).

---

## 📈 Lessons Learned

### What Worked Well ✅

- **Flat `std::vector<double>`** is all you need for $N \le 128$
- **Simple loops** beat expression templates for tiny matrices
- **Moving `Arr` to ellalgo-cpp** avoided circular dependency issues
- **Static cache** of cosine matrix reduced spectral factorization time by 5×
- **Direct FFTW C API** is simpler than going through xtensor-fftw

### What Was Tricky ⚠️

- **Constructor ambiguity**: `Arr(size_t, double)` vs `Arr(size_t, size_t)` — had to remove one
- **`xt::where` semantics**: Returns indices differently than simple `where()`
- **`xt::row`/`xt::col` assignments**: Need mutable view support
- **xtensor-blas operations**: Cholesky, inv, meshgrid — had to write from scratch for corr-solver

---

## 🧪 Test Results Summary

```text
ellalgo-cpp:
  Test cases: 128 | 128 passed | 0 failed
  RapidCheck: 1500/1500 passed

multiplierless-cpp:
  Test cases: 13 | 13 passed | 0 failed  
  RapidCheck: 1000/1000 passed

corr-solver-cpp:
  Test cases: 3 | 3 passed | 0 failed
```

All existing tests pass with identical numerical results (within machine epsilon).

---

## 🔮 Future Work

### For `arr.hpp` Improvements

- **Span-based views** (no data copy) for zero-cost slicing
- **Fixed-size small-buffer optimization** for $N \le 32$
- **SSE/AVX intrinsics** for dot product hot paths

---

## Q&A 🎤

---

# 🙏 Thank You

### Links

- 📦 [multiplierless-cpp](https://github.com/luk036/multiplierless-cpp)
- 📦 [ellalgo-cpp](https://github.com/luk036/ellalgo-cpp)
- 📦 [corr-solver-cpp](https://github.com/luk036/corr-solver-cpp)

**Key takeaway**: For small arrays, a **240-line header** can replace a **7-library dependency stack**.