전산물리학:압축_센싱 [statphys]

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
======개요======
섀넌-나이키스트 표본화 정리(Shannon-Nyquist sampling theorem)에 의하면, 신호의 최대 주파수가 B일 때 이를 온전히 재구성하려면 2B 이상의 주파수로 표본을 추출해야 한다.그런데 우리가 보는 대부분의 신호는 큰 패턴이 중요한 역할을 하므로 주파수 밀도 측면에서 희박(sparse)하고, 이를 이용하면 섀넌-나이키스트 표본화 정리가 요구하는 것보다 적은 수의 표본으로도 신호를 재구성할 수 있다. 이러한 아이디어는 Candès, Romberg, and Tao (2006)에 의해 정식화되었고 압축 센싱(compressed sensing)이라 불리고 있다.

======기본적인 설명======
길이 $N$인 '희박한' 신호 $\mathbf{x}$가 미지수로서 있다. $k\times N$의 '관찰' 행렬 $A$를 곱하자. 이 행렬의 원소는 난수로 주어져 있고 $k \ll N$이라고 가정한다. 그 곱의 결과로 $k$ 차원 벡터 $\mathbf{y} = A \mathbf{x}$를 관찰 결과로 가지고 있다.

$A$와 $\mathbf{y}$만을 가지고 미지수 $\mathbf{x}$를 복원하는 것은 불가능해 보인다. $\mathbf{y}=A\mathbf{x}$를 만족하는 $\mathbf{x}$가 무수히 많기 때문이다. 하지만 그러한 $\mathbf{x}$ 중에서 가장 '희박한' 것, 구체적으로는 $\sum_i |x_i|$를 최소화하는 $\mathbf{x}$를 택하면 원래의 $\mathbf{x}$를 상당히 잘 복원해낼 수 있다.

=====코드 구현=====
아래는 참고문헌 항목 중 codeproject.com에 있는 Matlab 코드를 python으로 옮긴 것이다.
<code:python>
from __future__ import print_function,division
from numpy import zeros,dot,transpose
from numpy.random import normal,choice
from pylab import plot,show
from l1eq_pd import l1eq_pd

# Initialize constants and variables
N = 256             # length of signal
P = 5               # number of non-zero peaks
K = 64              # number of measurements to take
x = zeros(N,float)  # original signal (P-sparse)

# Generate signal with P randomly spread values
peaks = choice(N, P, replace=False)
x[peaks] = normal(0.0, 1.0, P)

# Obtain K measurements
A = normal(0.0, 1.0, (K,N))
y = dot(A,x)

# Perform Compressed Sensing recovery
x0 = dot(transpose(A), y)
xp = l1eq_pd(x0, A, [], y)
plot(x)
show()
plot(xp-x)
show()
</code>
신호의 길이는 $N=256$인데 이 중 $P=5$개 지점에서만 0이 아닌 값을 가진다. 이 신호는 그런 의미에서 '희박'하다. $K=64$는 관찰의 횟수로서 $N$의 1/4에 불과하다. choice 명령을 써서 5개의 점을 임의로 골라내고 그 곳들의 값을 정규(normal) 분포로부터 결정한다.

관찰 행렬 $A$도 정규 분포로부터 만든 다음, 관찰 결과인 $\mathbf{y}=A\mathbf{x}$를 만들어낸다. 이제 $\mathbf{x}$를 모른다고 가정한 채, $A$와 $\mathbf{y}$로부터 신호를 복원해낸다. 여기에서 절대값 합의 최소화를 위해 쓰인 l1eq_pd라는 함수는 아래 페이지에서 내려받을 수 있다:
https://github.com/nikcleju/pyCSalgos

실행 결과는 다음과 같다. 첫 번째 그림은 원래 신호를, 두 번째 그림은 복원된 신호와 원래 신호의 차이를 보여준다. 그 차이는 매우 작다.
{{::전산물리학::sensing_x.png?400|}}{{::전산물리학::sensing_d.png?400|}}

=====푸리에 공간에서 희박한 신호=====
이 예는 참고문헌 항목 중 pyrunner.com의 내용으로부터 일부를 수정하여 가져온 것이다.
실공간 상에서는 조밀하지만 실제로는 몇 개의 주파수만이 중요한 신호를 생각해보자. 예컨대 다음과 같은 신호는 단 2개의 주파수 성분만을 지닌다:
$$x = \sin(10 \pi t) + \sin (54 \pi t).$$
이산 사인 변환(discrete sine transform, DCT)를 통해 이 신호를 푸리에 공간으로 옮기고 압축 센싱을 실행한 다음 역변환을 할 것이다. 그런데 앞의 코드를 재활용하기 위해 푸리에 변환을 행렬 곱 형태로 표현하면 편리하다 (계산 시간을 생각하면 비효율적인 선택일 것이다):
$$\tilde{x} \equiv \mathcal{F}[x] = \Psi \mathbf{x},$$
$$x \equiv \mathcal{F}^{-1}[\tilde{x}] = \Psi^{-1} \tilde{\mathbf{x}}.$$
이 떄 $x$는 원래 신호, $\tilde{x}$는 그것의 푸리에 변환($\mathcal{F}$) 결과이고 굵은 글씨로 적은 것은 신호를 벡터로서 적은 것을 의미한다. $\Psi$가 푸리에 변환의 역할을 해주는 행렬이다.

코드는 다음처럼 작성할 수 있다.
<code:python>
from __future__ import print_function,division
from math import pi
from numpy import sin,dot,transpose,identity,linspace
from numpy.random import choice
from pylab import plot,show
from l1eq_pd import l1eq_pd
from scipy.fftpack import dst,idst

# Initialize constants and variables
n = 600            # length of signal
m = 50             # number of measurements to take

# Generate signal with P randomly spread sinusoids
t = linspace(0,1,n)
x = sin(10*pi*t) + sin(54*pi*t)

# Orthonormal basis matrix
Psi = dst(identity(n), norm='ortho', axis=0)
Psi_inv = idst(identity(n), norm='ortho', axis=0)
xt = dot(Psi,x)

# Extract small sample of signal
ri = choice(n, m, replace=False)
ri.sort()
t2, x2 = t[ri], x[ri]
plot(t,x,'b')
plot(t2,x2, 'r.')
show()

A = Psi_inv[ri]      # sensing matrix
y = dot(A,xt)        # measured values

# Perform Compressed Sensing recovery
x0 = dot(transpose(A), y)
X_hat = l1eq_pd(x0, A, [], y)

x_rec = dot(Psi_inv,X_hat)     # IFFT of X_hat(f)
plot(x)
plot(x_rec, ".")
show()
</code>
코드의 실행 결과는 아래 그림들처럼 나온다.

{{::전산물리학::sensing_fourier1.png?400|}}
{{::전산물리학::sensing_fourier2.png?400|}}

첫 번째 그림은 원래 신호와 그 신호를 추출하기 위해서 무작위로 고른 위치들을 보여준다. 두 번째 그림은 원래 신호(실선)와 복원된 신호(점)를 보여주고 있다. 새넌-나이키스트 정리에서 요구하는 것보다 적은 수의 점을 추출했음에도 불구하고 복원된 신호는 원래 신호와 상당히 가깝다.

=====통계물리적 접근=====

2012년, L1 norm에 기반한 방법보다 더 적은 데이터를 가지고도 신호를 복원하는 방법이 개발되었다 [Krzakala et al., PRX 2, 021005 (2012)].


======참고문헌======
  * https://www.codeproject.com/Articles/852910/Compressed-Sensing-Intro-Tutorial-w-Matlab
  * http://www.pyrunner.com/weblog/2016/05/26/compressed-sensing-python/
  * https://dilawarnotes.wordpress.com/2017/09/06/compressed-sensing-a-python-demo/
  * https://doi.org/10.1103/PhysRevX.2.021005