Свободная библиотека линейной алгебры Eigen

Библиотека Eigen настолько хорошо документирована, что любое описание работы с ней будет неизбежно дублировать документацию. Я попытаюсь дать практическое представление о наиболее часто используемых возможностях библиотеки и «подводных камнях» при ее использовании.

Первая программа

Начнем с практического примера, в котором будем считывать из файла матрицу чисел и искать для нее собственные числа и собственные векторы:

#include <iostream>
 
#include <fstream>
 
#include <Eigen/Core>
 
#include <Eigen/Dense>
 
#include <boost/lexical_cast.hpp>
 
using namespace std; using namespace Eigen;
 
int main(int argc, char** argv){
 
int N = boost::lexical_cast<int>(argv[2]);
 
MatrixXf m(N,N);
 
ifstream f(argv[1]);
 
for(int i=0;i<N;++i)
 
for(int j=0;j<N;++j) f >> m(i,j);
 
cout << "Исходная матрица:" << endl << m << endl;
 
EigenSolver<MatrixXf> solver(m);
 
cout << "Собственные числа:" << solver.eigenvalues().transpose() << endl; cout << "Матрица собственных векторов (колонки):" << endl << solver.eigenvectors() << endl;
 
f.close();
 
}

Скомпилировать нашу программу можно так:

$ g++ prog1.cpp -I/<путь_к_Eigen>

Заголовок является базовым для Eigen и включает объявления классов для матриц и векторов. В нашем случае его можно было бы не включать явно, поскольку нужен еще и заголовок для операций линейной алгебры над плотными матрицами, который автоматически его подключает. Все объекты и функции Eigen расположены в одноименном пространстве имен.

Мы будем передавать в программу два параметра: имя файла и размер матрицы, которая в нем содержится. Объявление матрицы интуитивно понятно:

MatrixXf m(N,N);

Это матрица переменного размера («Х») типа float («f») размерностью NxN. Аналогично MatrixXd — матрица типа double и т.п. Если нужно создать нестандартную матрицу (скажем, с фиксированным числом строк, но переменным числом колонок), то можно использовать развернутую форму:

Matrix m; // 20 строк и изменяемое число столбцов

// хранение построчное (по умолчанию — по колонкам)

В случае использования хотя бы одной размерности Dynamic матрица создается динамически в куче, но если обе размерности фиксированы, то она создается в стеке, что выгоднее в плане производительности.

Матрицы с одной из размерностей 1 являются векторами-колонками или векторами строками. Разница между ними такая же, как в стандартной линейной алгебре. По умолчанию создаются вектора-колонки:

Vector3f v; // то же, что Matrix — вектора-колонка размера 3

Размер динамической матрицы можно изменить в любой момент:

m.resize(2*N, 3*N); // теперь размер 2Nx3N

Доступ к элементам матрицы производится с помощью оператора «()»:
m(1,1) = m(2,3)+m(1,3);

Можно работать с целыми колонками и строками:

VectorXf c = m.col(0); // Первая колонка

Матрицы и векторы поддерживают вывод с помощью оператора «<<», так что аккуратно напечатать матрицу можно одной строкой: cout << m << endl; После считывания матрицы из файла мы создаем объект EigenSolver для вычисления собственных чисел и векторов. В качестве шаблонного параметра он принимает тип нашей матрицы, а в конструкторе передается сама матрица. Сам конструктор и проводит нужные вычисления, так что можно сразу работать с полученными собственными значениями через метод eigenvalues() и векторами через eigenvectors(). Собственные числа выводятся в виде вектора-колонки, так что при выводе на экран мы превращаем его в вектор-строку методом transpose(). EigenSolver работает с произвольными матрицами и может выдавать комплексные результаты (типа complex или complex). Часто приходится иметь дело с симметричными матрицами, собственные числа которых гарантированно действительны. Для них вместо EigenSolver можно использовать SelfAdjointEigenSolver.

Приятные мелочи

Все библиотеки, интенсивно использующие шаблоны, как правило, страдают от очень большого времени компиляции. В случае Eigen разработчики сделали все возможное, чтобы его уменьшить, и добились в этом впечатляющих результатов — даже в сложных проектах Eigen увеличивает время компиляции незначительно. Кроме того, Eigen — чуть ли не единственная шаблонная библиотека, которая выводит вменяемые сообщения об ошибках. Например, простой, но неправильный код:

Vector3f a;

Vector3d b;

cout << a+b << endl; выдает такую ошибку:

‘YOU_MIXED_DIFFERENT_NUMERIC_TYPES_YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_

CAST_NUMERIC_TYPES_EXPLICITLY’

in ‘Eigen::internal::static_assertion

EIGEN_CHECK_BINARY_COMPATIBILIY(BinaryOp,typename Lhs::Scalar,typename Rhs::Scalar);

В отличие от обычной для шаблонных функций неудобочитаемой «тарабарщины», в этом сообщении четко указан источник ошибки — несовпадение числовых типов операндов (векторы типа float и double).

Еще одна приятная мелочь — инициализация «новорожденных» матриц и векторов нулями. Это позволяет избавиться от целого вороха неприятных ошибок и не заботиться специально об инициализации.

Линейная алгебра

В Eigen все операции линейной алгебры выполняются предельно интуитивным образом. Например, умножение матриц:

MatrixXf a,b,c; c = a*b;
Производится по обычным математическим правилам умножения матриц. Сумма или разность матриц или векторов — это обычная поэлементная сумма. Умножение или деление матрицы или вектора на число — это поэлементное скалярное умножение.

Умножение матрицы на вектор, как и в обычной математике, дает разные результаты в зависимости от типа вектора (колонка или строка) и от того, справа или слева он расположен в выражении:

 
 
 
Matrix2d mat;
 
mat << 1, 2, 3, 4; // простая инициализация
 
Vector2d u(-1,1), v(2,0);
 
std::cout << "Произведение матриц mat*mat:\n" << mat*mat << std::endl;
 
std::cout << "Умножение матрицы на вектор-столбец mat*u:\n" << mat*u << std::endl;
 
std::cout << "Умножение вектора-строки на матрицу uAT*mat:\n" J << u.transpose()*mat << std::endl;
 
std::cout << "Умножение матрицы на вектор-строку u*vAT:\n" J << u*v.transpose() << std::endl;

Скалярные и векторные произведения векторов реализованы отдельными методами dot() и cross():

Vector3d 3) 2 >
Vector3d w(0,1,2);
cout << ” 'Скалярное произведение: ” << v.dot(w) << endl; cout << ” 'Векторное произведение: ” << v.cross(w) << endl;

Поэлементные операции
Очень часто возникает необходимость обрабатывать матрицы и векторы не как математические объекты, а как обычные массивы чисел, т.е. поэлементно. Для этого в Eigen предусмотрен тип Array. По сути он ничем не отличается от типа Matrix, кроме специфики поэлементного доступа. Умножение двух массивов уже не является векторным произведением, а проводится строго поэлементно:

ArrayXXf a(2,2);
ArrayXXf b(2,2);
3, 2 < < a b << 5,6,7,8; cout << ”a * b = ’ ' << endl << a * b << endl; // выводит: a * b = 5 12 21 32

Любую матрицу можно «на лету» преобразовать в массив методом array() без потери производительности (все происходит на этапе компиляции):

VectorXf 2 2 a(
VectorXf 2 2 b(
a << 1,2, 3,4; b << 5,6, ,7,8; cout << ” 'a * b = ’ 1 << endl << a.array() * b.array() << endl;

C массивами можно производить очень сложные операции «в стиле Fortran 90/95». Например, в следующем примере перемножаются только те элементы двух массивов, для которых элемент первого массива меньше, а все остальные элементы складываются:

ArrayXf a(6); a << 1, 2, 3, 43, 15, 7; ArrayXf b(6); b << 51, 5, 0, 4, 5, 17; cout << (a < b).select(a*b, a+b).transpose() << endl; Подводные камни
Как и в любой сложной библиотеке, в Eigen есть подводные камни. Они являются расплатой за очень высокую производительность.

Главная проблема заключается в том, что для эффективной векторизации операций все объекты Eigen должны иметь специфическое 128-битное выравнивание. Для динамических матриц (например, MatrixXf) это не проблема, поскольку они создаются в куче с помощью специального распределителя.

А вот матрицы фиксированного размера создаются в стеке, и для них нельзя контролировать выравнивание. Если такой объект имеет размер, кратный 128 битам (16 байт), то его можно успешно векторизовать. Он называется «fixed-size vectorizable Eigen type» (их полный список приведен в документации: http://eigen.tuxfamily.org/dox/TopicFixedSizeVectorizable.html). Eigen «видит» объект, кратный 128 битам, но проверить его реальное выравнивание во время выполнения никак не может, поэтому генерирует код «в надежде на лучшее». Такой код будет «падать», если выравнивание в реальности окажется не 128-битным.

Такое может случиться в нескольких случаях — например, когда объекты Eigen являются членами классов:

class Foo{
 
Eigen::Vector2d v;
 
};
 
//...
 
Foo *foo = new Foo;
 
// Программа падает, поскольку Foo создается с неправильным выравниванием полей

Специальный макрос позволяет исправить ситуацию:

class Foo{
 
Eigen::Vector2d v;
 
public:
 
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
 
};
 
Foo *foo = new Foo; // Теперь работает!

Другой «тяжелый случай» — использование контейнеров STL с векторизуемыми объектами фиксированного размера — например, std::map. Чтобы обеспечить правильное выравнивание, нужно вручную указывать специальный распределитель Eigen::aligned_ allocator:

std::map,

Eigen::aligned_allocator > > cont;
В случае std::vector, из-за специфики реализации не помогает даже это, и нужно дополнительно использовать особую специализацию контейнера vector:

#include

std::vector > vec;

Другое неудобство связано с передачей объектов Eigen в функции. Рассмотрим такую простую функцию:

void add(VectorXf& a, VectorXf& b, VectorXf& res){ res = a+b;

}

На первый взгляд, все в порядке:
 
VectorXf a(3); a << 1, 2, 3;
 
VectorXf b(3); b << 4, 5, 6;
 
VectorXf r(3);
 
add(a,b,r);
 
cout << r << endl;
 
// Печатает 5 7 9

Но теперь попробуем немного видоизменить вызов нашей функции:
add(a,b+a,r);

При компиляции мы совершенно неожиданно увидим ошибку:

prog1.cpp:42:13: error: invalid initialization of non-const reference of type ‘Eigen::VectorXf& {aka Eigen::Matrix&}’ from an rvalue of type ‘const Eigen::CwiseBinaryOp, const Eigen::Matrix, const Eigen::Matrix >’

Всмотревшись в этот текст, можно уловить его смысл.
Выражение «a+b» — это не VectorXf, а некий внутренний шаблон выражения, имеющий тип Eigen::internal::scalar_sum_op. Естественно, что наша функция о нем не знает… Не вдаваясь в детали (их можно почерпнуть из документации: http://eigen.tuxfamily.org/dox/ TopicFunctionTakingEigenTypes.html), приведу пример правильно написанной функции add:

template void add( const MatrixBase& a, const MatrixBase& b, const MatrixBase& res){ const_cast< MatrixBase& >(res) = a + b;

}

Во-первых, функция стала шаблонной и каждый ее аргумент принимает значения типа MatrixBase<...>. Шаблонные аргументы для a,b и res разные. Это позволяет передавать в нашу функцию любые объекты Eigen, производные от базового типа MatrixBase — лишь бы их размерности совпадали. В частности, рассмотренное выше выражение типа Eigen::internal::scalar_ sum_op отвечает сигнатуре MatrixBase<...> и спокойно проходит компиляцию.

Во-вторых, «выходной» параметр res объявлен как const, а потом в теле функции с него вручную снимается константность с помощью const_cast. Это «шаманство» необходимо для того, чтобы в выходном параметре можно было также использовать произвольные выражения (например, отдельную колонку матрицы). Такой подход гарантированно избавляет от ненужных временных переменных при передаче аргументов и дает максимальную производительность кода. Расплатой является не самая красивая сигнатура функции.

Выводы

На сегодняшний день Eigen по праву считается лучшей свободной библиотекой линейной алгебры для C++, которая сочетает прекрасную производительность с предельно простым и интуитивно понятным синтаксисом. Знакомство с Eigen будет полезно любому программисту, кто так или иначе сталкивается с вычислительными задачами. eof|

Семен Есилевский [email protected]


http://blog.wel.org.ua

работаю админом, прогером сеошнегом :)

Leave a Comment

Ваш e-mail не будет опубликован. Обязательные поля помечены *

Загрузка...
Menu Title