01 — Eigen: Linear Algebra in C++

Eigen is the linear algebra library you reach for when you want NumPy-like ergonomics with C++ speed and no runtime dependencies. It is header-only, template-heavy, and uses expression templates so that an expression like y = A * x + b * c compiles to a single fused loop with zero temporary allocations. It ships with, or wraps around, BLAS/LAPACK when you need them, and it is what TensorFlow, Ceres, MediaPipe, g2o, and countless robotics/perception codebases use internally. In 2026, Eigen 3.4+ is still the dominant C++ template linear algebra library — Armadillo and xtensor exist and are worth knowing exist, but Eigen has the community, the docs, the interop, and the perf.

You will use Eigen for two things in Phase 5: (1) as the numerical core of your from-scratch MLP (P5.1), and (2) as a stepping stone to LibTorch — the tensor model is similar enough that fluency in Eigen makes LibTorch’s Tensor feel like a familiar object with more machinery around it.

1. Install

# macOS (Homebrew):
brew install eigen                          # header-only, installs to /opt/homebrew/include/eigen3

# Ubuntu:
sudo apt install libeigen3-dev

# vcpkg (recommended for your Phase 4 template):
# In vcpkg.json:
#   "dependencies": ["eigen3"]

# Conan:
# In conanfile.txt: eigen/3.4.0

CMake:

find_package(Eigen3 3.4 REQUIRED CONFIG)
target_link_libraries(mynn PUBLIC Eigen3::Eigen)

That’s it. Header-only, so no libraries to link and no ABI to worry about.

2. The core types

Eigen has two type families:

• Matrix<Scalar, Rows, Cols> — linear algebra semantics. * is matrix multiply, .transpose(), .inverse(), etc.

• Array<Scalar, Rows, Cols> — coefficient-wise semantics. * is element-wise. Use for NumPy-style broadcasting.

Convert between them with .array() and .matrix() — same underlying storage, different view.

Typedefs you’ll use daily

using Eigen::MatrixXd;   // dynamic × dynamic, double
using Eigen::MatrixXf;   // dynamic × dynamic, float
using Eigen::VectorXd;   // column vector, double, dynamic length
using Eigen::VectorXf;   // column vector, float
using Eigen::RowVectorXf;
using Eigen::Matrix3d;   // 3×3, double, stack-allocated (fixed size)
using Eigen::Vector4f;   // length 4, stack-allocated
using Eigen::ArrayXf;    // 1D array, float, coefficient-wise
using Eigen::ArrayXXf;   // 2D array

Fixed-size vs dynamic: if you know a dimension at compile time and it is small (≤ 16), use the fixed-size type. Eigen unrolls loops and skips the heap allocation. For ML shapes (batch × features), dimensions are dynamic — use MatrixXf.

First examples

#include <Eigen/Dense>
#include <iostream>

int main() {
    Eigen::MatrixXf W = Eigen::MatrixXf::Random(4, 3);   // 4×3 random in [-1, 1]
    Eigen::VectorXf x = Eigen::VectorXf::Ones(3);        // length-3 all ones
    Eigen::VectorXf b = Eigen::VectorXf::Zero(4);

    Eigen::VectorXf y = W * x + b;                       // matrix-vector + broadcast add
    std::cout << "y =\n" << y << "\n";

    Eigen::MatrixXf A = Eigen::MatrixXf::Random(3, 3);
    Eigen::VectorXf rhs = Eigen::VectorXf::Random(3);
    Eigen::VectorXf sol = A.colPivHouseholderQr().solve(rhs);   // solve A x = rhs
}

3. Expression templates — the reason Eigen is fast

When you write y = W * x + b, a naive library would:

1. Compute W * x into a temporary vector.

2. Add b into a second temporary vector.

3. Copy that into y.

Two heap allocations, two loops, unnecessary memory traffic.

Eigen instead builds an expression tree at compile time:

Assign(
  y,
  Add(
    MatMul(W, x),
    b
  )
)

When you finally assign to y, Eigen evaluates the whole tree in one fused loop, writing directly into y’s storage, with zero temporaries. This is the expression template pattern. It is why Ax + b + c + d runs at hand-written speed.

You get this automatically. You do not need to think about it. But two consequences:

1. auto is a footgun. auto z = W * x + b; binds z to the expression object, not to a materialized vector. If you then modify W, z changes with it. Rule: assign to a concrete type (VectorXf) unless you know what you’re doing.

2. Aliasing bugs. A = A.transpose(); is a bug — Eigen tries to write to A while still reading from it. Use A.transposeInPlace() or A = A.transpose().eval().

4. Aliasing and noalias()

For matrix multiplication C = A * B, Eigen assumes potential aliasing and creates a hidden temporary for safety. If you know the LHS doesn’t overlap with the RHS:

C.noalias() = A * B;    // skip the temporary — direct write into C

This is a real speedup on tight inner loops (10–30% on matmul-dominated code). Use it when you can guarantee no aliasing. If you’re wrong, you get silent numerical corruption — so be sure.

5. Storage order — Eigen is column-major by default

Eigen defaults to column-major storage. NumPy is row-major. This matters when:

• You memory-map data from a NumPy binary file.

• You interop with LibTorch, ONNX Runtime, or Arrow — all row-major.

• You read/write to raw pointer buffers.

Override per-type:

using RowMatrixXf = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;

For ML interop, define row-major typedefs project-wide and use those. This avoids a whole category of “the shapes look right but the numbers are wrong” bugs.

6. Map<> — zero-copy interop with raw pointers

When PyTorch, NumPy, Arrow, or ONNX Runtime hands you a float* and a shape, you don’t want to copy the data into Eigen. Use Map:

float* data = /* from pybind11 numpy buffer, torch tensor, or arrow array */;
int rows = 32, cols = 128;

// Zero-copy view onto external memory. No allocation.
Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
    view(data, rows, cols);

// Now use `view` like any MatrixXf — .rows(), .cols(), multiplication, everything works.
Eigen::VectorXf col_sums = view.colwise().sum();

This is the mechanism for hooking Eigen into the ML stack. You will use it every time you write pybind11 code that accepts a numpy array.

7. Interop with BLAS / LAPACK

Eigen ships built-in implementations of most operations. For large matrices (roughly n > 500), a hand-tuned BLAS (OpenBLAS, Intel MKL, Apple’s Accelerate) beats Eigen’s built-in. To route Eigen through BLAS:

#define EIGEN_USE_BLAS
#define EIGEN_USE_LAPACKE     // for LAPACK routines: SVD, QR, etc.
#include <Eigen/Dense>

Then link the BLAS lib. On macOS you get Apple Accelerate for free:

if(APPLE)
    target_link_libraries(mynn PUBLIC "-framework Accelerate")
    target_compile_definitions(mynn PUBLIC EIGEN_USE_BLAS EIGEN_USE_LAPACKE)
endif()

Benchmark before and after. On Apple Silicon with Accelerate, expect a 2–5x speedup on MatrixXf * MatrixXf at n=1024 vs Eigen’s default kernels. If your matrices are ML-shaped (batch × 128 × 768 kinda thing), the crossover is worth measuring for your workload.

8. Common operations cheat sheet

// Creation
MatrixXf::Zero(3, 4);              MatrixXf::Ones(3, 4);          MatrixXf::Identity(3, 3);
MatrixXf::Random(3, 4);            MatrixXf::Constant(3, 4, 2.5f);

// Shape
A.rows();  A.cols();  A.size();     A.resize(m, n);  // .size() = rows*cols

// Slicing (Eigen 3.4+)
A.row(i);                          A.col(j);
A.block(row0, col0, nrows, ncols); A(seq(0, 5), all);   // Eigen 3.4 seq/all
A.leftCols(2);  A.topRows(3);

// Reductions
A.sum();   A.mean();   A.maxCoeff();   A.minCoeff();
A.rowwise().sum();     A.colwise().mean();

// Coefficient-wise (use .array())
(A.array() * B.array()).matrix();     // element-wise product
(A.array().exp()).matrix();           // element-wise exp
(A.array() > 0).cast<float>();        // ReLU mask

// Matrix ops
A * B;                               A.transpose();
A.inverse();                         // 4x4 max — for bigger use a decomposition
A.determinant();                     A.trace();

// Decompositions
Eigen::LLT<MatrixXf> llt(A);         // Cholesky, positive-definite
Eigen::LDLT<MatrixXf> ldlt(A);       // Cholesky, semi-definite
Eigen::PartialPivLU<MatrixXf> lu(A); // General
Eigen::JacobiSVD<MatrixXf> svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);

9. Building your MLP’s building blocks (preview for P5.1)

A single linear layer forward pass in Eigen:

struct Linear {
    Eigen::MatrixXf W;   // [out_features, in_features], row-major recommended for ML
    Eigen::VectorXf b;   // [out_features]

    Linear(int in_f, int out_f)
        : W(Eigen::MatrixXf::Random(out_f, in_f) * std::sqrt(2.f / in_f)),  // He init
          b(Eigen::VectorXf::Zero(out_f)) {}

    // x: [batch, in_features]
    // returns: [batch, out_features]
    Eigen::MatrixXf forward(const Eigen::MatrixXf& x) const {
        // x @ W.T + b (broadcast)
        return (x * W.transpose()).rowwise() + b.transpose();
    }
};

Backprop:

// grad_out: [batch, out_features]
struct LinearGrad { Eigen::MatrixXf dW; Eigen::VectorXf db; Eigen::MatrixXf dx; };

LinearGrad backward(const Linear& l, const Eigen::MatrixXf& x, const Eigen::MatrixXf& grad_out) {
    LinearGrad g;
    g.dW = grad_out.transpose() * x;         // [out, in]
    g.db = grad_out.colwise().sum();         // [out]
    g.dx = grad_out * l.W;                   // [batch, in]
    return g;
}

Read these until the shapes are obvious. If you can’t derive them on paper, you can’t implement them in Eigen.

10. Eigen vs alternatives (2026)

Library	Style	When to use
Eigen 3.4+	Template, expression-tree, header-only	Default choice. Robotics, small-to-mid ML, general numerical work.
Armadillo 12+	MATLAB-like syntax (A * B, .t())	Migrating MATLAB code to C++ verbatim. Cleaner for stats; smaller ecosystem.
xtensor / xtensor-python	NumPy-like syntax and broadcasting	You want NumPy-in-C++ literally. Fewer optimizations than Eigen; good pybind11 integration.
Blaze	Fastest of the C++ template libraries in some benchmarks	Small community. Consider only if benchmarks show a clear win for your workload.
BLAS / LAPACK directly (OpenBLAS / MKL / Accelerate)	C API, raw pointers	You are writing performance-critical kernels and don’t want any abstraction. Verbose.
cuBLAS / rocBLAS	GPU BLAS	You are on NVIDIA / AMD GPUs and everything above is CPU.

Verdict: Learn Eigen deeply. Know that Armadillo and xtensor exist for reading unfamiliar code. Use BLAS through Eigen (EIGEN_USE_BLAS), not directly, unless you have a strong reason.

11. What most people get wrong

• auto on Eigen expressions. As above — you bind to the tree, not the value. Eigen::VectorXf y = W * x + b; — always concrete on the LHS.

• Ignoring storage order for interop. They copy data instead of using Map<>, or they use Map<> with the wrong storage order and get garbage. Always match the source’s row/column-major layout.

• A.inverse() * b to solve Ax = b. Slow and numerically bad. Use A.colPivHouseholderQr().solve(b) or A.llt().solve(b).

• Fixed-size types for ML. Matrix<float, 784, 128> allocates 784×128×4 = 400 KB on the stack. Stack overflow. Use MatrixXf for anything bigger than ~16×16.

• Not using noalias() on tight inner loops. Free 10–30%.

• Rewriting Eigen ops in raw loops because “loops must be faster.” No. Eigen’s inner loops are SIMD-vectorized. Your hand loop is not.

• Not enabling -O3 -march=native -DNDEBUG. In debug mode, Eigen is 10–50x slower because bounds checks and asserts fire on every access. Never benchmark in debug mode.

12. Practice exercises

1. Solve Ax = b for a random 500×500 A using colPivHouseholderQr. Compare timing to A.inverse() * b. Report both.

2. Implement one linear layer forward pass. Verify it matches a torch.nn.Linear with the same weights on the same input to within 1e-5.

3. Take a numpy array via a raw float* (mock it: std::vector<float> flat(rows*cols);) and build a Map<RowMatrixXf> view. Compute .colwise().sum(). Verify against a hand loop.

4. Enable Apple Accelerate via EIGEN_USE_BLAS. Benchmark 1024×1024 matmul before/after. Report the ratio.

5. Read the “Aliasing” and “Eigen expressions” sections of the Eigen docs end-to-end. This is 45 minutes and prevents an entire category of bug.

---

Nav: ← Phase 5 README · Next: 02 PyTorch C++ / LibTorch →