Pet-project: мини-библиотека по линейной алгебре

Введение

Однажды меня попросили рассказать о своем опыте разработки математических алгоритмов. Так как коммерческий опыт у меня был преимущественно в веб-разработке, то рассказать я мог только об университетском опыте, либо реализовать собственный pet-проект.Я выбрал тему линейной алгебры.

Существовало два варианта реализации проекта: с интерфейсом на Qt либо в виде решения, которое можно использовать в backend-разработке. Я выбрал второй вариант и реализовал небольшую библиотеку.

В ходе разработки мне пришлось ответить на следующие вопросы:

- Как в принципе реализуются библиотеки
  и как сделать её подключение и использование удобной для пользователя?

- Сколько прироста производительности даёт
  переход от 2D векторов к плоским матрицам?

- Какие возможности оптимизации есть?
  Какие подходят для данного проекта?

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

Проект не является промышленной библиотекой и не ставит своей целью создание аналога Eigen или Armadillo.

Архитектура

На данный момент проект состоит из двух основных компонентов:

-   Хранение (MatrixBase.hpp,FlatMatrix.hpp)

-   Декомпозиция(LU.hpp)

    - разложение PA = LU

    - решение систем Ax = b

    - вычисление определителя

    - обращение матрицы

LU-декомпозиция.

Реализована LU-декомпозиция с частичным выбором главного элемента (partial pivoting). В ходе разложения вычисляются верхне-треугольная матрица U и нижне-треугольная матрица L. Они хранятся совместно в одной матрице (in-place). Проверка корректности происходит в тестах, где проверяется равенство PA=LU. В классе есть ограничение это работа с типами с плавающей запятой, поскольку при вычислении элементов матриц U и L возникают операции деления, как правило в следствии деления будут дробные значения. Использование целочисленных типов приводит к потере дробной части из-за целочисленного деления и, как следствие, к искажению промежуточных результатов и ошибкам вычислений.

Если кратко:

1. Находим опорный элемент. Самый большой элемент в текущем столбце k. Метод:pivoting

2. Меняем опорную строку с текущей строкой

3. Далее заполняем треугольники L(нижняя часть матрицы) и U(верхняя часть матрицы), метод: elimination

4. Там же вычисляется знак определителя signP Определитель вычисляется как произведение диагональных элементов матрицы U с учётом знака перестановок где signP принимает значения 1 или -1 в зависимости от чётности числа перестановок строк.

  template<typename T, typename MatrixType>
    void LU<T, MatrixType>::decomposition() {
        if (matrix.getRows() != matrix.getCols())
        throw std::runtime_error("Matrix must be square for LU decomposition");
        int n = matrix.getRows();
        initP();
        for (int k = 0; k < n; ++k) {
            int pivot = pivoting(k);
            if(pivot != k) {
                for (int j = 0; j < n; j++) {
                    std::swap(matrix(k,j), matrix(pivot,j));
                }
                
                std::swap(P[k], P[pivot]);
                signP = -signP;
            }
            if (std::abs(matrix(k,k)) <= eps) {
                throw std::runtime_error("Matrix is singular or nearly singular at pivot " + std::to_string(k));
            }
            elimination(k);
        }
    }

Проверка корректности разложения

Для проверки корректности разложения нам необходимо проверить равенство. Изначально формируется матрица PA. Каждая строка результирующей матрицы берётся из соответствующей строки исходной матрицы согласно вектору перестановок:

Таким образом, матрица PA эквивалентна произведению перестановочной матрицы P на исходную матрицу A. Далее мы извлекаем матрицы L и U из основной матрицы. И произведение LU сравнивается с PA.

TEST_F(LU_Inplace_Test,decomposition) {
    FlatMatrix<double>exp = {
        {6, 3, 4},
        {5.0/6, -4.5, -19.0/3},
        {1.0/3, -8.0/9, 1.0/27}
    };

    std::vector<int> P = lu->getP();

    int n = A.getRows();
    FlatMatrix<double> PA = A;
    for (int i = 0; i < n; i++) {
        for(int j = 0; j < n; j++) {
            PA(i,j) = A(P[i],j);
        }
    }
    auto m = std::move(lu->getMatrix());
    auto L = extractL<double>(m);
    auto U = extractU<double>(m);
    auto LUproduct = L*U;
    compare<double, FlatMatrix<double>>(PA,LUproduct);
}

Обращение матриц

Обращение матрицы можно свести к решению линейной системы AX=I. После декомпозиции LU вида PA=LU вычисляется обратная матрица. Каждый столбец обратной матрицы вычисляется методами прямой() и обратной() подстановки. Деление на нет, так как диагональ в треугольнике равна единице.

Прямая подстановка

L в данной реализации это нижний треугольник матрицы. В данной реализации это вектор единичной матрицы где единица находится в позиции поэтому вместо того чтобы выделять память под новый вектор b на каждом шаге, мы передаём в метод позицию с единицей в векторе b.

Метод прямой подстановки для обратной матрицы:

    void forwardSubstitution(std::vector<T>& y, int b_pos) const;

Метод прямой подстановки для решения системы уравнений:

    void forwardSubstitution(std::vector<T>& y, const std::vector<T>& b) const;

Обратная подстановка:

U в данной реализации это верхний треугольник матрицы. Метод обратной подстановки:

void backwardSubstitution(std::vector<T>& x, const std::vector<T>& y, int n) const;

Метод обращения матрицы на основе LU декомпозиции:

   template<typename T, typename MatrixType>
    MatrixType LU<T, MatrixType>::inv() const {
        int n = matrix.getRows();
        MatrixType X(n,n);
        std::vector<int> invP = initInvP();
        std::vector<T>  y(n), x(n);
        for(int i = 0; i < n; ++i) {
            int b_pos = invP[i];

            forwardSubstitution(y, b_pos, n);
            backwardSubstitution(x, y, n);

            for (int k = 0; k < n; ++k) {
                X(k,i) = x[k];
            }
        }
        return X;
    }

Оптимизация LU

Изначально была реализация была на 2D векторах с раздельными матрицами L U. Цель выбора 2D матрицы и раздельных матриц L и U была в попытке сделать решение наглядным с явными шагами в решении. Но в после замеров было решено оптимизировать код. В ходе оптимизации был переход к хранению в плоской матрице и объединению L и U в одной матрице.

Бенчмарки были проведены в следующих условиях:

Флаги: -O2 -march=native Процессор: 12th Gen Intel® Core™ i5-12450H (2.00 GHz)

Расчёт обратной матрицы

Выигрыш в производительности небольшой. В дальнейшем планируется применить блочную оптимизацию LU.

Применение блочной оптимизации в умножении матриц

Для наглядности эффективности данного подхода я решил продемонстрировать на самой простой операции это матричное умножение.

Как видно эффективность данного подхода значительная.

Пример использования

Чтобы показать практическое применение библиотеки, был рассмотрен пример решения задачи линейной регрессии.

Подключение библиотеки в проекте

cmake_minimum_required(VERSION 3.14)
project(RandomMathApp LANGUAGES CXX)

find_package(LinearAlgebra REQUIRED)

add_executable(random_math_app main.cpp)

target_link_libraries(random_math_app PRIVATE LinearAlgebra::LinearAlgebra)
target_compile_features(random_math_app PRIVATE cxx_std_17)

Рассмотрим задачу нахождения коэффициентов линейной регрессии методом наименьших квадратов. Дана матрица признаков X и вектор наблюдений. Коэффициенты вычисляются по формуле: В данном примере производится обращение матрицы XTX с декомпозицией LU. Пример кода:

#include <iostream>
#include <decompose/LU/LU.hpp>
#include <matrix/FlatMatrix.hpp>
#include <decompose/LU/Inversion.hpp>
#include <memory>

int main() {
    using namespace LinearAlgebra;
    
    FlatMatrix<double> X = {
        {1,1},
        {1,2},
        {1,3},
        {1,4}
    };

    ColumnVector<double> y = {2,4,5,7};

    auto XT = ~X;
    auto XTX = XT*X;
    auto XTy = XT*y;
    auto lu = LU<double,FlatMatrix<double>>(std::move(XTX));

    auto Xinv = lu.inv();
    
    auto B = Xinv * XTy;
   
    std::cout << "Input data:n";
    std::cout << "X: " << X << "n";
    std::cout << "y:" << std::endl;
    for (int i = 0; i < y.size(); i++)
    {
        std::cout<<y[i]<<"n";
    }
    std::cout << std::endl;

    std::cout << "Computing linear regression coefficients "
            << "B0 and B1 using the least squares method B = (X^T X)^-1 X^T y:n";

    std::cout << "B0 = " << B[0] << "n";
    std::cout << "B1 = " << B[1] << std::endl;
}

Что можно улучшить

- Добавить методы декомпозиций(QR,SVD)
- Провести оптимизацию
- Сделать более масштабируемую архитектуру

Вывод

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

Проведённые бенчмарки показали, что переход на плоское хранение данных и объединение матриц L и U в одной структуре даёт лишь умеренный прирост производительности. В то же время измерения блочного умножения демонстрируют, что дальнейшие оптимизации, учитывающие особенности кэш-памяти процессора, могут дать значительно больший эффект.

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

Ссылки

Репозиторий: https://github.com/web-dev137/linear-library

Пример использования: https://github.com/web-dev137/random-math-app