Линейная алгебра в C++ с Eigen
Привет, Хабр!
Кто хоть раз пытался работать с матрицами в C++, знает, что это удовольствие сродни написанию своего STL — возможно, но зачем? Eigen — это библиотека, которая избавит вас от ручного управления памятью, оптимизирует вычисления и позволит писать код, похожий на чистую математику. Поэтому в этой статье мы разберем эту прекрасную библиотеку.
Основы работы с Eigen: матрицы, векторы и выражения
Матрицы и векторы
Всё в Eigen основано на шаблонном классе Eigen::Matrix. Будь то динамическая матрица или матрица фиксированного размера. Вот базовый пример:
#include
#include
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
#include
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
#include
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
#include
int main() {
Eigen::Matrix3d A;
A << 3, -1, 2,
1, 0, -1,
2, 1, 1;
Eigen::Vector3d b(1, 0, 4);
Eigen::PartialPivLU 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
#include
int main() {
Eigen::Matrix3d A;
A << 4, -2, 1,
-2, 4, -2,
1, -2, 3;
Eigen::Vector3d b(1, 2, 3);
Eigen::LLT 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
#include
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 qr(A);
Eigen::VectorXd x = qr.solve(b);
std::cout << "Решение системы Ax = b (QR):\n" << x << std::endl;
return 0;
} При решении переопределённых систем такой метод показывает отличные результаты в плане устойчивости.
Всё о разработке на С++ можно изучить на специализации C++ Developer.
Собственные значения и собственные векторы
Вычисление собственных значений — ещё одна фича Eigen. В зависимости от симметрии матрицы используйте разные методы.
Для симметричных матриц
Если матрица симметричная, используйте SelfAdjointEigenSolver — это быстро и стабильно:
#include
#include
int main() {
Eigen::Matrix3d M;
M << 4, 1, 2,
1, 3, 0,
2, 0, 5;
Eigen::SelfAdjointEigenSolver 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
#include
int main() {
Eigen::Matrix3d M;
M << 1, 2, 3,
0, 4, 5,
0, 0, 6;
Eigen::EigenSolver 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
#include
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 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
#include
int main() {
// Сырой массив данных (например, считанный из файла)
double rawData[] = {1, 2, 3, 4, 5, 6};
// Маппинг массива в матрицу 2x3 (без копирования данных)
Eigen::Map> mat(rawData);
std::cout << "Матрица из сырого массива:\n" << mat << std::endl;
// Можно менять данные напрямую
mat(0, 0) = 42;
std::cout << "\nПосле изменения первого элемента:\n" << mat << std::endl;
return 0;
} Работа с разреженными матрицами
Eigen умеет работать и с разреженными матрицами. Мастхев, когда размерность системы огромна, а ненулевых элементов мало.
#include
#include
int main() {
// Определяем разреженную матрицу размером 5x5
Eigen::SparseMatrix 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::InnerIterator it(spMat, k); it; ++it) {
std::cout << "(" << it.row() << "," << it.col() << ") = " << it.value() << "\n";
}
}
return 0;
} Работая с разреженными матрицами, всегда помните о том, что правильный выбор структуры данных может сократить время вычислений.
Оптимизация кода с Eigen
Без оптимизации даже самый красивый код может тупить. Пару советов:
Фиксированные размеры, когда возможно.
Если матрица имеет размер, известный на этапе компиляции, объявляйте её какEigen::Matrix, а неEigen::MatrixXd. Компилятор тут может провести множество оптимизаций.Используйте
.eval(), чтобы избежать лишних копирований.
В сложных цепочках вычислений результат может не материализоваться до тех пор, пока вы явно не вызовете.eval().Соблюдайте выравнивание памяти.
Eigen использует векторизацию, так что корректное выравнивание данных может дать солидный прирост производительности. Обычно это делается автоматически, но если вы работаете с внешними данными черезEigen::Map, следите за этим.
Пример, объединяющий несколько советов:
#include
#include
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 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
#include
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
#include
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
#include
#include
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
#include
#include
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(image.at(i, j));
// Применим простейшую нормализацию (можно заменить на более продвинутую операцию)
eigenImage = eigenImage.array() / 255.0;
std::cout << "Нормализованное изображение (первые 5 строк):\n" << eigenImage.topRows(5) << std::endl;
return 0;
} Это лишь малая часть возможностей.
Eigen — это не просто набор функций для работы с матрицами, это целый экосистемный подход к математическим вычислениям в C++. А как вы применяете эту библиотеку в своих проектах?
Пользуясь случаем, рекомендую обратить внимание на открытые уроки по C++:
5 марта в 20:00. Готовим рабочее место: C++ + VSCode. Подробнее
20 марта в 20:00. Пакетные менеджеры для C++ проектов. Подробнее
Habrahabr.ru прочитано 19041 раз
