{'fraction_non_alphanumeric': 0.03611674730763228, 'fraction_numerical': 0.014203215233338536, 'mean_word_length': 4.3048517966550754, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solving large linear least squares problems with linear equality constraints

We consider the problem of efficiently solving large-scale linear least squares problems that have one or more linear constraints that must be satisfied exactly. Whilst some classical approaches are theoretically well founded, they can face difficulties when the matrix of constraints contains dense rows or if an algorithmic transformation used in the solution process results in a modified problem that is much denser than the original one. To address this, we propose modifications and new ideas, with an emphasis on requiring the constraints are satisfied with a small residual. We examine combining the null-space method with our recently developed algorithm for computing a null space basis matrix for a"wide"matrix. We further show that a direct elimination approach enhanced by careful pivoting can be effective in transforming the problem to an unconstrained sparse-dense least squares problem that can be solved with existing direct or iterative methods. We also present a number of solution variants that employ an augmented system formulation, which can be attractive when solving a sequence of related problems. Numerical experiments using problems coming from practical applications are used throughout to demonstrate the effectiveness of the different approaches.


Introduction
Our interest lies in efficient and robust methods for solving large-scale linear least squares problems with linear equality constraints. We assume that A ∈ IR m×n and C ∈ IR p×n , with m > n ≫ p. We further assume that A is large and sparse and C represents a few, possibly dense, linear constraints. Given b ∈ IR m and d ∈ IR p , the least squares problem with equality constraints (the LSE problem) is having full column rank. In the case of non-uniqueness, there is a unique minimumnorm solution. LSE problems arise in a variety of practical applications, including scattered data approximation [12], fitting curves to data [15], surface fitting problems [33], real-time signal processing, and control and communication leading to recursive problems [48], as well as when solving nonlinear least squares problems and least squares problems with inequality constraints. For example, in fitting curves to data, equality constraints may arise from the need to interpolate some data or from a requirement for adjacent fitted curves to match with continuity of the curves and possibly of some derivatives. Motivations for LSE problems together with solution strategies are summarized in the research monographs [7,8,28].
Classical approaches for solving LSE problems derive an equivalent unconstrained linear least squares (LS) problem of lower dimension. There are two standard ways to perform this reduction: the null-space approach [20,28] and the method of direct elimination [10], both of which, with suitable implementation, offer good numerical stability. These methods, termed constraint substitution methods, consider the constraints (1.2) as the primary data and substitute from them into the LS problem (1.1). The former performs a substitution using a null-space basis of C obtained from a QR factorization, while the latter is based on substituting an expression for selected solution components from the constraints into (1.1). This can be done using either a pivoted LU factorization [10] or a QR factorization of C [7]. Other solution methods, which may be regarded as complementary to the constraint substitution approaches, reverse the direction of the substitution, substituting from the LS problem into the constraints. This involves the use of an augmented system and include a Lagrange multiplier formulation [21], updating procedures that force the constraints to be satisfied a posteriori [6,7], and a weighting approach [5,34,44], Solving large-scale LS problems is typically much harder than solving systems of linear algebraic equations, in part because key issues such as ill-conditioning or dense structures within an otherwise sparse problem can vary significantly between different problem classes. Consequently, we do not expect that there will be a single method that is optimal for all LSE problems, and having a range of approaches available that target different problems is important. Our main objective is to revisit classical solution strategies and to propose new ideas and modifications that enable large-scale systems to be solved, with an emphasis first on the possibility that the constraints may be dense, and second on requiring that the constraints be tightly satisfied. In Sections 2 and 3, we consider the null-space method and the direct elimination approach, respectively. We review the methods and show how they can be used for large-scale problems. In Section 4, we present complementary solution approaches within an augmented system framework. This allows us to treat the constraints and the least squares part of the problem using a single extended system of equations or via a global updating scheme. Both direct and iterative methods are discussed.
Much of the published literature related to LSE problems lacks numerical results. For instance, Björck [6] remarks "no attempt has yet been made to implement the (general updating LSE) algorithm", and as far as we are aware, attempts remain absent. We assume this is because implementing the algorithms is far from straightforward. While it is not the intention here to offer a full general comparison of the different approaches, throughout our study we use numerical experiments on problems arising from real applications to highlight key features that may make a method attractive (or unsuitable) for particular problems and to illustrate the effectiveness of the different approaches. Our key findings and recommendations are summarized in Section 5.
We end this introduction by describing our test environment. The test matrices are taken from the SuiteSparse Matrix Collection [14] and comprise a subset of those used by Gould and Scott in their study of numerical methods for solving large-scale LS problems [19]. If necessary, the matrix is transposed to give an overdetermined system. Basic information on our test set is given in Table 1.1. The problems in the top half of the table contain rows that are identified as dense by Algorithm 1 of [42] (with the density parameter set to 0.05). These rows are taken to form the constraint matrix C and all other rows form A. For the other problems, we form A by removing the 20 densest rows of the SuiteSparse matrix; some or all of these rows are used to form C (and the rest are discarded). In our experiments, we prescale the extended matrix A given by (1.3) by normalizing each of its columns. That is, we replace A by A D, where D is the diagonal matrix with entries D ii satisfying D ii = 1/ A e i 2 (e i denotes the i-th unit vector). The entries of A D are at most one in absolute value. The vectors b and d are set to be vectors of 1's (so that b 2 and d 2 are O(1)).
For the substitution approaches described in Sections 2 and 3, we have developed prototype Fortran codes; in Section 4, the augmented system methods are implemented using the SuiteSparseQR package of Davis [13] and Fortran software from the HSL mathematical software library [25]. The prototype codes are not optimised for efficiency and so computational times are not reported. Developing library quality implementations is far from trivial and is outside the scope of the current study, which focuses rather on determining which approaches are sufficiently promising for sophisticated implementations to be considered in the future.
Notation All norms are 2-norms and in the rest of the paper, to simply the notation, . 2 is denoted by . . I is used to denote the identity matrix of appropriate dimension. The entries of any matrix B are (B) i, j and its columns are denoted by b 1 , b 2 , . . .. The null space of B is N (B) and Z is used to denote a matrix whose columns form a basis for the null space (i.e., Z satisfies BZ = 0). Permutation matrices are denoted by P (possibly with a subscript). The normal matrix for (1.1) is H = A T A.

The null-space approach
The null-space approach is a standard technique for solving least squares problems. It is based on constructing a matrix Z ∈ IR n×(n−p) such that its columns form a basis for N (C). Any x ∈ IR n satisfying the constraints can be written in the form where x 1 ∈ IR n is a particular solution of the underdetermined system Cx 1 = d. The minimum norm solution can be obtained from the QR factorization of C, that is, CP C = Q C R C 0 , where the permutation P C ∈ IR n×n represents the pivoting, R C ∈ IR p×p is an upper triangular matrix and Q C ∈ IR p×p is an orthogonal matrix. x 1 is then given by Substituting (2.1) into (1.1) gives the transformed LS problem The method is summarized as Algorithm 1.
Algorithm 1 Null-space method for solving the LSE problem (1.1)-(1. 2) with C of full row rank.
In the 1970s, the null-space method was developed and discussed by a number of authors, including in relation to quadratic programming [20,28,37,43]. These and subsequent contributions formulate the approach via the orthogonal null-space basis obtained, for example, from the QR factorization of C T given by where P ∈ IR n×n is a column permutation of C and Q ∈ IR n×n is an orthogonal matrix. Z is equal to the last (n − p) columns of PQ and consequently is dense. Note that although it is possible to store Q implicitly using, for example, Householder transformations, the memory demands and implied operation counts are generally too high. Our interest is in large LS problems and therefore it may not be practical to solve the (n − p) × (n − p) system in Step 3 if Z is dense. To make the approach feasible for large problems we can exploit our recent work [41] on constructing sparse null-space bases of "wide" matrices such as C that have many more columns than rows and may include some dense rows. Scott and Tůma [41] propose a number of ways to construct sparse Z. In our experiments, we employ Algorithm 3 from Section 3 of [41]. This computes each column of Z independently using its own QR factorization that incorporates a threshold pivoting strategy. When selecting candidate pivots, the aim is to combine using a threshold parameter θ ∈ [0, 1] to ensure stability with locality (that is, column interchanges are between columns that are as close as possible) to maintain sparsity in Z. Small values of θ lead to Z having a narrow bandwidth and (assuming the normal matrix H is sparse) to Z T HZ being sparse but this is potentially at the expense of numerical stability. Our first results are for problems deter3 and gemat1. As discussed in the Introduction, we form the constraint matrix C by taking the p = 2, 5, 10, 20 densest rows of A . The sparse block A is the same for each case. In Figure 2.1, we plot the number of entries nnz(Z T HZ) in Z T HZ and the norm of the constraints residual r c = d − Cx . As expected, nnz(Z T HZ) increases with θ , and this increase grows with p. This is illustrated further by the results in Table 2.1. We can see that, independently of the choice of θ , for some problems (including lp fit2p and sctap1-2r) the constraints are not tightly satisfied. This demonstrates an inherent limitation of the null-space ap- Table 2.1 The density of Z T HZ (that is, nnz(Z T HZ)/(n − p) 2 ) and constraint residual r c for two values of the threshold pivoting parameter θ used in the computation of the null-space basis. ‡ indicates insufficient memory for HSL MA87. proach of [41] that focuses on constructing the columns of Z so as to keep Z T HZ sparse but does not result in Z having orthogonal columns.
The linear system in Step 3 of Algorithm 1 is symmetric positive definite. In the above experiments, we employ the sparse direct solver HSL MA87 [22] (combined with an approximate minimum degree ordering). However, for large problems, the memory demands mean it may not be possible to use a direct method; this is illustrated by problem south31 with θ = 1. If a preconditioned iterative solver is used instead, not only are the solver memory requirements much less but explicitly forming the potentially ill-conditioned normal matrix H can be avoided and because Z only needs to be applied implicitly, the need for sparsity can potentially be relaxed. Currently, finding a good preconditioner for use in this case remains an open problem [30].
If a sequence of LSE problems is to be solved with the same set of constraints but different A, the null-space basis can be reused, substantially reducing the work required. But if the constraints are changed, then Z will also change. In [41], we present a strategy that allows Z to be updated when a row (or block of rows) is added to C.

The method of direct elimination
The second method we look at is direct elimination [28]. The basic idea is to express the dependency of p selected components of the vector x on the remaining n − p components and to substitute this into the LS problem (1.1). Here we propose how to choose the p components so as to retain sparsity in the transformed problem.
Consider the constraints (1.2). The method starts by permuting and splitting the solution components as follows: where P c ∈ IR n×n is a permutation matrix chosen so that C 1 ∈ IR p×p is nonsingular. Let AP c = A 1 A 2 be a conformal partitioning of AP c . Substituting the expression with the transformed matrix Note that if C 1 is irreducible, the transformation combines all the rows of C 2 . If C is composed of dense rows then A T has more dense rows than A. We thus seek to add as few row patterns as possible replicating the (possibly) dense pattern of C within A T . If both A and C are sparse, the substitution leads to a sparse LS problem. We have the following straightforward result.

Lemma 3.1 Let A ∈ IR m×n be sparse. Let m > n > p and assume a conformal column splitting induced by the permutation P c is such that CP
Then the number of dense rows in A T given by (3.3) is at most the number of entries in Occupied.
Proof The result follows directly from the transformation. Assuming the rows of C −1 1 C 2 are dense, the substitution step (3.1) of the direct elimination implies a dense row k in A T only if there is a nonzero in the k-th row of A 1 .
A simple example is given in Figure 3.1. Here we ignore cancellation of nonzeros during arithmetic operations. We see that the pattern of A T satisfies Lemma 3.1. Note that, although in this example C −1 1 C 2 is shown as dense, it need not be fully dense and the number of entries in Occupied represents an upper bound on the number of dense rows in A T . Lemma 3.1 implies that the LSE problem is transformed to a LS problem (3.2) that has some dense rows, which we refer to as a sparse-dense LS problem. Consequently, existing methods for sparse-dense LS problems can be used, including those recently proposed in [38,39,40] (see also the recent direct LS solver HSL MA85 from the HSL library). A straightforward algorithmic implication of the lemma is that the  * * * * * * * * * * * * * * * * * * * * *  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * . The matrix C −1 1 C 2 ∈ IR p×n is given as fully dense.
permuting and splitting of C cannot be separated from considering the sparsity pattern of A because the splitting also determines A 1 and A 2 . Thus we want to permute the columns of C to allow a sufficiently well-conditioned factorization of C 1 while limiting the number of entries in Occupied and hence the number of dense rows in A T . The approach outlined in Algorithm 2 is one way to achieve this. There is an important difference between the pivoting used in Algorithm 3 of [41] (which we used in the previous section) and that of Algorithm 2 below. The former incorporates standard threshold pivoting into a QR factorization algorithm that aims to compute the columns of Z with very few nonzero entries so that the transformed normal matrix Z T HZ is as sparse as possible. Its threshold parameter θ controls the size of the pivots and thus the stability of the factorization. The threshold parameter τ ∈ (0, 1] used in Algorithm 2 also guarantees the pivots in the QR factorization of C are not small but the selection of the candidate pivots is very different. In this case, the choice is based on limiting the fill-in in the transformed matrix A T . Here it is the set of rows held in Occupied that potentially fill-in in A T that plays a crucial role. The use of different notation for the threshold parameters emphasises the difference between the two QR-based approaches. Observe that the pivoting strategy in Algorithm 2 considers C and A simultaneously and will not select a column as the pivot column if this column in A is dense (as it would lead to A T being dense). While we do not discuss the implementation details, we remark that care is needed to ensure efficiency. For example, QR with pivoting for a wide matrix is relatively cheap but it may be necessary to store the squares of the column norms using a heap, which is why we emphasize their role in the algorithm by using the explicit notation w i for these norms.
The effects of increasing the pivoting parameter τ on the number of dense rows in A T are illustrated in Figure 3.2 for problems deter3 and gemat1; results for the full test set are given in Table 3.1. The dense rows of the transformed matrix A T are determined using Algorithm 1 of [42] and to solve the transformed LS problem (3.2) we use the sparse-dense preconditioned iterative approach of [38]. This computes a Cholesky factorization of the normal matrix corresponding to the sparse part of A T and uses it as a preconditioner within a conjugate gradient (CG) method; the CG convergence tolerance that measures relative decrease of the transformed residual Algorithm 2 Assume C = (c 1 , . . . , c n ) ∈ IR p×n (p < n) has full row rank. Compute C 1 ∈ IR p×p and the column permutation P c ∈ IR n×n for the direct elimination method for solving the LSE problem (1.1)-(1.2). P c is determined by a QR factorization with threshold pivoting; τ ∈ (0, 1] is the threshold pivoting parameter. Find Update Occupied ← Occupied ∪ { j |(A) j,i l = 0}. 10: end for 11: Set P c to permute the columns of C with indices in S to obtain C 1 . ||A T T r|| 2 /||r|| 2 is set to 10 −11 . For the problems in the top half of the table for which the rows of C are much denser than those of A (recall Table 1.1), reducing τ leads to only a small reduction in the number ndense of dense rows in A T . However, when the constraints are not dense (the problems in the lower half of the table), ndense can be significantly decreased by choosing τ < 1, although if τ is too small, the matrix C 1 computed by Algorithm 2 can become highly ill-conditioned and A T close to being singular. In our experiments we occasionally observed this for τ < 10 −5 .
By comparing the pairs of problems in the lower half of the table (such as de-ter3 5 and deter3 20) and considering the plots in Figure 3.2, we see that increasing the number p of constraints can lead to a sharp increase in ndense (even if these constraints are relatively sparse), which can result in the transformed problem being hard to solve. The constraints are very well satisfied in all the successful test cases, making this an attractive approach if a good sparse-dense LS solver is available and the number of dense rows in the transformed problem is not too large. Furthermore, it can be used, without modification, if the matrix A contains a (small) number of dense rows. However, for a sequence of problems, if A and/or C changes then, because direct elimination couples the two matrices, the computation must be completely restarted.

Approaches described via augmented systems
We now focus on complementary approaches that are based on substitution from the unconstrained least squares problem into the constraints. A useful way to describe this is via the augmented (or saddle-point) system Here λ ∈ IR p is a vector of additional variables that are often called Lagrange multipliers [17,21]. The solution x of (4.1) solves the LSE problem. Using (4.1) can be particularly useful if the constraint matrix C is dense and p is small. As we shall see in the following discussions, this is because the work involved in the proposed algorithms that depends upon p is effectively independent of the density of C. Observe that because (4.1) has a zero (2, 2) block, the augmented system can be also used to give an alternative derivation of the null-space approach of Algorithm 1. For if Z is such that CZ = 0 and x 1 is a particular solution of the second equation of (4.1) so that Cx 1 = d (steps 1 and 2 of Algorithm 1), then if x = x 1 +x, (4.1) becomes The second equation in this system is equivalent to finding x 2 such thatx = Zx 2 . Substituting this into the first equation we have Hence (see steps 3 and 4 of Algorithm 1).

Direct use of Lagrange multipliers
Algorithm 3 presents a straightforward updating scheme for solving the LSE problem using Lagrange multipliers and (4.1). Any appropriate direct or iterative method can be used for Step 1, which is usually the most expensive part of the computation. There is no dependence on C so the solution y does not need to be recomputed when C changes. The method used to solve the system with a block of p right-hand sides in Step 2 can be chosen to exploit Step 1. For example, a sparse Cholesky factorization of H may be computed in Step 1 and the factors reused in Step 2. Using existing sparse LS solvers (and a dense linear solver for the p × p at Step 5), Algorithm 3 is straightforward to implement and, from Step 6, the solution y of the unconstrained LS problem can be compared with that of the LSE.
As discussed by Golub [16] and Heath [21], a numerically superior direct method that avoids both forming the potentially ill-conditioned normal matrix H and computing the multipliers λ can be derived using the QR factorization of A. Following [40], we obtain Algorithm 4. Here P is a permutation matrix chosen to ensure sparsity of the R factor. Note that, unless b (and hence f ) changes, the Q factor need not be retained and the R factor can be reused if the constraints change but A is fixed. Results for Algorithm 4 presented in Table 4.1 confirm that the computed solution  is such that the norm of the constraints residual r c is small. We omit results for problems such as deter 5 that have p = 5 constraints because they are similar (with r c typically smaller than for the corresponding problems with p = 20).

An extended augmented system approach
An equivalent formulation of (4.1) is given by the 3-block saddle-point system (the first order optimality conditions) Applying the analysis of Section 5 of [41] to this problem yields Algorithm 5. In exact arithmetic, the main difference between the work required by Algorithms 4 and 5 is that the former involves an additional solve with RP T . For both algorithms, K is independent of b and d.

Augmented regularized normal equations
The next approach weights the constraints and uses a regularization parameter within an augmented system formulation and then aims to balance these two modifications.
Consider the weighted least squares problem (WLS) for some very large γ (γ ≫ 1). Let x LSE be the solution of the LSE problem (1.1)-(1.2). Then because lim the WLS problem can be used to approximately solve the LSE problem [27]. An obvious solution method is to solve the normal equations for (4.3) The appeal of this is that no special methods are required: software for solving standard normal equations can be used. However, for very large values of the parameter γ, the normal matrix H γ becomes extremely ill-conditioned; this is discussed in Section 4 of [9], where it is shown that the method of normal equations can break down if γ > ε −1/2 (ε is the machine precision). Furthermore, if C contains dense rows then H γ will be dense. Another possibility is to use the regularized normal equations where ω > 0 is a regularization parameter [47]. Solving (4.4) is equivalent to solving the (m + p + n) × (m + p + n) augmented regularized normal equations The spectral condition number of (4.5) is cond(A (ω, γ)) = cond(H γ + ω 2 I) and Saunders [36] shows that cond(A (ω, γ)) ≈ A γ /ω regardless of the condition of A γ . Thus using (4.5) potentially gives a significantly more accurate approximation to the pseudo solution x = A + γ b γ (where (.) + denotes the Moore-Penrose pseudo inverse of a matrix) compared to the approximation provided by solving (4.4). In [46], the parameters are set to ω = 10 −q and γ = 10 q , where q = min{k : 10 −2k ≤ ν −t }.
Here t-bit floating-point arithmetic with base ν is used.
Rewriting (4.5) using (4.3) and a conformal partitioning of y gives (4.6) This system can be solved as in [40] using a modified version of Algorithm 5. Or, eliminating y s and setting ωγ = 1, yields We can solve this system using a QR factorization of A ωI and modifying Algorithm 4. Or, ignoring the block structure, we can treat it as a sparse symmetric indefinite linear system and compute an LDL T factorization (with L unit lower triangular and D block diagonal with blocks of size 1 and 2) using a sparse direct solver such as HSL MA97 [23] that incorporates pivoting for stability with a sparsity-preserving ordering. This factorization would have to be recomputed for each new set of constraints. Alternatively, a block signed Cholesky factorization of (4.7) can be used, that is, We then obtain Algorithm 6. Note that B need not be computed explicitly. Rather, the Schur complement S may be computed using ω 2 I + CL −T L −1 C T , and w = Bz may be computed by solving Lv = z and then setting w = −Cv, and w = −B T y c may be obtained by solving Lw = C T y c .

Algorithm 6
Given ω > 0, solve the augmented system (4.7) using Cholesky factorizations. Results for Algorithm 6 for three of our test problems using a range of values of ω are given in Table 4.2. Note that here r c is computed using r c = d − Cx (rather than using r c = ω * y c ). We see that, provided ω is sufficiently small, the values of x and r are consistent with those given in Table 1.1. By replacing the Cholesky factorization of H(ω) by an incomplete factorization H(ω) ≈LL T , we can obtain a preconditioner for solving (4.7). In particular, the rightpreconditioned system is withLB T = −C T andS = ω 2 I +BB T .
As the preconditioner (4.9) is indefinite, it needs to be used with a general nonsymmetric iterative method such as GMRES [35]. A positive definite preconditioner for use with MINRES [31] can be obtained by replacing −I in (4.9) by I. MINRES has the important advantage of only requiring three vectors of length equal to the size of the linear system. GMRES results are included in Table 4.2. The GMRES convergence tolerance is taken to be 10 −11 . We see that the GMRES iteration count is essentially independent of ω. We also ran MINRES with the same settings and, while the iteration counts were again insensitive to ω, they were significantly greater than for GMRES. For problems sctap1-2r, south31 and deter3 20 the counts were 17, 772 and 56, respectively (ω = 1.0 × 10 −5 ). Our findings in Section 4 suggest that, if we require the constraints to be solved with a small residual, then an augmented system based approach combined with a QR factorization performs better (in terms of r c ) than combining it with regularization and a Cholesky factorization. Unfortunately, QR factorizations are more expensive and while strategies for computing incomplete orthogonal factorizations for use in building preconditioners have been proposed (see, for instance, [2,3,4,26,29,32,45]), the only available software is the MIQR package of Li and Saad [29] (probably because developing high quality implementations is non-trivial). In their study of preconditioners for LS problems, Gould and Scott [18,19] found that MIQR generally performed less well than incomplete Cholesky factorization preconditioners and so is not considered here.
We have made the implicit assumption that A is sparse. However, it is straightforward to extend the augmented system-based approaches to the more general case that A contains rows that are dense. For example, if A is permuted and partitioned as where A 1 is sparse and A 2 is dense, then using a conformal partitioning of y s and of b, (4.7) can be replaced by the augmented system Finally, we remark that, if we use the 3-block form (4.6) then we can follow [40], which in turn generalises the work of Carson, Higham and Pranesh [11], and obtain an augmented system approach with multi-precision refinement. This has the potential to reduce the computational cost in terms of time and/or memory, thus allowing larger problems to be solved.

Conclusions
We have considered a number of approaches for solving large-scale LSE problems in which the constraints may be dense. Our main findings can be summarized as follows: -The classical null-space method relies on computing a null-space basis matrix Z for the "wide" constraint matrix C such that Z T A T AZ is sparse. In recent work [41], we proposed how this can be achieved using a method based on a QR factorization of C with threshold pivoting. This is not straightforward to implement. Furthermore, our numerical experiments show that, in some cases, the norm r c of the constraints residual can be larger than for other approaches considered in this study. Thus, although in some contexts null-space approaches are popular, we do not recommend the strategy of [41] for LSE problems. -The direct elimination approach couples the constraint matrix and the LS matrix, leading to a sparse-dense transformed least squares problem. Existing direct or iterative methods can be used to solve the transformed problem and our experiments found the computed constraint residuals are small. The approach can be used for problems for which A (as well as C) contains a small number of dense rows. A weakness is that, if solving a sequence of problems in which either A or C is fixed, the coupling of the two blocks in the solution process means that it must be restarted. Furthermore, the number of dense rows in the transformed problem can be relatively large, making it expensive to solve. -There are several options for using an augmented system formulation. This can be solved using standard building blocks, such as a sparse QR factorization, a sparse symmetric indefinite linear solver, or a block sparse Cholesky factorization. An attraction of each of these is that existing "black box" solvers can be exploited, thereby greatly reducing the effort required in developing robust and efficient implementations. The augmented system formulation can be generalised to handle dense rows in A and offers the potential for mixed-precision computation. Moreover, an incomplete Cholesky factorization can be used as a preconditioner with a Krylov subspace solver. -In the case of a series of LSE problems in which only the constraints change, both the null-space and direct elimination approaches have the disadvantage that the computation must be redone for each new set of constraints. For the augmented system approaches, a significant amount of work can be reused from the first problem in the sequence when solving subsequent problems.
Finally, we observe that there is a lack of iterative methods and preconditioners that can be used to extend the size of LSE problems that can be solved. We have shown that using an incomplete factorization within a block factorization of an augmented system can be effective, but most current incomplete factorizations that result in efficient preconditioners are serial in nature and not able to tackle extremely large problems (but see [1,24] for novel approaches that are designed to exploit parallelise). Addressing the lack of iterative approaches is a challenging subject for future work.