センサのスパース配置法: Sparse Sensor Placement

はじめに

本稿はこちらの参考書の3.8節、Sparse Sensor Placementの備忘。
後続のROMと関連が深い節なので、これも１投稿としたい。

概要

前回記事までは、ランダム測定値を用いてフーリエ変換/ウェーブレット変換などの一般基底から信号を復元することについて考えてきた。ところが、復元対象が例えば、「人の顔」であるという事前知識がある場合、SVDによる構築される特定の特徴量ライブラリ、あるいは縮約スパース表現行列$\boldsymbol{\Psi}_r=\tilde{\boldsymbol{U}}$に基づいてセンサを最適化することで、必要なセンサ数を大幅に減らすことができる（ここでいうセンサ配置は元信号からどうデータをサンプリングするか、と同義と捉えて良い）。これにより、汎用的な特徴量ライブラリ（同じ基底を使って、人の顔やコーヒーの画像を表現できたような）にランダムにセンサ配置する手法とは対照的に、特定のライブラリに合わせたセンサ配置を設計することが可能となる。これに類似した手法がROM（Reduced-Order Models）のサンプリング効率化としても用いられており、ハイパーリダクション（hyper reduction）と称されている。

本稿では、Pythonを用いてSparse Sensor Placementによる信号復元をデモしてみる。

Sparse Sensor Placement for Reconstruction

スパース表現行列$\boldsymbol{\Psi}_r\in\mathbb{R}^{n\times r}$によるセンサ配置最適化は、次の連立一次方程式の行列表現

\begin{align}
\boldsymbol{y}=
\boldsymbol{C\Psi}_r\boldsymbol{a}=
\boldsymbol{\Theta a}
\end{align}

の逆行列の条件数（condition number）を可能な限り減らすようにセンシング行列$\boldsymbol{C}\in\mathbb{R}^{p\times n}$を設計することである。そうすることで、ノイズを含んだ測定データ$\boldsymbol{y}$を元にして、元々の低ランク構造を持つ系の係数$\boldsymbol{a}$を特定するために逆行列を使うことができる。

条件数について

ここで、一旦条件数について少し掘り下げてみる。

行列$\boldsymbol{\Theta}$の条件数はその最大特異値と最小特異値の比率を指し、行列演算（乗算や逆行列）が入力誤差に対してどれだけ敏感かを示す。具体的には、行列$\boldsymbol{\Theta}$を特異値分解して得られた特異値$\lbrace\sigma_i\rbrace_i^n$で最大のものを$\sigma_{\rm max}$、最小のものを$\sigma_{\rm min}$とすれば、条件数$\kappa\left(\boldsymbol{\Theta}\right)$は以下のように表せられる。

\begin{align}
\kappa
\left(
    \boldsymbol{\Theta}
\right)
=
\cfrac{
    \sigma_{\rm max}
}{
    \sigma_{\rm min}
}
\end{align}

例えば、式$(1)$で入力信号$\boldsymbol{a}$にノイズ$\epsilon_{\boldsymbol{a}}$が加わったケースを考えてみると、

\begin{align}
\Delta \boldsymbol{y}
=
\boldsymbol{\Theta}\left(
\boldsymbol{a}+\epsilon_{\boldsymbol{a}}
\right)
-
\boldsymbol{\Theta a}
=\boldsymbol{\Theta}\epsilon_{\boldsymbol{a}}
\end{align}

行列$\boldsymbol{\Theta}$の特異値分解が$\boldsymbol{\Theta}=\boldsymbol{U\Sigma V}^*$という形で表されるとして、ノイズの出力変化$|\Delta\boldsymbol{y}|$は以下のように表現できる。

\begin{align}
\|\Delta\boldsymbol{y}\|
=\|\boldsymbol{U\Sigma V}^*\epsilon_{\boldsymbol{a}}\|
\end{align}

行列$\boldsymbol{U}$と$\boldsymbol{V}$はユニタリであることから、

\begin{align}
\|\Delta\boldsymbol{y}\|
&=\|\boldsymbol{\Sigma}\epsilon_{\boldsymbol{a}}\| \\
\Rightarrow
\sigma_{\rm min}\|\epsilon_{\boldsymbol{a}}\|
\leq
&\|\Delta \boldsymbol{y}\|
\leq
\sigma_{\rm max}\|\epsilon_{\boldsymbol{a}}\| \\
\Rightarrow
\sigma_{\rm min}\|\epsilon_{\boldsymbol{a}}\|
\leq
&\|\Delta \boldsymbol{y}\|
\leq
\kappa
\left(
    \boldsymbol{\Theta}
\right)
\cdot
\sigma_{\rm min}\|\epsilon_{\boldsymbol{a}}\| 
\end{align}

この式$(7)$の右側不等式に着目すると、入力ノイズに対する出力変化$|\Delta \boldsymbol{y}|$は条件数$\kappa\left(\boldsymbol{\Theta}\right)$と最小特異値$\sigma_{\rm min}$に依存することが分かる。つまり、条件数が小さい（=$1$に近づく）場合、出力変化は$\sigma_{\rm min}$のみに依存し、ノイズ影響を抑制することが出来る。

センサ配置の最適化と近似手法

センサ数$p$とスパース表現行列$\boldsymbol{\Psi}_r$の階数$r$が等しい時、行列$\boldsymbol{\Theta}=\boldsymbol{C\Psi}_r$は正方行列となり、この条件数$\kappa\left(\boldsymbol{\Theta}\right)$がなるべく小さくなるようにセンシング行列$\boldsymbol{C}$を設計する。
一方、$p>r$のとき、擬似逆行列で使われる$\boldsymbol{M}=\boldsymbol{\Theta}^T\boldsymbol{\Theta}$の条件数最適化を図るが、これはNP困難であり、可能なセンサ配置の全ての組み合わせに対して探索が必要となる。

この問題を解くヒューリスティックな手法として、凸最適化や半正定値計画法といった反復手法が存在するが、これらは毎イテレーション$n\times n$行列の因子分解をするため計算コストが高く付く。そのため、一般的には貪欲アルゴリズムが用いられる。

ランダムなスパースセンサ配置

一般的に、センサのランダム配置はモード係数$\boldsymbol{a}$を推定するために使用される。しかし、$p=r$のとき条件数$\kappa\left(\boldsymbol{\Theta}\right)$は非常に大きくなる傾向にある。こちらの記事で紹介したように、オーバーサンプリングを採用することで条件数を大幅に改善することが見込める。

QRピボット法によるスパースセンサ配置

行列$\boldsymbol{\Psi}^T_r$のQR分解を列ピボット選択を導入して行うことで、簡単かつ効果的なセンサ配置最適化を実現できる。具体的には次のような形で分解する。

\begin{align}
\boldsymbol{\Psi}^T_r\boldsymbol{C}^T=\boldsymbol{QR}
\end{align}

ここで、（同じ文字を使ってややこしいのだが）$\boldsymbol{C}$は列交換行列、$\boldsymbol{Q}$は直交行列、$\boldsymbol{R}$は上三角行列である。列交換行列$\boldsymbol{C}$による列順序の変更により、QR分解の数値的安定性向上を図る。

デモコード

Extended Yale Face Database Bを用いて、$p=r$時のQRピボット法およびランダム法によるセンサ配置を行い、それぞれで画像復元して結果を比較する。

大まかな流れとしては以下の通り。

• データ行列$\boldsymbol{X}$を訓練データとテストデータに分割（今回は38人目だけをテストデータに。残りは訓練用とする。）
• 事前知識（=人の顔の特徴量が欲しい）として、訓練データを特異値分解してスパース表現行列$\boldsymbol{\Psi}_r$を作成
• 作成したスパース表現行列$\boldsymbol{\Psi}_r$を用いて、テストデータを圧縮センシングで復元

先ずはライブラリのインポートから。

# Import libraries
import numpy as np
import matplotlib.pyplot as plt
import scipy.io
from tqdm import tqdm
plt.rcParams['figure.figsize'] = [7, 7]
plt.rcParams.update({'font.size': 18})

続いてデータの読み込み

# 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])

QRピボット法によるセンサ配置のクラスを作成する。
今回、スパース表現行列$\boldsymbol{\Psi}_r$は訓練データ$\boldsymbol{X}_{\rm train}$のSVDを実施して、階数$r=100$とした時の左特異ベクトル縮約行列$\tilde{\boldsymbol{U}}$から得る。

\begin{align}
\boldsymbol{X}_{\rm train}
&\simeq
\tilde{\boldsymbol{U}}\tilde{\boldsymbol{\Sigma}}\tilde{\boldsymbol{V}}^* \\
\Rightarrow
\boldsymbol{\Psi}_r&=\tilde{\boldsymbol{U}}
\end{align}

そして、この$\boldsymbol{\Psi}_r$に対して、QRピボット法を適用する。用いる関数はscipy.linalg.qrで、オプションpivoting=trueとすることで第3返り値に並び替えられたインデックスのリストが来る。それを元に列交換行列$\boldsymbol{C}$を生成

そして、テストデータの再構築$\tilde{\boldsymbol{X}}_{\rm test}$はスパースベクトル$\boldsymbol{a}$を用いて、

\begin{align}
\tilde{\boldsymbol{X}}_{\rm test}
&\simeq\boldsymbol{\Psi}_r\boldsymbol{a} \\
&=\boldsymbol{\Psi}_r
\left(
    \boldsymbol{\Theta}^\dagger\boldsymbol{y}
\right)
\end{align}

と表せられる。
測定値$\boldsymbol{y}$は選択したセンサとランクからX[pivot[:rank]]などと抽出できる。

# QR pivoting for Sparse Sensors
class QRplacement:

    def __init__(self, X, rank):
        self.X = X
        self.rank = rank
        self._SVD()
        self._QR()

    def _SVD(self):
        self.U, self.S, self.Vh = np.linalg.svd(X, full_matrices=False)
        self.Psi_r = self.U[:, :self.rank]

    def _QR(self):
        self.Q, self.R, self.piv = qr(self.Psi_r.T, pivoting=True)
        self.C = np.zeros_like(self.Psi_r.T)
        
        for col in range(self.rank):
            self.C[col, self.piv[col]] = 1
            
        self.Theta = np.dot(self.C, self.Psi_r)
        self.invTheta = np.linalg.pinv(self.Theta)
    
    def reconstruct(self, X):
        y = X[self.piv[:self.rank]]
        return self.Psi_r @ self.invTheta @ y

同様にして、ランダムセンサ配置を模倣したクラスを作成。QRピボット法で列交換行列$\boldsymbol{C}$を用いて上位rankを選択しているところを乱数にするだけで良い。

# Rondom Sparse Sensor placement
class Randplacement:

    def __init__(self, X, rank):
        self.X = X
        self.rank = rank
        self._SVD()
        self._random_selection()

    def _SVD(self):
        self.U, self.S, self.Vh = np.linalg.svd(X, full_matrices=False)
        self.Psi_r = self.U[:, :self.rank]

    def _random_selection(self):
        np.random.seed(0)
        self.selected_sensors = np.random.choice(self.X.shape[0], self.rank, replace=False)
        self.C = np.zeros_like(self.Psi_r.T)
        
        for col, sensor in enumerate(self.selected_sensors):
            self.C[col, sensor] = 1
        
        self.Theta = np.dot(self.C, self.Psi_r)
        self.invTheta = np.linalg.pinv(self.Theta)
    
    def reconstruct(self, X):
        y = X[self.selected_sensors]
        return self.Psi_r @ self.invTheta @ y

上記2手法を階数$r=100$を指定して、訓練データを食わせる。
テストデータの最初のインデックスをoffsetとして、

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

offset = get_offset(37) # ゼロオリジンなので、これが38人目の最初の画像

qr_placement = QRplacement(X[:,:offset], rank=100)
res_qr = qr_placement.reconstruct(X[:,offset:])

rand_placement = Randplacement(X[:,:offset],rank=100)
res_rand = rand_placement.reconstruct(X[:,offset:])

結果を出力

# Create a figure with 3 subplots
fig, axes = plt.subplots(1, 3, figsize=(13, 5))

# Plot the original image
axes[0].imshow(np.reshape(X[:, offset], (m, n)).T, cmap="gray")
axes[0].axis("off")
axes[0].set_title("Original")

# Plot the QRplacement result
axes[1].imshow(np.reshape(res_qr[:, 0], (m, n)).T, cmap="gray")
axes[1].axis("off")
axes[1].set_title("QR Pivoting")

# Plot the Randplacement result
axes[2].imshow(np.reshape(res_rand[:, 0], (m, n)).T, cmap="gray")
axes[2].axis("off")
axes[2].set_title("Random")

# Show the figure
plt.tight_layout()
plt.show()

乱数法よりもQRピボット法による再構築画像の方が綺麗になっている。

おわりに

参考にした書籍ではコードもなく、テキストだけだったので、理解やコーディングが誤っている部分があるかも知れない。
もし見つけた場合は指摘お願いします。