Как оптимизировать этот продукт из трех матриц в C++ для x86?

Question

У меня есть ключевой алгоритм, в котором большая часть времени выполнения тратится на вычисление продукта с плотной матрицей:

A*A'*Y, where: A is an m-by-n matrix, 
               A' is its conjugate transpose,
               Y is an m-by-k matrix

Typical characteristics:
    - k is much smaller than both m or n (k is typically < 10)
    - m in the range [500, 2000]
    - n in the range [100, 1000]

Исходя из этих измерений, в соответствии с уроками задачи умножения цепочки матриц, становится ясно, что в смысле числа операций оптимально структурировать вычисление как A*(A'*Y), Моя текущая реализация делает это, и заметно повышение производительности от простого навязывания этой ассоциативности выражению.

Мое приложение написано на C++ для платформы x86_64. Я использую библиотеку линейной алгебры Eigen с библиотекой Intel Math Kernel в качестве бэкэнда. Eigen может использовать интерфейс BLAS IMKL для выполнения умножения, и также важен переход от собственной реализации SSE2 Eigen к оптимизированной реализации Intel на основе AVX на моей машине Sandy Bridge.

Тем не менее, выражение A * (A.adjoint() * Y) (написано на собственном языке) раскладывается на два общих матрично-матричных произведения (вызовы xGEMM Подпрограмма BLAS), с временной матрицей, созданной между ними. Мне интересно, если, перейдя к специализированной реализации для оценки всего выражения сразу, я смогу получить реализацию, которая быстрее, чем общая, которая у меня есть сейчас. Пара наблюдений, которые заставляют меня поверить в это:

• Используя типичные размеры, описанные выше, входная матрица A обычно не помещается в кэш. Следовательно, конкретный шаблон доступа к памяти, используемый для вычисления трехматричного произведения, будет ключевым. Очевидно, что избежать создания временной матрицы для частичного продукта также было бы выгодно.

• A и его сопряженная транспонирование, очевидно, имеет очень связанную структуру, которую можно использовать для улучшения схемы доступа к памяти для всего выражения.

Существуют ли какие-либо стандартные методы для реализации такого рода выражений в дружественной кэш-памяти форме? Большинство методов оптимизации, которые я нашел для умножения матриц, предназначены для стандартных A*B регистр, не большие выражения. Меня устраивают аспекты этой проблемы, связанные с микрооптимизацией, такие как перевод в соответствующие наборы SIMD-инструкций, но я ищу какие-либо ссылки, чтобы разбить эту структуру наиболее удобным для памяти способом.

Изменить: Основываясь на ответах, которые пришли до сих пор, я думаю, что я был немного неясным выше. Тот факт, что я использую C++/Eigen, на самом деле является лишь подробностью реализации с моей точки зрения на эту проблему. Eigen отлично справляется с реализацией шаблонов выражений, но оценка этого типа проблемы как простого выражения просто не поддерживается (только продукты из 2 общих плотных матриц).

На более высоком уровне, чем то, как выражения будут оцениваться компилятором, я ищу более эффективную математическую разбивку операции сложного умножения, с целью избежать ненужных избыточных обращений к памяти из-за общей структуры A и его сопряженное транспонирование. Результат, вероятно, будет трудно реализовать эффективно в чистом Эйгене, поэтому я, скорее всего, просто реализую его в специализированной процедуре с использованием SIMD-функций.

Answer 1

Это не полный ответ (пока - и я не уверен, что он станет одним).

Давайте сначала немного подумаем о математике. Поскольку матричное умножение ассоциативно, мы можем сделать (A*A')Y или A(A' * Y).

Операции с плавающей точкой для (A*A')*Y

2*m*n*m + 2*m*m*k //the twos come from addition and multiplication

Операции с плавающей точкой для A*(A'*Y)

2*m*n*k + 2*m*n*k = 4*m*n*k

Поскольку k намного меньше, чем m и n, понятно, почему второй случай намного быстрее.

Но симметрично мы могли бы в принципе уменьшить количество вычислений для A*A'на два (хотя это может быть нелегко сделать с SIMD), поэтому мы могли бы уменьшить количество операций с плавающей запятой (A * A') * Y в

m*n*m + 2*m*m*k.

Мы знаем, что и m, и n больше, чем k. Давайте выберем новую переменную для m и n с именем z и выяснить, где случай один и два равны:

z*z*z + 2*z*z*k = 4*z*z*k  //now simplify
z = 2*k.

Таким образом, если m и n больше, чем дважды, во втором случае будет меньше операций с плавающей запятой. В вашем случае m и n оба больше 100, а k меньше 10, поэтому в втором случае используется гораздо меньше операций с плавающей запятой.

С точки зрения эффективного кода. Если код оптимизирован для эффективного использования кэша (как это делают MKL и Eigen), то умножение больших плотных матриц связано с вычислениями, а не с памятью, поэтому вам не нужно беспокоиться о кеше. MKL быстрее, чем Eigen, так как MKL использует AVX (а может быть, fma3 сейчас?).

Я не думаю, что вы сможете сделать это более эффективно, чем вы уже используете второй случай и MKL (через Eigen). Включите OpenMP, чтобы получить максимальные FLOPS.

Вы должны рассчитать эффективность, сравнивая FLOPS с пиковыми FLOPS вашего процессора. Предполагая, что у вас есть процессор Sandy Bridge/Ivy Bridge. Пик SP FLOPS составляет

frequency * number of physical cores * 8 (8-wide AVX SP) * 2 (addition + multiplication)

Для двойной прецессии делим на два. Если у вас есть Haswell и MKL использует FMA, тогда удвойте пиковые значения FLOPS. Чтобы получить правильную частоту, вы должны использовать значения Turbo Boost для всех ядер (это ниже, чем для одного ядра). Вы можете посмотреть это, если вы не разогнали свою систему или используете CPU-Z в Windows или Powertop в Linux, если у вас разогнанная система.

Answer 2

Используйте временную матрицу для вычисления A'*Y, но убедитесь, что вы говорите eigen, что псевдонимов не происходит: temp.noalias() = A.adjoint()*Y, Затем вычислите ваш результат, еще раз сказав, что объекты не являются псевдонимами: result.noalias() = A*temp,

Answer 3

Было бы избыточное вычисление, только если вы выполняете (A*A')*Y так как в этом случае (A*A') является симметричным, и требуется только половина вычислений. Однако, как вы заметили, это все еще намного быстрее A*(A'*Y) в этом случае нет избыточных вычислений. Я подтверждаю, что стоимость временного создания здесь совершенно незначительна.

Answer 4

Я думаю, что выполнить следующее

result = A * (A.adjoint() * Y)

будет так же, как это сделать

temp = A.adjoint() * Y
result = A * temp;

---

Если ваша матрица Y помещается в кэш, вы, вероятно, можете воспользоваться этим

result = A * (Y.adjoint() * A).adjoint()

или, если предыдущая запись не допускается, вот так

temp = Y.adjoint() * A
result = A * temp.adjoint();

Тогда вам не нужно делать сопряжение с матрицей A и хранить временную сопряженную матрицу для A, что будет намного дороже, чем для Y.

Если ваша матрица Y помещается в кэш, то она должна выполняться намного быстрее, выполняя цикл по столбцам A для первого умножения, а затем по строкам A для второго умножения (с Y.adjoint() в кэше для первое умножение и temp.adjoint() для второго), но я думаю, что внутренне собственный уже заботится об этих вещах.