# Orthogonal Matching Pursuit the QR Way

Continuing from the naive OMP and the Cholesky OMP, we now look at using QR decomposition to make OMP extremely efficient. I previously discussed this approach here, but below I make some additions.

Consider that in the $$k$$th step of OMP we have
$$\mathcal{I}_{0} = \{\langle \vu, \vphi_n \rangle/\|\vphi_n\| \}_{n \in \Omega}$$,
$$\mathcal{I}_{k-1} = \{\langle \vr_{k-1}, \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}$$,
as well as the QR decomposition of the subdictionary produced thus far,
i.e., $$\MPhi_{\Omega_{k-1}} = \MQ_{k-1}\MR_{k-1}$$.
The QR decomposition of a matrix will always exist,
but is not unique.

OMP searches through $$\mathcal{I}_{k-1}$$
to find a new atom index $$n_k$$.
It then defines $$\vw = \MQ_{k-1}^H\vphi_{n_k}$$
and updates the QR decomposition by
\begin{align} \MR_k & = \left [ \begin{array}{cc} \MR_{k-1} & \vw \\ \zerob^T & \sqrt{\|\vphi_{n_k}\|^2 – \|\vw\|^2} \end{array} \right ] \\ \MQ_{k} & = \left [ \begin{array}{cc} \MQ_{{k-1}} & \frac{\vphi_{n_k} – \MQ_{{k-1}}\vw}{ \sqrt{\|\vphi_{n_k}\|^2 – \|\vw\|^2}} \end{array} \right ]. \end{align}
With the search, updating these matrices has complexity
$$\mathcal{O}(N + Mk)$$.
OMP computes the residual energy recursively
\begin{align} \|\vr_k\|^2 & = \|\vu – \MQ_{{k-1}}\MQ_{{k-1}}^H \vu – \vq_k \vq_k^H\vu\|^2 \\ & = \|\vr_{k-1}\|^2 – |\langle \vu,\vq_k\rangle|^2 \end{align}
where $$\vq_k$$ is the last column of
the orthonormal basis in $$\MQ_k$$.
OMP can then update the projections
without computing the residual since
$$\vr_k = \vr_{k-1} – \vq_k \vq_k^H\vu$$:
$$\begin{multline} \mathcal{I}_{k} = \left \{\frac{\langle \vr_{k}, \vphi_n \rangle}{\|\vphi_n\|} \right\}_{n \in \Omega} = \left \{\frac{\langle \vr_{k-1},\vphi_n \rangle}{\|\vphi_n\|} – \frac{\langle\vq_k \vq_k^H\vu, \vphi_n \rangle}{\|\vphi_n\|} \right\}_{n \in \Omega} \\ = \left \{\mathcal{I}_{k-1}(n) – \langle \vu,\vq_k \rangle\frac{\langle\vq_k, \vphi_n \rangle}{\|\vphi_n\|} \right\}_{n \in \Omega}. \end{multline}$$
This has a complexity of $$\mathcal{O}(NM)$$.

Once finished, OMP finds the solution at iteration $$K$$ by
solving two triangular systems
\begin{align} \MR_K^H \vy & = \MPhi_{\Omega_K}^H \vu \\ \MR_K \vx_K & = \vy \end{align}
where the right hand side of the first system comes from $$\mathcal{I}_{0}$$.
Solving both of these with back-substitution
has complexity $$\mathcal{O}(K^2)$$.
In this way, the total complexity of OMP in its $$k$$th iteration
is $$\mathcal{O}(NM + Mk)$$,
which is linear in $$k$$.
Notice that in this implementation
OMP does not compute the solution or the residual
during decomposition.

If instead, OMP updates the set of projections using
$$\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}$$
then it might reduce the complexity at the cost of
computing the solution at each iteration.
After computing the solution at complexity $$\mathcal{O}(k^2)$$,
OMP can update the projections using the above,
and update the residual energy using
\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}
Thus, using this approach, the total complexity of its $$k$$th iteration is
$$\mathcal{O}((N+M)k + k^2)$$.

Here is my code: ompdQR.m