In linear algebra courses we learn that the solution to a linear system of equations in unknowns can be written , where is the matrix inverse. What is not always emphasized is that there are very few circumstances in which one should compute . Indeed one would not solve the scalar () system by computing , but rather would carry out a division . In the case, it is faster and more accurate to solve a linear system by LU factorization (Gaussian elimination) with partial pivoting than by inverting (which has, in any case, to be done by LU factorization).

Rare cases where is required are in statistics, where the diagonal elements of the inverse of the covariance matrix are relevant quantities, and in certain algorithms for computing matrix functions.

The solution to the linear least squares problem , where is a full-rank matrix with , satisfies the normal equations . It is therefore natural to form the symmetric positive definite matrix and solve the normal equations by Cholesky factorization. While fast, this method is numerically unstable when is ill conditioned. By contrast, solving the least squares problem via QR factorization is always numerically stable.

What is wrong with the cross-product matrix (also known as the Gram matrix)? It squares the data, which can cause a loss of information in floating-point arithmetic. For example, if

where is the unit roundoff of the floating point arithmetic, then

is positive definite but, since , in floating-point arithmetic rounds to and so

which is singular, and the information in has been lost.

Another problem with the cross product matrix is that the -norm condition number of is the square of that of , and this leads to numerical instability in algorithms that work with when the condition number is large.

The cost of evaluating a matrix product depends on the order in which the product is evaluated (assuming the matrices are not all ). More precisely, matrix multiplication is associative, so , and in general the cost of the evaluation of a product depends on where one puts the parentheses. One order may be much superior to others, so one should not simply evaluate the product in a fixed left-right or right-left order. For example, if , , and are -vectors then can be evaluated as

• : a vector outer product followed by a matrix–vector product, costing operations, or
• : a vector scalar product followed by a vector scaling, costing just operations.

In general. finding where to put the parentheses in a matrix product in order to minimize the operation count is a difficult problem, but for many cases that arise in practice it is easy to determine a good order.

Symmetric positive definite matrices (symmetric matrices with positive eigenvalues) are ubiquitous, not least because they arise in the solution of many minimization problems. However, a matrix that is supposed to be positive definite may fail to be so for a variety of reasons. Missing or inconsistent data in forming a covariance matrix or a correlation matrix can cause a loss of definiteness, and rounding errors can cause a tiny positive eigenvalue to go negative.

Definiteness implies that

• the diagonal entries are positive,
• ,
• for all ,

but none of these conditions, or even all taken together, guarantees that the matrix has positive eigenvalues.

The best way to check definiteness is to compute a Cholesky factorization, which is often needed anyway. The MATLAB function chol returns an error message if the factorization fails, and a second output argument can be requested, which is set to the number of the stage on which the factorization failed, or to zero if the factorization succeeded. In the case of failure, the partially computed factor is returned in the first argument, and it can be used to compute a direction of negative curvature (as needed in optimization), for example.

This sin takes the top spot in Schmelzer and Hauser’s Seven Sins in Portfolio Optimization, because in portfolio optimization a negative eigenvalue in the covariance matrix can identify a portfolio with negative variance, promising an arbitrarily large investment with no risk!

One of the fundamental tenets of numerical linear algebra is that one should try to exploit any matrix structure that might be present. Sparsity (a matrix having a large number of zeros) is particularly important to exploit, since algorithms intended for dense matrices may be impractical for sparse matrices because of extensive fill-in (zeros becoming nonzero). Here are two examples of structures that can be exploited.

Matrices from saddle point problems are symmetric indefinite and of the form

with symmetric positive definite. Much work has been done on developing numerical methods for solving that exploit the block structure and possible sparsity in and . A second example is a circulant matrix

Circulant matrices have the important property that they are diagonalized by a unitary matrix called the discrete Fourier transform matrix. Using this property one can solve in operations, rather than the operations required if the circulant structure is ignored.

Ideally, linear algebra software would detect structure in a matrix and call an algorithm that exploits that structure. A notable example of such a meta-algorithm is the MATLAB backslash function x = A\b for solving . Backslash checks whether the matrix is triangular (or a permutation of a triangular matrix), upper Hessenberg, symmetric, or symmetric positive definite, and applies an appropriate method. It also allows to be rectangular and solves the least squares problem if there are more rows than columns and the underdetermined system if there are more columns than rows.

An matrix is nonsingular if and only if its determinant is nonzero. One might therefore expect that a small value for indicates a matrix that is nearly singular. However, the size of tells us nothing about near singularity. Indeed, since we can achieve any value for the determinant by multiplying by a scalar , yet is no more or less nearly singular than for .

Another limitation of the determinant is shown by the two matrices

Both matrices have unit diagonal and off-diagonal elements bounded in modulus by . So , yet

So is ill conditioned for large . In fact, if we change the element of to then the matrix becomes singular! By contrast, is always very well conditioned. The determinant cannot distinguish between the ill-conditioned and the well-conditioned .

For any matrix and any consistent matrix norm it is true that for all , where the are the eigenvalue of . Since the eigenvalues of are , it follows that the matrix condition number is bounded below by the ratio of largest to smallest eigenvalue in absolute value, that is,

But as the matrix in (1) shows, this bound can be very weak.

It is singular values not eigenvalues that characterize the condition number for the 2-norm. Specifically,

where is a singular value decomposition (SVD), with and orthogonal and , . If is symmetric, for example, then the sets and are the same, but in general the eigenvalues and singular values can be very different.

• What Is a Cholesky Factorization? (2020)
• What Is a Circulant Matrix? (2022)
• What Is a Condition Number? (2020)
• What is a Sparse Matrix? (2020)
• What Is a Symmetric Positive Definite Matrix? (2020)
• What Is the Singular Value Decomposition? (2020)

Many algorithms are implemented in blocked form. For example, an inner product of two -vectors and can computed as

where and is the block size. The inner product has been broken into smaller inner products of size , which are computed independently then summed. Many linear algebra algorithms are blocked in an analogous way, where the blocking is into submatrices with rows or columns (or both). Careful analysis of the error analysis shows that a blocked algorithm has an error bound about a factor of smaller than that for the corresponding unblocked algorithm. Practical block sizes for matrix algorithms are typically or , so blocking brings a substantial reduction in the error bounds.

In fact, one can do even better than an error bound of order . By computing the sum with a more accurate summation method the error constant is further reduced to (this is the FABsum method of Blanchard et al. (2020)).

Intel x86 processors support an 80-bit extended precision format with a 64-bit significand, which is compatible with that specified in the IEEE standard. When a compiler uses this format with 80-bit registers to accumulate sums and inner products it is effectively working with a unit roundoff of rather than for double precision, giving error bounds smaller by a factor up to .

Some processors have a fused multiply–add (FMA) operation, which computes a combined multiplication and addition with one rounding error instead of two. This results in a reduction in error bounds by a factor .

Mixed precision block FMA operations , with matrices of fixed size, are available on Google tensor processing units, NVIDIA GPUs, and in the ARMv8-A architecture. For half precision inputs these devices can produce results of single precision quality, which can give a significant boost in accuracy when block FMAs are chained together to form a matrix product of arbitrary dimension.

Worst-case rounding error bounds suffer from the problem that they are not attainable for most specific sets of data and are unlikely to be nearly attained. Stewart (1990) noted that

To be realistic, we must prune away the unlikely. What is left is necessarily a probabilistic statement.

Theo Mary and I have recently developed probabilistic rounding error analysis, which makes probabilistic assumptions on the rounding errors and derives bounds that hold with a certain probability. The key feature of the bounds is that they are proportional to when a corresponding worst-case bound is proportional to . In the most general form of the analysis (Connolly, Higham, and Mary, 2021), the rounding errors are assumed to be mean independent and of mean zero, where mean independence is a weaker assumption than independence.

The different features we have described can be combined to obtain significantly smaller error bounds. If we use a blocked algorithm with block size then in an inner product the standard error bound of order reduces to a probabilistic bound of order , which is a significant reduction. Block FMAs and extended precision registers provide further reductions.

For example, for a linear system of order solved in single precision with a block size of , the probabilistic error bound is of order versus for the standard worst-case bound. If FABsum is used then the bound is further reduced.

Our conclusion is that we can successfully solve linear algebra problems of greater size and at lower precisions than the standard rounding error analysis suggests. A priori bounds will always be pessimistic, though. One should compute a posteriori residuals or backward errors (depending on the problem) in order to assess the quality of a numerical solution.

For full details of the work summarized here, see Higham (2021).

• Pierre Blanchard, Nicholas Higham, and Theo Mary, A Class of Fast and Accurate Summation Algorithms, SIAM J. Sci. Comput. 42(3), A1541–A1557, 2020.
• Pierre Blanchard, Nicholas Higham, Florent Lopez, Theo Mary, and Srikara Pranesh, Mixed Precision Block Fused Multiply-Add: Error Analysis and Application to GPU Tensor Cores, SIAM J. Sci. Comput. 42(3), C124-C141, 2020
• Michael P. Connolly, Nicholas J. Higham, and Theo Mary, Stochastic Rounding and Its Probabilistic Backward Error Analysis, SIAM J. Sci. Comput. 43(1), A566–A585, 2021.
• Nicholas J. Higham, Numerical Stability of Algorithms at Extreme Scale and Low Precisions, MIMS EPrint 2021.14, Manchester Institute for Mathematical Sciences, The University of Manchester, UK, September 2021.
• Nicholas Higham and Theo Mary, A New Approach to Probabilistic Rounding Error Analysis, SIAM J. Sci. Comput. 41(5), A2815–A2835, 2019.
• G. W. Stewart, Stochastic Perturbation Theory, SIAM Rev. 32(4), 579–610, 1990.

From the inequality , for any matrix norm, where is the spectral radius (the largest magnitude of any eigenvalue of ) we have

Fir the -norm, this bound is an equality for a normal matrix (one for which ), but it can be arbitrarily weak for nonnormal matrices.

Guggenheimer, Edelman, and Johnson (1995) obtain the bound

The proof of the bound applies the arithmetic–geometric mean inequality to the numbers , where the are the singular values of . This bound can be arbitrarily weak but it is an approximate equality when are of similar order of magnitude.

Merikoski, Urpala, Virtanen, Tam, and Uhlig (1997) obtain the bound

Their proof uses a more refined application of the arithmetic–geometric mean inequality, and they show that this bound is the smallest that can be obtained based on , , and only. Hence (3) is no larger than (2), and they show that it can be smaller by no more than 1.5. Equality holds in (3) if and only if .

As an example, for three random matrices with , generated by gallery('randsvd') with three different singular value dsitributions:

Mode	(2)	(3)
One large singular value	9.88e+07	9.88e+07
One small singular value	1.21e+01	1.20e+01
Geometrically distributed singular values	5.71e+04	5.71e+04

We note that for larger the formula (3) is prone to overflow, which can be avoided by evaluating it in higher precision arithmetic.

Merikoski et al. (1997) also give a version of (3) for Hermitian positive definite :

This is the smallest bound that can be obtained based on , , and only. Equality holds in (4) if and only if the eigenvalues of satisfy . We can rewrite this upper bound as

which gives the weaker bound

This bound is analogous to (2) and is up to a factor larger than (4), this factor being attained for .

If then (4) reduces to

These bounds hold for any positive definite matrix with unit diagonal, that is, any nonsingular correlation matrix.

We can sometimes get a sharper bound than (4) and (5) by writing , where and (thus is a correlation matrix), using

and bounding using (6). For example, for the Pascal matrix

the condition number is . The bounds from (4) and (5) are both , whereas combining (4) and (7) gives a bound of .

Many other condition number bounds are available in the literature. All have their pros and cons and any bound based on limited information such as traces of powers of and the determinant will be potentially very weak.

A drawback of the bounds (3)–(6) is that they require . Sometimes the determinant is easily computable, as for a Vandermonde matrix, or can be bounded: for example, for a matrix with integer entries. If a Cholesky, LU, or QR factorization of is available then is easily computable, but in this case a good order of magnitude estimate of the condition number can be cheaply computed using condition estimation techniques (Higham, 2002, Chapter 15).

The bounds (3) and (4) are used by Higham and Lettington (2021) in investigating the most ill conditioned symmetric matrices with integer elements bounded by ; see What Is the Wilson Matrix?

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

• Heinrich Guggenheimer, Alan Edelman, and Charles Johnson, A Simple Estimate of the Condition Number of a Linear System, College Math. J. 26, 2–5, 1995.
• Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms, second edition, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2002.
• Nicholas J. Higham and Matthew C. Lettington, Optimizing and Factorizing the Wilson Matrix, to appear in Amer. Math. Monthly, 2021.
• Jorma Kaarlo Merikoski, Uoti Urpala, Ari Virtanen, Tin-Yau Tam, and Frank Uhlig, A Best Upper Bound for the 2-Norm Condition Number of a Matrix, Linear Algebra Appl. 254, 355–365, 1997.

• Roger A. Horn and Charles R. Johnson, Matrix Analysis, second edition, Cambridge University Press, 2013. My review of the second edition.
• Ilse Ipsen, Relative Perturbation Results for Matrix Eigenvalues and Singular Values, Acta Numerica 7, 151–201, 1998.

• Eigenvalue Inequalities for Hermitian Matrices (2021)
• What Is the Singular Value Decomposition? (2020)

First, we consider a general matrix and let be an eigenvalue with (the spectral radius) and a corresponding eigenvector. With , where is the vector of ones, , so

which implies since . Hence

Let be a triangular matrix. Applying the latter bound to , whose eigenvalues are its diagonal entries , gives

Combining this bound with the analogous bound for gives

We note that commonly used norms satisfy , which yields another proof of (3) and (4).

For any and such that we have the lower bound . We can choose and then solve the triangular system for to obtain the lower bound. Condition number estimation techniques, which we will describe in another article, provide ways to choose that usually yield estimates of correct to within an order of magnitude.

For the -norm, we can choose and then compute by repeated triangular solves, obtaining the lower bound . This bound is simply the power method applied to .

Let be an upper triangular matrix. The upper bounds for that we will discuss depend only on the absolute values of the elements of . This limits the ability of the bounds to distinguish between well-conditioned and ill-conditioned matrices. For example, consider

The bounds for and will be the same, yet the inverses are of different sizes (the more so as the dimension increases).

Let and write

where is strictly upper triangular and hence nilpotent with . Then

Taking absolute values and using the triangle inequality gives

where the inequalities hold elementwise.

The comparison matrix associated with a general is the matrix with

It is not hard to see that is upper triangular with and so the bound (5) is

If we replace every element above the diagonal of by the most negative off-diagonal element in its row we obtain the upper triangular matrix with

Then , where , so

Finally, let , where is strictly upper triangular with every element above the diagonal equal to the maximum element of , that is,

Then

We note that , , and are all nonsingular -matrices. We summarize the bounds.

Theorem 1.

If is a nonsingular upper triangular matrix then

We make two remarks.

• The bounds (6) are equally valid for lower triangular matrices as long as the maxima in the definitions of and are taken over columns instead of rows.
• We could equally well have written . The comparison matrix is unchanged, and (6) continues to hold as long as the maxima in the definitions of and are taken over columns rather than rows.

It follows from the theorem that

for the 1-, 2-, and -norms and the Frobenius norm. Now , , and all have nonnegative inverses, and for a matrix with nonnegative inverse we have . Hence

where the big-Oh expressions show the asymptotic cost in flops of evaluating each term by solving the relevant triangular system. As the bounds become less expensive to compute they become weaker. The quantity can be explicitly evaluated for , using . It has the same value for , and since we have

This bound is an equality for for the matrix in (1).

For the Frobenius norm, evaluating , and using , gives

For the -norm, either of (7) and (8) can be the smaller bound depending on .

For the special case of a bidiagonal matrix it is easy to show that , and so can be computed exactly in flops.

These upper bounds can be arbitrarily weak, even for a fixed , as shown by the example

for which

As , . On the other hand, the overestimation is bounded as a function of for triangular matrices resulting from certain pivoting strategies.

Theorem 1.

Suppose the upper triangular matrix satisfies

Then, for the -, -, and -norms,

Proof. The first four inequalities are a combination of (3) and (6). The fifth inequality is obtained from the expression (7) for with .

Condition (9) is satisfied for the triangular factors from QR factorization with column pivoting and for the transpose of the unit lower triangular factors from LU factorization with any form of pivoting.

The upper bounds we have described have been derived independently by several authors, as explained by Higham (2002).

• Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms, second edition, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2002. (Chapter 8.)

The function is the quadratic form for evaluated on the unit sphere, since . As is Hermitian it has a spectral decomposition , where is unitary and . Then

from which is it clear that

with equality when is an eigenvector corresponding to and , respectively, This characterization of the extremal eigenvalues of as the extrema of is due to Lord Rayleigh (John William Strutt), and is called a Rayleigh quotient. The intermediate eigenvalues correspond to saddle points of .

The Courant–Fischer theorem (1905) states that every eigenvalue of a Hermitian matrix is the solution of both a min-max problem and a max-min problem over suitable subspaces of .

Theorem (Courant–Fischer).

For a Hermitian ,

Note that the equalities are special cases of these characterizations.

In general there is no useful formula for the eigenvalues of a sum of Hermitian matrices. However, the Courant–Fischer theorem yields the upper and lower bounds

from which it follows that

This inequality shows that the eigenvalues of a Hermitian matrix are well conditioned under perturbation. We can rewrite the inequality in the symmetric form

If is positive semidefinite then (1) gives

while if is positive definite then strict inequality holds for all . These bounds are known as the Weyl monotonicity theorem.

Weyl’s inequalities (1912) bound the eigenvalues of in terms of those of and .

Theorem (Weyl).

For Hermitian and ,

The Weyl inequalities yield much information about the effect of low rank perturbations. Consider a positive semidefinite rank- perturbation . Inequality (3) with gives

(which also follows from (1)). Inequality (3) with , combined with (2), gives

These inequalities confine each eigenvalue of to the interval between two adjacent eigenvalues of ; the eigenvalues of are said to interlace those of . The following figure illustrates the case , showing a possible configuration of the eigenvalues of and of .

A specific example, in MATLAB, is

>> n = 4; eig_orig = 5:5+n-1
>> D = diag(eig_orig); eig_pert = eig(D + ones(n))'
eig_orig =
     5     6     7     8
eig_pert =
   5.2961e+00   6.3923e+00   7.5077e+00   1.0804e+01

Since and the trace is the sum of the eigenvalues, we can write

where the are nonnegative and sum to . If we greatly increase , the norm of the perturbation, then most of the increase in the eigenvalues is concentrated in the largest, since (5) bounds how much the smaller eigenvalues can change:

>> eig_pert = eig(D + 100*ones(n))'
eig_pert =
   5.3810e+00   6.4989e+00   7.6170e+00   4.0650e+02

More generally, if has positive eigenvalues and negative eigenvalues then (3) with gives

while (4) with gives

So the inertia of (the number of negative, zero, and positive eigenvalues) determines how far the eigenvalues can move as measured relative to the indexes of the eigenvalues of .

An important implication of the last two inequalities is for the case , for which we have

Exactly eigenvalues appear in one of these inequalities and appear in both. Therefore of the eigenvalues are equal to and so only eigenvalues can differ from . So perturbing the identity matrix by a Hermitian matrix of rank changes at most of the eigenvalues. (In fact, it changes exactly eigenvalues, as can be seen from a spectral decomposition.)

Finally, if has rank then and and so taking in (3) and in (4) gives

The Cauchy interlace theorem relates the eigenvalues of successive leading principal submatrices of a Hermitian matrix. We denote the leading principal submatrix of of order by .

Theorem (Cauchy).

For a Hermitian ,

The theorem says that the eigenvalues of interlace those of for all . Two immediate implications are that (a) if is Hermitian positive definite then so are all its leading principal submatrices and (b) appending a row and a column to a Hermitian matrix does not decrease the largest eigenvalue or increase the smallest eigenvalue.

Since eigenvalues are unchanged under symmetric permutations of the matrix, the theorem can be reformulated to say that the eigenvalues of any principal submatrix of order interlace those of . A generalization to principal submatrices of order is given in the next result.

Theorem.

If is a principal submatrix of order of a Hermitian then

It follows by taking to be a unit vector in the formula that for all . And of course the trace of is the sum of the eigenvalues: . These relations are the first and last in a sequence of inequalities relating sums of eigenvalues to sums of diagonal elements obtained by Schur in 1923.

Theorem (Schur).

For a Hermitian ,

where is the set of diagonal elements of arranged in decreasing order: .

These inequalities say that the vector of eigenvalues majorizes the ordered vector of diagonal elements.

An interesting special case is a correlation matrix, a symmetric positive semidefinite matrix with unit diagonal, for which the inequalities are

and . Here is an illustration in MATLAB.

>> n = 5; rng(1); A = gallery('randcorr',n);
>> e = sort(eig(A)','descend'), partial_sums = cumsum(e)
e =
  2.2701e+00   1.3142e+00   9.5280e-01   4.6250e-01   3.6045e-04
partial_sums =
  2.2701e+00   3.5843e+00   4.5371e+00   4.9996e+00   5.0000e+00

Ky Fan (1949) proved a majorization relation between the eigenvalues of , , and :

For , the inequality is the same as the upper bound of (1), and for it is an equality: .

For a Hermitian and a nonsingular , the transformation is a congruence transformation. Sylvester’s law of inertia says that congruence transformations preserve the inertia. A result of Ostrowski (1959) goes further by providing bounds on the ratios of the eigenvalues of the original and transformed matrices.

Theorem (Ostrowski).

For a Hermitian and ,

where .

If is unitary then and so Ostrowski’s theorem reduces to the fact that a congruence with a unitary matrix is a similarity transformation and so preserves eigenvalues. The theorem shows that the further is from being unitary the greater the potential change in the eigenvalues.

Ostrowski’s theorem can be generalized to the situation where is rectangular (Higham and Cheng, 1998).

The results we have described are strongly interrelated. For example, the Courant–Fischer theorem and the Cauchy interlacing theorem can be derived from each other, and Ostrowski’s theorem can be proved using the Courant–Fischer Theorem.

• Rajendra Bhatia, Linear Algebra to Quantum Cohomology: The Story of Alfred Horn’s Inequalities, Amer. Math. Monthly 108(4), 289–318, 2001.
• Roger A. Horn and Charles R. Johnson, Matrix Analysis, second edition, Cambridge University Press, 2013. My review of the second edition.
• Nicholas J. Higham and Sheung Hun Cheng, Modifying the Inertia of Matrices Arising in Optimization, Linear Algebra Appl. 275–276, 261-279, 1998.
• Beresford Parlett, The Symmetric Eigenvalue Problem, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA 1998.

The rank-nullity theorem says that

where is the null space of .

The rank cannot exceed the number of columns, or, by (5) below, the number of rows:

For any and of the same dimension,

The upper bound follows from the fact that the dimension of the sum of two subspaces cannot exceed the sum of the dimensions of the subspaces. Interestingly, the upper bound is also a corollary of the bound (3) for the rank of a matrix product, because

For the lower bound, writing and applying the upper bound gives , and likewise with the roles of and interchanged.

For any ,

Indeed implies , and implies , which implies . Hence the null spaces of and are the same. The equality (2) follows from the rank-nullity theorem.

For any and for which the product is defined,

If then , so the columns of are linear combinations of those of and so cannot have more linearly independent columns than , that is, . Using (5) below, we then have

The latter inequality can be proved without using (5) (our proof of which uses (3)), as follows. Suppose . Let the columns of span , so that has columns and for some matrix with columns. Now by the first part, so for some nonzero . But then , which contradicts the linear independence of the columns of , so we must have .

We have

We note that and are both nonsingular matrices by (2), so their product has rank . Using (3),

and hence .

Another important relation is

This is a consequence of the equality for nonsingular and .

By (2) and (3) we have . Interchanging the roles of and gives and so

In other words, the rank of is equal to the maximum number of linearly independent rows as well as the maximum number of linearly independent columns.

has rank if and only if for some and , both of rank , and this is called a full-rank factorization. The existence of such a factorization implies that by (4). Conversely, suppose that has rank . Let the columns of form a basis for the range space of . Then there are -vectors such that , , and with we have . Finally, by (3), and since we have .

A characterization of rank that is sometimes used as the definition is that it is the size of the largest nonsingular square submatrix. Equivalently, the rank is the size of the largest nonzero minor, where a minor of size is the determinant of a submatrix.

Although and have some properties in common when both products are defined (notably they have the same nonzero eigenvalues), is not always equal to . A simple example is and with and orthogonal vectors: but . An example with square and is

Note that and , where has in the th position and zeros everywhere else. Such matrices are easy to manipulate in this form (e.g., ) and are useful for constructing examples.

If we have a full-rank factorization of then we can read off the rank from the dimensions of the factors. But finding a full-rank factorization is a nontrivial task. The ultimate full-rank factorization is the SVD

where and are orthogonal, , where , and . The rank of is , the number of nonzero singular values.

In floating-point arithmetic, the standard algorithms for computing the SVD are numerically stable, that is, the computed singular values are the exact singular values of a matrix with , where is a constant and is the unit roundoff. Unfortunately, will typically be full rank when is rank deficient. For example, consider this computation.

>> n = 4; A = zeros(n); A(:) = 1:n^2, svd(A)
A =
     1     5     9    13
     2     6    10    14
     3     7    11    15
     4     8    12    16
ans =
   3.8623e+01
   2.0713e+00
   1.5326e-15
   1.3459e-16

The matrix has rank and the two zero singular values are approximated by computed singular values of order . In general, we have no way to know whether tiny computed singular values signify exactly zero singular values. In practice, one typically defines a numerical rank based on a threshold and regards computed singular values less than the threshold as zero. Indeed the MATLAB rank function computes the rank as the number of singular values exceeding , where is the largest computed singular value. If the data from which the matrix is constructed is uncertain then the definition of numerical rank should take into account the level of uncertainty in the data. Dealing with rank deficiency in the presence of data errors and in finite precision arithmetic is a tricky business.

An excellent reference for further rank relations is Horn and Johnson. Stewart describes some of the issues associated with rank-deficient matrices in practical computation.

• Roger A. Horn and Charles R. Johnson, Matrix Analysis, second edition, Cambridge University Press, 2013. My review of the second edition.
• G. W. Stewart, Rank Degeneracy, SIAM J. Sci. Statist. Comput. 5 (2), 403–413, 1984

We take the Fiedler matrix

>> A = gallery('fiedler',5)
A =
     0     1     2     3     4
     1     0     1     2     3
     2     1     0     1     2
     3     2     1     0     1
     4     3     2     1     0

The smallest eigenvalue is , so is positive semidefinite for . The Gill–Murray–Wright method gives with diagonal elements

2.0000e+01   4.6159e+00   1.3194e+00   2.5327e+00   1.0600e+01

and has while the Schnabel–Eskow method gives with diagonal elements

6     6     6     6     6

and has . If we increase the diagonal elements of by 0.5 to give comparable smallest eigenvalues for the perturbed matrices then we have

Shift	5.2361	26.180
Gill-Murray-Wright	20.000	39.068
Schnabel–Eskow	6.000	30.000

• Roger Fletcher, Semi-Definite Matrix Constraints in Optimization, SIAM J. Control Optim. 23 (4), 493–513, 1985.
• Philip Gill, Walter Murray, and Margaret Wright, Practical Optimization, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2020. Republication of book first published by Academic Press, 1981.
• Robert Schnabel and Elizabeth Eskow, A Revised Modified Cholesky Factorization Algorithm, SIAM J. Optim. 9(4), 1135–1148, 1999.