We first consider bounds on for . It is possible to obtain a bound from first principles by using the relation , where is the adjugate matrix, along with the fact that since has integer entries. Higham and Lettington (2021) found that the smallest bound they could obtain came from a bound of Merikoski et al. (1997): for nonsingular ,

Applying this bound to , using the fact that is monotonically increasing for , gives

Another result from Merikoski et al. (1997) gives, for symmetric positive definite ,

For , since we have , and hence

Recall that Rutishauser’s bound is . The bounds (1) and (2) remain valid if we modify the definitions of and to allow zero elements (note that Rutishauser’s matrix has a zero element).

The sets and are large: has on the order of elements. Exhaustively searching over the sets in reasonable time is possible with a carefully optimized code. Higham and Lettington (2021) use a MATLAB code that loops over all symmetric matrices with integer elements between and and

• evaluates from an explicit expression (exactly computed for such matrices) and discards if the matrix is singular;
• computes the eigenvalues of and obtains the condition number as (since is symmetric); and
• for , checks whether is positive definite by checking whether the smallest eigenvalue is positive.

The code is available at https://github.com/higham/wilson-opt.

The maximum over is attained for

which has . and determinant . The maximum over is attained for

which has and determinant . Obviously, symmetric permutations of these matrices are also optimal.

The following table summarizes the condition numbers of the matrices discussed and how close they are to the bounds.

Matrix	Comment
	Wilson matrix		99.73	13.40
	Rutishauser’s matrix		8.31	1.12
	By random sampling		6.19	—
	Optimal matrices in		3.91	—
	Optimal matrices in		8.38	1.13

Clearly, the bounds are reasonably sharp.

We do not know how Wilson constructed his matrix or to what extent he tried to maximize the condition number subject to the matrix entries being small integers. One possibility is that he constructed it via the factorization in the next section.

The Cholesky factor of the Wilson matrix is

Apart from the zero element, it is unremarkable. If we factor out the diagonal then we obtain the factorization, which has rational elements:

Suppose we drop the requirement of triangularity and ask whether the Wilson matrix has a factorization with a matrix of integers. It is known that every symmetric positive definite matrix of integers with determinant has a factorization with an matrix of integers as long as , but examples are known for for which the factorization does not exist. This result is mentioned by Taussky (1961) and goes back to Hermite, Minkowski, and Mordell. Higham and Lettington (2021) found the integer factor

of , which is block upper triangular so can be thought of as a block Cholesky factor. Higham, Lettington, and Schmidt (2021) draw on recent research that links the existence of such factorizations to number-theoretic considerations of quadratic forms to show that for the existence of an integer solution to it is necessary that a certain quadratic equation in variables has an integer solution. In the case of the Wilson matrix the equation is

The authors solve this equation computationally and find and two rational factors:

They show that these matrices are the only factors of up to left multiplication by integer orthogonal matrices.

The Wilson matrix has provided sterling service throughout the digital computer era as a convenient symmetric positive definite matrix for use in textbook examples and for testing algorithms. The recent discovery of its integer factorization has led to the development of new theory on when general integer matrices can be factored as (when is symmetric positive definite) or (a problem also considered in Higham, Lettington, and Schmidt (2021)), with integer .

Olga Taussky Todd wrote in 1961 that “matrices with integral elements have been studied for a very long time and an enormous number of problems arise, both theoretical and practical.” We wonder what else can be learned from the Wilson matrix and other integer test matrices.

This is a minimal set of references, which contain further useful references within.

• Nicholas J. Higham and Matthew C. Lettington, Optimizing and Factorizing the Wilson Matrix, to appear in Amer. Math. Monthly, 2021.
• Nicholas J. Higham, Matthew C. Lettington, and Karl Michael Schmidt, Integer Matrix Factorisations, Superalgebras and the Quadratic Form Obstruction, Linear Algebra Appl. 622, 250–267, 2021.
• Cleve B. Moler, Reviving Wilson’s Matrix, 2018.
• H. Rutishauser, On Test Matrices, 7 (no. 165), 349-365, in: Programmation en Mathematiques Numeriques, Besancon, 1968.
• Olga Taussky, Some Computational Problems Involving Integral Matrices, Journal of Research of the National Bureau of Standards Section B—Mathematical Sciences 65(1), 15–17, 1961.

• What Is a Cholesky Factorization? (2020)
• What Is a Condition Number? (2020)
• What Is a Symmetric Positive Definite Matrix? (2020)

This article is part of the “What Is” series, available from https://nhigham.com/category/what-is and in PDF form from the GitHub repository https://github.com/higham/what-is.

An RRF helps to determine the numerical rank, which we now define.

Definition 2. For a given the numerical rank of is the largest integer such that .

By the Eckart–Young theorem, the numerical rank is the smallest rank attained over all with . For the numerical rank to be meaningful in the sense that it is unchanged if is perturbed slightly, we need not to be too close to or , which means that there must be a significant gap between these two singular values.

One might attempt to compute an RRF by using a QR factorization , where has orthonormal columns, is upper triangular, and we assume that . In Definition 1, we can take

However, it is easy to see that QR factorization in its basic form is flawed as a means for computing an RRF. Consider the matrix

which is a Jordan block with zero eigenvalue. This matrix is its own QR factorization (), and the prescription gives , so . The essential problem is that the diagonal of has no connection with the nonzero singular values of . What is needed are column permutations: for the permutation matrix that reorders to , and this is a perfect RRF with .

For a less trivial example, consider the matrix

Computing the QR factorization we obtain

R =
  -2.0000e+00   5.0000e-01  -2.5000e+00   5.0000e-01
            0   1.6583e+00  -1.6583e+00  -1.5076e-01
            0            0  -4.2640e-09   8.5280e-01
            0            0            0  -1.4142e+00

The element tells us that is within distance about of being rank deficient and so has a singular value bounded above by this quantity, but it does not provide any information about the next larger singular value. Moreover, in , is of order for this factorization. We need any small diagonal elements to be in the bottom right-hand corner, and to achieve this we need to introduce column permutations to move the “dependent columns” to the end.

A common method for computing an RRF is QR factorization with column pivoting, which for a matrix with computes a factorization , where is a permutation matrix, has orthonormal columns, and is upper triangular and satisfies the inequalities

In particular,

If with then we can write

with

Hence is within -norm distance of the rank- matrix . Note that if is partitioned conformally with in (4) then

so , which means that provides an approximation to the range of .

To assess how good an RRF this factorization is (with ) we write it as

Applying (1) gives

where , since has orthonormal columns and so has unit singular values. Now has unit diagonal and, in view of (2), its off-diagonal elements are bounded by . Therefore . On the other hand, by Theorem 1 in Bounds for the Norm of the Inverse of a Triangular Matrix. Therefore

The lower bound is an approximate equality for small for the triangular matrix

devised by Kahan, which is invariant under QR factorization with column pivoting. Therefore QR factorization with column pivoting is not guaranteed to reveal the rank, and indeed it can fail to do so by an exponentially large factor.

For the matrix , QR with column pivoting reorders to and yields

R =
  -3.0000e+00   3.3333e-01   1.3333e+00  -1.6667e+00
            0  -1.6997e+00   2.6149e-01   2.6149e-01
            0            0   1.0742e+00   1.0742e+00
            0            0            0   3.6515e-09

This suggests a numerical rank of for (say). In fact, this factorization provides a very good RRF, as in we have .

Consider a factorization with triangular factor partitioned as

We have

where (7) is from singular value interlacing inequalities and (8) follows from the Eckart-Young theorem, since setting to zero gives a rank- matrix. Suppose has numerical rank and . We would like to be able to detect this situation from , so clearly we need

In view of the inequalities (7) and (8) this means that we wish to choose maximize and minimize .

Some theoretical results are available on the existence of such QR factorizations. First, we give a result that shows that for the approximations in (9) can hold to within a factor .

Theorem 1. For with there exists a permutation such that has the QR factorization with and , where .

Proof. Let , with and let be such that satisfies . Then if is a QR factorization,

since , which yields the result.

Next, we write , where , and partition

with . Then

implies . On the other hand, if is an SVD with , , and then

so

Finally, we note that we can partition the orthogonal matrix as

and the CS decomposition implies that

Hence , as required.

Theorem 1 is a special case of the next result of Hong and Pan (1992).

Theorem 2. For with and any there exists a permutation matrix such that has the QR factorization where, with partitioned as in (6),

where .

The proof of Theorem 2 is constructive and chooses to move a submatrix of maximal determinant of to the bottom of , where comprises the last columns of the matrix of right singular vectors.

Theorem 2 shows the existence of an RRF up to the factor , but it does not provide an efficient algorithm for computing one.

Much work has been done on algorithms that choose the permutation matrix in a different way to column pivoting or post-process a QR factorization with column pivoting, with the aim of satisfying (9) at reasonable cost. Typically, these algorithms involve estimating singular values and singular vectors. We are not aware of any algorithm that is guaranteed to satisfy (9) and requires only flops.

By applying Householder transformations on the right, a QR factorization with column pivoting can be turned into a complete orthogonal decomposition of , which has the form

where is upper triangular and and are orthogonal. Stewart (1998) calls (6) with upper triangular or lower triangular a UTV decomposition and he defines a rank-revealing UTV decomposition of numerical rank by

The UTV decomposition is easy to update (when a row is added) and downdate (when a row is removed) using Givens rotations and it is suitable for parallel implementation. Initial determination of the UTV decomposition can be done by applying the updating algorithm as the rows are brought in one at a time.

Instead of QR factorization we can build an RRF from an LU factorization with pivoting. For with , let

where and are permutation matrices, and are lower and upper triangular, respectively, and and are . Analogously to (7) and (8), we always have and . With a suitable pivoting strategy we can hope that and .

A result of Pan (2000) shows that an RRF based on LU factorization always exists up to a modest factor . This is analogue for LU factorization of Theorem 2.

Theorem 3 For with and any there exist permutation matrices and such that

where is unit lower triangular, is upper triangular, and

where .

Again the proof is constructive, but the permutations it chooses are too expensive to compute. In practice, complete pivoting often yields a good RRF.

In terms of Definition 1, an RRF has

For the matrix (), the factor for LU factorization without pivoting is

U =
   1.0000e+00   1.0000e+00   1.0000e-08            0
            0  -2.0000e+00   2.0000e+00   1.0000e+00
            0            0   5.0000e-09  -1.5000e+00
            0            0            0  -2.0000e+00

As for QR factorization without pivoting, an RRF is not obtained from .. However, with complete pivoting we obtain

U =
   2.0000e+00   1.0000e+00  -1.0000e+00   1.0000e+00
            0  -2.0000e+00            0            0
            0            0   1.0000e+00   1.0000e+00
            0            0            0  -5.0000e-09

which yields a very good RRF with and .

QR factorization with column pivoting is difficult to implement efficiently, as the criterion for choosing the pivots requires the norms of the active parts of the remaining columns and this requires a significant amount of data movement. In recent years, randomized RRF algorithms have been developed that use projections with random matrices to make pivot decisions based on small sample matrices and thereby reduce the amount of data movement. See, for example, Martinsson et al. (2019).

This is a minimal set of references, which contain further useful references within.

• Shivkumar Chandrasekaran and Ilse C. F. Ipsen, On Rank-Revealing Factorisations, SIAM J. Matrix Anal. Appl. 15 (2), 592–622, 1994
• Y. P. Hong and C.-T. Pan, Rank-Revealing QR Factorizations and the Singular Value Decomposition, Math. Comp. 58 (197), 213–232, 1992.
• P. G. Martinsson, G. Quintana-Orti, and N. Heavner, randUTV: A Blocked Randomized Algorithm for Computing a Rank-Revealing UTV Factorization, ACM Trans. Math. Software 45(1), 4:1–4:26, 2019.
• C.-T. Pan, On the Existence and Computation of Rank-Revealing LU Factorizations, Linear Algebra Appl. 316, 199–222, 2000.
• G. W. Stewart, Matrix Algorithms. Volume I: Basic Decompositions, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 1998.

• What Is the CS Decomposition? (2020)
• What Is an LU Factorization? (2021)
• What Is a QR Factorization? (2020)
• What Is the Singular Value Decomposition? (2020)

This article is part of the “What Is” series, available from https://nhigham.com/category/what-is and in PDF form from the GitHub repository https://github.com/higham/what-is.

Pseudo-orthogonal matrices arise in hyperbolic problems, that is, problems where there is an underlying indefinite scalar product or weight matrix. An example is the problem of downdating the Cholesky factorization, where in the simplest case we have the Cholesky factorization of a symmetric positive definite and want the Cholesky factorization of , which is assumed to be symmetric positive definite. A more general downdating problem is that we are given

and the Cholesky factorization and wish to obtain the Cholesky factor of . Note that and are . This task arises when we solve a regression problem after the observations corresponding to have been removed. The simple case above corresponds to removing one row (). Assuming that , we would like to obtain cheaply from , and numerical stability considerations dictate that we should avoid explicitly forming . If we can find a pseudo-orthogonal matrix such that

with given by (5) and upper triangular, then

so is the desired Cholesky factor.

The factorization (6) is called a hyperbolic QR factorization and it can be computed by using hyperbolic rotations to zero out the elements of . A hyperbolic rotation has the form (4), and an hyperbolic rotation is an identity matrix with a hyperbolic rotation embedded in it at the intersection of rows and columns and , for some and .

In general, a hyperbolic QR factorization of with and has the form with pseudo-orthogonal with respect to and upper triangular. The factorization exists if is positive definite.

Another hyperbolic problem is the indefinite least squares problem

where , , and are given, and with . For or we have the standard least squares (LS) problem and the quadratic form is definite, while for the problem is to minimize a genuinely indefinite quadratic form. This problem arises, for example, in the area of optimization known as smoothing.

The normal equations for (7) are , and since the Hessian matrix of the quadratic objective function in (7) is it follows that the indefinite least squares problem has a unique solution if and only if is positive definite. To solve the problem we can use a hyperbolic QR factorization to write

Solving the problem now reduces to solving the triangular system , where comprises the first components of . The same equation can also be obtained without using the normal equations by substituting the hyperbolic QR factorization into (7).

A simple technique exists for converting pseudo-orthogonal matrices into orthogonal matrices and vice versa. Let with , partition

and assume is nonsingular. The exchange operator is defined by

It is easy to see that the exchange operator is involutory, that is,

and moreover (recalling that is given by (5)) that

The next result gives a formula for the inverse of .

Lemma 1. Let with nonsingular. If is nonsingular and exists then is nonsingular and .

Proof. Consider the equation

By solving the first equation for and then eliminating from the second equation we obtain

By the same argument applied to , we have

Hence for any and there is a unique and , which implies by (10) that is nonsingular and that .

Now we will show that the exchange operator maps pseudo-orthogonal matrices to orthogonal matrices and vice versa.

Theorem 2. Let . If is pseudo-orthogonal then is orthogonal. If is orthogonal and is nonsingular then is pseudo-orthogonal.

Proof. If is pseudo-orthogonal then , which implies that is nonsingular. Since , it follows that also has a nonsingular block and so exists. Furthermore, using Lemma 1, . But (9) shows that , and we conclude that is orthogonal.

Assume now that is orthogonal with nonsingular. Then exists and Lemma 1 shows that is nonsingular and . Hence, using (9),

which shows that is pseudo-orthogonal.

This MATLAB example uses the exchange operator to convert an orthogonal matrix obtained from a Hadamard matrix into a pseudo-orthogonal matrix.

>> p = 2; n = 4;
>> A = hadamard(n)/sqrt(n), Sigma = blkdiag(eye(p),-eye(n-p))
A =
   5.0000e-01   5.0000e-01   5.0000e-01   5.0000e-01
   5.0000e-01  -5.0000e-01   5.0000e-01  -5.0000e-01
   5.0000e-01   5.0000e-01  -5.0000e-01  -5.0000e-01
   5.0000e-01  -5.0000e-01  -5.0000e-01   5.0000e-01
Sigma =
     1     0     0     0
     0     1     0     0
     0     0    -1     0
     0     0     0    -1
>> Q = exc(A,p), Q'*Sigma*Q
Q =
     1     1    -1     0
     1    -1     0    -1
     1     0    -1    -1
     0     1    -1     1
ans =
     1     0     0     0
     0     1     0     0
     0     0    -1     0
     0     0     0    -1

The code uses the function

function X = exc(A,p)
%EXC     Exchange operator.
%   EXC(A,p) is the result of applying the exchange operator to 
%   the square matrix A, which is regarded as a block 2-by-2 
%   matrix with leading block of dimension p.  
%   p defaults to floor(n)/2.

[m,n] = size(A);
if m ~= n, error('Matrix must be square.'), end
if nargin < 2, p = floor(n/2); end

A11 = A(1:p,1:p);
A12 = A(1:p,p+1:n);
A21 = A(p+1:n,1:p);
A22 = A(p+1:n,p+1:n);

X21 = A11\A12;
X = [inv(A11) -X21;
     A21/A11  A22-A21*X21];

For an orthogonal matrix expressed in block form there is a close relationship between the singular value decompositions (SVDs) of the blocks, as revealed by the CS decomposition (see What Is the CS Decomposition?). An analogous decomposition holds for a pseudo-orthogonal matrix. Let be pseudo-orthogonal with respect to in (5), and suppose that is partitioned as

Then there exist orthogonal matrices and such that

where , , and , with for all . This is the hyperbolic CS decomposition, and it can be proved by applying the CS decomposition of an orthogonal matrix to .

The leading principal submatrix in (11) generalizes the matrix (4), and in fact it can be permuted into a direct sum of such matrices.

Note that the matrix on the right in (11) is symmetric positive definite. Therefore the singular values of are the eigenvalues of that matrix, namely

Since for all , the first singular values occur in reciprocal pairs, hence the largest and smallest singular values satisfy (with strict inequality unless ). This gives another proof of (3).

While an orthogonal matrix is perfectly conditioned, a pseudo-orthogonal matrix can be arbitrarily ill conditioned, as follows from (3). For example, the MATLAB function gallery('randjorth') produces a random pseudo-orthogonal matrix with a default condition number of sqrt(1/eps).

>> rng(1); A = gallery('randjorth',2,2) % p = 2, n = 4
A =
   2.9984e+03  -4.2059e+02   1.5672e+03  -2.5907e+03
   1.9341e+03  -2.6055e+03   3.1565e+03  -7.5210e+02
   3.1441e+03  -6.2852e+02   1.8157e+03  -2.6427e+03
   1.6870e+03  -2.5633e+03   3.0204e+03  -5.4157e+02
>> cond(A)
ans =
   6.7109e+07

This means that algorithms that use pseudo-orthogonal matrices are potentially numerically unstable. Therefore algorithms need to be carefully constructed and rounding error analysis must be done to ensure that an appropriate form of numerical stability is obtained.

Pseudo-orthogonal matrices form the automorphism group of the scalar product defined by for . More results for pseudo-orthogonal matrices can be obtained as special cases of results for automorphism groups of general scalar products. See, for example, Mackey, Mackey, and Tisseur (2006).

For the set of pseudo-orthogonal matrices is known to have four connected components, a topological property that can be proved using the hyperbolic CS decomposition (Motlaghian, Armandnejad, and Hall, 2018).

One can define pseudo-unitary matrices in an analogous way, as such that . These correspond to the automorphism group of the scalar product for . The results we have discussed generalize in a straightforward way to pseudo-unitary matrices.

The exchange operator is also known as the principal pivot transform and as the sweep operator in statistics. Tsatsomeros (2000) gives a survey of its properties

The hyperbolic CS decomposition was derived by Lee (1948) and, according to Lee, was present in work of Autonne (1912).

This is a minimal set of references, which contain further useful references within.

• Adam Bojanczyk, Nicholas J. Higham and Harikrishna Patel, Solving the Indefinite Least Squares Problem by Hyperbolic QR Factorization, SIAM J. Matrix Anal. Appl. 24(3), 914–931, 2003
• Nicholas J. Higham, -Orthogonal Matrices: Properties and Generation, SIAM Rev. 45(3), 504–519, 2003
• H. C. Lee, Canonical Factorization of Pseudo-Unitary Matrices, Proc. London Math. Soc. s2-50, 230–241, 1948.
• D. Steven Mackey, Niloufer Mackey, and Françoise Tisseur, Structured Factorizations in Scalar Product Spaces, SIAM J. Matrix Anal. Appl. 27(3), 821–850, 2006.
• Sara M. Motlaghian, Ali Armandnejad, and Frank J. Hall, Topological Properties of -Orthogonal Matrices, Linear and Multilinear Algebra 66(12), 2524–2533, 2018.
• Michael Stewart and G. W. Stewart, On Hyperbolic Triangularization: Stability and Pivoting, SIAM J. Matrix Anal. Appl. 19(4), 847–860, 1998
• Michael J. Tsatsomeros, Principal Pivot Transforms: Properties and Applications, Linear Algebra Appl. 307, 151–165, 2000.

• What Is an Orthogonal Matrix? (2020)
• What Is the CS Decomposition? (2020)

This article is part of the “What Is” series, available from https://nhigham.com/category/what-is and in PDF form from the GitHub repository https://github.com/higham/what-is.

LU factorization is intimately connected with Gaussian elimination. Recall that Gaussian elimination transforms a matrix to upper triangular form in stages. At the th stage, multiples of row are added to the later rows to eliminate the elements below the diagonal in column , using the formulas

where the quantities are the multipliers. Of course each must be nonzero for these formulas to be defined, and this is connected with the conditions of Theorem 1, since . The final is the upper triangular LU factor, with for , and for , that is, the multipliers make up the factor (for a proof of these properties see any textbook on numerical linear algebra).

The matrix factorization viewpoint is well established as a powerful paradigm for thinking and computing. Separating the computation of LU factorization from its application is beneficial. For example, given we saw above how to solve . If we need to solve for another right-hand side we can just solve and , re-using the LU factorization. Similarly, solving reduces to solving the triangular systems and .

An LU factorization can be computed by directly solving for the components of and in the equation . Indeed because has unit diagonal the first row of is the same as the first row of , and then determines the first column of . One can go on to determine the th row of and the th row of , for . This leads to the Doolittle method, which involves inner products of partial rows of and partial columns of .

Given the equivalence between LU factorization and Gaussian elimination we can also employ the Gaussian elimination equations:

This ordering of the loops in the factorization is the basis of early Fortran implementations of LU factorization, such as that in LINPACK. The inner loop travels down the columns of , accessing contiguous elements of since Fortran stores arrays by column. Interchanging the two inner loops gives the ordering, which updates the matrix a row at a time, and is appropriate for a language such as C that stores arrays by row.

The and orderings correspond to the Doolittle method. The last two of the orderings are the and orderings, to which we will return later.

For with we can write

The matrix is called the Schur complement of in .

The first row and column of and have the correct forms for a unit lower triangular matrix and an upper triangular matrix, respectively. If we can find an LU factorization then

is an LU factorization of . Note that this is simply another way to express the algorithm above.

For several matrix structures it is immediate that . If we can show that the Schur complement inherits the same structure then it follows by induction that we can compute the factorization for , and so an LU factorization of exists. Classes of matrix for which and being in the class implies the Schur complement is also in the class include

• symmetric positive definite matrices,
• -matrices,
• matrices (block) diagonally dominant by rows or columns.

(The proofs of these properties are nontrivial.) Note that the matrix (1) is row diagonally dominant, as is its factor, as must be the case since its rows are contained in Schur complements.

We say that has upper bandwidth if for and lower bandwidth if for . Another use of (2) is to show that and inherit the bandwidths of .

Theorem 2. Let have lower bandwidth and upper bandwidth . If has an LU factorization then has lower bandwidth and has upper bandwidth .

Proof. In (2), the first column of and the first row of have the required structure and has upper bandwidth and lower bandwidth , since and have only and nonzero components, respectively. The result follows by induction.

In order to achieve high performance on modern computers with their hierarchical memories, LU factorization is implemented in a block form expressed in terms of matrix multiplication and the solution of multiple right-hand side triangular systems. We describe two block forms of LU factorization. First, consider a block form of (2) with block size , where is :

Here, is the Schur complement of in , given by . This leads to the following algorithm:

1. Factor .
2. Solve for .
3. Solve for .
4. Form .
5. Repeat steps 1–4 on to obtain .

The factorization on step 1 can be done by any form of LU factorization. This algorithm is known as a right-looking algorithm, since it accesses data to the right of the block being worked on (in particular, at each stage lines 2 and 4 access the last few columns of the matrix).

An alternative algorithm can derived by considering a block partitioning, in which we assume we have already computed the first block column of and :

We now compute the middle block column of and , comprising columns, by

1. Solve for .
2. Factor .
3. Solve for .
4. Repartition so that the first two block columns become a single block column and repeat steps 1–4.

This algorithm corresponds to the ordering. Note that the Schur complement is updated only a block column at a time. Because the algorithm accesses data only to the left of the block column being worked on, it is known as a left-looking algorithm.

Our description of these block algorithms emphasizes the mathematical ideas. The implementation details, especially for the left-looking algorithm, are not trivial. The optimal choice of block size will depend on the machine, but is typically in the range —.

An important point is that all these different forms of LU factorization, no matter which ordering or which value of , carry out the same operations. The only difference is the order in which the operations are performed (and the order in which the data is accessed). Even the rounding errors are the same for all versions (assuming the use of “plain vanilla” floating-point arithmetic).

Although it is most commonly used for square matrices, LU factorization is defined for rectangular matrices, too. If then the factorization has the form with lower triangular and upper trapezoidal. The conditions for existence and uniqueness of an LU factorization of are the same as those for , where .

Another form of LU factorization relaxes the structure of and from triangular to block triangular, with having identity matrices on the diagonal:

Note that is not, in general, upper triangular.

An example of a block LU factorization is

LU factorization fails on because of the zero pivot. This block LU factorization corresponds to using the leading principal submatrix of to eliminate the elements in the submatrix. In the context of a linear system , we have effectively solved for the variables and in terms of and and then substituted for and in the last two equations.

Conditions for the existence of a block LU factorization are analogous to, but less stringent than, those for LU factorization in Theorem 1.

Theorem 3. The matrix has a unique block LU factorization if and only if the first leading principal block submatrices of are nonsingular.

The conditions in Theorem 3 can be shown to be satisfied if is block diagonally dominant by rows or columns.

Note that to solve a linear system using a block LU factorization we need to solve and , but the latter system is not triangular and requires the solution of systems involving the diagonal blocks of , which would normally be done by (standard) LU factorization.

If has a unique LU factorization then for a small enough perturbation an LU factorization exists. To first order, this equation is , which gives

Since is strictly lower triangular and is upper triangular, we have, to first order,

where denotes the strictly lower triangular part and the strictly upper triangular part. Clearly, the sensitivity of the LU factors depends on the inverses of and . However, in most situations, such as when we are solving a linear system , it is the backward stability of the LU factors, not their sensitivity, that is relevant.

Since not all matrices have an LU factorization, we need the option of applying row and column interchanges to ensure that the pivots are nonzero unless the column in question is already in triangular form.

In finite precision computation it is important that computed LU factors and are numerically stable in the sense that with , where is a constant and is the unit roundoff. For certain matrix properties, such as diagonal dominance by rows or columns, numerical stability is guaranteed, but in general it is necessary to incorporate row interchanges, or row or column interchanges, in order to obtain a stable factorization.

See What Is the Growth Factor for Gaussian Elimination? for details of pivoting strategies and see Randsvd Matrices with Large Growth Factors for some recent research on growth factors.

This is a minimal set of references, which contain further useful references within.

• Jack J. Dongarra, Iain S. Duff, Danny C. Sorensen, and Henk A. Van der Vorst, Numerical Linear Algebra for High-Performance Computers, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 1998. (For different implementations of LU factorization.)
• Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms, second edition, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2002.

• Randsvd Matrices with Large Growth Factors (2020)
• What Is a Block Matrix? (2020)
• What is a Diagonally Dominant Matrix? (2021)
• What Is the Growth Factor for Gaussian Elimination? (2020)

This article is part of the “What Is” series, available from https://nhigham.com/category/what-is and in PDF form from the GitHub repository https://github.com/higham/what-is.

Theorem 1 can be used to obtain information about the location of the eigenvalues of a matrix. Indeed if is an eigenvalue of then is singular and hence cannot be strictly diagonally dominant, by Theorem 1. So cannot be true for all . Gershgorin’s theorem is simply a restatement of this fact.

Theorem 3 (Gershgorin’s theorem).

The eigenvalues of lie in the union of the discs in the complex plane

If is symmetric with positive diagonal elements and satisfies the conditions of Theorem 1 or Theorem 2 then it is positive definite. Indeed the eigenvalues are real and so in Gershgorin’s theorem the discs are intervals and , so , so the eigenvalues are nonnegative, and hence positive since nonzero. This provides another proof that the matrix in (5) is positive definite.

In some situations is not diagonally dominant but a row or column scaling of it is. For example, the matrix

is not diagonally dominant by rows or columns but

is strictly diagonally dominant by rows.

A matrix is generalized diagonally dominant by rows if is diagonally dominant by rows for some diagonal matrix with for all , that is, if

It is easy to see that if is irreducible and there is strictly inequality in (6) for some then is nonsingular by Theorem 2.

It can be shown that is generalized diagonally dominant by rows if and only if it is an -matrix, where an -matrix is a matrix for which the comparison matrix , defined by

is an -matrix (see What Is an M-Matrix?).

A matrix is block diagonally dominant by rows if, for a given norm and block partitioning , the diagonal blocks are all nonsingular and

is block diagonally dominant by columns if is block diagonally dominant by rows. If the blocks are all then block diagonal dominance reduces to the usual notion of diagonal dominance. Block diagonal dominance holds for certain block tridiagonal matrices arising in the discretization of PDEs.

Analogues of Theorems 1 and 2 giving conditions under which block diagonal dominance implies nonsingularity are given by Feingold and Varga (1962).

If a matrix is strictly diagonally dominant then we can bound its inverse in terms of the minimum amount of diagonal dominance. For full generality, we state the bound in terms of generalized diagonal dominance.

Theorem 4.

If and is strictly diagonally dominant by rows for a diagonal matrix with for all , then

where .

Proof. Assume first that . Let satisfy and let . Applying (3) gives . The result is obtained on applying this bound to and using . .

Another bound for when is strictly diagonally dominant by rows can be obtained by writing , where , , and for . It is easy to see that , which gives another proof that is nonsingular. Then

This bound implies that , so in view of its sign pattern is an -matrix, which essentially proves one direction of the -matrix equivalence in the previous section. The same bound holds if is diagonally dominant by columns, by writing .

An upper bound also holds for block diagonal dominance.

Theorem 5.

If is block diagonally dominant by rows then

where .

It is interesting to note that the inverse of a strictly row diagonally dominant matrix enjoys a form of diagonal dominance, namely that the largest element in each column is on the diagonal.

Theorem 6.

If is strictly diagonally dominant by rows then satisfies for all .

Proof. For we have . Let . Taking absolute values in gives

or , since . This inequality holds for all , so we must have , which gives the result.

Theorems 1 and 2 have a long history and have been rediscovered many times. Theorem 1 was first stated by Lévy (1881) with additional assumptions. In a short but influential paper, Taussky (1949) pointed out the recurring nature of the theorems and gave simple proofs (our proof of Theorem 2 is Taussky’s). Schneider (1977) attributes the surge in interest in matrix theory in the 1950s and 1960s to Taussky’s paper and a few others by her, Brauer, Ostrowski, and Wielandt. The history of Gershgorin’s theorem (published in 1931) is intertwined with that of Theorems 1 and 2; see Varga’s 2004 book for details.

Theorems 4 and 5 are from Varah (1975) and Theorem 6 is from Ostrowski (1952).

This is a minimal set of references, which contain further useful references within.

• David G. Feingold and Richard S. Varga, Block Diagonally Dominant Matrices and Generalizations of the Gerschgorin Circle Theorem, Pacific J. Math. 12(4), 1241–1250, 1962.
• A. M. Ostrowski, Note on Bounds for Determinants with Dominant Principal Diagonal, Proc. Amer. Math. Soc. 3, 260–30, 1952.
• Hans Schneider, Olga Taussky-Todd’s Influence on Matrix Theory and Matrix Theorists: A Discursive Personal Tribute, Linear and Multilinear Algebra 5, 197–224, 1977.
• Olga Taussky, A Recurring Theorem on Determinants, Amer. Math. Monthly 56(2), 672–676, 1949.
• J. M. Varah, A Lower Bound for the Smallest Singular Value of a Matrix, Linear Algebra Appl. 11, 3–5, 1975.
• Richard Varga, Geršgorin and His Circles, Springer-Verlag, Berlin, 2004.

• What Is a Symmetric Positive Definite Matrix? (2020)
• What Is an M-Matrix? (2021)

This article is part of the “What Is” series, available from https://nhigham.com/category/what-is and in PDF form from the GitHub repository https://github.com/higham/what-is.