Eigen: álgebra lineal en C++

Hasta ahora has escrito bucles para sumar y multiplicar. En álgebra lineal eso se vuelve insostenible: una multiplicación de matrices de 1000 por 1000 son mil millones de operaciones y nadie quiere escribirlas a mano. Eigen es una biblioteca de C++ que ofrece matrices, vectores y descomposiciones con una sintaxis parecida a la de NumPy, pero compilada.

¿Qué es Eigen?

Eigen es una biblioteca header-only: no se compila por separado ni se enlaza, solo se incluyen sus cabeceras. Está escrita en C++ con plantillas, de modo que el compilador genera código específico para cada tipo de matriz y cada expresión. Es la base de muchas bibliotecas de robótica y visión por computador, y una de las formas más directas de hacer álgebra lineal seria en C++.

• Tipos con tamaño en el tipo: Matrix3d es una 3 por 3 de double; MatrixXd es de tamaño dinámico.
• Expresiones diferidas: A + B no calcula nada hasta que se asigna a una matriz, y entonces el compilador fusiona las operaciones.
• Descomposiciones incluidas: LU, QR, Cholesky, SVD y autovalores, cada una con su método solve().
• Sin dependencias: solo cabeceras de la biblioteca estándar, así que se integra en cualquier proyecto.

¿Para qué sirve en Machine Learning?

Casi todo modelo lineal se reduce a una operación matricial: la predicción es X * w, la pérdida es una norma, el ajuste por mínimos cuadrados es un sistema sobredeterminado y la regularización añade una matriz identidad. Con Eigen escribes esas ideas tal como las piensas, y el compilador se encarga del bucle.

• Regresión lineal: resolver X w = y en el sentido de mínimos cuadrados.
• PCA: autovalores y autovectores de la matriz de covarianza.
• Sistemas de recomendación: factorización de matrices de baja dimensión.
• Robótica y visión: transformaciones rígidas, calibración de cámaras, filtros de Kalman.
• Prototipos antes de GPU: comprobar el algoritmo en CPU con dobles antes de portarlo a CUDA.

Cómo funciona por dentro: expresiones y almacenamiento

Eigen no interpreta tus expresiones en tiempo de ejecución: las convierte en tipos. Cuando escribes A * B + C, el compilador construye un objeto que representa esa operación y solo la evalúa al asignarla. Gracias a eso puede recorrer las matrices una vez y sumar el producto en el mismo bucle, sin crear temporales intermedios.

El almacenamiento por defecto es por columnas, igual que en Fortran y BLAS: los elementos de una columna son contiguos en memoria. Esa decisión explica por qué recorrer una matriz de Eigen por filas puede ser más lento que por columnas, el mismo fenómeno de localidad que viste en la lección de rendimiento. Para la mayoría de usos no importa, pero sí cuando escribes bucles propios sobre A(i, j).

Los índices empiezan en cero y se escriben con paréntesis: A(0, 1) es la primera fila, segunda columna. A.rows(), A.cols() y A.size() dan las dimensiones, y A.transpose() devuelve una vista, no una copia.

Sintaxis: instalar, incluir y declarar

Antes de escribir una línea, comprueba que las cabeceras están donde el compilador las busca. Eigen no se enlaza: basta con indicar la carpeta que contiene el directorio Eigen.

1234567891011

# Debian y Ubuntu
sudo apt install libeigen3-dev

# macOS con Homebrew (queda en /opt/homebrew/include/eigen3)
brew install eigen

# Compilar indicando la carpeta que contiene el directorio Eigen
g++ -std=c++20 -O2 -Wall -Wextra -I/usr/include/eigen3 algebra.cpp -o algebra

# En macOS con Homebrew
g++ -std=c++20 -O2 -Wall -Wextra -I/opt/homebrew/include/eigen3 algebra.cpp -o algebra

Con las cabeceras localizadas, un programa mínimo declara una matriz, la rellena y opera con ella. Fíjate en que el tamaño se pasa al constructor porque el tipo es dinámico.

12345678910111213141516

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

int main() {
    Eigen::MatrixXd A(2, 2);      // dinamica: el tamano se fija al construir
    A << 1.0, 2.0,
         3.0, 4.0;

    Eigen::VectorXd v(2);
    v << 5.0, 6.0;

    std::cout << "A * v =\n" << A * v << '\n';
    std::cout << "traspuesta =\n" << A.transpose() << '\n';
    std::cout << "determinante = " << A.determinant() << '\n';
    return 0;
}

El operador << sirve para rellenar: los valores se leen por filas, separados por comas. A * v da [17, 39] y el determinante es -2, porque 1·4 − 2·3. No hace falta ninguna función matmul: el producto es el operador *.

Ejemplo intermedio: multiplicar, trasponer y resolver

Tres operaciones que usarás a diario: un producto de matrices de distinta forma, una trasposición y la resolución de un sistema. Fíjate en que el producto comprueba las dimensiones en tiempo de compilación cuando los tamaños son fijos.

123456789101112131415161718192021222324252627

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

int main() {
    Eigen::MatrixXd A(2, 3);
    A << 1, 2, 3,
         4, 5, 6;

    Eigen::MatrixXd B(3, 2);
    B << 7, 8,
         9, 10,
         11, 12;

    const Eigen::MatrixXd C = A * B;   // (2x3)(3x2) = 2x2
    std::cout << "C =\n" << C << '\n';

    Eigen::Matrix2d S;                 // tamano fijo, sin reservar memoria
    S << 2, -1,
        -1,  3;

    Eigen::Vector2d b;
    b << 1, 2;

    const Eigen::Vector2d x = S.ldlt().solve(b);
    std::cout << "x = " << x.transpose() << '\n';
    return 0;
}

C es [[58, 64], [139, 154]]: cada elemento es el producto de una fila de A por una columna de B. El sistema 2x − y = 1, −x + 3y = 2 tiene solución exacta x = 1, y = 1, y ldlt() la obtiene sin invertir la matriz. La matriz es simétrica y definida positiva, así que Cholesky es la opción más rápida.

Visualización: una multiplicación paso a paso

Recorre los pasos y observa de dónde sale cada celda de C. La regla es siempre la misma: la fila de A que marca el primer índice se multiplica escalarmente por la columna de B que marca el segundo.

El coste de esta operación es el que mediste en la lección de complejidad: para una matriz m×k por otra k×n se hacen m·k·n multiplicaciones. Con 1000 por 1000 son mil millones, y ahí es donde Eigen marca la diferencia frente a un triple bucle escrito a mano: usa bloques que caben en caché y, si el procesador lo permite, instrucciones SIMD.

Caso práctico: mínimos cuadrados en tres líneas

Tienes cuatro muestras y sospechas que y depende linealmente de x. Añades una columna de unos para el término independiente, con lo que el modelo es y = b + w·x. El sistema tiene 4 ecuaciones y 2 incógnitas: no tiene solución exacta en general, pero sí la mejor aproximación.

123456789101112131415161718

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

int main() {
    Eigen::MatrixXd X(4, 2);   // columna 0: unos; columna 1: x
    X << 1, 1,
         1, 2,
         1, 3,
         1, 4;

    Eigen::VectorXd y(4);
    y << 2, 4, 6, 8;

    const Eigen::VectorXd w = X.colPivHouseholderQr().solve(y);
    std::cout << "b = " << w(0) << "  w = " << w(1) << '\n';
    std::cout << "residuo = " << (X * w - y).norm() << '\n';
    return 0;
}

El resultado es b = 0, w = 2 con un residuo del orden de 1e-15: en binario pueden aparecer valores como -4e-16 en lugar de cero, porque la solución se calcula con coma flotante. Los cuatro puntos están exactamente sobre la recta y = 2x. colPivHouseholderQr() resuelve el sistema en el sentido de mínimos cuadrados cuando hay más ecuaciones que incógnitas, y lo hace sin formar las ecuaciones normales, que empeoran el condicionamiento numérico.

Descomposiciones: cuál elegir

Todas las descomposiciones de Eigen tienen un método solve() con la misma firma, así que cambiar de estrategia es cambiar una palabra. La elección depende de lo que sepas de la matriz.

La regla práctica: si sabes que la matriz es simétrica y definida positiva, usa llt(). Si es cuadrada e invertible y no sabes más, partialPivLu(). Si hay más ecuaciones que incógnitas, QR con pivoteo o SVD. Y si la matriz puede ser casi singular, paga el precio del SVD: te dirá qué columnas son redundantes.

Errores comunes

• Olvidar el -I: el error de cabecera no encontrada casi siempre es una ruta mal escrita.
• Usar auto con una expresión: auto C = A + B guarda la expresión, no el resultado. Si A o B cambian después, C cambia contigo. Escribe Eigen::MatrixXd C = A + B;.
• Aliasing en la asignación: A = A.transpose() da un resultado incorrecto. Usa A.transposeInPlace() o evalúa antes con .eval().
• Confundir = con <<: A = 1, 2, 3, 4 no rellena la matriz; el relleno es A << 1, 2, 3, 4 y debe coincidir con el número de elementos.
• Mezclar tamaños: multiplicar una 2 por 3 con una 2 por 3 no compila si los tamaños son fijos, y aborta en tiempo de ejecución con aserción si son dinámicos.

Buenas prácticas y rendimiento

• Compila con -O2 o -O3: sin optimización, Eigen pierde casi toda su ventaja porque no vectoriza.
• Usa tipos fijos (Matrix3d, Vector4f) para tamaños de 4 o menos; el compilador los guarda en registros.
• Deja el tamaño dinámico (MatrixXd) para datos cuyo tamaño solo conoces al leerlos.
• En C = A * B con C ya dimensionada, escribe C.noalias() = A * B; para evitar un temporal.
• Prefiere solve() a inverse(); y si solo necesitas el determinante, determinant() es más barato que invertir.
• Convierte a arrays con .array() para operaciones elemento a elemento: A.array() * B.array().

Ejercicios

Eigen::Matrix3d A;
A << 2, 1, 0,
     1, 3, 1,
     0, 1, 4;
Eigen::Vector3d b;
b << 3, 5, 5;

const Eigen::Vector3d x = A.().solve(b);
std::cout << A.() << '\n';   // la matriz con filas y columnas intercambiadas

Soluciones explicadas

En el primer ejercicio, comprobar el residuo vale más que mirar la solución: si A * x - b es del orden de 1e-16, la solución es correcta aunque no sea exactamente [1, 1, 1] en binario. En coma flotante, «exacto» significa «con el error de redondeo esperado», y ese error se mide, no se supone.

En el segundo, la matriz X tiene más filas que columnas porque hay cinco observaciones y dos parámetros. Ese desajuste es la definición de un problema de mínimos cuadrados: no buscas una solución exacta, buscas la que deja el residuo más pequeño. Cuando los datos son perfectamente lineales, esa mejor solución resulta tener residuo cero y coincide con la recta que ya conocías.

Resumen

• Eigen es header-only: se incluye con #include <Eigen/Dense> y se compila con -I a su carpeta.
• MatrixXd es dinámica; Matrix3d y Vector3d tienen el tamaño en el tipo y son más rápidas.
• El producto es A * B, la traspuesta A.transpose() y los elementos se leen con A(i, j).
• Para resolver usa una descomposición con solve(), no inverse().
• llt() para matrices simétricas definidas positivas, LU para cuadradas, QR o SVD para mínimos cuadrados.
• El almacenamiento por defecto es por columnas, lo que afecta a los bucles que escribas a mano.
• auto con una expresión guarda la expresión; declara el tipo de la matriz para forzar el cálculo.

Evaluación

12345678

import numpy as np

A = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
B = np.array([[7.0, 8.0], [9.0, 10.0], [11.0, 12.0]])

C = A @ B
print(C)
print(np.linalg.solve(np.array([[2.0, -1.0], [-1.0, 3.0]]), np.array([1.0, 2.0])))

123456789

Eigen::MatrixXd A(2, 3);
A << 1, 2, 3,
     4, 5, 6;
Eigen::MatrixXd B(3, 2);
B << 7, 8,
     9, 10,
     11, 12;

const Eigen::MatrixXd C = A * B;