Матричное умножение собственной матрицы для подмножества столбцов

Решение проблемы

Использование симметрии

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

Mat sub_mat = subset_cols(mat, idx); // From your original post
Mat a = Mat::Zero(numRows, numRows);
a.selfadjointView<Eigen::Lower>().rankUpdate(sub_mat); // (1)
a.triangularView<Eigen::Upper>() = a.transpose(); // (2)


Линия (1)будет вычисляться a += sub_mat * sub_mat.transpose()только для нижней части. (2)затем запишет нижнюю часть в верхнюю часть. Также смотрите документацию ( здесь и здесь ). Конечно, если вы можете жить только с нижней частью, шаг (2) можно пропустить.

Для матрицы 100x100000 matя получаю ускорение примерно в 1 раз.

~ 1,1x при 10 столбцах,

~ 1,5x при 100 столбцах,

~1,7x при использовании 1000 столбцов

как в Windows с использованием MSVC, так и в Linux с использованием clang с полной оптимизацией и AVX.

Включение распараллеливания

Еще один способ ускорить вычисления — включить распараллеливание путем компиляции с OpenMP. Эйген позаботится обо всем остальном. Однако приведенный выше код, использующий симметрию, не выигрывает от этого. Но исходный код

Eigen::MatrixXd sub_mat = subset_cols(mat, idx);
Eigen::MatrixXd a = sub_mat * sub_mat.transpose();


делает.

Для матрицы 100x100000 matс помощью clang в Linux, работающей с 4 потоками (на 4 реальных ядрах) и по сравнению с одним потоком, я получаю ускорение примерно в 1 раз.

~1.0x при использовании 10 столбцов, т.е. никакого ускорения.

~ 1,8x при 100 столбцах

~ 2,0x при использовании 1000 столбцов

Другими словами, 4 ядра или более превосходят симметричный метод, показанный выше, за исключением очень небольшого количества столбцов. Использование только 2 ядер всегда было медленнее. Обратите внимание, что использование SMT ухудшало производительность в моих тестах, иногда заметно.

Другие примечания

Я уже писал это в комментарии, но для полноты картины:
Eigen::Mapне получится, потому что шаги не равноудалены. Использование нарезки дает мне примерно на 10% лучшую производительность, чем ваш метод копирования в Linux с clang и gcc, но несколько хуже на MSVC. Кроме того, как вы заметили, он недоступен в ветке Eigen 3.3. Существует собственный способ имитировать его, но в моих тестах он всегда работал хуже. Кроме того, в моих тестах он не экономил память по сравнению с методом копирования.

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

Кроме того, следует не только включить оптимизации, но и использовать максимально возможный набор инструкций. Включение AVX в моих тестах улучшило производительность примерно в 1,5 раза. К сожалению, я не могу протестировать AVX512.