【圧縮センシング】複素数に対応したOrthogonal Matching Pursuitの紹介

0. はじめに

圧縮センシングとは、観測したデータがスパースなベクトルの線形写像としてモデリングできる場合に、そのスパースなベクトルを観測データから探索する技術です。
Orthogonal Matching Pursuit(OMP)は、圧縮センシングの代表的なアルゴリズムです。
Pythonでは、sklearnにてOMPが提供されています[1]。しかし、データ型がRealである必要があります。
本稿では、データ型が複素数であっても動作するOMPのPythonプログラムを紹介します。
尚、今回実装したOMPのアルゴリズム詳細を知りたい方は、文献[2]の"The Naive Implementation of OMP"を参照ください。

1. OMPプログラム

複素数に対応したOMPのプログラムを示します。

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)))
	if DEBUG == 1:
		print(A1[p]) 
		print(A[:,p]) 
		print(np.sum(np.abs(A.T-A1)))
		print(np.dot(A,X))
	return MSE, X, L

2. テスト

2.1 テスト用のプログラム

テスト用のプログラムを示します。
プログラムでは、観測データの生成やOMPによる正解データの推定、元の正解データと推定したデータのグラフ表示を行います。尚、今回の正解データは、512要素の内、3要素のみが非ゼロの複素数であるスパースなデータとしました。

test.py

import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import make_sparse_coded_signal

n_components, n_features = 512, 100
n_nonzero_coefs = 3

# generate the data
# y = Xw
# |x|_0 = n_nonzero_coefs
y, X, w = make_sparse_coded_signal(
    n_samples=1,
    n_components=n_components,
    n_features=n_features,
    n_nonzero_coefs=n_nonzero_coefs,
    random_state=0,
    data_transposed=True,
)

w = np.zeros(w.shape[0],dtype=complex)
w[10] = 1.0+1j*0.5
w[100] = -0.5+1j*1.0
w[200] = 0.5-1j*0.5
(idx,) = w.nonzero()

y = np.dot(X,w)

# plot the sparse signal
plt.figure(figsize=(8, 12))
plt.subplot(2, 2, 1)
plt.xlim(0, 512)
plt.title("Real Part of Sparse signal")
plt.stem(idx, w[idx].real)
plt.subplot(2, 2, 2)
plt.xlim(0, 512)
plt.title("Imaginary part of Sparse signal")
plt.stem(idx, w[idx].imag)

from omp import omp as myomp
MSE, coef, idx_r = myomp(X,y,n_nonzero_coefs)

plt.subplot(2, 2, 3)
plt.xlim(0, 512)
plt.title("Real Part of Recovered signal")
plt.stem(idx_r, coef[idx_r].real)
plt.subplot(2, 2, 4)
plt.xlim(0, 512)
plt.title("Imaginary Part of Recovered signal")
plt.stem(idx_r, coef[idx_r].imag)
plt.show()

2.2 実行結果

プログラムを実行すると(python test.py)、以下のグラフが表示されます。
グラフの上段は、元の正解データのReal（左）とImaginary（右）をプロットしたものです。
グラフの下段は、ompで推定したデータのReal（左）とImaginary（右）をプロットしたものです。
ompで推定したデータが正しく復元できていることが確認できます。

3. まとめ

複素数のデータ型に対応したOMPのプログラムを紹介し、正しく動作することを示しました。

参考文献

[1] https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.OrthogonalMatchingPursuit.html
[2] https://wnt.sjtu.edu.cn/papers/conf/vtc13-hufeizhu-efficient.pdf