A comparative study of null-space factorizations for sparse symmetric saddle point systems

Summary Null-space methods for solving saddle point systems of equations have long been used to transform an indefinite system into a symmetric positive definite one of smaller dimension. A number of independent works in the literature have identi-fied that we can interpret a null-space method as a matrix factorization. We review these findings, highlight links between them, and bring them into a unified framework. We also investigate the suitability of using null-space factorizations to derive sparse direct methods and present numerical results for both practical and academic problems.


INTRODUCTION
A saddle point system is an indefinite linear system of equations of the form (1) Here, we will assume that B ∈ R m×n (n > m) has full rank and A ∈ R n×n is symmetric positive definite on the null space of B.
We are particularly interested in the case where A and B are large, sparse matrices, and our discussion will focus on A symmetric and positive semidefinite. Saddle point systems (also known as KKT systems or augmented systems) arise frequently in scientific applications, particularly when solving constrained optimization problems. Sometimes the optimization framework is explicit, for example, when applying an interior point method to solve a linear or quadratic program, or when solving a least squares problem. More often, the optimization context is not so obvious, as is the case with incompressible fluid flow, electronic circuit simulation, and structural mechanics. The survey paper of Benzi et al. 1 describes a wide range of applications that require the solution of systems of the form 1.
One approach for solving Equation 1 is to use a null-space method. [1, section 6] These methods, which we will describe in detail below, have been used in the fields of optimization (where they are known as reduced Hessian methods), structural mechanics (where they are known as a "force" method or "direct elimination"), fluid mechanics (where they are known as the "dual variable" method), and electrical engineering (where they are known as "loop analysis"). This approach remains particularly popular in the large-scale optimization literature. [2][3][4][5][6][7] There has been a sizable body of work-some historical but much recent-that has (directly or indirectly) revisited null-space methods and put them into the context of a matrix factorization, see, for example, the references. [8][9][10][11][12][13][14][15][16][17][18] Our main contribution is to bring these factorizations together within a unified framework, highlighting the relationships between them (some of which do not appear to be well known). In addition to drawing links between different null-space factorizations, we compare their relative merits. We focus on two issues that are vital for determining the potential of a sparse direct solver, namely the sparsity of the factors and the stability of the factorization.
Null-space methods require a basis for the null space of B. This typically comes from either finding an invertible subblock of B or performing a QR factorization of B; we discuss different choices in Section 2, drawing them into a unified framework. The choice of null-space basis affects the conditioning of the resulting factorization and thus its stability. In Section 3, we highlight stability results that have appeared in the literature.
We also observe that null-space factorizations allow us to predict the level of fill a priori, which is not the case for a general sparse symmetric indefinite solver that employs numerical pivoting for stability. However, as a null-space factorization predetermines an elimination ordering, it may or may not produce sparser factors than a general sparse solver. Section 4 presents our numerical experiments that explore the effectiveness of null-space factorizations and compare their performance with that of a general sparse solver. We consider a range of problems, both artificial and from real-world applications, and investigate the stability of the factorization and the sparsity of the factors. We give concluding remarks in Section 5.

NULL-SPACE METHODS AS A FACTORIZATION
Suppose we are given a matrix Z ∈ R n×(n−m) whose columns form a basis for the null space of B, that is, BZ = 0. Suppose additionally that we have a particular solution for the second equation, that is, a vectorx such that Bx = g.

Then solving Equation 1 is equivalent to solving
where x =x +x. The second equation in this system is equivalent to finding a vector z ∈ R (n−m) such thatx = Zz. Substituting this into the first equation we have AZz + B T y = f − Ax Therefore, by solving the reduced system (Equation 2), we can straightforwardly recover x =x + Zz. We can then obtain y by solving the overdetermined system Ax + B T y = f . This is a null-space method, which we summarise as Algorithm 1.
The approach outlined in Algorithm 1 is well known. What is less well known is that we can interpret this as a matrix factorization. Consider again Equation 1. The primal variable x can be expressed in the form where Y ∈ R n×m is chosen so that and hence, and so, It is clear that this is a matrix representation of Algorithm 1. First,x = Yx R is found by solving the linear system BYx R = g. Then, the component x N of x in the null space of B is found by solving the linear system Finally, y is recovered by solving Note that the matrix is square and nonsingular and so, using Equations 3 and 4, we can rewrite Equation 1 as Thus, the factorization is a block LTL T factorization, with L lower triangular and T reverse block triangular, that is equivalent to the null-space method.
Because there are infinitely many potential bases Y, Z, this factorization is nonunique and the main difficulty of the null-space method is choosing these bases. In the Sections 2.1 to 2.5, we discuss some special cases that have been proposed in the literature.

Factorizations based on the fundamental basis
One way of fixing Y in Equation 3 is to extend B T to an n × n nonsingular matrix then it is easy to see that BZ = 0 and BY = I. The factorization in Equation 5 reduces to This factorization was given by Fletcher et al., [19, eq. 3.6] who described it as "readily observed (but not well known)."

A first null-space factorization
Suppose we have a nonsingular subset of m columns of B. We may then write, without loss of generality, , where B 1 ∈ R m×m is nonsingular. If, as suggested by Fletcher et al, 19 we make the choice of V =

This gives us the bases
This choice for Z is often called the fundamental basis, [1, section 6] and we consequently label it Z f . Substituting Equation 6 into Equation 5 gives the factorization denotes the (n − m) × (n − m) null-space matrix and with A 11 ∈ R m×m . It is easy to see that Equation 7 is equivalent to Indeed, this LTL T factorization appeared in the survey paper by Benzi et al., [1, eq. 10.35] where it was attributed to a personal communication from Michael Saunders and was described as being "related to the null-space method." We will refer to this decomposition as Factorization 1. The factor  1 is well conditioned provided B 1 is chosen appropriately; this is discussed in Section 3. Factorization 1 is then a stable factorization of  and the diagonal blocks of  1 (that is, B 1 , N, and B T 1 ) accurately reflect the condition of the full system.
In practice, once we have found the m × (n − m) matrix  by solving the sparse system the only additional calculations required to generate Factorization 1 are matrix-matrix products with , some matrix additions, and a factorization of the null-space matrix N. This is due to the relationships X = − T A 11 + A 21 and N = −X −  T A 12 + A 22 .
The overall cost is summarised in Table 1. We point out that, because the subblocks involved are sparse, it is not possible to give explicit operation counts in Table 1, as these would depend on the density of the submatrices. Note further that, in some applications (for example, where N is large and dense) it may not be possible to explicitly compute N and form its factorization; in this case, it is necessary to use an iterative solver.
Having formed Factorization 1, Table 2 gives the costs of using it to solve Equation 1. Again, note that exact operation counts are not possible. Table 2 describes two variants: implicit and explicit. In the explicit version, we store the off-diagonal matrices (X, ) that are formed in the construction of N and apply them via matrix-vector products to solve Equation 1. In the implicit version, we discard the off-diagonal matrices and recompute them as needed. This is computationally more expensive but saves storing the potentially dense matrices X and  of size m × (n − m). It is clear that the implicit version is exactly equivalent to the null-space method as presented in Algorithm 1.
The costs in Tables 1 and 2 are upper bounds, and they may be reduced in certain circumstances. For example, as we discuss in Section 3 below, it is usual to find B 1 , B 2 by forming an LU factorization of B T where L a is lower triangular, U a is upper triangular, and P and Q are permutation matrices. Then B −1 1 B 2 = PL −T a L T b P T , and so this can be calculated without reference to U a , although L b is needed and is likely to be less sparse than B 2 . Furthermore, in this , so Z T AZ can also be formed efficiently. See, for example, Fletcher et al. 19 for more details.
As explicit a plus: Note. mat-vec denotes the product of a sparse matrix with a vector, and vec-add denotes the addition of two vectors. a The numbers of solves and matrix-vector products will be the same, but the matrices in the matrix-vector products will generally be sparser in the implicit case.

A factorization due to Lungten, Schilders, and Maubach
Assume now that A is symmetric and positive semidefinite so that N is symmetric positive definite and a Cholesky factorization of the form N = L 2 L T 2 exists, where L 2 is lower triangular. Then we can decompose the reverse triangular  1 matrix in Equation 9 as with L A a strictly lower triangular matrix and D A a diagonal matrix. Combining the outer matrices here with the outer matrices in Equation 9 yields the alternative, but equivalent, LTL T factorization where This factorization was recently proposed both for use as a direct method and as the basis of a preconditioner for an iterative method by Lungten et al. 15 We refer to it as Factorization 2 or the LSM factorization.
Note that forming Equation 12 is more expensive than Equation 9, as it requires one more matrix-matrix multiply of  T (recall Equation 10) with an m × m matrix and one more (n − m) × m matrix addition, both coming from the formation of K. In terms of applying Equation 12 explicitly, one matrix-vector multiply with A 11 is replaced by matrix-vector multiplies with its strictly upper, lower, and diagonal parts, and two extra m × m matrix additions; there is a similar increase in cost when applying the factorization implicitly. This suggests that, in terms of the computational cost, Factorization 1 is preferable; we perform tests with both versions in Section 3.
Lungten et al. 15 focus on problems for which the nonsingular matrix B 1 is also upper triangular (or it is easy to transform the problem into this form). In this case, if Factorization 2 is formed via an modified Cholesky algorithm, they show that, for the dense case, it takes flops to factorize the saddle point matrix .

The Schilders factorization
We next consider the relationship of Factorizations 1 and 2 to the decomposition known in the literature as the Schilders factorization. Dollar et al. 10 were interested in developing constraint preconditioners for symmetric systems as in Equation 1. Such preconditioners represent the blocks B exactly but approximate the (1,1) block A. Dollar et al. choose symmetric matrices and M ∈ R (n−m)×m . To obtain an inexpensive preconditioner, E 2 , F 2 are chosen so that they are easy to invert. The Schilders factorization is then given by where Note that the (1,1) block A is implicitly defined by the choices of E i , F i . Nevertheless, we can use this construction to give a factorization for a given A. One possible choice for E 1 , F 1 is . The matrices M, E 2 , and F 2 are then given by the relations With these choices, transferring a factor of the block diagonal matrix with diagonal blocks B T 1 , I, and B −1 1 from the left outer matrix to the central matrix in Equation 13 again gives Factorization 2.
The original Schilders factorization, 14 of which the formulation 13 is a generalization, was given only for matrices for which B 1 is upper triangular and used the choice where diag() and lower() denote the diagonal and strictly lower triangular parts of a matrix, respectively. Again, it is straightforward to show the equivalence of this factorization to Factorization 2. Generating this factorization is clearly significantly more work than Equation 9, not least because we are unable to reuse the matrix  in forming the subblocks. There are, of course, other ways of rearranging the decompostion given by Equation 9. Dollar et al. 20 give a list of forms the factorization can take. As already observed, their focus was on implicit factorizations of constraint preconditioners and only five of their factorizations are applicable in the case of an arbitrary symmetric (1,1) block.

Relationship to the Schur complement factorization
A commonly used block LDL T factorization for generalized saddle point systems (which have a negative definite (2,2) see, for example, Benzi et al. [1,Equation 3.9] Approximating the terms of this factorization has proved an effective strategy for developing preconditioners for saddle point systems. It can also be used to develop another factorization that is equivalent to the null-space method. First, note that where we are again assuming, without loss of generality, that B 1 is a nonsingular m × m subblock of B. Applying Equation 14 with C = A 22 , we obtain Here, S denotes the Schur complement, which satisfies It follows that the null-space matrix in Equation 8 is the Schur complement for an alternative blocking of the matrix, and we have the factorization Again, this can be derived from Equation 9 by simply noting that the reverse triangular matrix  1 is equal to the product . This factorization is therefore not "new"-and it is significantly more expensive to form than Equation 9-but it highlights the connection between Schur complement methods and null-space factorizations. In particular, given the success in finding approximations to the Schur complement using techniques from functional analysis (see, e.g., Mardel et al. 21 ) it is hoped that viewing the null-space matrix this way could yield alternative preconditioners for certain classes of saddle point systems.

Connection with Cholesky decomposition
Schilders developed his original factorization and the subsequent variant (Equation 12), by considering what he terms a microfactorization. In this formulation, the matrix in Equation 1 is reordered by pairing every entry on the diagonal of the (1,1) block A with a corresponding nonzero entry in the constraint block B so that (after permutations) the entries on the diagonal form micro saddle point systems. This is known as a tiling in the optimization community. Below is an illustrative example of this ordering for n = 3, m = 2: Note that there is no requirement for a ij to be combined with b ij and b ji ; instead, a suitable pairing that preserves sparsity and is numerically stable should be chosen-see Section 3 for further discussion. This is an example of a constrained ordering. [22][23][24] Because the entries on the (block) diagonal are now micro saddle point systems (chosen to be nonsingular) a modified sparse Cholesky code can be used to solve this system, and it is guaranteed (at least in exact arithmetic) that this will not break down. [14, section 3] By this process, a factorization = LDL T can be computed, where D has 1×1 and 2×2 blocks on the diagonal in the appropriate places. Furthermore, uniqueness results in, for example, Maubach et al., 12 show that the factors generated by the Cholesky process will be equivalent to those generated by the (macro) factorizations described earlier in this section.
In addition to the work by Schilders et al., 12,14,15 this approach has been considered by Forsgren et al. 25 (whose focus was on inertia control), Gould, 17 and de Niet et al. 9 ; each of these works, to varying degrees, made the connection between the microfactorization and the null-space method.

The antitriangular factorization
An alternative matrix factorization can be obtained by assuming we have a QR factorization where is an n × n orthogonal matrix and R is an m × m upper triangular and nonsingular matrix. Then Q 1 spans the range of B T and Q 2 spans the null space of B. We can, therefore, substitute Y = Q 1 and Z = Q 2 into the factorization given by Equation 5 to obtain We refer to this as Factorization 3 or the AT (for antitriangular) factorization. The reverse triangular matrix  3 in Equation 17 appeared in the proof of theorem 2.1 in Keller et al. 18 which established eigenvalue bounds for Equation 1 preconditioned with a constraint preconditioner, although the link to a factorization of  does not appear to have been made there. Recently, Pestana et al. 13 derived Factorization 3 in the case A is symmetric and B is of full rank and showed that it is (up to a permutation) the representation of the antitriangular factorization of Mastronardi et al. 11 applied to a saddle point system. They also gave a flop count for forming the factorization in the dense case. The work is dominated by performing the QR factorization (Equation 16), with additional work being two n × n matrix-matrix products and a subsequent factorization of Q T 2 AQ 2 . It is clear from the above formulation that-provided a rank-revealing QR is applied-Factorization 3 is well defined for a B with linearly dependent rows and also for nonsymmetric A.
An advantage of having an orthogonal null basis is that the reduced system (Equation 4) is guaranteed to be well conditioned if A is. However, even for sparse problems, Q 1 and Q 2 may be rather dense. The explicit version of the factorization requires storing Q 1 , Q 2 , Q T 1 AQ 2 and the lower triangular part of Q T 1 AQ 1 in addition to keeping the factors of Q T 2 AQ 2 . For large systems this incurs prohibitive storage costs, see Section 3. Sparse QR routines-for example, SuiteSparseQR 26 -normally allow the user to compute the action of Q i on a vector using the stored Householder reflectors, and this facility is the preferred option here.
Pestana et al. 13 give a complexity analysis of the factorization in Equation 17 in the dense case and show that the number of flops required is An alternative way of obtaining the factorization in Equation 17 which again focuses on the dense case, is given by Mastronardi et al., 27 and the number of flops there is dependent on the size of m relative to n; see Pestana et al. [13,

Other bases and converting between factorizations
Consider again the general factorization as given in Equation 5. As already observed, there is an infinite number of potential matrices Z (whose columns span the null space of B) and Y (whose columns are such that [ Y Z ] span R n ). In fact, for any given basis matrix Z and for any nonsingular matrix G ∈ R (n−m)×(n−m) , another factorization that is equivalent to Equation 5-and hence Algorithm 1-is whereŶ is a matrix such that To highlight this fact, we will show explicitly how we can transform Factorization 1 into Factorization 3. Consider the fundamental basis Z f (Equation 6), and assume that we have its "skinny" QR factorization Z f = Q 2 R 2 . Then, we see that the columns of which clearly span the null-space of B. Now, let B T = Q 1 R 1 be the skinny QR factorization of B T . Because B is of full rank, we have obtained two matrices, Q 1 ∈ R m×n that spans the range of B T and Q 2 ∈ R n×(n−m) that spans the null-space of B. We, therefore, know that the n × n matrix [ spans R n . Substituting G = R −1 2 andŶ = Q 1 into Equation 18, we see that these choices transform Factorization 1 (Equation 9) into Factorization 3 (Equation 17). This construction may not have any practical use, but it is theoretically satisfying to explicitly connect two seemingly different methods. In the same way, by an appropriate choice of G andŶ, we can transform any null-space factorization into any other.

STABILITY AND SPARSITY PROPERTIES
For a factorization to be effective for solving sparse linear systems, it must be stable in the classical sense [29, chapter 7] and the factors must be sparse. In this section, we bring together results from the literature that concern the viability of the factorizations described in Section 2 as the basis of a sparse direct solver.
Consider Factorization 3, the antitriangular factorization. This was shown to be backward stable by Mastronardi et al. 11 Stability is a direct consequence of the fact that the only calculation needed is a stable QR factorization, see also Gill et al., 16 and Fletcher et al. 19 The penalty for this stability is that the factors may not be sparse, even when using a sparse QR factorization routine, as the Q matrices can fill in significantly; this is illustrated in Section 4.1.
An alternative approach is to choose the permutations in Factorizations 1 and 2 in such a way that we guarantee a stable factorization. Forsgren et al. 25 -thinking in terms of a microfactorization-describe a pivoting strategy to do this. Their method was refined further by Gould. 17 Another way to ensure stability is to perform an LU factorization on B T (incorporating pivoting to ensure the L factor is well conditioned). If we obtain B 1 using this, then B T 2 B −T 1 is bounded, and the factorization can be proved to be backward stable, whereas the forward error in x and y are proportional to ([B T V T ]) (Z T AZ) and ([B T V T ]) 2 (Z T AZ), respectively, where denotes the condition number-see Fletcher et al. 19 for details; an alternative description is given by Gill et al. 30 This is the approach we adopt in the tests in Section 4.
A weakness of a general solver, such as an LDL T factorization, is that, due to the need for pivoting to preserve stability, the fill in the computed factors may be significantly higher than predicted on the basis of the sparsity pattern of the matrix alone. This use of numerical pivoting not only leads to a higher operation count and denser factors but also prohibits the exploitation of parallelism (and significantly complicates the code development). On the other hand, null-space factorizations reduce the solution of Equation 1 to solving a positive definite matrix, and as such, they give a direct method with a predictable level of fill without pivoting (other than to find B 1 ). This is a potentially attractive feature. However, by identifying a nonsingular subblock B 1 of B we are essentially predetermining an ordering that may or may not produce sparser factors than that used by a general sparse indefinite direct solver.
Instead of choosing B 1 to give stability at the expense of sparsity, we can alternatively choose B 1 so that the fundamental basis Z f is sparse, for example, by ensuring B −1 1 is sparse. Pinar et al. 31 describe several ways to do this using methods based on graph matchings and hypergraph partitioning. However, such approaches are reportedly time consuming, and no results about the stability of the resulting decomposition are given. Murray 32 describes a method for picking Z so that Z T AZ is diagonal. However, to quote Fletcher, [33, section 10.1] this "may be doubtful in terms of stability." Note also that B −1 1 can only be sparse when the graph of B 1 is disconnected, [34, section 12.6] which may not be possible to achieve in certain applications.
Favourable sparsity and stability properties are achievable for certain classes of matrices . In particular, de Niet et al. 9 prove some strong results for  matrices. An  matrix is a saddle point matrix (Equation 1) in which A is symmetric positive definite and B is a gradient matrix, that is, B has at most two entries per row, and if there are two entries, their sum is zero. Such matrices appear naturally in certain formulations of fluid flow problems (see, for example, Arioli et al. 35 ) and also in electronic simulations, where they arise as the incidence matrix of the graph representing the circuit. 36 The original Schilders factorization 14 was developed specifically for solving systems with this structure.
De Niet et al. find the basis B 1 implicitly by considering a microfactorization (recall Section 2.3). They pair entries in A with entries in B in such a way as to ensure the stability of the factorization. Moreover, de Niet et al. show that the number of nonzeros in the symbolic factorization of F(A) ∪ F(BB T ), where F(·) denotes the sparsity pattern, "provides a reasonable estimate" of the number of nonzeros in the factorization. We perform tests on this class of matrices in Section 4.3.

NUMERICAL RESULTS
With the exception of the special case of  matrices, the literature focuses on proving either stability, with little said about sparsity or vice versa. In this section, we present results using Factorizations 1, 2, and 3-that is, Equations 9, 12, and 17, respectively-as the basis of a direct method. Our intention is to demonstrate the stability and sparsity of these factorizations when applied to a range of problems, and we compare their performance with that of a general sparse symmetric indefinite solver. We concentrate on the solution of a single saddle-point system and do not exploit the potential advantages null-space factorizations can offer over a general indefinite solver in terms of updating the factors after a row is added to (or removed from) .
We perform tests using MATLAB Version R2014a. We factorize the positive definite null-space matrix N using cholmod 37 and use the SuiteSparseQR algorithm, 26 both from Tim Davis' SuiteSparse package and applied via a MATLAB interface. To find an invertible subset of B, we use the package LUSOL 38 with partial pivoting (with a pivot tolerance of 1.9, as suggested for basis repair in the LUSOL documentation) to perform the sparse factorization of B T . LUSOL is called via a MATLAB interface. 39 For comparison, results for the MATLAB command ldl (with default settings) are also given. ldl computes an LDL T factorization (where L is unit lower triangular and D is block diagonal with 1 × 1 and 2 × 2 blocks on the diagonal) of a sparse symmetric indefinite matrix by employing the direct solver 40 MA57 from the HSL mathematical software library. 41 MA57 first computes a sparsity-preserving ordering of  and a scaling that is designed to help with numerical stability. It does not explicitly exploit the block structure of  but is designed to efficiently and robustly solve general symmetric indefinite sparse linear systems. To ensure a fair comparison in our experiments, we always prescale  using the scaling for  calculated by ldl. We recognise, however, that, in practice, it may be more effective to scale the A and B blocks separately as x and y may represent very different types of quantity.
To measure stability, we report the scaled backwards error for Equation 1 given by , see Higham. [29,Equation 7.23] Again, we note that in certain circumstances it may be more appropriate to consider the residual of the first and second rows of Equation 1 individually. In particular, note that we can bound Equation 19 by ||Ax + B T y − f|| 2 2 ∕||b|| 2 2 + ||Bx − g|| 2 2 ∕||b|| 2 2 ⩽ ||w − b|| 2 2 ∕||b|| 2 2 . We highlight the fact that, due to the nature of null-space factorizations, ||Bx − g|| 2 2 is generally well behaved and close to machine precision. This is because the LU or QR factorization used to find the null space of B will be accurate. When comparing null-space factorizations, we can, therefore, assume that most of the error in ||w − b|| 2 actually comes from the term ||Ax + B T y − f || 2 . In contrast, MA57 does not treat the two types of variable differently. Because we assume no prior knowledge of the source of the matrix in the following examples, we present the results as residuals for the unreduced matrix only.
We require a measure of the amount of fill that a method uses; this will be different in the explicit and implicit cases. For the explicit methods, the measure of fill we use is where  (·) denotes the number of nonzeros needed to store the information required to solve with  i . For all the factorizations considered here,  ( i ) denotes the number of nonzeros in the leftmost matrix in Factorization i. The central  i is symmetric, and we only count the lower triangular part. Here,  ( i ) denotes the number of nonzeros of all blocks that only need to be multiplied, plus the number of nonzeros in the factors of any matrix that needs to be solved. where tril(·) denotes the lower triangular part of a matrix. For implicit methods, we only store the number of nonzeros needed by the factorizations of the blocks we must solve for; we assume that we have access to the original matrix, and do not count entries that can be obtained from there. For example, for Factorization 1, the measure of fill would be nnz(L B ) + nnz(U B ) + nnz(L N ), as all the other blocks can be obtained from the original matrix. Recall that, although the storage is lower for implicit methods, the operation count is higher, as certain quantities are computed "on the fly"; see Table 2.

An academic example
We start by considering some small academic matrices that illustrate some worst case behaviour. We form a system of the form given in Equation 1 using the MATLAB commands: 0.1,1e-1,2) We take n = 1024 and m = {100, 512, 900} and construct the matrices so that the leading m columns of B form a nonsingular submatrix; these can be used without permutation as B 1 . Figure 1 presents plots showing the backward error and the fill. We should expect the norm of the residual (Equation 19) to be close to machine precision if the method is stable. However, we observe that, if B 1 taken as being the leading columns of B, Factorizations 1 and 2 are numerically unstable, with the measure of the backwards error spanning four orders of magnitude, and, at worst, greater than 10 −9 . This indicates that the accuracy of the solution is variable and hence demonstrates that choosing an invertible subset for B 1 is not sufficient for stability. If, however, we choose B 1 using LUSOL with partial pivoting, as expected given the results quoted in Section 3, the norms of the residuals obtained using Factorizations 1 and 2 all lie below 10 −13 , confirming that this method is stable. Factorization 3 (the antitriangular factorization) also behaves as expected, being the most stable of the null-space factorizations, although it also results in the largest amount of fill (Figure 1b).

Optimization examples
We next consider examples from the optimization community; namely, a selection of relatively small quadratic programs from the CUTEst 42 test collection. In particular, the problems are taken from the Maros and Meszaros quadratic programming test set and are solved without consideration for any inequality constraints. These are of the form We chose the subset of problems for which A is symmetric positive semidefinite, has at least n − m nonzeros on the diagonal, and where n > m. Our test problems, and the relative sizes of their submatrices, are listed in Table 3. Our tests showed that the families {HS21, HS35, HS35MOD, HS51, HS52, HS53}, {AUG2D, AUG2DC, AUG2DCQP, AUG2DQP}, and {AUG3D, AUG3DC, AUG3DCQP} all exhibited very similar behaviour, and so we only report results for a representative problem from each of these sets: HS51, AUG2DC, and AUG3DC. We again select the basis B 1 using LUSOL with partial pivoting. We do not compare Factorization 3 here, and in the subsequent tests, as the cost of doing a QR factorization proved to be excessive (both in terms of timing and storage) for many of these problems.
The results are shown in Figures 2 and 3. The backwards errors for the null-space factorizations and the ldl factorization are generally similar, with the null-space factorizations giving slightly better accuracy on the CONT-xxx problems, which are known to be difficult for direct solvers. If we apply one step of iterative refinement (as is usual in practical applications) then, as we see in Figure 2b, all methods behave comparably.
In Figure 3, we compare the sparsity of the factors. The explicit versions of Factorizations 1 and 2 result in similar levels of fill, and, with a few exceptions, these are greater than for ldl, while the implicit version of Factorization 1 gives significant storage savings (at the cost of more computation). The latter needs the least storage for 38 of the 53 problems, whereas ldl performs best for the remaining 18 problems, in some cases, producing significantly sparser factors than the null-space factorizations. For example, about 100 times more storage is required by Factorization 1 (implicit) compared to ldl for the AUGxx problems, and over 1,500 times the storage is needed for HUES-MOD and HUESTIS.
We do not report timings, as we cannot give a fair comparison between our proof-of-concept MATLAB code and the compiled ldl command, but we observed the null-space factorizations to be significantly slower. They do, however, potentially offer . The legend gives the order of performance, running from least sparse to most sparse (on average) scope for parallelization, as the bulk of the computational work is from forming and factorizing N. The columns of N can be computed in parallel and, because the Cholesky factorization does not require pivoting, it can be parallelized more effectively than an indefinite sparse solver (see, e.g., Hogg et al. 43 ). However, it is still necessary to factorize the nonsquare matrix B T to obtain the basis B 1 .

 matrices
Recall that both stability and sparsity can be shown for  matrices. We give two examples of such problems below: one from resistor networks and one from fluid flow. In the former, the (1,1) block is a diagonal matrix, and in the latter, it is nondiagonal but symmetric positive definite.

Resistor networks
An important problem arising from the electronics industry is to find the voltage, V, and current, I, of a network of resistors. The current and voltages are related by the equation AI + BV = 0, where A is a diagonal matrix containing the values of the resistors, and B is the incidence matrix of the network. From Kirchhoff's law, we also have B T I = 0. One node is usually grounded (algebraically, we remove a column of B) so that B has full rank. It is clear that putting these two equations together we get a system of the form 1 that is also an  matrix. For further details, we refer the reader to Rommes et al. 44 In order to analyse the behaviour of null space factorizations on such problems, we run tests on artificial systems with n = 1, 024 resistors, joined at m = {100, 250, 512} nodes at random (while forming a complete network). The resistor values are chosen at random from a uniform distribution between 0 and 10 −2 . Plots showing backward errors and fill are given in Figure 4. A property of matrices of this type is that it is possible to permute B to make it upper trapezoidal, and so a well-conditioned, easy to invert, block B 1 is possible to achieve without arithmetic (see Benzi et al. [1, chapter 6] for a discussion and references).
The results illustrate that the null-space factorizations are stable (as predicted by the theory in de Niet et al. 9 ); indeed, they produce smaller backward errors than ldl (without iterative refinement). In terms of sparsity, the antitriangular factorization (Factorization 3) is again the poorest, whereas Factorization 1 gives the sparsest factors. As we do not need to factorize a block of B, the fill for the factorizations based on the fundamental basis is negligible.

Fluid flow problems
Fluid flow problems provide another source of  matrices. In particular, we use some matrices derived from the mixed hybrid finite element approximation of Darcy's law and the continuity equation, which describes fluid flow through porous media. 45 Our test matrices * are described in Table 4; the same examples were used as test cases by, for example, Tůma 46 and de Neit et al. 9 Sparsity and stability results are given in Figure 5. As with example 4.2, we do not test Factorization 3 here, as computing a QR factorization proved to be too expensive. *We thank Miroslav Tůma for providing us with these matrices.   For these real-word problems, the null-space factorizations perform markedly worse than ldl. The fill, in particular, is very high. Consider the problem dan2, which has 318,750 nonzeros in A and 127,054 nonzeros in B. For Factorization 1, if Z T AZ = LL T , then the number of nonzeros in L is 59,146,699. This is an order of magnitude larger than any of the other terms in the fill calculation and is due to the fact that the null-space matrix, Z T AZ, is significantly denser than A. This means that the null-space factorizations can need more than 20 times the storage of ldl, as ldl has greater freedom to choose appropriate pivots to preserve sparsity in the LDL T factorization. Furthermore, null-space factorizations also result in larger backwards errors (although a single step of iterative refinement reduces all the backwards errors to the order of machine precision).

CONCLUSION
Null-space methods for solving systems of the form given by Equation 1 can be thought of in terms of matrix factorizations. We have presented a number of such decompositions that have been proposed in the literature. By placing these in the unified framework of null-space factorizations we see that, although they may appear to be different, they are in fact fundamentally related.
By highlighting such relationships it is clear that, for example, stability results proved for one type of factorization are equally applicable to another. In particular, we can say that all factorizations relying on a fundamental basis enjoy a certain degree of stability and sparsity when applied to  matrices. Similarly, the stability of Factorization 2 can be shown if B 1 is chosen appropriately, a result that is only clear by putting it in the null-space factorization framework.
In the literature we find null-space factorizations that are known to be either stable or sparse, but (with the exception of  matrices) not both. However, a practical factorization must possess favourable properties in both these areas. In Section 4, we investigated numerically the stability and sparsity properties of various null-space factorizations using matrices derived from a number of academic and practical applications. Our results suggest that these factorizations-particularly the form in Equation 9-have the potential to be competitive with a sparse symmetric indefinite solver for some classes of problems.
However, we stress that any set of numerical tests can only be indicative of typical behaviour, and cannot predict worse case performance. In the future, we would like to see the development a null-space factorization for general saddle point systems that gives, as much as possible, a theoretical guarantee of a certain measure of both stability and sparsity. Our hope is that, by exposing connections between previously disparate solution methods, we have brought such an algorithm one step closer.