Bỏ qua điều hướng, tới nội dung chính
Học C
Bài 10.424 phút đọc

Nhân ma trận

Sau bài này bạn sẽ làm được

  • Cài phép nhân ma trận đúng công thức
  • Kiểm tra điều kiện số cột bằng số hàng
  • Giải thích vì sao đổi thứ tự ba vòng lặp làm chương trình nhanh gấp nhiều lần
  • Đo và so sánh hai thứ tự vòng lặp

Nhân ma trận là bài toán đầu tiên trong khóa học có ba vòng lặp lồng nhau, và cũng là bài toán đầu tiên mà chỉ đổi thứ tự các vòng lặp, không đổi một phép tính nào, đã làm chương trình nhanh lên nhiều lần.

#Công thức

Tích hai ma trận
Phần tử ở hàng i cột j của tích bằng tổng các tích của hàng i ma trận thứ nhất với cột j ma trận thứ hai.
C[i][j] = tổng theo k của ( A[i][k] * B[k][j] )
Ma trậnKích thướcGhi chú
Am x nm hàng, n cột
Bn x pSố hàng của B PHẢI bằng số cột của A
Cm x pSố hàng lấy từ A, số cột lấy từ B

Ví dụ tính tay

A = | 1  2 |     B = | 5  6 |
    | 3  4 |         | 7  8 |

C[0][0] = 1*5 + 2*7 = 19
C[0][1] = 1*6 + 2*8 = 22
C[1][0] = 3*5 + 4*7 = 43
C[1][1] = 3*6 + 4*8 = 50

C = | 19  22 |
    | 43  50 |

#Cài đặt

nhan.c
/* A là m x n, B là n x p, C là m x p. Mọi ma trận đều phẳng. */
void nhan_ma_tran(const double *A, const double *B, double *C,
                  size_t m, size_t n, size_t p)
{
    for (size_t i = 0; i < m; ++i)
        for (size_t j = 0; j < p; ++j) {
            double s = 0.0;

            for (size_t k = 0; k < n; ++k)
                s += A[i * n + k] * B[k * p + j];

            C[i * p + j] = s;
        }
}

Phiên bản có kiểm tra tham số

an-toan.c
int nhan_ma_tran(const double *A, const double *B, double *C,
                 size_t m, size_t n, size_t p)
{
    if (A == NULL || B == NULL || C == NULL) return -1;
    if (m == 0 || n == 0 || p == 0)          return -1;
    if (C == A || C == B)                    return -1;   /* chồng lấn */

    for (size_t i = 0; i < m; ++i)
        for (size_t j = 0; j < p; ++j) {
            double s = 0.0;

            for (size_t k = 0; k < n; ++k)
                s += A[i * n + k] * B[k * p + j];

            C[i * p + j] = s;
        }

    return 0;
}

Kiểm tra C == A chỉ bắt được trường hợp trùng hoàn toàn. Trường hợp hai vùng chồng lấn một phần thì không có cách nào kiểm tra di động trong C. Đây là lý do các thư viện thật ghi rõ trong tài liệu rằng vùng đích không được chồng lấn vùng nguồn.

#Độ phức tạp lập phương

Với ba ma trận vuông cạnh n, thuật toán thực hiện đúng n³ phép nhân và n³ phép cộng. Đó là con số tăng rất nhanh.

nSố phép nhânThời gian ước tínhBộ nhớ cho ba ma trận double
1001 triệuvài mili giây240 KB
1 0001 tỷkhoảng 1 giây24 MB
10 0001 nghìn tỷkhoảng 20 phút2.4 GB
100 0001 triệu tỷkhoảng 2 tuần240 GB

Có thuật toán nhanh hơn không

Thuật toánĐộ phức tạpThực tế
Ba vòng lặpn³Dùng cho mọi kích thước thông thường
Strassen, 1969n^2.807Có lợi từ khoảng n bằng 1000 trở lên, nhưng kém ổn định về số học
Coppersmith và Winograd cùng các cải tiếndưới n^2.372Chỉ có ý nghĩa lý thuyết, hằng số quá lớn

Trong thực tế, thư viện tính toán số dùng ba vòng lặp nhưng chia ma trận thành các khối vừa bộ nhớ đệm, và tận dụng lệnh véc tơ cùng nhiều nhân xử lý. Một thư viện tốt nhanh hơn bản ngây thơ của bạn khoảng năm mươi tới một trăm lần, mà vẫn là cùng độ phức tạp n³.

terminal
# Với ma trận 1024 nhân 1024, cùng máy
./so-sanh-nhan
ba vong lap i j k : 3.412 s
ba vong lap i k j : 0.782 s
thu vien BLAS     : 0.031 s

#Thứ tự ba vòng lặp

Ba vòng lặp i, j, k có sáu cách sắp xếp. Cả sáu đều cho kết quả toán học giống hệt nhau, nhưng tốc độ chênh nhau tới năm lần trên ma trận lớn.

hai-thu-tu.c
/* Thứ tự i j k: vòng trong chạy theo k */
for (size_t i = 0; i < n; ++i)
    for (size_t j = 0; j < n; ++j) {
        double s = 0.0;

        for (size_t k = 0; k < n; ++k)
            s += A[i * n + k] * B[k * n + j];   /* B nhảy n phần tử mỗi lượt */

        C[i * n + j] = s;
    }

/* Thứ tự i k j: vòng trong chạy theo j */
for (size_t k = 0; k < n * n; ++k) C[k] = 0.0;

for (size_t i = 0; i < n; ++i)
    for (size_t k = 0; k < n; ++k) {
        double a = A[i * n + k];                /* đọc một lần, giữ trong thanh ghi */

        for (size_t j = 0; j < n; ++j)
            C[i * n + j] += a * B[k * n + j];   /* B và C đều đi liền mạch */
    }
Thứ tựA truy cập thế nàoB truy cập thế nàoC truy cập thế nào
i j kLiền mạch theo hàngNhảy quãng theo cộtGhi một lần mỗi ô
i k jMột phần tử, giữ trong thanh ghiLiền mạch theo hàngLiền mạch theo hàng
j k iNhảy quãngMột phần tửNhảy quãng

Cách khác: chuyển vị B trước

/* Chuyển vị B một lần, tốn n² phép sao chép,
   sau đó cột của B trở thành hàng của Bt nên đọc liền mạch. */
chuyen_vi(B, Bt, n, n);

for (size_t i = 0; i < n; ++i)
    for (size_t j = 0; j < n; ++j) {
        double s = 0.0;

        for (size_t k = 0; k < n; ++k)
            s += A[i * n + k] * Bt[j * n + k];   /* cả hai đều liền mạch */

        C[i * n + j] = s;
    }

Chi phí chuyển vị là n², nhỏ hơn hẳn n³ của phép nhân, nên nó gần như miễn phí. Bài 10.5 sẽ cài phép chuyển vị.

#Ứng dụng: lũy thừa ma trận

Một ứng dụng đẹp của phép nhân ma trận là tính số Fibonacci thứ n trong thời gian log n, thay vì n như vòng lặp ở Bài 8.10.

| F(n+1)  F(n)   |     | 1  1 |^n
| F(n)    F(n-1) |  =  | 1  0 |
fib-ma-tran.c
typedef struct { long long m[2][2]; } M2;

static M2 nhan2(M2 a, M2 b)
{
    M2 c = { { { 0, 0 }, { 0, 0 } } };

    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            for (int k = 0; k < 2; ++k)
                c.m[i][j] += a.m[i][k] * b.m[k][j];

    return c;
}

/* Lũy thừa nhanh, cùng ý tưởng với Bài 7.2 */
static M2 luy_thua(M2 a, unsigned n)
{
    M2 kq = { { { 1, 0 }, { 0, 1 } } };      /* ma trận đơn vị */

    while (n > 0) {
        if (n & 1u) kq = nhan2(kq, a);

        a = nhan2(a, a);
        n >>= 1;
    }

    return kq;
}

long long fib(unsigned n)
{
    if (n == 0) return 0;

    M2 a = { { { 1, 1 }, { 1, 0 } } };

    return luy_thua(a, n - 1).m[0][0];
}
terminal
./fib-ma-tran 90
fib(90) = 2880067194370816120
so phep nhan ma tran: 11

Tự làm thử

  1. Cài phép nhân ma trận, kiểm tra bằng ví dụ tính tay trong bài, rồi thử với ma trận đơn vị.
  2. Xác nhận bằng thực nghiệm rằng A * B khác B * A với hai ma trận vuông ngẫu nhiên.
  3. Cài cả sáu thứ tự vòng lặp, đo thời gian với n bằng 512 và 1024, lập bảng xếp hạng.
  4. Cài phiên bản chuyển vị B trước rồi nhân, so sánh thời gian với hai thứ tự nhanh nhất.
  5. Đo thời gian với n bằng 128, 256, 512, 1024 và kiểm tra xem tỷ lệ có đúng gấp tám lần mỗi khi n gấp đôi không.
  6. Cài lũy thừa ma trận để tính Fibonacci, so sánh số phép tính với vòng lặp thường khi n bằng một tỷ.
  7. Dùng lũy thừa ma trận kề để đếm số đường đi độ dài 5 giữa hai đỉnh trong một đồ thị nhỏ.

Trình chấm điểm tự động sẽ được bổ sung ở giai đoạn sau. Hiện tại bạn tự chạy thử trên máy.

Tóm tắt

  • Tích hai ma trận đòi hỏi số cột của A bằng số hàng của B, và phép nhân không có tính giao hoán.
  • Ba vòng lặp lồng nhau cho độ phức tạp lập phương: gấp đôi kích thước là gấp tám lần thời gian.
  • Tích lũy vào biến cục bộ rồi ghi một lần, đừng cộng dồn thẳng vào ô của ma trận kết quả.
  • Ma trận đích không được trùng với ma trận nguồn, khác với phép cộng theo từng phần tử.
  • Thứ tự i k j nhanh hơn i j k nhiều lần vì cả ba mảng đều được duyệt liền mạch, dù số phép tính không đổi.
  • Lũy thừa ma trận tính số hạng thứ n của mọi dãy truy hồi tuyến tính trong log n phép nhân.