Вопрос по linear-algebra, algorithm, matrix-multiplication, c++ – Матрица умножения Ладермана 3х3 с 23 умножениями, стоит ли это того?

18

Возьмите произведение двух матриц 3х3A*B=C, Наивно это требует 27 умножений, используястандартный алгоритм, Если бы вы были умны, вы могли бы сделать это, используя только 23 умножения,результат, найденный в 1973 году Ладерманом, Техника предполагает сохранение промежуточных шагов и их правильное объединение.

Теперь давайте исправим язык и тип, скажем, C ++ с элементамиdouble, Если бы алгоритм Ладермана был жестко запрограммирован по сравнению с простым двойным циклом, можем ли мы ожидать, что производительность современного компилятора уменьшит различия в алгоритмах?

Notes about this question: Этоprogramming сайт, и вопрос задается в контексте наилучшей практики для критичного ко времени внутреннего цикла; преждевременной оптимизации это не так. Советы по реализации приветствуются в качестве комментариев.

Общий ответ на этот вопрос - «нет». Умные алгоритмы все еще нужны в мире. phs
Поздний комментарий. Преимущества алгоритма Ладермана становятся видимыми, только если он используется рекурсивно для больших матриц. Это было экспериментально проанализировано Паоло Д. Альберто в егоFastMMW site, С точки зрения сложности умножения алгоритм Ладермана O (n ^ 2.854) лучше, чем наивный O (n ^ 3), но уступает алгоритму O Штрассена (n ^ 2.807). Алгоритм 3х3 с 21 умножением O (n ^ 2,77) будет лучше, чем у Штрассена. Axel Kemper
@phs: Ответ на вопрос в заголовке или на вопрос чуть выше заметки? Они противоположны. MSalters
& quot; ... можно ли ожидать, что производительность современного компилятора устранит различия в алгоритмах? & quot; Почему бы не попробовать это? Кодируйте два, запускайте их каждые 1000 раз и сравните время выполнения. AndyPerfect
Существует также метод Макарова, который требует только 22 умножения (original in Russian, PDF; English translation, behind paywall) Jitse Niesen

Ваш Ответ

4   ответа
10

Timing Tests:

Я сам проверил временные тесты, и результаты меня удивили (поэтому я и задал вопрос в первую очередь). Суть в том, что при стандартной компиляцииladerman на ~ 225% быстрее, но с-03 оптимизирующий флаг это50% slower! Я должен был добавить случайный элемент в матрицу каждый раз во время-O3 флаг или компиляторcompletely optimized away простое умножение, занимающее нулевое время с точностью до часов. Посколькуladerman Алгоритм был трудным проверить / перепроверить. Я опубликую полный код ниже для потомков.

Спецификации: Ubuntu 12.04, Dell Prevision T1600, gcc. Процентная разница во времени:

  • g++ [2.22, 2.23, 2.27]
  • g++ -O3 [-0.48, -0.49, -0.48]
  • g++ -funroll-loops -O3 [-0.48, -0.48, -0.47] 

Код бенчмаркинга вместе с реализацией Ладермана:

#include <iostream>
#include <ctime>
#include <cstdio>
#include <cstdlib>
using namespace std;

void simple_mul(const double a[3][3], 
        const double b[3][3],
        double c[3][3]) {
  int i,j,m,n;
  for(i=0;i<3;i++) {
    for(j=0;j<3;j++) {
      c[i][j] = 0;
      for(m=0;m<3;m++) 
    c[i][j] += a[i][m]*b[m][j];
    }
  }
}

void laderman_mul(const double a[3][3], 
           const double b[3][3],
           double c[3][3]) {

   double m[24]; // not off by one, just wanted to match the index from the paper

   m[1 ]= (a[0][0]+a[0][1]+a[0][2]-a[1][0]-a[1][1]-a[2][1]-a[2][2])*b[1][1];
   m[2 ]= (a[0][0]-a[1][0])*(-b[0][1]+b[1][1]);
   m[3 ]= a[1][1]*(-b[0][0]+b[0][1]+b[1][0]-b[1][1]-b[1][2]-b[2][0]+b[2][2]);
   m[4 ]= (-a[0][0]+a[1][0]+a[1][1])*(b[0][0]-b[0][1]+b[1][1]);
   m[5 ]= (a[1][0]+a[1][1])*(-b[0][0]+b[0][1]);
   m[6 ]= a[0][0]*b[0][0];
   m[7 ]= (-a[0][0]+a[2][0]+a[2][1])*(b[0][0]-b[0][2]+b[1][2]);
   m[8 ]= (-a[0][0]+a[2][0])*(b[0][2]-b[1][2]);
   m[9 ]= (a[2][0]+a[2][1])*(-b[0][0]+b[0][2]);
   m[10]= (a[0][0]+a[0][1]+a[0][2]-a[1][1]-a[1][2]-a[2][0]-a[2][1])*b[1][2];
   m[11]= a[2][1]*(-b[0][0]+b[0][2]+b[1][0]-b[1][1]-b[1][2]-b[2][0]+b[2][1]);
   m[12]= (-a[0][2]+a[2][1]+a[2][2])*(b[1][1]+b[2][0]-b[2][1]);
   m[13]= (a[0][2]-a[2][2])*(b[1][1]-b[2][1]);
   m[14]= a[0][2]*b[2][0];
   m[15]= (a[2][1]+a[2][2])*(-b[2][0]+b[2][1]);
   m[16]= (-a[0][2]+a[1][1]+a[1][2])*(b[1][2]+b[2][0]-b[2][2]);
   m[17]= (a[0][2]-a[1][2])*(b[1][2]-b[2][2]);
   m[18]= (a[1][1]+a[1][2])*(-b[2][0]+b[2][2]);
   m[19]= a[0][1]*b[1][0];
   m[20]= a[1][2]*b[2][1];
   m[21]= a[1][0]*b[0][2];
   m[22]= a[2][0]*b[0][1];
   m[23]= a[2][2]*b[2][2];

  c[0][0] = m[6]+m[14]+m[19];
  c[0][1] = m[1]+m[4]+m[5]+m[6]+m[12]+m[14]+m[15];
  c[0][2] = m[6]+m[7]+m[9]+m[10]+m[14]+m[16]+m[18];
  c[1][0] = m[2]+m[3]+m[4]+m[6]+m[14]+m[16]+m[17];
  c[1][1] = m[2]+m[4]+m[5]+m[6]+m[20];
  c[1][2] = m[14]+m[16]+m[17]+m[18]+m[21];
  c[2][0] = m[6]+m[7]+m[8]+m[11]+m[12]+m[13]+m[14];
  c[2][1] = m[12]+m[13]+m[14]+m[15]+m[22];
  c[2][2] = m[6]+m[7]+m[8]+m[9]+m[23];    
}

int main() {
  int N = 1000000000;
  double A[3][3], C[3][3];
  std::clock_t t0,t1;
  timespec tm0, tm1;

  A[0][0] = 3/5.; A[0][1] = 1/5.; A[0][2] = 2/5.;
  A[1][0] = 3/7.; A[1][1] = 1/7.; A[1][2] = 3/7.;
  A[2][0] = 1/3.; A[2][1] = 1/3.; A[2][2] = 1/3.;

  t0 = std::clock();
  for(int i=0;i<N;i++) {
    // A[0][0] = double(rand())/RAND_MAX; // Keep this in for -O3
    simple_mul(A,A,C);
  }
  t1 = std::clock();
  double tdiff_simple = (t1-t0)/1000.;

  cout << C[0][0] << ' ' << C[0][1] << ' ' << C[0][2] << endl;
  cout << C[1][0] << ' ' << C[1][1] << ' ' << C[1][2] << endl;
  cout << C[2][0] << ' ' << C[2][1] << ' ' << C[2][2] << endl;
  cout << tdiff_simple << endl;
  cout << endl;

  t0 = std::clock();
  for(int i=0;i<N;i++) {
    // A[0][0] = double(rand())/RAND_MAX; // Keep this in for -O3
    laderman_mul(A,A,C);
  }
  t1 = std::clock();
  double tdiff_laderman = (t1-t0)/1000.;

  cout << C[0][0] << ' ' << C[0][1] << ' ' << C[0][2] << endl;
  cout << C[1][0] << ' ' << C[1][1] << ' ' << C[1][2] << endl;
  cout << C[2][0] << ' ' << C[2][1] << ' ' << C[2][2] << endl;
  cout << tdiff_laderman << endl;
  cout << endl;

  double speedup = (tdiff_simple-tdiff_laderman)/tdiff_laderman;
  cout << "Approximate speedup: " << speedup << endl;

  return 0;
}
Нет, но я дам им шанс (я новичок в оптимизации), вы еще что-нибудь предлагаете? Hooked
Играя с различными вещами, такими как параметры строки cmd (просто прочитайте man-страницу и попробуйте. Также у gcc есть атрибут для установки флагов оптимизации, так что вы можете взломать что-нибудь, что будет проверять несколько таких флагов в одном двоичном файле) и __restrict gcc расширение наgcc.godbolt.org часто дает хорошее представление о том, что происходит.
@PlasmaHH Я отредактировал ответ, чтобы отразить время с развертыванием, без реальных изменений. Я подозреваю, что какие бы флаги я ни реализовывал, простая версия может оказаться быстрее, хотя я немного поиграюсь с этим. Hooked
Вы пытались использовать больше опций оптимизации, таких как -funroll-loops или sse?
2

Я ожидаю, что основной проблемой производительности будет задержка памяти.double[9] обычно составляет 72 байта. Это уже нетривиальная сумма, и вы используете три из них.

72 байта слишком велики, чтобы поместиться в кэш L1? У меня сложилось впечатление, что вы можете помещать килобайты данных в кеш, и это очень быстро. Я мало знаю о сроках различных операций в процессоре; я что-то здесь упускаю?
@Dan: это будет соответствовать легко, но займет несколько строк кэша. Основное беспокойство вызывает то, что происходит, когда матрица не находится в кэше. В этом случае вы, вероятно, загрузите 2 * 64 байта из основной памяти (поскольку кэши читают целые строки).
2

Хотя в вопросе упоминался C ++, я реализовал матричное умножение 3x3 C = A * B в C # (.NET 4.5) и провел некоторые базовые тесты синхронизации на моей 64-битной машине Windows 7 с оптимизацией. 10 000 000 умножений заняло около

  1. 0.556 seconds with a naive implementation and
  2. 0.874 seconds with the laderman code from the other answer.

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

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

    public static Matrix3D operator *(Matrix3D a, Matrix3D b)
    {
        double c11 = a.M11 * b.M11 + a.M12 * b.M21 + a.M13 * b.M31;
        double c12 = a.M11 * b.M12 + a.M12 * b.M22 + a.M13 * b.M32;
        double c13 = a.M11 * b.M13 + a.M12 * b.M23 + a.M13 * b.M33;
        double c21 = a.M21 * b.M11 + a.M22 * b.M21 + a.M23 * b.M31;
        double c22 = a.M21 * b.M12 + a.M22 * b.M22 + a.M23 * b.M32;
        double c23 = a.M21 * b.M13 + a.M22 * b.M23 + a.M23 * b.M33;
        double c31 = a.M31 * b.M11 + a.M32 * b.M21 + a.M33 * b.M31;
        double c32 = a.M31 * b.M12 + a.M32 * b.M22 + a.M33 * b.M32;
        double c33 = a.M31 * b.M13 + a.M32 * b.M23 + a.M33 * b.M33;
        return new Matrix3D(
            c11, c12, c13,
            c21, c22, c23,
            c31, c32, c33);
    }

где Matrix3D является неизменяемой структурой (только для чтения двойных полей).

Самое сложное, чтобы придуматьvalid бенчмарк, где вы измеряете ваш код, а не то, что компилятор сделал с вашим кодом (отладчик с кучей лишних вещей или оптимизирован без вашего реального кода, так как результат никогда не использовался). Я обычно пытаюсь «дотронуться» результат, такой, что компилятор не может удалить тестируемый код (например, проверить матричные элементы на равенство с 89038.8989384 и throw, если он равен). Однако, в конце концов, я даже не уверен, что компилятор взломает это сравнение.

15

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

Одновременное указание нескольких команд данных выполняет одну и ту же операцию с несколькими операндами параллельно. Так что вы можете взять

double a,b,c,d;
double w = d + a; 
double x = a + b;
double y = b + c;
double z = c + d;

и заменить его на

double256 dabc = pack256(d, a, b, c);
double256 abcd = pack256(a, b, c, d);
double256 wxyz = dabc + abcd;

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

Важным фактором является число с плавающей запятой, некоторые DSP могут значительно быстрее работать с целыми числами. Графические процессоры, как правило, отлично работают с плавающей запятой, хотя некоторые из них в 2 раза быстрее с одинарной точностью. В случае 3x3 этой проблемы может поместиться в один поток CUDA, так что вы можете одновременно выполнить 256 этих вычислений.

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

Похожие вопросы