пʼятницю, 20 липня 2012 р.

Робота з матрицями в С++ -- бібліотека eigen

(c) Eigen FAQ
Багато років я шукав зручну і швидку бібліотеку лінійної алгебри для використання з C++ (*1). Вимоги ніби скромні -- щоб вміла створювати, додавати, віднімати матриці, множити їх на скаляр, множити між собою, шукати обернені та детермінанти. Все. При чому, для зовсім маленьких матриць -- розміром до 8х8, а в основному взагалі 2х2. (*2) Ну, ще Open Source ліцензія була б бажаною. :-) Але ця вимога не критична.

Відволікаючись поки від теми продуктивності, С++ має всі засоби, щоб зробити роботу з матрицями зручною. Однак більшість авторів матричних бібліотек вдавалися в пояснення, чому писати: add_matrix(A,B) багато крутіше, ніж A+B. Завершилося все тим, що я, на базі простенької бібліотеки Slate, склепав свою матричну бібліотеку -- Slate2, зручну для моїх потреб, хоча й доволі повільну (*3). Подивитися на останню версію цього чуда можна тут (ніяким іншим чином не розповсюджував). Щастя наступило лише кілька років тому -- натрапив на eigen!

Виникла eigen всередині проекту KDE. Зараз вона надзвичайно зручна, швидка, та вміє доволі багато. Зокрема, використовує шаблонне метапрограмування, щоб обійти деякі проблеми матричних бібліотек в C++ (*4). Вміє самостійно векторизувати, використовуючи SIMD-інструкції, такі як SSE, AltiVec і ARM NEON. Eigen3 -- третя версія, векторизує навіть роботу з матрицями комплексних чисел. Працює із фіксованими, динамічними та розрідженими матрицями. В ролі елементів матриць можуть бути числа з плаваючою крапкою, комплексні числа, цілі числа, користувацькі типи. Що вміє робити над матрицями? Багато, перераховувати довго, див. краще сюди: 1, 2, 3. Надійна -- містить більше півтисячі власних тестів, проходить стандартні тести BLAS та LAPACK , використовується багатьма проектами.

З того часу, як я натрапив на eigen ситуація покращилася -- з'явилися Armadillo, MTL4 (її попередня версія, MTL2 -- якраз наглядний приклад згадуваного вище багатослівного підходу), тощо. Як на мене, вони безперечно вартують уваги, але поступаються eigen, хоча фахово не скажу -- багато тестів не писав, багато часу на ознайомлення з ними не витрачав. Щоб відобразити своє, хай не до кінця науково перевірене, ставлення до eigen, зацитую її розробників:

Some other libraries do satisfy very well certain specialized needs, but none is as versatile as Eigen, has such a nice API, etc.

*** skipped ***

The state of existing matrix libraries before Eigen is that:
 some are Free Software
 some are fast
 some have a decent API
 some handle fixed-size matrices, some handle dynamic-size dense matrices, some handle sparse matrices
 some provide linear algebra algorithms (LU, QR, ...)
 some provide a geometry framework (quaternions, rotations...)

However Eigen is the first library to satisfy all these criteria.

Ну і ліцензія хороша -- LGPL3+.


Як почати використовувати eigen


Бібліотека цілком складається з include-файлів, для свого використання не вимагає компіляції. Достатньо розархівувати її кудись та вказати це місце компілятору. Як її можна отримати, звідки завантажити, описано тут. Ми вважатимемо, що використовується eigen3, орієнтуючись, коли це гратиме роль, на версію 3.1.

Базується бібліотека на С++ стандарту 1998-го року, тому мала б працювати з будь яким сумісним компілятором. Перевірено з наступними компіляторами:
  • gcc -- новіші за версію 3.4, зокрема -- MinGW під Win32,
  • MSVC (Microsoft Visual Studio) -- 2005 та новіші,
  • Intel C++ Compiler (icc),
  • LLVM/CLang++ (2.8 і новіші),
  • QNX's QCC compiler.
Тести коректності та продуктивності (benchmarks) компілюються за допомогою CMake, але про них іншим разом.

Надалі вважатиметься, що використовується компілятор gcc (зокрема, це може бути MinGW).

Найпростіша програмка виглядатиме якось так:

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

using Eigen::MatrixXd;

int main()
{
  MatrixXd m(2,2);
  m << 3, 2.5, -1, 0;
  m(1,1) = m(1,0) + m(0,1);
  std::cout << m << std::endl;
}

Увага: приклади часто базуються на прикладах з офіційної документації.

Скомпілювати її можна так:

g++.exe -O3 -I<PATH-TO-EIGEN-HEADERS> program.cpp -o program

Де опція -O3 задає максимальну оптимізацію, -I -- шлях до include-файлів бібліотеки, а -o -- назву вихідного файлу (під Windows MinGW сам додасть розширення exe).

Результат виконання буде таким:

  3 2.5
 -1 1.5

Header-файл, який ми включили -- Eigen/Dense, підключає базові типи векторів, матриць, багато операцій лінійної алгебри (обчислення детермінантів, LU-декомпозиція, розклад Холецького, сингулярний розклад матриці (SVD drcomposition), QR розклад,  обчислення власних значень -- див. "Modules and Header files").

Всі типи оголошено в просторі імен Eigen. MatrixXd -- матриця довільного розміру (X) з елементами типу double (d). Конструктору MatrixXd передається розмір матриці, 2x2 в цьому прикладі. Потім за допомогою оператора << та оператора коми ініціалізуємо її -- зліва направо, по рядках, зверху вниз.

Показано, як звертатися до окремих елементів матриці за їх індексами:

m(1,1) = m(1,0) + m(0,1);

Приклад матричних обчислень з використанням eigen


Подивимося, як із допомогою цієї бібліотеки можна виконувати різноманітні матричні операції. Щоб кожен раз не копіювати, вважатимемо що всі програми починаються таким заголовком:

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

using namespace Eigen;
using namespace std;

Додавання та віднімання матриць:


    // Матриці фіксованого розміру 2х2
    Matrix2d A,B,C;
    A << 1, 2,
         3, 4;

    B << 4, 3,
         2, 1;

    //Виведемо їх
    cout << "A=\n" << A << endl;
    cout << "B=\n" << B << endl;
    cout << endl;
    //Матриці змінного розміру, встановлено ті ж 2х2
    MatrixXd Ar(2,2), Br(2,2), Cr(2,2);

    Ar << 10, 20,
          30, 40;

    Br << 40, 30,
          20, 10;

    //Виведемо їх
    cout << "A-resizable=\n" << Ar << endl;
    cout << "B-resizable=\n" << Br << endl;
    cout << endl;

    // Додамо та віднімемо матриці фіксованого розміру
    C = A+B;
    cout << "A+B=\n" << C << endl;
    C = A-B;
    cout << "A-B=\n" << C << endl;

    // Додамо та віднімемо матриці змінного розміру
    Cr = Ar+Br;
    cout << "A-resizable+B-resizable=\n" << Cr << endl;
    Cr = Ar-Br;
    cout << "A-resizable-B-resizable=\n" << Cr << endl;

    cout << endl;

    // Додамо сумісні за розміром матриці фіксованого та змінного розміру
    C = Ar+A;
    cout << "A-resizable+A-fixed=\n" << C << endl;

Результат роботи цього фрагменту буде таким:

A=
1 2
3 4
B=
4 3
2 1

A-resizable=

10 20
30 40
B-resizable=
40 30
20 10

A+B=

5 5
5 5
A-B=
-3 -1
1  3
A-resizable+B-resizable=
50 50
50 50
A-resizable-B-resizable=
-30 -10
10  30

A-resizable+A-fixed=

11 22
33 44

Як бачимо, можна легко змішувати матриці різного типу -- змінного і фіксованого розміру, у тих же виразах. Eigen намагатиметься перевірити коректність такого змішування. Якщо це можливо - на стадії компіляції, і під час виконання у відладочному режимі, якщо ні.

Додавання та віднімання векторів:


Все вищесказане стосується і векторів:

    Vector2d V, U, R;
    V << 1, 5;
    U << 5, 1;
    cout << "V = \n" << V << endl;
    cout << "U = \n" << U << endl  << endl;

    R = V+U;
    cout << "V+U = \n" << R << endl << endl;

    VectorXd Vr(2), Ur(2), Rr(2);
    Vr << 10, 50;
    Ur << 50, 10;
    cout << "V-resizable = \n" << Vr << endl;
    cout << "U-resizable = \n" << Ur << endl << endl;

    Rr = Ur-U+Vr;
    cout << "U-resizable-U+V-resizable = \n" << Rr << endl;

Результатом роботи буде:

V =
1
5
U =
5
1

V+U =
6
6

V-resizable =
10
50
U-resizable =
50
10

U-resizable-U+V-resizable =
55
59


Множення:


Матриці та вектори можна множити як між собою, так і на скаляри. Нотація залишається природною. На скаляр можна також і ділити :-)

    C = A*B; // Множимо матрицю на матрицю

    Cr = A*Br - 2*Ar + B/5; // Множимо матриці різних типів, 
                            // множимо та ділимо матриці на скаляр

    R = A*V+5*U; //Множимо матрицю на вектор, вектор на скаляр

(Матриці та вектори використано ті ж, що й у прикладах вище).

Eigen коректно обробляє також потенційно небезпечні (для бібліотеки, що використовує lazy evaluation) операції, типу A=A*A -- за потреби тимчасові змінні створюються автоматично. Детальніше див. тут.

Також реалізовано скалярний і векторний (для трьохмірних векторів) добуток:

    Vector3d V3(1,2,3);
    Vector3d W3(0,1,2);
    Vector3d R3;

    double dp = V3.dot(W3); // скалярний добуток
    R3 = V3.cross(W3); // векторний добуток

    dp = V.dot(Ur); // скалярний добуток 2-мірних векторів


Редукція:


Бібліотека вміє знаходити певні "редуковані" характеристики матриць та векторів -- суму ( A.sum() ) , добуток ( A.prod() ) всіх елементів, максимальний і мінімальний елемент (також і його координати в матриці) -- A.minCoeff(), A.maxCoeff() ( A.minCoeff(&i,&j), A.maxCoeff(&i,&j), V.maxCoeff(&i), V.minCoeff(&i) ), слід (шпур) матриці - A.trace(), середнє компонент - A.mean(). Можна скористатися частковою редукцією та іншими "суміжними" засобами -- див. "Tutorial page 7 - Reductions, visitors and broadcasting". 

Масиви:


На додачу до спеціалізованих матричних/векторних типів, eigen містить і загальний тип -- масив, Array, із поелементною обробкою.

    ArrayXXd  M1(2,2), M2(2,2);

    M1 << 1,2,3,4;
    M2 << 4,3,2,1;

    cout << "M1*M2=\n" << M1*M2 << endl; // Добуток - поелементний! 
    cout << "M1 + =\n" << M1 + 4 << endl; // Добуток - поелементний! 

Результатом виконання цього коду буде щось таке:

M1*M2=
4 6
6 4 
M1 + 4=
5 6
7 8 

Крім арифметичних операцій, та вже згаданих для матриць/векторів функцій редукції, існують поелементні функції, такі як обчислення квадратного кореня, абсолютне значення числа.

Змішані операції над матрицями та масивами вимагають мінімальних зусиль.

Перемножити матриці поелементно можна, наприклад, так :

    C =  A.array() * B.array(); //Поелементне множення матриць
    C =  A.array() * Br.array(); // Поелементне множення матриць різних типів

    C = A.cwiseProduct(B); // Альтернативни спосіб
    C = A.cwiseProduct(Br);

Тобто, можна скористатися як API масивів так і матричним API. І навпаки, матрично перемножити два масиви чи масив і матрицю, можна так:

    C = M1.matrix()*M2.matrix(); // Матричне множення масивів M1 і M2
    C = A*M2.matrix(); // Матричне множення матриці A і масиву M1

Лінійна алгебра, розкладання матриць (LU, QR, SVD), власні значення, тощо:


Без всього, перерахованого в назві цього пункту, цінність бібліотеки була б значно меншою. Почнемо з найпростішої задачі -- розв'язати систему лінійних рівнянь: \[\hat{A} \vec{x}=\vec{b}\]. Розв'яжемо її за допомогою LU-розкладу:

   Matrix3d A3;
   Vector3d b3;
   A3 << 1,2,3,  4,5,6,  7,8,10;
   b3 << 3, 3, 4;
   cout << "A3=\n" << A3 << endl;
   cout << "b=\n" << b3 << endl;
   Vector3d x = A3.lu().solve(b3);
   cout << "x=\n" << x << endl;

Результат буде наступним:

A3=
 1  2  3
 4  5  6
 7  8 10
b=
3
3
4
x=
-2
1
1

Легко можна перевірити його правильність. :-)

Якщо постає задача розв'язати багато таких систем, із тією ж матрицею A, але різними b, результат розкладу матриці можна зберегти та використовувати скільки завгодно. Зробимо це на прикладі іншого доступного в eigen методу розкладу -- QR. (Зберегти результат LU-розкладу так само просто). Використаємо ту ж A3, b3, що й в попередньому прикладі, а також ще одну праву частину:

    Vector3d b3_2;
    b3_2 << 3, 3, 3;

    ColPivHouseholderQR<Matrix3d> decomposition(A3);

    x = decomposition.solve(b3);
    cout << "x=\n" << x << endl << endl;
    x = decomposition.solve(b3_2);
    cout << "x=\n" << x << endl << endl;

Результат першого обчислення, як і можна сподіватися, співпадає із обчисленим методом LU-розкладу.

x=
-2
1
1

x=
-3
3
2.95037e-016

Список розкладів, які можна використовувати -- "Catalogue of decompositions offered by Eigen".

Вміє бібліотека знаходити також власні вектори і власні значення:

     SelfAdjointEigenSolver<Matrix3d> eigensolver(A3);
     if (eigensolver.info() != Success)
        cout << "Error!" << endl;
     else
     {
        cout << "The eigenvalues of A3 are:\n" << eigensolver.eigenvalues() << endl;
        cout << "Here's a matrix whose columns are eigenvectors of A3 \n"
             << "corresponding to these eigenvalues:\n"
             << eigensolver.eigenvectors() << endl;
     }

Даючи наступний результат:

The eigenvalues of A3 are:
-2.85563
-0.524024
19.3796
Here's a matrix whose columns are eigenvectors of A3
corresponding to these eigenvalues:
 0.844671  0.355502  0.400186
 0.104531 -0.842767  0.528032
 -0.52498  0.404182  0.749022

Не є проблемою обчислення детермінанту та оберненої матриці:

    cout << "The determinant of A3 is " << A3.determinant() << endl;
    cout << "The inverse of A3 is:\n" << A3.inverse() << endl;

The determinant of A3 is -2
The inverse of A3 is:
  -2    1
 1.5 -0.5

Є засоби для вирішення задачі методу найменших квадратів, обчислення рангу матриць і т.д. -- див. "Tutorial page 6 - Linear algebra and decompositions".


Решта:


Вміє eigen працювати і з блоками, вирізаними з матриць. Деталі поки див. тут.
Володіє розвинутою системою ініціалізації матриць -- "Tutorial page 5 - Advanced initialization".
Працює з розрідженими матрицями. Включає невеликий, але зручний, геометричиний модуль.
Цікавими є "unsupported" -- написані користувачами, модулі: "Unsupported modules". Зокрема, поміж них - FFT (швидке перетворення Фур'є), сплайни, підтримка чисел із плаваючою крапкою довільної точності (MPFRC++ Support module), нелінійна оптимізація, чисельне диференціювання, поліноми, взаємодія з OpenGL, тощо.
Багато важливих, але не згаданих тут питань, обговорюється в документації: взаємодія з STL, використання MKL, багатопоточність і т.д. і т.п.

Післямова


Щоб цей пост теж не роздувся до нереальних розмірів, на цьому і завершу. Нехай він буде більше рекламним за своєю природою. ;-)

Звичайно, скоро заплановано наступні пости, присвячені eigen! На загал, орієнтуватимуся в них на офіційний підручник, однак не ставлячи собі цілі відтворити його повністю.


Виноски


(*1) Тому BLAS і LAPACK та їх прив'язки (bindings) до C/C++ використовувати можна, але краще не треба. Не зважаючи на регулярні спроби мене переконати, ніяк не можу повірити, що щось типу (чесно - писав згідно мануалів, компілятором не перевіряв, вірю, що можна краще):

    double A[9];
    double b[3];
    double b_copy[3];
    double relative_error;
    double dummy;

// Ініціалізація A i b, та чи інша

    int N = 3;
    int nrhs = 1;
    int lda = 3;
    int ipiv[3];
    int ldb = 3;
    int info;

    dgesv_(&N, &nrhs, A, &lda, ipiv, b, &ldb, &info);
    for(int i=0;i<N;++i)
         b_copy[i]=b[i];

    dgemm_('N', 'N', &N, &nrhs, &N, 1.0, A, &lda, x, &ldb, -1.0, b, &ldb);
    //now b == A*x-b -- spoiled
   relative_error=dlange_('F', &N, &nrhs, b,  &lda, &dummy)/dlange_('F' , &N, &nrhs, b_copy, &lda, &dummy);

є більш зручним, ніж:

    Matrix<double, 3, 3> A;
    Matrix<double, 3, 1> b,x;

// Ініціалізація A i b, та чи інша

    x = A.fullPivLu().solve(b)
    double relative_error = (A*x - b).norm() / b.norm();

Особливо, якщо різниця в швидкодії - лічені проценти, та ще й не завжди на користь "класики". Див., наприклад, сюди.

(*2) Так, матиці 2х2 можна взагалі вручну. Але, навіть не говорячи, що це відверте велосипедобудування, яке забирає сили і час, воно ще й збільшує ймовірність помилок.

(*3) Забігаючи наперед, eigen2 вона поступалася процентів на 30, і це для комплексних чисел, які остання не вміла векторизувати.

(*4) Приклад таких проблем: якщо матриці оголошені об'єктами, при операціях із ними регулярно виникатиме потреба робити їх копії. Скажімо, єдиний коректний спосіб повернення результату додавання матриць наступний:

SomeMatrixType operator+(const SomeMatrixType& A, const SomeMatrixType& B)
{
    SomeMatrixType result;
    result=... Певний алгоритм додавання ...
    return result;
}


Повертати посилання чи ще якось хитрувати, як показали дослідження класиків C++, Меєрса, Саттера, Александреску та інших -- небезпечно і призводить до помилок.

Компілятор намагається їх "відоптимізувати", але не всі такі марні копіювання може виявити.

Ще одна потенційна проблема з продуктивністю може виникати через так-званий aliasing. Грубо кажучи, виконуючи операції над послідовністю елементів, наприклад додаючи два масиви і зберігаючи результат в третьому, компілятор далеко не завжди може визначити, чи ці масиви не перекриваються між собою та з ще якимись даними, що призводить до необхідності повторного завантаження даних із пам'яті -- для забезпечення коректності. Див. також "Understanding Strict Aliasing".

Детальніше про те, як eigen вирішує ці проблеми можна почитати, на простому прикладі, тут.

Немає коментарів:

Дописати коментар