ResQPASS: an algorithm for bounded variable linear least squares with asymptotic Krylov convergence

TL;DR

ResQPASS is a novel algorithm for large-scale bounded linear least squares problems that combines active-set methods with Krylov subspace techniques, offering efficient convergence and practical implementation strategies.

Contribution

The paper introduces ResQPASS, a new method that integrates active-set quadratic programming with Krylov subspace methods for improved large-scale constrained least squares solving.

Findings

Converges like CG/LSQR when constraints are inactive

Efficient implementation with QR and Cholesky updates

Links convergence to asymptotic Krylov subspace

Abstract

We present the Residual Quadratic Programming Active-Set Subspace (ResQPASS) method that solves large-scale linear least-squares problems with bound constraints on the variables. The problem is solved by creating a series of small problems of increasing size by projecting onto the basis of residuals. Each projected problem is solved by the active-set method for convex quadratic programming, warm-started with a working set and solution from the previous problem. The method coincides with conjugate gradients (CG) or, equivalently, LSQR when none of the constraints is active. When only a few constraints are active the method converges, after a few initial iterations, like CG and LSQR. An analysis links the convergence to an asymptotic Krylov subspace. We also present an efficient implementation where QR factorizations of the projected problems are updated over the inner iterations and Cholesky or Gram-Schmidt over the outer iterations.

Equations

$$\underset{l\le x\le u}{\min }\frac{1}{2}\Vert Ax-b{\Vert }_2^2,$$

$$A^T(Ax-b)-\lambda +\mu$$

$${\lambda }_i(x_i-{\ell }_i)$$

$${\mu }_i(u_i-x_i)$$

$${\ell }_i\le x_i$$

$$(\lambda ,\mu )$$

$${\mathcal{K}}_k(A^{T}A,A^{T}b)=\text{span}\left\{A^{T}b,(A^{T}A)A^{T}b,\dots ,(A^{T}A{)}^{k-1}{A}^{T}b\right\}.$$

$$x_k=\underset{x\in \mathrm{span}\{r_0,\dots ,r_{k-1}\}}{\mathrm{argmin}}$$

$$l\le x\le u$$

$$r_k:=A^T(A{x}_k-b)-{\lambda }_k+{\mu }_k$$

$$\underset{y_k\in {\mathbb{R}}^k}{\min }$$

$${\ell }_i\le [V_k{y}_k{]}_i\le u_{i}i\in \{1,\dots ,m\}.$$

$$V_k^T{A}^T(A{V}_k{y}_k-b)-V_k^T{\lambda }_k+V_k^T{\mu }_k$$

$$({\lambda }_k{)}_i([V_k{y}_k{]}_i-{\ell }_i)$$

$$({\mu }_k{)}_i(u_i-[V_k{y}_k{]}_i)$$

$${\ell }_i\le [V_k{y}_k{]}_i$$

$$({\lambda }_k,{\mu }_k)$$

$$V_k^T\left(A^T(A{V}_k{y}_k-b)-{\lambda }_k+{\mu }_k\right)=V_k^T{r}_k=0.$$

$${\mathcal{V}}_{k+1}:=\mathrm{span}\left\{r_0,r_1,\dots ,r_k\right\}\cong {\mathcal{K}}_{k+1}(A^{T}A,r_0):=\mathrm{span}\left\{r_0,A^{T}A{r}_0,\dots ,{\left(A^{T}A\right)}^k{r}_0\right\}.$$

$$r_k=A^T(A{V}_k{y}_k-b).$$

$$r_k$$

$$=A^{T}A\left(y_0{r}_0+\cdots +y_{k-1}{r}_{k-1}\right)+r_0$$

$$=r_0+y_0{A}^{T}A{r}_0+\cdots +y_{k-1}{A}^{T}A{r}_{k-1}.$$

$$-\frac{1}{2}\Vert x_i^{\ast }\Vert -0.01\le x_i\le \frac{1}{2}\Vert x_i^{\ast }\Vert +0.01i\in \{1,\dots ,m_{\max }\}$$

$$\text{span}\left\{-{\lambda }_1+{\mu }_1,-{\lambda }_2+{\mu }_2,\dots ,-{\lambda }_k+{\mu }_k\right\},$$

$$\text{span}\left\{r_1-A^T(A{V}_1{y}_1-b),r_2-A^T(A{V}_2{y}_2-b),\dots ,r_k-A^T(A{V}_k{y}_k-b)\right\}.$$

$$-{\lambda }_k+{\mu }_k=r_k-A^{T}A{V}_k{y}_k=\sum _{l=1}^{k-1}(r_l-A^T(A{V}_l{y}_l-b))a_l\in \text{span}[V_k,A^{T}A{V}_{k-1}],$$

$$r_k$$

$$=A^{T}A{V}_k{y}_k+V_{k_0}{\alpha }^{(k)}+A^{T}A{V}_{k_0-1}{\beta }^{(k)}$$

$$=A^{T}A{V}_k(y_k+({\beta }^{(k)},0{)}^T)+V_{k_0}{\alpha }^{(k)}\forall k\ge k_0.$$

$$r_k=A^T(A{V}_k{y}_k-b)-{\lambda }_k+{\mu }_k=\sum _{m=0}^{k-k_0+1}(A^{T}A{)}^m{V}_{k_0}{\gamma }_m^{(k)}\forall k\ge k_0.$$

$$r_{k_0}$$

$$=A^{T}A{V}_{k_0}\left(y_k+({\beta }^{(k_0)}0{)}^T\right)+V_{k_0-1}{\alpha }^{(k_0)}=\sum _{m=0}^1(A^{T}A{)}^m{V}_{k_0}{\gamma }_m^{(k_0)},$$

$$r_{(k_0+1)}=(A^{T}A{)}^2{V}_{k_0}{\gamma }_2^{(k_0+1)}+(A^{T}A)V_{k_0}{\gamma }_1^{(k_0+1)}+V_{k_0}{\gamma }_0^{(k_0+1)}.$$

$$r_k$$

$$=\sum _{m=0}^{k-k_0+1}(A^{T}A{)}^m{V}_{k_0}{\gamma }_m^{(k)}.$$

$$W_{k_0,l}:=\text{span}\left\{v_{k_0},v_{k_0+l-1},\dots ,v_k\right\}=\text{span}\left\{r_{k_0-1},\dots ,r_{k_0+l-2}\right\}$$

$$B_{k_0,l}{V}_{k_0}\eta$$

Full text

{MSCcodes}

49N45, 65K05, 65F10, 90C20

1 Introduction

Inverse problems reconstruct an unknown object from measurements based on a mathematical model that describes the data acquisition. It inverts the measured signal into information about an unknown object. A well known inverse problem is computed tomography (CT), where an unknown object is reconstructed from a projection sinogram. In algebraic reconstruction the problem is formulated as a linear least squares problem, $\min_{x\in\mathbb{R}^{n}}\|Ax-b\|^{2}_{2}$, that minimizes the mismatch between the data and model. Here, the matrix $A\in\mathbb{R}^{m\times n}$ models the propagation of the X-rays through the object [8], $x\in\mathbb{R}^{n}$ describes the unknown pixel values and $b\in\mathbb{R}^{m}$ is the vector with noisy measurements.

However, since $A$ is ill-conditioned, e.g. a vector with oscillating $+1$,$-1$ pixel values is in the near null space, the straightforward least-squares solution, $x=(A^{T}A)^{-1}A^{T}b$ will be contaminated due to the noise in the measurements [6]. To avoid overfitting to the noise, regularization techniques are used. For example, in Tikhonov regularization [2] the reconstruction problem is replaced by $\min_{x\in\mathbb{R}^{n}}\|Ax-b\|_{2}^{2}+\lambda\|Lx\|_{2}$ where a regularization term with a regularization operator $L$ is added. The regularization parameter $\lambda>0$ can be tuned such that the Morozov discrepancy principle is satisfied [12]. In LASSO [17] the problem becomes $\min_{x\in\mathbb{R}^{n}}\|Ax-b\|_{2}^{2}+\lambda\|Lx\|_{1}$. The regularization term now uses a $\ell_{1}$-norm that promotes sparsity in the solution. Elastic net [20] combines the $\ell_{1}$-norm and the $\ell_{2}$-norm in the regularization.

The focus of this paper is on an alternative regularization technique where bounds on the variables limit the sensitivity of the solution to measurement noise. We solve

[TABLE]

where the lower and bounds are vectors $l,u\in\mathbb{R}^{n}$. This problem is known as the bounded-variable least squares (BVLS) that can, for example, be solved by the Stark-Parker algorithm [16]. When only lower bounds are used, it is a non-negative least squares (NNLS) problem [9]. Note that also the Lasso problem can be formulated as a non-negative QP problem.

This results in a constrained quadratic programming (QP) problem. We will call the objective $f(x)\coloneqq\frac{1}{2}\|Ax-b\|_{2}^{2}$. In tomography, for example, the bounds limit the pixel values of the reconstructed object.

The Karush-Kuhn-Tucker (KKT) optimality conditions [14] for Eq. 1 are

[TABLE]

where $\lambda\in\mathbb{R}^{n}$ and $\mu\in\mathbb{R}^{n}$ are vectors with Lagrange multipliers associated, respectively, for the lower and upper bounds.

The state-of-the-art method to solve a general large-scale QP or linear programming (LP) problem are interior-point methods [14, 3], where each iteration a weighted sparse-matrix needs to be factorized by a Cholesky factorization. For large-scale LP problems column generation [11] is a technique where the block structure of the matrix is exploited by writing the solution as a convex combination of solutions of subproblems. This results in an iterative method that alternates between solving a master problem, which finds the coefficients of the convex combination, and a new subproblem that expands the vectors for the convex combination.

In this paper we use a linear combination of residuals as an approximation to the solution. Each iteration we solve a projected problem, similar to the master problem in column generation. The use of residuals as a basis for the expansion is common in Krylov methods.

The use of Krylov methods is widespread in scientific computing. They are attractive for applications with sparse matrices. Since a sparse-matrix vector product, $Av$, is computationally inexpensive and paralellisable, it is cheap to generate a Krylov subspace. The convergence of Krylov methods is well understood and various preconditioning techniques such as multigrid or incomplete factorisation techniques can accelerate the convergence. For a review of Krylov methods and various preconditioning techniques we refer to [10, 15].

For inverse problems that typically involve a non-square matrix Golub-Kahan bi-diagonalization is used [1]. It searches a solution in the subspace formed by the matrix powers,

[TABLE]

So far as we know Krylov methods have not been extended to include lower and upper bounds on the variables. There is a literature on enforcing non-negative constraints in [13].

In this paper we make the following contributions. We propose a new subspace method that uses the residuals from the KKT condition, (2). Each iteration we solve a small projected QP problem with dense matrices using the active-set QP algorithm (QPAS). When only a few bounds are active in the solution we observe superlinear convergence behavior. We explain this superlinear behavior by making the link to Krylov convergence theory. From a certain point on, the residuals can be written as polynomials of the normal matrix applied to a vector and the convergence is determined by the Ritz values of the projected matrix.

In addition, we contribute an efficient implementation that uses warm starting, as the basis is expanded, and updates the factorization of the matrices that need to be solved each iteration. Additional recurrence relations gives further improvements.

The analysis of the propagation of the rounding errors and the final attainable accuracy is not part of this study and will be the subject of a future paper.

The outline of the paper is as follows. In Section 2 we derive the algorithm and prove some properties of the basis. Section 3 discusses several ways to speed up the algorithm and Section 4 discusses some synthetic numerical experiments representative for different application areas. We conclude in Section 5.

2 The Residual QPAS Subspace method and its convergence

In this section we introduce the methodology and analyse the convergence.

2.1 Residual QPAS Subspace Method

We propose to solve Eq. 1 by iteratively solving a projected version of the problem. We write the solution $x_{k}=V_{k}y_{k}$ with a zero initial guess.

Definition 2.1**.**

The residual QPAS subspace iteration for a $A\in\mathbb{R}^{m\times n}$, $b\in\mathbb{R}^{m}$ and $l,u\in\mathbb{R}^{n}$, lower and upper bounds such that $l\leq 0\leq u$, generates a series of guesses $\{x_{k}\}_{k\in\mathbb{N}}$ that are solution of

[TABLE]

where

[TABLE]

*and $\lambda_{k},\mu_{k}\in\mathbb{R}^{n}$ are the Lagrange multipliers associated with the lower and upper bound. The feasible initial guess is $x_{0}=0$ and $r_{0}:=-A^{T}b$. *

Remark 2.2**.**

*The condition $l\leq 0\leq u$, to ensure feasibility of $x_{0}$, does not imply any restrictions on the problems that can be solved. The problem with arbitrary $l<u$, can always be shifted, such that $l\leq 0\leq u$ holds. *

A high level implementation of this algorithm is given in algorithm 1 and a more detailed implementation in algorithm 3.

If $V_{k}\in\mathbb{R}^{n\times k}$ with $r_{k}$ is defined in (5), is a basis for $\text{span}\{r_{0},\ldots,r_{k-1}\}$, the optimization problem at iteration $k$ is

[TABLE]

The coefficients $y_{k}\in\mathbb{R}^{k}$ are the projections of the unknown solution on $V_{k}$ and $k$ is the dimension of the subspace.

The corresponding KKT conditions for (6) are now:

[TABLE]

We call the active set the set $\mathcal{A}_{k}\subset\{1,\ldots,n\}$ where the bound constraint, (7d), become equalities.

As is done in Krylov-type methods (see [10]) we will use residuals $r_{k}$ to expand our basis. For every iteration $k$, the basis of the subspace becomes $V_{k}=\begin{bmatrix}r_{0}/\|r_{0}\|&\cdots&r_{k-1}/\|r_{k-1}\|\end{bmatrix}$. This choice of residual basis is natural for Krylov methods such as conjugate gradients (CG). We will prove their pairwise orthogonality in lemma 2.3, how this is a generalization of Krylov in lemma 2.7 and how there is asymptotic Krylov convergence in lemma 2.19.

Lemma 2.3** (Orthogonality of $r_{k}$ for optimal solution of projected system).**

*If $y^{*}_{k}\in\mathbb{R}^{k}$ is an optimal solution of Eq. 6 then the $\{r_{k}\}_{k\in\{1,2\ldots\}}$ as defined in Eq. 5 are pairwise orthogonal. *

Proof 2.4**.**

If $y^{*}_{k}$ is optimal, the KKT conditions (7) hold. The stationarity conditions Eq. 7a can be rewritten as follows:

[TABLE]

*Since $V_{k}$ is a basis for $\text{span}[r_{0},r_{1},\ldots r_{k-1}]$, this proofs the orthogonality of the $r_{k}$’s to $V_{k}$. *

Lemma 2.5** (Orthogonality of $r_{k}$ for feasible solution of projected system).**

*If $y_{k},\lambda_{k},\mu_{k}$ is a feasible solution of 7a then the $\{r_{k}\}_{k\in\{1,2\ldots\}}$ as defined in Eq. 5 are pairwise orthogonal. *

Proof 2.6**.**

*Again, Eq. 7a can be rewritten as as a orthogonality condition. So a feasible solution of the stationarity condition is orthogonal to the previous residuals. *

There exists a natural link between the residuals in Eq. 5 and the Krylov-subspace.

Lemma 2.7** (Generalization of Krylov-subspace).**

In the unconstrained case none of the bounds of Eq. 6 are ever active, we get:

[TABLE]

*With $r_{k}$ defined as in Eq. 5. *

Proof 2.8**.**

If the problem is unconstrained, there are no Lagrange multipliers $\lambda$ and $\mu$. If none of the constraints are ever active, the complementarity condition Eqs. 7b and 7c ensures that $\lambda=0=\mu$. Thus our residuals $r_{k}$ simplify to

[TABLE]

For $k=1$, the statement is trivial. Assume Eq. 9 holds for $k=j-1$. Then

[TABLE]

Example 2.9**.**

*If we now compare our method applied to the unconstrained problem, with the implementation of CG in matlab (pcg) applied to $A^{T}Ax=A^{T}b$, we notice in figure 1 that they converge similarly. This also suggests the possibility of preconditioning. *

Example 2.10**.**

Let us look at an example with a limited number of active constraints, where we have control over the maximal number of active constraints. Consider the least-squares problem $\min\|Ax-b\|_{2}^{2}$ with $A\in\mathbb{R}^{1000\times 600}$ where 4% of the entries are 1 and the rest are 0. Let half of the entries of the exact solution $x^{*}$ be 0 and the others be $\pm 1$. Finally, let $b=Ax$. Let $m_{\max}$ be the maximal number of active constraints, then we add the following constraints to the problem:

[TABLE]

The $0.01$ offset is to ensure that the lower- and upper bounds are never equal. The experiment is performed for $m_{\max}=\{0,1,2,4,8,16,32,64,128\}$. From the results in figure 2 we can conjecture that the method has two steps: discovery of the active set and Krylov convergence. Notice that the discovery phase takes a number of iterations, roughly equal to $m_{\max}$, or the number of active constraints. The number of iterations for the Krylov convergence is always more or less the same as for the unconstrained problem (where there is no discovery phase).

2.2 Convergence theory

Based in these observations we develop a convergence theory.

Lemma 2.11**.**

*Let $\lambda_{k},\mu_{k}\in\mathbb{R}^{n}$ be the solution for the Lagrange multipliers of the projected KKT conditions, (7) for iteration $k$. Let $V_{k}\in\mathbb{R}^{n\times k}$ be the subspace generated by the residuals $r_{k}$, as in Definition 2.1. Then there exists an iteration $0<k_{0}\leq n$ such that for all iterations $k$, where $k_{0}\leq k$ holds, that $-\lambda_{k}+\mu_{k}\in\operatorname{span}([V_{k_{0}},A^{T}AV_{k_{0}-1}])$. *

Proof 2.12**.**

$V_{k}$* can maximally grow to $n$ independent vectors and then ${\text{span}\{V_{n}\}=\mathbb{R}^{n}}$. In that case, the $-\lambda+\mu\in\mathbb{R}^{n}$, solutions of (2), are in $\text{span}(V_{n})$.*

The number of active constraints is determined by the number of non-zeros in $-\lambda+\mu$, solutions of (2) and it is usually much smaller then $n$.

Because we solve (7) with a basis, the number of active bounds in $l\leq V_{k}y_{k}\leq u$ we have at least $k$ equalities. So the number of non-zeros in $-\lambda_{k}+\mu_{k}$, at iteration $k$, is at least $k$.

We can construct the following subspace of $\mathbb{R}^{n}$

[TABLE]

where the first vector has, at least, one non-zero, the second at least two, and so on. At some point, this space spans the exact active set of the problem and also $-\lambda^{*}+\mu^{*}$, the solution of the full problem (2), can be written as a linear combination of the vectors in this subspace. Once these vectors span the space of the active set, this basis does not expand anymore and linear dependence appears.

By construction, see definition of the residual vector (5), these vectors in (15) are equivalent with the following vectors

[TABLE]

As soon as linear dependence appears in (15), the last vector, $-\lambda_{k}+\mu_{k}$, can be written as a linear combination of the previous vectors. This means that there are coefficients $a_{l}\in\mathbb{R}$ such that

[TABLE]

where we use that $-A^{T}b=V_{k_{0}}\|A^{T}b\|e_{1}$, since we start with $x_{0}=0$.

*At the first iteration $k$, where this linear dependence appears, we set $k_{0}:=k-1$. *

The lemma 2.11 means that, for an iteration $k_{0}\leq k$, there is a $\alpha^{(k)}\in\mathbb{R}^{k_{0}}$ and a $\beta^{(k)}\in\mathbb{R}^{k_{0}-1}$ such that we can write the residual as

[TABLE]

Here we use again that $A^{T}b=V_{k_{0}}\|A^{T}b\|e_{1}$, since we start with $x_{0}=0$. Since $\lambda_{k}$ and $\mu_{k}$ can still change with the iteration, the coefficients $\alpha^{(k)}$ and $\beta^{(k)}$ depend on $k$.

Lemma 2.13**.**

If $y_{k}\in\mathbb{R}^{k}$ and $\lambda_{k},\mu_{k}\in\mathbb{R}^{n}$ are the solution of (7) and the iteration $k\geq k_{0}$ is such that 2.11 holds, then there are $\gamma_{i}^{(k)}\in\mathbb{R}^{k_{0}}$ with $i\in\{0,\ldots,k-k_{0}\}$ such that the residual is

[TABLE]

Proof 2.14**.**

Since lemma 2.11 holds, the residual ${r_{k}=A^{T}(AV_{k}y_{k}-b)-\lambda_{k}+\mu_{k}}$ can be rewritten for iteration $k_{0}$, as

[TABLE]

where we can use $\gamma^{k_{0}}_{1}=y_{k_{0}}+(\alpha_{2}^{(k_{0})}0)^{T}$ and $\gamma^{(k_{0})}_{0}=\alpha_{1}^{(k_{0})}$. So the lemma holds for $k=k_{0}$.

The residual $r_{k_{0}}/\|r_{k_{0}}\|$ will be added to the basis and it becomes now $V_{k_{0}+1}=[V_{k}\,\,r_{k_{0}}/\|r_{k_{0}}\|]$.

The solution of the next iteration, $k_{0}+1$, is $V_{k_{0}+1}y_{k_{0}+1}=V_{k_{0}}(y_{k_{0}+1})_{1:k_{0}}+r_{k_{0}}/\|r_{k_{0}}\|(y_{k_{0}+1})_{k_{0}+1}$ and can be rewritten as a linear combination of the vectors $A^{T}AV_{k_{0}}$ and $V_{k_{0}}$. As a result, the residual for iteration $k_{0}+1$ is now

[TABLE]

So for each next iteration $k>k_{0}$ this now becomes.

[TABLE]

We now define the subspace $W_{k_{0},l}$ as the space spanned by $l$ vectors since the linear dependence appeared in the Lagrange multipliers for iterations $k_{0}<k$, including the vector $v_{k_{0}}$. It is

[TABLE]

where $k=k_{0}+l-1$.

We introduce the operator $B_{k_{0},l}$ on this subspace $W_{k_{0},l}$. The action of $B_{k_{0},l}$ on $W_{k_{0},l}$, restricted to $W_{k_{0},l}$ corresponds with the action of $A^{T}A$ on $W_{k_{0},l}$. The operator $B_{k_{0},l}$ is fully determined by its action on the basis vector of $W_{k_{0},l}$. But instead of the basis vectors we use the vectors $(A^{T}A)^{m}V_{k_{0}}\eta$ with $m\in\{1,\ldots,l\}$ with an arbitrary choice of $\eta\in\mathbb{R}^{k_{0}}$.

We then have the equalities, that should hold for any choice of $\eta\in\mathbb{R}^{k_{0}}$

[TABLE]

where, in the last equality, we project again on the subspace $W_{k_{0},l}$. Since these equalities must hold for any $\eta$, we can replace them with matrix equalities. This is similar to the Vorobyev moment problem, see [10].

Lemma 2.15**.**

Let $B_{k_{0},l}=W_{k_{0},l}\,T_{l,l}\,W_{k_{0},l}^{T}$ with $T_{l,l}\in\mathbb{B}^{k_{0}\times k_{0}}$ the operator defined by the equations (24). Let $\{r_{0},r_{1},\ldots,r_{k}\}$ the series of residual vectors from system (7). Then there is a polynomial $P_{l}(t)$ of order $l$ with $P_{l}(0)=1$ such that the projection of $r_{k}$ on $W_{k_{0},l}$ satisfies

[TABLE]

Proof 2.16**.**

Because $r_{k}\perp V_{k}$ it is also holds that also $r_{k}\perp W_{k_{0},l}$, since $W_{k_{0},l}\subset V_{k}$. We can then write, for $k\geq k_{0}$

[TABLE]

*where we use the properties of (24) *

Lemma 2.17**.**

Let $k>k_{0}$ such that lemma 2.11 holds. Let $W_{k_{0},l}=\operatorname{span}\{v_{k_{0}},\ldots,v_{k}\}$ then

[TABLE]

*is tridiagonal. *

Proof 2.18**.**

The subspace $W_{k_{0},l}$ is spanned by the residuals $\text{span}\{r_{k_{0}-1},\ldots,r_{k_{0}+l-2}\}$, this equivalent to $\text{span}\{v_{k_{0}},\ldots,v_{k_{0}+l-1}\}$. Based on lemma 2.13, each of these can be written as

[TABLE]

for $t\in\{0,\ldots,l-2\}$.

When we apply that $A^{T}A$ on these residuals, we have

[TABLE]

Since $\text{span}\left\{r_{k_{0}-1},\ldots,r_{k_{0}+t}\right\}=\text{span}\left\{v_{k_{0}},\ldots,v_{k_{0}+t+1}\right\}$ we can write

[TABLE]

where $H_{l+1,l}\in\mathbb{R}^{(l+1)\times l}$ has a Hessenberg structure, i.e all its elements below the first sub-diagonal are zero. Now, if we project back on $W_{k_{0},l}$ we get

[TABLE]

*Since $A^{T}A$ is symmetric, the left hand side is symmetric, hence the matrix $H_{l,l}$ is also symmetric and becomes tridiagonal. We denote the matrix with $T_{l,l}$. *

Lemma 2.19**.**

*The polynomial $P_{l}(t)$ has zeros in the Ritz values $\theta_{j}^{(l)}$ of $T_{l,l}$. *

Proof 2.20**.**

The proof follows section 5.3 of [18] or theorem 3.4.3 of [10]. The eigenvectors $y_{j}^{(l)}$ of the tridiagonal matrix $T_{l,l}$ span $\mathbb{R}^{l}$. We can then write $e_{1}=\sum_{j}\gamma_{j}y_{j}^{(l)}$.

When we assume that $\gamma_{j}=(y_{j}^{(l)},e_{1})=0$ then

[TABLE]

Since $(T_{l,l}y^{(l)}_{j},e_{1})=(y^{(l)}_{j},T_{l,l}e_{1})=0$ this implies that also $(y_{j}^{(i)},e_{2})=0$. We can then see that $y_{j}^{(l)}\perp\{e_{1},e_{2},\ldots,e_{l}\}$. This leads to a contradiction. Hence, the linear system

[TABLE]

*Since $\gamma_{j}$ are non-zero, the linear system determines the coefficients of $P_{l}(t)$ such that the Ritz values $\theta_{j}^{(l)}$ are zeros $P_{l}(t)$. *

Remark 2.21**.**

*Note that the objective of the minimization problem $f(V_{k}y_{k}^{*})$ decreases monotonically. The objective of our initial guess $f(V_{k+1}[(y_{k}^{*})^{T},0]^{T})$ is equal to the optimal objective of the previous iteration $f(V_{k}y_{k}^{*})$. At worst, QPAS will not make a step towards an improved solution. This does, however, not imply that we remain stuck in this point. Bounds could be removed from the working set and the basis is expanded. *

As soon as the subset spanned by the Lagrange multipliers is invariant, it is clear which variables will be bound by the lower and the upper bounds. But the subspace might not be large enough to accurately solve the full KKT system. If the Lagrange multipliers at iteration $k>k_{0}$ can be written as $-\lambda_{k}+\mu_{k}=V_{k_{0}}d^{(k)}$ for some $d^{(k)}\in\mathbb{R}^{k_{0}}$. The eq. (7a) becomes

[TABLE]

Depending on the sign of $(V_{k_{0}}d^{(k)})_{i}$ for $i\in\mathcal{W}_{k}$ the corresponding $x_{i}=(V_{k}y_{k})_{i}$ is equal to either $l_{i}$ or $u_{i}$.

Let $\mathcal{N}=\{1,\ldots,n\}\setminus\mathcal{A}$ the non-active set. We can write

[TABLE]

when we insert the lower and upper bound we can replace the second part by there respective bounds and we can effectively summarize the terms by a right hand side $\tilde{b}$. In fact, we have eliminated the non-linear complementarity conditions. Using the notation $\tilde{A}={A}_{j\in\{1,2,\ldots,m\},i\in\mathcal{N}}$ we can rewrite the first equation from the KKT condition as

[TABLE]

The method then optimizes the remaining variables. This corresponds to solving a linear system. The Krylov convergence then determines the $y_{k}$ and the corresponding $d_{k}$ that determines the precise value of the Lagrange multipliers.

2.3 Algorithm

Until now we have used the $y_{k}\in\mathbb{R}^{k}$ to denote the solution of the projected problem (6). We will now introduce an inner iteration to solve for the optimal coefficients. We will use $y_{k}^{*}$ to denote the optimal solution of subspace $V_{k}$ and use $y_{k}$ to denote a guess for this solution, in the inner iteration. But this is not necessarily the optimal solution. The working set of an intermediate solution is denoted as $\mathcal{W}_{k}$. In similar way we will use $\mathcal{A}_{k}:=\mathcal{W}_{k}^{*}$ to denote the active set and the optimal working set.

Each outer iteration $k$, the projected optimization problems is solved with the QPAS algorithm [14]. The rationale behind the choice of QPAS is that it has warm-start capabilities that drastically improve the runtime of the algorithm. This is discussed in detail in Section 3.1.

We follow the same strategy as in [19] where a simplex method solves the minimization of $\|Ax-b\|_{1}$ and $\|Ax-b\|_{\infty}$ by creating a sequence of small projected problems.

The algorithm 1 describes the proposed method. The initial basis is the normalized initial residual. We then project the solution to the subspace and solve the small constraint optimization problem with QPAS. We start with $\mathcal{W}_{0}=\emptyset$. For $V_{1}$, the initial guess is $y_{1}=0$ and the working set $\mathcal{W}_{1}=\mathcal{W}_{0}$ is empty. After we have solved the projected problem, we calculate the residual and expand the basis. Because of lemma 2.3 the basis vector is orthogonal to the previous ones.

2.4 Stopping criteria

It is useful to have a stopping criterion to stop the method and return the solution. One possibility is to look at the residual $\|r_{k}\|_{2}$ and stop when it is small. Intuitively, this makes sense because $r_{k}$ is a measure of the distance to the solution that accounts for the bounds using the Lagrangian multipliers. Alternatively, but closely related, we can use the loss of positive definiteness in $V^{T}_{k}A^{T}AV_{k}$. This occurs when a basis vector is added that is linearly dependent on the previous basis. The detection of this loss is relatively straightforward, because a Cholesky factorisation is used and updated (see Section 3.2.2). Positive definiteness is necessary for Cholesky factorisation, and thus our update will fail. In experiments, both methods seem to work fine, as long as the threshold is not taken too small (because of rounding errors). Using both of these methods in tandem worked the best in our experiments, as sometimes semi-convergence would appear. When the semi-convergence sets in the residual would rise again.

3 Implementation

3.1 Warm start

By using the QPAS algorithm for the inner iterations, we can employ warm-starting. Let $\mathcal{W}^{*}_{k}$ be the active set, i.e the optimal working set, corresponding to the solution $y_{k}^{*}$, obtained from QPAS (see algorithm 4). Next, in iteration $k+1$, we expand the basis and warm start with the previous solution. We start QPAS with initial guess $\begin{bmatrix}(y_{k}^{*})^{T}&0\end{bmatrix}^{T}$. Because of this choice, the previous working set $\mathcal{W}_{k}$ is a valid initial working set for outer iteration $k+1$.

This approach yields a significant runtime improvement. In figure 3 we compare the number of QPAS inner iteration for each outer iteration $k$. It shows that once the active set is more or less discovered and superlinear Krylov convergence sets in (see Section 2.2) warm-starting bring significant benefit.

A final advantage of the warm starting is that when the number of QPAS iterations is limited (see also Section 3.3) and a cold start is used, the maximal size of the active set is equal to this limit. Warm starting starts from the previous working set and can then change, add or remove a number of constraints equal to this limit.

3.2 Factorisation update

3.2.1 Solving the linear system in QPAS with QR

In the QPAS algorithm (see algorithm 4), a system of equation is solved each iteration to obtain the Lagrange multiplier $\lambda_{k}$. The linear system is

[TABLE]

Where $G$ is the Hessian of the QP-problem to be solved, $f$ is the linear term and $C$ is the matrix representing the inequality constraints. $C_{\mathcal{W}_{k}}$ is the matrix containing rows of $C$ with indices from the set $\mathcal{W}$. Block elimination leads to the following linear system that needs to be solved, each iteration,

[TABLE]

Every iteration, the working set $\mathcal{W}_{k}$ is expanded, reduced by one index or remains unchanged. This then means that the QR-factorization of $\left(C_{\mathcal{W}_{k}}G^{-1}C^{T}_{\mathcal{W}_{k}}\right)$ can be easily expanded or reduced. There exist routines to efficiently add or remove rows and columns from a QR factorization, [5].

3.2.2 Cholesky

The only linear system that we need to solve is the projected Hessian $G=V_{k}^{T}A^{T}AV_{k}^{T}\in\mathbb{R}^{k\times k}$. Because of the positive definiteness, a natural choice for a factorization is the (lower) Cholesky decomposition. The lower Cholesky decomposition of some Hermitian, positive definite matrix $H$ is $LL^{T}$ with $L$ a lower triangular matrix. As mentioned before, $G$ would only lose positive definiteness once the columns of $V$ become linearly dependent (and we thus have a solution). When this factorization (or its update) fails, the algorithm can be terminated. This allows us to properly terminate in cases where the residual does not shrink below the tolerance (semi-convergence).

When the basis expands $H$ becomes larger. If $H$ is expanded as follows

[TABLE]

with $c\in\mathbb{R}^{k}$ a column vector and $d\in\mathbb{R}$ a scalar. Then the Cholesky factorisation can be updated as follows:

[TABLE]

with

[TABLE]

Remark 3.1**.**

*The forward substitution in line 2 in 2 can be started before the Cholesky update has completed. *

3.3 Limiting the inner iterations

In the early iterations it is often better to expand the basis first and then find the optimal coefficients. Indeed, a small basis cannot fully represent the solution and finding the optimal solution within this subspace is premature. An alternative is to limit the number of inner iterations and calculate the residual when this maximum is reached. In principle, we only require optimality at the final subspace. Limiting the number of inner iterations requires care. If we stop the inner iteration without feasibility then the orthogonality between the residuals is lost.

As can be seen in figure 4 that shows the number of inner iterations required to find the optimal solution for subspace $V_{k}$, there is a sweet spot for what the maximal number of inner iterations should be. If the number of inner iterations is high, we find the optimal solution in each intermediate subspace and we converge in the fewest number of outer iterations. However, each outer iteration takes longer because of the large number of inner iterations. Especially at the halfway point, it appears to be not useful to solve the projected problem exactly. An approximate (feasible) solution is sufficient up to a point. If the number of inner iterations is limited harshly, additional basis vectors are needed in the outer iterations leading to a larger subspace. This then has the effect of taking longer to converge. The sweet spot, in this example, seems to be 5 inner iterations, where we almost get a 33% speed-up.

Note that there is no effect on the orthogonality of the residuals if we stop with a feasible solution. This can be guaranteed if we stop the inner QPAS iterations when the maximum number of inner iteration is reached and $p=0$.

The final algorithm, with all improvements is written down in algorithms 3 and 4

3.4 Stability of QPAS

When a nonzero step $p$ is computed, a stepsize $\alpha_{k}$ is determined. The definition of $\alpha_{k}$ can be found in line 31 in 31, $C_{i}$ is the $i$-th row of matrix $C$. Two matrix vector products are computed in this step: $Cx_{k}$ and $Cp_{k}$. However, the definition of $x_{k+1}=x_{k}+\alpha_{k}p_{k}$ implies the possibility of a recursive relationship for $Cx_{k+1}$:

[TABLE]

where $Cp_{k}$ was already computed to find $\alpha_{k}$. We expect this recursion to reduces runtime. This doesn’t cause instability in our experiment.

In figure 5 we compare the recursive case with 5 and 10 inner iterations to the direct case for a problem from example 2.10 with $A\in\mathbb{R}^{10,000\times 6,000}$ and $m_{\max}=256$. The norm of the recursion error is also shown. The direct method is the slowest, as expected. The case with 10 inner iterations behaves exactly the same as the direct case, only 5% faster. The case with 5 inner iterations is slightly slower because of the delayed convergence caused by the stringent bound on inner iterations. It is however still faster than the direct method.

Note that a similar recurrence is used in CG for the residual. It also avoids a sparse-matrix multiplication and its effect on the rounding error is well studied, see for example Chapter 5 in [10].

4 Numerical experiments

The paper is accademic in nature and mainly explores the possibilities that projections and Krylov-like methods still offer. Therefore the problems are synthetic for now, to show that ResQPASS is promissing. Detailed error analysis and final attainable accuracy are future work along with in depth studies of specialised applications.

4.1 Conjugate Gradients (CG)

As mentioned in example 2.9, the unconstrained ResQPASS applied to $\|Ax-b\|$ converges like CG applied to $A^{T}Ax=A^{T}b$. For this experiment, the entries of $A\in\mathbb{R}^{m\times n}$ and $x_{ex}\in\mathbb{R}^{n}$ are drawn from a normal distribution with mean 0 and standard deviation 1 and all values smaller than 0 and bigger than 0.1 are set to 0. This gives us a sparse $A$ with $3.98\%$ fill. The right hand side is defined as $b=Ax_{ex}$. For our experiment we look at overdetermined and underdetermined systems with an aspect ratio of 10:6. The matrix $A^{T}A$ is close to a full matrix, especially for big $A$. In figure 6 we can see that ResQPASS is still able to exploit the sparsity. If we perform dense matrix-vector multiplication, the performance of ResQPASS and CG is similar for the overdetermined case, but ResQPASS is faster in the underdetermined case. CG can be sped up in this case by not computing $A^{T}A$ explicitly and performing the dense matrix-vector product $(A^{T}A)p$. Instead, we perform two sparse multiplications $A^{T}(Ap)$. In this case ResQPASS performs similarly but with a slight overhead.

4.2 Bounded Variable Least Squares (BVLS)

In this experiment, we use the formulation of example 2.10 with $A\in\mathbb{R}^{10,000\times 6,000}$. First of all, the numerical improvements done to the QPAS method give already a significant advantage compared to MATLAB’s active-set based solver and a naive QPAS implementation. Projection still seems beneficial, as ResQPASS is faster than the improved QPAS. In figure 7 it is also clear that for this problem with a small number of active constraints, ResQPASS outperforms IP. For sufficiently large $m_{\max}$ ResQPASS would however cross over IP (like MATLAB’s QPAS does). This is because the runtime of interior-point (IP) methods is known to be independent of the number of active constraints, whereas for ResQPASS (and other active-set algorithms) it does depend on the number of active constraints [14].

4.3 Nonnegative Matrix Factorisation (NMF)

The $n\times m$ matrix $A$ is given and we look for the unknown $U\in\mathbb{R}^{n\times p}$ and $V\in\mathbb{R}^{p\times m}$, with $X,Y\geq 0$ (elementwise), such that $\|A-XY\|$ is minimal. This can be achieved with the Alternating Least Squares (ALS) algorithm [7]. In the ALS algorithm nonnegative least squares problems need to be solved, for which ResQPASS could be used. As mentioned in remark 2.2 the problem should be translated such that 0 is feasible.

Let the function $\operatorname{vec}(X)$ be the vectorization of the matrix $X$, where columns are stacked. We can rewrite $\|A-UV\|^{2}$ as $\|\operatorname{vec}(A)-(I_{m}\otimes U)\operatorname{vec}(V)\|^{2}$ and $\|\operatorname{vec}(A)-(V^{T}\otimes I_{n})\operatorname{vec}(U)\|^{2}$ to get a classical least squares formulation.

We use the synthetic example of [4] where the matrix $A$ is constructed as follows. $X,Y$ are $n\times p$ and $p\times m$ matrices respectively, with uniformly drawn random entries between 0 and 1. Gaussian noise with mean 0 and standard deviation 0.1 is added to the product $XY$ and all negative values are set to 0 to obtain $A$.

The number of ALS iterations is fixed to 10 (20 minimizations) and the maximum number of inner QPAS iterations for ResQPASS is set to 10. In figure 8 we can see that ResQPASS has a big speed advantage (one or two orders of magnitude), without sacrifying accuracy of the solution.

4.4 Contact problem

In contact problems the freedom of the solution is limited due to spatial limitations. These bounds limit the solution space and we require a subspace that is adapted to these circumstances. We study the convergence of a cartoon model for an inflated balloon with limited space for inflation. In a balloon the elastic forces are balanced by the pressure inside the balloon. This leads to a equilibrium between the curvature and the magnitude of the pressure. So the Laplacian applied to the solution is equal to the pressure vector, which is constant through the domain. This leads to a system of equations $Lx=p$, where $L\in\mathbb{R}^{n\times n}$ is the Laplacian, $x\in\mathbb{R}^{n}$ is the solution that models the surface deformation of the balloon, $p\in\mathbb{R}^{n}$ is the pressure vector.

However, if the solution is limited by lower and upper bounds, this results in a bounded variables QP problem. In this case we minimize the energy with bounds on the solution. This is

[TABLE]

For the cartoon example, we are using $l=0$ and $u=0.1$ and a constant pressure vector. The domain is $[0,1]^{2}$.

The matrix is the 3-point stencil approximation of the Laplacian. The matrix can be written as $A_{\text{2D}}=A_{1D}\otimes I+I\otimes A_{1D}$, where $A_{1D}$ is the 1D finite-difference Laplacian matrix, a tridiagonal matrix with constant diagonals. We use a mesh with 50x50 grid points and a pressure vector with a constant value of 4.

In figure 10 we show the convergence of the ResQPASS method for the contact problem, both for the preconditioned and the unpreconditioned problem. The corresponding solutions after 200 iterations are shown in figure 9.

The preconditioning matrix is based on the incomplete LU factorisation of the Hessian $A^{T}A$. We use a threshold of 0.1 for the ilut factorisation. Instead of using the residual (5) we use a preconditioned residual, solution of

[TABLE]

where $M=LU$ is built out of the incomplete factorisation. When none of the bounds are active this corresponds to the ILU preconditioned CG.

The unpreconditioned iteration does not reach the bounds in the first 200 iterations. It follows the CG convergence. However, for the preconditioned solution we reach the bounds after 15 iterations. From then on the subspace is adapted to the bounds. We see that the solution is adapted to these bounds, see left hand side of figure 9.

5 Discussion and conclusion

In this paper we present the Residual QPAS subspace method (ResQPASS). The method solves box-constrained linear least squa-res problems with a sparse matrix, where the variables can take values between a lower and an upper bound. In this paper we have proposed the subspace method, analysed the convergence and documented how it is efficiently implemented. It uses warm-starting and updates the factorisation.

Similar as in Krylov methods the sequence of residual vectors are pairwise orthogonal. The residual vectors that we use include the current guess for the Lagrange multipliers. The inner problem only needs to be solved for feasibility to obtain an orthogonal basis.

The ResQPASS works fast if the vector of Lagrange multipliers is sparse, i.e. when only a few of the box constraints are active and most Lagrange multipliers are zero. As soon as the active set is discovered, Krylov convergence sets in. We then converge fast to the solution. As in classical Krylov methods, when you restart, it is hard to benefit from superlinear convergence.

Note that it is possible to further accelerate the Krylov convergence with the help of a preconditioner. We then use residuals that are the solution of $Mr_{k}=A^{T}(Ax_{k}-b)-\lambda_{k}+\mu_{k}$ for some non-singular matrix $M$ that is cheap to invert. The residuals are then $M$-orthogonal.

However, there are many inverse problems where many unknowns hit the boundary. For example, a deblurring of a satellite picture against a black background. In that case the Lagrange multipliers are not sparse. However, the deviation from the bounds in $x$ is sparse, hence a dual ResQPASS algorithm might be applicable. Initial experiments confirm this intuition and the dual algorithm exhibits similar superlinear convergence. The analysis of the dual algorithm is the subject of a next paper.

Other applications can be found in contact problems in mechanics, where the subproblems touch each other in only a few places. These points of contact are often unknown upfront.

Other future work is a thorough analysis of the propagation of rounding errors and effects of the loss of orthogonality.

Acknowledgments

We thank Jeffrey Cornelis for fruitful discussions during the initial phase of the research. Bas Symoens acknowledges financial support from the Department of Mathematics, U. Antwerpen.