Python
機械学習
DeepLearning
可視化
Keras

最新手法!「Grad-CAM++」のレビューと実装


はじめに

現象数理学科4年Dasonです。Qiita8記事目。

この記事は「明治大学 Advent Calendar 2018」第8日目の記事です。

本当はこの日に書く気はなかったのですが、PayPayでNintendo Switchを買ってテンションが上がった...兼ねてから実装を試みていたGrad-CAM++の実装ができたので、勢いで書いちゃうことにしました。

簡単になんでこの内容を選んだのか、どういう仕組みなのかを説明した後、実装&実験といった流れで話を進めて行きます。

では、いきましょう。


背景

私自身、普段はDeepLearningによる医療画像の疾病分類など4種類の研究をやってますが、それらに共通するのは画像を機械学習によって学習して、分類精度や回帰結果を見て一喜一憂する、といったことになります。しかし、機械学習はなぜそのように判断したのかの根拠が示されない、いわばブラックボックス となっていて、研究者を悩ませる、といった状況が続いていました。そこに登場したのが、CAMやGrad-CAMといった、どこに注目して判断したのかを可視化できるツールです。(それらについては、また後日)

あとで出てきますが、こんな感じで可視化できます。

both.png

cat_dog.pngFri Dec  7 19:31:23 2018_GCAM++_%s.jpg

近年Grad-CAMは世界中に普及して、よく使われるツールに仲間入りしたのですが、なぜかGrad-CAMの進化版とされる「Grad-CAM++」はどこにも記事がない!「実装してみた」もほとんどない!困った!

といったわけで、僕が今回Grad-CAM++を実装して使ってみたをやってみるわけです。


仕組み

まずは、Grad-CAM++がどのような仕組みで動いているのか、ちょっと数式を追ってみます。


概要

さて、Grad-CAM++は何をしているのか、確認していきましょう。

なお、元の論文はこちらにあるので、適宜参照してください。また、論文中に登場しますが、TensorFlowで書かれた公式Githubはこちらです。

そもそも、CAMやGrad-CAMがどのようにして可視化していたか、簡単に振り返ります。まず「Global Average Pooling layer (GAP層)」を持つCNNに対して、機械学習が入力画像が「あるクラス$c$である」と判断した、その最後の分類スコア$Y^c$を、GAP層の特徴量マップ$A^k_{ij}$を使って

Y^c = \sum_kw^{c}_k\sum_i\sum_jA_{ij}^k \tag{1}\\

と定義しましょう。これは、$A_{ij}^k$:$k$番目のチャンネルでの位置$(i,j)$での活性であることに着目すれば、$\sum_i\sum_j A_{ij}^k$は、その全体で割ってあげると「$k$番目のチャンネルでの活性の平均」となることから、ここがGAP層そのものに相当します。さらに、$Y^c$が加重平均をとっていると捉えることができるので、$w^c_k$は画像がクラス$c$であると返すための$\sum_i\sum_jA_{ij}^k$の重要度であると解釈できますね。さて、有限和の順序交換を施してあげて

Y^c = \sum_i\sum_j(\sum_kw^{c}_kA_{ij}^k)

としてあげましょう。ここで、

L^c_{ij} = \sum_{k}w_k^cA_{ij}^k \tag{2}

としてあげると

Y^c = \sum_i\sum_jL_{ij}^c

とできます。この$L^c_{ij}$が画像中のピクセル位置($i$,$j$)がその画像がクラス$c$であると判断された時の、重要度であると定義していました。(つまりは、この$L^c_{ij}$が求まれば、画像中の位置$(i,j)$での重要度のヒートマップが描けるということになりますね。)あとは$A^k_{ij}$はGAP層を見ればわかるので、$w_k^c$をどう定義するか、といったところに話が遷移していきます。

さて、CAMでは$w^c_k$を最後のconv層からfc層への重みとして定義、Grad-CAMでは

w_k^c = \dfrac{1}{Z}\sum_i\sum_j\dfrac{\partial Y^c}{\partial A_{ij}^k} \tag

{3}

として定義していました(但し、Zはactivation mapのピクセル数なので、Zで割るということはつまり平均とってるんですね)。今回の主役、Grad-CAM++では

w_k^c = \sum_i\sum_j\alpha_{ij}^{kc}relu\left(\dfrac{\partial Y^c}{\partial A_{ij}^k}\right) \tag{4}

と定義します。単純にピクセルの平均を見ていたGrad-CAMと違って、何やら不思議な重み$\alpha$が登場しています。ここが今回の要点です。

ちなみに、忙しい人は以下を読み飛ばせるように$\alpha$だけ記しておきます。

\alpha_{ij}^{kc} = \dfrac{\dfrac{\partial^2Y^c}{(\partial A_{ij}^k)^2}}{2\dfrac{\partial^2Y^c}{(\partial A_{ij}^k)^2}+\sum_a\sum_b A_{ab}^k\left\{\dfrac{\partial^3Y^c}{(\partial A_{ij}^k)^3}\right\}}


Grad-CAM++

さて、あるクラスcに対するスコア$Y^c$は、上の(3)式を用いて次のように書き換えられます。

Y^c = \sum_k\left\{\sum_a\sum_b\alpha^{kc}_{ab}relu\left(\dfrac{\partial Y^c}{\partial A_{ab}^k}\right)\right\}\sum_i\sum_jA_{ij}^k  \tag{5}

ここで、($i,j$)と($a,b$)は同じactivation map上の点を表すが、混乱が生じないように別の文字を使っていることに注意しておきます。さぁ、このReLU関数は

ReLU(x) = max(0,x) \tag{6}

で表される活性化関数です。ここで、この論文では「勾配を逆流させるための閾値としてReLU関数は機能するので、一般性を失うことなくReLU関数を無視してもよい」としています。その結果$\frac{\partial Y^c}{\partial A_{ij}^k}$は次のようにできます。

\dfrac{\partial Y^c}{\partial A_{ij}^k}=\sum_a\sum_b\alpha_{ab}^{kc}\dfrac{\partial Y^c}{\partial A_{ab}^k}+\sum_a\sum_bA_{ab}^{k}\left\{\alpha_{ij}^{kc}\dfrac{\partial^2Y^c}{(\partial A_{ij}^k)^2}\right\}\tag{7}

引き続き2階微分もしましょう。

\dfrac{\partial^2Y^c}{(\partial A_{ij}^k)^2}=2\alpha_{ij}^{kc}\dfrac{\partial^2Y^c}{(\partial A_{ij}^k)^2}+\sum_a\sum_bA_{ab}^k\left\{\alpha_{ij}^{kc}\dfrac{\partial^3Y^c}{(\partial A_{ij}^k)^3}\right\}\tag{8}

さて、上の2つの式から$\alpha$を表現できて、

\alpha_{ij}^{kc}=\dfrac{\dfrac{\partial^2Y^c}{(\partial A_{ij}^k)^2}}{2\dfrac{\partial^2Y^c}{(\partial A_{ij}^k)^2}+\sum_a\sum_bA_{ab}^k\left\{\dfrac{\partial^3Y^c}{(\partial A_{ij}^k)^3}\right\}}\tag{9}

以上から、Grad-CAM++で使う重み$w$が求まりました。

w_k^c=\sum_i\sum_j\dfrac{\dfrac{\partial^2Y^c}{(\partial A_{ij}^k)^2}}{2\dfrac{\partial^2Y^c}{(\partial A_{ij}^k)^2}+\sum_a\sum_bA_{ab}^k\left\{\dfrac{\partial^3Y^c}{(\partial A_{ij}^k)^3}\right\}}relu \left(\dfrac{\partial Y^c}{\partial A_{ij}^k}\right)\tag{10}

さて、ここで1つ注目してみましょう。

もし、$\alpha_{ij}^{kc}=\frac{1}{Z}$だとどうでしょうか。

w_k^c=\sum_i\sum_j\dfrac{1}{Z}relu \left(\dfrac{\partial Y^c}{\partial A_{ij}^k}\right)\tag{11}

こうなりますね。ここで、reluが先ほどの閾値の機能である、という捉え方をすると、ざっくり言ってGrad-CAMにも見えます。そうです。Grad-CAM++は、$\frac{1}{Z}$をかけることによって、ピクセル平均を出していたGrad-CAMの平均について2階微分、3階微分で重み付けをした、いわばGrad-CAMの一般化版と言えるものだったんですね。

最後に顕著量マップを作っておきましょう。

次の節の冒頭に出て来ますが、最後の層がReLU関数or線形関数となっていると、計算が楽になるのがこのGrad-CAM++の特徴だったりします。もし、ReLU関数が活性化関数となっている最終層をもつネットワークならば顕著量マップは

L_{ij}^c = relu(\sum_k w_{k}^c A_{ij}^k) \tag{12}

となります。


高階微分に立ち向かう

理論的になんかすごそうなのは(多分)分かりました。ネックなのは、2階微分や3階微分をどうするかです。ここで、論文の著者さんは次のように考えました。


もし、最後から2番目のLayerのスコアを指数関数を通して得られて、最後の層が線形orReLU関数からなる層ならば高階微分の計算結果は自明じゃん。


賢いですね。賢すぎて僕には言ってる意味がよくわかりません。実際に計算してみましょう。

$S^c$を最後から2番目の層のクラスcに対するスコアとします。このとき、

Y^c = \exp(S^c)\tag{13}

とせよ。と言われているわけですね。なるほど。微分します。

\dfrac{\partial Y^c}{\partial A_{ij}^k} = \exp(S^c)\dfrac{\partial S^c}{\partial A_{ij}^k}\tag{14}

いやいや、微分ありますやんけ。って思う方。実はTensorFlowやPyTorchには自動微分がライブラリに実装されてますから、すぐに計算できちゃうのです。

さて、2階微分してみましょう。

\dfrac{\partial^2Y^c}{(\partial A_{ij}^k)^2} = \exp(S^c)\left[\left(\dfrac{\partial S^c}{\partial A_{ij}^k}\right)^2+\dfrac{\partial^2S^c}{(\partial A_{ij}^k)^2}\right]\tag{15}

ここで、ReLU関数$f(x)=max(0,x)$を微分してみます。

\begin{align*}

\dfrac{\partial}{\partial x}f(x) &= 1\ \tt{for}\ x>0 \tag{16}\\
&=0\ \tt{for}\ x\leq 0
\end{align*}

さらに微分して

\dfrac{\partial^2 }

{\partial x^2}f(x) = 0 \tag{17}

を得ます。ここで、仮にReLU関数でなくて、線形な関数であったとしても2階微分は0になることに注意すれば、先ほどの式は次のようになります。

\dfrac{\partial^2Y^c}{(\partial A_{ij}^k)^2} = \exp(S^c)\left(\dfrac{\partial S^c}{\partial A_{ij}^k}\right)^2 \tag{18}

全く同様にして

\dfrac{\partial^3Y^c}{(\partial A_{ij}^k)^3} = \exp(S^c)\left(\dfrac{\partial S^c}{\partial A_{ij}^k}\right)^3 \tag{19}

を得ます。先ほど言った通り1階微分については、各種機械学習のフレームワークに搭載された自動微分で対処できますので、2階微分、3階微分も簡単に計算できることがわかりました。


Grad-CAM++よ!Softmax関数はいかが?

ここまで、一見無敵感溢れて来そうですが、いま使った「ReLU関数or線形な関数」という仮定が実は結構辛い。例えば、機械学習のクラス分類のときによく使われるsoftmax関数は線形ではありません。そのときも、似たようなことをすればいけますが、複雑になってしまいます。その件について論文でも述べてあり、計算結果が載ってますので、ここでもその結果を見せるだけに留めておきます。

まず、最終的なクラス分類におけるスコア$Y^c$を

Y^c = \dfrac{\exp(S^c)}{\sum_k \exp(S^c)}\tag{20}

とします。このとき$S^k$は最後から2番目の層でクラスkと出力されるように学習したスコアです。

このとき、各階微分は以下のようになります。

\dfrac{\partial Y^c}{\partial A_{ij}^k}= Y^c\left[\dfrac{\partial S^c}{\partial A_{ij}^k}-\sum_kY^k\dfrac{\partial S^k}{\partial A_{ij}^k}\right]\tag{21}

\dfrac{\partial^2 Y^c}{(\partial A_{ij}^k)^2}=\dfrac{\partial Y^c}{\partial A_{ij}^k}\left[\dfrac{\partial S^c}{\partial A_{ij}^k}-\sum_kY^k\dfrac{\partial S^k}{\partial A_{ij}^k}\right]-Y^c\left(\sum_k\dfrac{\partial Y^k}{\partial A_{ij}^k}\dfrac{\partial S^k}{\partial A_{ij}^k}\right)\tag{22}

\dfrac{\partial^3 Y^c}{(\partial A_{ij}^k)^3}=\dfrac{\partial^2 Y^c}{(\partial A_{ij}^k)^2}\left[\dfrac{\partial S^c}{\partial A_{ij}^k}-\sum_kY^k\dfrac{\partial S^k}{\partial A_{ij}^k}\right]-2\dfrac{\partial Y^c}{\partial A_{ij}^k}\left(\sum_k\dfrac{\partial Y^k}{\partial A_{ij}^k}\dfrac{\partial S^k}{\partial A_{ij}^k}\right)-Y^c\left(\sum_k\dfrac{\partial^2 Y^k}{(\partial A_{ij}^k)^2}\dfrac{\partial S^k}{\partial A_{ij}^k}\right)\tag{23}

はい。大変重い。これを先ほどの$\alpha$の式に入れてくれ、というのです。まぁ、計算はCPUくんが頑張ってくれるので、これを実装すりゃいいわけですね。もちろん、これは先ほどのReLU関数や線形関数の場合を含んでいるという点について注意をしておきます。


実装

いよいよ実装していきます。その前に、私が実装した環境をまとめておきます。

MacOSX 10.13.6

Python3.6.5
Keras2.2.4

余談ですが、Kerasのバージョンって

import keras

print(keras.__version__)

でわかるんですね。

本題に参ります。最初に必要なライブラリをインポートしましょう。

import numpy as np

import cv2
import argparse
import keras
import sys
from keras import backend as K
from keras.preprocessing.image import array_to_img, img_to_array, load_img
from keras.applications.resnet50 import ResNet50

今回は、使用するモデル、着目するレイヤー名、画像、画像サイズが与えられたもとで、Grad-CAM++によって描かれた画像を返す関数を作ることにします。雰囲気はこんな感じ。

def Grad_Cam_plus_plus(input_model, layer_name, img, target_size):

いろいろやって
return ヒートマップ画像(array)

さて、まずは必要な前処理から行きます。次元数をnp.expand_dimsで3次元から4次元に拡張した上で正規化します。

    #モデルを指定。

model = input_model
#画像の前処理
X = np.expand_dims(x, axis=0)

X = X.astype('float32')
preprocessed_input = X / 255.0

次に、予測クラスの算出をしましょう。今回はその画像がどのクラスである確率が高いか判定した上で、そのクラスであると判断したときに、画像中のどこが決め手であったかをヒートマップによって描くことにします。

    predictions = model.predict(preprocessed_input)

class_idx = np.argmax(predictions[0])

それぞれ必要な重みをとってきます。S_cが$S^c$に、conv_outputが$A_{ij}^k$に相当します。また、Kerasの場合だと、そのバックエンドとなるTensorFlowの関数を一部使うことができて、ここではK.gradients(A,B)という関数を使用しています。これは、BのAに関する勾配を返す関数です。(Kerasのバックエンドについてはこちら)

    S_c = model.layers[-1].output

conv_output = model.get_layer(layer_name).output
grads = K.gradients(S_c, conv_output)[0]

今回使用するモデルはkeras.applications.resnet50にて実装されているResNet50です。ここでは、最後の層がReLU関数を使ってあるので簡単な方の高階微分が使用できます。

    #first_derivative:1階微分

first_derivative = K.exp(class_output)[0][class_idx] * grads
#second_derivative:2階微分
second_derivative = K.exp(class_output)[0][class_idx] * grads * grads
#third_derivative:3階微分
third_derivative = K.exp(class_output)[0][class_idx] * grads * grads * grads

ここで、新しいTensorFlowバックエンド関数gradient_functionを定義します(つまりはgradient_functionという関数のインスタンスを作成します)。ここでは、ResNet50に対するinputが入力されると上に出て来たconv_outputfirst_derivativesecond_derivativethird_derivativeが出力されるものとしましょう。

    gradient_function = K.function([model.input], [conv_output, first_derivative, second_derivative, third_derivative])  

conv_output, conv_first_grad, conv_second_grad, conv_third_grad = gradient_function([preprocessed_input])
conv_output, conv_first_grad, conv_second_grad, conv_third_grad = conv_output[0], conv_first_grad[0], conv_second_grad[0], conv_third_grad[0]

ここからは、$\alpha_{ij}^{kc}$を計算していきます。まずは、分母を計算しましょう。1行目で$\sum_a\sum_bA_{ab}^k$、2行目で $2\frac{\partial^2Y^c}{(\partial A_{ij}^k)^2}+\sum_a\sum_bA_{ab}^k \left( \frac{\partial^3Y^c}{(\partial A_{ij}^k)^3}\right)$ を計算しています。なお、これを分母として割り算をするため、0になっては困ります。ここでは、0になったら、全て1とする処理を噛ませておきます。

    global_sum = np.sum(conv_output.reshape((-1, conv_first_grad.shape[2])), axis=0)

alpha_denom = conv_second_grad * 2.0 + conv_third_grad * global_sum.reshape(( 1, 1, conv_first_grad.shape[2]))
alpha_denom = np.where(alpha_denom != 0.0, alpha_denom, np.ones(alpha_denom.shape))

分子については、2階微分をそのまま適用して良いですね。では、$\alpha_{ij}^k$を算出します。

    alpha_num = conv_second_grad

alphas = alpha_num / alpha_denom

以下で、$\alpha_{ij}^k$を正規化します。


alpha_normalization_constant = np.sum(np.sum(alphas_thresholding, axis = 0), axis = 0)
alpha_normalization_constant_processed = np.where(alpha_normalization_constant != 0.0, alpha_normalization_constant, np.ones(alpha_normalization_constant.shape))

alphas /= alpha_normalization_constant_processed.reshape((1,1,conv_first_grad.shape[2]))

次に、

w_k^c=\sum_i\sum_j\dfrac{\dfrac{\partial^2Y^c}{(\partial A_{ij}^k)^2}}{2\dfrac{\partial^2Y^c}{(\partial A_{ij}^k)^2}+\sum_a\sum_bA_{ab}^k\left\{\dfrac{\partial^3Y^c}{(\partial A_{ij}^k)^3}\right\}}relu \left(\dfrac{\partial Y^c}{\partial A_{ij}^k}\right)

を計算します。

    #ReLU

weights = np.maximum(conv_first_grad, 0.0)
deep_linearization_weights = np.sum((weights * alphas).reshape((-1, conv_first_grad.shape[2])),axis = 0)

最後に

L_{ij}^c = relu(\sum_kw_k^cA_{ij}^k)

を計算しましょう。正規化もしておきます。

    grad_CAM_map = np.sum(deep_linearization_weights * conv_output, axis=2)

grad_CAM_map = np.maximum(grad_CAM_map, 0)
grad_CAM_map = grad_CAM_map / np.max(grad_CAM_map)

あとは、描画の設定をしてあげましょう。

    grad_CAM_map = cv2.resize(grad_CAM_map, (target_size, target_size), cv2.INTER_LINEAR)

jetcam = cv2.applyColorMap(np.uint8(255 * grad_CAM_map), cv2.COLORMAP_JET)
jetcam = (np.float32(jetcam) + x / 2)

以上まとめて、1つの関数にしてあげましょう。


Grad_Cam_plus_plus.py

def Grad_Cam_plus_plus(input_model, layer_name, x, row, col):

model = input_model

# 前処理
X = np.expand_dims(x, axis=0)
X = X.astype('float32')
preprocessed_input = X / 255.0

# 予測クラスの算出
predictions = model.predict(preprocessed_input)
class_idx = np.argmax(predictions[0])

# 使用する重みの抽出、高階微分の計算
class_output = model.layers[-1].output
conv_output = model.get_layer(layer_name).output
grads = K.gradients(class_output, conv_output)[0]
#first_derivative:1階微分
first_derivative = K.exp(class_output)[0][class_idx] * grads
#second_derivative:2階微分
second_derivative = K.exp(class_output)[0][class_idx] * grads * grads
#third_derivative:3階微分
third_derivative = K.exp(class_output)[0][class_idx] * grads * grads * grads

#関数の定義
gradient_function = K.function([model.input], [conv_output, first_derivative, second_derivative, third_derivative]) # model.inputを入力すると、conv_outputとgradsを出力する関数

conv_output, conv_first_grad, conv_second_grad, conv_third_grad = gradient_function([preprocessed_input])
conv_output, conv_first_grad, conv_second_grad, conv_third_grad = conv_output[0], conv_first_grad[0], conv_second_grad[0], conv_third_grad[0]

#alphaを求める
global_sum = np.sum(conv_output.reshape((-1, conv_first_grad.shape[2])), axis=0)
alpha_num = conv_second_grad
alpha_denom = conv_second_grad*2.0 + conv_third_grad*global_sum.reshape((1,1,conv_first_grad.shape[2]))
alpha_denom = np.where(alpha_denom!=0.0, alpha_denom, np.ones(alpha_denom.shape))
alphas = alpha_num / alpha_denom

#alphaの正規化
alpha_normalization_constant = np.sum(np.sum(alphas, axis = 0), axis = 0)
alpha_normalization_constant_processed = np.where(alpha_normalization_constant != 0.0, alpha_normalization_constant, np.ones(alpha_normalization_constant.shape))
alphas /= alpha_normalization_constant_processed.reshape((1,1,conv_first_grad.shape[2]))

#wの計算
weights = np.maximum(conv_first_grad, 0.0)
deep_linearization_weights = np.sum((weights * alphas).reshape((-1, conv_first_grad.shape[2])))

#Lの計算
grad_CAM_map = np.sum(deep_linearization_weights * conv_output, axis=2)
grad_CAM_map = np.maximum(grad_CAM_map, 0)
grad_CAM_map = grad_CAM_map / np.max(grad_CAM_map)

#ヒートマップを描く
grad_CAM_map = cv2.resize(grad_CAM_map, (row, col), cv2.INTER_LINEAR)
jetcam = cv2.applyColorMap(np.uint8(255 * grad_CAM_map), cv2.COLORMAP_JET) # モノクロ画像に疑似的に色をつける
jetcam = (np.float32(jetcam) + x / 2) # もとの画像に合成

return jetcam



実験

前節で作成したGrad-CAM++による可視化関数を使って実際にヒートマップを描いてみましょう。

今回使う画像はこの子です。可愛い。

both.png

描くのに使用したサンプルプログラムはこちらです。

    model = ResNet50(weights = 'imagenet')

target_layer = 'activation_49'
image_path = 'cat_dog.png'
row = 224
col = 224

img = img_to_array(load_img(image_path,target_size=(row,col)))
img_GCAMplusplus = Grad_Cam_plus_plus(model, img, target_layer, args.size)
img_Gplusplusname = args.image_path+"_GCAM++_%s.jpg"%args.model
cv2.imwrite(img_Gplusplusname, img_GCAMplusplus)

その結果はこちら。

cat_dog.pngFri Dec  7 19:31:23 2018_GCAM++_%s.jpg

まず、この写真は「犬」の画像だとクラス分けされたようですね。1000種類のラベルがついたImageNetのデータベースによって学習されたモデルを使用しているので、「犬」と判断されたのは間違ってないでしょう。さて、機械学習さんはどこを見て「犬」と判断したのか。それはヒートマップの赤いところを見れば良いのです。

うん。鼻のあたりを見たんですね。

可愛い。

なぜ可愛いのかとかはまた次回のadvent calenderでお話します。


さいごに

今回はKerasを用いたGrad-CAM++の実装を行いました。

ちなみに実装してみたものはこちらのgithubにもありますので、ご覧ください。

次回、CAMやGrad-CAMと今回のGrad-CAM++の比較をしていこうと思うので、興味ある方はお楽しみに。


参考