# Orthogonal Matching Pursuit the Cholesky Way

Continuing from my previous post, we now look at using Cholesky decomposition to make OMP extremely efficient.

Consider that in its $$k$$th step, OMP has the two sets of projections
$$\mathcal{I}_{k-1} = \{\langle \vr_{k-1}, \vphi_n \rangle/\|\vphi_n\| \}_{n \in \Omega}$$
and $$\mathcal{I}_{0} = \{\langle \vu, \vphi_n \rangle/\|\vphi_n\|\}_{n \in \Omega}$$,
the set of dictionary inner products
$$\Gamma_l = \{\langle \vphi_l, \vphi_n \rangle/\|\vphi_n\|\}_{n \in \Omega}$$
and norms $$\{\|\vphi_n\|\}_{n \in \Omega}$$,
and the Cholesky decomposition of the Gramian of the
subdictionary $$\MPhi_{\Omega_{k-1}}^H \MPhi_{\Omega_{k-1}} = \ML_{k-1}\ML_{k-1}^H$$.
Every Gramian can be expressed as a Cholesky decomposition.

OMP selects the new atom index $$n_k$$ by searching $$\mathcal{I}_{k-1}$$,
and updates the Cholesky factorization by
first solving for $$\vv$$ in the triangular system
$$\ML_{k-1}\vv = \MPhi_{\Omega_{k-1}}^H\vphi_{n_k}$$
(the righthand side coming from $$\Gamma_{n_k}$$),
and then updating the triangular matrix by
$$\ML_{k} = \left [ \begin{array}{cc} \ML_{k-1} & \zerob \\ \vv^H & \sqrt{\|\vphi_{n_k}\|^2 – \|\vv\|^2} \end{array} \right ].$$
The complexity of the search and update is $$\mathcal{O}(N + k^2 + k)$$.
The minimum $$\ell_2$$-norm of the solution error satisfies
$$\MPhi_{\Omega_{k}}^H\MPhi_{\Omega_{k}} \vx_k = \ML_{k}\ML_{k}^H \vx_k = \MPhi_{\Omega_{k}}^H \vu.$$
where the righthand side comes from $$\mathcal{I}_{0}$$.
OMP only needs to solve two triangular systems by backsubstitution
\begin{align} \ML_k \vy & = \MPhi_{\Omega_k}^H \vu \\ \ML_k^H \vx_k & = \vy \end{align}
which has complexity $$\mathcal{O}(k^2)$$.
The new residual energy is just
\begin{align} \|\vr_k\|^2 & = \|\vu – \MPhi_{\Omega_k} \vx_k\|^2 \\ & = \|\vu\|^2 + \|\MPhi_{\Omega_k} \vx_k\|^2 – 2\text{Re}\{\vx_k^H \MPhi_{\Omega_k}^H\vu\} \\ & = \|\vu\|^2 + \vx_k^H \ML_{k}\ML_{k}^H \vx_k – 2\text{Re}\{\vx_k^H \ML_{k} \vy\} \\ & = \|\vu\|^2 – \|\vy\|^2. \end{align}
Finally, OMP updates the inner products in one of two ways.
For the first,
OMP explicitly computes the residual signal
$$\vr_k = \vu – \MPhi_{\Omega_k} \vx_k$$
with a cost of $$\mathcal{O}(Mk)$$,
and then directly compute the inner products
$$\mathcal{I}_{k} = \{\langle \vr_k, \vphi_n\rangle/\|\vphi_n\|\}_{n \in \Omega}$$
with a cost of $$\mathcal{O}(MN)$$.
In the second case, OMP computes
$$\begin{multline} \mathcal{I}_{k} = \left \{\frac{\langle \vr_{k}, \vphi_n \rangle}{\|\vphi_n\|} \right\}_{n \in \Omega} = \left \{\frac{\langle \vu,\vphi_n \rangle}{\|\vphi_n\|} – \frac{\langle\MPhi_{\Omega_k} \vx_k, \vphi_n \rangle}{\|\vphi_n\|} \right\}_{n \in \Omega} \\ = \left \{\mathcal{I}_{0}(n) – \sum_{l=1}^k x_{l} \frac{\langle\vphi_{n_l}, \vphi_n \rangle}{\|\vphi_n\|} \right\}_{n \in \Omega} \\ = \left \{\mathcal{I}_{0}(n) – \sum_{l=1}^k x_{l} \Gamma_{n_l}(n) \right \}_{n \in \Omega} \end{multline}$$
where $$x_l$$ is the $$l$$th element of $$\vx_k$$.
This entire update has a complexity of $$\mathcal{O}(Nk)$$.
Thus, the total complexity of its $$k$$th iteration is either
$$\mathcal{O}(MN + k^2 + Mk)$$ in the first case,
or $$\mathcal{O}(Nk + k^2)$$ by using the latter approach.

Here is my MATLAB implementation of the Cholesky OMP: ompdChol.m