Python Matrix Multiplication and Caching
<p>I am investigating caching behavior in different languages. I create two matrices in python using lists (yes, I know it's a linked list but bear with me here). I then multiply these matrices together in three ways:</p>

<pre><code>def baseline_matrix_multiply(a, b, n):
    '''
    baseline multiply
    '''
    c = zero_matrix(n)
    for i in range(n):
        for j in range(n):
            for k in range(n):
                c[i][j] += a[i][k] * b[k][j]
    return c


def baseline_matrix_multiply_flipjk(a, b, n):
    '''
    same as baseline but switch j and k loops
    '''
    c = zero_matrix(n)
    for i in range(n):
        for k in range(n):
            for j in range(n):
                c[i][j] += a[i][k] * b[k][j]
    return c


def fast_matrix_multiply_blocking(a, b, n):
    '''
    use blocking
    '''
    c = zero_matrix(n)

    block = 25;
    en = int(block * n/block)

    for kk in range(0, en, block):
        for jj in range(0, en, block):
            for i in range(n):
                for j in range(jj, jj + block):
                    sum = c[i][j]
                    for k in range(kk, kk + block):
                        sum += a[i][k] * b[k][j]
                    c[i][j] = sum
    return c
</code></pre>

<p>My timings are as follows:</p>

<pre><code>Baseline:
3.440004294627216

Flip j and k:
3.4685347505603144
100.83% of baseline

Blocking:
2.729924394035205
79.36% of baseline
</code></pre>

<p>Some things to note:</p>

<ol>
<li><p>I am familiar with CPU caching behavior. To see my experiment in C see <a href="http://codereview.stackexchange.com/questions/143061/cache-conscious-simd-matrix-multiply-of-unsigned-integers">here</a> though I havent gotten any reviews for it.</p></li>
<li><p>I've done this in Javascript and C# and the flip-j-k function provides significant performance gains using arrays (JS run on chrome browser)</p></li>
<li><p>Python implementation is Python 3.5 by way of Anaconda</p></li>
<li><p>Please dont tell me about numpy. My experiment is not about absolute performance but rather understanding caching behavior.</p></li>
</ol>

<p><strong>Question</strong>:
Anyone know what is going on here? Why does flip-j-k not provide speedup? Is it because it's a linked list? But then why does blocking provide a non-marginal improvement in performance?</p>

<p>The latter version (block multiplication) is only faster by 30% because you save 30% of index accessing by using a local variable in the inner loop!</p>

<p>(and FYI this is not C++: <code>list</code> type is like C++ <code>vector</code> otherwise people would waste their time trying to access elements by index. So this is the fastest standard random access container available in python)</p>

<p>I just made some test program based on your code to prove my point and what I was suspecting (I had to complete your code, sorry for the big block but at least it is minimal complete &amp; verifiable for all).</p>

<pre><code>def zero_matrix(sz):
    return [[0]*sz for i in range(sz)]

def baseline_matrix_multiply(a, b, n):
    '''
    baseline multiply
    '''
    c = zero_matrix(n)
    for i in range(n):
        for j in range(n):
            for k in range(n):
                c[i][j] += a[i][k] * b[k][j]
    return c


def baseline_matrix_multiply_flipjk(a, b, n):
    '''
    same as baseline but switch j and k loops
    '''
    c = zero_matrix(n)
    for i in range(n):
        for k in range(n):
            for j in range(n):
                c[i][j] += a[i][k] * b[k][j]
    return c

def baseline_matrix_multiply_flipjk_faster(a, b, n):
    '''
    same as baseline but switch j and k loops
    '''
    c = zero_matrix(n)
    for i in range(n):
        ci = c[i]
        for k in range(n):
            bk = b[k]
            aik = a[i][k]
            for j in range(n):
                ci[j] += aik * bk[j]
    return c

def fast_matrix_multiply_blocking(a, b, n):
    '''
    use blocking
    '''
    c = zero_matrix(n)

    block = 25;
    en = int(block * n/block)

    for kk in range(0, en, block):
        for jj in range(0, en, block):
            for i in range(n):
                for j in range(jj, jj + block):
                    sum = c[i][j]
                    for k in range(kk, kk + block):
                        sum += a[i][k] * b[k][j]
                    c[i][j] = sum
    return c

def fast_matrix_multiply_blocking_faster(a, b, n):
    '''
    use blocking
    '''
    c = zero_matrix(n)

    block = 25;
    en = int(block * n/block)

    for kk in range(0, en, block):
        for jj in range(0, en, block):
            for i in range(n):
                ai = a[i]
                ci = c[i]
                for j in range(jj, jj + block):
                    s = ci[j]
                    for k in range(kk, kk + block):
                        s += ai[k] * b[k][j]
                    ci[j] = s
    return c

def init_ab(sz):
    return [list(range(sz)) for i in range(sz)],[[3]*sz for i in range(sz)]

sz=200


import time

a,b = init_ab(sz)

start_time=time.time()
baseline_matrix_multiply(a,b,sz)
print("baseline_matrix_multiply: "+str(time.time()-start_time))

a,b = init_ab(sz)

start_time=time.time()
baseline_matrix_multiply_flipjk(a,b,sz)
print("baseline_matrix_multiply_flipjk: "+str(time.time()-start_time))

a,b = init_ab(sz)


start_time=time.time()
fast_matrix_multiply_blocking(a,b,sz)
print("fast_matrix_multiply_blocking: "+str(time.time()-start_time))

a,b = init_ab(sz)

start_time=time.time()
baseline_matrix_multiply_flipjk_faster(a,b,sz)
print("**baseline_matrix_multiply_flipjk_faster**: "+str(time.time()-start_time))

a,b = init_ab(sz)

start_time=time.time()
fast_matrix_multiply_blocking_faster(a,b,sz)
print("**fast_matrix_multiply_blocking_faster**: "+str(time.time()-start_time))
</code></pre>

<p>results on my PC (last results surrounded with stars are my versions):</p>

<pre><code>baseline_matrix_multiply: 2.578160047531128
baseline_matrix_multiply_flipjk: 2.5315518379211426
fast_matrix_multiply_blocking: 1.9359750747680664
**baseline_matrix_multiply_flipjk_faster**: 1.4532990455627441
**fast_matrix_multiply_blocking_faster**: 1.7031919956207275
</code></pre>

<p>As you can see, my version of your <code>baseline_matrix_multiply_flipjk</code> (the fourth one) is faster than even the block multiply, meaning that index checking &amp; accessing dwarves the cache effect that you can experience in compiled languages &amp; using direct pointers like C or C++.</p>

<p>I just stored the values that were not changing in the inner loop (the one to optimize most) to avoid index access.</p>

<p>Note that I tried to apply the same recipe to the block multiply and I gained some time compared to your version, but is still beaten by the flipjk_faster version because of the unability to avoid index access.</p>

<p>Maybe compiling the code using Cython and drop the checks would get the result you want. But pre-computing indexes never hurts.</p>

python caching matrix