Разреженное матрично-плотное векторное умножение с матрицей, известной во время компиляции

Однако я не уверен, почему - настолько ли размер кода замедляет его?

Проблема в том, что исполняемый файл большой, ОС будет загружать много страниц с вашего устройства хранения. Этот процесс очень медленный. Процессор часто зависает в ожидании загрузки данных. И даже если бы код был уже загружен в ОЗУ (кеширование ОС), это было бы неэффективно, потому что скорость ОЗУ (латентность + пропускная способность) совсем плохая. Основная проблема здесь в том, что все инструкции выполняются только один раз. Если вы повторно используете функцию, то код необходимо перезагрузить из кеша, а если он слишком велик, чтобы поместиться в кеш, он будет загружен из медленной оперативной памяти. Таким образом, накладные расходы на загрузку кода очень высоки по сравнению с его фактическим выполнением. Чтобы решить эту проблему, вам нужно использовать довольно небольшой код с циклами, повторяющимися на довольно небольшом количестве данных.

Исключает ли такой код оптимизацию компилятора?

Это зависит от компилятора, но большинство основных компиляторов (например, GCC или Clang) оптимизируют код таким же образом (отсюда и медленное время компиляции).

Как вы думаете, стоит ли это попробовать (думаю, я все равно попробую в понедельник)?

Да, это решение явно лучше, особенно если индексы хранить компактно. В вашем случае вы можете хранить их, используя тип uint16_t. Все индексы можно поместить в большой буфер. Начальное/конечное положение индексов для каждой строки можно указать в другом буфере, ссылаясь на первый (или используя указатели). Этот буфер можно загрузить один раз в память в начале вашего приложения из специального файла, чтобы уменьшить размер результирующей программы (и избежать выборки с устройства хранения в критическом цикле). С вероятностью 0,01% наличия ненулевых значений результирующая структура данных займет менее 500 КиБ ОЗУ. На среднестатистическом процессоре для настольных ПК он может поместиться в кэше L3 (что довольно быстро), и я думаю, что ваши вычисления не должны занимать более 1 мс, если предположить, что код multiply тщательно оптимизирован.

Другая идея состоит в том, чтобы попытаться сохранить входной (и, возможно, также выходной?) вектор в виде упакованных битов и выполнить вычисления таким образом.

Битовая упаковка хороша только в том случае, если ваша матрица не слишком разрежена. С матрицей, заполненной 50% ненулевыми значениями, метод битовой упаковки великолепен. С 0,01% ненулевых значений метод битовой упаковки явно плох, так как занимает слишком много места.

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

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

Вот примерные тайминги для различных операций, которые может выполнять компьютер:

Я реализовал второй метод (массивы constexpr, хранящие матрицу в формате хранения сжатых столбцов), и он намного лучше. Для компиляции с -O3 требуется (для двоичной матрицы 64 000 x 22 000, содержащей 35 000 единиц) 1 минута, а на моем ноутбуке выполняется одно умножение за ‹300 микросекунд (у Джулии на то же вычисление уходит около 350 микросекунд). Общий размер исполняемого файла составляет ~1 Мбайт.

Вероятно, можно еще сделать намного лучше. Если у кого-то есть идея, дайте мне знать!

Ниже приведен пример кода (показывающий матрицу 5x10), иллюстрирующий то, что я сделал.

#include <iostream>
#include <array>

// Compressed sparse column storage for binary matrix
constexpr std::size_t M = 5;
constexpr std::size_t N = 10;
constexpr std::size_t num_nz = 5;
constexpr std::array<std::uint16_t, N + 1> colptr = {
0x0,0x1,0x2,0x3,0x4,0x5,0x5,0x5,0x5,0x5,0x5
};
constexpr std::array<std::uint16_t, num_nz> row_idx = {
0x0,0x1,0x2,0x3,0x4
};

template<typename Bit>
constexpr void encode(const std::array<Bit, N>& in, std::array<Bit, M>& out) {

    for (std::size_t col = 0; col < N; col++) {
        for (std::size_t j = colptr[col]; j < colptr[col + 1]; j++) {
            out[row_idx[j]] = (static_cast<bool>(out[row_idx[j]]) != static_cast<bool>(in[col]));
        }
    }
}

int main() {
    using Bit = bool;
    std::array<Bit, N> input{1, 0, 1, 0, 1, 1, 0, 1, 0, 1};
    std::array<Bit, M> output{};
    
    for (auto i : input) std::cout << i;
    std::cout << std::endl;

    encode(input, output);

    for (auto i : output) std::cout << i;
}