Привет, Хабр!
Кто хоть раз пытался работать с матрицами в C++, знает, что это удовольствие сродни написанию своего STL — возможно, но зачем? Eigen — это библиотека, которая избавит вас от ручного управления памятью, оптимизирует вычисления и позволит писать код, похожий на чистую математику. Поэтому в этой статье мы разберем эту прекрасную библиотеку.
Основы работы с Eigen: матрицы, векторы и выражения
Матрицы и векторы
Всё в Eigen основано на шаблонном классе Eigen::Matrix. Будь то динамическая матрица или матрица фиксированного размера. Вот базовый пример:
#include <iostream>
#include <Eigen/Dense>
int main() {
// Фиксированная матрица 2x2 (double)
Eigen::Matrix2d mat;
mat << 1, 2,
3, 4;
// Фиксированный вектор из 2-х элементов
Eigen::Vector2d vec(5, 6);
// Матрично-векторное умножение
Eigen::Vector2d result = mat * vec;
std::cout << "Результат: " << result.transpose() << std::endl;
return 0;
}На первый взгляд, всё выглядит просто, но именно такая лаконичность и заложена в основу производительности.
Выражения и ленивые вычисления
Одной из фич Eigen является ленивость вычислений. Что это значит? Допустим, есть выражение:
Eigen::Matrix3f C = A + B * D;В этом случае библиотека не создаёт лишних временных матриц — все вычисления оптимизируются на этапе компиляции. Если вдруг вам нужно форсировать вычисление, вызывайте .eval(), чтобы материализовать результат.
Eigen::Matrix3f result = (A * B + C).eval();Операции с матрицами
Блоки и подматрицы
Иногда нужно работать не со всей матрицей, а с её частью. Для этого есть методы для выделения блоков, строк и столбцов. Пример:
#include <iostream>
#include <Eigen/Dense>
int main() {
Eigen::Matrix4d mat;
mat << 1, 2, 3, 4,
5, 6, 7, 8,
9, 10, 11, 12,
13, 14, 15, 16;
// Извлекаем блок 2x2 из центральной части матрицы
Eigen::Matrix2d block = mat.block<2,2>(1,1);
std::cout << "Центральный блок:\n" << block << "\n\n";
// Получаем первую строку
Eigen::RowVector4d firstRow = mat.row(0);
std::cout << "Первая строка: " << firstRow << "\n\n";
// Получаем последний столбец
Eigen::Vector4d lastCol = mat.col(3);
std::cout << "Последний столбец:\n" << lastCol << std::endl;
return 0;
}Такое разделение матрицы на блоки позволяет создавать алгоритмы, где можно обрабатывать части матрицы независимо.
Элементные операции с .array()
Иногда требуется не матричное, а элементное умножение. Для этого Eigen имеет метод .array(), который переводит матрицу в массивный режим:
#include <iostream>
#include <Eigen/Dense>
int main() {
Eigen::Matrix3f A;
A << 1, 2, 3,
4, 5, 6,
7, 8, 9;
// Поэлементное умножение на константу
Eigen::Matrix3f B = A.array() * 2.0f;
std::cout << "Поэлементное умножение:\n" << B << std::endl;
// Элементное умножение двух матриц
Eigen::Matrix3f C = A.array() * A.array();
std::cout << "\nКвадраты элементов матрицы A:\n" << C << std::endl;
return 0;
}Удобно, когда нужно работать с математическими функциями над каждым элементом: .sin(), .cos(), .exp() и т. д.
Решение систем линейных уравнений
Теперь перейдём к одной из самых распространённых задач — решению системы линейных уравнений вида Ax = b. В зависимости от свойств матрицы AA мы можем использовать различные разложения.
LU-разложение
Если матрица не обладает специальными свойствами, можно использовать LU‑разложение с частичной перестановкой:
#include <Eigen/Dense>
#include <iostream>
int main() {
Eigen::Matrix3d A;
A << 3, -1, 2,
1, 0, -1,
2, 1, 1;
Eigen::Vector3d b(1, 0, 4);
Eigen::PartialPivLU<Eigen::Matrix3d> lu(A);
if(lu.isInvertible()) {
Eigen::Vector3d x = lu.solve(b);
std::cout << "Решение системы Ax = b:\n" << x << std::endl;
} else {
std::cerr << "Матрица A необратима!" << std::endl;
}
return 0;
}Всегда проверяйте, что матрица обратима.
Cholesky-разложение
Если AA симметричная и положительно определённая, можно применить Cholesky‑разложение, которое работает быстрее:
#include <Eigen/Dense>
#include <iostream>
int main() {
Eigen::Matrix3d A;
A << 4, -2, 1,
-2, 4, -2,
1, -2, 3;
Eigen::Vector3d b(1, 2, 3);
Eigen::LLT<Eigen::Matrix3d> llt(A);
if(llt.info() == Eigen::Success) {
Eigen::Vector3d x = llt.solve(b);
std::cout << "Решение системы Ax = b (Cholesky):\n" << x << std::endl;
} else {
std::cerr << "Ошибка: матрица A не является симметричной положительно определённой!" << std::endl;
}
return 0;
}QR-разложение
Иногда матрица может быть прямоугольной. Для таких случаев QR‑разложение — отличное решение:
#include <Eigen/Dense>
#include <iostream>
int main() {
// Прямоугольная матрица 4x3
Eigen::MatrixXd A(4, 3);
A << 1, 2, 3,
4, 5, 6,
7, 8, 9,
10, 11, 12;
// Вектор b с 4 элементами
Eigen::VectorXd b(4);
b << 1, 2, 3, 4;
// Решаем систему, используя QR-разложение
Eigen::ColPivHouseholderQR<Eigen::MatrixXd> qr(A);
Eigen::VectorXd x = qr.solve(b);
std::cout << "Решение системы Ax = b (QR):\n" << x << std::endl;
return 0;
}При решении переопределённых систем такой метод показывает отличные результаты в плане устойчивости.
Всё о разработке на С++ можно изучить на специализации C++ Developer.
Собственные значения и собственные векторы
Вычисление собственных значений — ещё одна фича Eigen. В зависимости от симметрии матрицы используйте разные методы.
Для симметричных матриц
Если матрица симметричная, используйте SelfAdjointEigenSolver — это быстро и стабильно:
#include <Eigen/Dense>
#include <iostream>
int main() {
Eigen::Matrix3d M;
M << 4, 1, 2,
1, 3, 0,
2, 0, 5;
Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigensolver(M);
if (eigensolver.info() == Eigen::Success) {
std::cout << "Собственные значения:\n" << eigensolver.eigenvalues() << "\n\n";
std::cout << "Собственные векторы:\n" << eigensolver.eigenvectors() << std::endl;
} else {
std::cerr << "Ошибка вычисления собственных значений!" << std::endl;
}
return 0;
}Для общих матриц
Если матрица не симметричная, есть обычный EigenSolver:
#include <Eigen/Dense>
#include <iostream>
int main() {
Eigen::Matrix3d M;
M << 1, 2, 3,
0, 4, 5,
0, 0, 6;
Eigen::EigenSolver<Eigen::Matrix3d> solver(M);
if (solver.info() == Eigen::Success) {
std::cout << "Собственные значения:\n" << solver.eigenvalues() << "\n\n";
std::cout << "Собственные векторы:\n" << solver.eigenvectors() << std::endl;
} else {
std::cerr << "Ошибка вычисления собственных значений!" << std::endl;
}
return 0;
}Вычисление собственных значений может быть вычислительно затратным, так что профилируйте свой код, если система большая.
SVD-разложение
Сингулярное разложение — универсальный инструмент для анализа матриц, будь то обработка изображений, рекомендации или уменьшение размерности данных. В Eigen для этого есть класс JacobiSVD.
#include <Eigen/Dense>
#include <iostream>
int main() {
// Создаём прямоугольную матрицу 4x3
Eigen::MatrixXd M(4, 3);
M << 1, 0, 0,
0, 2, 0,
0, 0, 3,
1, 1, 1;
// Вычисляем SVD с тонкими матрицами U и V
Eigen::JacobiSVD<Eigen::MatrixXd> svd(M, Eigen::ComputeThinU | Eigen::ComputeThinV);
std::cout << "Сингулярные значения:\n" << svd.singularValues() << "\n\n";
std::cout << "Матрица U:\n" << svd.matrixU() << "\n\n";
std::cout << "Матрица V:\n" << svd.matrixV() << std::endl;
return 0;
}Флаги ComputeThinU и ComputeThinV позволяют получить компактное представление разложения.
Работа с маппингом и внешними данными
Иногда данные приходят из сторонних источников, например, из С‑структур или файлов, и хочется обернуть их в удобный интерфейс Eigen. Для этого есть класс Eigen::Map:
#include <Eigen/Dense>
#include <iostream>
int main() {
// Сырой массив данных (например, считанный из файла)
double rawData[] = {1, 2, 3, 4, 5, 6};
// Маппинг массива в матрицу 2x3 (без копирования данных)
Eigen::Map<Eigen::Matrix<double, 2, 3>> mat(rawData);
std::cout << "Матрица из сырого массива:\n" << mat << std::endl;
// Можно менять данные напрямую
mat(0, 0) = 42;
std::cout << "\nПосле изменения первого элемента:\n" << mat << std::endl;
return 0;
}Работа с разреженными матрицами
Eigen умеет работать и с разреженными матрицами. Мастхев, когда размерность системы огромна, а ненулевых элементов мало.
#include <Eigen/Sparse>
#include <iostream>
int main() {
// Определяем разреженную матрицу размером 5x5
Eigen::SparseMatrix<double> spMat(5, 5);
// Заполняем разреженную матрицу: только диагональные элементы ненулевые
for (int i = 0; i < 5; ++i) {
spMat.insert(i, i) = i + 1;
}
// Преобразуем в формат Compressed Sparse Row (CSR) для ускорения операций
spMat.makeCompressed();
// Выводим ненулевые элементы
std::cout << "Ненулевые элементы разреженной матрицы:\n";
for (int k = 0; k < spMat.outerSize(); ++k) {
for (Eigen::SparseMatrix<double>::InnerIterator it(spMat, k); it; ++it) {
std::cout << "(" << it.row() << "," << it.col() << ") = " << it.value() << "\n";
}
}
return 0;
}Работая с разреженными матрицами, всегда помните о том, что правильный выбор структуры данных может сократить время вычислений.
Оптимизация кода с Eigen
Без оптимизации даже самый красивый код может тупить. Пару советов:
Фиксированные размеры, когда возможно.
Если матрица имеет размер, известный на этапе компиляции, объявляйте её какEigen::Matrix<double, 3, 3>, а неEigen::MatrixXd. Компилятор тут может провести множество оптимизаций.Используйте
.eval(), чтобы избежать лишних копирований.
В сложных цепочках вычислений результат может не материализоваться до тех пор, пока вы явно не вызовете.eval().Соблюдайте выравнивание памяти.
Eigen использует векторизацию, так что корректное выравнивание данных может дать солидный прирост производительности. Обычно это делается автоматически, но если вы работаете с внешними данными черезEigen::Map, следите за этим.
Пример, объединяющий несколько советов:
#include <Eigen/Dense>
#include <iostream>
int main() {
// Фиксированная матрица 3x3, оптимизированная под векторизацию
Eigen::Matrix3d A;
A << 4, -2, 1,
-2, 4, -2,
1, -2, 3;
// Вектор правой части
Eigen::Vector3d b(1, 2, 3);
// Решаем систему, используя Cholesky (быстро и эффективно)
Eigen::LLT<Eigen::Matrix3d> llt(A);
if(llt.info() == Eigen::Success) {
// Используем eval() для форсированного вычисления результата
Eigen::Vector3d x = llt.solve(b).eval();
std::cout << "Оптимизированное решение системы:\n" << x << std::endl;
} else {
std::cerr << "Ошибка: матрица не подходит для Cholesky-разложения!" << std::endl;
}
return 0;
}Полезные мелочи
Транспонирование и обращение матрицы
Простые операции, но часто используемые:
#include <Eigen/Dense>
#include <iostream>
int main() {
Eigen::Matrix2d mat;
mat << 1, 2,
3, 4;
// Транспонирование матрицы
Eigen::Matrix2d matT = mat.transpose();
std::cout << "Транспонированная матрица:\n" << matT << "\n\n";
// Обращение матрицы (если обратима)
if(mat.determinant() != 0) {
Eigen::Matrix2d invMat = mat.inverse();
std::cout << "Обратная матрица:\n" << invMat << std::endl;
} else {
std::cout << "Матрица необратима!" << std::endl;
}
return 0;
}Вычисление детерминанта
Бывает, что детерминант нужен для оценки устойчивости системы:
#include <Eigen/Dense>
#include <iostream>
int main() {
Eigen::Matrix3d mat;
mat << 1, 2, 3,
0, 1, 4,
5, 6, 0;
double det = mat.determinant();
std::cout << "Детерминант матрицы: " << det << std::endl;
return 0;
}Использование функций высшего порядка
Иногда хочется применять нестандартные операции ко всем элементам матрицы. Eigen позволяет это делать через метод .unaryExpr(), который принимает лямбда‑функцию:
#include <Eigen/Dense>
#include <iostream>
#include <cmath>
int main() {
Eigen::MatrixXd mat(3,3);
mat << 1, 4, 9,
16, 25, 36,
49, 64, 81;
// Применяем квадратный корень к каждому элементу
Eigen::MatrixXd sqrtMat = mat.unaryExpr([](double x) { return std::sqrt(x); });
std::cout << "Матрица после sqrt:\n" << sqrtMat << std::endl;
return 0;
}Интеграция с другими библиотеками
Eigen можно интегрировать с другими библиотеками, такими как OpenCV, Boost и TensorFlow (да‑да, иногда приходится связывать C++ и Python). Например, работа с изображениями в OpenCV может быть усилена за счёт Eigen при выполнении математических операций:
#include <opencv2/opencv.hpp>
#include <Eigen/Dense>
#include <iostream>
int main() {
// Загрузим изображение в OpenCV (grayscale)
cv::Mat image = cv::imread("example.jpg", cv::IMREAD_GRAYSCALE);
if(image.empty()) {
std::cerr << "Не удалось загрузить изображение!" << std::endl;
return -1;
}
// Преобразуем cv::Mat в Eigen::MatrixXd
Eigen::MatrixXd eigenImage(image.rows, image.cols);
for(int i = 0; i < image.rows; ++i)
for(int j = 0; j < image.cols; ++j)
eigenImage(i, j) = static_cast<double>(image.at<uchar>(i, j));
// Применим простейшую нормализацию (можно заменить на более продвинутую операцию)
eigenImage = eigenImage.array() / 255.0;
std::cout << "Нормализованное изображение (первые 5 строк):\n" << eigenImage.topRows(5) << std::endl;
return 0;
}Это лишь малая часть возможностей.
Eigen — это не просто набор функций для работы с матрицами, это целый экосистемный подход к математическим вычислениям в C++. А как вы применяете эту библиотеку в своих проектах?
Пользуясь случаем, рекомендую обратить внимание на открытые уроки по C++: