Почему умножение матриц 4x4 в Eigen происходит более чем в два раза быстрее, чем 3x3?

TL;DR: одновременно происходит несколько проблем, но главная из них заключается в том, что 3D-версия имеет более высокую задержку (из-за проблем с планированием и файлом регистров), которую нельзя перекрыть с вычислениями из-за зависимости, переносимой в цикле. Изменение порядка операций помогает ЦП скрыть задержку 3D-версии, поэтому она может быть быстрее, чем 4D. Другие проблемы, влияющие на стабильность теста, связаны с масштабированием частоты (и, в более общем плане, управлением питанием) и выравниванием памяти.

Установка хорошего эталона

Я могу воспроизвести этот процессор AMD Ryzen 7 5700U (Zen2) с GCC 13.2.0 на Ubuntu 24.04. Результаты производительности за один прогон очень нестабильны (поскольку бенчмарк слишком короткий). При этом при многократном запуске программы (например, 1000 раз) начинают происходить интересные вещи и неожиданно появляется странный эффект!

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

Обратите внимание, что вертикальная ось — это сообщаемое время (в мкс), а горизонтальная ось — количество итераций.

Мы видим, что 3D-версия менее стабильна, чем 4D-версия, но шум настолько велик, что мы не можем ничего сделать (может просто 3D менее стабильно, чем 4D и есть какой-то странный порог в таймингах) 3D-версии).

Если увеличить количество запусков до 1000, то получим удивительные результаты:

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

Теперь мы можем сказать, что 3D-версия действительно в среднем медленнее и менее стабильна, чем 4D-версия. Интересно, что мы видим, что 4D-версия через определенное время работает на удивление быстрее, чем 3D-версия.

Несвязанные параметры

Обратите внимание, что использование -march=native ускоряет вычисления (благодаря AVX/AVX2), но не влияет на проблемы в целом. Также обратите внимание, что -march=native и -mavx2 -mfma дают одинаковые результаты производительности.

Обратите внимание, что аналогичный эффект можно увидеть и с векторами, в 20 раз большими. Эффект также не зависит от привязки потока, поэтому это определенно не связано с переключением контекста или эффектами NUMA. Изменение порядка вызова функций не привело к существенному влиянию на производительность. Регулятор масштабирования установлен на performance на всех ядрах (хотя частота может быть переменной). Может ли это быть из-за масштабирования частоты?

Влияние частоты (масштабирования) и управления питанием

Что ж, похоже, масштабирование частоты каким-то образом играет роль в начальном пороге для 4D-версии. Я установил минимальную и максимальную частоту на одно и то же значение (800 МГц на моей машине из-за ошибки в регуляторе масштабирования на моей машине, который не заботится о минимальной частоте, которую я действительно могу установить, поэтому минимальная частота имеет тенденцию всегда составлять 800 МГц для активные ядра). Проверил, чтобы частота была стабильной во время прогона (по крайней мере, с грубой детализацией). Вот результат:

Мы видим, что для 4D-версии порога больше нет! В более общем плане я заметил два эффекта при попытке использовать разные частоты:

• более низкая частота уменьшает разрыв между высокими таймингами (первые) и низкими таймингами (последние после порога) 4D-версии;
• более низкая частота уменьшает пороговую задержку (т.е. количество итераций перед порогом);
• Большее количество элементов в векторах также приводит к меньшей пороговой задержке.

Что касается двух последних пунктов, то на самом деле это связано с тем, что пороговая задержка представляет собой фиксированную пороговую задержку, независимую от частоты (около 2,5 секунды), поэтому более низкая частота делает каждую итерацию медленнее, вызывая меньшую задержку. Все это доказывает, что частота влияет на порог, хотя я не могу объяснить, почему. Я предполагаю, что это связано с управлением питанием процессоров AMD Zen. Я ожидаю, что эта константа будет либо установлена ​​в ядре (возможно, в файлах, зависящих от поставщика), либо константа оборудования/прошивки, основанная на частоте вне ядра (независимо от частоты ядра).

Я также обнаружил (недавно), что эффекта не происходит, когда ноутбук заряжается, хотя производительность немного менее стабильна! Это подтверждает, что эта проблема напрямую связана с управлением питанием ЦП.

Влияние выравнивания памяти

Когда время выполнения функции, выполняющей базовые числовые вычисления, квантуется (и эффект не связан с частотой), квантование часто вызывается проблемами выравнивания памяти. Это то, что происходит с 3D-версией, обычно потому, что матрица двойной точности 3x3 занимает 3*3*8=72 байт (как указано Raildex в комментариях). На практике компиляторы, такие как GCC, могут добавлять некоторые дополнения в конце структуры, чтобы доступ к памяти обычно был быстрее (что приводит к 80-байтовым матрицам 3x3). Обратите внимание, что это не степень двойки. Более того, матрица 3x3 обычно пересекает две строки кэша, что приводит к дополнительной загрузке кэша L1 во время невыровненного доступа. Такая дополнительная нагрузка приводит к увеличению задержки. Хуже того: накладные расходы на задержку не являются постоянными, поэтому я думаю, что планировщику инструкций ЦП может быть сложнее эффективно планировать целевые инструкции. Чтобы проверить гипотезу выравнивания памяти, мы можем просто выровнять каждую матрицу по строке кэша. Это не очень эффективно с точки зрения объема памяти (поэтому может снизить производительность кодов, привязанных к памяти), но это должно, по крайней мере, решить проблему выравнивания. Вот код для принудительного выравнивания матриц в памяти:

// Align each matrix to 64 bytes (i.e. a cache line)
class alignas(64) Matrix3d : public Eigen::Matrix3d
{
public:
    Matrix3d() = default;
    Matrix3d(const Matrix3d&) = default;
    Matrix3d(Matrix3d&&) = default;
    template<typename MatrixType>
    Matrix3d(const MatrixType& other) : Eigen::Matrix3d(other) {}
    template<typename MatrixType>
    Matrix3d& operator=(const MatrixType& other) { Eigen::Matrix3d::operator=(other); return *this; }
};

class alignas(64) Matrix4d : public Eigen::Matrix4d
{
public:
    Matrix4d() = default;
    Matrix4d(const Matrix4d&) = default;
    Matrix4d(Matrix4d&&) = default;
    template<typename MatrixType>
    Matrix4d(const MatrixType& other) : Eigen::Matrix4d(other) {}
    template<typename MatrixType>
    Matrix4d& operator=(const MatrixType& other) { Eigen::Matrix4d::operator=(other); return *this; }
};

// Replace all occurrence of Eigen::Matrix by Matrix in the rest of the code
// [...]

Вот результаты производительности:

Мы видим, что тайминги 3D-версии больше не квантуются. Обратите внимание, что, к сожалению, время соответствует верхнему порогу в предыдущих тестах. Также можно отметить, что 4D-версия немного медленнее (возможно, из-за менее эффективного генерируемого кода). Как ни странно, синхронизация 4D-версии менее стабильна. Я ожидаю, что выравнивание массива, выделенного std::vector, повлияет на такое время.

В сочетании с фиксированной (низкой) частотой мы получаем следующие результаты производительности:

Результаты сейчас стабильны, но остается одна большая проблема: почему 4D-версия по-прежнему быстрее 3D-версии?

Планирование инструкций, задержка и зависимости

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

Планирование инструкций может оказать сильное влияние на производительность. Это особенно актуально для процессора Zen2, который AFAIK использует несколько отдельных планировщиков FP (а не единый планировщик). Более того, некоторые инструкции могут быть запланированы только для определенных портов. Таким образом, если планировщик примет неправильное решение, не поместив инструкции в лучшую очередь для данной итерации, задержка может немного увеличиться (вероятно, на несколько циклов), поэтому общее выполнение будет медленнее. Несколько циклов кажутся небольшими, но несколько циклов на итерацию, когда каждая итерация длится десятки циклов, — это много. Здесь, в последнем тесте на частоте 800 МГц, мы видим, что 4D-версии требуется около 4800 мкс для выполнения 100_000 итераций. Это означает ~38 циклов/итераций (при этом необходимо выполнить около 140-150 инструкций/итераций)!

На практике счетчики производительности оборудования, как правило, указывают на плохое использование внутренних ресурсов (в модифицированном примере, предназначенном для измерения каждой версии отдельно):

Use of the 4 FP SIMD ports (i.e. execution units):
  - 3D version:
       467 703 544      fpu_pipe_assignment.total0                                            
       940 583 576      fpu_pipe_assignment.total1                                            
        22 580 911      fpu_pipe_assignment.total2                                            
        10 002 224      fpu_pipe_assignment.total3  
  - 4D version:
       334 752 204      fpu_pipe_assignment.total0                                            
       280 487 188      fpu_pipe_assignment.total1                                            
       269 800 454      fpu_pipe_assignment.total2                                            
       268 277 318      fpu_pipe_assignment.total3     

В 3D-версии два порта практически не используются вообще. Из двух оставшихся портов один насыщен в два раза больше другого! Это действительно плохо. Между тем, 4D-версия работает хорошо. Это имеет тенденцию подтверждать проблемы с планированием инструкций.

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

Source of stalls (in cycles):
  - 3D version:
     7 936 000 400      cycles                                                                
     1 609 456 587      de_dis_dispatch_token_stalls1.fp_sch_rsrc_stall                                      
     3 998 358 801      de_dis_dispatch_token_stalls1.fp_reg_file_rsrc_stall                                       
  - 4D version:
     6 950 481 062      cycles                                                                
       194 523 072      de_dis_dispatch_token_stalls1.fp_sch_rsrc_stall                                      
     2 858 781 759      de_dis_dispatch_token_stalls1.fp_reg_file_rsrc_stall                                      

fp_sch_rsrc_stall укажите «циклы, в которых группа отправки действительна, но не отправляется из-за задержки токена. Задержка ресурсов планировщика FP.». Хотя fp_reg_file_rsrc_stall аналогичен, но означает «остановка ресурса файла регистра с плавающей запятой». Таким образом, в 3D-версии, хотя некоторые задержки (20%) происходят из-за проблем с планированием, основной источник задержек (50%), по-видимому, исходит из файла регистров FP. Обычно это происходит из-за отсутствия записей в файле регистров FP, но здесь это было бы странно, поскольку предполагается, что Zen2 имеет 160 регистров FP (что кажется достаточно большим для этого кода). Интересно, может ли это быть каким-то образом связано со спекулятивным выполнением инструкций FP. В версии 4D также наблюдается большое количество зависаний регистров FP (41%) и почти полное отсутствие зависаний из-за планирования. Хотя я не до конца понимаю причину, по которой так много зависаний в файлах регистров FP, последний момент подтверждает приведенную выше гипотезу.

Более быстрый/лучший код

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

Более того, зависимость между итерациями обычно снижает производительность. Ключом к ускорению этих вычислений является вычисление попарного сокращения (то есть что-то вроде (A*A)*(A*A)). Это помогает ЦП выполнять несколько умножений матриц одновременно (возможно, даже параллельно), уменьшая задержки в планировании. Кроме того, попарная редукция гораздо более устойчива численно. Это улучшение делает 3D-версию быстрее по сравнению с 4D-версией! Вот результаты производительности при вычислении матриц попарно (для обеих версий):

Мы видим, что 3D-версия теперь работает быстрее, а 4D-версия не сильно затронута этим улучшением. Обратите внимание, что вычисление матриц по группе из 4 снижает производительность, поскольку сгенерированный код кажется менее эффективным для моего процессора (из-за планирования инструкций и слишком большого количества требуемых регистров).

Инструкции FMA (т. е. инструкции, объединяющие умножение и сложение) могут повысить точность вычислений, одновременно уменьшая задержку кода. Вы можете указать -mfma компилятору, чтобы сгенерировать инструкции FMA. Однако, как указывалось ранее, это не влияет на разницу между версией 3D и 4D.

Добавление заполнения в 3D-матрицу может помочь выполнить возможно меньшую (выровненную) загрузку/сохранение, избегая некоторых накладных расходов на задержку (AFAIK, хранящий данные в строке кэша, может повлиять на скорость чтения в одной и той же строке кэша, но в разных местах на некоторых процессорах) . Компиляторы, как правило, не делают этого, поскольку им необходимо доказать, что доступ безопасен и операции не вызывают каких-либо побочных эффектов для незамаскированных значений (не говоря уже о том, что маскирование может быть дорогостоящим на некоторых платформах).

Для базовых матриц 3x3 и 4x4 лучшим решением может быть написание вручную ассемблерного кода для целевой архитектуры, несмотря на множество недостатков (непростой, утомительный в написании, не очень удобный в обслуживании, непереносимый).