行列の乗算の高速化

よく知られた方法に行列の片方を転置して掛け合わせるという方法があります。
通常のように

A_ikB_kj

としてkで回すと、Bのほうは飛び飛びのアドレスを見ることになり、キャッシュのヒットミスが起こりやすくなります。わざわざ転置してから、

A_ikB_jk

とするとそれがなくなります。

#include <iostream>
#include <algorithm>
#include <windows.h>

#pragma comment(lib, "winmm.lib")

const int   N = 512;

void mat_mul_trans(const int A[][N], const int B[][N], int C[][N]) {
    int (*tB)[N] = new int[N][N];
    for(int i = 0; i < N; i++) {
        for(int j = 0; j < N; j++)
            tB[j][i] = B[i][j];
    }
    
    for(int i = 0; i < N; i++) {
        for(int j = 0; j < N; j++) {
            C[i][j] = 0;
            for(int k = 0; k < N; k++)
                C[i][j] += A[i][k] * tB[j][k];
        }
    }
    delete [] tB;
}

void mat_mul(const int A[][N], const int B[][N], int C[][N]) {
    for(int i = 0; i < N; i++) {
        for(int j = 0; j < N; j++) {
            C[i][j] = 0;
            for(int k = 0; k < N; k++)
                C[i][j] += A[i][k] * B[k][j];
        }
    }
}

int main() {
    int (*A)[N] = new int[N][N];
    int (*B)[N] = new int[N][N];
    for(int i = 0; i < N; i++)
        std::fill(A[i], A[i] + N, 1);
    
    const int   t0 = timeGetTime();
    mat_mul_trans(A, A, B);
    const int   t1 = timeGetTime();
    std::cout << B[0][0] << std::endl;
    std::cout << (t1 - t0) * 1e-3 << "s" << std::endl;
    mat_mul(A, A, B);
    const int   t2 = timeGetTime();
    std::cout << B[0][0] << std::endl;
    std::cout << (t2 - t1) * 1e-3 << "s" << std::endl;
    
    delete [] A;
    delete [] B;
}

転置あり 0.359s 転置なし 7.375sとなりました。ふつうはこんなに差がつかないと思いますが。
では、Pythonではどうでしょう。

import time

N = 512

def mat_mul_trans(A, B, C):
    tB = [ [ B[j][i] for j in xrange(N) ] for i in xrange(N) ]
    for i in xrange(N):
        for j in xrange(N):
            C[i][j] = sum(A[i][k] * tB[j][k] for k in xrange(N))

def mat_mul(A, B, C):
    for i in xrange(N):
        for j in xrange(N):
            C[i][j] = sum(A[i][k] * B[k][j] for k in xrange(N))

A = [ [ 1 ] * N for k in xrange(N) ]
B = [ [ 0 ] * N for k in xrange(N) ]
t0 = time.clock()
mat_mul_trans(A, A, B)
t1 = time.clock()
print B[0][0]
print t1 - t0
mat_mul(A, A, B)
t2 = time.clock()
print B[0][0]
print t2 - t1

転置あり 75.6s 転置なし 78.9sとなりました。予想通りです。Pythonはインデックス参照が非常に遅くて、C++より2桁遅いという印象があります。こんなところは効いてこないです。