Оптимизация свёртки эрмитовской матрицы с вектором

Лапша на плюсах для бенча:

#include <Eigen/Core>
#include <random>
#include <vector>
#include <complex>
#include <numeric>
#include <iostream>
#include <iomanip>

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/hermitian.hpp>
#include <boost/numeric/ublas/io.hpp>

#include <chrono>
using std::chrono::high_resolution_clock;
using std::chrono::duration;
using std::chrono::milliseconds;

// Number of partial waves
constexpr unsigned n_PWs = 10;
// Actual number of partial waves, always <= n_PWs
constexpr unsigned n_compressed_default = 7;
constexpr unsigned n_test_ev_default = 100000;

using data_t =
std::vector <std::array <std::array <std::complex <double>, n_PWs>, n_PWs>>;
using params_t = std::array <std::complex <double>, n_PWs>;

template <typename real_t> void report (real_t res, double fcn_lat)
{
    std::cout << "\t" << std::setw (15) << "- ln L = "
        << std::setw (15) << res << " [n. u.]" << std::endl;
    std::cout << "\t" << std::setw(15) << "Float size: " 
        << std::setw (15) << sizeof (res) << " [Bytes] " << std::endl;
    std::cout << "\t" << std::setw (15) << "FCN duration:"
        << std::setw (15) << fcn_lat << " [ms]" << std::endl;
}

    template <typename real_t>
void eigen_fixed (const data_t &f_data, const params_t &f_params,
        int n_compressed, unsigned long n_test_ev)
{
    using namespace Eigen;
    using matrix_t = Matrix<std::complex<real_t>, n_PWs, n_PWs>;
    using vector_t = Matrix<std::complex<real_t>, n_PWs, 1>;
    std::vector <matrix_t> data (f_data.size ());
    for (auto i = 0U; i < n_test_ev; i++) {
        for (auto w = 0U; w < n_compressed; w++) {
            for (auto v = 0U; v < n_compressed; v++) {
                data [i] (w, v) = f_data [i] [w][v];
            }
        }
    }
    vector_t params;
    for (auto w = 0U; w < n_compressed; w++) {
        params (w) = std::conj (f_params [w]);
    }
    // Start benchmark
    auto t1 = high_resolution_clock::now ();
    real_t res = std::accumulate (data.begin (), data.end (),
            static_cast <real_t>(0.),
            [&] (const auto &x, const auto &y) {
            return x - std::log (std::abs (std::real (
                            params.dot(y*params))));});
    auto t2 = high_resolution_clock::now ();
    auto fcn_lat = duration<real_t, std::milli>(t2-t1).count();
    std::cout << "=== Eigen with fixed size matrix" << std::endl;
    report (res, fcn_lat);
    data.clear ();
}

    template <typename real_t>
void eigen_dynamic (const data_t &f_data, const params_t &f_params,
        int n_compressed, unsigned long n_test_ev)
{
    using namespace Eigen;
    using matrix_t = Matrix<std::complex<real_t>, Dynamic, Dynamic>;
    using vector_t = Matrix<std::complex<real_t>, Dynamic, 1>;
    std::vector <matrix_t> data (f_data.size ());
    for (auto i = 0U; i < n_test_ev; i++) {
        data [i].resize (n_compressed, n_compressed);
        for (auto w = 0U; w < n_compressed; w++) {
            for (auto v = 0U; v < n_compressed; v++) {
                data [i] (w, v) = f_data [i] [w][v];
            }
        }
    }
    vector_t params (n_compressed);
    for (auto w = 0U; w < n_compressed; w++) {
        params (w) = std::conj (f_params [w]);
    }
    // Start benchmark
    auto t1 = high_resolution_clock::now ();
    real_t res = std::accumulate (data.begin (), data.end (),
            static_cast <real_t>(0.),
            [&] (const auto &x, const auto &y) {
            return x - std::log (std::abs (std::real (
                            params.dot(y*params))));});
    auto t2 = high_resolution_clock::now ();
    auto fcn_lat = duration<real_t, std::milli>(t2-t1).count();
    // Report
    std::cout << "=== Eigen with dynamic size matrix" << std::endl;
    report (res, fcn_lat);
    data.clear ();
}

    template <typename real_t>
void boost_herm (const data_t &f_data, const params_t &f_params,
        int n_compressed, unsigned long n_test_ev)
{
    namespace bnu = boost::numeric::ublas;
    using matrix_t = bnu::hermitian_matrix<std::complex<real_t>,
          bnu::upper>;
    using vector_t = bnu::vector<std::complex<real_t>>; 
    std::vector <matrix_t> data (f_data.size ());
    for (auto i = 0U; i < n_test_ev; i++) {
        matrix_t m (n_compressed);
        for (auto w = 0U; w < n_compressed; w++) {
            for (auto v = 0U; v < n_compressed; v++) {
                m (w, v) = f_data [i] [w][v];
            }
        }
        data [i] = m;
    }
    vector_t params (n_compressed);
    for (auto w = 0U; w < n_compressed; w++) {
        params (w) = f_params [w];
    }
    // Start benchmark
    auto t1 = high_resolution_clock::now ();
    real_t res = std::accumulate (data.begin (), data.end (),
            static_cast <real_t>(0.),
            [&] (const auto &x, const auto &y) {
            return x - std::log (std::abs (std::real (
                            bnu::inner_prod (params, bnu::prod (y, bnu::conj (
                                        params))))));});
    auto t2 = high_resolution_clock::now ();
    auto fcn_lat = duration<real_t, std::milli>(t2-t1).count();
    // Report
    std::cout << "=== Boost hermitian matrix" << std::endl;
    report (res, fcn_lat);
    data.clear ();
}

int parse_args (int argc, char **argv, unsigned long &n_test_ev,
        unsigned &n_compressed)
{
    if (argc == 2) {
        n_test_ev = atoi (argv[1]);
    }
    if (argc == 3) {
        n_test_ev = atoi (argv[1]);
        n_compressed = atoi (argv[2]);
        if (n_compressed > n_PWs) {
            std::cout << "Input error: second argument must be in [0, "
                << n_PWs << "] range" << std::endl;
            return 1;
        }
    }
    return 0;
}

int main (int argc, char **argv)
{
    unsigned long n_test_ev = n_test_ev_default;
    unsigned n_compressed = n_compressed_default;
    if (parse_args (argc, argv, n_test_ev, n_compressed)) {
        return 1;
    }
    std::mt19937 gen;
    std::random_device rd;
    gen.seed (rd ());
    std::uniform_real_distribution <> rand_flat_PW (-1.0, 1.0);
    std::uniform_real_distribution <> rand_flat_diag (1.0, 5.0);
    data_t data (n_test_ev);
    params_t params;
    for (auto i = 0U; i < n_test_ev; i++) {
        for (auto w = 0U; w < n_compressed; w++) {
            for (auto v = w+1; v < n_compressed; v++) {
                data [i][w][v] = std::complex <double> (
                        rand_flat_PW (gen),
                        rand_flat_PW (gen));
            }
            data [i][w][w] = std::complex <double> (
                        (rand_flat_diag (gen)), 0.);
        }
        for (auto w = 0U; w < n_compressed; w++) {
            for (auto v = 0; v < w; v++) {
                data [i][w][v] = std::conj (data [i][v][w]);
            }
        }
    }
    for (auto w = 0U; w < n_compressed; w++) {
        params [w] = std::complex <double> (
                rand_flat_PW (gen), rand_flat_PW (gen));
    }
    eigen_fixed<float> (data, params, n_compressed, n_test_ev);
    eigen_fixed<double> (data, params, n_compressed, n_test_ev);
    eigen_dynamic<float> (data, params, n_compressed, n_test_ev);
    eigen_dynamic<double> (data, params, n_compressed, n_test_ev);
    boost_herm<float> (data, params, n_compressed, n_test_ev);
    boost_herm<double> (data, params, n_compressed, n_test_ev);
    std::cout.flush ();
    return 0;
}

Клей на баше:

#!/bin/bash
FLAGS="-O3 -ffast-math $(pkg-config --cflags eigen3)"
echo
echo Compiler flags: $FLAGS
echo
g++ mll-bench.cpp $FLAGS
N_EVENTS=1000000
echo Number of events: $N_EVENTS
echo "************************* Matrix size 10x10 ************************"
./a.out $N_EVENTS 10
echo "*************************  Matrix size 5x5 ************************"
./a.out $N_EVENTS 5
rm a.out

Ничего не понимаю, но

params.dot(y*params)

эквивалентно

bnu::inner_prod (params, bnu::prod (y, bnu::conj (params)))

?

В последнем выражении берется сопряженное (conj?) от params.

Свернуть матрицу с параметрами это в данном случае умножить матрицу на вектор?:-)

Взять gpgpu потолще и не мучаться. Можно конечно че то из эрмитовости матриц наоптимизировать, но скорее всего уже все написано, тем более матрицы то небольшие…

как лучше всего сворачивать матрицы с параметрами?

Свёрточной сетью с параметрическими слоями. // К.О.
Подбор слоёв иногда можно автоматизировать библиотекой AutoKeras ©.

Как делать это быстро?

Параллельной версией умножения матриц по Штрассену (PDF).

умножения матриц по Штрассену

Судя по всему, надо не матрицы умножать, а именно сворачивать a*H*a’, где вектор a’ - это сопряженный к a. И скорее всего правильно так a’*H*a

надо ещё не забыть комплексное сопряжение в первом выражении

Хотя судя по результату (совпадает) оно уже учтено.

А в однопотоке?

А ещё там ранг 4, а матрицы могут быть и 20x20 в перспективе.

А для каждого вектора параметров своя матрица?

Приходите к нам на семинар. Только пока что все в разъездах, я чо то вообще болею после отпуска, башка чугунная… Ближе к осени.

Свернуть матрицу с параметрами это в данном случае умножить матрицу на вектор?:-)

Не знаю как правильно назвать и при этом коротко. В квантовой механике это бы называлось средним значением оператора M по состоянию V, и в Дираковской нотации бы записывалось как < V | M | V >. Но аналогие неточная, ибо M не является оператором.

А для каждого вектора параметров своя матрица?

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

А так-то да, разложить бы этот логарифм как-нибудь. Большой прорыв в науке будет.

А в однопотоке?

Многопоточить данные для однопоточных инструкций (Single instruction, multiple data - SIMD).

В рамках стандартного с++ смотреть в сторону Execution policy https://en.cppreference.com/w/cpp/algorithm/execution_policy_tag_t и алгоритмов использующих эти policy: reduce вместо accumulate, transform_reduce вместо inner_product и тд. Хотя все равно плохо получается.

reduce вместо accumulate, transform_reduce вместо inner_product

Хорошо, надо попробовать.

Я-то надеялся boost и Eigen компилироваться в код, который компилятор потом засунет в SIMD регистры.

Eigen (вроде) должен поддерживать. А boost::ublas - написан, чтобы было.

Хорошо. Эх, придётся код перелопачивать на Eigen, строчек-то не много с этой арифметикой, но чтобы оно считало правильно — это надо постараться.

Хотя вот даже по бенчу видно, что Eigen лучше работает. А пользует ли он при этом SIMD — так это надо в ассемблер глядеть.

А пользует ли он при этом SIMD

На словах https://eigen.tuxfamily.org/index.php?title=FAQ#Which_SIMD_instruction_sets_are_supported_by_Eigen.3F

так это надо в ассемблер глядеть.

Использовать-то использует, а насколько эффективно - я не эксперт.

SSE2 is generally enabled by default, but you can enable AVX and FMA for better performance

На некоторых машинках кластера есть AVX, так что придётся наверное ещё и ключики подсовывать разные при компиляции на кластере… зато может быть прирост!

что ты называешь сверткой? что такое ^dag? псевдообращение или сопряжение? что такое «комплексные параметры»? это матрицы? какого размера? вопрос то сформулируй.

а так, общо, чтобы что-либо численное делать быстро, надо помнить что ты никогда сам не напишешь хоть сколько-то эффективный код для линейной алгебры. только библиотеки, только blas/lapack/mkl

dag

dagger. Транспонированный и комплексно-сопряжённый вектор.

«Свёртка» — умножение V^dagger на матрицу и затем на V, то бишь V^dagger M V.

что такое «комплексные параметры»?

параметры функции правдоподобия, они же собираются вместе в вектор V.

какого размера?

У матриц размер порядка 10x10, может быть в двойку больше в перспективе. При этом ранг матриц — 4, и собственные числа неотрицательные.

только библиотеки, только blas/lapack/mkl

Дык я в тесте и сравниваю Eigen с boost::ublas. Можно ещё конечно лапак попробовать…

«Свёртка» — умножение V^dagger на матрицу и затем на V, то бишь V^dagger M V.

дружище, это называется квадратичной формой.

У матриц размер порядка 10x10

ну айген вроде как раз выбирает оптимальный алгоритм для константных размеров? Можешь для симметричных\эрмитовых матриц использовать соответствующие процедуры (почему ты так сделал в бусте но не в айгене?) http://eigen.tuxfamily.org/dox/group__QuickRefPage.html#title12 или дернуть бласовский dsymm.

Вообще, советую тебе проверить что ты слинковал все нужные либы https://eigen.tuxfamily.org/dox/TopicUsingIntelMKL.html , просто бывает х10 прирост от того что профессиональные библиотеки используешь (на твоих матрицах такого не будет канешн).

Далее, почему такой трепет перед скоростью умножения матричек 10-10? Кого и почему это должно ипать? Если у тебя их надо миллионы на каждой итерации умножать, то посмотри в сторону batched-операций. https://software.intel.com/content/www/us/en/develop/articles/introducing-batch-gemm-operations.html

При этом ранг матриц — 4

ну если ты так сильно угараешь про задротство и эффективность, то храни эти 4 вектора (можешь их через eigendecomposition/svd получить), которые дают тебе матрицу

M=sum_{i<=4} v^H_i * v_i

и матричные умножения можешь реализовать через векторные умножения (очень не советую таким заниматься – там жизнь за окном идет:)) ).

Можно ещё конечно лапак попробовать

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

ПС. Ты этого не спрашивал, но хочу тебе искренне посоветовать (так как сам в области твоего ОП работаю): не дрочи ты такие штуки на крестах, бери матлаб\питон\джулиу, там все это удобнее делать, а на джулии может конкретно быть быстрее (она сыроватая и непростая, но я начинаю верить что она перевернет все численное программирование) и у тебя будет больше времени думать о методе, чем о микросекундах и прочем крестовом задротстве Ну и если припечет, то просто напишешь свои ботлнеки на крестах. Да и на матлабе попрогать стоит хотя бы из-за тамошней документации :))) Думаю ты сам видишь уже, что численное программирование мало общего имеет с нормальным программированием, так что всему свой инструмент.

дружище, это называется квадратичной формой.

Хоть горшком назови, только в печку не суй. В тензорной форме так и так свёртка.

почему ты так сделал в бусте но не в айгене?

Потому что я прочитал только беглое введение в айген, решил вот затестить как он работает по сравнению с бустом.

Далее, почему такой трепет перед скоростью умножения матричек 10-10?

Ну ты видал же, 200 мсек на одну итерацию в таких тепличных условиях если юзать eigen, если буст — то 500, а в реальной задаче конечно хуже, потому что там минуит на каждую итерацию чего-то там ещё своего насчитывает и думает куда наступать дальше.

то храни эти 4 вектора

Придётся скорее всего так делать.

не дрочи ты такие штуки на крестах

Я долбанусь переписывать парциальные волны ещё раз на какой-нибудь джулии, у меня уже всё работает и насчитывается. Если только в helicity-формализм переходить, но этого я делать в ближайшие пару лет не собираюсь. Тем более это значит что надо изучать какой-то другой языг погроммирования, я уж лучше просто Eigen подучу.

матлабе

Проприетарщина же.

Хоть горшком назови, только в печку не суй.

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

В тензорной форме так и так свёртка.

Слово свертка гвоздями намертво прибито к операции свертки сигналов и соответствующим матрицам. Недаром на английском то что ты упоминаешь называется tensor contraction, а не convolution. А на русский, видимо, какой-то тракторист переводил.

Про остальное, у меня, конечно, тоже бывает нереальный стояк когда я своими числодробилками утилизирую каждый такт процессора, но надо понимать что тут нету предела совершенства. Если ты для дела этим занимаешься, то пока ты там на 20% сделаешь что-то эффективнее, другие ребята просто купят гпу или проц попиже. Ну или дунут, пока их алгоритм считается. Ты же не какой-то риалтайм пилишь?

Проприетарщина же.

У тебя это жизненная позиция? В матлабе много здравого, его точно стоит хотя бы потрогать)

ПС. Еще раз советую посмотреть batched linear algebra. Если повезет с типом задачи ты рили можешь до 10 раз на современных процах ускорить свои числодробилки. Я же правильно понимаю что тебе надо на каждой итерации посчитать просто 100_000 квадратичных форм (каждой матрице свой вектор да? просто я не разбираюсь в современных крестах уже совсем) от матрицы 10х10? Попробую завтра в пути глянуть как оно в МКЛе сделано.

ППС. Добавь сейфгард в логарифм log(говно + EPS), ведь у тебя там близко к нулю может быть, все будет численно неустойчиво.

Вектор для одной итерации один, а матрицы разные.

Да, про билинейную форму я чёт подзабыл, это было на первом семестре, а потом нигде не надо было особо помнить название, а уж тем более применять на практике.

ГПУ конечно хочется купить, но в нашей лабе деньги уходят точно не на компьютинг, дескать и двух сотен оптеронов хватит.

batched linear algebra постараюсь посмотреть.

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

Вектор для одной итерации один, а матрицы разные.

Необычный вариант. Я сейчас сходу не прикину в голове, можно ли это как-то преобразовать, но я бы попробовал бенчмаркнуть einsum из нампи. Только пиши операцию сразу для массива из 100_000 матриц. Там он вроде должен хитрые оптимизации таких сумм делать. В этом, кстати, сила джулии, она сильно строже пистонов всяких, но юзабельнее крестов. И из-за тайп стабилити там ллвм оптимизирует как бес иногда.

Да, про билинейную форму я чёт подзабыл,

да я сам линал понял по-настоящему только когда ученым мужем стал)

Близко к нулю у меня не бывает, точнее если бывает, то это значит модель косячная.

Из своего опыта, там столько тонкостей с флоатами, что советую всегда пользоваться сейфгардами. Искать место откуда просираются наны или инфы – то еще веселье.

ГПУ конечно хочется купить, но в нашей лабе деньги уходят точно не на компьютинг, дескать и двух сотен оптеронов хватит.

да с этим у всех у нас проблемы. Мы у гугла стандартно клянчим денюжки на их GCP и гпушки в клауде. В россии, к сожалению эти программы не действуют (по крайней мере у гугла). Можешь посмотреть google co-lab, хотя думаю и так знаешь.

Твою задачу рили гляну плотнее на недельке, у меня самого давно была идея batched операции из МКЛ потрогать и заюзать у себя.

Искать место откуда просираются наны или инфы – то еще веселье

Дык я делаю по 50 фитов для одного только бина по массе (а всего их 40), ну и кто улетел в бесконечность — того и проблемы. На самом деле когда я отладил матрицы, чтобы они правильно считались, проблемы с нанами ещё ни разу не всплыли. Ну если выпадет один фит из 50 в бесконечность — пофиг. Главное чтобы хотя бы несколько минимум нашли.

но я бы попробовал бенчмаркнуть einsum из нампи

У меня обычно мозг от нумпая сворачивается, постоянно лезу гуглить как и что там устроено. Но это потому что не очень часто им пользуюсь, только лишь пока для отрисовки в matplotlib ну и сопутствующих посчитать среднее или построить гистограмму.

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

Ты все напутал, на плюсах числодробят зеленые и красные, причем никаких бенчей для определения победителя там не надо.

Чёрт, unsequential_policy только в 20 плюсах, а на кластере только gcc 9.