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

はじめに

現象数理学科4年Dasonです。Qiita８記事目。
この記事は「明治大学 Advent Calendar 2018」第８日目の記事です。
本当はこの日に書く気はなかったのですが、PayPayでNintendo Switchを買ってテンションが上がった...兼ねてから実装を試みていたGrad-CAM++の実装ができたので、勢いで書いちゃうことにしました。

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

では、いきましょう。

※注意 (2019年5月追記)
この記事を最初に執筆して半年近く経過しますが、半年もすると「最新手法！」とは言えなくなってきました。ごめんなさい。釣りじゃないです。
ただ、Grad-CAM++はCAMやGrad-CAM同様一目で分かりやすい結果の描画ができたり、Grad-CAMのもついくつかの問題点を改善した可視化手法とも言えますので、十分有意義なものであると思っています。

背景

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

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


近年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}

さて、上の２つの式から$\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}

さて、ここで１つ注目してみましょう。
もし、$\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の平均について２階微分、３階微分で重み付けをした、いわばGrad-CAMの一般化版と言えるものだったんですね。

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

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

となります。

高階微分に立ち向かう

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

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

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

$S^c$を最後から２番目の層のクラス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には自動微分がライブラリに実装されてますから、すぐに計算できちゃうのです。
さて、２階微分してみましょう。

\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関数でなくて、線形な関数であったとしても２階微分は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}

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

Grad-CAM++よ！Softmax関数はいかが？

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

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

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

とします。このとき$S^k$は最後から２番目の層でクラス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で３次元から４次元に拡張した上で正規化します。

    #モデルを指定。
    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：１階微分
    first_derivative = K.exp(class_output)[0][class_idx] * grads
    #second_derivative：２階微分
    second_derivative = K.exp(class_output)[0][class_idx] * grads * grads
    #third_derivative：３階微分
    third_derivative = K.exp(class_output)[0][class_idx] * grads * grads * grads

ここで、新しいTensorFlowバックエンド関数gradient_functionを定義します(つまりはgradient_functionという関数のインスタンスを作成します)。ここでは、ResNet50に対するinputが入力されると上に出て来たconv_outputとfirst_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}$を計算していきます。まずは、分母を計算しましょう。１行目で$\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になったら、全て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))

分子については、２階微分をそのまま適用して良いですね。では、$\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)  

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

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：１階微分
    first_derivative = K.exp(class_output)[0][class_idx] * grads
    #second_derivative：２階微分
    second_derivative = K.exp(class_output)[0][class_idx] * grads * grads
    #third_derivative：３階微分
    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++による可視化関数を使って実際にヒートマップを描いてみましょう。
今回使う画像はこの子です。可愛い。

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

    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)

その結果はこちら。

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

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

さいごに

今回はKerasを用いたGrad-CAM++の実装を行いました。
ちなみに実装してみたものはこちらのgithubにもありますので、ご覧ください。

参考

• A.Chattopadhyay et al., " Grad-CAM++: Improved Visual Explanations forDeep Convolutional Networks," arXiv preprint arXiv:1710.11063v3, 2018https://arxiv.org/pdf/1710.11063.pdf
• abityac94/Grad_CAM_plus_plus https://github.com/adityac94/Grad_CAM_plus_plus
• Keras Documentation https://keras.io/ja/