Orthogonal Matching Pursuit(omp)

17 Aug 2014

Recently, I have been learning to use omp methods to find sparse representation of the new data via an over complete dictionary. I want to put my summary here as a look-up page for my futher research. This is an incomplete version, since I also want to include several improved way such as Cholesky, QR, Batch and some example code. I will add some reference papers and webpages to those material later.

• Sparse Learning
• OMP Approach
  • Naive OMP
  • Cholesky way OMP
  • QR way OMP
  • Batch OMP

Sparse Learning

We are interested in modeling a column signal $\underline{x}_{m\times 1}$ (for example a new data observation) by a linear combination of atoms (columns) selected from a dictionary $\mathbf{D_{m\times n}}$ where $m$ is number of features and $n$ for # of observations. We want to find a sparse representation $\underline{\gamma}_{n\times 1}$ of $\underline{x}_{m\times 1}$ such that the difference between $\mathbf{D}\underline{\gamma}$ and $\underline{x}$ becomes negligible. Formally, the problem can be described by the following optimization problem:

---

$$\begin{equation} \begin{aligned} & \hat{\gamma}=\underset{\underline{\gamma}}{argmin}\left\Vert \underline{\gamma}\right\Vert _{0} & \text{subject to} && \left\Vert \underline{x} - \mathbf{D}\underline{r}\right\Vert^{2}_{2} \leq \epsilon \end{aligned} \end{equation}$$

---

In formula $\underline{\gamma}$ is the sparse representation of $\underline{x}$, $\epsilon$ is the error tolerance, and $\left\Vert \cdotp \right\Vert_{0}$ is a pseudo-norm which counts the number of non-zero entries. However, the above optimization problem is NP-hard which can be efficiently solved using several approximation methods. And OMP (Orthogonal Matching Pursuit) is one of them.

OMP Approach

OMP is an approximation to the above optimization by solving one of the following two problems:

---

sparsity-constrained sparse coding problems given by

$$\begin{equation} \begin{array} & \underline{\hat{\gamma}} = \underset{\underline{\gamma}}{argmin} \left\Vert \underline{x} - \mathbf{D}\underline{\gamma} \right\Vert_{2}^{2} & \text{subject to} && \left\Vert\underline{\gamma}\right\Vert_{0} \leq K \end{array} \end{equation}$$

error-constrained sparse coding problem given by

$$\begin{equation} \begin{array} & \underline{\hat{\gamma}} = \underset{\underline{\gamma}}{argmin}\left\Vert\underline{\gamma}\right\Vert_{0} & \text{subject to} && \left\Vert\underline{x} - \mathbf{D}\underline{\gamma}\right\Vert_{2}^{2} \leq \epsilon \end{array} \end{equation}$$

---

Note that the columns of dictionary $\mathbf{D_{m\times n}}$ should be normalized to unit $\mathcal{L}^{2}$-length.

Naive OMP

The OMP use method conceptually similar to Gram-Schmidt, since the Gram-Schmidt can also find up to $m$ orthogonal vectors to express the new observations $\underline{x}$. You see that for a $m$-dimensional vector space, you can find at most $m$ orthogonal directions.

The greedy OMP method selects each column atom from dictionary $\mathbf{D_{m\times n}}$ with the highest correlation to the current residual. Once the atom is selected, the signal $\underline{x}$ is orthogonally projected to the span of all selected atoms (means all linear combinations of the selected atoms set), the residual is computed, and process repeats.

The detail algorithm is the table below:

Step	Description
1	Input: Dictionary $\mathbf{D_{m \times n}}$, signal $\underline{x}_{m\times 1}$, target sparsity $K$ or target error $\epsilon$
2	Ouput: Sparse representation $\underline{\gamma}$ such that $\underline{x} \approx \mathbf{D}\underline{\gamma}$
3	Initialize: Set $I:=()$, $\underline{r}:=\underline{x}$, $\underline{\gamma}:=\underline{0}$
4	while (stopping criterion not meet) do
5	$\hat{k}:=\underset{k}{argmax}\left|\underline{d}_{k}^{T}\underline{r}\right|$
6	$I:=(I,\hat{k})$
7	$\underline{\gamma}_{I}:=(\mathbf{D}_{I})^{\dagger}\underline{x}$
8	$\underline{\gamma}:= \underline{x} - \mathbf{D}_{I} \underline{\gamma}_{I}$
9	end while

where $\underline{d}_{k}$ is one column element from $\mathbf{D_{m\times n}}=\left[\underline{d}_{1}\mid\underline{d}_{2}\mid\cdots\mid\underline{d}_{n}\right]$, $I$ is an ordered list of selected column indexes, $\mathbf{D}_{I}$ and $\underline{\gamma}$ are columns or entries selected by index $I$.

Cholesky way OMP

The naive OMP requires pseudo inverse of $\mathbf{D}_{I}^{\dagger} = (\mathbf{D}_{I}^{T}\mathbf{D}_{I})^{-1}\mathbf{D}_{I}^{T}$ which brings high computation cost with each iteration. However, $\mathbf{D}_{I}^{T}\mathbf{D}_{I}$ is updated by simply adding a single row an column to it every iteration (details below). We can use methods such as Cholesky or QR which only requires the update of its last row.

---

Let us denote $\mathbf{D}_{I}^{T}\mathbf{D}_{I}$ as $A$. $A$ is a symmetric positive definite matrix and can be used to decompose into the form of $A=LL^{*}$. $L$ is a lower triangular matrix and $L^{*}$ denotes the conjugate transpose of $L$. For iteration $(n-1)$ th, we have: (again $\mathbf{D}_{I}$ denotes a subdictionary of $\mathbf{D}$ selected by ordered set $I$ and $I_{n}$ is the ordered set in the $(n)$th iteration).

$$\begin{equation} \mathbf{D}_{I_{n-1}}^{T}\mathbf{D}_{I_{n-1}} = A_{n-1}=L_{n-1}L_{n-1}^{*}\in\mathbb{R}^{(n-1)\times(n-1)} \end{equation}$$

After selecting the maximum $\hat{k}_{n-1}$ in step 5, $\mathbf{D}_{I_{n}}$ adds a new column $\mathbf{D}_{I_{n}}=\left[\mathbf{D}_{I_{n-1}} \mid \mathbf{D}_{\{\hat{k}_{n-1 }\}} \right]$ and $A_{n}$ becomes

$$\begin{equation} A_{n}=\left[ \begin{array}{cc} A_{n-1} & \underline{v}\\ \underline{v}^{T} & c \end{array} \right] =L_{n}L_{n}^{*} \end{equation}\in\mathbb{R}^{(n)\times(n)}$$

where each entry in $\underline{v}$ is the dot product between every column of $\mathbf{D}_{I_{n-1}}$ and $\mathbf{D}_{\{\hat{k}_{n-1 }\}}$ and $c$ is the self dot product of $\mathbf{D}_{\{\hat{k}_{n-1 }\}}$ or $\mathbf{D}_{\{\hat{k}_{n-1 }\}}^{T}\mathbf{D}_{\{\hat{k}_{n-1 }\}}$.

$L_{n}$ is updated by

$$\begin{equation} \begin{array}{ccc} L_{n}=\left[\begin{array}{cc} L_{n-1} & \underline{0}\\ \underline{\omega}^{T} & \sqrt{c-\underline{\omega}^{T}\underline{\omega}} \end{array}\right] & \text{,} & \underline{\omega}=L_{n-1}^{-1}\underline{v}\end{array} \end{equation}$$

QR way OMP

Batch OMP