C++ Eigen

  • Description: Eigen linear algebra — matrices, storage order, decompositions and solvers, geometry (rotations/quaternions/axis-angle for SLAM, wxyz vs xyzw, homogeneous coords), aliasing and the auto-with-expression-templates pitfall
  • My Notion Note ID: K2A-B2-5
  • Created: 2020-01-17
  • Updated: 2026-04-30
  • License: Reuse is very welcome. Please credit Yu Zhang and link back to the original on yuzhang.io

Table of Contents

1. Setup and Expression Templates

  • Header-only. #include <Eigen/Dense> (or <Eigen/Geometry>, <Eigen/Sparse>). Link nothing.
  • Provides dense + sparse matrices, decompositions, geometric primitives.
  • Internally uses expression templates — A + B * C produces expression tree, compiler evaluates lazily, fuses ops, skips intermediates.
  • Makes Eigen fast + creates the auto and aliasing pitfalls in § 8 and § 9.

Build flags:

  • -O2 or -O3 + vector-ISA target (-march=native, -mavx2, -mfma) for SIMD.
  • -DEIGEN_DONT_PARALLELIZE if you parallelize higher up (TBB, OpenMP, custom pool) and don't want Eigen using OpenMP internally.

2. Matrix and Vector Types

  • Fundamental template: Matrix<Scalar, Rows, Cols>.
Type Equivalent
Vector3d Matrix<double, 3, 1>
Vector3f Matrix<float, 3, 1>
RowVector3d Matrix<double, 1, 3>
Matrix3d Matrix<double, 3, 3>
Matrix4f Matrix<float, 4, 4>
MatrixXd Matrix<double, Dynamic, Dynamic>
VectorXd Matrix<double, Dynamic, 1>
Matrix<double, 3, Dynamic> 3 rows, runtime cols

  • Fixed-size (Matrix3d, Vector4f) — stack-allocated; major source of Eigen's speed (compiler unrolls).

  • Dynamic-size (MatrixXd) — heap-allocated storage.

  • Pick fixed when you can. Vector3d is several × faster to construct, copy, operate than VectorXd of length 3. Dynamic only when size genuinely unknown at compile time.

  • Alignment. SIMD-friendly fixed-size types (Vector4d, Vector2d, Matrix4f) require 16/32-byte alignment.

  • Storing inside your classes — class needs EIGEN_MAKE_ALIGNED_OPERATOR_NEW (pre-C++17) or recent compiler with C++17 over-aligned new (GCC ≥ 7, clang ≥ 5, MSVC ≥ 19.12).

  • Containers — std::vector<Vector4d, Eigen::aligned_allocator<Vector4d>> on older toolchains; modern C++17 handles plain std::vector<Vector4d> correctly.

3. Construction and Access

#include <Eigen/Dense>
using namespace Eigen;

Vector3d v(1, 2, 3);
Matrix3d M;
M << 1, 2, 3,
     4, 5, 6,
     7, 8, 9;

MatrixXd A(3, 4);     // uninitialized
A.setZero();
A.setIdentity();      // for square: identity; otherwise diag of 1s
A.setRandom();
A.setConstant(7.0);

double x = M(0, 1);   // 2.0  -- (row, col)
v(2) = 30.0;          // mutate

// Block access
M.row(0);             // first row as a 1x3 expression
M.col(1);             // second col as a 3x1 expression
M.block<2, 2>(0, 0);  // top-left 2x2 (compile-time size)
M.block(0, 0, 2, 2);  // top-left 2x2 (runtime size)
M.topLeftCorner(2, 2);
v.head<2>();          // first 2 elements
v.tail(2);
v.segment(1, 2);      // 2 elements starting at index 1
  • Indexing always (row, col), regardless of storage order.

4. Arithmetic and Common Operations

Matrix3d A, B;
Vector3d u, v;

Matrix3d C = A + B;
Matrix3d D = A * B;            // matrix product
Vector3d w = A * u;            // matrix-vector product
double dot = u.dot(v);
Vector3d cross = u.cross(v);   // only for size-3 vectors

double n  = v.norm();
double n2 = v.squaredNorm();
Vector3d hat = v.normalized();
v.normalize();                  // in place

// Component-wise (use .array())
Matrix3d E = A.array() * B.array();           // Hadamard product
Matrix3d F = A.array().square();
Matrix3d G = (A.array() > 0).select(A, -A);   // |A|

// Reductions
double s    = A.sum();
double mn   = A.minCoeff();
double mx   = A.maxCoeff();
double tr   = A.trace();
double det  = A.determinant();
Matrix3d In = A.inverse();   // small fixed-size only; use a solver otherwise

// Transpose / adjoint
Matrix3d At = A.transpose();
A.transposeInPlace();        // ok
A = A.transpose();           // BUG -- aliasing; see § 8
  • .array() view → coefficient-wise context (* element-wise, scalar functions work). .matrix() returns. Zero-cost views.

5. Storage Order: Row-Major vs Column-Major

  • Default = column-major. Override per type:
Matrix<double, 3, 3, RowMajor> Mrow;
Matrix<double, 3, 3, ColMajor> Mcol;        // explicit form of the default
typedef Matrix<float, Dynamic, Dynamic, RowMajor> RowMatrixXf;
  • Indexing, transposes, decompositions — behave identically regardless of storage order. Difference matters in 3 places:
  1. Interop with external memory. Eigen::Map wraps raw pointer; layout in memory must match declared storage order. Most numpy buffers, OpenCV Mats, image data = row-major. Map<MatrixXd, 0, Stride<Outer, Inner>> for arbitrary strides.
  2. Slice performance. M.row(i) cheap on row-major, crosses strides on column-major (vice versa for M.col(j)). Pick layout matching access pattern.
  3. Vector typedefs. VectorXd column-major (column vector). RowVectorXd row-major (row vector). Not interchangeable in shape-aware expressions.
// Wrap an existing row-major buffer (e.g., from a Python/C library)
double buf[12];
Map<Matrix<double, 3, 4, RowMajor>> M(buf);
M(1, 2) = 5.0;       // writes through to buf

6. Solving Linear Systems and Decompositions

  • Don't compute A.inverse() * b to solve A x = b. Use decomposition:
Decomposition Matrix shape and properties Use when
PartialPivLU Square, invertible General fast solve, stable enough
FullPivLU Square More stable, also gives rank/null space
HouseholderQR Any Least-squares, fast
ColPivHouseholderQR Any Least-squares, more stable, gives rank
LLT (Cholesky) Square, symmetric positive-definite Fast SPD solver
LDLT Square, symmetric (positive or negative semidefinite) Symmetric semidefinite — handles PSD/NSD cases that strict-PD LLT rejects, no sqrt
JacobiSVD Any Most stable, computes full SVD
BDCSVD Any (large) Faster SVD for large matrices 
MatrixXd A = MatrixXd::Random(3, 3);
VectorXd b = VectorXd::Random(3);

VectorXd x = A.colPivHouseholderQr().solve(b);
double residual = (A * x - b).norm();

// SPD case
MatrixXd K = A.transpose() * A;     // SPD
VectorXd y = K.llt().solve(A.transpose() * b);

// Eigenvalues (self-adjoint case is much faster)
SelfAdjointEigenSolver<Matrix3d> es(K);
Vector3d evals = es.eigenvalues();
Matrix3d evecs = es.eigenvectors();

// SVD
JacobiSVD<MatrixXd> svd(A, ComputeThinU | ComputeThinV);
  • Small fixed-size matrices → A.inverse() OK (special-cases sizes ≤ 4).
  • Bigger or runtime-sized → solver. A.inverse() * b slower (full inverse) + less numerically stable than solver.solve(b).

7. The Geometry Module

  • #include <Eigen/Geometry> adds rotation, translation, rigid-transform types for 3D.
  • Bread-and-butter of SLAM, robotics, CV, graphics.
  • Conventions: Hamilton (right-handed) quaternions + right-handed coordinate frames. What Eigen, ROS, GTSAM, Sophus, most of robotics use.

7.1 Rotation Representations

Type Storage Best for
Matrix3d 9 doubles Composition, rotating vectors (one matvec), interop with linear-algebra code
Quaterniond 4 doubles Smooth interpolation (slerp), no gimbal lock, compact, the right thing to store on disk
AngleAxisd 4 doubles (angle + axis) Human-readable construction, small rotations as a tangent-space step
Euler angles (Vector3d of eulerAngles) 3 doubles Intuition, configs, RPY for humans — avoid internally (gimbal lock, convention chaos) 
// Constructors
AngleAxisd  aa(M_PI / 4, Vector3d::UnitZ());           // 45° about +Z
Quaterniond q(aa);                                     // from axis-angle
Matrix3d    R = q.toRotationMatrix();

// Apply to a point
Vector3d p(1, 0, 0);
Vector3d p_rotated = q * p;                            // overloaded * for unit quats
Vector3d p_rotated_alt = R * p;

// Interpolation between two orientations (slerp)
Quaterniond q0 = Quaterniond::Identity();
Quaterniond q1(aa);
Quaterniond q_mid = q0.slerp(0.5, q1);

7.2 Quaternion Storage: wxyz vs xyzw

  • Catches everyone exactly once. 2 conventions for the component order; Eigen uses one in constructor, the other in internal storage:
Where Order
Quaterniond(w, x, y, z) constructor (and .w(), .x(), .y(), .z() accessors) w, x, y, z — scalar first
Quaterniond::coeffs() (4-vector), and the underlying contiguous memory x, y, z, w — scalar last 
Quaterniond q(0.7071, 0.0, 0.0, 0.7071);   // w=0.7071, x=0, y=0, z=0.7071
std::cout << q.w() << " " << q.x() << " "
          << q.y() << " " << q.z() << '\n'; // 0.7071 0 0 0.7071

Vector4d v = q.coeffs();
// v == (x=0, y=0, z=0.7071, w=0.7071) -- xyzw order!

Matters when serializing / interop:

  • ROS geometry_msgs/Quaternion — serialized (x, y, z, w). Same byte order as q.coeffs(). Build with Quaterniond q(msg.w, msg.x, msg.y, msg.z).

  • glm — glm::quat(w, x, y, z) matches Eigen's constructor.

  • Unity / DirectX / many graphics APIs — (x, y, z, w). Match coeffs().

  • JSON / YAML — write down which order; never pass a 4-tuple unnamed.

  • Axis-angle convention universal: angle first, axis second.

  • Common-sense check after deserializing: unit quaternion has q.norm() == 1; applying to known vector produces expected rotation. If off → almost always order convention.

7.3 Axis-Angle

  • AngleAxisd(angle, axis) — rotation by angle (rad) around unit axis. Right-hand rule (fingers curl in positive rotation direction).
AngleAxisd aa(0.1, Vector3d::UnitZ());     // 0.1 rad ≈ 5.73° about +Z

// Convert to/from
Quaterniond q(aa);
Matrix3d    R = aa.toRotationMatrix();
AngleAxisd  aa_back(R);                    // canonical form: angle in [0, π]

// Compose three axis-aligned rotations explicitly
Quaterniond q_rpy =
    AngleAxisd(yaw,   Vector3d::UnitZ()) *
    AngleAxisd(pitch, Vector3d::UnitY()) *
    AngleAxisd(roll,  Vector3d::UnitX());   // Z-Y-X intrinsic
  • Natural for small rotations (tangent space of SO(3) at identity = axis * angle) → IMU integration, gradient steps in pose optimization, exponential maps:
Vector3d omega_dt = ...;                    // small rotation as 3-vector (axis * angle)
double angle = omega_dt.norm();
Quaterniond dq = (angle < 1e-9)
    ? Quaterniond::Identity()
    : Quaterniond(AngleAxisd(angle, omega_dt / angle));

q = q * dq;                                 // integrate rotation
q.normalize();                              // periodically re-normalize to combat drift

7.4 Conversions Between Representations

Matrix3d    R  = q.toRotationMatrix();
Quaterniond q2(R);                          // matrix → quat
AngleAxisd  aa2(q);                         // quat → axis-angle
AngleAxisd  aa3(R);                         // matrix → axis-angle

// Euler (use sparingly; read the convention!)
Vector3d ypr = R.eulerAngles(2, 1, 0);      // intrinsic Z-Y-X
  • eulerAngles(a, b, c) — (a, b, c) is the axis order. (0, 1, 2) = XYZ, (2, 1, 0) = ZYX. Result is intrinsic, in radians; 2nd-axis angle in [0, π] (so XZX-style "first and third axis the same" decompositions are well-defined).
  • Document the convention or Friday-afternoon debugging awaits.

7.5 Composing Rotations

  • Rotations compose by multiplication, right-to-left (rightmost applied first):
// "rotate first by a, then by b" applied to vector v:
Quaterniond ab  = b * a;
Vector3d    out = ab * v;                   // == b * (a * v)
  • For unit quaternions, inverse = conjugate:
Quaterniond q_inv = q.conjugate();          // == q.inverse() for unit quats
  • q.conjugate() faster than q.inverse() (no division). Agree as long as quaternion is unit. q.normalize() if unsure.

7.6 Rigid Transforms and SLAM Patterns

  • For rotation + translation → Isometry3d (or parent Affine3d). Isometry3d = rigid-body transform group SE(3) — what SLAM calls a "pose."
Isometry3d T_world_body = Isometry3d::Identity();
T_world_body.linear()      = q.toRotationMatrix();    // 3×3 rotation block
T_world_body.translation() = Vector3d(1, 2, 3);

Vector3d p_body  = ...;
Vector3d p_world = T_world_body * p_body;             // applies R then t
Isometry3d T_body_world = T_world_body.inverse();     // exact, not LU

Patterns common in SLAM:

Compose poses along a chain. Common bug: left-to-right composition.

// T_W_C: pose of camera in world; T_C_L: pose of LiDAR relative to camera
Isometry3d T_W_L = T_W_C * T_C_L;             // correct
Vector3d   p_W   = T_W_L * p_lidar;
  • "Frame-on-frame" naming (T_A_B reads "B in A") makes chain rule mechanical — subscripts cancel inwards. T_W_C * T_C_L = T_W_L.

Transforming a point cloud. N points as 3×N matrix:

Matrix<double, 3, Eigen::Dynamic> P_body  = ...;            // 3×N
Matrix<double, 3, Eigen::Dynamic> P_world =
    (T_world_body.linear() * P_body).colwise() + T_world_body.translation();
  • Or with Isometry3d directly per column (clearer, equally fast):
Matrix<double, 3, Eigen::Dynamic> P_world(3, N);
for (int i = 0; i < N; ++i) P_world.col(i) = T_world_body * P_body.col(i);

Homogeneous coordinates — 3D point as 4-vector (x, y, z, 1); affine transforms (rotation + translation) → single 4×4 matrix multiply. Eigen helpers:

Vector3d p(1, 2, 3);
Vector4d p_h = p.homogeneous();              // (1, 2, 3, 1)
Vector3d p_back = p_h.hnormalized();         // (1, 2, 3) (divides by last component)

// Apply a 4x4 transform stored as Matrix4d
Matrix4d M = T_world_body.matrix();
Vector4d p_world_h = M * p.homogeneous();
Vector3d p_world   = p_world_h.hnormalized();
  • Homogeneous form mostly useful for projective transforms (camera projection, NDC) — perspective divide in hnormalized() does work.
  • Pure rigid motion → Isometry3d multiplication cleaner, avoids the divide.

Storing + re-normalizing rotations. Repeated composition accumulates FP drift — result slowly stops being orthonormal (matrices) or unit (quaternions).

q.normalize();                                 // unit quaternion
// or for a rotation matrix: re-orthogonalize via SVD
JacobiSVD<Matrix3d> svd(R, ComputeFullU | ComputeFullV);
Matrix3d R_clean = svd.matrixU() * svd.matrixV().transpose();
  • Long-running SLAM — normalize every few hundred frames or after every optimization step.
  • Quaternions drift much more slowly than matrices → SLAM front-ends keep state as (translation, quaternion) even when computations use matrices.

8. Aliasing and .eval()

  • Eigen evaluates expression trees lazily, writes result coefficient-by-coefficient.
  • Destination overlaps source → can read coefficient already overwritten. Aliasing bug.
// BUG: reading transposed entries from M while writing to M
MatrixXd M = MatrixXd::Random(3, 3);
M = M.transpose();          // wrong result, no warning
  • Most ops — Eigen detects aliasing at compile time + inserts temporary automatically (matrix products do this).
  • Transpose, blocks, reverse — user must opt in:
M = M.transpose().eval();   // forces evaluation into a temporary first
M.transposeInPlace();       // in-place version that handles aliasing

  • Other *InPlace helpers for reverse, adjoint, etc.

  • Run with -DEIGEN_RUNTIME_NO_MALLOC + eigen_assert in tests → catch aliasing issues early.

  • Flip side: when you know no aliasing (most matrix products into fresh dest) → noalias() to skip safety temporary:

C.noalias() = A * B;        // no internal temporary
  • Measurable win for hot loops doing many matrix-matrix mults.

9. The auto Pitfall with Expression Templates

  • Eigen returns expression objects, not matrices. auto captures the expression, not value.
  • Inputs go out of scope or mutated → cached expression points at garbage.
// BUG: auto deduces an expression, not a Matrix.
auto E = A + B;             // E is an Eigen::CwiseBinaryOp<...>, not a Matrix
A.setZero();                // mutating A invalidates anything E referenced
std::cout << E;             // wrong values
  • Same trap with temporaries:
// BUG: temporary B*c is gone after the line, but E refers into it.
auto E = (B * c).transpose();
use(E);                     // dangling expression
  • Fix: never use auto for Eigen expression unless you really want the lazy expression. Store concrete types:
MatrixXd E = A + B;                        // forces evaluation
Vector3d  e = (B * c).transpose();         // ditto
auto      E_lazy = (A + B).eval();         // explicit eval -> Matrix
  • Eigen documents this as the auto keyword pitfall. Treat as project-wide rule.

10. References