スパース性と圧縮センシング その②

はじめに

前回記事では圧縮センシングの概要と二音信号を用いた復元のデモを行った。本稿では、スパース表現行列$\boldsymbol{\Psi}$に焦点を置いてさらに掘り下げた実験を行ってみる。

スパース表現の概要

自然データにみられる固有の構造は、データがある適当な座標系でスパースに表現できる。つまり、高次元の自然データが低次元の構造を持つ場合、適切な基底や辞書を使ってスパースな表現が可能になる、という前提に基づく。例えば、信号がSVDやフーリエ規程のような特定の基底に対してスパースである場合に加えて、訓練データそのものから構成されるオーバーコンプリート辞書（overcomplete dictionary）に対してもスパースで表現できることがある。

Wrightらの研究では、このスパース表現を用いた、分類のためのスパース表現（SRC: Sparse Representation for Classification）が、ノイズや遮蔽（=マスキング）がある場合でも顔認識の精度を向上させることを示した。以下から、例としてYale B databaseを用いてオーバーコンプリート辞書の構築とノイズや遮蔽があるケースでの画像認識の精度検証を行う。

Yale B databaseのデータの中身についてはこちらの記事で少し触れている。

SRCの概略

上記図のように、SRCによる分類手法はスパースベクトル$\boldsymbol{s}$が特定の人物に対応（faces行列の列添字がそれに該当）するように大きな係数をもつ持つことになる。即ち、以下の凸$\ell_1$-最小化問題を解けば良い。

\begin{align}
\begin{cases}
{\rm 目的関数:}\quad&\hat{\boldsymbol{s}}
=\mathop{\rm arg~min}\limits_{\boldsymbol{s}}\|\boldsymbol{s}\|_1 \\
{\rm 制約条件:}\quad&\boldsymbol{y}
=\boldsymbol{\Theta s}
\end{cases}
\end{align}

ここで、$\boldsymbol{\Theta}$は全ての訓練データから構成されるオーバーコンプリート辞書である。

デモコード

ライブラリのインポート

まずはライブラリのインポート

# Import libraries
import random

import numpy as np
import matplotlib.pyplot as plt
import scipy.io
from scipy.optimize import minimize
from skimage.transform import resize
from matplotlib.image import imread
plt.rcParams['figure.figsize'] = [7, 7]
plt.rcParams.update({'font.size': 18})

データの読み込み

Yale B databaseを読み込む。これはallFaces.matのfacesというキーに1つの行列として格納されている。
この行列は二次配列がflattenされたものが列ベクトルとして格納されており、そのサイズは$(192_{(高)}\times168_{(横)})\times 2,410_{(サンプル数)}$となっていることに留意

# Read data
contents = scipy.io.loadmat('../../data/allFaces.mat')
X = contents['faces']
nfaces = contents['nfaces'].reshape(-1)
m = int(contents['m'][0,0])
n = int(contents['n'][0,0])

上記コードでcfacesには38名のデータ数がリストとして格納されている。
例えば、$n$人目の最初のインデックスを取得する関数を作成して、実際に6人目（=インデックスは5）の最初のデータを出力してみると下記のようになる。

def get_offset(n, _nfaces=nfaces):
    return sum(_nfaces[:n])


plt.imshow(
    np.reshape(contents['faces'][:, get_offset(5)], (m, n)).T,
    cmap='gray'
)
plt.axis('off')
plt.show()

出力結果は以下の通り。

データセット作成

続いて、訓練データセットとテストデータセットを作成する。
バラつきはあるが、同一人物の画像データ（光の当たり具合や表情が微妙に変わっている）は60枚強あるので、ひとり当たり訓練データ数をnTrain=30、テストデータ数をnTest=20として38名中ランダムで20人（=nPeoples）で選択し抽出する。
当然、訓練データとテストデータには重複が生じ無いように配慮して、random_splitなど適当な関数を実装

# Build Training and Test sets
nTrain = 30
nTest = 20
nPeople = 20

peoples = random.sample(range(len(nfaces)), nPeople)
Train = np.zeros((m*n, nTrain*len(peoples)))
Test = np.zeros((m*n, nTest*len(peoples)))


def random_split(n, _nTrain=nTrain, _nTest=nTest):
    train_idxs = random.sample(n, _nTrain)
    remainings = [num for num in n if num not in train_idxs]
    test_idxs = random.sample(remainings, _nTest)
    return train_idxs, test_idxs


for i, p in enumerate(peoples):
    offset = get_offset(p)
    idxs = range(offset, offset+nfaces[p])
    train_idxs, test_idxs = random_split(idxs)
    Train[:, i*nTrain:(i+1)*nTrain] = X[:, train_idxs]
    Test[:, i*nTest:(i+1)*nTest] = X[:, test_idxs]

ダウンサンプリングとオーバーコンプリート辞書の作成

さて、訓練データセットとテストデータセットを作成できたところで、本節冒頭で記載したSRC概略図のダウンサンプリングを実施する。
scikit-imageのskimage.transformを用いたが、正直画像リサイズのアルゴリズムは良くわかっていない。が、どうやら局所平均を取る（厳密には何らかの方法で新しい数値を補完してるらしい）ことによるダウンサンプリングを行なっているようだ。引数anti_aliasingをTrueにすることで、前処理にガウシアンフィルタを適用し、局所平均化によるノイズを軽減している。
今回は幅を$168\rightarrow10$、高さを$192\rightarrow12$へと変更し、リサイズした訓練データを全て束ねることによって、オーバーコンプリート辞書$\boldsymbol{\Theta}$を構築する。

# Downsample Training Images (Build Theta)
M = nTrain*nPeople
sm = 10
sn = 12
Theta = np.zeros((sm*sn, M))

for k in range(M):
    temp = np.reshape(np.copy(Train[:, k]), (m, n))
    sampled = resize(temp, (sm, sn), anti_aliasing=True)
    Theta[:, k] = np.reshape(sampled, (sm*sn,))

これで作成された$\boldsymbol{\Theta}$のサイズは行$120=10_{\rm (sm)}\times12_{\rm (sn)}$、列$600=30_{\rm (nTrain)}\times20_{\rm (nPeople)}$となっている。

画像データの再構築と分類

では、6人目の最初のテストデータを復元してみよう。因みに、ランダムで選択した20名のリストの6人目は以下のような画像データである。

plt.imshow(
    np.reshape(contents['faces'][:, get_offset(peoples[5])], (m, n)).T,
    cmap='gray'
)
plt.axis('off')
plt.show()

訓練データと同様に、一度リサイズしてからflattenし、それを観測データ$\boldsymbol{y}_1$とする。

# Downsample Test Images
sample_test = nTest*5
t1 = np.copy(Test[:, sample_test])  # Clean image
temp = np.reshape(np.copy(t1), (m, n))
sampled = resize(temp, (sm, sn), anti_aliasing=True)
y1 = np.reshape(sampled, (sm*sn,))

これを用いて、凸$\ell_1$-最適化を解く。

# L1 Search, Testclean
eps = 0.01

# L1 Minimum norm solution s_L1
def L1_norm(x):
    return np.linalg.norm(x, ord=1)


constr = ({
    'type': 'ineq',
    'fun': lambda x:  eps - np.linalg.norm(Theta @ x - y1, 2)
})
x0 = np.linalg.pinv(Theta) @ y1  # initialize with L2 solution
res = minimize(
    L1_norm, x0,
    method='SLSQP',
    constraints=constr,
    options={'maxiter': 4000, 'disp': True}
)
s1 = res.x

この処理は10分程度掛かってしまうが、うまく動作すると以下のような出力を得る。

Optimization terminated successfully    (Exit mode 0)
            Current function value: 2.7742555744517174
            Iterations: 2656
            Function evaluations: 1601941
            Gradient evaluations: 2656

ここで、Current function valueが経験的に小さい値である程良い結果が得られている。これで得たスパースベクトル$\boldsymbol{s}$を描画してみる。

plt.plot(s1)
plt.show()

可視化してみると、スパースとは言えないが、インデックス=169付近でピークがあり、特定の人物を特徴付けられていることが分かる。

次に、Thetaの代わりにTrainを用いて、$\boldsymbol{s}$からオリジンサイズのデータを復元してみる。

plt.imshow(np.reshape(Train @ (s1),(m,n)).T,cmap='gray')
plt.axis('off')
plt.show()

見た目では比較できないので、元データと復元データの差分を描画

plt.imshow(np.reshape(t1 - Train @ (s1),(m,n)).T,cmap='gray')
plt.axis('off')
plt.show()

ガウシアンフィルタを前処理として加えたとは言え、やはりエッジでの差分が目立つ。とは言え、人が観察する分には十分に復元できており、特定の人物であることが判定できる精度だと感じている。

最後に誤差の定量評価を行う。選択された20名の人物に対して、復元データとの相対誤差を計算して、棒グラフで可視化してみる。

binErr = np.zeros(nPeople)
for k in range(nPeople):
    L = range(k*nTrain, (k+1)*nTrain)
    binErr[k] = np.linalg.norm(t1-Train[:, L] @ (s1[L]))/np.linalg.norm(t1)

plt.bar(range(nPeople), binErr)
plt.show()

結果は以下の通り。

インデックス=5（つまり、6人目）の箇所で誤差が明らかに小さくなっており、正解である、特定の人物である可能性が高いことが示唆できている。

マスキングデータからの再構築と分類

次に、画像に一部マスキングがある場合の元画像の復元および人物の分類を行ってみる。マスキングは口髭のデータを使って単純に直積を取ることで実現する。

# Read occulated data (fake mustache)
mustache = imread('../../data/mustache.jpg')
mustache = np.mean(mustache, -1)
mustache = (mustache/255).astype(int)
mustache = mustache.T

実際にどんな画像かを可視化してみる。

plt.imshow(mustache.T,cmap='gray')
plt.axis('off')
plt.show()

このバイナリの巻きひげデータとテストデータの直積を取る。

# Downsample Test Images with Mustache mask
t2 = np.copy(Test[:, sample_test]) * mustache.reshape(n*m)
temp = np.reshape(np.copy(t2), (m, n))
sampled = resize(temp, (sm, sn), anti_aliasing=True)
y2 = np.reshape(sampled, (sm*sn,))

plt.imshow(np.reshape(t2,(m,n)).T,cmap='gray')
plt.axis('off')
plt.show()

このテストケースを用いて、先ほど同様に$\ell_1$-最適化を行う。

## L1 Search, Mustache
eps = 100
constr = ({
    'type': 'ineq',
    'fun': lambda x:  eps - np.linalg.norm(Theta @ x - y2, 2)
})
x0 = np.linalg.pinv(Theta) @ y2  # initialize with L2 solution
res = minimize(
    L1_norm, x0,
    method='SLSQP',
    constraints=constr,
    options={'maxiter': 4000, 'disp': True}
)
s2 = res.x

この出力結果は以下の通り。

Optimization terminated successfully    (Exit mode 0)
            Current function value: 6.845300326650376
            Iterations: 3434
            Function evaluations: 2070360
            Gradient evaluations: 3434

さて、まずは得られたスパースベクトル$\boldsymbol{s}$の分布を見てみる。

plt.plot(s2)
plt.show()

この結果を見ると、顕著なピークが2つあり、その他の要素でもバラつきがあることが確認できる。
次に復元されたデータを描画する。

plt.imshow(np.reshape(Train @ (s2),(m,n)).T,cmap='gray')
plt.axis('off')
plt.show()

かなり不気味なAIチックな仕上がりになってしまった。元ニキとAIニキの差分を見てみる。

plt.imshow(np.reshape(t2 - Train @ (s2),(m,n)).T,cmap='gray')
plt.axis('off')
plt.show()

髭は勿論、ノイズなので明らかだが顔全体が比較的誤差として現れてしまっているように感じる。精度の定性評価としてはイマイチといったところだろうか。
最後に定量的な評価も観察する。

binErr = np.zeros(nPeople)
for k in range(nPeople):
    L = range(k*nTrain, (k+1)*nTrain)
    binErr[k] = np.linalg.norm(t2-Train[:, L] @ (s2[L]))/np.linalg.norm(t2)

plt.bar(range(nPeople), binErr)
plt.show()

このグラフだけ観察して、エラー最小の人物のラベルはインデックス=4と正解であるインデックス=5と異なっているが、5の方も次いで小さい値を取っている。
その他では、インデックス=14、16あたりが特徴として目立っている。この巻きひげのようなマスキングのケースは既存の辞書だけで表現するのは少し分が悪いといった所感である。

オクルージョンデータからの再構築と分類

今度は、画像データの一部のピクセルを選択してランダムな値に置き換える、オクルージョンと呼ばれる欠損処理を施したテストデータを用いる。

# Occlude Test Image
def occlude(x, ratio=0.3):
    randvec = np.random.permutation(n*m)
    first_ratio = randvec[:int(np.floor(ratio*len(randvec)))]
    vals_ratio = (255*np.random.rand(*first_ratio.shape)).astype(int)
    _x = np.copy(x)
    _x[first_ratio] = vals_ratio
    return _x


t3 = occlude(t1)
temp = np.reshape(np.copy(t3), (m, n))
sampled = resize(temp, (sm, sn), anti_aliasing=True)
y3 = np.reshape(sampled, (sm*sn,))

このように画像の30%のピクセルがランダムに選択され、0~255の適当の値に置換される。

## L1 Search, Occlusion
eps = 100
constr = ({
    'type': 'ineq',
    'fun': lambda x:  eps - np.linalg.norm(Theta @ x - y3, 2)
})
x0 = np.linalg.pinv(Theta) @ y3  # initialize with L2 solution
res = minimize(
    L1_norm, x0,
    method='SLSQP',
    constraints=constr,
    options={'maxiter': 4000, 'disp': True}
)
s3 = res.x

Optimization terminated successfully    (Exit mode 0)
            Current function value: 0.8741597687516052
            Iterations: 2739
            Function evaluations: 1652772
            Gradient evaluations: 2739

得られたスパースベクトル$\boldsymbol{s}$の分布を見てみる。

plt.plot(s3)
plt.show()

スパース性は良好であるが、複数の要素でピークが立っており、特定の個人だけで特徴付けられているとは言い難い。
次に、復元データおよび元データとの差分データを可視化してみる。

plt.imshow(np.reshape(Train @ (s3),(m,n)).T,cmap='gray')
plt.axis('off')
plt.show()

plt.imshow(np.reshape(t3 - Train @ (s3),(m,n)).T,cmap='gray')
plt.axis('off')
plt.show()

出力結果はそれぞれ以下の通り。

差分の可視化結果を見ると、加えたノイズとエッジは目立つが総じてエラーは小さそうである。
最後に定量評価を観察

binErr = np.zeros(nPeople)
for k in range(nPeople):
    L = range(k*nTrain, (k+1)*nTrain)
    binErr[k] = np.linalg.norm(t1-Train[:, L] @ (s3[L]))/np.linalg.norm(t3)

plt.bar(range(nPeople), binErr)
plt.show()

スパースベクトル$\boldsymbol{s}$が複数のピークを持っていたので、エラーもまちまちだが、インデックス=5の人物が最小となっているので、結果としては悪くは無いようである。

おわりに

今回は圧縮センシングとオーバーコンプリート辞書を用いた画像の分類デモンストレーションを行った。ダウンサンプリングと、取得した情報はかなり少ないものの、人間が観察する分には復元精度は割と良さそうである。
また、ノイズやマスキングを加えたテストデータによる検証も行ったが、ノイズにはある程度頑健だが、マスキングには向かないことが分かった。