*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