Sparse Linear Regression Via Generalized Orthogonal Least-Squares
The Orthogonal Least Squares (OLS) algorithm sequentially selects columns of the coefficient matrix to greedily find an approximate sparse solution to an underdetermined system of linear equations. In this paper, conditions under which OLS recovers sparse signals from a low number of random linear measurements with probability arbitrarily close to one are stated. Moreover, a computationally efficient generalization of Orthogonal Least-Squares which relies on a recursive relation between the components of the optimal solution to select $L$ columns at each step and solve the resulting overdetermined system of equations is proposed. This generalized OLS algorithm is empirically shown to outperform existing greedy algorithms broadly used in literature. % such as orthogonal Matching pursuits (OMP), OLS, and thus reduces the performance gap between the class of greedy algorithms and convex optimization methods for sparse linear regression. %other greedy solvers from literature.
main.tex
\documentclass{article}
\usepackage{spconf,amsmath,graphicx,multicol}
\usepackage{amsfonts}
\usepackage{bbm}
\usepackage{cite}
\usepackage{amssymb}
\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{tabularx}
\usepackage[utf8]{inputenc}
\usepackage[english]{babel}
\usepackage{amsthm}
\usepackage{subcaption}
\theoremstyle{definition}
\newtheorem{definition}{Definition}[section]
\newtheorem{proposition}{Proposition}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{corollary}{Corollary}[theorem]
\newtheorem{lemma}[theorem]{Lemma}
\let\oldemptyset\emptyset
\let\emptyset\varnothing
\newcommand\norm[1]{\left\lVert#1\right\rVert}
\newcommand{\argmax}{\arg\!\max}
\newcommand{\argmin}{\arg\!\min}
% \captionsetup[subfigure]{labelformat=simple,labelsep=colon}
% \renewcommand{\thesubfigure}{\arabic{subfigure}}

\def\x{{\mathbf x}}
\def\s{{\mathbf s}}
\def\y{{\mathbf y}}
\def\e{{\mathbf e}}
\def\r{{\mathbf r}}
\def\A{{\mathbf A}}
\def\B{{\mathbf B}}
\def\P{{\mathbf P}}
\def\L{{\cal L}}
\def\R{{\mathbb{R}}}
\def\N{{\mathbb{N}}}
\def\I{{\mathbf I}}
\def\a{{\mathbf a}}
\def\r{{\mathbf r}}
\usepackage{mathtools}
\DeclarePairedDelimiter{\ceil}{\lceil}{\rceil}



\title{Sparse Linear Regression Via Generalized Orthogonal Least-Squares}
%
% Single address.
% ---------------
\name{Abolfazl Hashemi and Haris Vikalo}
\address{Department of Electrical and Computer Engineering\\
		University of Texas at Austin, Austin, TX, USA}
%
% For example:
% ------------
%\address{School\\
%	Department\\
%	Address}
%
% Two addresses (uncomment and modify for two-address case).
% ----------------------------------------------------------
%\twoauthors
%  {A. Author-one, B. Author-two\sthanks{Thanks to XYZ agency for funding.}}
%	{School A-B\\
%	Department A-B\\
%	Address A-B}
%  {C. Author-three, D. Author-four\sthanks{The fourth author performed the work
%	while at ...}}
%	{School C-D\\
%	Department C-D\\
%	Address C-D}
%
\begin{document}
%\ninept
%
\maketitle
%
\begin{abstract}
The Orthogonal Least Squares (OLS) algorithm sequentially selects columns of the coefficient 
matrix to greedily find an approximate sparse solution to an underdetermined system of linear equations. 
In this paper, conditions under which OLS recovers sparse signals from a low number of random linear 
measurements with probability arbitrarily close to one are stated. Moreover, a computationally efficient 
generalization of Orthogonal Least-Squares which relies on a recursive relation between the components 
of the optimal solution to select $L$ columns at each step and solve the resulting overdetermined 
system of equations is proposed. This generalized OLS algorithm is empirically shown to outperform 
existing greedy algorithms broadly used in literature. % such as orthogonal Matching pursuits (OMP), OLS, and thus reduces the performance gap between the class of greedy  algorithms  and  convex  optimization  methods  for sparse linear regression.
%other greedy solvers from literature.
\end{abstract}
%
\begin{keywords}
linear regression, compressed sensing, greedy algorithm, orthogonal least-squares
\end{keywords}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{sec:intro}
The problem of finding sparse solution to an underdetermined system of linear equations arises in a number 
of practical scenarios. Examples include compressed sensing \cite{donoho2006compressed}, sparse channel
estimation in communication systems \cite{mitr07,barik2014sparsity}, compressive DNA microarrays 
\cite{parvaresh2008recovering} as well as a number of other applications in signal processing and machine 
learning \cite{lustig2007sparse,elad2010role,mishali2010theory,elhamifar2009sparse}. Consider the linear 
measurement model
\begin{equation} \label{eq:1}
{\y=\A\x+{\bf e}},
\end{equation}
where $\y \in\R^{n}$ denotes the vector of observations, $\A \in\R^{n \times m}$ is the coefficient matrix (i.e., a collection of features) assumed to be full rank, ${\bf e} \in\R^{n}$ is the additive observation noise vector, and $\x \in\R^{m}$ is a vector known to have at most $k$ non-zero components (i.e., $k$ is the sparsity level of $\x$).  We are interested in finding a sparse approximation to $\x$; in particular, we would like to solve the so-called $l_0$-constrained least-squares
\begin{equation}  \label{eq:2}
\begin{aligned}
& \underset{\x}{\text{minimize}}
&  \norm{\y-\A\x}^{2}_{2}
& &\text{subject to}
& & \norm{\x}_{0} \leq k.
\end{aligned}
\end{equation}
%\subsection{Prior Work}
The number of possible locations of non-zero entries in $\x$ scales combinatorially with $n$ which renders
(\ref{eq:2}) computationally challenging; in fact, the problem is NP-hard. To enable computationally efficient 
search for sparse $\x$ approximating (\ref{eq:1}), the non-convex $l_0$-norm-constrained optimization
(\ref{eq:2}) can be replaced by a sparsity-promoting $l_1$-norm optimization
\begin{equation} 
\begin{aligned}
& \underset{\x}{\text{minimize}}
& \norm{\x}_{1}
& & \text{subject to}
& & \norm{\y - \A \x}_2 \leq \varepsilon,
\end{aligned}
\label{eq:L1}
\end{equation}
where $\varepsilon > 0$ is a predetermined measure of noise power. In the noise-free scenario where
$\e$ in (\ref{eq:1}) and $\varepsilon$ in (\ref{eq:L1}) are both zero and where $\A$ satisfies certain properties, 
it is known that a sufficiently sparse $\x$ can be reconstructed exactly \cite{candes2006robust}. 
However, while the convexity of $l_{1}$-norm enables finding 
the optimal solution to the reformulated sparse vector recovery problem, the complexity of doing so (by 
means of, e.g., iterative shrinkage-thresholding algorithms such as \cite{beck2009fast}, or alternating 
direction method of multipliers \cite{boyd2011distributed}) is often prohibitive when one
deals with high dimensional data. For this reason, a number of fast greedy heuristics that attempt to solve 
(\ref{eq:2}) directly by successively identifying columns of $\A$ which correspond to non-zero components 
of $\x$ have been proposed \cite{tropp2004greed,needell2010signal}. Among those,
particular attention has been paid to the orthogonal matching pursuit (OMP) algorithm 
\cite{tropp2007signal,pati1993orthogonal} which has an intuitive geometric interpretation and is characterized 
by high speed and competitive performance; numerous modifications of OMP that explore the trade-off 
between accuracy and speed have been proposed in literature 
\cite{donoho2012sparse,needell2009cosamp}. A related
Orthogonal Least-Squares (OLS) method \cite{chen1989orthogonal}, proposed as an identification 
algorithm for parameter estimation of generally multivariable non-linear systems which are linear in 
parameters, has recently been employed in compressed sensing \cite{soussen2013joint}. In general, OLS 
outperforms OMS in settings where the columns of $\A$ are non-orthogonal but it does so at a moderate
increase in complexity. The existing analysis and performance guarantees for OLS are limited to the
case of non-random measurements \cite{soussen2013joint,herzet2013exact,herzet2016relaxed}.

In this paper, we provide a result establishing that in the noiseless scenario where the coefficient matrix 
is drawn at random from a Gaussian or a Bernoulli distribution, with ${\cal O}\left(k\log(m)\right)$ linear random 
measurements OLS guarantees recovery of $\x$ with high probability. This result is comparable to those previously 
provided for OMP \cite{tropp2007signal,rangan2009orthogonal}. Moreover, we propose a generalization of 
OLS, the Generalized Orthogonal Least-Squares (GOLS), an efficient algorithm which relies on a recursive 
relation between the components of the optimal solution to (\ref{eq:1}) to select a pre-determined number of 
columns and provide performance superior to existing methods.



%The rest of the paper is organized as follows. Section 2 reviews OLS and its connection to OMP, and 
%states the results on sparse recovery via OLS from random measurements. In Section 3 
%we develop the proposed algorithm. Section 4 and Section 5 present simulation results and 
%concluding remarks, respectively.

\section{Performance Guarantee for Orthogonal Least-Squares}

The OLS algorithm sequentially projects columns of $\A$ onto a residual vector and each time selects the 
column that leads to the smallest residual norm. Specifically, OLS chooses a new index $j_s$ as
\begin{equation*}
{j}_{s}=\argmin_{j \in {\cal I}}{\norm{\y - \A_{{\cal S}_{i-1}\cup\{j\}}\A_{{\cal S}_{i-1}\cup\{j\}}^{\dagger}\y}_2},
\end{equation*}
where $\cal I$ is the set of indices that are not yet selected. This procedure is computationally more expensive 
than OMP since in addition to solving the least-square problem to update the residual vector, orthogonal 
projection of each column needs to be found at each step of OLS. Perhaps in part due to this increase in 
complexity, OLS has not played as prominent role in sparse signal recovery literature as OMP did.

Note that the performances of OLS and OMP are identical when the columns of $\A$ are 
orthogonal.\footnote{In fact, orthogonality of the columns of $\A$ leads to a modular objective function in
(\ref{eq:2}), implying optimality of both methods.} It is beneficial to further clarify the difference between 
OMP and OLS. In each iteration of OMP, an element that best correlates with the current residual is chosen. 
OLS, on the other hand, selects a column that has the largest portion which is inexpressible by previously 
selected columns which, in turn, minimizes related approximation error. 

The following theorem states that for Gaussian and Bernoulli matrices with normalized columns, which are
often considered in compressed sensing problems, in the noiseless scenario OLS is with high probability 
capable of the exact recovery of sparse signals if the number of measurements grows linearly with the 
sparsity level and logarithmically with the dimension of the unknown signal.
\begin{theorem}
\textit{Suppose that $\x \in \R^m$ is an arbitrary sparse vector with sparsity level $k$. Consider a random 
matrix $\A \in \R^{n\times m}$ such that its entries are drawn uniformly and independently from either 
${\cal N}(0,1/n)$ or $\{+1\big\slash\sqrt{n},-1\big\slash\sqrt{n}\}$. Given the noiseless observation 
$\y=\A\x$, the OLS algorithm can recover $\x$ in $k$ iterations with probability of success exceeding 
$1-\delta$ if $n={\cal O}\left(k\log(m/\delta)\right)$ for some $0<\delta<\gamma$, where $\gamma$ is a 
positive constant which is independent of $n$, $m$, and $k$.}
\end{theorem}

The proof, which exploits the fact that the columns of A are spherically symmetric random vectors and
relies on Johnson-Lindenstrauss lemma, is omitted for brevity.

%A more detailed geometric comparison is given in \cite{blumensath2007difference}. %The OLS 
%algorithm is formalized as Algorithm 1. 
% \begin{algorithm}
% \vspace{0.1cm}
% \begin{tabularx}{\textwidth}{l>{$}c<{$}X}
% \textbf{Input:}  \hspace{0.85cm}observation vector $\y$ \\ \hspace{1.95cm}coefficient matrix $\A$ \\ \hspace{1.95cm}sparsity level $k$ \vspace{0.1cm}\\
% \textbf{Output:} \hspace{0.5cm} recovered support ${\cal S}_k$\\ \hspace{1.92cm}estimated signal $\hat{\x}_{k}$\vspace{0.1cm}\\
% \textbf{Initialize:} \hspace{0.3cm} ${\cal S}_0=\oldemptyset$, ${\cal I} = \{1,2,\dots,m\}$%, ${\bf r}=\y$, f\left(S_{0}\right)=0$
% \end{tabularx}
% \begin{algorithmic}
% \FOR  {\hspace{0.40cm}$i=1$ to $k$ \hspace{0.40cm}}\vspace{0.1cm}
% \STATE 1. ${j}_{s}=\argmin_{j \in {\cal I}}{\norm{\y - \A_{{\cal S}_{i-1}\cup\{j\}}\A_{{\cal S}_{i-1}\cup\{j\}}^{\dagger}\y}_2}$ \vspace{0.1cm}\\
% 2. ${\cal S}_{i}={\cal S}_{i-1}\cup\{j_{s}\}$ , ${\cal I} = {\cal I} \backslash {\cal S}_{i}$ \vspace{0.1cm}\vspace{0.1cm}\\
% 3. $\hat{\x}_{i}=\A_{{\cal S}_{i}}^{\dagger}\y$\vspace{0.1cm} \\% , ${\bf r}=\y - \A_{{\cal S}_{i}}\hat{\x}_{i}$ 
% \ENDFOR \vspace{0.1cm}
% \end{algorithmic}
% \caption{Orthogonal Least-Squares (OLS)}
% \label{algo:1}
% \end{algorithm}
% \label{sec:algo}

\section{Generalized Orthogonal Least-Squares}

To formulate generalized OLS, we start by establishing a recursive relation between the components of the 
optimal solution to the 
$l_0$-constrained least-squares problem. Let ${\bf B}_i$ denote the sub-matrix of $\bf A$ constructed 
by selecting $i$ of its columns and let ${\P}_{i}={\bf B}_{i}{\bf B}^{\dagger}_{i}$ denote the projection 
matrix onto the span of the columns of $\B_i$, where 
${\bf B}^{\dagger}_{i} = \left({\bf B}^{T}_{i}{\bf B}_{i}\right)^{-1}{\bf B}^{T}_{i}$  is the Moore-Penrose 
pseudo-inverse of $\B_i$. Then, after appending $\B_i$ with another column vector $\a$ to form
$\B_{i+1} = [\B_i \;\;\; \a]$, we can write 
\begin{equation} \label{eq:8}
\begin{aligned}
{\P}_{i+1} & = \B_{i+1} \left({\bf B}^{T}_{i+1}{\bf B}_{i+1}\right)^{-1}{\bf B}^{T}_{i+1} \\
&=
\begin{bmatrix}
    {\B}_i & \a 
\end{bmatrix}
\begin{bmatrix}
    {\B}^{T}_i{\B}_i   &  {\bf B}^{T}_i\a  \\
    \a^{T}{\B}_i       & \a^{T}\a
\end{bmatrix}^{-1}
\begin{bmatrix}
    {\B}^{T}_i        \\
    \a^{T}     
\end{bmatrix}\\
&\stackrel{(a)}{=}
\begin{bmatrix}
    {\B}_i & \P^\bot_i\a
\end{bmatrix}
\begin{bmatrix}
    \left({\bf B}^{T}_i{\B}_i\right)^{-1}      &  0  \\
    0      & \left(\a^{T}\P^\bot_i\a\right)^{-1}
\end{bmatrix}
\begin{bmatrix}
    {\B}^{T}_i        \\
    \a^{T}\P^\bot_i
\end{bmatrix}\\
&\stackrel{(b)}{=} \P_i+\frac{\P^\bot_i\a\a^T\P^\bot_i}{\norm{\P^\bot_i\a}^2_2}
\end{aligned}
\end{equation}
where ($a$) follows after the LDU decomposition of the intermediate matrix inverse, 
i.e., we use the identity \cite{kailath2000linear}
\begin{equation}
\begin{aligned}
\hspace{-0.4cm}\begin{bmatrix}
    {\bf A}       &  {\bf  E}  \\
    {\bf C}         & {\bf D}   
\end{bmatrix}^{-1} = 
\begin{bmatrix}
    {\bf I}       &  {\bf A}^{-1} {\bf E}   \\
    {\bf 0}      & {\bf I}   
\end{bmatrix} 
 \begin{bmatrix}
    {\bf A}^{-1}       &  {\bf 0}   \\
    {\bf 0}        & {\bf \Delta}^{-1} 
\end{bmatrix}
\begin{bmatrix}
    {\bf I}       &  {\bf 0}  \\
    {\bf C}{\bf A}^{-1}        & {\bf I}   
\end{bmatrix}
\end{aligned}
\nonumber
\end{equation}
where we identify ${\bf A}={\B}^{T}_i{\B}_i$, ${\bf E}={\bf B}^{T}_i\a$, ${\bf C}=\a^{T}{\B}_i $, 
${\bf D}=\a^{T}\a$, and ${\bf \Delta}={\bf D}-{\bf C}{\bf A}^{-1}{\bf E}$, and introduce 
${\P}_i^{\bot}=\I-{\P}_i$. The identity ($b$) follows from the idempotent property of the projection matrix.
%, i.e.,
% \[
% {\P}_i^{\bot}={{\P}_i^{\bot}}^T={\P_i^{\bot}}^2. 
% \]
Alternatively, we write (\ref{eq:8}) as
\begin{equation}
{\P}_{i+1}^{\bot}={\P}_i^{\bot}-\frac{{\P}_i^{\bot}\a\a^T {\P}_i^{\bot}}{\norm{{\P}_i^{\bot}\a}_2^2}.
\label{eq:perp}
\end{equation}
Note that, (\ref{eq:perp}) is related to order-recursive least-squares \cite{kay2013fundamentals}. However, this specific derivation makes it suitable for iterative sparse reconstruction applications.

Now, the OLS algorithm in each step selects a column with index $j_s$ from the set ${\cal I}$ of the 
previously non-selected columns according to
\begin{equation}
\begin{aligned}
{j}_{s}&=\argmin_{j \in {\cal I}}{\norm{\y - \A_{{\cal S}_{i-1}\cup\{j\}}\A_{{\cal S}_{i-1}\cup\{j\}}^{\dagger}\y}_2}\\
&\stackrel{(a)}{=}\argmin_{j \in {\cal I}}{\norm{\left(\I-\P_i\right)\y}_2^2}\\
&=\argmin_{j \in {\cal I}}{\y^T\y-\y^T\P_i\y-\y^T\P_i^T\y+\y^T\P_i^T\P_i\y}\\
&\stackrel{(b)}{=}\argmin_{j \in {\cal I}}{\y^T\y-\y^T\P_i\y}\\
&\stackrel{(c)}{=}\argmax_{j \in {\cal I}}{\y^T\P_{i-1}\y+\y^T\frac{\P^\bot_{i-1}\a_j\a^T_j\P^\bot_{i-1}}{\norm{\P^\bot_{i-1}\a_j}^2_2}\y}\\
&\stackrel{(d)}{=} \argmax_{j \in {\cal I}}{\frac{\norm{\y^T\P^\bot_{i-1}\a_j}^2_2}{\norm{\P^\bot_{i-1}\a_j}^2_2}}
=\argmax_{j \in {\cal I}}{\left|\y^T\frac{{\bf P}_{i-1}^{\bot}{\bf a}_{j}}{\norm{{\bf P}_{i-1}^{\bot}{\bf a}_{j}}_2}\right|}
\end{aligned}
\end{equation}
where ($a$) follows from the definition of $\P_i$, ($b$) is due to $P_i^\perp$ being an idempotent projection 
matrix, ($c$) follows from eq. (\ref{eq:8}), and ($d$) is due to the fact that $\y^T\P_{i-1}\y$ is not a function 
of the optimization variable. 

%The recursion at the core of the procedure is given by the next proposition.
%\begin{proposition}\label{prop1}
%If $(i+1)^{th}$ iteration of OLS algorithm identifies column $\a$, then
%\begin{equation}
%{\P}_{i+1}^{\bot}={\P}_i^{\bot}-\frac{{\P}_i^{\bot}\a\a^T {\P}_i^{\bot}}{\norm{{\P}_i^{\bot}\a}_2^2}
%\end{equation}
%where ${\P}_i^{\bot}$ is the projection matrix onto the subspace orthogonal to span of previously selected columns.
%\end{proposition}
%The proof is given in the Appendix. Proposition 1 provides a computationally efficient way of computing OLS update with no requirement of finding projections in the first step of OLS algorithm. Next proposition summarizes this result.
%\begin{proposition}\label{prop2}
%For OLS algorithm, selection criterion can be equivalently written as
%\begin{equation}
%j_{max}=\argmax_{j \in {\cal I}}{\left|\y^T\frac{{\bf P}_{i-1}^{\bot}{\bf a}_{j}}{\norm{{\bf P}_{i-1}^{\bot}{\bf a}_{j}}_2}\right|}
%\end{equation}
%\end{proposition}

We propose a straightforward extension of OLS which selects multiple (say, $L$) columns of $\A$ 
in each step rather than choosing a single column, ultimately replacing the underdetermined 
$n \times m$ system of equations by an overdetermined $Lk \times k$ one. This strategy is 
motivated by the observation that the candidate columns whose projection onto the space orthogonal
to that spanned by the previously selected columns is strongly correlated with the observation
vector but not chosen in the current step of OLS will likely be selected in subsequent steps of the 
algorithm; therefore, selecting several ``good" candidates in each step accelerates the selection 
procedure and enables sparse reconstruction with fewer steps (and, therefore, fewer calculations
of the mutual correlations needed to perform the selection). More specifically, the proposed 
generalized OLS algorithm performs the following: in each step, the algorithm select $L$ columns 
of matrix $\A$ such that their normalized projection onto the orthogonal complement of the 
subspace spanned by previously chosen columns have the highest correlation with the observation 
vector among the non-selected columns. After such columns are identified, we update the orthogonal 
projection matrix by repeatedly applying (\ref{eq:perp}) $L$ times. We continue until a stopping 
criterion is met. Generalized orthogonal least-squares algorithm is formalized as Algorithm 1.

\subsection{Computational complexity} 

To analyze the computational complexity of GOLS, note that its first step involves a matrix-vector 
multiplication and a vector inner product; the computational cost of these two is dominated by the 
former and thus requires ${\cal O}\left(mn^2\right)$ operations. The other operations are those in 
Step 3 where a matrix-vector multiplication and a matrix addition, needing ${\cal O}\left(n^2\right)$ 
calculation, need to be repeated $\min\{k,\lfloor \frac{n}{L}\rfloor\}$ times. Therefore, the aggregate 
cost of operations in this step is ${\cal O}\left(mkn^2\right)$. Finally, finding the estimate $\hat{\x}$ 
entails solving a least-squares problem which can be implemented with a small cost 
${\cal O}\left(kn\right)$ by relying on a QR factorization of $\A_{{\cal S}_{k}}$. 
% Here we apply the modified Gram-Schmidt (MGS) algorithm because the coefficient matrix is assumed unstructured and dense. 
Therefore, the total complexity of the algorithm is ${\cal O}\left(mkn^2+kn\right)$.
%and is lower than the computational cost of OLS.

%\begin{figure}[h]
%\begin{minipage}[b]{.5\textwidth}
%\begin{subfigure}[b]{0.30\textwidth}
%  \centering
%    \includegraphics[width=0.99\textwidth]{ex_gaus}\quad\caption{ERR}
%    \end{subfigure}
%    \begin{subfigure}[b]{.30\textwidth}
%  \centering
%    \includegraphics[width=0.99\textwidth]{mse_gaus_log}\quad\caption{MSE}
%    \end{subfigure}
%    \begin{subfigure}[b]{.30\textwidth}
%  \centering
%    \includegraphics[width=0.99\textwidth]{run_gaus1}\quad\caption{Running time}
%\end{subfigure}
%\caption{Performance comparison of GOLS, OLS, OMP, $l_1$-norm minimization and LASSO
%for $n=64$, $m=128$, $\A$ having Gaussian ${\cal N}(0,1/n)$ entries, and the $k$ non-zero
%components of $\x$ drawn from ${\cal N}(0,1)$ distribution.}
%\end{minipage}
%\end{figure}
%\begin{figure}[h]
%\begin{minipage}[b]{.5\textwidth}
%\begin{subfigure}{.32\textwidth}
%  \centering
%    \includegraphics[width=0.99\textwidth]{ex_ber}\quad\caption{ERR}
%    \end{subfigure}
%    \begin{subfigure}{.32\textwidth}
%  \centering
%    \includegraphics[width=0.99\textwidth]{mse_ber_log}\quad\caption{MSE}
%        \end{subfigure}
%        \begin{subfigure}{.32\textwidth}
%  \centering
%    \includegraphics[width=0.99\textwidth]{run_ber1}\quad\caption{Running time}
%\end{subfigure}
%\caption{Performance comparison of GOLS, OLS, OMP, $l_1$-norm minimization and LASSO
%for $n=64$, $m=128$, $\A$ having Gaussian ${\cal N}(0,1/n)$ entries, and the $k$ non-zero
%components of $\x$ randomly and equally likely set to $1$ or $-1$.}
%\end{minipage}
%\end{figure}
%\begin{figure}[h]
%\begin{minipage}[b]{.5\textwidth}
%\begin{subfigure}{0.32\textwidth}
%  \centering
%    \includegraphics[width=0.99\textwidth]{ex_Aber}\quad\caption{ERR}
%    \end{subfigure}
%    \begin{subfigure}{.32\textwidth}
%  \centering
%    \includegraphics[width=0.99\textwidth]{mse_Aber_log}\quad\caption{MSE}
%        \end{subfigure}
%        \begin{subfigure}{.32\textwidth}
%  \centering
%    \includegraphics[width=0.99\textwidth]{run_Aber1}\quad\caption{Running time}
%\end{subfigure}
%\caption{Performance comparison of GOLS, OLS, OMP, $l_1$-norm minimization and LASSO
%for $n=64$, $m=128$, $\A$ having uniformly i.i.d entries from $\{+1\big\slash\sqrt{n},-1\big\slash\sqrt{n}\}$, and the $k$ non-zero
%components of $\x$ drawn from ${\cal N}(0,1)$ distribution.}
%\end{minipage}
%\end{figure}

\begin{figure*}[t]
%\begin{minipage}[b]{.5\textwidth}
\begin{subfigure}[b]{0.32\textwidth}
  \centering
    \includegraphics[width=1\textwidth]{ex_gaus}\quad\caption{ERR}
    \end{subfigure}
    \begin{subfigure}[b]{.32\textwidth}
  \centering
    \includegraphics[width=1\textwidth]{mse_gaus_log}\quad\caption{MSE}
    \end{subfigure}
    \begin{subfigure}[b]{.32\textwidth}
  \centering
    \includegraphics[width=1\textwidth]{run_gaus1}\quad\caption{Running time}
\end{subfigure}
\caption{Performance comparison of GOLS, OLS, OMP, $l_1$-norm minimization and LASSO
for $n=64$, $m=128$, $\A$ having Gaussian ${\cal N}(0,1/n)$ entries, and the $k$ non-zero
components of $\x$ drawn from ${\cal N}(0,1)$ distribution.}
%\end{minipage}
\end{figure*}
\begin{figure*}[t]
%\begin{minipage}[b]{.5\textwidth}
\begin{subfigure}{.32\textwidth}
  \centering
    \includegraphics[width=1\textwidth]{ex_ber}\quad\caption{ERR}
    \end{subfigure}
    \begin{subfigure}{.32\textwidth}
  \centering
    \includegraphics[width=1\textwidth]{mse_ber_log}\quad\caption{MSE}
        \end{subfigure}
        \begin{subfigure}{.32\textwidth}
  \centering
    \includegraphics[width=1\textwidth]{run_ber1}\quad\caption{Running time}
\end{subfigure}
\caption{Performance comparison of GOLS, OLS, OMP, $l_1$-norm minimization and LASSO
for $n=64$, $m=128$, $\A$ having Gaussian ${\cal N}(0,1/n)$ entries, and the $k$ non-zero
components of $\x$ randomly and equally likely set to $1$ or $-1$.}
%\end{minipage}
\end{figure*}
\begin{figure*}[t]
%\begin{minipage}[b]{.5\textwidth}
\begin{subfigure}{0.32\textwidth}
  \centering
    \includegraphics[width=1\textwidth]{ex_Aber}\quad\caption{ERR}
    \end{subfigure}
    \begin{subfigure}{.32\textwidth}
  \centering
    \includegraphics[width=1\textwidth]{mse_Aber_log}\quad\caption{MSE}
        \end{subfigure}
        \begin{subfigure}{.32\textwidth}
  \centering
    \includegraphics[width=1\textwidth]{run_Aber1}\quad\caption{Running time}
\end{subfigure}
\caption{Performance comparison of GOLS, OLS, OMP, $l_1$-norm minimization and LASSO
for $n=64$, $m=128$, $\A$ having uniformly i.i.d entries from $\{+1\big\slash\sqrt{n},-1\big\slash\sqrt{n}\}$, and the $k$ non-zero
components of $\x$ drawn from ${\cal N}(0,1)$ distribution.}
%\end{minipage}
\end{figure*}


\begin{algorithm}
\vspace{0.1cm}
\begin{tabularx}{\textwidth}{l>{$}c<{$}X}
\textbf{Input:}  \hspace{0.85cm}observation $\y$, coefficient matrix ${\bf A}$, sparsity $k$ 
\vspace{0.1cm}\\
\textbf{Output:} \hspace{0.5cm} recovered support ${\cal S}_k$, estimated signal $\hat{\x}_{k}$\vspace{0.1cm}\\
\textbf{Initialize:} \hspace{0.3cm}  ${\cal S}_0=\oldemptyset$, ${\bf P}_0^{\bot}={\bf I}$, ${\cal I} = \{1,2,\dots,m\}$
\end{tabularx}
\begin{algorithmic}
\FOR  { $i=1$ to $\min\{k,\lfloor \frac{n}{L}\rfloor\}$  }\vspace{0.1cm}
\STATE 1. Select $\{i_{s_1},\dots,i_{s_L}\}$ corresponding to $L$ largest terms:
$\left|\y^T\frac{{\bf P}_{i-1}^{\bot}{\bf a}_{j}}{\norm{{\bf P}_{i-1}^{\bot}{\bf a}_{j}}_2}\right|$ for $j \in {\cal I}$\vspace{0.1cm}\\
2. ${\cal S}_{i}={\cal S}_{i-1}\cup\{i_{s_1},\dots,i_{s_L}\}$ , ${\cal I} = {\cal I} \backslash {\cal S}_{i}$\vspace{0.05cm}\vspace{0.1cm}\\
3. ${\bf D}={\bf P^{\bot}}_{i-1}$ \vspace{0.1cm}\\
\FOR { $l=1$ to $L$ }\vspace{0.1cm}
\STATE  ${\bf d}=\frac{{\bf D}{\bf a}_{j_{s_l}}}{\norm{{\bf D}{\bf a}_{j_{s_l}}}_2}, {\bf D}={\bf D}-{\bf d}{\bf d}^T$\vspace{0.1cm}
\ENDFOR\vspace{0.1cm}\\
${\bf P}_{i}^{\bot}={\bf D}$ \vspace{0.1cm}\\
\ENDFOR \vspace{0.1cm}
\STATE $\hat{\x}_k=\A_{{\cal S}_{k}}^{\dagger}\y$
\end{algorithmic}
\caption{Generalized Orthogonal Least-Squares}
\label{algo:2}
\end{algorithm}
\section{SIMULATION RESULTS} \label{sec:sim}
To evaluate the algorithm, we compared its performance with four other sparse recovery algorithms as 
a function of the sparsity level $k$. In particular, we considered OMP, OLS, $l_1$-norm minimization 
\cite{candes2006robust}, and Least Absolute Shrinkage and Selection Operator (LASSO) 
\cite{tibshirani1996regression}. As typically done in benchmarking tests \cite{dai2009subspace}, we used 
CVX \cite{grant2008cvx} to implement the $l_1$ minimization and LASSO. The tuning parameter in LASSO 
is found by means of 10-fold cross validation. We draw entries of the coefficient matrix $\A$ from two
distributions. First, we generate entries of $\A$ by drawing independently from a Gaussian distribution with 
zero-mean and variance $n^{-1}$. We then consider two different scenarios for this choice of coefficient 
matrix: (1) the non-zero elements of $\x$ are independent and identically distributed normal random variables, 
and (2) the non-zero components of $\x$ are drawn uniformly from alphabet $\{+1,-1\}$. Second, the entries 
of $\A$ are drawn independently and uniformly from $\{+1\big\slash\sqrt{n},-1\big\slash\sqrt{n}\}$. In all 
settings, the locations of non-zero entries of $\x$ are drawn uniformly at random. The number of equations 
is $n=64$, the dimension of $\x$ is $m=128$; the experiment is repeated 1000 times. Performance of each 
algorithm is characterized by three metrics: (i) exact recovery rate (ERR), defined as the fraction of the 
correctly recovered signal components, (ii) mean-square error (MSE), measuring the distance between the 
unknown signal and its estimate, and (iii) the running time of the algorithm. Results for the Gaussian 
coefficient matrices are illustrated in Fig. 1 and Fig. 2. Fig. 1 shows the performance of the algorithms for 
non-zero values of $\x$ being normally distributed while Fig. 2 corresponds to the second scenario. Fig. 3 
shows the performance of the methods for $\A$ being constructed according to the second option while 
the non-zero values of $\x$ are normally distributed. As can be seen from Fig. 1 and Fig. 3, the generalized 
OLS (GOLS) outperforms all the competing methods in terms of the exact recovery rate, and is better than 
OLS and OMP in terms of the MSE. Moreover, the runtimes of GOLS is 2nd only to OMP but the accuracy 
of the latter is significantly worse than that of GOLS. In the case of $\{+1,-1\}$ non-zero entries of $\x$ 
studied in Fig. 2, $l_1$-norm/LASSO methods perform the best (and are the slowest) while the GOLS 
offers reasonably accurate performance at relatively high speed.


% Below is an example of how to insert images. Delete the ``\vspace'' line,
% uncomment the preceding line ``\centerline...'' and replace ``imageX.ps''
% with a suitable PostScript file name.
% -------------------------------------------------------------------------
% \begin{figure}[htb]

% \begin{minipage}[b]{1.0\linewidth}
%   \centering
%   \centerline{\includegraphics[width=8.5cm]{image1}}
% %  \vspace{2.0cm}
%   \centerline{(a) Result 1}\medskip
% \end{minipage}
% %
% \begin{minipage}[b]{.48\linewidth}
%   \centering
%   \centerline{\includegraphics[width=4.0cm]{image3}}
% %  \vspace{1.5cm}
%   \centerline{(b) Results 3}\medskip
% \end{minipage}
% \hfill
% \begin{minipage}[b]{0.48\linewidth}
%   \centering
%   \centerline{\includegraphics[width=4.0cm]{image4}}
% %  \vspace{1.5cm}
%   \centerline{(c) Result 4}\medskip
% \end{minipage}
% %
% \caption{Example of placing a figure with experimental results.}
% \label{fig:res}
% %
% \end{figure}
% To start a new column (but not a new page) and help balance the last-page
% column length use \vfill\pagebreak.

\section{CONCLUSION} \label{sec:concl}
We show that for Gaussian and Bernoulli coefficient matrices, Orthogonal Least-Squares (OLS) is with 
high probability guaranteed to recover any sparse signal from a low number of random linear measurements. 
Moreover, we introduced a greedy algorithm for sparse linear regression that generalizes OLS and forms 
the subset of features (i.e., columns of a coefficient matrix in an underdetermined system of equations) by 
sequentially selecting multiple candidate columns. Since multiple indices are selected without additional cost, 
the running time of the algorithm is reduced compared to OLS. Thus, generalized OLS is more favorable than 
convex optimization based methods whose complexity grows faster with the dimension of the problem, i.e., 
$n$ and $m$. Simulation studies demonstrate that the generalized orthogonal least-squares algorithm 
outperforms competing greedy methods, OLS and OMP, while being computationally more efficient than 
$l_1$-norm minimization and LASSO. 
%Future work will include search for more efficient techniques for column selection.
\bibliographystyle{IEEEbib}\small
\bibliography{IEEEabrv,refs}

\end{document}

main.bbl
\begin{thebibliography}{10}

\bibitem{donoho2006compressed}
David~L Donoho,
\newblock ``Compressed sensing,''
\newblock {\em Information Theory, IEEE Transactions on}, vol. 52, no. 4, pp.
  1289--1306, 2006.

\bibitem{barik2014sparsity}
Somsubhra Barik and Haris Vikalo,
\newblock ``Sparsity-aware sphere decoding: algorithms and complexity
  analysis,''
\newblock {\em Signal Processing, IEEE Transactions on}, vol. 62, no. 9, pp.
  2212--2225, 2014.

\bibitem{parvaresh2008recovering}
Farzad Parvaresh, Haris Vikalo, Sidhant Misra, and Babak Hassibi,
\newblock ``Recovering sparse signals using sparse measurement matrices in
  compressed dna microarrays,''
\newblock {\em Selected Topics in Signal Processing, IEEE Journal of}, vol. 2,
  no. 3, pp. 275--285, 2008.

\bibitem{lustig2007sparse}
Michael Lustig, David Donoho, and John~M Pauly,
\newblock ``Sparse mri: The application of compressed sensing for rapid mr
  imaging,''
\newblock {\em Magnetic resonance in medicine}, vol. 58, no. 6, pp. 1182--1195,
  2007.

\bibitem{elad2010role}
Michael Elad, Mario~AT Figueiredo, and Yi~Ma,
\newblock ``On the role of sparse and redundant representations in image
  processing,''
\newblock {\em Proceedings of the IEEE}, vol. 98, no. 6, pp. 972--982, 2010.

\bibitem{mishali2010theory}
Moshe Mishali and Yonina~C Eldar,
\newblock ``From theory to practice: Sub-nyquist sampling of sparse wideband
  analog signals,''
\newblock {\em Selected Topics in Signal Processing, IEEE Journal of}, vol. 4,
  no. 2, pp. 375--391, 2010.

\bibitem{elhamifar2009sparse}
Ehsan Elhamifar and Ren{\'e} Vidal,
\newblock ``Sparse subspace clustering,''
\newblock in {\em Computer Vision and Pattern Recognition, 2009. CVPR 2009.
  IEEE Conference on}. IEEE, 2009, pp. 2790--2797.

\bibitem{candes2006robust}
Emmanuel~J Cand{\`e}s, Justin Romberg, and Terence Tao,
\newblock ``Robust uncertainty principles: Exact signal reconstruction from
  highly incomplete frequency information,''
\newblock {\em Information Theory, IEEE Transactions on}, vol. 52, no. 2, pp.
  489--509, 2006.

\bibitem{beck2009fast}
Amir Beck and Marc Teboulle,
\newblock ``A fast iterative shrinkage-thresholding algorithm for linear
  inverse problems,''
\newblock {\em SIAM journal on imaging sciences}, vol. 2, no. 1, pp. 183--202,
  2009.

\bibitem{boyd2011distributed}
Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein,
\newblock ``Distributed optimization and statistical learning via the
  alternating direction method of multipliers,''
\newblock {\em Foundations and Trends{\textregistered} in Machine Learning},
  vol. 3, no. 1, pp. 1--122, 2011.

\bibitem{tropp2004greed}
Joel~A Tropp,
\newblock ``Greed is good: Algorithmic results for sparse approximation,''
\newblock {\em Information Theory, IEEE Transactions on}, vol. 50, no. 10, pp.
  2231--2242, 2004.

\bibitem{needell2010signal}
Deanna Needell and Roman Vershynin,
\newblock ``Signal recovery from incomplete and inaccurate measurements via
  regularized orthogonal matching pursuit,''
\newblock {\em Selected Topics in Signal Processing, IEEE Journal of}, vol. 4,
  no. 2, pp. 310--316, 2010.

\bibitem{tropp2007signal}
Joel~A Tropp and Anna~C Gilbert,
\newblock ``Signal recovery from random measurements via orthogonal matching
  pursuit,''
\newblock {\em Information Theory, IEEE Transactions on}, vol. 53, no. 12, pp.
  4655--4666, 2007.

\bibitem{pati1993orthogonal}
Yagyensh~Chandra Pati, Ramin Rezaiifar, and PS~Krishnaprasad,
\newblock ``Orthogonal matching pursuit: Recursive function approximation with
  applications to wavelet decomposition,''
\newblock in {\em Signals, Systems and Computers, 1993. 1993 Conference Record
  of The Twenty-Seventh Asilomar Conference on}. IEEE, 1993, pp. 40--44.

\bibitem{donoho2012sparse}
David~L Donoho, Yaakov Tsaig, Iddo Drori, and Jean-Luc Starck,
\newblock ``Sparse solution of underdetermined systems of linear equations by
  stagewise orthogonal matching pursuit,''
\newblock {\em Information Theory, IEEE Transactions on}, vol. 58, no. 2, pp.
  1094--1121, 2012.

\bibitem{needell2009cosamp}
Deanna Needell and Joel~A Tropp,
\newblock ``Cosamp: Iterative signal recovery from incomplete and inaccurate
  samples,''
\newblock {\em Applied and Computational Harmonic Analysis}, vol. 26, no. 3,
  pp. 301--321, 2009.

\bibitem{chen1989orthogonal}
Sheng Chen, Stephen~A Billings, and Wan Luo,
\newblock ``Orthogonal least squares methods and their application to
  non-linear system identification,''
\newblock {\em International Journal of control}, vol. 50, no. 5, pp.
  1873--1896, 1989.

\bibitem{soussen2013joint}
Charles Soussen, R{\'e}mi Gribonval, J{\'e}r{\^o}me Idier, and C{\'e}dric
  Herzet,
\newblock ``Joint k-step analysis of orthogonal matching pursuit and orthogonal
  least squares,''
\newblock {\em Information Theory, IEEE Transactions on}, vol. 59, no. 5, pp.
  3158--3174, 2013.

\bibitem{herzet2013exact}
C{\'e}dric Herzet, Charles Soussen, J{\'e}r{\^o}me Idier, and R{\'e}mi
  Gribonval,
\newblock ``Exact recovery conditions for sparse representations with partial
  support information,''
\newblock {\em Information Theory, IEEE Transactions on}, vol. 59, no. 11, pp.
  7509--7524, 2013.

\bibitem{herzet2016relaxed}
C{\'e}dric Herzet, Ang{\'e}lique Dr{\'e}meau, and Charles Soussen,
\newblock ``Relaxed recovery conditions for omp/ols by exploiting both
  coherence and decay,''
\newblock {\em Information Theory, IEEE Transactions on}, vol. 62, no. 1, pp.
  459--470, 2016.

\bibitem{rangan2009orthogonal}
Sundeep Rangan and Alyson~K Fletcher,
\newblock ``Orthogonal matching pursuit from noisy random measurements: A new
  analysis,''
\newblock in {\em Advances in Neural Information Processing Systems}, 2009, pp.
  540--548.

\bibitem{kailath2000linear}
Thomas Kailath, Ali~H Sayed, and Babak Hassibi,
\newblock {\em Linear estimation}, vol.~1,
\newblock Prentice Hall Upper Saddle River, NJ, 2000.

\bibitem{kay2013fundamentals}
Steven~M Kay,
\newblock {\em Fundamentals of Statistical Signal Processing: Practical
  Algorithm Development}, vol.~3,
\newblock Pearson Education, 2013.

\bibitem{tibshirani1996regression}
Robert Tibshirani,
\newblock ``Regression shrinkage and selection via the lasso,''
\newblock {\em Journal of the Royal Statistical Society. Series B
  (Methodological)}, pp. 267--288, 1996.

\bibitem{dai2009subspace}
Wei Dai and Olgica Milenkovic,
\newblock ``Subspace pursuit for compressive sensing signal reconstruction,''
\newblock {\em Information Theory, IEEE Transactions on}, vol. 55, no. 5, pp.
  2230--2249, 2009.

\bibitem{grant2008cvx}
Michael Grant and Stephen Boyd,
\newblock ``Cvx: Matlab software for disciplined convex programming,'' .

\end{thebibliography}
