# The Lawson-Hanson Algorithm with Deviation Maximization:
Finite Convergence and Sparse Recovery

Monica Dessole  Department of Mathematics “Tullio Levi Civita”, University of
Padova, Via Trieste 63, 35121 Padova, Italy  Marco Dell’Orto  Department of
Mathematics “Tullio Levi Civita”, University of Padova, Via Trieste 63, 35121
Padova, Italy  Fabio Marcuzzi  Department of Mathematics “Tullio Levi Civita”,
University of Padova, Via Trieste 63, 35121 Padova, Italy

(August 2021)

###### Abstract

In this work we devise a new version of the Lawson-Hanson with Deviation
Maximization (LHDM) algorithm, a block method for the solution of NonNegative
Least Squares (NNLS). It is shown that the new algorithm here presented
convergences in a finite number of steps under suitable conditions. An
extensive campaign of experiments is performed in order to evaluate the
performance gain with respect to the standard Lawson-Hanson algorithm. We also
explore the sparse recovery ability of LHDM, comparing it against several
$\ell_{1}$-minimization solvers in terms of solution quality and time-to-
solution on a large set of dense instances.

Keywords: nonnegative least squares, sparse recovery, block pivoting

AMS: 65F25, 65K99

##  1 Introduction

NonNegative Least Squares (NNLS) problems arise in many applications where
data points can be represented as nonnegative linear combinations of some
meaningful components. Such problems are frequently encountered in signal and
image processing and they are core problems in more complex computations, such
as nonnegative matrix and tensor decompositions, see e.g. [20, 18]. Moreover,
when dealing with underdetermined systems of equations, the nonnegativity
constraint is known to naturally enhance sparsity of the solution, that is the
solution attained has few nonzeros, see e.g. [6, 17, 31, 32]. An important
outcome of this body of work is that nonnegativity alone may attain a
satisfactory sparse recovery. Moreover, NNLS solvers can be adopted to solve
unconstrained least squares problems with minor adjustments. Over the last two
decades, sparsity has become one of the most relevant topics in signal
processing. In general, sparsity in signals describes the phenomenon where
high dimensional data can be expressed with few measurements and it results in
finding a sparse solution to underdetermined systems of equations. This
problem is highly nonconvex and its solution is NP-hard in general, although
it is well known that $\ell_{1}$-minimization leads to the sparsest solution
for a restricted class of matrices. This fact is known as $\ell_{0}-\ell_{1}$
equivalence and it has been found empirically [9] and theoretically [13, 16].

The first algorithm devised for NNLS is the Lawson-Hanson algorithm [19].
Since this seminal work, many modifications have been proposed in order to
improve the standard Lawson-Hanson algorithm: Bro and Jong [5] have proposed a
variation specifically designed for use in nonnegative tensor decompositions;
their algorithm, called “fast NNLS” (FNNLS), reduces the execution time by
avoiding redundant computations in nonnegative matrix factorization problems
arising in tensor decompositions and performs well with multiple right-hand-
sides problems, which is not the case here discussed, thus we omit a
comparison. Van Benthem and Keenan [30] presented a different NNLS solution
algorithm, namely “fast combinatorial NNLS” (FCNNLS), also designed for the
specific case of a large number of right-hand sides. Principal block pivoting
method introduced by Portugal et al. [22] is an active set type algorithm
which differs from the standard Lawson-Hanson algorithm, since the sequence of
iterates produced does not necessarily fall into the feasible region. The
convergence is ensured provided the objective function is strictly convex,
while when we deal with underdetermined matrices it is simply convex,
therefore this algorithm may fail in sparse recovering. Surprisingly, the
Lawson-Hanson algorithm [19] does not suffer from this drawback.

The Lawson-Hanson algorithm is an active set method which iteratively builds
an optimal solution by solving at each iteration an unconstrained least
squares subproblem, e.g. by computing a QR decomposition, on a subset of
columns of the initial matrix, in such a way that the objective value strictly
decreases while the iterates are kept within the feasible region. The
algorithm performs a column pivoting by selecting one column at a time and
therefore it relies only on BLAS–2 operations, resulting in poor performance.
In an actual implementation, the columns are not permuted, rather a vector of
indices is used to keep trace of columns ordering. In general, column pivoting
makes it more difficult to achieve high performances in QR computation, see
[1, 2, 3, 23, 4].

In this work, we address the NNLS problem in which the objective function
matrix is dense and we present a modified Lawson-Hanson algorithm in which
multiple column selection is performed in order to exploit BLAS–3 operations
for efficiency. The column selection is based on the “deviation maximization”
pivoting presented in [10] to obtain rank-revealing QR factorizations. The
main difference is that the column selection here does not aim at revealing
the rank of the matrix, but is driven from the active-set method to solve the
non-negativity constrained problem. The proposed algorithm is named Lawson-
Hanson with Deviation Maximization (LHDM). A first try of combining the two
techniques was presented in [11], where we observed good sparsity recovery
capacity and a significant performance gain with respect to the standard
Lawson-Hanson algorithm when dealing with underdetermined systems. However, in
this preliminary version the algorithm may fall in an infinite loop without
converging to a solution: indeed, it is known that manipulating multiple
columns at once may allow cycling. In the present paper, we introduce a new
mechanism which allows a similar performance gain, but the termination in a
finite number of steps is ensured. Moreover, we explore its sparsity recovery
ability against several $\ell_{1}$-minimization solvers, extending a
comparison that has been done in the recent literature [21].

The rest of the paper is organized as follows. In Section 2, we formally
introduce the NNLS problem and the Lawson-Hanson algorithm. We present the
deviation maximization algorithm for column selection and the new version of
the Lawson-Hanson with Deviation Maximization algorithm. We provide a
theoretical result about finite convergence of LHDM by generalizing the
analogous result for the standard Lawson-Hanson. In Section 3, we introduce
the $\ell_{0}$-minimization problem as sparse recovery and the convex
$\ell_{1}$-minimization problem. We recall results that ensure the so-called
$\ell_{0}-\ell_{1}$ equivalence and conditions under which the underdetermined
NNLS problem has a unique solution. Last, we recall how arbitrary signed
sparse recovery can be attained by NNLS solvers. In Section 5, we present a
comparison of LHDM against several $\ell_{1}$-minimization solvers in terms of
performance and solution found. The results are shown with an extensive
campaign of experiments. Finally, Section 6 concludes the paper.

##  2 Active set algorithms for nonnegative least squares

Consider the nonnegative least squares problem

$\min_{\mathbf{x}}\|A\mathbf{x}-\mathbf{b}\|_{2},\qquad\text{s.t.
}{\mathbf{x}\geq 0},$ (1)

where $A=(\mathbf{a}_{1}\ \dots\ \mathbf{a}_{n})$ is a matrix of size $m\times
n$, $\mathbf{b}$ is a vector of length $m$. In this work, we address the case
in which $A$ is a dense matrix, i.e. it is stored as a full array with $mn$
entries. The case in which $A$ is sparse would require the design of a
specific implementation in order to exploit the presence of many null entries.
Trivially, a solution ${\mathbf{x}}^{\star}$ of (1) also solves

$\min_{\mathbf{x}}\frac{1}{2}\|A\mathbf{x}-\mathbf{b}\|_{2}^{2},\qquad\text{s.t.
}{\mathbf{x}\geq 0}.$

Define $\phi(\mathbf{x}):=\|A\mathbf{x}-\mathbf{b}\|_{2}^{2}$. Note that the
gradient of $\phi$ at a point $\mathbf{x}$, i.e. steepest descend direction
with opposite sign, is

$\nabla\phi(\mathbf{x})=A^{T}A\mathbf{x}-A^{T}\mathbf{b}=-A^{T}(\mathbf{b}-A\mathbf{x})=-A^{T}\mathbf{r}(\mathbf{x}),$

where $\mathbf{r}$ is the residual function
$\mathbf{r}(\mathbf{x})=\mathbf{b}-A\mathbf{x}$. The set of points satisfying
the constraints is called feasible region, here
$\Omega=\left\\{\mathbf{x}\in\mathbb{R}^{n}:\ \mathbf{x}\geq 0\right\\}$, and
a point $\mathbf{x}$ is said to be feasible if it belongs to the feasible
region. Let us recall the following basic definition.

######  Definition 1.

Let $\mathbf{x}\in\mathbb{R}^{n}$ be a feasible point, i.e. $\mathbf{x}\geq
0$. A vector $\mathbf{s}\in\mathbb{R}^{n}$ is a feasible direction at
$\mathbf{x}$ if there exists $\varepsilon>0$ such that

$\mathbf{x}+\varepsilon\mathbf{s}\geq 0.$

Notice that if $\mathbf{s}$ is a feasible direction at $\mathbf{x}$ and
$s_{i}<0$, then we must have $x_{i}>0$. The following theorem characterizes
the solutions of the convex problem (1). For a proof see e.g. [19].

######  Theorem 1.

(Karush-Kuhn-Tucker conditions for NNLS) A point
$\mathbf{x}^{\star}\in\mathbb{R}^{n}$ is a solution of problem (1) if and only
if there exists $\mathbf{w}^{\star}\in\mathbb{R}^{n}$ and a partition
$Z^{\star}\cup P^{\star}=\left\\{1,\dots,n\right\\}$ such that

$\displaystyle\mathbf{w}^{\star}=A^{T}(\mathbf{b}-A\mathbf{x}^{\star}),$ (2)
$\displaystyle x^{\star}_{i}=0,\ i\in Z^{\star},\quad x^{\star}_{i}>0,\ i\in
P^{\star},$ (3) $\displaystyle w^{\star}_{i}\leq 0,\ i\in Z^{\star},\quad
w^{\star}_{i}=0,\ i\in P^{\star}.$ (4)

Notice that equations (3-4) imply

$x^{\star}_{i}w^{\star}_{i}=0,\ i=1,\dots,n,$ (5)

which are the well known complementarity conditions. The points $\mathbf{x}$
and $\mathbf{w}$ are referred as primal and dual variables, respectively.
Notice that the primal and dual variables are related as follows
$\mathbf{w}=A^{T}\mathbf{r}=-\nabla\phi(\mathbf{x})$, i.e. the dual vector
$\mathbf{w}$ is the opposite of the gradient direction, i.e. the steepest
descend direction, at $\mathbf{x}$.

Active set methods are in the family of descend algorithms, that is they look
for the solution by decreasing the objective function value, in particular at
each iteration the new solution is a feasible point in a feasible descend
direction with respect to the current solution. Recall that the $i$-th
constraint is said to be active if it holds with equality, here $x_{i}=0$, it
is said to be passive when it holds with strict inequality, here $x_{i}>0$,
otherwise it is violated, namely $x_{i}<0$.

###  2.1 The Lawson-Hanson algorithm

The first algorithm to solve (1) is due to Lawson and Hanson [19], and it is
presented in Algorithm 1, in an efficient formulation that uses Householder
and Givens orthogonal transformations, which are then detailed here below and
in sec. 4. It is an active set method and a particular case of the algorithm
introduced in [26] for the least squares problem with linear inequality
constraints.

Algorithm 1 LH$(A,\mathbf{b})$

1:${P}_{0}=\emptyset$, ${Z}_{0}=\left\\{1,\dots,M\right\\}$,
$\mathbf{x}^{(0)}=0$, $\mathbf{w}^{(0)}=A^{T}\mathbf{b}$, $s=0$

2:while ${Z}_{s}\neq\emptyset$ and $\max(\mathbf{w}^{(s)}_{Z})>0$ do
$\triangleright$ outer loop

3: $j_{s}\in\operatorname*{argmax}_{i}{w}^{(s)}_{i}$

4: set ${P}_{s+1}={P}_{s}\cup\left\\{j_{s}\right\\}$,
${Z}_{s+1}={Z}_{s}\setminus\left\\{j_{s}\right\\}$

5: if $s=0$ then

6: set $R^{(1)}$ as the Householder reflector related to $\mathbf{a}_{j_{s}}$
and compute $\mathbf{b}^{(1)}$

7: else

8: use rank-1 modification to update $R^{(s+1)}$ and update
$\mathbf{b}^{(s+1)}$

9: end if

10: compute $\mathbf{y}^{(s+1)}_{P}$ as the solution of
$R^{(s+1)}\mathbf{y}=\mathbf{b}^{(s+1)}$, set $\mathbf{y}^{(s+1)}_{Z}=0$

11: while $\min\left(\mathbf{y}^{(s+1)}_{P}\right)\leq 0$ do $\triangleright$
inner loop

12: $s=s+1$

13:
$\alpha=\min\left\\{{x^{(s-1)}_{i}}\left({x^{(s-1)}_{i}-y^{(s)}_{i}}\right)^{-1}:\
i\in P_{s},\ y^{(s)}_{i}\leq 0\right\\}$

14:
$\mathbf{x}^{(s)}=\mathbf{x}^{(s-1)}+\alpha\left(\mathbf{y}^{(s)}-\mathbf{x}^{(s-1)}\right)$

15: ${Q}=\left\\{i\ :\ x^{(s)}_{i}\leq 0\right\\}$

16: ${P}_{s+1}={P}_{s}\setminus Q$

17: ${Z}_{s+1}={Z}_{s}\cup Q$

18: $R^{(s+1)}=R^{(s)}$

19: for $j\in Q$ do

20: delete column corresponding to index $j$ from $R^{(s+1)}$

21: reduce $R^{(s+1)}$ to triangular form by Givens transformation and update
$\mathbf{b}^{(s+1)}$

22: end for

23: compute $\mathbf{y}^{(s+1)}_{P}$ as the solution of
$R^{(s+1)}\mathbf{y}=\mathbf{b}^{(s+1)}$, set $\mathbf{y}^{(s+1)}_{Z}=0$

24: end while

25: $\mathbf{x}^{(s+1)}=\mathbf{y}^{(s+1)}$

26: $\mathbf{w}^{(s+1)}_{Z}=\left(A^{(s+1)}_{Z}\right)^{T}\mathbf{b}$,
$\mathbf{w}^{(s+1)}_{P}=0$

27: $s=s+1$

28:end while

29:return ${P}^{\star}={P}_{s}$, ${Z}^{\star}={Z}_{s}$,
$\mathbf{x}^{\star}=\mathbf{x}^{(s)}$

Let us now describe the main features of the Lawson-Hanson algorithm. Let
${P}_{s}$ denote the passive set and ${Z}_{s}$ denote the active set at $k$-th
algorithmic step, with cardinality respectively $n_{s}$ and $n-n_{s}$. Let us
define the following submatrices

$\displaystyle A^{(s)}_{P}$
$\displaystyle=\left(\mathbf{a}_{i_{1}}\dots\mathbf{a}_{i_{n_{s}}}\right),$
$\displaystyle i_{1},\dots i_{n_{s}}\in{P}_{s},$ $\displaystyle{A}^{(s)}_{Z}$
$\displaystyle=\left(\mathbf{a}_{j_{1}}\dots\mathbf{a}_{j_{n-n_{s}}}\right),$
$\displaystyle j_{1},\dots j_{n-n_{s}}\in{Z}_{s}.$

We define the $k$-th iterate as
$\mathbf{x}^{(s)}=\left(x^{(s)}_{1},\dots,x^{(s)}_{n}\right)^{T}\in\mathbb{R}^{n}$.
With an analogous notation to the one introduced above for matrices, we have

$\displaystyle\mathbf{x}^{(s)}_{P}$
$\displaystyle=\left(x^{(s)}_{i_{1}},\dots,x^{(s)}_{i_{n_{s}}}\right)^{T}\in\mathbb{R}^{n_{s}},$
$\displaystyle i_{1},\dots,i_{n_{s}}\in{P}_{s},$
$\displaystyle{\mathbf{x}}^{(s)}_{Z}$
$\displaystyle=\left(x^{(s)}_{j_{1}},\dots,x^{(s)}_{j_{n-n_{s}}}\right)^{T}\in\mathbb{R}^{n-n_{s}},$
$\displaystyle j_{1},\dots,j_{n-n_{s}}\in{Z}_{s}.$

Since the iterates produced by Lawson-Hanson algorithm do not leave the
feasible region, we always have ${\mathbf{x}}^{(s)}_{Z}=0$. The vector
$\mathbf{y}^{(s)}_{P}$ is the least squares solution of the following
unconstrained subproblem

$\min_{\mathbf{y}}\left\|A^{(s)}_{P}\mathbf{y}-\mathbf{b}\right\|,$ (6)

and $\mathbf{y}^{(s)}_{Z}=0$. Notice that these vectors change size at each
iteration.

With the notation introduced, we have

$\displaystyle\mathbf{r}^{(s)}$
$\displaystyle=\mathbf{b}-A\mathbf{x}^{(s)}=\mathbf{b}-A^{(s)}_{P}\mathbf{x}^{(s)}_{P},$
$\displaystyle\mathbf{w}^{(s)}$
$\displaystyle=A^{T}\left(\mathbf{b}-A\mathbf{x}^{(s)}\right)=A^{T}\left(\mathbf{b}-A^{(s)}_{P}\mathbf{x}^{(s)}_{P}\right)=A^{T}\mathbf{r}^{(s)},$
$\displaystyle\mathbf{w}^{(s)}_{Z}$
$\displaystyle=\left(A^{(s)}_{Z}\right)^{T}\left(\mathbf{b}-A^{(s)}_{P}\mathbf{x}^{(s)}_{P}\right)=\left(A^{(s)}_{Z}\right)^{T}\mathbf{r}^{(s)},$
$\displaystyle\mathbf{w}^{(s)}_{P}$
$\displaystyle=\left(A^{(s)}_{P}\right)^{T}\left(\mathbf{b}-A^{(s)}_{P}\mathbf{x}^{(s)}_{P}\right)=0.$

where the last identity is a consequence of normal equations for (6), since at
the end of the outer loop we have $\mathbf{x}^{(s)}_{P}=\mathbf{y}^{(s)}_{P}$.
Clearly we have $\mathbf{w}^{(s)}\cdot\mathbf{x}^{(s)}=0$ for every $k$. As
shown in [19], Algorithm 1 terminates in a finite number of steps and
$\mathbf{x}^{(s)},\mathbf{w}^{(s)},P_{s},Z_{s}$ satisfy KKT conditions (3-4)
on termination.

Since the columns of $A^{(s)}_{P}$ are linearly independent, the vector
$\mathbf{y}^{(s)}_{P}$ is the unique solution of the linear system
$R^{(s)}\mathbf{y}^{(s)}_{P}=\mathbf{b}^{(s)}$, where $R^{(s)}$ is the square
upper triangular matrix of order $n_{s}=|P_{s}|$ in the QR decomposition of
$A^{(s)}_{P}$, that is

$A^{(s)}_{P}=Q^{(s)}\left(\begin{array}[]{c}R^{(s)}\\\ 0\end{array}\right).$

Here the passive set $P_{s}$ is an ordered set (the order in which the column
indices have been inserted) and therefore $P_{s}$ also imposes the column
ordering in QR computation. Similarly, the updated right-hand side is given by

$\mathbf{b}^{(s)}=\left(Q^{(s)}\right)^{T}\mathbf{b}.$

With respect to the original matrix $A$, we obtain the following partial QR
decomposition with column pivoting

$\left(Q^{(s)}\right)^{T}A\
\Pi^{(s)}=\left(\begin{array}[]{cc}R^{(s)}&B^{(s)}\\\
0&C^{(s)}\end{array}\right),$

where $\Pi^{(s)}$ is the permutation matrix that permuted the elements in
$P_{s}$ in the leftmost positions. In an actual implementation [19], the
columns are not physically permuted, rather a permutation vector is used to
keep track of the column ordering. Here, $B^{(s)}$ has size
$n_{s}\times(n-n_{s})$ and $C^{(s)}$ has size $(m-n_{s})\times(n-n_{s})$.

Let us spend the rest of this section by giving some details about the
implementation of Algorithm 1. Notice that the matrix $R^{(s)}$ is never
computed from scratch, instead it is modified by means of rank-$1$ updates or
downdated by means of Givens rotations. Let us briefly detail these
strategies. At the beginning of the outer loop, the new passive set is
incremented of exactly one index, therefore $n_{s+1}=|P_{s+1}|=n_{s}+1$. Then,
the corresponding Householder vector is computed and applied to the right-hand
side and to the trailing matrix, i.e.

$\left(I-\mathbf{v}_{s}^{T}\mathbf{v}_{s}\right)C^{(s)}=\left(\begin{array}[]{cc}r_{n_{s}+1}&b_{n_{s}+1,2}\
\dots\ b_{n_{s}+1,n-n_{s}}\\\ 0&C^{(s+1)}\end{array}\right),$ (7)

where we supposed that the selected column was in the leading position. The
new block $C^{(s+1)}$ has size $(m-n_{s}-1)\times(n-n_{s}-1)$. The partial
factorization is updated as follows

$R^{(s+1)}=\left(\begin{array}[]{cc}R^{(s)}&\mathbf{b}_{1}^{(s)}\\\
&r_{n_{s}+1}\end{array}\right),\quad
B^{(s+1)}=\left(\begin{array}[]{ccc}\mathbf{b}_{2}^{(s)}&\dots&\mathbf{b}_{n-n_{s}}^{(s)}\\\
b_{n_{s}+1,2}&\dots&b_{n_{s}+1,n-n_{s}}\end{array}\right).$

On the other hand, when a column is removed from the passive set, the R factor
is almost an Hessenberg matrix, i.e. it has some nonzero entries in its first
lower diagonal, and it has to be reduced to triangular form. This can be done
by computing a sequence of Givens rotations, which have to be applied to the
whole matrix $A$ and to the right-hand side. In this way, the orthogonal
factor $Q^{(s)}$ is never computed nor stored. We detail the modified QR
decomposition procedure in Section 4.

###  2.2 The Lawson-Hanson algorithm with Deviation Maximization

In this section we propose a new active set method in which multiple columns
are manipulated at a time and for which finite termination can be established
in exact arithmetic. Before coming to the description of the main algorithm,
let us discuss a strategy, for columns selection, based upon the deviation
maximization [10].

Recall that for an $m\times k$ matrix $C=(\mathbf{c}_{1}\ \dots\
\mathbf{c}_{k})$ whose columns $\mathbf{c}_{j}$ are non-null, the correlation
matrix $\Theta$ has entries

$\theta_{ij}=\frac{\mathbf{c}_{i}^{T}\mathbf{c}_{j}}{\|\mathbf{c}_{i}\|\|\mathbf{c}_{j}\|},\quad
1\leq i,j\leq k.$ (8)

In particular, we have
$\Theta=\left(CD^{-1}\right)^{T}CD^{-1}=D^{-1}C^{T}CD^{-1}$, where $D$ is the
diagonal matrix with entries $d_{i}=\|\mathbf{c}_{i}\|\ ,\ 1\leq i\leq k$. It
is immediate to see that $\Theta$ is symmetric positive semidefinite, it has
only ones on the diagonal, and its entries range from $-1$ to $1$. Notice that
$\theta_{ij}$ is the cosine of
$\alpha_{ij}=\alpha(\mathbf{c}_{i},\mathbf{c}_{j})\in[0,\pi)$, the angle
(modulo $\pi$) between $\mathbf{c}_{i}$ and $\mathbf{c}_{j}$. In order to
emphasize this geometric interpretation, from now on we refer to $\Theta$ as
the cosine matrix.

Algorithm 2 presents the deviation maximization algorithm for column
selection. Here, the procedure for computing the candidate set $I$ has been
adapted to fit the Lawson-Hanson algorithm and it is slightly different from
the one presented in [10].

Algorithm 2 Deviation Maximization:
$[J]$=DM$(A,\mathbf{u}_{1},\mathbf{u}_{2},\tau_{1},\tau_{2},\delta,k_{\max})$

1:${J}=\left\\{j\ :\ j\in\operatorname*{argmax}\mathbf{u}_{1}\right\\}$

2:${I}=\left\\{i\ :\ \mathbf{u}_{l}(i)\geq\tau_{l}\max\mathbf{u}_{l},\ i\neq
j,l=1,2\right\\}$

3:if $|I|>k_{\max}$ then

4: $I=\left\\{i_{s}\in I\ :\ k=1,\dots,k_{\max}\right\\}$

5:end if

6:compute the cosine matrix $\Theta$ associated to $[A]_{:,{I}}$

7:for $i\in{I}$ do

8: if $|\theta_{i,j}|<\delta,\forall j\in{J}$ then

9: ${J}={J}\cup\left\\{i\right\\}$

10: end if

11:end for

12:return $J$

The resulting procedure, namely the Lawson-Hanson algorithm with Deviation
Maximization (LHDM), is presented in Algorithm 3, and it terminates in a
finite number of steps as shown in Theorem 3 below.

######  Theorem 2.

Let $A^{(s)}_{P}$ be a full column rank matrix, and let $R^{(s+1)}$ be the R
factor of the QR decomposition of $A^{(s+1)}_{P}=\left(A^{(s)}_{P}\
A^{(s)}_{J}\right)$, where $A^{(s)}_{J}$ is the submatrix of $A$ that is
obtained by taking columns with indices in the set ${J_{s}}$, which is chosen
by the deviation maximization with $0<\tau_{1},\tau_{2}<1$, $k_{\max}>1$, and

$\delta\leq\frac{{\tau}_{2}}{k_{\max}-1}.$ (9)

Then

$\bar{\sigma}^{(s+1)}\geq\sigma_{n_{s+1}}(A)\frac{\bar{\sigma}^{(s)}}{\sigma_{1}(A)}\frac{1}{\sqrt{2(n-n_{s+1})n_{s+1}}}\frac{\tau_{2}\sqrt{1-\tau_{2}}}{k_{s}^{2}n_{s}},$
(10)

where $k_{s}=|J_{s}|$ and $\bar{\sigma}^{(s+1)}$ is the smallest singular
value of $R^{(s+1)}$, which is in particular nonsingular.

###### Proof.

Apply Theorem 2 of [10] to conclude. ∎

It is known, it has been showed in [10], that even if the matrix $A^{(s)}_{P}$
is well conditioned, there could be a potentially dramatic increase in the
condition number of $A^{(s+1)}_{P}$ when using the deviation maximization
strategy, as well as for the other column selection strategy like the standard
column pivoting. A famous example is the Kahan matrix, however it is a very
unlikely occurrence in practice.

The choice (9) is not the only one possible: in fact, the values of $\delta$
and $\tau_{2}$ can be set independently, as shown in [10].

Algorithm 3 LHDM$(A,\mathbf{b},\tau_{1},\tau_{2},\tau_{c},k_{\max})$

1:${P}_{0}=\emptyset$, ${Z}_{0}=\left\\{1,\dots,M\right\\}$,
$\mathbf{x}^{(0)}=0$, $\mathbf{w}^{(0)}=A^{T}\mathbf{b}$,
$\mathbf{y}^{(0)}=0$, $s=0$

2:$R^{(0)}=0$, $B^{(0)}=0$, $C^{(0)}=A$

3:initialize $\mathbf{u}^{(0)}$ as the vector of column norms of $A$

4:while ${Z}_{s}\neq\emptyset$ and $\max\left(\mathbf{w}^{(s)}\right)>0$ do
$\triangleright$ outer loop

5:
$[J_{s}]$=DM$(C^{(s)},\mathbf{w}^{(s)},\mathbf{u}^{(s)},\tau_{1},\tau_{2},\delta,k_{\max})$

6: ${P}_{s+1}={P}_{s}\cup J_{s}$

7: ${Z}_{s+1}={Z}_{s}\setminus J_{s}$

8: if $s=0$ then

9: set $R^{(1)}$ as the R factor of the QR decomposition of $A^{(1)}_{J}$ and
compute $\mathbf{b}^{(1)}$, $B^{(1)}$ and $C^{(1)}$

10: else

11: use rank-$k$ update to get $R^{(s+1)}$, $B^{(s+1)}$, $C^{(s+1)}$, and
update $\mathbf{b}^{(s+1)}$

12: end if

13: compute $\mathbf{y}^{(s+1)}_{P}$ as the solution of
$R^{(s+1)}\mathbf{y}=\mathbf{b}^{(s+1)}$, set $\mathbf{y}^{(s+1)}_{Z}=0$

14: while $\min\left(\mathbf{y}^{(s+1)}_{J}\right)\leq 0$ do $\triangleright$
ensure a feasible descend direction

15: delete the last element from $J_{s}$

16: ${P}_{s+1}={P}_{s}\cup J_{s}$

17: ${Z}_{s+1}={Z}_{s}\setminus J_{s}$

18: delete the last column of $R^{(s+1)}$ and update $\mathbf{b}^{(s+1)}$,
$B^{(s+1)}$, $C^{(s+1)}$,

19: compute $\mathbf{y}^{(s+1)}_{P}$ as the solution of
$R^{(s+1)}\mathbf{y}=\mathbf{b}^{(s+1)}$, set $\mathbf{y}^{(s+1)}_{Z}=0$

20: end while

21: while $\min\left(\mathbf{y}^{(s+1)}_{P}\right)\leq 0$ do $\triangleright$
inner loop

22: $s=s+1$

23:
$\alpha=\min\left\\{{x^{(s-1)}_{i}}\left({x^{(s-1)}_{i}-y^{(s)}_{i}}\right)^{-1}:\
i\in P_{s},\ y^{(s)}_{i}\leq 0\right\\}$

24:
$\mathbf{x}^{(s)}=\mathbf{x}^{(s-1)}+\alpha\left(\mathbf{y}^{(s)}-\mathbf{x}^{(s-1)}\right)$

25: ${Q}=\left\\{i\ :\ x^{(s)}_{i}\leq 0\right\\}$

26: ${P}_{s+1}={P}_{s}\setminus Q$

27: ${Z}_{s+1}={Z}_{s}\cup Q$

28: $R^{(s+1)}=R^{(s)}$

29: for $j\in Q$ do

30: delete column corresponding to index $j$ from $R^{(s+1)}$

31: reduce $R^{(s+1)}$ to triangular form by Givens transformation and update
$\mathbf{b}^{(s+1)}$

32: end for

33: compute $\mathbf{y}^{(s+1)}_{P}$ as the solution of
$R^{(s+1)}\mathbf{y}=\mathbf{b}^{(s+1)}$, set $\mathbf{y}^{(s+1)}_{Z}=0$

34: end while

35: $\mathbf{x}^{(s+1)}=\mathbf{y}^{(s+1)}$

36: $\mathbf{w}^{(s+1)}=A^{T}\left(\mathbf{b}-A\mathbf{x}^{(s+1)}\right)$

37: update $\mathbf{u}^{(s)}$ as the vector of column norms of $C^{(s+1)}$

38: $s=s+1$

39:end while

40:return ${P}^{\star}={P}_{s}$, ${Z}^{\star}={Z}_{s}$,
$\mathbf{x}^{\star}=\mathbf{x}^{(s)}$

Let us now prove the main result of this work.

######  Theorem 3.

Consider the NNLS problem

$\min\|A\mathbf{x}-\mathbf{b}\|^{2}_{2},\qquad\text{s.t. }{\mathbf{x}\geq 0}.$

$k_{\max}>1$, Algorithm 3 terminates in a finite number of steps. On
termination, $\mathbf{x}$ will be a solution vector and $\mathbf{w}$ will be a
dual solution vector.

###### Proof.

In order to ensure finite convergence of Algorithm 3, we have to ensure the
following

  1. 1.

termination of the outer loop, by showing that
$\mathbf{s}^{(s+1)}=\mathbf{y}^{(s+1)}-\mathbf{x}^{(s)}$ is a feasible descent
direction;

  2. 2.

termination of the inner loop with a nonempty passive set.

Let us address the first task. Theorem 2 ensures that
$A^{(s+1)}_{P}=\left(A^{(s)}_{P}\ A^{(s)}_{J}\right)$ has numerically linearly
independent columns. Let us claim the following statement: if
$\mathbf{y}^{(s+1)}_{J}>0$, then
$\mathbf{s}^{(s+1)}=\mathbf{y}^{(s+1)}-\mathbf{x}^{(s)}$ is a feasible descend
direction at $\mathbf{x}^{(s)}$. Then, if $\mathbf{y}^{(s+1)}_{j}\leq 0$ for
some $j\in J_{s}$, the last column of $A^{(s)}_{J}$ is dropped and the
corresponding index $j_{k_{s}}$ is moved back to the active set $Z_{s+1}$. We
then look for a new value of $\mathbf{y}^{(s+1)}$ such that

$R^{(s+1)}_{1}\mathbf{y}^{(s+1)}_{P}=\mathbf{b}^{(s+1)}_{1},$

where $R^{(s+1)}_{1}$ is equal to $R^{(s+1)}$ where the last column is
missing, and the right hand side $\mathbf{b}^{(s+1)}_{1}$ has been properly
updated, e.g. by applying the Givens rotation that zeros out the last diagonal
element of $R^{(s+1)}$, corresponding to the deleted column. We continue
dropping a column index at a time from ${J_{s}}$ until
$\mathbf{y}^{(s+1)}_{j}>0$ for all indices $j$ left in $J_{s}$. It should be
noticed that this procedure terminates at most when only the first column of
$A^{(s)}_{J}$ is left, that is the same column index that would be selected by
the standard Lawson-Hanson algorithm, i.e. it would guarantee the positivity
of the inserted component.

For what concerns the second task, the inner loop terminates each time in at
most $n_{s+1}-1$ steps. Set $s=s+1$ and notice that the first time the inner
loop is entered, we have $\mathbf{y}^{(s)}_{J}>0$, implying
$\mathbf{x}^{(s)}_{J}=\alpha\mathbf{y}^{(s)}_{J}>0$ and thus $J_{s}$ and $Q$
do not intersect. Therefore, indices in $J_{s}$ are not removed before
computing $\mathbf{y}^{(s+1)}$. Moreover, at least the first index $j_{1}$ of
$J_{s}$ is left, because of Lemma 23.17 in [19], and thus the inner loop must
terminate at most in $n_{s+1}-1$ step with a nonempty passive set
$P_{s+1}\supseteq\\{j_{1}\\}$.

Let us now prove the claim, i.e. if $\mathbf{y}^{(s+1)}_{J}>0$, then
$\mathbf{s}^{(s+1)}=\mathbf{y}^{(s+1)}-\mathbf{x}^{(s)}$ is a feasible descend
direction at $\mathbf{x}^{(s)}$. Consider the least squares solution
$\mathbf{y}$ computed at step 13 of Algorithm 1 with deviation maximization
pivoting, that is $\mathbf{y}$ solves

$\mathbf{y}_{P}=\operatorname*{argmin}_{\mathbf{y}}\left\|A^{(s)}_{P}\mathbf{y}-\mathbf{b}\right\|\quad\mathbf{y}_{Z}=0.$
(11)

where $Z=Z_{s}$ and $P=P_{s}$. The additional inner loop ensures that
$y^{(s)}_{j}>0$ for $j\in J_{s}$ on termination. This, together with
$\mathbf{x}^{(s-1)}\geq 0$, means that
$\mathbf{s}^{(s)}:=\mathbf{y}^{(s)}-\mathbf{x}^{(s-1)}$ is a feasible
direction. Let us now prove that $\mathbf{s}^{(s)}$ is also a descend
direction. Said $\mathbf{r}^{(s-1)}=\mathbf{b}-A\mathbf{x}^{(s-1)}$, let us
define $\hat{\mathbf{z}}$ as follows

$\hat{\mathbf{z}}=\operatorname*{argmin}_{\mathbf{z}_{Z}=0}\|A\mathbf{z}-\mathbf{r}^{(s-1)}\|=\operatorname*{argmin}_{\mathbf{z}_{Z}=0}\|A(\mathbf{x}^{(s-1)}+\mathbf{z})-\mathbf{b}\|,$
(12)

where $Z=Z_{s}$ and $P=P_{s}$. Notice that $\hat{\mathbf{z}}$ can be obtained
as

$\hat{\mathbf{z}}_{P}=\operatorname*{argmin}_{\mathbf{z}}\|A^{(s)}_{P}\mathbf{z}-\mathbf{r}^{(s-1)}\|,\quad\hat{\mathbf{z}}_{Z}=0.$
(13)

Then we have

$\displaystyle\|A\hat{\mathbf{z}}-\mathbf{r}^{(s-1)}\|$
$\displaystyle=\|A\hat{\mathbf{z}}-(\mathbf{b}-A\mathbf{x}^{(s-1)})\|=\|A(\hat{\mathbf{z}}+\mathbf{x}^{(s-1)})-\mathbf{b}\|$
(14) $\displaystyle\geq\|A\mathbf{y}^{(s)}-\mathbf{b}\|,$

since $\hat{\mathbf{z}}+\mathbf{x}^{(s-1)}$ vanishes on $Z_{s}$ and the
minimum is reached at $\mathbf{y}^{(s)}$. On the other hand

$\displaystyle\|A\mathbf{y}^{(s)}-\mathbf{b}\|$
$\displaystyle=\|A\mathbf{y}^{(s)}-A\mathbf{x}^{(s-1)}-\mathbf{b}+A\mathbf{x}^{(s-1)}\|=\|A(\mathbf{y}^{(s)}-\mathbf{x}^{(s-1)})-\mathbf{r}^{(s-1)}\|$
(15) $\displaystyle\geq\|A\hat{\mathbf{z}}-\mathbf{r}^{(s-1)}\|,$

since $\mathbf{y}^{(s)}-\mathbf{x}^{(s-1)}$ vanishes on $Z_{s}$ and the
minimum is reached at $\hat{\mathbf{z}}$. Therefore

$\|A\hat{\mathbf{z}}-\mathbf{r}^{(s-1)}\|=\|A\mathbf{y}^{(s)}-\mathbf{b}\|,$

and the inequality chain (15) collapses in an equality chain. As a consequence
of the uniqueness of minimizers of (11) and (12), we have

$\hat{\mathbf{z}}=\mathbf{s}^{(s)}=\mathbf{y}^{(s)}-\mathbf{x}^{(s-1)}.$ (16)

Consider the function

$\varphi(\varepsilon)=\frac{1}{2}\|A(\mathbf{x}^{(s-1)}+\varepsilon\mathbf{s}^{(s)})-\mathbf{b}\|^{2}.$

We have

$\displaystyle\varphi^{\prime}(0)$
$\displaystyle=(\mathbf{s}^{(s)})^{T}A^{T}(A\mathbf{x}^{(s-1)}-\mathbf{b})=-(\mathbf{s}^{(s)})^{T}\mathbf{w}^{(s-1)}=-(\mathbf{y}^{(s)})^{T}\mathbf{w}^{(s-1)}$
(17) $\displaystyle=\sum_{i\in J}y^{(s)}_{i}(-w^{(s-1)}_{i})<0,$

since $w^{(s-1)}_{i}>0$ for all ${i\in J}$ and $y^{(s)}_{i}>0$ for all ${i\in
J}$. Therefore there exists $\varepsilon>0$ such that

$\|A(\mathbf{x}^{(s-1)}+\varepsilon\mathbf{s}^{(s)})-\mathbf{b}\|^{2}<\|A\mathbf{x}^{(s-1)}-\mathbf{b}\|^{2}=\phi(\mathbf{x}^{(s-1)}).$

By definition (12) and relation (16), we have

$\|A(\mathbf{x}^{(s-1)}+\mathbf{s}^{(s)})-\mathbf{b}\|^{2}\leq\|A(\mathbf{x}^{(s-1)}+\varepsilon\mathbf{s}^{(s)})-\mathbf{b}\|^{2}<\|A\mathbf{x}^{(s-1)}-\mathbf{b}\|^{2}.$
(18)

Thus, we take $\varepsilon=1$ (i.e. $\mathbf{x}^{(s)}=\mathbf{y}^{(s)}$), if
$\mathbf{y}^{(s)}$ is feasible, otherwise the algorithm enters the inner loop
and we choose a smaller step length $\alpha$, more precisely the largest
possible to keep the new iterate within the feasible region. It’s easy to
check that

$\displaystyle\|A(\mathbf{x}^{(s-1)}+\alpha\mathbf{s}^{(s)})-\mathbf{b}\|^{2}$
$\displaystyle=\min_{\varepsilon}\left\\{\|A(\mathbf{x}^{(s-1)}+\varepsilon\mathbf{s}^{(s)})-\mathbf{b}\|^{2}:\
\mathbf{x}^{(s-1)}+\varepsilon\mathbf{s}^{(s)}\geq 0\right\\}$
$\displaystyle\leq\|A\mathbf{x}^{(s-1)}-\mathbf{b}\|^{2}.$

We deduce that the objective function is always non increasing. Equation (18)
also implies that every time step 35 is reached, we have a strictly smaller
value of the objective function. Hence at every iteration of the outer loop,
the current objective value is strictly smaller than the preceding one. Since
the solution solves a least squares subproblem like those of steps 13 and 33,
the corresponding passive set must be different from the preceding one. As the
number of possible combinations of active and passive sets is finite, namely
it is equal to $2^{n}$, the algorithm must terminate in at most $2^{n}$ outer
loop repetitions. ∎

The proposed strategy may look naive. In order to give an idea about the
reasons why this should work well, let us prove the following.

######  Theorem 4.

If $\tau_{1}=1$ and $\tau_{2}\leq 1/2$, then in the outer loop we have
$\mathbf{y}^{(s)}_{J}>0$ for all $s>0$.

###### Proof.

The proof is an extension of Lemma 23.17 in [19], in which we end up showing

$\mathbf{y}^{(s)}_{J}=(C^{T}C)^{-1}\mathbf{w}^{(s+1)}_{J},$

where $C=A[n_{s}+1:m,n_{s}+k_{s}]$ and
$\mathbf{w}^{(s+1)}_{J}=(\omega,\dots,\omega)$, where
$\omega=\max\mathbf{w}^{(s+1)}>0$, because $\tau_{1}=1$. The matrix $C^{T}C$
is $\tau_{2}$-scaled diagonally dominant because it has been selected by
deviation maximization, with $\tau_{2}\leq 1/2$, hence its inverse is also
diagonally dominant with a positive diagonal [24], yielding

$\mathbf{y}^{(s)}_{J}>0.$

∎

This result shows that it is possible to analytically prove that cycling is
avoided without the need of an additional loop to ensure a feasible descend
direction. However, $\tau_{1}=1$ may be a too demanding choice and the
algorithm may loose in performance. In Section 5, we give experimental
evidence of the significant performance gain that this strategy can lead.

##  3 Sparse recovery

In this section we present the problem of sparse recovery and its connections
with nonnegative least squares and with the deviation maximization technique.
Sparse recovery is a fundamental problem in compressed sensing, signal
denoising, statistical model selection and related fields, and it consists in
finding the sparsest solution to an underdetermined system of equations, see
e.g. [14, 12, 15]. There is a wide literature in the field of signal
processing, in which a matrix is usually referred as a dictionary and its
columns are called atoms. Here, we abandon this terminology in favour of the
classical linear algebra one.

Given a vector $\mathbf{x}\in\mathbb{R}^{n}$, we define its support
${S}={S}(\mathbf{x})=\left\\{i:\ x_{i}\neq 0\right\\}$. We aim at identifying
the solution $\mathbf{x}^{\star}$ of the linear system $Ax=b$ with the
sparsest support ${S}^{\star}$. Formally, we consider the optimization problem

$\min\|\mathbf{x}\|_{0},\qquad\text{s.t. }\ A\mathbf{x}=\mathbf{b}.$ (19)

Recall the “$\ell_{0}$-norm” of a vector $\mathbf{x}$ is defined as

$\|\mathbf{x}\|_{0}=\left|{S}(\mathbf{x})\right|=\left|\left\\{i:x_{i}\neq
0\right\\}\right|,$

and it is not an actual norm because it does not hold that for any scalar $s$
we have $\|s\mathbf{x}\|=|s|\ \|\mathbf{x}\|$. Actually, problem (19) of
combinatorial search is NP-hard, and in general we seek an approximate
solution. Two well known approaches are: “greedy algorithms”, that focus on
the determination of the support and then a simple LS solution on it, and
“relaxation methods”, that focus on smoothing the $l_{0}$ norm and then adopt
smooth optimization. A famous greedy algorithm is the Orthogonal-Matching-
Pursuit (OMP), while a popular convex relaxation technique is the so-called
Basis Pursuit (BP) [8] problem

$\min\|\mathbf{x}\|_{1},\qquad\text{s.t. }\ A\mathbf{x}=\mathbf{b},$ (20)

which is a convex optimization problem, since $\|\mathbf{x}\|_{1}$ is a convex
function of $x$ and the solutions of $A\mathbf{x}=\mathbf{b}$ form a convex
set. We will concentrate on this second approach. Indeed, the two problems
(19) and (20) yield the same solution provided it is sparse enough. This
result is known as $\ell_{0}-\ell_{1}$ equivalence, and it has been found
empirically [9] and theoretically [13, 16].

###  3.1 Exact recovery

There are different conditions under which we have $\ell_{0}-\ell_{1}$
equivalence, i.e. an exact recovery. In the literature, these conditions are
usually found by showing the uniqueness of solution to problems (19) and (20)
and their coincidence. As a consequence of these results, efficient algorithms
to solve (20) can be used to find the sparsest solution to an underdetermined
system of equations.

####  3.1.1 Uniqueness based on RIP

Let us recall the following definition.

######  Definition 2.

Consider a matrix $A\in\mathbb{R}^{m\times n}$, and suppose there exists
$\delta_{s}\in(0,1)$ such that for every $m\times s$ submatrix $A_{s}$ of $A$
and for every $y\in\mathbb{R}^{s}$ we have

$(1-\delta_{s})\|y\|_{2}^{2}\leq\|A_{s}y\|_{2}^{2}\leq(1-\delta_{s})\|y\|_{2}^{2}.$
(21)

We say that $A$ has the $s$-order Restricted Isometry Property (RIP) with
constant $\delta_{s}$.

For a given $s$, one is usually interested in the smallest constant
$\delta_{s}$ for which (21) holds, which is referred as the $s$-order
restricted isometry constant.

A popular result [7] states that if $\delta_{2s}<1$, then problem (19) has a
unique solution $\mathbf{x}$ with support size obeying $\|\mathbf{x}\|_{0}\leq
s$. Moreover, if $\delta_{2s}<\sqrt{2}-1$, then the solution to problem (20)
also solves problem (19). However, the problem of establishing whether a given
matrix $A$ fulfills the $s$-order RIP is NP-hard in general [28].

####  3.1.2 Uniqueness based on ERC

Let us state the following result from [29].

######  Theorem 5.

Consider a linear system $A\mathbf{x}=\mathbf{b}$. If there exists a solution
$\bar{\mathbf{x}}$ with support $\bar{{S}}={S}(\bar{\mathbf{x}})$ such that

$\max_{i\not\in\bar{{S}}}\left\|A^{\dagger}_{\bar{{S}}}\mathbf{a}_{i}\right\|_{1}<1,$
(22)

then $\bar{\mathbf{x}}$ is the unique solution to the minimum $\ell_{0}$
problem (19), which can be recovered by solving the minimum $\ell_{1}$ problem
(20).

We refer to (22) as Exact Recovery Condition (ERC). This provides an easy
sufficient optimality check for a given support ${S}$, but finding a support
${S}$ satisfying (22) is a combinatorial problem.

###  3.2 Sparse recovery by NNLS

The nonnegativity constraint is known to naturally produce sparse solutions,
see e.g. [6, 17, 31, 32]. An important outcome of this body of work is that
nonnegativity alone may attain a satisfactory sparse recovery. There are also
important similarities. For example, [17] shows that the before mentioned
Orthogonal-Matching-Pursuit (OMP) resembles the Lawson-Hanson algorithm except
for the internal loop, which has indeed an evident benefit on nonnegative
recovery. They note also that this internal loop acts as a shrinkage: it
brings to zero the component which is relatively more negative and brings to
positive with a reduced magnitude the other components that have become
negative. So, it modifies all the components and, more importantly, reduces
the $l_{1}$-norm of the solution. In synthesis, the Lawson-Hanson algorithm
finds monotonically the solution if the growing support is correct, otherwise
enters in a inner loop that reduces the support and diminishes the
$l_{1}$-norm of the solution.

Notice also that arbitrary signed sparse recovery is easily achievable. Given
$\bar{\mathbf{x}}\in\mathbb{R}^{n}$, decompose it as
$\bar{\mathbf{x}}=\bar{\mathbf{x}}^{+}-\bar{\mathbf{x}}^{-}$, where
$\bar{\mathbf{x}}^{+}\geq 0$ and $\bar{\mathbf{x}}^{-}\geq 0$. Then the
solutions of the linear system $A\bar{\mathbf{x}}=\mathbf{b}$ can be attained
as the solutions of the nonnegative least squares problem

$\min_{\mathbf{x}}\left\|(A\
-A){\mathbf{x}}-\mathbf{b}\right\|_{2}^{2},\quad{\mathbf{x}}\geq 0,$ (23)

where ${\mathbf{x}}=\left(\begin{array}[]{l}\bar{\mathbf{x}}^{+}\\\
\bar{\mathbf{x}}^{-}\end{array}\right)\in\mathbb{R}^{2n}$. This has been shown
e.g. in [17] and it is sometimes referred as “positivity trick”.

####  3.2.1 Uniqueness for NNLS problem

This topic has been treated in [6, 25, 17].

Let us define the following class of matrices

${\cal M}^{+}=\left\\{A\in\mathbb{R}^{m\times n}:\ \exists\
\mathbf{h}\in\mathbb{R}^{m}\text{ s.t.
}A^{T}\mathbf{h}=\mathbf{w}>0\right\\}.$ (24)

A necessary condition in order to have a unique solution to the NNLS problem
(1) is that $A$ belongs to ${\cal M}^{+}$ [31].

######  Definition 3.

The coherence or mutual coherence or two-sided coherence of a matrix $A$ is
defined as

$\mu(A)=\max_{i<j}|\theta_{ij}|\left(=\max_{i>j}|\theta_{ij}|\right).$ (25)

In [29] it is demonstrated the following uniqueness result, under the
necessary condition (24).

######  Theorem 6.

Consider a linear system $A\mathbf{x}=\mathbf{b}$ such that $A\in{\cal
M}^{+}$. If there exists a nonnegative solution $\bar{\mathbf{x}}$ with
$\|\bar{\mathbf{x}}\|_{0}\leq s$ and

$\mu(A)<\frac{1}{2s-1},$ (26)

then $\bar{\mathbf{x}}$ is the unique nonnegative solution and the minimum
$\ell_{1}$-norm solution to $A\mathbf{x}=\mathbf{b}$.

###  3.3 Sparsity enhancing methods and approximate measurements

Let us now consider the more general situation in which the linear system does
not have to be exactly solved, i.e. we rather look for an approximate solution
satisfying $A\mathbf{x}=\mathbf{b}+{\boldsymbol{\varepsilon}}$, where
${\boldsymbol{\varepsilon}}$ is the vector containing the error in the
measurements. A famous method addressing this problem is the so-called LASSO
[27]. Following [17], we address to the solution of the following closely
related problem

$\min\|\mathbf{x}\|_{1}^{2}+\lambda\|A\mathbf{x}-\mathbf{b}\|^{2}_{2}.$ (27)

As we already pointed out that no theoretical result states that NNLS solvers
ensure sparse recovery, problem (27) can be interpreted as a sparsity
enhancing or squared $\ell_{1}$-regularization method consisting in adding a
penalization term to the objective function. This problem can be thought as a
weighted-sum form of a multiobjective problem where the parameter $\lambda$
controls the tradeoff between the two objectives. This technique suffers from
an evident drawback that is the choice of $\lambda$, which is not explicitly
related norm to the sparsity, nor to the linear system.

In [17], problem (27) is recast as a nonnegative least squares problem. By
making use of the positivity trick, problem (27) becomes

$\min\left\|\left(\begin{array}[]{c}1\ \dots\ 1\\\ \lambda{A}\ -\lambda
A\end{array}\right)\left(\begin{array}[]{l}\mathbf{x}^{+}\\\
\mathbf{x}^{-}\end{array}\right)-\left(\begin{array}[]{c}0\\\
\lambda\mathbf{b}\end{array}\right)\right\|^{2}_{2},\qquad\text{s.t.
}\mathbf{x}^{+},\mathbf{x}^{-}\geq 0$ (28)

which will be referred as $\ell_{1}$-NNLS. A nice feature of this form is that
the matrix involved clearly belongs to ${\cal M}^{+}$ and uniqueness of the
solution is possible for any $\lambda$. The equivalence between problem (20)
and problem (28) above can be established for $\lambda\rightarrow\infty$,
meaning that the squared residual is much more important then the
$\ell_{1}$-norm, which becomes more and more vacuous as $\lambda$ gets larger.
This suggests it could be sufficient to solve a simple NNLS problem in order
to get a sparse solution.

##  4 Upgrading the QR decomposition

We have already mentioned that in the methods here presented it is not
necessary to compute the QR decomposition from scratch, but instead one can
upgrade the current QR decomposition. In this section we detail the “updating”
and “downdating” procedures to be adopted in the outer and inner loop,
respectively. In the next subsection we see the general case, while in the
following we see a simplified downdating applicable with the positivity trick.

###  4.1 The general case

Let us first consider the problem of updating the QR decomposition in the
outer loop. Let $A$ be an $m\times n$ matrix and $B$ an $m\times k$ matrix,
both full rank by columns, such that $k+n\leq m$, and suppose moreover that
the columns in $B$ are linearly independent from those in $A$. We want to
exploit the knowledge of the decomposition $A=QR$ in order to compute
$\widetilde{Q},\widetilde{R}$ such that

$\widetilde{A}\coloneqq(A\ B)=\widetilde{Q}\widetilde{R}.$

The following procedure is based on Householder transformations, and it is
simply the Block Recursive Householder decomposition (without computing the
initial QR decomposition $A=QR$, that has already been computed in this
context).

Algorithm 4 Updating

1:$A=QR$, $B$

2:$B=Q^{T}B$

3:$\widetilde{Q}=Q$, $\widetilde{R}=(R\ B)$

4:for $i=1,\dots,k$ do

5: $[\beta,v]=$house$(B(n+i:m,i))$ $\triangleright$ compute Householder
reflector

6: $\widetilde{R}(n+i:m,n+i:m)=(I_{m-(n+i)+1}-\beta v\
v^{T})\widetilde{R}(n+i:m,n+i:m)$

7: ${H}_{n+i}=I_{m}$

8: ${H}_{n+i}(n+i:m,n+i:m)=(I_{m-(n+i)+1}-\beta v\ v^{T})$

9: $\widetilde{Q}={H}_{n+i}\widetilde{Q}$

10:end for

Note that explicit computation of $\widetilde{Q}$ is not necessary if we
express it as a product of Householder matrices

$\widetilde{Q}=Q{H}_{n+1}\cdots{H}_{n+k}.$

Operation count for the above algorithm (without taking into account the
explicit computation of $\widetilde{Q}$ in steps 7–9) if that of an
Householder $QR$ decomposition of a matrix of size $m-(n+1)\times k$, i.e.
about $2k^{2}(m-(n+1)-k/3)$ flops, plus the operation count in step 2. If we
suppose that $Q$ itself is expressed as a product of Householder matrices,
i.e.

${Q}={H}_{n}\cdots{H}_{n},$

we get about $kn(2m-n)$ flops. The whole procedure involves about
$2k^{2}(m-(n+1)-k/3)+kn(2m-n)=\mathcal{O}(k^{2}m+knm)$ flops.

Let us now consider the downdating problem. Consider a set of $k$ column
indices $\left\\{j_{1},\dots,j_{k}\right\\}$, with $1\leq
j_{1}<\dots<j_{k}<n$, and let $\widetilde{A}$ be the $m\times(n-k)$ matrix
whose columns are the columns of $A$ with index not in
$\left\\{j_{1},\dots,j_{k}\right\\}$. Again, we want to exploit the knowledge
of the decomposition $A=QR$ in order to compute $\widetilde{Q},\widetilde{R}$
such that

$\widetilde{A}=\widetilde{Q}\widetilde{R}.$

The matrix $\widetilde{A}$ has a structure close to the Hessenberg form.
Suppose $m=8$, $n=7$, $k=2$, $j_{1}=2$ and $j_{2}=5$. Then we have

$\widetilde{A}=\begin{bmatrix}\begin{array}[]{ccccc}\times&\times&\times&\times&\times\\\
&\times&\times&\times&\times\\\ &\times&\times&\times&\times\\\
&&\times&\times&\times\\\ &&&\times&\times\\\ &&&\times&\times\\\
&&&&\times\\\ &&&&\\\ \end{array}\end{bmatrix}.$

The following procedure is based on Givens rotations. We do this in a column-
major fashion, using Givens rotation to annihilate the highlighted entries one
at a time:

$\displaystyle\begin{bmatrix}\begin{array}[]{ccccc}\times&\times&\times&\times&\times\\\
&\times&\times&\times&\times\\\ &\boldsymbol{\times}&\times&\times&\times\\\
&&\times&\times&\times\\\ &&&\times&\times\\\ &&&\times&\times\\\
&&&&\times\\\ &&&&\\\
\end{array}\end{bmatrix}\rightarrow\begin{bmatrix}\begin{array}[]{ccccc}\times&\times&\times&\times&\times\\\
&\times&\times&\times&\times\\\ &&\times&\times&\times\\\
&&\boldsymbol{\times}&\times&\times\\\ &&&\times&\times\\\ &&&\times&\times\\\
&&&&\times\\\ &&&&\\\
\end{array}\end{bmatrix}\rightarrow\begin{bmatrix}\begin{array}[]{ccccc}\times&\times&\times&\times&\times\\\
&\times&\times&\times&\times\\\ &&\times&\times&\times\\\ &&&\times&\times\\\
&&&\times&\times\\\ &&&\boldsymbol{\times}&\times\\\ &&&&\times\\\ &&&&\\\
\end{array}\end{bmatrix}\rightarrow$
$\displaystyle\begin{bmatrix}\begin{array}[]{ccccc}\times&\times&\times&\times&\times\\\
&\times&\times&\times&\times\\\ &&\times&\times&\times\\\ &&&\times&\times\\\
&&&\boldsymbol{\times}&\times\\\ &&&&\times\\\ &&&&\times\\\ &&&&\\\
\end{array}\end{bmatrix}\rightarrow\begin{bmatrix}\begin{array}[]{ccccc}\times&\times&\times&\times&\times\\\
&\times&\times&\times&\times\\\ &&\times&\times&\times\\\ &&&\times&\times\\\
&&&&\times\\\ &&&&\times\\\ &&&&\boldsymbol{\times}\\\ &&&&\\\
\end{array}\end{bmatrix}\rightarrow\begin{bmatrix}\begin{array}[]{ccccc}\times&\times&\times&\times&\times\\\
&\times&\times&\times&\times\\\ &&\times&\times&\times\\\ &&&\times&\times\\\
&&&&\times\\\ &&&&\boldsymbol{\times}\\\ &&&&\\\ &&&&\\\
\end{array}\end{bmatrix}.$

Thus, columns with indices $j_{l}<j<j_{l+1}$ have exactly $l$ subdiagonal
entries to annihilate. Note that one could also zero out one subdiagonal at a
time, with exactly the same number of floating point operations. The following
algorithm describes the resulting procedure (in column-major fashion).

Algorithm 5 Downdating

1:$A=QR$, $j_{1}<\dots<j_{k}$

2:$\widetilde{Q}=Q$, $\widetilde{R}=R$, $j_{k+1}=n+1$

3:for $l=1,\dots,k$ do

4: for $j=j_{l}+1,\dots,j_{l+1}-1$ do

5: for $p=1,\dots,l$ do

6: $i=j-p+1$

7: $\widetilde{G}=$givens$(\widetilde{A}(i-1:i,j))$ $\triangleright$ compute
Givens rotation

8: $\widetilde{R}(i-1:i,j+1:m)=\widetilde{G}\widetilde{R}(i-1:i,j+1:m)$

9: $\widetilde{Q}(:,i-1:i)=\widetilde{Q}(:,i-1:i)\widetilde{G}^{T}$

10: end for

11: end for

12:end for

13:$\widetilde{R}(:,j_{1},\dots,j_{k})=[\ ]$ $\triangleright$ delete columns
$j_{1},\dots,j_{k}$

Again, the explicit computation of $\widetilde{Q}$ is not necessary if we
express it as a product of Givens rotation matrices

$\widetilde{Q}=Q{G}_{1}^{T}\cdots{G}_{t}^{T},$

where $t$ is the total number of computed rotations and $G_{l}$ is the
$m\times m$ matrix obtained enlarging the corresponding $\widetilde{G}_{l}$.

The operation count can be computed by observing that this procedure boils
down to the annihilation of $k$ subdiagonals of different length, as in the
second exemplification. After removing column $j_{l}$, we have to zero out the
last element of each column of the remaining $m\times(n-j_{l}-(k-l))$ matrix:
in fact the columns from $j_{l}$-th (not included) until the last one are
$n-j_{l}$, but there are also $k-l$ columns we don’t have to downdate, those
with indices $j_{l+1},\dots,j_{k}$, that have already been removed. Define
$n_{l}=(n-j_{l}-(k-l))$. Then, this operation takes
$\sum_{i=0}^{n_{l}}6(n_{l}-i)=\mathcal{O}(3n_{l}^{2})=\mathcal{O}(3(n-j_{l}-(k-l))^{2})$.
We compute the worst case of the total operation count, in which the first $k$
columns are removed: we have $j_{l}=l$ with $1\leq l\leq k$, thus we have to
zero out the last $k$ subdiagonals of the remaining $m\times(n-k)$ matrix.
This involves
$\mathcal{O}(\sum_{l=1}^{k}3(n-k)^{2})=\mathcal{O}(3n^{2}k+3k^{3}-6nk^{2})$
flops.

The updating and downdating procedures allow to modify the current QR
decomposition with a remarkable computational saving. In a practical
implementation, one can apply each time the orthogonal transformations to the
whole matrix, or only to the columns whose indices are in the current passive
set. In the first case, it is not necessary to keep trace of the
transformations themselves. In the second case, we need to keep in memory the
sequence of orthogonal transformations computed so far; moreover, an
additional vector is needed to keep trace of the number of orthogonal
transformations that have been applied to each column.

###  4.2 The positivity trick case

Let us now consider the particular case in which we deal with the NNLS problem
(23). In such a case, the objective function’s matrix is $[A\ ;\ -A]$. The
updating procedure by Householder triangularization is unchanged, while the
downdating procedure by Givens rotations in the inner loop can be avoided with
a simple “sign flip”, as we show here.

Suppose that $\mathbf{y}^{(s)}$ is not a feasible point, thus
$\mathbf{y}^{(s)}_{k}<0$, for some $k\in P_{s}$. Since
$\mathbf{s}^{(s)}:=\mathbf{y}^{(s+1)}-\mathbf{x}^{(s-1)}$ is a feasible
descend direction, we have $k\not\in J_{s}$. If $k\leq n$, the index $k$
corresponds to a column of $A$, otherwise if $n<k\leq 2n$, $k$ corresponds to
a column of $-A$. Let us denote by $\bar{k}$ the twin index of $k$, i.e.
$\bar{k}=k+n$ if $k\leq n$, and $\bar{k}=k-n$ if $n<k\leq 2n$. Let us show
that if we let
$P^{(s+1)}=(P^{(s)}\cup\left\\{\bar{k}\right\\})\setminus\left\\{{k}\right\\}$,
then $\mathbf{y}^{(s+1)}$ is a feasible point and
$\mathbf{s}^{(s+1)}:=\mathbf{y}^{(s+1)}-\mathbf{x}^{(s-1)}$ is a feasible
descend direction.

We first show that $\mathbf{y}^{(s+1)}$ is a feasible point. With a little
abuse of notation, we will denote by $k$ both the index itself and its
position within the passive set $P_{s}$. If we write by columns
$A^{(s)}_{P}=(\mathbf{a}_{1}\ \dots\mathbf{a}_{k-1}\ \mathbf{a}_{k}\
\mathbf{a}_{k+1}\dots\mathbf{a}_{n_{s}})$, then we have
$A^{(s+1)}_{P}=(\mathbf{a}_{1}\ \dots\mathbf{a}_{k-1}\ -\mathbf{a}_{k}\
\mathbf{a}_{k+1}\dots\mathbf{a}_{n_{s}})$. Consider the QR decomposition

$A^{(s)}_{P}=Q^{(s)}\ R^{(s)}=Q^{(s)}\ (\mathbf{r}_{1}\ \dots\mathbf{r}_{k-1}\
\mathbf{r}_{k}\ \mathbf{r}_{k+1}\dots\mathbf{r}_{n_{s}}),$

then we can write

$A^{(s+1)}_{P}=Q^{(s)}\ R^{(s+1)}=Q^{(s)}\ (\mathbf{r}_{1}\
\dots\mathbf{r}_{k-1}\ -\mathbf{r}_{k}\
\mathbf{r}_{k+1}\dots\mathbf{r}_{n_{s}}).$

Indeed, if $H$ is the Householder reflector of a vector $\mathbf{x}$, i.e.
$H\mathbf{x}=\|\mathbf{x}\|\ \mathbf{e}_{1}$, then
$H(-\mathbf{x})=-\|\mathbf{x}\|\ \mathbf{e}_{1}$. Now, $\mathbf{y}^{(s)}$ is
the solution of the least squares problem

$\mathbf{y}^{(s)}_{P}=\operatorname*{argmin}_{\mathbf{y}}\left\|A^{(s)}_{P}\mathbf{y}-\mathbf{b}\right\|,\quad\mathbf{y}^{(s)}_{Z}=0,$
(29)

which can be obtained as the solution of

$R^{(s)}\mathbf{y}^{(s)}_{P}=\mathbf{b}^{(s)}_{P}:=(Q^{(s)})^{T}\mathbf{b},\quad\mathbf{y}^{(s)}_{Z}=0.$
(30)

Let us split
$\mathbf{y}^{(s)}_{P}=(\mathbf{y}^{(s)}_{1},y^{(s)}_{k},\mathbf{y}^{(s)}_{2})^{T}$,
where $\mathbf{y}^{(s)}_{1}\in\mathbb{R}^{k-1}$,
$\mathbf{y}^{(s)}_{2}\in\mathbb{R}^{n_{s}-k}$, and $y^{(s)}_{k}<0$. Similarly,
we write
$\mathbf{b}^{(s)}_{P}=(\mathbf{b}^{(s)}_{1},b^{(s)}_{k},\mathbf{b}^{(s)}_{2})^{T}$,
where $\mathbf{b}^{(s)}_{1}\in\mathbb{R}^{k-1}$,
$\mathbf{b}^{(s)}_{2}\in\mathbb{R}^{m-k}$, and
$\mathbf{y}^{(s+1)}_{P}=(\mathbf{y}^{(s+1)}_{1},y^{(s+1)}_{\bar{k}},\mathbf{y}^{(s+1)}_{2})^{T}$.
The triangular matrices can be split as

$R^{(s)}=\left(\begin{array}[]{lll}R_{11}&\mathbf{r}_{k}(1:k-1)&R_{12}(1:k-1,:)\\\
0^{T}&\mathbf{r}_{k}(k)&R_{12}(k,:)\\\ 0&0&R_{22}\end{array}\right),$ (31)

and

$R^{(s+1)}=\left(\begin{array}[]{lll}R_{11}&\mathbf{-r}_{k}(1:k-1)&R_{12}(1:k-1,:)\\\
0^{T}&\mathbf{-r}_{k}(k)&R_{12}(k,:)\\\ 0&0&R_{22}\end{array}\right),$ (32)

where $R_{11}$ is upper triangular of order $k-1$, $R_{12}$ has size
$k\times(n_{s}-k)$, and $R_{22}$ is upper triangular of size
$(m-k)\times(n_{s}-k)$.

Since $\mathbf{y}^{(s)}_{2}$ solves $\min\|R_{22}\mathbf{y}-\mathbf{b}_{2}\|$,
we have $\mathbf{y}^{(s)}_{2}=\mathbf{y}^{(s+1)}_{2}$. By backsubstitution, we
obtain

$y^{(s)}_{k}=\frac{1}{\mathbf{r}_{k}(k)}\left(b^{(s)}_{k}-R_{12}(k,1:n_{s}-k)^{T}\mathbf{y}^{(s)}_{2}\right)<0,$
(33)

and similarly

$y^{(s+1)}_{\bar{k}}=-\frac{1}{\mathbf{r}_{k}(k)}\left(b^{(s)}_{k}-R_{12}(k,1:n_{s}-k)^{T}\mathbf{y}^{(s+1)}_{2}\right)=-y^{(s)}_{k}>0,$
(34)

where we used the fact that $\mathbf{y}^{(s+1)}_{2}=\mathbf{y}^{(s)}_{2}$.
Last, we have

$\mathbf{y}^{(s)}_{1}=R_{11}^{-1}\left(\mathbf{b}^{(s)}_{1}-y^{(s)}_{k}\mathbf{r}_{k}(1:k-1)-R_{12}\mathbf{y}^{(s)}_{2}\right),$
(35)

and

$\displaystyle\mathbf{y}^{(s+1)}_{1}$
$\displaystyle=R_{11}^{-1}\left(\mathbf{b}^{(s)}_{1}+y^{(s+1)}_{\bar{k}}\mathbf{r}_{k}(1:k-1)-R_{12}(1:k-1,:)\mathbf{y}^{(s+1)}_{2}\right)$
(36)
$\displaystyle=R_{11}^{-1}\left(\mathbf{b}^{(s)}_{1}-y^{(s)}_{k}\mathbf{r}_{k}(1:k-1)-R_{12}(1:k-1,:)\mathbf{y}^{(s)}_{2}\right)=\mathbf{y}^{(s)}_{1}.$

Therefore,
$\mathbf{y}^{(s+1)}_{P}=(\mathbf{y}^{(s)}_{1},-y^{(s)}_{k},\mathbf{y}^{(s)}_{2})^{T}>0$,
meaning that $\mathbf{y}^{(s+1)}$ is a feasible point and
$\mathbf{s}^{(s+1)}=(\mathbf{y}^{(s+1)}-\mathbf{x}^{(s-1)})$ is a feasible
direction.

Let us show that $\mathbf{s}^{(s+1)}$ is also a descend direction. Suppose now
without loss of generality that $k\leq n$ ($\bar{k}>n$), i.e. $k$ corresponds
to a column of $A$, and let $A=(\mathbf{a}_{1}\ \dots\ \mathbf{a}_{n})$.
Recall that, at the beginning of the inner loop, the vector
$\mathbf{s}^{(s)}=(\mathbf{y}^{(s)}-\mathbf{x}^{(s-1)})$ is a descend
direction, i.e.

$\|\mathbf{r}^{(s)}\|=\|(A\ -A)\mathbf{y}^{(s)}-\mathbf{b}\|<\|(A\
-A)\mathbf{x}^{(s-1)}-\mathbf{b}\|.$

We have

$\displaystyle\mathbf{r}^{(s+1)}$ $\displaystyle=\mathbf{b}-(A\
-A)\mathbf{y}^{(s+1)}$ $\displaystyle=\mathbf{b}-\sum_{\begin{subarray}{c}i\in
P\\\ i\leq
n\end{subarray}}\mathbf{a}_{i}y^{(s+1)}_{i}+\sum_{\begin{subarray}{c}i\in P\\\
i>n,i\neq\bar{k}\end{subarray}}\mathbf{a}_{i}y^{(s+1)}_{i}+\mathbf{a}_{i}y^{(s+1)}_{\bar{k}}\mathbf{a}_{\bar{k}}$
$\displaystyle=\mathbf{b}-\sum_{\begin{subarray}{c}i\in P\\\ i\leq n,i\neq
k\end{subarray}}\mathbf{a}_{i}y^{(s)}_{i}+\sum_{\begin{subarray}{c}i\in P\\\
i>n\end{subarray}}\mathbf{a}_{i}y^{(s)}_{i}+\mathbf{a}_{i}y^{(s)}_{k}\mathbf{a}_{k}$
$\displaystyle=\mathbf{r}^{(s)},$

therefore $\mathbf{s}^{(s+1)}$ is a feasible descend direction with
$\mathbf{y}^{(s+1)}$ feasible. We can set
$\mathbf{x}^{(s+1)}=\mathbf{y}^{(s+1)}$. Notice that once computed
$\mathbf{y}^{(s)}$, no further computation is needed to obtain
$\mathbf{y}^{(s+1)}$, except for the sign flip.

##  5 Numerical experiments

In this section, we present a few results obtained by applying LHDM and other
methods for sparse recovery, to a quite large set of linear systems.

We compared LHDM with LH, both written by us in Matlab, with the other convex
relaxation techniques tested in [21] and with the Orthogonal Matching Pursuit
algorithm (OMP), whose implementation is publicly available online:

  * •

“$l1\\_magic$”, a primal-dual interior-point method for a LP reformulation of
Basis-Pursuit (20), from `www.l1-magic.org`;

  * •

“$l1\\_homotopy$”, from users.ece.gatech.edu/~sasif/ homotopy;

  * •

“$solveBP$”, from SparseLab (sparselab.stanford.edu);

  * •

“$spgl1$”, from www.cs.ubc.ca/labs/scl/spgl1;

  * •

“$isal1$”, from wwwopt.mathematik.tu-darmstadt.de/spear;

  * •

“$solveOMP$”, from SparseLab (sparselab.stanford.edu).

As done in [21], we used the default setting for all optional parameters of
these solvers.

The comparisons are made in terms of solution quality, i.e. the $\ell_{2}$
norm of the error (“distance to optimum”), and execution times. Note that the
distance to optimum is computable because the unicity of the solution is
guaranteed, under the ERC property (22).

Our purpose is not to provide an exhaustive experimental comparison between
the methods but, instead, to assess the substantial performance gain attained
by the Deviation Maximization to the Lawson-Hanson algorithm, i.e. by the LHDM
algorithm, and its competitiveness with the other methods publicly available
for sparse recovery of dense, underdetermined systems. Actually, numerical
experiments reveal that LHDM is a good choice for sparse recovery, for a wide
class of instances.

Numerical tests have been carried out on a dataset freely available online and
fully described in [21]. For brevity, we refer to the cited article for
details about the $548$ matrices that make up the database. It contains many
instances satisfying the ERC (22), thus representing a favorable situation for
compressed sensing in terms of recoverability of sparse solutions via
$\ell_{1}$ minimization. It also includes instances not satisfying the ERC in
order to assess solvers’ performance outside the so-called $\ell_{0}-\ell_{1}$
equivalence. In our test we have used $444$ instances over $548$, since we
restrict our investigation to dense matrices from this dataset. In fact,
dealing with sparse matrices would require a dedicated implementation of LHDM
which can be addressed in the future. In order to simplify the comparison of
the results here presented with the ones in [21], we keep the same indexing
from $1$ to $548$ to show the results. Results concerning sparse instances are
left blank. Instances satisfying the ERC are indexed in two intervals
$[1,200]$ and $[275,474]$. This helps the understanding of the Figures.

We used the positivity trick presented in (23) in order to achieve an
arbitrary signed solution by means of LHDM, which succeeded in exact recovery
for all tests verifying the ERC, meaning that the solution with smallest
support was found.

We used a fixed choice of the thresholds used in LHDM, that is
$\tau_{1}=0.6,\delta=0.9,\tau_{2}=0.15$ and $k_{\max}=32$, that are the same
values used in [10], hence giving evidence of the wide applicability of the
method proposed without the need of choosing problem dependent parameters.
However, a fine tuning of the parameters can give slight improvements in
specific problem families. For example, for very well-conditioned matrices,
$\delta$ can be substantially smaller, to get a more efficient choice of
columns. This fine tuning is out of the scope of this paper since here we want
to exploit the general applicability of the method.

The results are grouped into subsections, each devoted to a single
experimental aspect.

###  5.1 Efficiency of the “sign flip”

We have seen in sec. 4.2 that, in the case of the positivity trick, the inner
loop of LH can be substituted by a much simpler “sign flip” of the involved R
columns. In Figure 1 we see in LH the difference between running the inner
loop and the sign-flip: there is an interesting speedup, around $1.5\times$,
for the more expensive tests.

|  |
---|---|---
Figure 1: Execution times (left), speedup of LH with sign flip vs LH with
internal loop (center) and distance to optimum
$||\mathbf{x}-\mathbf{x}_{calc}||_{2}$ (right), for the Lawson-Hanson (LH)
algorithm, with the classic internal loop and with its replacement of sec. 4.2
(“sign flip”). The vector $\mathbf{x}_{calc}$ is the solution computed by the
algorithm and $\mathbf{x}$ the true solution (unique, in ERC tests).

###  5.2 Speedup of LHDM versus LH

In Figure 2 and 3 we see the speedup of LHDM with respect to LH, for solution
vectors at different levels of sparsity. The original dataset of [21] has
instances with very sparse solution vectors; we created a few modified
datasets, starting from the original, with random generated solutions at
prescribed sparsity levels, and perturbing matrices to satisfy the ERC
property at each instance. As we can notice, the speedup varies a lot with the
sparsity of the solution: the more the solution is less sparse, i.e. the
number of nonzeros increases, the more LHDM accelerates on LH. The Figures
shows the two extreme situations: the original database with very sparse
solutions and a modified database with almost dense solution vectors.

|  |
---|---|---
Figure 2: Execution times (left), speedup of LHDM vs LH (center) and distance to optimum $||\mathbf{x}-\mathbf{x}_{calc}||_{2}$ (right), for the original dataset, whose instances have very sparse solution vectors. The vector $\mathbf{x}_{calc}$ is the solution computed by the algorithm and $\mathbf{x}$ the true solution (unique, in ERC tests). |  |
---|---|---
Figure 3: Execution times (left), speedup of LHDM vs LH (center) and distance
to optimum $||\mathbf{x}-\mathbf{x}_{calc}||_{2}$ (right), for the modified
dataset, whose instances have almost dense solution vectors. The vector
$\mathbf{x}_{calc}$ is the solution computed by the algorithm and $\mathbf{x}$
the true solution (unique, in ERC tests).

###  5.3 Comparison with the other convex relaxation solvers

In Figure 4 we see the comparison of LHDM with the other convex relaxation
solvers listed at the beginning of this section.

|  |
---|---|---
|  |
|  |
|  |
|  |
|  |
Figure 4: Execution times (left), speedup of LHDM vs the other convex
relaxation solvers listed at the beginning of this section (center) and
distance to optimum $||\mathbf{x}-\mathbf{x}_{calc}||_{2}$ (right). Each row
of plots corresponds to a convex relaxation solver, in this order (from the
top): $l1\\_magic$, $l1\\_homotopy$, $solveBP$, $spgl1$, $isal1$ and
$solveOMP$. The vector $\mathbf{x}_{calc}$ is the solution computed by the
algorithm and $\mathbf{x}$ the true solution (unique, in ERC tests).

We can notice that LHDM is a good performer except versus solveOMP (an
implementation of the Orthogonal Matching Pursuit), that surprisingly has not
been considered in [21] for theoretical reasons, but actually outperforms all
the other methods in the whole database considered. In the next subsection we
will therefore go into a detailed comparison between LHDM and solveOMP.

###  5.4 LHDM vs SolveOMP and the DM role

To compare LHDM with SolveOMP we will perturb the database in three main
directions: sparsity level, condition number of the system matrix, and
nonnegativity of the solution.

Let us start with the first database modification. As already did in sec. 5.2,
we created a few modified datasets, starting from the original, with random
generated solutions at prescribed sparsity levels, and perturbing matrices to
satisfy the ERC property at each instance. As we can notice in Figure 5, where
a modified database with almost dense solution vectors is used, LHDM behaves
quite better than SolveOMP now, with a substantial speedup (compare with
Figure 4, where the original database with very sparse solution vectors were
used). If we consider also the speedup gain of LHDM vs LH with the same
modified dataset, shown in Figure 3, we can conclude that the performance gain
of LHDM vs SolveOMP is due to the Deviation Maximization (DM) and not to LH
alone. This is not surprising, since SolveOMP is based on Cholesky
factorization, which is computationally lighter that the QR factorization on
which LH is based. However, the QR is a more stable algorithm and it allows
the DM to improve the performances considerably in a wide range of instances.

|  |
---|---|---
Figure 5: Execution times (left), speedup of LHDM vs SolveOMP (center) and
distance to optimum $||\mathbf{x}-\mathbf{x}_{calc}||_{2}$ (right), for the
modified dataset, whose instances have almost dense solution vectors. The
vector $\mathbf{x}_{calc}$ is the solution computed by the algorithm and
$\mathbf{x}$ the true solution (unique, in ERC tests).

This observation suggests us to try a modified dataset along the second
direction: the matrix condition number. The matrices of the original dataset
in [21] are very well conditioned (condition number $<10^{2}$). We created a
few modified datasets, starting from the original, with random generated
matrices at prescribed condition numbers, and then perturbing the resulting
matrices to satisfy the ERC property at each instance. If the matrix condition
number is moderate, say between $10^{5}$ and $10^{6}$, SolveOMP fails $\approx
70\%$ of instances, while LHDM solves all of them. Moreover, as we can notice
in Figure 6, LHDM keeps its performance gain vs LH, independently from the
conditioning of the problem (compare with Figure 2).

|  |
---|---|---
Figure 6: Execution times (left), speedup of LHDM vs LH (center) and distance
to optimum $||\mathbf{x}-\mathbf{x}_{calc}||_{2}$ (right), for the modified
dataset, whose matrices have a moderate condition number (between $10^{5}$ and
$10^{6}$). The vector $\mathbf{x}_{calc}$ is the solution computed by the
algorithm and $\mathbf{x}$ the true solution (unique, in ERC tests).

The third direction in which we modified the database is the nonnegativity of
the solution. We created one modified dataset, starting from the original, by
taking componentwise the absolute value of each solution vector, and then
perturbing the matrices to satisfy the ERC property at each instance. When a
nonnegative solution is requested, SolveOMP does not satisfy this requirement,
while LHDM does. Moreover, as we can notice in Figure 7, LHDM keeps its
performance gain vs LH.

|  |
---|---|---
Figure 7: Execution times (left), speedup of LHDM vs LH (center) and distance
to optimum $||\mathbf{x}-\mathbf{x}_{calc}||_{2}$ (right), for the modified
dataset, whose solution vectors are componentwise the absolute value of each
solution vector of the original dataset. The vector $\mathbf{x}_{calc}$ is the
solution computed by the algorithm and $\mathbf{x}$ the true solution (unique,
in ERC tests).

##  6 Conclusions

In this work we have presented a new NNLS solver, namely the Lawson-Hanson
with Deviation Maximization algorithm (LHDM), based on the “deviation
maximization” columns selection strategy. Our method relies on cosine
evaluation in order to select a subset of sufficiently linearly independent
columns to compute a feasible descend direction and we have used it as column
pivoting strategy to devise a new active set method. We have provided a
theoretical analysis proving the finite convergence of LHDM, which is shown to
terminate in at most $2^{n}$ steps, just like the standard Lawson-Hanson. This
is a worst case bound, and in practice a polynomial rate of convergence is
observed.

We have given also the details for an efficient implementation of LHDM,
through Householder and Givens orthogonal transformations. In particular, the
algorithmic clue given for the positivity-trick case in sec. 4.2 is by itself
a novel contribution, as far as we know.

Extensive numerical experiments have been carried out over a wide set of
instances, confirming that LHDM yields a significant performance gain over LH
with an average speedup of $3\times$ with peaks up to $6\times$. Moreover, we
compared it with several $\ell_{1}$ solvers whose implementation is publicly
available online. Numerical testing has confirmed that LHDM is competitive
with $\ell_{1}$ solvers for sparse recovery in terms of solution quality and
execution times, revealing it as a good choice for sparse recovery for a wide
class of instances.

The proposed LHDM algorithm exploits BLAS-3 operations for efficiency. A
thorough experimental evidence is not possible in this Matlab implementation
and it requires an implementation in a compiled programming language, e.g. the
C language, which is an ongoing project, for an open-source delivery. With a C
implementation, the performances seen in sec. 5 should sensibly improve, like
we did experience with QRDM [10], included the speedup with respect to LH.

We did not cope with sparse instances, since it would require a dedicated
implementation based on sparse QR with column pivoting and would be
interesting to compare with randomization techniques. Finally, it would be
also interesting to extend the deviation maximization as pivoting strategy to
other problems which require column selection, e.g. more general constrained
optimization problems such as least squares problems with linear inequality
constraints.

## Acknowledgments

The authors gratefully acknowledge the doctoral grant funded by BeanTech
s.r.l. “GPU computing for modeling, nonlinear optimization and machine
learning”. This work was partially supported by the Project BIRD192932 of the
University of Padova.

## References

  * Bischof and Hansen [1992] C. Bischof and P. Hansen.  A block algorithm for computing rank-revealing QR factorizations.  _Numerical Algorithms_ , 2:371–391, 10 1992.  doi: 10.1007/BF02139475.
  * Bischof and Quintana-Ortí [1998a] C. Bischof and G. Quintana-Ortí.  Computing Rank-Revealing QR Factorizations of Dense Matrices.  _ACM Trans. Math. Softw._ , 24:226–253, 06 1998a.  doi: 10.1145/290200.287637.
  * Bischof and Quintana-Ortí [1998b] C. Bischof and G. Quintana-Ortí.  Algorithm 782: Codes for Rank-Revealing QR Factorizations of Dense Matrices.  _ACM Transactions on Mathematical Software_ , 24:254–257, 07 1998b.  doi: 10.1145/290200.287638.
  * Bischof [1989] J. R. Bischof.  A block QR factorization algorithm using restricted pivoting.  In _Supercomputing ’89:Proceedings of the 1989 ACM/IEEE Conference on Supercomputing_ , pages 248–256, 1989.  doi: 10.1145/76263.76290.
  * Bro and Jong [1997] R. Bro and S. Jong.  A Fast Non-negativity-constrained Least Squares Algorithm.  _Journal of Chemometrics_ , 11:393–401, 09 1997.  doi: 10.1002/(SICI)1099-128X(199709/10)11:53.0.CO;2-L.
  * Bruckstein et al. [2008] A. M. Bruckstein, M. Elad, and M. Zibulevsky.  On the uniqueness of nonnegative sparse solutions to underdetermined systems of equations.  _IEEE Transactions on Information Theory_ , 54(11):4813–4820, 2008.  doi: 10.1109/TIT.2008.929920.
  * Candès [2008] E. J. Candès.  The restricted isometry property and its implications for compressed sensing.  _Comptes Rendus Mathematique_ , 346(9):589–592, 2008.  ISSN 1631-073X.  doi: https://doi.org/10.1016/j.crma.2008.03.014.  URL https://www.sciencedirect.com/science/article/pii/S1631073X08000964.
  * Chen and Donoho [1994] S. Chen and D. Donoho.  Basis pursuit.  In _Proceedings of 1994 28th Asilomar Conference on Signals, Systems and Computers_ , volume 1, pages 41–44 vol.1, 1994.  doi: 10.1109/ACSSC.1994.471413.
  * Chen et al. [2001] S. S. Chen, D. L. Donoho, and M. A. Saunders.  Atomic decomposition by basis pursuit.  _SIAM Rev._ , 43(1):129–159, Jan. 2001.  ISSN 0036-1445.  doi: 10.1137/S003614450037906X.  URL https://doi.org/10.1137/S003614450037906X.
  * Dessole and Marcuzzi [2022] M. Dessole and F. Marcuzzi.  Deviation maximization for rank-revealing qr factorizations.  _Numerical Algorithms_ , 2022.  doi: 10.1007/s11075-022-01291-1.  URL https://doi.org/10.1007/s11075-022-01291-1.
  * Dessole et al. [2020] M. Dessole, F. Marcuzzi, and M. Vianello.  Accelerating the Lawson-Hanson NNLS solver for large-scale Tchakaloff regression designs.  _Dolomites Research Notes on Approximation_ , 13:20 – 29, 2020.  ISSN 2035-6803.  doi: http://dx.doi.org/10.14658/PUPJ-DRNA-2020-1-3.
  * Donoho [2006] D. Donoho.  Compressed sensing.  _Information Theory, IEEE Transactions on_ , 52:1289 – 1306, 05 2006.  doi: 10.1109/TIT.2006.871582.
  * Donoho and Huo [2001] D. Donoho and X. Huo.  Uncertainty principles and ideal atomic decomposition.  _IEEE Transactions on Information Theory_ , 47(7):2845–2862, 2001.  doi: 10.1109/18.959265.
  * Donoho and Elad [2003] D. L. Donoho and M. Elad.  Optimally sparse representation in general (nonorthogonal) dictionaries via $\ell_{1}$ minimization.  _Proceedings of the National Academy of Sciences_ , 100(5):2197–2202, 2003.
  * Elad [2010] M. Elad.  _Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing_.  Springer Publishing Company, Incorporated, 1st edition, 2010.  ISBN 144197010X.
  * Elad and Bruckstein [2002] M. Elad and A. M. Bruckstein.  A generalized uncertainty principle and sparse representation in pairs of bases.  _IEEE Trans. Inform. Theory_ , 48:2558–2567, 2002.
  * Foucart and Koslicki [2014] S. Foucart and D. Koslicki.  Sparse Recovery by Means of Nonnegative Least Squares.  _Signal Processing Letters, IEEE_ , 21:498–502, 04 2014\.  doi: 10.1109/LSP.2014.2307064.
  * Gillis [2020] N. Gillis.  _Nonnegative Matrix Factorization_.  SIAM, 2020.
  * Lawson and Hanson [1995] C. L. Lawson and R. J. Hanson.  _Solving least squares problems_ , volume 15.  SIAM, 1995.
  * Lee and Seung [1999] D. D. Lee and H. S. Seung.  Learning the parts of objects by non-negative matrix factorization.  _Nature_ , 401(6755):788–791, 1999.  doi: 10.1038/44565.  URL https://doi.org/10.1038/44565.
  * Lorenz et al. [2015] D. A. Lorenz, M. E. Pfetsch, and A. M. Tillmann.  Solving basis pursuit: Heuristic optimality check and solver comparison.  _ACM Trans. Math. Softw._ , 41(2), Feb. 2015.  ISSN 0098-3500.  doi: 10.1145/2689662.  URL https://doi.org/10.1145/2689662.
  * Portugal et al. [1994] L. F. Portugal, J. J. Júdice, and L. N. Vicente.  A Comparison of Block Pivoting and Interior-Point Algorithms for Linear Least Squares Problems with Nonnegative Variables.  _Mathematics of Computation_ , 63(208):625–643, 1994.  ISSN 00255718, 10886842.
  * Quintana-Ortí et al. [1998] G. Quintana-Ortí, X. Sun, and C. H. Bischof.  A BLAS-3 Version of the QR Factorization with Column Pivoting.  _SIAM Journal on Scientific Computing_ , 19(5):1486–1494, 1998.  doi: 10.1137/S1064827595296732.
  * Radons [2016] M. Radons.  Direct solution of piecewise linear systems.  _Theoretical Computer Science_ , 626:97–109, 2016.  ISSN 0304-3975.  doi: https://doi.org/10.1016/j.tcs.2016.02.009.
  * Slawski and Hein [2013] M. Slawski and M. Hein.  Non-negative least squares for high-dimensional linear models: Consistency and sparse recovery without regularization.  _Electronic Journal of Statistics_ , 7(1):3004–3056, 2013.  URL www.scopus.com.  Cited By :81.
  * Stoer [1971] J. Stoer.  On the Numerical Solution of Constrained Least-Squares Problems.  _SIAM Journal on Numerical Analysis_ , 8(2):382–411, 1971.  ISSN 00361429.
  * Tibshirani [1996] R. Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 58(1):267–288, 1996.  ISSN 00359246.  URL http://www.jstor.org/stable/2346178.
  * Tillmann and Pfetsch [2014] A. M. Tillmann and M. E. Pfetsch.  The computational complexity of the restricted isometry property, the nullspace property, and related concepts in compressed sensing.  _IEEE Transactions on Information Theory_ , 60(2):1248–1259, 2014.  doi: 10.1109/TIT.2013.2290112.
  * Tropp [2004] J. Tropp.  Greed is good: algorithmic results for sparse approximation.  _IEEE Transactions on Information Theory_ , 50(10):2231–2242, 2004.  doi: 10.1109/TIT.2004.834793.
  * Van Benthem and Keenan [2004] M. H. Van Benthem and M. R. Keenan.  Fast algorithm for the solution of large-scale non-negativity-constrained least squares problems.  _Journal of Chemometrics_ , 18(10):441–450, 2004\.  doi: https://doi.org/10.1002/cem.889.  URL https://analyticalsciencejournals.onlinelibrary.wiley.com/doi/abs/10.1002/cem.889.
  * Wang and Tang [2009] M. Wang and A. Tang.  Conditions for a unique non-negative solution to an underdetermined system.  In _Proceedings of the 47th Annual Allerton Conference on Communication, Control, and Computing_ , Allerton’09, page 301–307. IEEE Press, 2009.  ISBN 9781424458707.
  * Wang et al. [2011] M. Wang, W. Xu, and A. Tang.  A Unique “Nonnegative” Solution to an Underdetermined System: From Vectors to Matrices.  _IEEE Transactions on Signal Processing_ , 59(3):1007–1016, Mar. 2011.  doi: 10.1109/TSP.2010.2089624.