La multiplicación de matrices 3x3 de Laderman con solo 23 multiplicaciones, ¿vale la pena?

Question

Tome el producto de dos matrices de 3x3 A*B=C. Naively esto requiere 27 multiplicaciones usando el standard algorithm. Si uno fuera inteligente, podría hacer esto usando solo 23 multiplicaciones, a result found in 1973 by Laderman. La técnica implica salvar pasos intermedios y combinarlos de la manera correcta.

Ahora vamos a corregir un idioma y un tipo, por ejemplo C++ con elementos de double. Si el algoritmo de Laderman estaba codificado de forma rígida versus el doble bucle simple, ¿podríamos esperar que el rendimiento de un compilador moderno supere las diferencias de los algoritmos?

Notas sobre esta pregunta: Esta es una programación sitio, y se hace la pregunta en el contexto de la mejor práctica para un bucle interno de tiempo crítico; optimización prematura, esto no es así. Los consejos sobre la implementación son bienvenidos como comentarios.

Answer 1

pruebas de tiempo:

que corrían el tiempo pone a prueba a mí mismo y los resultados me sorprendieron (de ahí hice la pregunta en el primer lugar). En resumen, bajo una compilación estándar, el laderman es ~ 225% más rápido, pero con el indicador de optimización -03 es 50% más lento! Tuve que agregar un elemento aleatorio en la matriz cada vez durante la bandera -O3 o el compilador completamente optimizado la multiplicación simple, tomando un tiempo de cero dentro de la precisión del reloj. Dado que el algoritmo laderman fue difícil de verificar/revisar, publicaré el código completo a continuación para la posteridad.

Especificaciones: Ubuntu 12.04, Dell Prevision T1600, gcc.por ciento de diferencia en los tiempos:

• 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]

código de evaluación comparativa, junto con la aplicación Laderman:

#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; 
} 

Answer 2

La clave está en dominar el conjunto de instrucciones en su plataforma. Depende de tu plataforma. Hay varias técnicas, y cuando tiende a necesitar el máximo rendimiento posible, su compilador vendrá con herramientas de creación de perfiles, algunas de las cuales incorporan sugerencias de optimización. Para obtener las mejores operaciones de grano, mire la salida del ensamblador y vea si hay alguna mejora en ese nivel también.

Instrucción simultánea Los comandos de datos múltiples realizan la misma operación en varios operandos en paralelo. De modo que usted puede tomar

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

y reemplazarlo con

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

de manera que cuando los valores se cargan en registros, se cargan en un registro único de ancho de 256 bits por alguna plataforma de ficción con 256 -bit registros de ancho.

El punto flotante es una consideración importante, algunos DSP pueden ser mucho más rápidos en números enteros. Las GPU tienden a ser excelentes en coma flotante, aunque algunas son 2 veces más rápidas con una sola precisión. El caso 3x3 de este problema podría caber en un solo hilo CUDA, por lo que podría transmitir 256 de estos cálculos simultáneamente.

Elija su plataforma, lea la documentación, implemente varios métodos diferentes y perfílleselos.

Answer 3

Espero que el principal problema de rendimiento sea la latencia de la memoria. A double[9] es típicamente de 72 bytes. Eso ya es una cantidad no trivial, y estás usando tres de ellos.

Answer 4

Aunque la cuestión mencionada C++, que implementa la matriz 3x3 multiplicación C = A * B en C# (.NET 4.5) y ejecutó algunas pruebas básicas de tiempo en mi máquina Windows 7 de 64 bits con optimizaciones. 10.000.000 multiplicaciones tomó cerca de

1. 0,556 segundos con una aplicación ingenua y
2. 0,874 segundos con el código Laderman de la otra respuesta.

Curiosamente, el código de laderman era más lento que el ingenuo. No investigué con un generador de perfiles, pero creo que las asignaciones adicionales son más costosas que algunas multiplicaciones adicionales.

Parece que los compiladores actuales son lo suficientemente inteligentes como para hacer esas optimizaciones para nosotros, lo cual es bueno. Aquí está el código ingenua he usado, por su interés:

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); 
    } 

donde Matrix3D es una estructura inmutable (de sólo lectura campos dobles).

Lo difícil es llegar a una válida referencia , donde se mide el código y no es así, lo que el compilador hizo con su código (depurador con un montón de material extra, o optimizado sin su código real ya que el resultado nunca fue usado). Por lo general, trato de "tocar" el resultado, de manera que el compilador no puede eliminar el código bajo prueba (por ejemplo, verifique los elementos de la matriz para la igualdad con 89038.8989384 y throw, si es igual). Sin embargo, al final ni siquiera estoy seguro de si el compilador hackeó esta comparación :)