Safe Feature Elimination for Non-Negativity Constrained Convex Optimization

TL;DR

This paper introduces a safe feature elimination method for non-negativity constrained convex optimization, enabling the removal of guaranteed zero features using primal-dual feasible pairs, improving efficiency and robustness especially in ill-conditioned problems.

Contribution

The paper develops a novel safe feature elimination strategy for non-negativity constrained convex problems, utilizing primal-dual feasible pairs to enhance robustness and applicability.

Findings

Successfully eliminates zero features in ill-conditioned problems

Certifies the uniqueness of non-negative least-squares solutions

Demonstrates effectiveness on synthetic and microscopy datasets

Abstract

Inspired by recent work on safe feature elimination for $1$-norm regularized least-squares, we develop strategies to eliminate features from convex optimization problems with non-negativity constraints. Our strategy is safe in the sense that it will only remove features/coordinates from the problem when they are guaranteed to be zero at a solution. To perform feature elimination we use an accurate, but not optimal, primal-dual feasible pair, making our methods robust and able to be used on ill-conditioned problems. We supplement our feature elimination problem with a method to construct an accurate dual feasible point from an accurate primal feasible point; this allows us to use a first-order method to find an accurate primal feasible point, then use that point to construct an accurate dual feasible point and perform feature elimination. Under reasonable conditions, our feature elimination strategy will eventually eliminate all zero features from the problem. As an application of our methods we show how safe feature elimination can be used to robustly certify the uniqueness of non-negative least-squares (NNLS) problems. We give numerical examples on a well-conditioned synthetic NNLS problem and a on set of 40000 extremely ill-conditioned NNLS problems arising in a microscopy application.

Tables

Table 1. Table 1: Number of problems (40000 total) certified to have unique solutions using feature elimination with the strong concavity subproblem ( 10 ) and with the strong concavity plus partial dual feasibility subproblem ( 19 ). Adding partial dual feasibility constraints in ( 19 ) can eliminate sufficiently many features to certify solution uniqueness at fewer iterations than using strong concavity alone in ( 10 ).

	Solutions Certified Unique		Features Eliminated
Iterations	using (10)	using (19)	using (10)	using (19)
$10,000$	$7237$	$9510$	$20.0\%$	$24.8\%$
$50,000$	$18,339$	$23,174$	$46.2\%$	$56.9\%$
$100,000$	$23,512$	$28,826$	$57.9\%$	$68.3\%$
$500,000$	$34,462$	$38,094$	$81.6\%$	$87.8\%$
$500,000$ + lsqnonneg	$-$	$40,000$	$-$	$91.0\%$

Equations

$$f^{\ast }(y)\stackrel{\mathrm{def}}{=}\underset{x}{\sup }⟨y,x⟩-f(x).$$

$$\begin{array}[]{ll}\min_{x}&f(Ax)\\ \text{s.t.}&x\geq 0.\end{array}$$

$$p^{\ast }=\underset{x}{\min }f(Ax)+{\iota }_{{\mathbb{R}}_{+}^n}(x)$$

$$\underset{\nu }{\max }-f^{\ast }(\nu )-({\iota }_{{\mathbb{R}}_{+}^n}{)}^{\ast }(-A^T\nu ).$$

$$\begin{array}[]{lll}d^{\ast}=&\max_{\nu}&g(\nu)\\ &\text{s.t.}&A^{T}\nu\geq 0,\end{array}$$

$$\frac{1}{2L}\Vert \nu -{\nu }^{\ast }{\Vert }^2\le g({\nu }^{\ast })-g(\nu )\forall \nu \in \mathrm{dom}g.$$

$$A^T\nabla f(Ax)-A^T\nu$$

$$x$$

$$A^T\nu$$

$$x_i\{A^T\nu {\}}_i$$

$$\begin{array}[]{ll}\min_{\nu}&\langle a_{i},\nu\rangle\\ \text{s.t.}&\nu\in N.\end{array}$$

$$\frac{1}{2L}\Vert \hat{\nu }-{\nu }^{\ast }{\Vert }^2\le g({\nu }^{\ast })-g(\hat{\nu }).$$

$$g({\nu }^{\ast })-g(\hat{\nu })=f(A{x}^{\ast })-g(\hat{\nu })\le f(A\hat{x})-g(\hat{\nu })=ϵ.$$

$$\begin{array}[]{ll}\min_{\nu}&\langle a_{i},\nu\rangle\\ \text{s.t.}&\|\nu-\hat{\nu}\|^{2}\leq 2L\epsilon.\end{array}$$

$$\begin{array}[]{ll}\text{min}_{\hat{\nu}}&\dfrac{1}{2}\|\hat{\nu}-\nu^{\prime}\|^{2}\\ \text{s.t.}&A^{T}\hat{\nu}\geq 0,\end{array}$$

$$\begin{array}[]{ll}\min&t\\ \text{s.t.}&A^{T}((1-t)\nu^{\prime}+t\nu^{\text{strict}})\geq 0\\ &0\leq t\leq 1.\end{array}$$

$$\begin{array}[]{ll}\max_{\nu,t}&t\\ \text{s.t.}&A^{T}\nu\geq t\\ &\mathbbm{1}^{T}A^{T}\nu=1.\end{array}$$

$$t(\lambda;\lambda_{0})\stackrel{{\scriptstyle\operatorname{def}}}{{=}}\left\{\begin{array}[]{ll}0&\lambda\geq 0\\ \dfrac{\lambda}{\lambda-\lambda_{0}}&\lambda<0,\end{array}\right.$$

$$t^{\ast }=\underset{i}{\max }t(a_i^T{\nu }^{\prime };a_i^T{\nu }^{\text{strict}}).$$

$$\hat{\nu }=(1-t^{\ast }){\nu }^{\prime }+t^{\ast }{\nu }^{\text{strict}}.$$

$$\Vert {\hat{\nu }}^k-{\nu }^{\ast }\Vert \le \Vert {\hat{\nu }}^k-{\nu }^{\prime k}\Vert +\Vert {\nu }^{\prime k}-{\nu }^{\ast }\Vert .$$

$$\Vert {\hat{\nu }}^k-{\nu }^{\prime k}\Vert =\Vert (1-t^k){\nu }^{\prime k}+t^k{\nu }^{\text{strict}}-{\nu }^{\prime k}\Vert =t^k\Vert {\nu }^{\prime k}-{\nu }^{\text{strict}}\Vert ,$$

$$\nu^{\prime}\mapsto t(a_{i}^{T}\nu^{\prime};a_{i}^{T}\nu^{\text{strict}})=\left\{\begin{array}[]{ll}0&a_{i}^{T}\nu^{\prime}\geq 0\\ \dfrac{a_{i}^{T}\nu^{\prime}}{a_{i}^{T}\nu^{\prime}-a_{i}^{T}\nu^{\text{strict}}}&a_{i}^{T}\nu^{\prime}<0,\end{array}\right.$$

$$\sqrt{ϵ}<\frac{1}{\sqrt{2L}}\underset{i\in \mathcal{I}}{\min }\frac{⟨a_i,\hat{\nu }⟩}{\Vert a_i\Vert }.$$

$$⟨a_i,\hat{\nu }⟩-\sqrt{2Lϵ}\Vert a_i\Vert >0⟺\sqrt{ϵ}<\frac{1}{\sqrt{2L}}\frac{⟨a_i,\hat{\nu }⟩}{\Vert a_i\Vert }.$$

$$\begin{array}[]{ll}\min_{x}&\frac{1}{2}\|Ax-b\|^{2}\\ \text{s.t.}&x\geq 0.\end{array}$$

$$\begin{array}[]{ll}\max_{\nu}&g(\nu)=-\frac{1}{2}\|\nu+b\|^{2}+\frac{1}{2}\|b\|^{2}\\ \text{s.t.}&A^{T}\nu\geq 0.\end{array}$$

$$\begin{array}[]{ll}\min_{\nu}&\langle a_{i},\nu\rangle\\ \text{s.t.}&\|\nu-\hat{\nu}\|^{2}\leq 2\epsilon,\end{array}$$

$$A=\left[\begin{array}{ccccc}1& 6& -1& 8& 0\\ -2& 7& 1& 8& 2\\ 3& 1& 4& 1& -5\end{array}\right],$$

$$x^{+}\leftarrow x-t{A}^T(Ax-b)$$

$$A_{\text{red}}=\left[\begin{array}{ccc}1& -1& 0\\ -2& 1& 2\\ 3& 4& -5\end{array}\right],$$

$$\Vert \hat{x}-x^{\ast }{\Vert }^2\le \frac{2}{{\sigma }_{\text{min}}(A_{\text{red}}{)}^2}\left(f(A\hat{x})-f(A{x}^{\ast })\right)\le \frac{2}{{\sigma }_{\text{min}}(A_{\text{red}}{)}^2}ϵ\doteq \mathrm{0.066.}$$

Full text

This is a pre-print of an article published in Journal of Optimization Theory and Applications. The final authenticated version is available online at: https://doi.org/10.1007/s10957-019-01612-w.

1 Introduction

There is an expanding body of work on safe feature elimination for $1$-norm regularized optimization problems, particularly for $1$-norm regularized least-squares (the lasso). Safe feature elimination removes features/columns of the dictionary/observation matrix when they are guaranteed not to be present in a solution. El Ghaoui et al.’s influential work in this direction [GVR12] is based on using complementary slackness between primal and dual optimization problems to identify zero coordinates in a solution to the primal problem. Complementary slackness implies that if the dual optimal point satisfies an inequality constraint strictly, then the corresponding primal optimal coordinate must be equal to zero in any primal optimal point (we will make this statement precise shortly).

Using duality to identify zero coordinates has been used before, of course; for instance duality has been used to eliminate features in linear programs [TTZ66]. What is novel is that safe feature elimination (SAFE) strategies are designed to avoid the use of the exact dual optimal point, which may be very expensive to compute in practice. SAFE strategies instead use an auxiliary dual feasible point to construct a compact set that is guaranteed to contain the dual optimal point. If all points in this compact set satisfy a dual inequality constraint strictly, then the exact dual optimal point also satisfies the inequality strictly and we can safely eliminate the corresponding primal coordinate. We refer the reader to [FGS15] for a discussion of numerous such safe sets for the lasso and [XWR17] for a survey on both safe and unsafe feature elimination strategies for lasso problems.

In this paper we develop a SAFE strategy for non-negativity constrained convex optimization problems which uses an accurate, but non-optimal, primal-dual feasible pair. This is similar to the SAFE strategy for the lasso proposed in [FGS15] which is more robust than El Ghaoui et al.’s original in [GVR12]. We show that under reasonable conditions, a sufficiently accurate primal-dual pair will eliminate all zero coordinates from the problem.

A recent technique in super-resolution fluorescence microscopy uses many tens of thousands of non-negative least-squares (NNLS) problems to form a super-resolved image. Motivated by these problems we focus our efforts on the case where only an accurate primal feasible point is known, as is usually the case when using first-order methods to solve the primal. To enable the use of SAFE we propose an efficient method to construct an accurate dual feasible point from a given primal feasible point. We also show that the construction depends continuously on the given primal feasible point, meaning that as the primal feasible point converges to an optimal point (e.g., as one iterates a first-order method) so too does the dual, and hence SAFE will eliminate all zero features.

We apply our SAFE strategy to the task of robustly certifying the uniqueness of solutions to NNLS problems. In a small synthetic numerical example we compare our method with an existing uniqueness sufficient condition that relies on a strong assumption on the structure of the data matrix. The strong assumption is that the columns are in general linear position (18), which can be checked only for very small matrices or if the matrix has a generating model of a certain form. In a real-data numerical example of a much larger size we use our SAFE strategy to certify the uniqueness of reconstructed images from a microscopy application. It is infeasible to check if the columns of the data matrix are in general linear position, so the existing uniqueness condition cannot be used. We instead find an approximate solution (a reconstructed image) to the NNLS problem using an efficient gradient method and use our SAFE strategy to confirm that the exact reconstructed image is unique.

Constructing the dual feasible point and performing feature elimination costs about as much as a primal gradient evaluation, which is to say that it is not expensive. Although we do not explore this direction in this work, the inexpensiveness of our SAFE strategy likely allows it to be used to decrease the cost of solving the primal problem with a first-order method, as has already been demonstrated for the lasso [XWR17].

The rest of this paper is organized as follows. In Section 2 we state a general non-negativity constrained primal problem, develop a dual problem, and state KKT optimality conditions. Section 3 derives the general structure of a SAFE strategy using an accurate primal-dual feasible pair. We also give a simple, but effective, instantiation of this strategy. To enable SAFE to work with first-order methods, Section 4 derives and analyses a dual line search that allows us to construct an accurate dual feasible point from an accurate primal feasible point. Section 5 gives the proof that our SAFE strategy eventually eliminates all zero features. Sections 6 and 7 discuss robustly certifying solution uniqueness for NNLS problems.

2 Preliminaries

Let $\overline{\mathbb{R}}\stackrel{{\scriptstyle\operatorname{def}}}{{=}}\mathbb{R}\cup\{\pm\infty\}$ be the set of extended real values and $\mathbb{R}^{n}_{+}=\{x\in\mathbb{R}^{n}\,:\,x\geq 0\}$ be the non-negative orthant. We denote the standard inner product on $\mathbb{R}^{n}$ by both $\langle x,y\rangle$ and $x^{T}y$, and the induced norm by $\|x\|=\sqrt{x^{T}x}$. The convex conjugate of a function $f$ is defined via

[TABLE]

We denote domain of a function $f$ by $\operatorname{dom}{f}$ and the $i$th coordinate of a vector $x$ by $x_{i}$ or $\{x\}_{i}$. Except for the convex conjugate, we use a superscript $\ast$ to denote the value of a quantity at an optimum, e.g., $p^{\ast}$ for the optimal value of a primal optimization problem.

We consider a general optimization problem involving a convex objective $f$ subject to a non-negativity constraint on the optimization variables:

[TABLE]

For example, this generic problem structure captures non-negative least-squares (NNLS) with ${f(z)=\frac{1}{2}\|z-b\|^{2}}$. We assume the following throughout the paper:

• •

$f:\mathbb{R}^{m}\to\overline{\mathbb{R}}$ is a proper, convex, extended value function

• •

$f$ has globally $L$-Lipschitz continuous gradient with $L>0$

• •

$A$ is a real ${m\times n}$ matrix with full rank

• •

$\operatorname{int}{\operatorname{dom}f}\cap A\mathbb{R}^{n}_{+}$ is non-empty

We begin by deriving a dual problem to (1). The problem (1) fits naturally into the framework of Fenchel-Rockafellar duality (though one can use Lagrange duality to find the same dual problem; see [BV04]). See the text [Bec17] for a self-contained introduction or [BC17] for a thorough treatment. We can write (1) as

[TABLE]

where $\iota_{\mathbb{R}^{n}_{+}}$ is the indicator function for the non-negative orthant, i.e. $\iota_{\mathbb{R}^{n}_{+}}(x)$ is [math] if $x\geq 0$ and ${+\infty}$ otherwise.

We can directly write the dual problem as

[TABLE]

The conjugate of $\iota_{\mathbb{R}^{n}_{+}}$ is readily found to be $\iota_{\mathbb{R}^{n}_{-}}$, the indicator for the non-positive orthant. Writing this as a constraint on $A^{T}\nu$, we have the dual problem

[TABLE]

where we have defined the dual objective ${g(\nu)\stackrel{{\scriptstyle\operatorname{def}}}{{=}}-f^{\ast}(\nu)}$. Since we assumed $f$ to be proper and convex, $f^{\ast}$ is proper and convex. Further, by assuming $f$ has $L$-Lipschitz continuous gradient, the “conjugate correspondence theorem” (Theorem 5.26 of [Bec17]) implies that $g=-f^{\ast}$ is $1/L$-strongly concave. The strong concavity of $g$ implies that the dual optimal point $\nu^{\ast}$ exists and is unique. It also provides us with the bound

[TABLE]

We will use this bound as a fundamental building block for our feature elimination procedure in Section 3.

Observe that Slater’s condition holds for the primal problem (1). This implies that strong duality holds, so that the primal optimal value $p^{\ast}$ and the dual optimal value $d^{\ast}$ are equal. Note that Slater’s condition holding for the primal problem also shows that the dual optimal value is attained (a $\nu^{\ast}$ exists that achieves $d^{\ast}=g(\nu^{\ast})$) [BV04]. Furthermore, Slater’s condition holds for the dual problem (2), which implies that the primal optimal value is attained and that the KKT conditions are necessary and sufficient for primal and dual optimal points.

The KKT conditions for the primal problem (1) can be written as

[TABLE]

3 Safe Feature Elimination

Let $x^{\ast}$ be a (primal) optimal point of (1) and $\nu^{\ast}$ the (dual) optimal point of (2). Let $a_{i}$ be the $i$th column/feature of the matrix $A$. The complementary slackness condition (7) implies that if ${\{A^{T}\nu^{\ast}\}_{i}>0}$, then $x_{i}^{\ast}=0$. The key idea of safe feature elimination is that if we can certify that ${\{A^{T}\nu^{\ast}\}_{i}=\langle a_{i},\nu^{\ast}\rangle>0}$, then we can guarantee that $x_{i}^{\ast}=0$. This allows us to eliminate the $i$th column of $A$, $a_{i}$, from the problem with a guarantee that it will not be present in a solution. What remains is to robustly determine for each column $a_{i}$ if ${\langle a_{i},\nu^{\ast}\rangle>0}$ without knowledge of the exact solutions $\nu^{\ast}$ or $x^{\ast}$.

Observe that we do not require the precise value of ${\langle a_{i},\nu^{\ast}\rangle}$; we merely need to certify that ${\langle a_{i},\nu^{\ast}\rangle}$ is strictly positive to certify $x_{i}^{\ast}=0$. This allows us to avoid the apparent need for the exact solution $\nu^{\ast}$. Suppose we have a set of dual points $N$ that is guaranteed to contain $\nu^{\ast}$. We then find a lower bound for ${\langle a_{i},\nu^{\ast}\rangle}$ by solving the “feature elimination subproblem”

[TABLE]

That $N$ contains $\nu^{\ast}$ makes (8) safe. The optimal value of (8) is guaranteed to be no larger than ${\langle a_{i},\nu^{\ast}\rangle}$, so that if the optimal value is strictly positive we can certify that ${\langle a_{i},\nu^{\ast}\rangle>0}$. The feature elimination subproblem (8) tests for elimination of the single feature $a_{i}$, so to test for feature elimination on all of $A$ we simply solve (8) for each column of $A$.

We now construct a simple, but effective, search set $N$ without the use of any exact solutions. Let us assume that we have access to both a primal feasible point $\hat{x}$ and a dual feasible point $\hat{\nu}$, neither of which are assumed to be optimal. This gives the duality gap ${\epsilon=f(A\hat{x})-g(\hat{\nu})}$. Since strong duality holds, the duality gap $\epsilon$ will shrink to zero as $\hat{x}$ and $\hat{\nu}$ become increasingly accurate. Using the strong-concavity bound (3), we have

[TABLE]

Since strong duality holds, ${g(\nu^{\ast})=f(Ax^{\ast})}$, from which we see

[TABLE]

Combining these gives us a bound on the distance from $\hat{\nu}$ to $\nu^{\ast}$ in terms of the duality gap $\epsilon$. We therefore define the search set $N$ to be the set of all points satisfying this bound: $N\stackrel{{\scriptstyle\operatorname{def}}}{{=}}{\{\nu\,:\,\|\hat{\nu}-\nu\|^{2}\leq 2L\epsilon\}}$. As desired, $N$ is guaranteed to contain $\nu^{\ast}$, but is constructed using only the feasible points $\hat{x}$ and $\hat{\nu}$.

The associated feature elimination subproblem is

[TABLE]

The problem (10) has a linear objective and the constraint set is a ball of radius $\sqrt{2L\epsilon}$ centered at $\hat{\nu}$. See Figure 1 for a diagram of the dual geometry for this problem. Taking the search set $N$ to be a ball as we have done is very similar to the GAP SAFE sphere test of [FGS15].

The optimal value is easily found in closed-form to be ${\langle a_{i},\hat{\nu}\rangle-\sqrt{2L\epsilon}\|a_{i}\|}$. As $\hat{x}$ and $\hat{\nu}$ become more accurate, $\hat{\nu}$ approaches $\nu^{\ast}$ and $\epsilon$ shrinks to zero; as this occurs the optimal value of the subproblem approaches $\langle a_{i},\nu^{\ast}\rangle$, giving more precise lower bounds and thus increasing the strength of the subproblem to eliminate features. It is therefore crucial to have accurate, feasible $\hat{x}$ and $\hat{\nu}$ in order to apply (10) effectively.

4 Dual Line Search

In order to use feature elimination subproblem (10), we must have a primal feasible $\hat{x}$ and a dual feasible $\hat{\nu}$ that achieve a reasonably tight duality gap $\epsilon=f(\hat{x})-g(\hat{\nu})$. When using a first-order method on the primal we have access to an accurate primal feasible point $\hat{x}$ simply by taking one of the iterates. But we typically do not have access to an accurate dual feasible point $\hat{\nu}$. Hence we derive an inexpensive method to find an accurate dual feasible $\hat{\nu}$ from an accurate primal feasible $\hat{x}$.

4.1 Finding an Accurate Dual Feasible $\hat{\nu}$ From an Accurate Primal Feasible $\hat{x}$

To leverage the accuracy of $\hat{x}$, we form ${\nu^{\prime}=\nabla f(A\hat{x})}$, since if $\hat{x}$ were optimal, then $\nabla f(A\hat{x})$ would be the dual optimal point (see Lemma 4.2). But note that $\nu^{\prime}$ is not guaranteed to be dual feasible since $\hat{x}$ is not necessarily optimal (i.e., $A^{T}\nu^{\prime}\not\geq 0$ is possible). To fix this, perhaps the “best” approach is to solve the orthogonal projection problem

[TABLE]

which finds the closest dual feasible point to $\nu^{\prime}$. But this projection subproblem is expensive to solve, as it is closely related to the dual problem (2). We do not need the “best” possible $\hat{\nu}$, however. We only need $\hat{\nu}$ that does not spoil the accuracy provided by $\hat{x}$, thereby providing a small duality gap $\epsilon$.

To that end let us assume we have access to a strictly dual feasible point $\nu^{\text{strict}}$. We can use $\nu^{\text{strict}}$ to construct $\hat{\nu}$ nearby $\nu^{\prime}$ that is also dual feasible using a simple line search: we find the closest dual feasible point to $\nu^{\prime}$ along the line segment between $\nu^{\prime}$ and $\nu^{\text{strict}}$ via

[TABLE]

Once we have solved the line search for $t^{\ast}$, we form $\hat{\nu}=(1-t^{\ast})\nu^{\prime}+t^{\ast}\nu^{\text{strict}}$. We can view this line search as a not-necessarily-orthogonal projection onto the dual feasible set. We will give a few simple methods to find a strictly dual feasible $\nu^{\text{strict}}$ in Subsection 4.1.1, and in Subsection 4.1.2 we will show that $\nu^{\text{strict}}$ being strictly dual feasible (instead of just dual feasible) is necessary for $\hat{\nu}$ from the line search to converge to $\nu^{\ast}$ as $\hat{x}$ converges to $x^{\ast}$.

See Figure 2 for a diagram of this line search in two dimensions. The boundary of the dual feasible set is given by two hyperplanes ${\langle a_{1},\nu\rangle=0}$ and ${\langle a_{2},\nu\rangle=0}$. We can see ${\langle a_{2},\nu^{\prime}\rangle<0}$, so $\nu^{\prime}$ is not dual feasible.

The constraint $t\geq 0$ is used in the line search only so that $\hat{\nu}=\nu^{\prime}$ in the case when $\nu^{\prime}$ is already dual feasible. Additionally, the optimal value is never greater than $1$, since the point $\nu^{\text{strict}}$ is assumed to be dual feasible. By precomputing $A^{T}\nu^{\prime}$ and $A^{T}\nu^{\text{strict}}$, the optimal value of the line search and the resulting dual feasible point $\hat{\nu}$ can be found in closed-form, which is given in Subsection 4.1.2.

4.1.1 Finding a Strictly Dual Feasible $\nu^{\text{strict}}$

For a given $A$, we can search for a strictly dual feasible point $\nu^{\text{strict}}$ via the linear program (LP)

[TABLE]

This problem maximizes a lower bound of $A^{T}\nu$ while the constraint ${\mathbbm{1}^{T}A^{T}\nu=1}$ serves to keep $\nu$ bounded. Note that this linear program may fail to find a strictly dual feasible $\nu^{\text{strict}}$. But it will do so only when $A$ does not admit any strictly feasible points (e.g., $A=0$ does not does not admit any strictly dual feasible points). The cost of this LP is not much of a concern (unless $A$ is huge), as it only needs to be solved once to find $\nu^{\text{strict}}$; once we have a strictly dual feasible $\nu^{\text{strict}}$, we can use it for any primal problem of the form (1) with the same $A$.

There are other methods to find a suitable $\nu^{\text{strict}}$ without solving the LP (13). If the sum of each row of $A$ is positive, the point ${\nu^{\text{strict}}=\mathbbm{1}}$ is strictly dual feasible. In particular, if $A$ is elementwise positive (as is the case in our microscopy example in Subsection 7.2), ${\nu^{\text{strict}}=\mathbbm{1}}$ is strictly dual feasible. We can also take $\nu^{\text{strict}}=\max\{0,\nu^{\prime}\}$, where $\nu^{\prime}=\nabla f(A\hat{x})$; if $\nu^{\text{strict}}$ has at least one positive entry, then it is strictly dual feasible. In our microscopy example $A$ is elementwise positive and we find that using ${\nu^{\text{strict}}=\max\{0,\nu^{\prime}\}}$ reliably produces strictly dual feasible points (and avoids the need for solving the LP (13)).

4.1.2 The Dual Line Search is a Continuous Mapping

In Subsection 4.2 we will show how $\hat{\nu}$ from the dual line search converges to the dual optimal point $\nu^{\ast}$ as $\hat{x}$ converges to a primal optimal point. This will then be used to show that, under reasonable conditions, our dual line search and feature elimination strategy will eventually eliminate all zero features from the problem. This means that if we perform sufficiently many iterations of a first-order method, we can eliminate all zero features from the problem. To enable that analysis, we find a closed-form solution to the line search and prove a lemma on the continuity of the mapping from $\nu^{\prime}=\nabla f(A\hat{x})$ to $\hat{\nu}$ found via the line search.

We find the closed-form solution to the line search by identifying two cases:

If $\nu^{\prime}$ is dual feasible, $t=0$ is the minimum feasible value, which leads to $\hat{\nu}=\nu^{\prime}$. 2. 2.

Otherwise, there is at least one index $i$ such that ${\langle a_{i},\nu^{\prime}\rangle=\{A^{T}\nu^{\prime}\}_{i}<0}$. In this case, we must increase $t$ until the all coordinates of $A^{T}((1-t)\nu^{\prime}+t\nu^{\text{strict}})$ are non-negative.

We define the scalar-valued function

[TABLE]

where $\lambda$ is the scalar independent variable and $\lambda_{0}$ is a fixed parameter. We can write the dual feasible point returned from the line search as $\hat{\nu}=(1-t^{\ast})\nu^{\prime}+t^{\ast}\nu^{\text{strict}}$ where

[TABLE]

Lemma 4.1**.**

If $\nu^{\text{strict}}$ is strictly dual feasible (i.e., $A^{T}\nu^{\text{strict}}>0$), then the dual line search (12) mapping $\nu^{\prime}$ to $\hat{\nu}$ is continuous in $\nu^{\prime}$.

Proof The strict dual feasibility assumption states that $a_{i}^{T}\nu^{\text{strict}}>0$ for each $i$. The dual line search produces the point

[TABLE]

To show continuity of the mapping $\nu^{\prime}\mapsto\hat{\nu}$, it is sufficient to show $t^{\ast}$ is continuous in $\nu^{\prime}$.

Observe that if $\lambda_{0}>0$ the function $t(\lambda;\lambda_{0})$ is continuous for all $\lambda$. Since we take $\nu^{\text{strict}}$ strictly dual feasible, $a_{i}^{T}\nu^{\text{strict}}>0$ for each $i$, meaning that $t(a_{i}^{T}\nu^{\prime};a_{i}^{T}\nu^{\text{strict}})$ depends continuously on $\nu^{\prime}$ for each $i$. Since $t^{\ast}$ is the pointwise maximum of continuous functions of $\nu^{\prime}$, it is continuous in $\nu^{\prime}$, completing the proof.∎

The strict dual feasibility assumption in Lemma 4.1 is necessary for the continuity of the mapping. Let us look at the dual geometry when $\nu^{\text{strict}}$ is not strictly dual feasible. The non-strictly dual feasible point $\nu^{\text{strict}}$ is on the boundary of the dual feasible set, since it satisfies $a_{i}^{T}\nu^{\text{strict}}=0$ for some $a_{i}$. When $\nu^{\prime}$ is not dual feasible, the only dual feasible point on the line segment between $\nu^{\prime}$ and $\nu^{\text{strict}}$ is $\nu^{\text{strict}}$, so the dual line search returns $\hat{\nu}=\nu^{\text{strict}}$. Thus, when $\nu^{\text{strict}}$ is not strictly dual feasible, the dual line search will return one of two values: when $\nu^{\prime}$ is not feasible, the line search returns $\nu^{\text{strict}}$; when $\nu^{\prime}$ is feasible, the line search returns $\nu^{\prime}$.

Figure 3 illustrates what would happen if $\nu^{\prime}$ converges to $\nu^{\ast}$ while remaining dual infeasible. The line search is “stuck”, returning $\hat{\nu}=\nu^{\text{strict}}$ for all $\nu^{\prime}$. So as $\nu^{\prime}\to\nu^{\ast}$, $\hat{\nu}=\nu^{\text{strict}}$ is “stuck” and does not converge to $\nu^{\ast}$. Picking $\nu^{\text{strict}}$ to be strictly dual feasible “unsticks” the dual line search, allowing the returned value $\hat{\nu}$ to converge to $\nu^{\ast}$ as $\nu^{\prime}\to\nu^{\ast}$.

4.2 Convergence of Dual Sequence Given Primal Sequence

Using a first-order method to solve (1) typically provides a sequence $x^{k}$ of primal feasible points that converges to a primal optimal point $x^{\ast}$ (which may not be unique). For example, the projected gradient method produces such a sequence (see Theorem 10.24 of [Bec17] for a proof of this). In Subsection 4.1, we discussed a dual line search that allows us to produce a dual feasible point $\hat{\nu}^{k}$ from $\nu^{\prime k}\stackrel{{\scriptstyle\operatorname{def}}}{{=}}\nabla f(Ax^{k})$ and a separate dual feasible point $\nu^{\text{strict}}$.

In Subsection 3 we saw a simple feature elimination subproblem (10) that utilized an accurate, but not necessarily optimal, primal-dual pair. The strength of the subproblem depends on the size of the duality gap $\epsilon$. In other words, as the duality gap shrinks, the lower bound produced by the subproblem increases, possibly eliminating the feature. In order for the duality gap $\epsilon$ to shrink to zero as $x^{k}$ converges, we must also have that the dual line search produces $\hat{\nu}^{k}$ that converges to the dual optimal point $\nu^{\ast}$ (recall that strong duality holds for (1) and (2)). If $\nu^{\text{strict}}$ is strictly dual feasible, we will see that $\hat{\nu}^{k}\to\nu^{\ast}$ as $x^{k}\to x^{\ast}$, thus giving $\epsilon\to 0$ as desired.

First, we give a lemma that states ${\nu^{\ast}=\nabla f(Ax^{\ast})}$ for any primal optimal $x^{\ast}$. This comes somewhat directly from the KKT conditions for the dual problem.

Lemma 4.2**.**

Let $x^{\ast}$ be a primal optimal point of (1). Then ${\nu=\nabla f(Ax^{\ast})}$ is the unique dual optimal point of (2).

Proof.

This is a known result coming from optimality conditions for the dual problem. For example, this is a consequence of Theorem 19.1 of [BC17]. ∎

Now we show that the dual line search produces a sequence of points that converges to the dual optimal point.

Theorem 4.3**.**

Let $x^{k}$ be a sequence of primal feasible points that converge to a primal optimal point $x^{\ast}$, and let $\nu^{\text{strict}}$ be a strictly dual feasible point (i.e., $A^{T}\nu^{\text{strict}}>0$). For each $k$, define the dual feasible point $\hat{\nu}^{k}$ by performing the dual line search (12) using ${\nu^{\prime k}=\nabla f(Ax^{k})}$ and $\nu^{\text{strict}}$. Then the sequence $\hat{\nu}^{k}$ of dual feasible points converges to the unique dual optimal point $\nu^{\ast}$.

Proof Note that the map ${x^{k}\mapsto\nu^{\prime k}=\nabla f(Ax^{k})}$ is continuous by assumption. By continuity and Lemma 4.2, $\nu^{\prime k}\to\nu^{\ast}$ as $x^{k}\to x^{\ast}$. But $\nu^{\prime k}$ is not guaranteed to be dual feasible, hence our use of the dual line search. We seek to show that the dual feasible sequence $\hat{\nu}^{k}\to\nu^{\ast}$ as $x^{k}\to x^{\ast}$. By the triangle inequality,

[TABLE]

We already have that $\|\nu^{\prime k}-\nu^{\ast}\|\to 0$, so to complete the proof it remains to show $\|\hat{\nu}^{k}-\nu^{\prime k}\|\to 0$.

Since $\hat{\nu}^{k}$ is computed using the dual line search using $\nu^{\prime k}$ and $\nu^{\text{strict}}$, we have that

[TABLE]

where $t^{k}$ is the optimal value of $t$ from the dual line search for that particular $\nu^{\prime k}$. Since $\nu^{\prime k}\to\nu^{\ast}$ and $\nu^{\ast}$ is unique, we know that $\|\nu^{\prime k}-\nu^{\text{strict}}\|$ is eventually bounded above by a constant. Therefore, we just need to show $t^{k}\to 0$ to imply that $\|\hat{\nu}^{k}-\nu^{\prime k}\|\to 0$.

Recalling the proof of Lemma 4.1, the function

[TABLE]

is continuous precisely when $a_{i}^{T}\nu^{\text{strict}}>0$ (i.e., when $\nu^{\text{strict}}$ is strictly dual feasible). Since $\nu^{\prime k}$ converges to a dual feasible point, $a_{i}^{T}\nu^{\prime k}$ converges to a nonnegative value for each $a_{i}$. Therefore the limiting value of $t(a_{i}^{T}\nu^{\prime k};a_{i}^{T}\nu^{\text{strict}})$ is [math] for each $a_{i}$, and thus $t^{k}=\max_{i}t(a_{i}^{T}\nu^{\prime k};a_{i}^{T}\nu^{\text{strict}})\to 0$. This then shows that $\|\hat{\nu}^{k}-\nu^{\prime k}\|\to 0$, which, using the triangle inequality above, implies that $\hat{\nu}^{k}\to\nu^{\ast}$ as desired. ∎

5 When Will the Screening Rule Eliminate all Zero Features?

Here we show that the feature elimination subproblem (10) will eliminate all zero features from the solution, under reasonable assumptions and with sufficiently small duality gap $\epsilon$. Coupled with the use of a first-order method and the dual line search of Subsection 4.1, this means that our feature elimination strategy will eventually eliminate all features that can be eliminated.

Theorem 5.1**.**

Assume that strict complementary slackness holds (i.e., $x_{i}^{\ast}=0$ iff $\langle a_{i},\nu^{\ast}\rangle>0$). Define $\mathcal{I}\stackrel{{\scriptstyle\operatorname{def}}}{{=}}{\{i\,:\,\langle a_{i},\nu^{\ast}\rangle>0\}}$ to be the index set of zero features. Let the pair $\hat{x},\hat{\nu}$ produce a duality gap estimate $\epsilon$ such that

[TABLE]

Then the feature elimination problem (10) will eliminate all zero features.

Proof.

Since we have assumed strict complementary slackness, a zero feature $x_{i}^{\ast}=0$ always corresponds to ${\langle a_{i},\nu^{\ast}\rangle>0}$. Since $\nu^{\ast}$ is unique, the set $\mathcal{I}$ uniquely determines the indexes of zero features. If we did not assume strict complementary slackness, then a zero feature $x_{i}^{\ast}=0$ may be associated with ${\langle a_{i},\nu^{\ast}\rangle=0}$, which cannot be eliminated by the subproblem (10). So we see that strict complementary slackness implies that each zero feature corresponds to a strictly satisfied dual inequality. We must show that the subproblem (10) produces a strictly positive lower bound for $\langle a_{i},\nu^{\ast}\rangle$ for every $i\in\mathcal{I}$.

For the subproblem to produce a strictly positive lower bound for the $i$th feature, the search set $N$ must be contained strictly in the interior of the halfspace $\langle a_{i},\nu\rangle\geq 0$. The search set $N$ is a (closed) ball of radius $\sqrt{2L\epsilon}$ centered at $\hat{\nu}$. Recall from Section 3 that the optimal value of (10) is $\langle a_{i},\hat{\nu}\rangle-\sqrt{2L\epsilon}\|a_{i}\|$ for the $i$th feature. Up to a scaling factor, this is the minimum distance between the search set $N$ and the hyperplane $\langle a_{i},\nu\rangle=0$. Therefore the search set $N$ is strictly separated from the hyperplane precisely when the optimal value is strictly positive:

[TABLE]

By assumption on the size of $\epsilon$, this condition is satisfied for each $i\in\mathcal{I}$, so feature elimination will eliminate all zero features. ∎

Theorem 5.1 shows that under strict complementarity and if the duality gap $\epsilon$ is sufficiently small, then feature elimination will eliminate all zero features from the problem. We used knowledge of the the exact dual optimal point only to assist in quantifying how small $\epsilon$ must be in order to imply that feature elimination will work. But even without knowledge of the dual optimal point, we still know that if $\epsilon$ is sufficiently small, then feature elimination will have worked. Indeed, by combining Theorems 4.3 and 5.1, we have the following:

Corollary 5.2**.**

Assume that strict complementary slackness holds (i.e., $x_{i}^{\ast}=0$ iff $\langle a_{i},\nu^{\ast}\rangle>0$). Let $x^{k}$ be a sequence of primal feasible points that converges to $x^{\ast}$ (e.g., from a first-order method) and let $\hat{\nu}^{k}$ be the sequence of dual feasible points produced as in Theorem 4.3. Then the duality gap ${\epsilon=f(x^{k})-g(\hat{\nu}^{k})}\to 0$ as $k\to\infty$ and Theorem 5.1 will eventually apply. This means that the feature elimination subproblem (10) will eventually eliminate all zero features.

This tells us that if we do enough iterations of a first-order method, perform the dual line search, and then do feature elimination, we will eliminate all possible features. But we don’t know how many iterations are sufficient (without knowledge of the dual optimal point, that is). Furthermore, if strict complementarity does not hold, we can only eliminate zero features that correspond to $\langle a_{i},\nu^{\ast}\rangle>0$; a zero feature that corresponds to $\langle a_{i},\nu^{\ast}\rangle=0$ cannot be eliminated. These issues notwithstanding, we can still use feature elimination very effectively in practice, including certifying that underdetermined NNLS problems have unique solutions.

6 Certifying NNLS Solution Uniqueness

Here we consider applying safe feature elimination to the problem of certifying the uniqueness of the primal solution. We consider the case of NNLS, where $f(z)=\frac{1}{2}\|z-b\|^{2}$, which reduces the primal problem (1) to

[TABLE]

We have assumed throughout that the ${m\times n}$ matrix $A$ is full-rank. But we have not yet assumed anything about the shape of $A$, which may be “overdetermined” ($m\geq n$) or “underdetermined” ($m<n$). In the overdetermined case, the primal objective is $\sigma_{\text{min}}(A^{T}A)$-strongly convex, where $\sigma_{\text{min}}(A^{T}A)$ is the minimum singular value of $A^{T}A$. Since $A$ is overdetermined, the Hessian $A^{T}A$ of the NNLS problem is non-singular, so $\sigma_{\text{min}}(A^{T}A)>0$. The NNLS primal problem (14) therefore has a unique optimal point [Bec17]. In the underdetermined case the Hessian $A^{T}A$ is singular and the primal objective is neither strongly nor strictly convex, so there is no such uniqueness guarantee.

To attempt to certify the uniqueness of solutions to underdetermined problems we use our feature elimination strategy to reduce the problem to an overdetermined, full-rank NNLS problem. Suppose we eliminate $r$ features/columns of $A$. This allows us to form the reduced matrix $A_{\text{red}}$ with those $r$ columns removed and with the guarantee that the removed columns are not used by a solution of the original problem. If $r\geq n-m$, so that the reduced matrix $A_{\text{red}}$ is overdetermined, and if $A_{\text{red}}$ is full-rank, then the reduced NNLS problem has a strongly convex objective and has a unique solution. Since our feature elimination strategy is safe, the solution to the original NNLS problem is guaranteed to be the same as the solution to the reduced problem (with appropriate zero padding), meaning that the original NNLS problem has a unique solution. Thus we have a procedure to robustly certify the uniqueness of NNLS problems via our safe feature elimination strategy, which requires an accurate, but not optimal, primal-dual pair.

Note that to certify uniqueness we need not eliminate all zero features, as was the goal of Theorem 5.1. In a sense, certifying uniqueness is an easier problem than eliminating all zero features; indeed, to certify uniqueness we need only eliminate sufficiently many ($r\geq n-m$) features. We therefore do not require (full) strict complementary slackness, as was assumed in Theorem 5.1. There may be some zero features $x_{i}^{\ast}=0$ paired with $\langle a_{i},\nu^{\ast}\rangle=0$ but that does not concern us as long as there are sufficiently many $x_{i}^{\ast}=0$ such that $\langle a_{i},\nu^{\ast}\rangle>0$, enabling us to certify the uniqueness of $x^{\ast}$.

Observe that we can generalize this uniqueness certification procedure to a general objective $f(Ax)$ where $f$ is strictly convex. When $A$ is underdetermined the objective $f(Ax)$ is no longer strictly convex. But if we use SAFE to eliminate sufficiently many features such that the reduced matrix is overdetermined and full-rank, the reduced objective is strictly convex and therefore a minimizer is unique. One can also modify this uniqueness certification technique to work with $\ell_{1}$-regularized problems (like lasso), for instance using the GAP SAFE rules of [FGS15], which are similar to (10).

6.1 A Small NNLS Example

Let us illustrate our procedure with a small example. For NNLS, the dual problem (2) reduces to

[TABLE]

Suppose we have a primal feasible point $\hat{x}$ and dual feasible point $\hat{\nu}$. This allows us to compute the duality gap $\epsilon=f(A\hat{x})-g(\hat{\nu})$. The basic feature elimination subproblem (10) reduces to

[TABLE]

where we have used that the dual objective $g$ is 1-strongly convex (since $f$ has $1$-Lipschitz continuous gradient).

Consider the following matrix with randomly chosen entries

[TABLE]

and right-hand side (RHS) $b=\begin{bmatrix}-1&2&1\end{bmatrix}^{T}$. We note that there is nothing special about these entries; the entries are the first few digits of the golden ratio, the base of the natural logarithm, and pi, with some negative signs added. Projected gradient descent (PGD) for NNLS produces the iteration

[TABLE]

where we pick step size $t=1/\|A\|^{2}$. Starting with $x=0$ and iterating $250$ times yields the primal feasible point $\hat{x}\doteq\begin{bmatrix}0&0&0.9282&0&0.5409\end{bmatrix}^{T}$.111We use $\doteq$ to denote equality up the number of digits shown.

We now find a dual feasible point $\hat{\nu}$ with the dual line search of Section 4. First we find a strictly dual feasible point $\nu^{\text{strict}}$ via the LP (13). Because it makes the numbers more presentable on paper, we opt to rescale the solution $\nu^{\text{strict}}$ so that $\|\nu^{\text{strict}}\|_{1}=1$, which gives $\nu^{\text{strict}}\doteq\begin{bmatrix}0.56&0.34&0.1\end{bmatrix}^{T}$. Then we perform the dual line search (12) with $\nu^{\text{strict}}$ and $\nu^{\prime}=A\hat{x}-b$, giving us the dual feasible point $\hat{\nu}\doteq\begin{bmatrix}0.1387&0.0552&0.0209\end{bmatrix}^{T}$. Together $\hat{x}$ and $\hat{\nu}$ produce the duality gap $\epsilon\doteq 0.0069$.

If we instead found $\hat{\nu}$ via the orthogonal projection subproblem (11), we would find the improved duality gap $\epsilon\doteq 0.0013$. But recall that the orthogonal projection subproblem is closely related to the NNLS dual problem and is computationally expensive to solve. Avoiding this expense is precisely the motivation for the dual line search, and we see for this example that the dual line search is not terribly worse than the orthogonal projection.

We are now ready to solve the feature elimination subproblem (15) once for each of the five columns of $A$. Using the closed-form solution given in Section 3, we find the following lower bounds on $A^{T}\nu^{\ast}$: $\begin{bmatrix}-0.34&0.17&-0.49&0.26&-0.61\end{bmatrix}^{T}$. The lower bounds for $\langle a_{2},\nu^{\ast}\rangle$ and $\langle a_{4},\nu^{\ast}\rangle$ are strictly positive, so we can eliminate them from the problem; the lower bounds for the remaining columns are non-positive, so the test is inconclusive. The reduced matrix is

[TABLE]

which is overdetermined and full-rank. We can therefore certify that the original NNLS problem has a unique solution. In fact, only $206$ iterations of PGD are required to certify uniqueness, though this is unknown a priori.

With the solution certified to be unique, we can bound the distance from $\hat{x}$ to the unique solution $x^{\ast}$ via the strong convexity of the reduced primal problem:

[TABLE]

6.2 An Alternative Method to Certify Uniqueness

Slawski and Hein, as part of their analysis of NNLS problems in [SH*+*13], prove a lemma on the uniqueness of NNLS solutions. Their result is very similar to existing results for $\ell_{1}$-regularized least-squares and related problems [Tib13, ZYC15]. The lemma relies on a strong assumption on the columns of $A$, but provides a simple condition to certify the uniqueness of a solution. We discuss this condition first, state their lemma, and finally discuss how to use their lemma to certify uniqueness in practice.

For an index set $\mathcal{J}\subseteq\{1,\ldots,n\}$, we denote by $A_{\mathcal{J}}$ the submatrix of $A$ formed by taking column $j$ for $j\in\mathcal{J}$. The columns of the matrix $A\in\mathbb{R}^{m\times n}$ are said to be in general linear position (GLP) in $\mathbb{R}^{m}$ if the following condition holds:

[TABLE]

In other words, every subset of $\min\{m,n\}$ columns is linearly-independent. For brevity, we will say “$A$ is in GLP” to mean “the columns of $A$ are in GLP”. It is easy to see that $A$ in GLP implies that $A$ is full-rank, but the converse is not true: GLP is strictly stronger than full-rank. $A$ being in GLP is also related to the spark of $A$, where $\operatorname{spark}(A)$ is defined in [DE03] to be minimum number of columns that form a linearly dependent set. If $m<n$, then $A$ is in GLP iff $\operatorname{spark}(A)=m+1$.

Unlike computing the rank of a matrix, computing the spark of $A$ and determining if $A$ is in GLP may be prohibitively difficult in the worst case. The straightforward computation to determine if $A$ is in GLP requires computing combinatorially many determinants. Indeed, determining if $A$ is in GLP (equivalently, if $\operatorname{spark}(A)=m+1$) is coNP-complete [ACM12, TP14]; computing $\operatorname{spark}(A)$ is NP-hard in general [TP14]. So numerically verifying that $A$ is in GLP is likely intractable except for very small $A$.

But these are worst-case results, when we know nothing about the matrix $A$; there are matrices that are known to be in GLP or have known spark. For example, if the entries of $A\in\mathbb{R}^{m\times n}$ are drawn i.i.d. from an absolutely continuous distribution, then $A$ is in GLP with probability one [Tib13]. Though it is complex, another example is ${\textstyle A=\begin{bmatrix}I_{n}&F_{n}\end{bmatrix}}$ where $I_{n}$ is the $n\times n$ identity matrix and $F_{n}$ is the $n\times n$ discrete Fourier transform matrix. When $n$ is a perfect square, the spark is known to be exactly $2\sqrt{n}$, and hence it is not in GLP for $n>1$ [DE03].

There are also lower bounds for $\operatorname{spark}(A)$ [DE03, Tro04]. One such bound is $\operatorname{spark}(A)>1/\mu(A)$, where

[TABLE]

is called the coherence parameter of $A$. For the $A\in\mathbb{R}^{1681\times 2822}$ from our microscopy example in Subsection 7.2, we have $\mu(A)\approx 0.99$ which gives the uninformative bound $\operatorname{spark}(A)\geq 2$.

Assuming we know that $A$ is in GLP (e.g., if $A$ is drawn with entries from a continuous distribution like the standard normal distribution), the following lemma from [SH*+*13] gives a simple condition implying the uniqueness of the NNLS solution.

Lemma 6.1** (Lemma 5 from [SH*+*13]).**

Let the columns of $A\in\mathbb{R}^{m\times n}$, $m<n$, be in GLP. If the NNLS optimal value is strictly positive,

[TABLE]

then the NNLS problem has a unique solution. Furthermore there are at most $m-1$ non-zero values in the solution.

For underdetermined NNLS problems with $A$ in GLP, we can certify uniqueness simply by certifying $p^{\ast}>0$. Assuming we know that $A$ is in GLP, this is simple to check and certify in practice, including when using a first-order method that produces only primal points. We can produce a dual feasible point $\hat{\nu}$ from a primal feasible point $\hat{x}$ using the dual line search from Section 4. If we have that $g(\hat{\nu})>0$, then $p^{\ast}>0$ by weak duality and the solution is certified to be unique. But of course if $A$ is not known to be in GLP we cannot invoke Lemma 6.1.

The small example problem in the previous subsection has $A$ in GLP, which can be checked directly since it is so small. It takes $286$ iterations of PGD to certify that $p^{\ast}>0$, which is slightly more than the $206$ iterations needed for SAFE to certify uniqueness.

To certify uniqueness using safe feature elimination, there must be at least ${n-m}$ zero features with strict complementarity. If this condition does not hold, safe feature elimination will never certify uniqueness. If $A$ is in GLP and ${p^{\ast}>0}$, then one can certify uniqueness using Lemma 6.1. But notice from Lemma 6.1 that under such conditions, there are at least $n-m+1$ zero features in the solution; so safe feature elimination will also certify uniqueness, provided the solution exhibits enough strict complementarity. Even if feature elimination fails to certify uniqueness, it still provides certificates that features are not present in the solution. This is a positive result, whereas Lemma 6.1 provides no additional benefit when it fails to certify uniqueness.

7 Certifying NNLS Solution Uniqueness - Examples

7.1 Synthetic Data Examples

Let us now see how safe feature elimination performs on a larger synthetic example. We construct a random NNLS problem by drawing a random $50\times 100$ matrix $A$ and $50\times 1$ RHS $b$ each with entries drawn i.i.d. from the standard normal distribution $\mathcal{N}(0,1)$. Such an NNLS problem may not necessarily have a unique solution. To check if it does we find a high-accuracy solution using matlab’s lsqnonneg, which implements an active set method from [LH95]. Using the numerically optimal solution and noting that $A$ is in GLP with probability one, we check if $p^{\ast}>0$ to certify that the solution is unique. If it not certified to be unique, we draw another random NNLS problem until we have a problem with a unique solution.

Although it is outside the scope of the present work, it is interesting to note that the uniqueness of the solution to random NNLS problems of this form appears to depend sharply on the shape of $A$. If $m>n/2$, the solution appears to be unique with high probability for large $m,n$; if $m<n/2$, the solution appears to be non-unique with high probability for large $m,n$. “Phase transitions” of a similar form are analyzed in [ALMT14] and it seems quite possible to extend their results to random NNLS problems of the form used here.

In preparation for the dual line search, we find a strictly dual feasible point $\nu^{\text{strict}}$ by solving the LP (13) once and precomputing $A^{T}\nu^{\text{strict}}$. We run $7500$ iterations of projected gradient descent starting with $x=0$, with each iteration giving a primal feasible point $\hat{x}$. At each iteration we use the dual line search (12) to construct a dual feasible point $\hat{\nu}$. Performing the dual line search requires computing $A^{T}\hat{\nu}$ plus $\mathcal{O}(n)$ work, which is on the order of a single gradient evaluation. Using $\hat{x}$ and $\hat{\nu}$, we find the duality gap and use the feature elimination subproblem (15) to eliminate features.

Figure 4 shows the result of using feature elimination to certify the uniqueness of the random NNLS problem. After about $1500$ iterations, the duality gap is small enough that feature elimination has started to eliminate features. Just after $3000$ iterations, sufficiently many features are eliminated to certify uniqueness ($A$ in GLP implies that the reduced matrix $A_{\text{red}}$ is full-rank, but we can also verify this numerically). In accordance with Corollary 5.2 we find that SAFE eventually eliminates all zero features.

For comparison, we also use SAFE with the orthogonal projection (11) instead of the dual line search. Solving the orthogonal projection subproblem at each step of projected gradient descent is tractable for this small problem, but is impractical for larger problems. We see that the dual line search, which scales well to large problems, performs only a bit worse than the orthogonal projection.

Figure 4 also shows that Lemma 6.1 certifies uniqueness for quite a large duality gap. For the problem used for Figure 4, $\nu^{\text{strict}}$ is sufficient to certify $p^{\ast}>0$, which certifies uniqueness before even the first iteration. While impressive, this is not “for free” since we still solve the LP (13) to find $\nu^{\text{strict}}$. For other instances more PGD iterations are required, but it is typical for these problems that $p^{\ast}>0$ is certified before SAFE has eliminated sufficiently many features. After SAFE has certified uniqueness, the reduced primal problem is strongly convex, which allows us to bound the distance from the primal iterate $\hat{x}$ to the true solution $x^{\ast}$ à la (17). Though not always possible (unlike a bound on the duality gap, which we can always find), this provides quite strong information relating the iterate $\hat{x}$ to the optimal point $x^{\ast}$.

7.2 Microscopy Uniqueness Example

We now consider a challenging set of NNLS problems arising from a new technique in super-resolution fluorescence microscopy. In this instance, the image formation process involves solving $40000$ NNLS problems each using the same matrix $A$ but different RHS $b$. The solutions of the NNLS problems are then assembled into the final image. We refer the reader to [YBF*+*18] for details on the microscope and NNLS problem setup.

It is natural to ask if the final, super-resolved image is uniquely determined given the data. If each of the $40000$ NNLS problems has a unique solution, then the final image is unique. We answer that question in the affirmative by using feature elimination to certify the uniqueness of each NNLS problem. In fact these microscopy problems motivated our development of feature elimination, in particular developing them to work with just a primal feasible point $\hat{x}$ coming from a first-order method. Note that we do not know a priori, and cannot verify numerically, if $A$ is in GLP, so we cannot use Lemma 6.1.

The matrix $A$ is $1681\times 2822$ and has $2$-norm condition number $\kappa_{2}(A)=2.4\times 10^{20}$ (computed using dgesvj compiled to use quadruple precision [ABB*+*99]). Even though these NNLS problems are extremely ill-conditioned, first-order methods are well-suited to solve them since each NNLS problem uses the same $A$. This structure allows us to combine gradient computations for many RHS into matrix-matrix products with $A$ and $A^{T}$, instead of repeated matrix-vector products with $A$ and $A^{T}$. High-performance matrix-matrix product implementations take advantage of modern hierarchical memory computers to achieve higher performance than repeated matrix-vector products [GVL98, GG08]. The result is an order of magnitude speedup in the gradient evaluation time (throughput, specifically). For further improved speed, we implement an optimal/accelerated first-order method from [AT06] instead of using PGD; we will refer to this method as AT. We include many of the implementation tricks from TFOCS [BCG11], including an adaptive step size selection method. Our implementation uses a GPU for the matrix-matrix products and array operations in the iteration, leading to further improved runtime.

Like the example in Subsection 7.1, we will iterate AT for some number of iterations, then stop and perform feature elimination. We perform the dual line search using $\nu^{\text{strict}}=\max\{0,\nu^{\prime}\}$, as mentioned in Subsection 4.1.1. Using the closed-form solution to the dual line search requires computing $A^{T}\hat{\nu}$ and $A^{T}\nu^{\text{strict}}$, which is on the order of the cost of a gradient evaluation. That is to say the dual line search is not terribly expensive, though we do not generally want to do it after each iteration of AT.

Table 1 shows the results of using the strong concavity subproblem (10) (we also show results for subproblem (19), which we discuss shortly). We show the number of iterations of AT, the total number of solutions certified to be unique, and the number of features eliminated across all NNLS problems. There are $40000$ NNLS problems, each with $2822$ features, giving approximately $113$ million features total. As the accuracy of the primal feasible points $\hat{x}$ increases, the duality gap closes and more features are eliminated from the problem. But even at 500K iterations, not all problems are certified to have a unique solution.

We know from Corollary 5.2 that we could simply perform more iterations of AT to shrink the duality gap. But instead let us construct a stronger feature elimination subproblem, allowing us to expend a little more work in solving the new subproblem to avoid computing more iterations of AT. We do this by introducing another constraint on (10) to shrink the search set $N$ while still ensuring that $\nu^{\ast}\in N$. One of many ways to do this is by adding the dual feasibility constraint: $A^{T}\nu\geq 0$. This has the possibility to shrink $N$, thereby increasing the feature elimination lower bound, while guaranteeing that $\nu^{\ast}\in N$. Thus this leads to a stronger but still safe subproblem.

But the resulting feature elimination subproblem is too difficult to solve for our purposes. We relax the subproblem by enforcing dual feasibility for only a single column at a time with $a_{j}^{T}\nu\geq 0$. Since we can pick any $j\in\{1,\ldots,n\}$, we solve the subproblem for each column and take the largest lower bound:

[TABLE]

Each inner problem is a “dome subproblem”, since the feasible set is the intersection of a ball and a halfspace. A closed-form solution to the dome subproblem exists (see [GVR12] for instance), allowing us to compute the optimal value cheaply and accurately. The dome subproblem optimal value uses the inner products $\langle a_{j},\hat{\nu}\rangle$ and $\langle a_{i},a_{j}\rangle$. With all required inner products computed, evaluation of the optimal value requires $\mathcal{O}(1)$ work.

Observe that we can compute the all the required inner products for all dome subproblems as $A^{T}\hat{\nu}$ and $A^{T}A$. Since $A$ is fixed we precompute $A^{T}A$ and discard this one-time cost. This brings the cost of computing the optimal value of (19) to about the cost of a gradient evaluation plus $\mathcal{O}(n)$ work for $n$ evaluations of the dome subproblem optimal value.

Table 1 shows the analogous results when using the strong concavity and partial dual feasibility subproblem (19). We see a marked improvement in the number of solutions certified to be unique, though we still fall a bit short of certifying uniqueness for all 40000 problems. This appears to be due to a few particularly slow-to-converge problems where the accuracy of $\hat{x}$ is still quite low. We fix this by computing a high-accuracy solution for the remaining 1996 problems using matlab’s lsqnonneg (note that we do not use lsqnonneg for the 5538 remaining problems when using subproblem (10)). This results in a sufficiently accurate $\hat{x}$ and we certify the remaining problems as having unique solutions. The final image is constructed by assembling the individual NNLS solutions, so by certifying that all NNLS solutions are unique we also guarantee that the final image is uniquely determined from the data.

8 Conclusions

We have developed a safe feature elimination strategy for non-negativity constrained convex optimization problems which uses an accurate, but non-optimal, primal-dual feasible pair. We show that under reasonable conditions, a sufficiently accurate primal-dual pair will eliminate all zero coordinates from the problem. To enable our methods to work with optimization algorithms that produce only primal points we also developed a dual line search to construct an accurate dual feasible point from an accurate primal feasible point. This allows us to use a first-order method to solve the primal, use the dual line search to cheaply construct a dual feasible point, and then use SAFE to eliminate features. We demonstrate the use of SAFE to robustly certify the uniqueness of a non-negative least-squares solution in a small synthetic data example and also for a large-scale, extremely ill-conditioned problem set arising from a microscopy application. Once an NNLS solution has been certified unique, safe feature elimination also provides a bound on the distance to the unique optimal point. Possible future directions of this work include strengthening the feature elimination subproblems and dual line search, and extending the uniqueness certification technique to $1$-norm regularized problems like lasso.

9 Acknowledgements

Stephen Becker acknowledges the donation of a Tesla K40c GPU from NVIDIA.