卷積圖層矩陣計算方法 — Convolution using Matrix to calculation

本文介紹如何比較有效率的計算 Convolution ( 卷積 )方法，不使用雙迴圈，改用矩陣運算

只要你懂矩陣，矩陣就會幫助你

前言 — 卷積( Convolution ) 介紹

卷積(Convolution) 如果有聽過深度學習( Deep Learning )的人都略有所知
其概念在影像處理上是非常有幫助且行之有年，不只適用於 Deep Learning。

本文需要有矩陣運算與 numpy 相關背景知識，著重在如何用比較有效率的計算方式來計算卷積影像 ( 我們這邊為了方便講解，只說明長寬相同的 image / kernel 以及不做 padding, stride 另外處理 )

卷積的做法是透過類似遮罩 ( Filter/Kernel) 和圖片局部運算並且掃過所有像素得出新的圖片，然後剛好這種做法可以產生特徵值 ( 對 Machine Learning 很有效)，所以就對 Deep Learning 有很大的幫助

套用不同的 Kernel / Feature 效果

如上圖，左上是原圖(灰階) 然後套用不同的效果後得出不同的圖片
這些圖片各自的數值再餵給 Machine Learning 讓他自己去學習與分析就可以更準確的提升辨識率，關於 Kernel 如何產生與更好應用這部分因為礙於篇幅，超出本文範圍，待日後有機會再說明。 Kernel 有些資料可以參考 wiki
https://en.wikipedia.org/wiki/Kernel_(image_processing)

傳統算法 — 雙雙迴圈

最傳統的 Convolution 算法就是兩個迴圈，從左到右從上到下，然後一路與 kernel 相加，我們先來看看傳統演算法，我們這邊先用 1 至 25 的 5x5 矩陣當作“影像”，然後 3x3 皆為 1 的矩陣為 kernel

img = np.arange(1, 26, 1).reshape(5, 5)
kernel = np.ones((3,3))

左 img, 右 kernel

很直白的程式碼

def convolution(img, kernel):
    img_size = img.shape[0]
    kernel_size = kernel.shape[0]
    new_img_size = img_size - kernel_size + 1
    # 準備新圖片
    new_img = np.zeros((new_img_size, new_img_size))
    
    # 經典雙迴圈
    for i in range(new_img_size):
        for j in range(new_img_size):
            value = 0
            # 經典雙雙迴圈 加大不加價(誤)
            for ki in range(kernel_size):
                for kj in range(kernel_size):
                    value += img[i+ki][j+kj] * kernel[ki][kj]
            new_img[i][j] = value
    return new_img

來驗證一下數值

[[ 63.  72.  81.]
 [108. 117. 126.]
 [153. 162. 171.]]

簡單手算一下

print((1 + 2 + 3) + (6 + 7 + 8) + (11 + 12 + 13)) # 位置 0x0  = 63
print((2 + 3 + 4) + (7 + 8 + 9) + (12 + 13 + 14)) # 位置 0x1 = 72
...
print((13 + 14 + 15) + (18 + 19 + 20) + (23 + 24 + 25)) # 位置 2x2 = 171

看來沒錯，接下來就是改善演算法了

傳統算法 — 還是雙迴圈

只要你懂矩陣，矩陣就會幫助你

矩陣運算對於 CPU / GPU 其實都有特別的加速指令，所以能夠將這種大量運算變成矩陣的話，對於效能提昇會非常有幫助
觀察一下程式碼可以發現內迴圈其實就是

a1 * b1 + a2 * b2 + a3 * b3 …

這種行為，剛好就是矩陣的 dot 運算 ( np / pytorch 可以使用 @ 表示 )
而 a 我們假設為 kernel, b 假設為 image
我們只要用 numpy 的 flatten 就可以得到 row vector
但是 column vector 依然需要靠迴圈取得 ( 後面會說更有效率的做法 )

程式碼如下

def convolution(img, kernel):
    img_size = img.shape[0]
    kernel_size = kernel.shape[0]
    new_img_size = img_size - kernel_size + 1
    # 準備新圖片
    new_img = np.zeros((new_img_size, new_img_size))
    # -> 將 kernel 打平變成 row vector
    row_vector = kernel.flatten() 
    
    # 經典雙迴圈
    for i in range(new_img_size):
        for j in range(new_img_size):
            value = 0
            column_vector = []
            # 經典雙雙迴圈 但是快了一點
            for ki in range(kernel_size):
                for kj in range(kernel_size):
                    # 從 img 中取出數值放入 column_vector
                    column_vector.append(img[i+ki][j+kj])
            # 這邊就是使用矩陣運算
            new_img[i][j] = row_vector @ column_vector
    return new_img

現在剩下的問題就是如何取得大量的 column_vector ( 從 img 中轉換出來 )

另外取得 column_vector 其實有辦法讓內回圈只要一次 ( 用 array np.concatenate ) ，但是不好閱讀所以上文才用雙迴圈示意
而且經典雙迴圈問題依然存在

# 經典雙迴圈
for i in range(new_img_size):
    for j in range(new_img_size):
        value = 0
        column_vector = np.array([], dtype=img.dtype)
        # 一個迴圈, 用 np.concatenate
        for ki in range(kernel_size):
            column_vector = np.concatenate([column_vector, img[i+ki, j:j+kernel_size]])

Convolution — 矩陣運算法

觀察一下 column_vector 爬過的軌跡

依此類推爬完整個 image

可以發現是從左到右，從上到下的去跑
然後再觀察一下裡面的內容
[
[ 1,2,3...] <- 小矩陣,
[ 2,3,4...] <- 小矩陣,
[ 3,4,5...] <- 小矩陣
…
] <- 大矩陣
會發現其實是兩個矩陣在跑

Numpy 有個工具
https://numpy.org/doc/stable/reference/generated/numpy.lib.stride_tricks.as_strided.html

在計算機概論中，陣列的位置存取其實就是對應記憶體位置中的數值
所以我們只要改變對於記憶體位置的看法，就可以讓 numpy 按照我們的意思去抓取數值

其中 numpy 對應記憶體使用 “stride” (跨步) 方式，也和 convolution 中的 stride 概念一樣，下一步要跳多少

理論說來複雜，我們先來看範例

a = np.arange(1, 10, 1)  # 這是一維陣列
stride = a.strides[0] # 我們要改變 numpy 對於資料的看法, 要先確定"跨步"基本單位
# 我們改變 numpy 對於 array 的看法
# 告訴他形狀會是 (3,3), 
# stride(跨步) 規則是 (3單位, 1單位)
b = np.lib.stride_tricks.as_strided(a, shape=(3,3), strides=(stride * 3, stride))
print('a:', a)
print('b:', b)

可以看到

a: [1 2 3 4 5 6 7 8 9]
b: [[1 2 3]
 [4 5 6]
 [7 8 9]]

我們改變了對於 a 的看法確實得到 b 矩陣
再來玩看看 stride(跨步) 規則是 ( 1單位, 3單位)
理論上:
第一行會被選中 [1], 2, 3, [4], 5, 6, [7], 8, 9
第二行會被選中 1, [2], 3, 4, [5], 6, 7, [8], 9
第三行會被選中 1, 2, [3], 4, 5, [6], 7, 8, [9]

np.lib.stride_tricks.as_strided(a, shape=(3,3), strides=(stride, stride*3))
# 輸出
 [[1 4 7]
 [2 5 8]
 [3 6 9]]

果然沒錯

有了以上概念我們就可以透過這個方法來取得 column_vector 了
(這個是雙雙矩陣在內概念… 各位可以自己拿紙筆算一下 )

測試一下概念
我們假設 kernel size 為 2, 測試的 image size 為 5

img_size = 5
kernel_size = 2
new_img_size = img_size - kernel_size + 1
test_img = np.arange(0, img_size **2, 1).reshape(img_size, img_size) + 1
def test_stride():
    b = np.lib.stride_tricks.as_strided(test_img,
        shape=(
            new_img_size,
            new_img_size,
            kernel_size,
            kernel_size,
        ),
        strides=(
        stride * img_size,
        stride,
        stride * img_size,
        stride,
    ))
return b.reshape(-1, kernel_size * kernel_size)

我們視覺化一下資料

左為原矩陣, 右為轉化過的

看起來正確，以上就是我們要與 kernel 相乘的矩陣
然後我們只要 transpose 矩陣就可以了，包裝成 function 的程式碼如下範例

隨便套點 kernel 試看看效果

結語

numpy 的 stride 這部分其實不好透過文字說明，但是讀者應該可以透過調整參數去玩看看，當然實務上 Tensorflow / Pytorch / … 都已經有打包好的方法可以使用，但是有辦法透過自己實作並且了解原理也是一種加強實力的方法

只要你懂矩陣，矩陣就會幫助你

祝大家新年快樂 : )

最後這是本文用來顯示矩陣數值(圖表)的方法

def show_number(img):
    shape = img.shape
    for i in range(shape[0]):
        for j in range(shape[1]):
            c = img[i, j]
            plt.text(j, i, str(c), va='center', ha='center')
            
def show_matrix(img):
    plt.imshow(img, cmap='gray')
    show_number(img)
show_matrix(img)
plt.clim(-10, 90)