匹配追蹤

**匹配追蹤**(matching pursuit, MP)最早是時頻分析的分析工具，目的是要將一已知訊號拆解成由許多被稱作為原子訊號的加權總和，而且企圖找到與原來訊號最接近的解。其中原子訊號為一極大的原子庫中的元素。以數學式子表示可以得到：

:  f(t) = \sum_{n=0}^{+\infty} a_n g_{\gamma_n}(t)

其中，a_n是權重，g_{\gamma_n}是由字典D中獲得的原子訊號。

如同傅立葉級數將一訊號拆解成一系列的正弦波的相加，其中每個成分擁有不同的係數作為權重，其數學式子如下：

:f(x) = \sum_{n=-\infty}^ \infty c_n e^{inx}

而匹配追蹤也具有將訊號拆解成一系列原子相加的意涵，甚至可以使用匹配追蹤去描述傅立葉級數，也就是原子庫對應到的所有正弦函數的集合。

## 演算法
為了找到最符合原訊號的一組原子加權總合，如果對原子庫進行所有組合的嘗試過於耗費時間。在1993年由Mallat S和Zhang Z的論文[1]中，提出了一個貪婪演算法（Greedy Algorithm），並大幅降低找出近似解的時間。其作法首先在原子庫中尋找與原訊號內積結果最大的原子g_{\gamma_n}，找到此訊號以及其內積結果a_n之後再將原訊號減掉a_n g_{\gamma_n}作為下一次重複運算的原始訊號，如此反覆做下去即可得到一系列的a_n以及原子g_{\gamma_n}，直到達到停止條件為止，其詳細的演算法如下：

:  輸入：Signal: f(t), dictionary D.
:  輸出：List of coefficients:   \left( a_n, g_{\gamma_n}\right) .
:  初始化：
:   R_1\,\leftarrow\,f(t);
:   n\,\leftarrow\,1;
:  重複：
:    寻找 g_{\gamma_n} \in D 具有最大内积  | \langle R_n, g_{\gamma_n} \rangle | ;
:     a_n\,\leftarrow\,\langle R_n, g_{\gamma_n}\rangle ;
:     R_{n+1}\,\leftarrow\,R_n - a_n g_{\gamma_n};
:     n\,\leftarrow\,n + 1;
:直到達到停止條件（例如：\|R_n\| ）

## 時頻原子分解(time-frequency atom decomposition)
在信號處理的許多應用中，需要將信號分解為一群在時域和頻域都具有良好局部性(集中在某一範圍)的函數，這些函數稱為時頻原子(time-frequency atom)。

選擇不同的時頻原子時，分解方式的特性會有很大的差異。窗函數傅立葉轉換(window Fourier transform)和小波轉換(wavelet transform)都是時頻信號分解的方法。

通常一個時頻原子群可以由單一的窗函數g(t)\in L^2(R)經過scale、translation和modulation產生，令g(t)\in O(\frac{1}{t^2+1})為一個實數的連續可微函數，且限制\|g\|=1、g(t)的積分不為零、g(0)\neq0。

以\gamma=(s, u, \xi)表示scale參數s (s>0)、translation參數u和modulation參數\xi，定義g_\gamma(t)為

g_\gamma(t)=\frac{1}{\sqrt{s}}g(\frac{t-u}{s})e^{i\xi t}

其中\gamma是集合\Gamma=R^+\times R^2中的元素，\frac{1}{\sqrt{s}}使得\|g\|=1。

事實上，函數群D={(g_\gamma(t))}_{\gamma\in\Gamma}含有許多冗餘的元素，對於任何函數f(t)，更有效率的表示方法是，在原子{(g_{\gamma_n}(t))}_{n\in N}中，只選取適當數量的子集合，其中\gamma_n=(s_n, u_n, \xi_n)，則f(t)可以表示為

f(t)=\sum_{n=-\infty}^{+\infty}a_ng_{\gamma_n}(t)

在窗函數傅立葉轉換中，所有原子g_{\gamma_n}具有相同的scale參數s_n=s_0，因此主要分布在一個大小為s_0倍數的區間內，由於上述特性，窗函數傅立葉轉換無法準確地描述比s_0大許多或小許多的函數結構。

小波轉換將信號分解為不同尺度的時頻原子，稱為小波(wavelet)，小波群{(g_{\gamma_n}(t))}_{n\in N}的建構方法是令\xi_n=\xi_0/s_n ，其中\xi_0是一個常數。小波轉換可以分析不同尺寸的信號成分，然而，受限於參數\xi_n和s_n必須成反比的條件，小波轉換的係數無法精準估計傅立葉轉換後具有良好局部性的頻率成分。

## 希爾伯特空間(Hilbert space)的匹配追蹤
自適應時頻分解(adaptive time-frequency decomposition)的目的是將信號展開到一組波形(waveform)上，這些波形選自一個數量龐大的冗餘字典，而匹配追蹤是能達到自適應分解的一種方法。

一個希爾伯特空間可表示為L^2(R)，其組成的複數函數f必須滿足

\|f\|=\int_{-\infty}^{+\infty}^2dt

令H代表一個希爾伯特空間，則將「字典」定義為H中的一個向量群D={(g_\gamma)}_{\gamma\in\Gamma}，滿足\|g_\gamma\|=1，其中\gamma是集合\Gamma=R^+\times R^2中的元素。V代表字典向量的封閉線性生成空間(closed linear span)，在空間V中，集合D之向量的有限線性展開(finite linear expansion)是稠密(dense)的，如果V=H，則稱此字典具有完備性(completeness)。對於「時頻原子分解」段落所描述的字典，H=L^2(R)，在空間L^2(R)中，時頻分子的有限線性展開是稠密的，因此該字典具有完備性。

假設有一信號f\in H，欲將其線性展開到由集合D中選出的一組向量上，使得結果最匹配原來的信號結構。匹配追蹤的方法是連續地將f以其在集合D中元素的正交投影(orthogonal projection)近似。

令g_{\gamma_0}\in D，向量f可以被分解為

f=\langle f, g_{\gamma_0}\rangle g_{\gamma_0}+Rf

其中Rf是將f以g_{\gamma_0}的方向近似後的剩餘向量(residual vector)，由於g_{\gamma_0}和Rf正交，可得下式

{\|f\|}^2=^2+{\|Rf\|}^2

為了最小化\|Rf\|，必須選取g_{\gamma_0}\in D使得|\langle f, g_{\gamma_0}\rangle|最大化。在某些情況下，只能找到近似最佳的向量g_{\gamma_0}，符合

|\langle f, g_{\gamma_0}\rangle|\geq\alpha\ \underset{\gamma\in\Gamma}{\sup}|\langle f, g_\gamma\rangle|

其中0，在選擇向量g_{\gamma_0}時，並非隨機選擇，而是由一個選擇函數C決定。

重複上述步驟，疊代地將剩餘向量Rf投影到集合D中最匹配Rf的向量，並將Rf分解。

匹配追蹤的步驟可以由數學歸納法來表示
# 令R^0f=f
# 假設已經計算第n次的剩餘向量R^nf，n\geq0
# 根據選擇函數C，選取一個最匹配R^nf的元素g_{\gamma_n}\in D
|\langle R^nf, g_{\gamma_n}\rangle|\geq\alpha\ \underset{\gamma\in\Gamma}{\sup}|\langle R^nf, g_\gamma\rangle|剩餘向量R^nf被分解為R^nf=\langle R^nf, g_{\gamma_n}\rangle g_{\gamma_n}+R^{n+1}f由於g_{\gamma_n}和R^{n+1}f正交，可得下式{\|R^nf\|}^2=^2+{\|R^{n+1}f\|}^2當分解到第m次時，f被分解為\begin{align} f & = \sum_{n=0}^{m-1}(R^nf-R^{n+1}f)+R^mf \\ & = \sum_{n=0}^{m-1}\langle R^nf, g_{\gamma_n}\rangle g_{\gamma_n}+R^mf \\ \end{align}{\|f\|}^2被分解為\begin{align} {\|f\|}^2 & = \sum_{n=0}^{m-1}({\|R^nf\|}^2-{\|R^{n+1}f\|}^2)+{\|R^mf\|}^2 \\ & = \sum_{n=0}^{m-1}^2+{\|R^mf\|}^2 \\ \end{align}此公式具有能量守恆的意義，原來的向量f被分解為許多字典中元素的總和。

## 有限空間的匹配追蹤
當信號存在的空間H具有有限的維度N時，匹配追蹤方法會有特殊的特性。在字典D中，可能含有無限多的元素，假設此字典具有完備性，此時可以用一種有效率的匹配追蹤方法，剩餘向量的範數(norm)會以指數方式下降。

當字典含有非常多冗餘的元素時，要尋找和剩餘向量最匹配的向量，通常可以只限制在一個子字典D_\alpha={(g_\gamma)}_{\gamma\in{\Gamma_\alpha}}\subset D中尋找，假設\Gamma_\alpha是一個包含於\Gamma的有限索引集，使得對於所有信號f\in H，滿足

\underset{\gamma\in\Gamma_\alpha}{\sup}|\langle f, g_\gamma\rangle|\geq\alpha\ \underset{\gamma\in\Gamma}{\sup}|\langle f, g_\gamma\rangle|

依據\Gamma_\alpha的大小和字典D的冗餘程度，集合\Gamma_\alpha可以比\Gamma小許多。

以數學歸納法表示此處的匹配追蹤方法
# 計算內積{(\langle f, g_{\gamma}\rangle)}_{\gamma\in\Gamma_\alpha}
# 假設已經計算{(\langle R^nf, g_{\gamma}\rangle)}_{\gamma\in\Gamma_\alpha}，n\geq 0
# 從子字典D_\alpha中找出一個元素g_{\tilde{\gamma_n}}，使得
|\langle R^nf, g_\tilde{\gamma_n}\rangle|=\underset{\gamma\in\Gamma_\alpha}{\sup}|\langle R^nf, g_\gamma\rangle|為了從字典中找到一個比g_{\tilde{\gamma_n}}更匹配f的元素，可以利用牛頓法(Newton’s method)，在\Gamma中尋找\tilde{\gamma_n}的鄰近索引\gamma_n，使得內積達到局部最大值，在此情況下，可以得出下式|\langle R^nf, g_{\gamma_n}\rangle|\geq|\langle R^nf, g_\tilde{\gamma_n}\rangle|\geq \alpha\ \underset{\gamma\in\Gamma}{\sup}|\langle R^nf, g_\gamma\rangle|在此處，選擇函數C與希爾伯特空間中的匹配追蹤不同，必須進行二次搜尋。在選出一個g_{\tilde{\gamma_n}}後，必須計算新的剩餘向量R^{n+1}f和任何g_\gamma\in D_\alpha的內積，更新公式如下\langle R^{n+1}f, g_\gamma\rangle=\langle R^nf, g_\gamma\rangle-\langle R^nf, g_{\gamma_n}\rangle\langle g_{\gamma_n}, g_\gamma\rangle由於先前的計算已經得到\langle R^nf, g_\gamma\rangle和\langle R^nf, g_{\gamma_n}\rangle，因此上式的更新只需要計算\langle g_{\gamma_n}, g_\gamma\rangle。對於一個給定的信號f，要對其剩餘向量分解多少次，決定於要求的精準度\epsilon，重複的次數為能夠滿足下式的最小值p\|R^pf\|=\|f-\sum_{n=0}^{p-1}\langle R^nf, g_{\gamma_n}\rangle g_{\gamma_n}\|\leq\epsilon\|f\|根據能量守恆，此公式等價於\|f\|-\sum_{n=0}^{p-1}^2\leq\epsilon^2\|f\|由於在此方法中，每一次重覆計算時，並沒有計算剩餘向量R^nf，因此只根據上式來判斷是否達到停止分解的條件。

重複次數p決定於\|R^nf\|的下降速率，依據信號的不同，p可以有很大的變化，但在一般情況下，p會比空間H的維度N小很多。

## 性質

* 任何訊號f都會在由原子庫所張的空間中找到收斂的解。
* 稀疏性：當原子庫很大的時候，MP演算法找出來的最佳吻合解，其中的大部分原子訊號的係數可能都是0，只有少部分的係數不為0，此性質稱為稀疏代表性，而此特性對於影像或視訊編碼和壓縮很有幫助。

## 應用

匹配追蹤演算法的靈活性和效率在訊號處理領域中越來越重要，尤其在以下幾種領域中更有其重要的應用：

在視訊編碼和影像壓縮上，對於運動的影像估計和補償，在提出新的原子庫或是擴展的演算法之後，有相當的改良。在影像辨識和形狀辨認上，匹配追蹤演算法的稀疏性對於同樣具有稀疏性的圖像提供新的研究方向。另外在音樂、語音方面，最早即在時頻分析上作為MP演算法研究對象。

## 參見
壓縮感知

## 參考文獻
[1] S. G. Mallat and Z. Zhang, "Matching pursuits with time-frequency dictionaries," in IEEE Transactions on Signal Processing, vol. 41, no. 12, pp. 3397-3415, Dec 1993.

[2] Jian-Jiun Ding, Time frequency analysis and wavelet transform class note,the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2012.

[3] B. Torresani, "Wavelets associated with representations of the affine Weyl-Heisenberg group," 1. Math. Physics, vol. 32, pp. 1273-1279, May 1991.