Lesson Explainer: LU Decomposition: Doolittle’s Method Mathematics

In this explainer, we will learn how to find the LU decomposition (factorization) of a matrix using Doolittle’s method.

Despite all of the elegance and simplicity of the algebraic rules which govern many of the associated operations, there is no escaping the fact that linear algebra is a computationally heavy discipline and that this is especially true for those who are new to the subject. In order to properly understand and appreciate the benefits of linear algebra, it is necessary for a new student to undertake many of the common calculations manually, before these are either streamlined or replaced with more efficient methods. However, there are some operations such as matrix multiplication where there is little room for maneuver in terms of simplifying a calculation and, ultimately, there is no escape from having to perform a large number of arithmetic operations.

Due to the number of computations involved in linear algebra, we are constantly motivated to find better, less burdensome methods for completing calculations. There are few advantages in completing more calculations than is strictly necessary, especially in terms of the time taken and the possibility of introducing mistakes with each superfluous stage of the working. Even in the modern age of computing, where calculations can be completed flawlessly and autonomously at a rate of billions per second, we are still motivated to decrease the number of calculations that are needed, given that there are so few reasons not to!

In linear algebra, one of the most common tricks when working with a matrix is to see if it can be shifted or split up into a form that will be easier for subsequent calculations (especially matrix multiplication, exponentiation, and inversion). In many of these methods, it is often desirable to shift a matrix of interest into specific forms, with one of the most common being the triangular forms. In this explainer, we will show, where possible, how to take a general matrix and split it into two separate triangular matrices that can be combined together with matrix multiplication to obtain the original matrix. Although the process for completing this is actually less complicated than it often appears, it will help us to first clarify several definitions before demonstrating this method. We will begin with the notion of upper and lower-triangular matrices, which are common objects of interest in linear algebra.

We will demonstrate this concept with several examples. Notice that all of the following matrices are upper triangular, because all entries below the (highlighted) diagonal entries are zero:

As shown in the middle matrix, it is perfectly valid for any of the diagonal entries to be equal to zero. The definition only requires that the entries below the diagonal entries all have a value of zero, and there is no mention of the value of the diagonal entries themselves or of the entries that are above the diagonal entries. Any matrix is not in upper-triangular form if any of the entries below the diagonal are nonzero. Taking the matrices above, we could easily add in some errant entries that would ruin the upper-triangular form:

Lower-triangular matrices are defined in a very similar way, except we require that all entries above the diagonal entries be zero. For example, all of the following are lower-triangular matrices:

This property can be removed by adding some nonzero entries above the diagonal entries. For example,

The aim of this explainer will be to take a general matrix and to split it into its “LU decomposition,” meaning that , where L is a lower-triangular matrix and U is an upper-triangular matrix. We should begin by saying that this will not always be possible and that we may occasionally need the “PLU” decomposition, for which there is a separate explainer. We will continue undaunted by this fact. We will begin to describe the general method that we will use.

As we begin to develop our method for finding the LU decomposition of a matrix, it will be easiest for us to think in terms of elementary row operations, as this is the method that we would normally choose if we were aiming to manipulate the entries of a matrix, either as part of Gauss–Jordan elimination or as one of the many other methods that rely on row operations. More specifically, we will show how each type of row operation can be written in terms of a simple, invertible, lower-triangular matrix, hence allowing a series of row operations to be encoded exactly and in full by an accompanying series of matrices. There are three types of elementary row operations, with the simplest two of the three being the row swap operation . This is the only one of the three types of elementary row operations that we do not cover in this explainer, as it pertains to PLU decomposition rather than LU decomposition. For the other two types of row operations, we will define the equivalent elementary matrix. Please be aware that the literature on linear algebra has never reached a consensus as to how the elementary matrices should be labeled, so we are using our own terminology.

This result is simple enough to demonstrate by example. Suppose that we had the square matrix and that we wished to perform the row operation in order to obtain the row-equivalent matrix . Then, we would expect to find that

Using the definition above, we can construct the corresponding elementary matrix. This will just be a copy of the identity matrix but with the diagonal entry in the third row being equal to 5 instead of 1. The matrix is

From this, we can check that has produced the row-equivalent matrix that we were expecting. Practically speaking, the second type of row operation is not as useful as the third type and we will find ourselves making less use of this type of row operation. However, given that we may choose to use this row operation, we will still explain how to calculate the inverse of . Fortunately, the result is simple to both state and verify, so we will summarize it in the following theorem without proof.

If we had scaled a row by a constant and wanted to undo this operation, then we would just multiply the same row in the new matrix by the constant , thus returning the original matrix. This rationale is respected by the theorem above and can be demonstrated by example. We can check that the theorem is correct in the case of the example above, by noting that the inverse matrix of is the matrix . We can then apply this to the row-equivalent matrix to recover the original matrix :

We could equivalently have checked directly that , where is the identity matrix.

To complete the LU decomposition of a matrix, the most important type of row operation is the third type defined by . For example, it is possible to complete the LU decomposition of a matrix without using the second type of row operation that we described above. Generally speaking, however, it would be impossible to achieve this without the third type of row operation. This importance is easy to appreciate, given that the third type of row operation is the one which is used most commonly when completing any process that involves the elementary row operations. For example, there would generally be little hope of completing Gauss–Jordan elimination if we did not have license to use the third type of row operation. With this in mind, we now define the third type of elementary matrix.

We will demonstrate this in relation to the example above, where we had the matrix

If we wanted to apply the row operation , then we would expect the resulting row-equivalent matrix

Because of the definition above, we know that we can encode the elementary row operation by the elementary matrix . This matrix will be identical to the identity matrix , except for an entry in the second row and first column having a value of 3 instead of 0. This matrix therefore has the form

We can verify that this matrix can be used to recover the row-equivalent matrix by matrix multiplication:

With the second type of row operation, we said that it was important to understand both the matrix itself and the corresponding inverse matrix, so that we can complete the LU decomposition. This is also the case of the third type of row operation that we were describing in the definition above. Even though the third type of row operation is more complicated than the second type, the inverse operation is little more difficult to express.

Suppose that we had taken a matrix and performed the elementary row operation . If we wanted to undo this effect, then we would use the row operation on the new matrix, hence returning the original matrix. This understanding is totally encapsulated by the statement of the theorem above, for which we now provide a check. To test this, we can produce the inverse matrix corresponding to the example above. Since we used , we have that the inverse matrix is . We can then check that

Now that we have discussed the necessary ingredients for LU decomposition, we will both define this concept and then show how it can be achieved. We will do so initially by way of example, wherein we will briefly justify all of the steps that we are taking. Many of the techniques are ostensibly similar to those used in Gauss–Jordan elimination and other related methods. The completion of the LU decomposition focuses on achieving an upper-triangular form, whereas Gauss–Jordan elimination focuses on achieving an echelon form. These two concepts are similar but are not identical, and in fact when attempting to complete the LU decomposition we will have fewer options than we do generally for Gauss–Jordan.

We will demonstrate this with the simplest possible example of a matrix. We take the matrix and task ourselves with finding a lower-triangular matrix L and an upper-triangular matrix U such that . Suppose that we begin with the familiar method of first highlighting all of the pivot entries of :

After having practiced Gauss–Jordan elimination, our first instinct might be to simplify , which in this case would involve removing the factor of 3 that appears in every entry in the top row. To do this, we would use the elementary row operation , which has the corresponding elementary matrix

We will also state the inverse of this matrix:

Given that these two matrices are inverses of each other, they must obey the relationship where is the identity matrix. This is certainly the case, which means that we can write this equation in full:

It is this statement that allows us to begin finding the LU decomposition of . Given that for any matrix of order , we can write

Now the reason for writing this equation is explained by equation (1), which we duly substitute into the equation immediately above, giving

Matrix multiplication is associative, which means that we can combine the second and third matrices using the operation of matrix multiplication. This gives precisely the desired effect, as the new middle matrix is the same as the original matrix but with the constant factor of 3 now removed from the first row:

The leftmost matrix is a lower-triangular matrix, as all of the entries above diagonal entries are zero. In the case of a matrix, this leaves only the entry in the top right, which is shown to be zero. Unfortunately, the middle matrix is not yet in upper-triangular form and so we have not yet found our LU decomposition. Now, once again inspecting the pivot entries shows that we can remove the pivot in the second row to turn the middle matrix into an upper-triangular matrix. This is completed with the row operation . The corresponding elementary matrix is which has the inverse elementary matrix

Once again, these two elementary matrices are multiplicative inverses of each other, so we have the relationship . Written out in long form, this is

We can take equation (2) and insert the identity matrix as shown, without changing the validity of the equation:

Substituting equation (3) then gives

We once again rely on the property of matrix multiplication being associative. Given that this is always the case, we can order the calculations in any way we would prefer. In this situation, it is best to order the calculations as shown:

Then, completing the two matrix multiplications gives the required decomposition:

We have therefore obtained the two matrices

such that . It is easily checked that L is a lower-triangular matrix and that U is an upper-triangular matrix and the final confirmation of the result can be completed by checking that equation (4) is true.

The calculation of an LU decomposition, once practiced, would normally be completed without much of the supporting explanation that was given in the above example. Much as elementary row operations can be used to manipulate a matrix into any type of possible decomposition, obtaining an LU decomposition requires that we only ever produce lower-triangular matrices on the leftmost side of our equations. This means that our elementary matrices (and hence our elementary row operations) are chosen accordingly, as it will be counterproductive to use elementary row operations that incur a non-lower-triangular form in the leftmost terms. This is not as challenging as it might seem at first, and the criteria are essentially those described by Doolittle’s algorithm, which we will discuss in more detail later in this explainer.

The main point of understanding before the statement of Doolittle’s algorithm is that, for each elementary symmetric matrix, we are also calculating the corresponding inverse of this matrix. When a matrix is multiplied by its inverse, the result is the identity matrix. Also, when an identity matrix is multiplied by any other matrix, the result is this matrix. Doolittle’s algorithm requires both of these properties as it is normally phrased in a way that is more technical than the way we present the result below.

There are two key points to this theorem. Firstly, the matrices are lower-triangular, which means that the inverses will also be lower-triangular and so will the product of these matrices. This guarantees that the matrix is in lower-triangular form, as needed. The second key point is that this decomposition, if possible, does actually achieve the matrix . This is a simple matter of checking that

The most useful property in this situation is that matrix multiplication is associative, meaning that we can bracket the terms in any order that we would find to be most helpful. In this case, we will pair together the matrices, working from the central point outward. Given that is the inverse of , we would find that , the identity matrix. This allows equation (9) to be written in the following way: where the final line has been achieved by recalling that for any matrix . It is now clear that because every elementary matrix will become paired with its own, which will then return the identity matrix, the next step follows exactly the same pattern:

Continuing to pair up these matrices until the final pair gives

It is now clear why Doolittle’s algorithm works. The only reason why this approach would not work would be if one of the elementary matrices was not lower triangular. If it is not possible to create the LU decomposition without the use of non-lower-triangular matrices, then a minor adjustment is needed, which is referred to as the PLU decomposition, wherein non-lower-triangular matrices of a certain type are temporarily permitted in the calculations. For the remainder of this explainer, we will give examples where it is possible to find the LU decomposition, with the PLU decomposition being treated in a separate explainer.

As ever with linear algebra, once a technique has been mastered for a matrix of a particular order, it is normally little extra conceptual effort to apply the same approach to matrices with different orders. This will again show to be the case with the application of Doolittle’s algorithm to a matrix with order . Much as the calculations for a matrix are more numerous, the overall method is no more difficult than the previous examples, relying on exactly the same principles to decompose the matrix into LU form.

The advantages of the LU decomposition are both conceptual and computational. If a matrix can be written in the form , then we are essentially able to think of this matrix as two different matrices, both of which are triangular. Triangular matrices are often easier to work with, especially when trying to understand the determinant of a square matrix. If a square matrix is in either of the two triangular forms, then the determinant of this matrix is the product of the diagonal entries. Suppose that we had the a square matrix with order , having the form

Ordinarily, finding the determinant of might be a computationally expensive process, which we would rather avoid if an alternative method was available. Let us further suppose that it is possible to achieve an LU decomposition for , meaning that we can write , where L and U are, respectively, the lower and triangular matrices:

Due to the fact that L and U are triangular matrices, we can calculate their determinants as the products of the diagonal entries. Therefore,

One of the key properties of the determinant is that it is multiplicative, meaning that if then . The purpose of the above representations is now revealed: we can use the multiplicative property of the determinant on the LU decomposition , giving the final result

Therefore, if the LU decomposition of is either known or easily obtainable, the determinant can be calculated in a way that is essentially trivial.