0. はじめに

圧縮センシングは、観測ベクトルをy、あるスパースなベクトルをxとし、y = Axとモデリングできる場合に、xがスパースであることを条件にyからxを復元する技術です。
圧縮センシングを利用するときの注意点は、観測ベクトルyの長さが十分でないと、xの復元に必要な情報を得ることができずに復元に失敗することです。
観測ベクトルyの長さ(M)に関しては、以下のように理論的な目安（復元成功の理論閾値）が示されています[1]。

M ≧ C*S*log(S/N) \:\:\:\:\:\:\:\:・・・ (1)

ここで、Nはベクトルxの長さ、Sはxの非ゼロ成分の数、Cは固定値です。
復元成功の理論閾値は示されていますが、行列Aの構成と圧縮センシングのアルゴリズムによって、復元に失敗する場合があります。
本稿では、圧縮センシングのアルゴリズムの一つである、Orthogonal Matching Pursuit (OMP) を対象に、復元に失敗する確率をシミュレーションで確認します。OMPによるスパースなベクトルxの復元はxとAの構成に依存するため、シミュレータでは、双方をランダムな複素数の値で構成することを繰り返し、確率的に復元の成否を求めます。

【コラム】OMPにおいてxの復元に失敗することの感覚的な説明

OMPは、yとAの相関が強いAの列を順番に抽出していくことで、元のスパースなベクトルxを復元するアルゴリズムです。
OMPで復元に失敗する場合を説明します。例えば、yがAの第一列と第二列で構成されることが正解であった場合に、yとAの第三列の相関が運悪く高いと、OMPの初期段階でAの第三列が抽出されてしまい、正解の列の組み合わせを抽出できず、復元に失敗してしまいます。
一方で、yの長さが長くなると、上記のような誤った相関が発生しにくいため、復元の失敗を避けられます。

1. シミュレータの仕様

シミュレータのプログラムを示す前に、シミュレータの仕様を以下に示します。

y, x, Aは複素数のベクトルまたは行列とする
圧縮センシングのアルゴリズムとしてOMPを利用する
xやAの要素のRealとImaginary成分は、(-1,1]の範囲でランダムな値とする
観測ベクトルyの長さ(M)をSから100まで変えながら、OMPによるxの復元成功確率を計算する
あるMの値におけるxの復元成功確率を求める際には、Aをランダムな値で置き換えることを100回行う
引数として、ベクトルxの長さ、ベクトルxの非ゼロ成分の数を設定可能にする
シミュレーションの結果と共に、復元成功の理論閾値もグラフとしてプロットする

2. シミュレータの実装

上記に示した仕様に基づき、Pythonで実装したプログラムを示します。

main.py

import sys
import numpy as np
import matplotlib.pyplot as plt
from omp import omp

np.random.seed(1159)

def main(N,S):
	mse_threshold = 1e-10 # Threshold for MSE
	Nitr = 100 # Iteration 
	MSEL = []
	ML = []

	# M : Length of output vector
	for M in range(S,100):
		vote = 0
		for itr in range(Nitr):
			Weight = 2.* (np.random.random([M,N]) + 1j *np.random.random([M,N]))
			Weight -= np.ones([M,N]) + 1j* np.ones([M,N])
			Xs = np.zeros(N,dtype=complex)
			for i in range(S):
				Xs[i*3] = 2.*(np.random.random(1) + 1j * np.random.random(1)) - 1 - 1j
			Yr = np.dot(Weight,Xs)
			MSE,Xe,idx = omp(Weight,Yr,S)
			if MSE > mse_threshold:
				vote = vote + 1
		ML.append(M)
		MSEL.append(vote/Nitr*100)
	MSEth = (S * np.log(N/S))
	print(MSEth)
	print(MSEL)

	# Theoretical 
	TH = np.int32(np.round(MSEth,0))
	MSETHL = []
	for k in range(S,TH):
		MSETHL.append(100.)

	for k in range(TH,100):
		MSETHL.append(0.)
	
	plt.plot(ML,MSEL,color="b",label="simulation")
	plt.plot(ML,MSETHL,color="r",label="theoretical bound")
	plt.xlim(S,100)
	plt.ylim(0.,100.)
	plt.xlabel("Length of Received Signals")
	plt.ylabel("Prob. in Failure  [%]")
	plt.legend()
	plt.show()


if __name__ == '__main__':
	np.set_printoptions(suppress=True)
	args = sys.argv
	if len(args) == 3:
		N  = int(args[1])
		S  = int(args[2])
	else:
		print("Set  Arguments as int")
		print("#1 : the number of input vector")
		print("#2 : the number of sparsity in input vector")
		print("example :  python3 main.py 256 8")
		sys.exit(0)
	main(N,S)

omp.py

import numpy as np

# y = A * x
def omp(A,y,n_sparse):
	r = y.copy()
	L = []
	L0 = A.shape[0]
	L1 = A.shape[1]
	X = np.zeros(L1,dtype=complex)
	Vacancy = np.zeros(L0,dtype=complex) 
	A1 = (A.T).copy()
	
	for k in range(n_sparse):
		p = np.argmax(np.abs(np.dot(np.conjugate(A1),r)))
		L.append(p)
		if k == 0:
			A0 = A1[p].copy()
			IAA = 1./(np.dot(np.conjugate(A0),A0.T))
		else:
			A0 = np.vstack((A0,A1[p]))
			IAA = np.linalg.inv(np.dot(np.conjugate(A0),A0.T))

		IA = np.dot(IAA,np.conjugate(A0))
		z = np.dot(IA,y)
		r = y - np.dot(A0.T,z)
		A1[p] = Vacancy 
	i = 0
	for idx in L:
		X[idx] = z[i]
		i = i + 1
	MSE = np.sum(np.abs(y - np.dot(A,X)))
	return MSE, X, L

3. シミュレータの実行結果

シミュレータの第一引数はベクトルxの長さ、第二引数はベクトルxの非ゼロ成分の数です。
以下では、2パターンの引数でシミュレーションした結果をグラフ表示しています。

python main.py 100 3

python main.py 256 8

グラフの横軸は観測ベクトルyの長さ、グラフの縦軸はxの復元に失敗する確率です。
グラフの青線はシミュレーションの結果です。グラフの赤線は、復元成功の理論閾値をC=1の条件で計算したものです。
グラフの横軸の左端で復元に必ず成功する理由は、行列Aのランクとxの非ゼロ成分の数が一致しているために、一意に復元できるからです。
観測ベクトルyの長さが短いと、OMPの性質により復元に失敗しますが、ある程度の長さになると、復元に成功し始める様子が確認できます。この復元に成功し始める境として、赤線で示した復元成功の理論閾値がよく対応していることが確認できます。従って、観測ベクトルの長さを設定する際には、期待する復元失敗率よりも低くなるように、復元成功の理論閾値やシミュレーション結果を参考にしながら決めると良いと言えます。

4. まとめ

本稿では、圧縮センシングのアルゴリズムの一つであるOMPを対象に、観測ベクトルの長さと復元に失敗する確率をシミュレーションや理論値で評価しました。

参考文献

[1] https://authors.library.caltech.edu/10092/1/CANieeespm08.pdf

前回の記事