Nhân ma trận 3x3 của Laderman chỉ với 23 phép nhân, có xứng đáng không?

Question

Lấy sản phẩm của hai ma trận 3x3 A*B=C. Điều này đòi hỏi 27 phép nhân sử dụng standard algorithm. Nếu một người thông minh, bạn có thể làm điều này chỉ bằng 23 phép nhân, a result found in 1973 by Laderman. Kỹ thuật này liên quan đến việc lưu các bước trung gian và kết hợp chúng một cách đúng đắn.

Bây giờ, hãy khắc phục ngôn ngữ và loại, giả sử C++ với các thành phần của double. Nếu thuật toán Laderman được mã hóa cứng so với vòng lặp đôi đơn giản, chúng ta có thể mong đợi hiệu suất của một trình biên dịch hiện đại để vượt qua sự khác biệt của các thuật toán không?

Ghi chú về câu hỏi này: Đây là một lập trình trang web, và câu hỏi được hỏi trong bối cảnh của việc thực hành tốt nhất cho một vòng lặp bên trong thời gian rất quan trọng; sớm tối ưu hóa điều này là không. Lời khuyên về thực hiện được hoan nghênh rất nhiều như ý kiến.

Answer 1

Thời gian thử nghiệm:

Tôi chạy thời gian thử nghiệm bản thân mình và kết quả làm tôi ngạc nhiên (vì thế tại sao tôi hỏi những câu hỏi ở nơi đầu tiên). Tóm lại, dưới chuẩn biên dịch, laderman nhanh hơn ~ 225%, nhưng với cờ tối ưu hóa -03, nó là 50% chậm hơn! Tôi đã phải thêm một phần tử ngẫu nhiên vào ma trận mỗi lần trong cờ -O3 hoặc trình biên dịch được tối ưu hóa hoàn toàn phép nhân đơn giản, lấy một thời gian bằng không trong độ chính xác của đồng hồ. Vì thuật toán laderman là một nỗi đau để kiểm tra/kiểm tra lại, tôi sẽ đăng mã hoàn chỉnh dưới đây cho hậu thế.

Thông số: Ubuntu 12.04, Dell Prevision T1600, gcc.chênh lệch phần trăm trong timings:

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

Benchmarking đang cùng với thực hiệ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

Điều quan trọng là làm chủ tập lệnh trên nền tảng của bạn. Nó phụ thuộc vào nền tảng của bạn. Có một số kỹ thuật, và khi bạn có xu hướng cần hiệu suất tối đa, trình biên dịch của bạn sẽ đi kèm với các công cụ lược tả, một số công cụ có tối ưu hóa gợi ý được tích hợp sẵn. ở cấp độ đó.

Lệnh cùng lúc nhiều lệnh dữ liệu thực hiện cùng thao tác trên một số toán hạng song song. Vì vậy mà bạn có thể mang

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

và thay thế bằng

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

Vì vậy mà khi các giá trị được nạp vào thanh ghi, họ được nạp vào một thanh ghi rộng 256-bit duy nhất cho một số nền tảng hư cấu với 256 -bit đăng ký rộng.

Điểm nổi là một yếu tố quan trọng, một số DSP có thể hoạt động nhanh hơn đáng kể trên các số nguyên. GPU có xu hướng tuyệt vời trên điểm nổi, mặc dù một số có tốc độ nhanh gấp 2 lần ở độ chính xác đơn. Trường hợp 3x3 của vấn đề này có thể phù hợp với một chuỗi CUDA duy nhất, vì vậy bạn có thể phát 256 đồng thời các tính toán này cùng một lúc.

Chọn nền tảng của bạn, đọc tài liệu, triển khai một số phương pháp khác nhau và cấu hình chúng.

Answer 3

Tôi hy vọng mối quan tâm hiệu suất chính sẽ là độ trễ của bộ nhớ. A double[9] thường là 72 byte. Đó là một số tiền không đáng kể, và bạn đang sử dụng ba trong số đó.

Answer 4

Mặc dù câu hỏi nêu C++, tôi thực hiện ma trận 3x3 phép nhân C = A * B trong C# (.NET 4.5) và chạy một số kiểm tra thời gian cơ bản trên máy tính Windows 64 bit 7 của tôi với tối ưu hóa. 10.000.000 phép nhân mất khoảng

1. 0,556 giây với một thực hiện ngây thơ và
2. 0,874 giây với mã laderman từ câu trả lời khác.

Điều thú vị là mã laderman chậm hơn con đường ngây thơ. Tôi đã không điều tra với một hồ sơ, nhưng tôi đoán phân bổ thêm là tốn kém hơn so với một vài phép nhân thêm.

Dường như các trình biên dịch hiện tại đủ thông minh để thực hiện các tối ưu hóa đó cho chúng tôi, điều này là tốt. Đây là mã ngây thơ tôi đã sử dụng, vì sở thích của bạn:

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

nơi Matrix3D là cấu trúc không thay đổi (trường đôi chỉ đọc).

Điều khó khăn là đưa ra tiêu chuẩn hợp lệ điểm chuẩn, nơi bạn đo mã của mình chứ không phải, trình biên dịch đã làm với mã của bạn (trình gỡ rối với tấn công cụ bổ sung hoặc được tối ưu hóa mà không có mã thực tế của bạn không bao giờ được sử dụng). Tôi thường cố gắng "chạm vào" kết quả, sao cho trình biên dịch không thể loại bỏ mã đang được kiểm tra (ví dụ: kiểm tra các phần tử ma trận cho sự bình đẳng với 89038.8989384 và ném, nếu bằng nhau). Tuy nhiên, cuối cùng tôi thậm chí không chắc chắn nếu trình biên dịch hacks so sánh này ra khỏi con đường :)