# Orthogonal Matching Pursuit with the Matrix Inversion Lemma

Continuing from the naive OMP, the Cholesky OMP, and the QR OMP, we now look at using the matrix inversion lemma (MIL) to make OMP extremely efficient. This is one I haven’t ever seen, but it might be the same as the block inversion approach discussed by Pati et al. 1993. I still have to take a closer look.

Note: If you see some LaTeX error messages about macros not being defined, shift reload the page.

Consider that in its $$k$$th iteration, OMP has
the set of dictionary inner products
$$\Gamma_l = \{\langle \vphi_l, \vphi_n \rangle\}_{n \in \Omega}$$
and squared norms $$\{\|\vphi_n\|^2\}_{n \in \Omega}$$,
the sets of original and current inner products,
$$\mathcal{I}_{0} = \{\langle \vu, \vphi_n \rangle \}_{n \in \Omega}$$ and
$$\mathcal{I}_{k} = \{\langle \vr_{k}, \vphi_n \rangle \}_{n \in \Omega}$$,
the current solution $$\vx_k$$,
and the basis dual to $$\MPhi_{\Omega_k}$$, i.e.,
$$\MPsi_{\Omega_k} = \MPhi_{\Omega_k}(\MPhi_{\Omega_k}^H\MPhi_{\Omega_k})^{-1}$$.
OMP guarantees that $$\MPhi_{\Omega_k}$$
has full column rank when $$k \le M$$,
and so $$\MPsi_{\Omega_k}$$ always exists when $$k \le M$$.

When it selects the new atom index $$n_{k+1}$$ from $$\mathcal{I}_{k}$$,
OMP must update the dual basis and solution.
With the matrix inversion lemma, we see
$$\begin{multline} (\MPhi_{\Omega_{k+1}}^H\MPhi_{\Omega_{k+1}})^{-1} = \left [ \begin{array}{cc} \MPhi_{\Omega_k}^H\MPhi_{\Omega_k} & \MPhi_{\Omega_k}^H\vphi_{n_{k+1}} \\ \vphi_{n_{k+1}}^H \MPhi_{\Omega_k} & \|\vphi_{n_{k+1}}\|^2 \end{array} \right ]^{-1} \\ = \left [ \begin{array}{cc} (\MPhi_{\Omega_k}^H\MPhi_{\Omega_k})^{-1} + \lambda \vh \vh^H & -\lambda \vh\\ -\lambda \vh^H & \lambda \end{array} \right ] \end{multline}$$
where $$\vh = \MPsi_{\Omega_k}^H\vphi_{n_{k+1}}$$, and
$$1/\lambda = \|\vphi_{n_{k+1}}\|^2 – \|\MPhi_{\Omega_k}\vh\|^2$$.
OMP thus computes the new dual basis by
$$\begin{multline} \MPsi_{\Omega_{k+1}} = \MPhi_{\Omega_{k+1}}(\MPhi_{\Omega_{k+1}}^H\MPhi_{\Omega_{k+1}})^{-1} \\ = \left [ \MPhi_{\Omega_{k}} | \vphi_{n_{k+1}} \right ] \left [ \begin{array}{cc} (\MPhi_{\Omega_k}^H\MPhi_{\Omega_k})^{-1} + (\lambda \vh) \vh^H & -\lambda \vh\\ -\lambda \vh^H & \lambda \end{array} \right ] \\ = \left [ \MPsi_{\Omega_k} – (\lambda\vv)\vh^H | \lambda\vv \right ] \end{multline}$$
where $$\vv = \vphi_{n_{k+1}} – \MPhi_{\Omega_{k}} \vh$$.
To update the dual basis,
OMP only computes $$\vh$$ and $$\lambda$$
at a complexity of $$\mathcal{O}(N + Mk)$$, including the search.
The new solution is then given by
$$\vx_{k+1} = \MPsi_{\Omega_{k+1}}^H \vu = \left [ \begin{array}{c} \vx_k \\ 0\end{array} \right ] – \Delta \left [ \begin{array}{c} \vh \\ -1 \end{array} \right ]$$
where
\begin{align} \Delta & = \lambda \vv^H\vu = \lambda (\vphi_{n_{k+1}} – \MPhi_{\Omega_{k}} \vh)^H\vu \\ & = \lambda [\langle \vu, \vphi_{n_{k+1}} \rangle – \langle\vu, \MPhi_{\Omega_{k}} \vh \rangle] \\ & = \lambda \left [\mathcal{I}_{0}(n_{k+1}) – \sum_{l=1}^k h_l^* \mathcal{I}_{0}(n_{l}) \right ]. \end{align}
OMP can update the set of projections without computing the solution:
\begin{align} \mathcal{I}_{k+1} & = \{\langle \vr_{k+1}, \vphi_n \rangle \}_{n \in \Omega} = \{\langle \vu – \MPhi_{\Omega_{k+1}}\vx_{k+1}, \vphi_n \rangle \}_{n \in \Omega} \\ & = \left \{\langle \vu – \MPhi_{\Omega_{k+1}} \left (\left [ \begin{array}{c} \vx_k \\ 0\end{array} \right ] – \Delta \left [ \begin{array}{c} \vh \\ -1 \end{array} \right ] \right ), \vphi_n \rangle \right \}_{n \in \Omega} \\ & = \left \{ \mathcal{I}_{k}(n) + \Delta \langle \MPhi_{\Omega_{k}}\vh – \vphi_{n_{k+1}}, \vphi_n \rangle \right \}_{n \in \Omega} \\ & = \left \{ \mathcal{I}_{k}(n) – \Delta\langle \vphi_{n_{k+1}}, \vphi_n \rangle + \Delta \sum_{l=1}^k h_l \Gamma_{n_l}(n) \right \}_{n \in \Omega}. \end{align}
At the same complexity, OMP can also use
$$\begin{multline} \mathcal{I}_{k+1} = \left \{\langle \vr_{k+1}, \vphi_n \rangle\right\}_{n \in \Omega} = \left \{ \langle \vu,\vphi_n \rangle – \langle\MPhi_{\Omega_{k+1}} \vx_{k+1}, \vphi_n \rangle \right\}_{n \in \Omega} \\ = \left \{\mathcal{I}_{0}(n) – \sum_{l=1}^{k+1} x_{l} \langle\vphi_{n_l}, \vphi_n \rangle\right\}_{n \in \Omega} \\ = \left \{\mathcal{I}_{0}(n) – \sum_{l=1}^{k+1} x_{l} \Gamma_{n_l}(n) \right \}_{n \in \Omega} \end{multline}.$$
OMP can easily update the residual energy.
First, we see
\begin{align} \|\vr_{k+1}\|^2 & = \|\vu – \MPhi_{\Omega_{k+1}}\vx_{k+1} \|^2 \\ & = \left \| \vu – \MPhi_{\Omega_{k+1}} \left ( \left [ \begin{array}{c} \vx_k \\ 0\end{array} \right ] – \Delta \left [ \begin{array}{c} \vh \\ -1 \end{array} \right ]\right ) \right \| \\ & = \left \|\vr_k + \Delta \MPhi_{\Omega_{k+1}}\left [ \begin{array}{c} \vh \\ -1 \end{array} \right ] \right \|^2 \\ & = \|\vr_k\|^2 + \left \| \Delta \MPhi_{\Omega_{k+1}}\left [ \begin{array}{c} \MPsi_{\Omega_k}^H\vphi_{n_{k+1}} \\ -1 \end{array} \right ]\right \|^2 + \\ & \qquad\qquad 2\text{Re}\left \{ \Delta \vr_k^H \MPhi_{\Omega_{k+1}}\left [ \begin{array}{c} \MPsi_{\Omega_k}^H\vphi_{n_{k+1}} \\ -1 \end{array} \right ]\right \}. \end{align}
Since $$\vr_k$$ is orthogonal to all selected atoms except the last,
$$\vr_k^H \MPhi_{\Omega_{k+1}} = [\zerob^T \mathcal{I}_k(n_{k+1})]$$,
and so this finally becomes
\begin{align} \|\vr_{k+1}\|^2 & = \|\vr_k\|^2 + |\Delta|^2 \| \vv \|^2 – 2\text{Re} \{ \Delta \mathcal{I}_k(n_{k+1}) \} \end{align}
since $$\vv = \vphi_{n_k} – \MPhi_{\Omega_{k}} \MPsi_{\Omega_k}^H\vphi_{n_{k+1}}$$.
In total, OMP with the matrix inversion lemma form,
has complexity of $$\mathcal{O}((N + M)k)$$.

Here is my code: ompdMIL.m