Section TD Triangular Decomposition

Our next decomposition will break a square matrix into a product of two matrices, one lower triangular and the other upper triangular. So we will write A = LU, and hence many refer to this as LU decomposition. We will see that this decomposition is very easy to compute and that it has a direct application to solving systems of equations. Since this section is about triangular matrices you might want to review the definitions and a couple of basic theorems back in Subsection OD.TM.

Subsection TD: Triangular Decomposition

With a slight condition on the nonsingularity of certain submatrices, we can split a matrix into a product of two triangular matrices.

Theorem TD
Triangular Decomposition
Suppose A is a square matrix of size n. Let {A}_{k} be the k × k matrix formed from A by taking the first k rows and the first k columns. Suppose that {A}_{k} is nonsingular for all 1 ≤ k ≤ n. Then there is a lower triangular matrix L with all of its diagonal entries equal to 1 and an upper triangular matrix U such that A = LU. Furthermore, this decomposition is unique. □

Proof We will row reduce A to a row-equivalent upper triangular matrix through a series of row operations, forming intermediate matrices {A}_{j}^{′}, 1 ≤ j ≤ n, that denote the state of the conversion after working on column j. First, the lone entry of {A}_{1} is {\left [A\right ]}_{11} and this scalar must be nonzero if {A}_{1} is nonsingular (Theorem SMZD). We can use row operations Definition RO of the form α{R}_{1} + {R}_{k}, 2 ≤ k ≤ n, where α = −{\left [A\right ]}_{1k}∕{\left [A\right ]}_{11} to place zeros in the first column below the diagonal. The first two rows and columns of {A}_{1}^{′} are a 2 × 2 upper triangular matrix whose determinant is equal to the determinant of {A}_{2}, since the matrices are row-equivalent through a sequence of row operations strictly of the third type (Theorem DRCMA). As such the diagonal entries of this 2 × 2 submatrix of {A}_{1}^{′} are nonzero. We can employ this nonzero diagonal element with row operations of the form α{R}_{2} + {R}_{k}, 3 ≤ k ≤ n to place zeros below the diagonal in the second column. We can continue this process, column by column. The key observations are that our hypothesis on the nonsingularity of the {A}_{k} will guarantee a nonzero diagonal entry for each column when we need it, that the row operations employed are always of the third type using a multiple of a row to transform another row with a greater row index, and that the final result will be a nonsingular upper triangular matrix. This is the desired matrix U.

Each row operation described in the previous paragraph can be accomplished with matrix multiplication by the appropriate elementary matrix (Theorem EMDRO). Since every row operation employed is adding a multiple of a row to a subsequent row these elementary matrices are of the form {E}_{j,k}\left (α\right ) with j < k. By Definition ELEM, these matrices are lower triangular with every diagonal entry equal to 1. We know that the product of two such matrices will again be lower triangular (Theorem PTMT), but also, as you can also easily check using a proof with a style similar to one above, that the product maintains all 1’s on the diagonal. Let {E}_{1},\kern 1.95872pt {E}_{2},\kern 1.95872pt {E}_{3},\kern 1.95872pt \mathop{\mathop{…}},\kern 1.95872pt {E}_{m} denote the elementary matrices for this sequence of row operations. Then

\eqalignno{ U & = {E}_{m}{E}_{m−1}\mathop{\mathop{…}}{E}_{3}{E}_{2}{E}_{1}A = {L}^{′}A & & }

where {L}^{′} is the product of the elementary matrices, and we know {L}^{′} is lower triangular with all 1’s on the diagonal. Our desired matrix L is then L ={ \left ({L}^{′}\right )}^{−1}. By Theorem ITMT, L is lower triangular with all 1’s on the diagonal and A = LU, as desired.

The process just described is deterministic. That is, the proof is constructive, with no freedom for each of us to walk through it differently. But could there be other matrices with the same properties as L and U that give such a decomposition of A. In other words, is the decomposition unique (Technique U)? Suppose that we have two triangular decompositions, A = {L}_{1}{U}_{1} and A = {L}_{2}{U}_{2}. Since A is nonsingular, two applications of Theorem NPNT imply that {L}_{1},\kern 1.95872pt {L}_{2},\kern 1.95872pt {U}_{1},\kern 1.95872pt {U}_{2} are all nonsingular. We have

\eqalignno{ {L}_{2}^{−1}{L}_{ 1} & = {L}_{2}^{−1}{I}_{ n}{L}_{1} & &\text{@(a href="fcla-jsmath-latestli31.html#theorem.MMIM")Theorem MMIM@(/a)} & & & & \cr & = {L}_{2}^{−1}A{A}^{−1}{L}_{ 1} & &\text{@(a href="fcla-jsmath-latestli32.html#definition.MI")Definition MI@(/a)} & & & & \cr & = {L}_{2}^{−1}{L}_{ 2}{U}_{2}{\left ({L}_{1}{U}_{1}\right )}^{−1}{L}_{ 1} & & & & \cr & = {L}_{2}^{−1}{L}_{ 2}{U}_{2}{U}_{1}^{−1}{L}_{ 1}^{−1}{L}_{ 1} & &\text{@(a href="fcla-jsmath-latestli32.html#theorem.SS")Theorem SS@(/a)} & & & & \cr & = {I}_{n}{U}_{2}{U}_{1}^{−1}{I}_{ n} & &\text{@(a href="fcla-jsmath-latestli32.html#definition.MI")Definition MI@(/a)} & & & & \cr & = {U}_{2}{U}_{1}^{−1} & &\text{@(a href="fcla-jsmath-latestli31.html#theorem.MMIM")Theorem MMIM@(/a)} & & & & }

Theorem ITMT tells us that {L}_{2}^{−1} is lower triangular and has 1’s as the diagonal entries. By Theorem PTMT, the product {L}_{2}^{−1}{L}_{ 1} is again lower triangular, and it is simple to check (as before) that the diagonal entries of the product are again all 1’s. By the entirely similar process we can conclude that the product {U}_{2}{U}_{1}^{−1} is upper triangular. Because these two products are equal, their common value is a matrix that is both lower triangular and upper triangular, with all 1’s on the diagonal. The only matrix meeting these three requirements is the identity matrix (Definition IM). So, we have,

\eqalignno{ {I}_{n} = {L}_{2}^{−1}{L}_{ 1} & ⇒ {L}_{2} = {L}_{1} &{I}_{n} = {U}_{2}{U}_{1}^{−1} & ⇒ {U}_{ 1} = {U}_{2} & & & & }

Studying the proofs of some previous theorems will perhaps give you an idea for an approach to computing a triangular decomposition. In the proof of Theorem CINM we augmented a nonsingular matrix with an identity matrix of the same size, and row-reduced until the original matrix became the identity matrix (as we knew in advance would happen, since we knew Theorem NMRRI). Theorem PEEF tells us about properties of extended echelon form, and in particular, that B = JA, where A is the matrix that begins on the left, and B is the reduced row-echelon form of A. The matrix J is the result on the right side of the augmented matrix, which is the result of applying the same row operations to the identity matrix. We should recognize now that J is just the product of the elementary matrices (Subsection DM.EM) that perform these row operations. Theorem ITMT used the extended echelon form to discern properties of the inverse of a triangular matrix. Theorem TD proves the existence of a triangular decomposition by applying specific row operations, and tracking the relevant elementary row operations. It is not a great leap to combine these observations into a computational procedure.

To find the triangular decomposition of A, augment A with the identity matrix of the same size and call this new 2n × n matrix, M. Perform row operations on M that convert the first n columns to an upper triangular matrix. Do this using only row operations that add a scalar multiple of one row to another row with higher index (i.e. lower down). In this way, the last n columns of M will be converted into a lower triangular matrix with 1’s on the diagonal (since M has 1’s in these locations initially). We could think of this process as doing about half of the work required to compute the inverse of A. Take the first n columns of the row-equivalent version of M and call this matrix U. Take the final n columns of the row-equivalent version of M and call this matrix {L}^{′}. Then by a proof employing elementary matrices, or a proof similar in spirit to the one used to prove Theorem PEEF, we arrive at a result similar to the second assertion of Theorem PEEF. Namely, U = {L}^{′}A. Multiplication on the left, by the inverse of {L}^{′}, will give us a decomposition of A (which we know to be unique). Ready? Lets try it.

Example TD4
Triangular decomposition, size 4
In this example, we will illustrate the process for computing a triangular decomposition, as described in the previous paragraphs. Consider the nonsingular square matrix A of size 4,

\eqalignno{ A & = \left [\array{ −2& 6 & −8 & 7\cr −4 &16 &−14 &15 \cr −6&22&−23&26\cr −6 &26 &−18 &17 } \right ] & & }

We form M by augmenting A with the size 4 identity matrix {I}_{4}. We will perform the allowed operations, column by column, only reporting intermediate results as we finish converting each column. It is easy to determine exactly which row operations we perform, since the final four columns contain a record of each such operation. We will not verify our hypotheses about the nonsingularity of the {A}_{k}, since if we do not have these conditions, we will reach a stage where a diagonal entry is zero and we cannot create the row operations we need to zero out the bottom portion of the associated column. In other words, we can boldly proceed and the necessity of our hypotheses will become apparent.

\eqalignno{ M = &\left [\array{ −2& 6 & −8 & 7 &1&0&0&0\cr −4 &16 &−14 &15 &0 &1 &0 &0 \cr −6&22&−23&26&0&0&1&0\cr −6 &26 &−18 &17 &0 &0 &0 &1 } \right ] & & \cr → &\left [\array{ −2&6&−8& 7 & 1 &0&0&0\cr 0 &4 & 2 & 1 &−2 &1 &0 &0 \cr 0 &4& 1 & 5 &−3&0&1&0\cr 0 &8 & 6 &−4 &−3 &0 &0 &1 } \right ] & & \cr → &\left [\array{ −2&6&−8& 7 & 1 & 0 &0&0\cr 0 &4 & 2 & 1 &−2 & 1 &0 &0 \cr 0 &0&−1& 4 &−1&−1&1&0\cr 0 &0 & 2 &−6 & 1 &−2 &0 &1 } \right ] & & \cr → &\left [\array{ −2&6&−8&7& 1 & 0 &0&0\cr 0 &4 & 2 &1 &−2 & 1 &0 &0 \cr 0 &0&−1&4&−1&−1&1&0\cr 0 &0 & 0 &2 &−1 &−4 &2 &1 } \right ] & & }

\eqalignno{ U & = \left [\array{ −2&6&−8&7\cr 0 &4 & 2 &1 \cr 0 &0&−1&4\cr 0 &0 & 0 &2 } \right ] &{L}^{′} & = \left [\array{ 1 & 0 &0&0\cr −2 & 1 &0 &0 \cr −1&−1&1&0\cr −1 &−4 &2 &1 } \right ] & & & & }

\eqalignno{ L & ={ \left ({L}^{′}\right )}^{−1} = \left [\array{ 1&0& 0 &0\cr 2&1 & 0 &0 \cr 3&1& 1 &0\cr 3&2 &−2 &1 } \right ] & & }

Subsection TDSSE: Triangular Decomposition and Solving Systems of Equations

In this section we give an explanation of why you might be interested in a triangular decomposition for a matrix. Many of the computational problems in linear algebra revolve around solving large systems of equations, or nearly equivalently, finding inverses of large matrices. Suppose we have a system of equations with coefficient matrix A and vector of constants b, and suppose further that A has the triangular decomposition A = LU.

Let y be the solution to the linear system ℒS\kern -1.95872pt \left (L,\kern 1.95872pt b\right ), so that by Theorem SLEMM, we have Ly = b. Notice that since L is nonsingular, this solution is unique, and the form of L makes it trivial to solve the system. The first component of y is determined easily, and we can continue on through determining the components of y, without even ever dividing. Now, with y in hand, consider the linear system, ℒS\kern -1.95872pt \left (U,\kern 1.95872pt y\right ). Let x be the unique solution to this system, so by Theorem SLEMM we have Ux = y. Notice that a system of equations with U as a coefficient matrix is also straightforward to solve, though we will compute the bottom entries of x first, and we will need to divide. The upshot of all this is that x is a solution to ℒS\kern -1.95872pt \left (A,\kern 1.95872pt b\right ), as we now show,

\eqalignno{ Ax & = LUx = L\left (Ux\right ) = Ly = b & & }

An application of Theorem SLEMM demonstrates that x is a solution to ℒS\kern -1.95872pt \left (A,\kern 1.95872pt b\right ).

Example TDSSE
Triangular decomposition solves a system of equations
Here we illustrate the previous discussion, recycling the decomposition found previously in Example TD4. Consider the linear system ℒS\kern -1.95872pt \left (A,\kern 1.95872pt b\right ) with

\eqalignno{ A & = \left [\array{ −2& 6 & −8 & 7\cr −4 &16 &−14 &15 \cr −6&22&−23&26\cr −6 &26 &−18 &17 } \right ] &b & = \left [\array{ −10\cr −2 \cr −1\cr −8 } \right ] & & & & }

First we solve the system ℒS\kern -1.95872pt \left (L,\kern 1.95872pt b\right ) (see Example TD4 for L),

\eqalignno{ {y}_{1} & = −10 & & \cr 2{y}_{1} + {y}_{2} & = −2 & & \cr 3{y}_{1} + {y}_{2} + {y}_{3} & = −1 & & \cr 3{y}_{1} + 2{y}_{2} − 2{y}_{3} + {y}_{4} & = −8 & & }

\eqalignno{ {y}_{1} & = −10 & & \cr {y}_{2} & = −2 − 2{y}_{1} = −2 − 2(−10) = 18 & & \cr {y}_{3} & = −1 − 3{y}_{1} − {y}_{2} = −1 − 3(−10) − 18 = 11 & & \cr {y}_{4} & = −8 − 3{y}_{1} − 2{y}_{2} + 2{y}_{3} = −8 − 3(−10) − 2(18) + 2(11) = 8 & & }

\eqalignno{ y & = \left [\array{ −10\cr 18 \cr 11\cr 8} \right ] & & }

Then we solve the system ℒS\kern -1.95872pt \left (U,\kern 1.95872pt y\right ) (see Example TD4 for U),

\eqalignno{ − 2{x}_{1} + 6{x}_{2} − 8{x}_{3} + 7{x}_{4} & = −10 & & \cr 4{x}_{2} + 2{x}_{3} + {x}_{4} & = 18 & & \cr − {x}_{3} + 4{x}_{4} & = 11 & & \cr 2{x}_{4} & = 8 & & }

\eqalignno{ {x}_{4}& = 8∕2 = 4 && \cr {x}_{3}& = \left (11 − 4{x}_{4}\right )∕(−1) = \left (11 − 4(4)\right )∕(−1) = 5 && \cr {x}_{2}& = \left (18 − 2{x}_{3} − {x}_{4}\right )∕4 = \left (18 − 2(5) − 4\right )∕4 = 1 && \cr {x}_{1}& = \left (−10 − 6{x}_{2} + 8{x}_{3} − 7{x}_{4}\right )∕(−2) = \left (−10 − 6(1) + 8(5) − 7(4)\right )∕(−2) = 2&& }

\eqalignno{ x & = \left [\array{ 4\cr 5 \cr 1\cr 2 } \right ] & & }

is the solution to ℒS\kern -1.95872pt \left (U,\kern 1.95872pt y\right ) and consequently is the unique solution to ℒS\kern -1.95872pt \left (A,\kern 1.95872pt b\right ), as you can easily verify. ⊠

Subsection CTD: Computing Triangular Decompositions

It would be a simple matter to adjust the algorithm for converting a matrix to reduced row-echelon form and obtain an algorithm to compute the triangular decomposition of the matrix, along the lines of Example TD4 and the discussion preceding this example. However, it is possible to obtain relatively simple formulas for the entries of the decomposition, and if computed in the proper order, an implementation will be straightforward. We will state the result as a theorem and then give an example of its use.

Theorem TDEE
Triangular Decomposition, Entry by Entry
Suppose that A is a squarematrix of size n with a triangular decomposition A = LU, where L is lower triangular with diagonal entries all equal to 1, and U is upper triangular. Then

\eqalignno{ {\left [U\right ]}_{ij} & ={ \left [A\right ]}_{ij} −{\mathop{∑ }}_{k=1}^{i−1}{\left [L\right ]}_{ ik}{\left [U\right ]}_{kj} & &1 ≤ i ≤ j ≤ n & & & & \cr {\left [L\right ]}_{ij} & = {1\over { \left [U\right ]}_{jj}}\left ({\left [A\right ]}_{ij} −{\mathop{∑ }}_{k=1}^{j−1}{\left [L\right ]}_{ ik}{\left [U\right ]}_{kj}\right ) & &1 ≤ j < i ≤ n & & & & }

Proof Consider a single scalar product of an entry of L with an entry of U of the form {\left [L\right ]}_{ik}{\left [U\right ]}_{kj}. By Definition LTM, if k > i then {\left [L\right ]}_{ik} = 0, while Definition UTM, says that if k > j then {\left [U\right ]}_{kj} = 0. So we can combine these two facts to assert that if k >\mathop{ min}(i,\kern 1.95872pt j), {\left [L\right ]}_{ik}{\left [U\right ]}_{kj} = 0 since at least one term of the product will be zero. Employing this observation,

\eqalignno{ {\left [A\right ]}_{ij} & ={ \mathop{∑ }}_{k=1}^{n}{\left [L\right ]}_{ ik}{\left [U\right ]}_{kj} & &\text{@(a href="fcla-jsmath-latestli31.html#theorem.EMP")Theorem EMP@(/a)} & & & & \cr & ={ \mathop{∑ }}_{k=1}^{\mathop{min}(i,\kern 1.95872pt j)}{\left [L\right ]}_{ ik}{\left [U\right ]}_{kj} & & & & }

\eqalignno{ {\left [U\right ]}_{ij} & ={ \left [A\right ]}_{ij} −{\left [A\right ]}_{ij} +{ \left [U\right ]}_{ij} & & \cr & ={ \left [A\right ]}_{ij} −{\mathop{∑ }}_{k=1}^{\mathop{min}(i,\kern 1.95872pt j)}{\left [L\right ]}_{ ik}{\left [U\right ]}_{kj} +{ \left [U\right ]}_{ij} & & \cr & ={ \left [A\right ]}_{ij} −{\mathop{∑ }}_{k=1}^{i}{\left [L\right ]}_{ ik}{\left [U\right ]}_{kj} +{ \left [U\right ]}_{ij} & & \cr & ={ \left [A\right ]}_{ij} −{\mathop{∑ }}_{k=1}^{i−1}{\left [L\right ]}_{ ik}{\left [U\right ]}_{kj} −{\left [L\right ]}_{ii}{\left [U\right ]}_{ij} +{ \left [U\right ]}_{ij} & & \cr & ={ \left [A\right ]}_{ij} −{\mathop{∑ }}_{k=1}^{i−1}{\left [L\right ]}_{ ik}{\left [U\right ]}_{kj} −{\left [U\right ]}_{ij} +{ \left [U\right ]}_{ij} & & \cr & ={ \left [A\right ]}_{ij} −{\mathop{∑ }}_{k=1}^{i−1}{\left [L\right ]}_{ ik}{\left [U\right ]}_{kj} & & }

\eqalignno{ {\left [L\right ]}_{ij} & = {1\over { \left [U\right ]}_{jj}}\left ({\left [L\right ]}_{ij}{\left [U\right ]}_{jj}\right ) & & \cr & = {1\over { \left [U\right ]}_{jj}}\left ({\left [A\right ]}_{ij} −{\left [A\right ]}_{ij} +{ \left [L\right ]}_{ij}{\left [U\right ]}_{jj}\right ) & & \cr & = {1\over { \left [U\right ]}_{jj}}\left ({\left [A\right ]}_{ij} −{\mathop{∑ }}_{k=1}^{\mathop{min}(i,\kern 1.95872pt j)}{\left [L\right ]}_{ ik}{\left [U\right ]}_{kj} +{ \left [L\right ]}_{ij}{\left [U\right ]}_{jj}\right ) & & \cr & = {1\over { \left [U\right ]}_{jj}}\left ({\left [A\right ]}_{ij} −{\mathop{∑ }}_{k=1}^{j}{\left [L\right ]}_{ ik}{\left [U\right ]}_{kj} +{ \left [L\right ]}_{ij}{\left [U\right ]}_{jj}\right ) & & \cr & = {1\over { \left [U\right ]}_{jj}}\left ({\left [A\right ]}_{ij} −{\mathop{∑ }}_{k=1}^{j−1}{\left [L\right ]}_{ ik}{\left [U\right ]}_{kj} −{\left [L\right ]}_{ij}{\left [U\right ]}_{jj} +{ \left [L\right ]}_{ij}{\left [U\right ]}_{jj}\right ) & & \cr & = {1\over { \left [U\right ]}_{jj}}\left ({\left [A\right ]}_{ij} −{\mathop{∑ }}_{k=1}^{j−1}{\left [L\right ]}_{ ik}{\left [U\right ]}_{kj}\right ) & & }

At first glance, these formulas may look exceedingly complex. Upon closer examination, it looks even worse. We have expressions for entries of U that depend on other entries of U and also on entries of L. But then the formula for entries of L depend on entries from L and entries from U. Do these formula have circular dependencies? Or perhaps equivalently, how do we get started? The key is to be organized about the computations and employ these two (similar) formulas in a specific order. First compute the first row of L, followed by the first column of U. Then the second row of L, followed by the second column of U. And so on. In this way, all of the values required for each new entry will have already been computed previously.

Of course, the formula for entries of L require division by diagonal entries of U. These entries might be zero, but in this case A is nonsingular and does not have a triangular decomposition. So we need not check the hypothesis carefully and can launch into the arithmetic dictated by the formulas, confident that we will be reminded when a decomposition is not possible. Note that these formula give us all of the values that we need for the decomposition, since we require that L has 1’s on the diagonal. If we replace the 1’s on the diagonal of L by zeros, and add the matrix U, we get an n × n matrix containing all the information we need to resurrect the triangular decomposition. This is mostly a notational convenience, but it is a frequent way of presenting the information. We’ll employ it in the next example.

\eqalignno{ A & = \left [\array{ 3 & 3 & −3 &−2& −1 & 0\cr −6 & −4 & 5 & 2 & 4 & 2 \cr 9 & 9 & −7 &−7& 0 & 1\cr −6 &−10 & 8 & 10 & −1 &−7 \cr 6 & 4 & −9 &−2&−10& 1\cr 9 & 3 &−12 &−3 &−21 &−2 } \right ] & & }

Using the notational convenience of packaging the two triangular matrices into one matrix, and using the ordering of the computations mentioned above, we display the results after computing a single row and column of each of the two triangular matrices.

\eqalignno{ &\left [\array{ 3 &3&−3&−2&−1&0\cr −2 \cr 3\cr −2 \cr 2\cr 3 } \right ] & &\left [\array{ 3 & 3 &−3&−2&−1&0\cr −2 & 2 &−1 &−2 & 2 &2 \cr 3 & 0\cr −2 &−2 \cr 2 &−1\cr 3 &−3 } \right ] & & & & \cr &\left [\array{ 3 & 3 &−3&−2&−1&0\cr −2 & 2 &−1 &−2 & 2 &2 \cr 3 & 0 & 2 &−1& 3 &1\cr −2 &−2 & 0 \cr 2 &−1&−2\cr 3 &−3 &−3 } \right ] & &\left [\array{ 3 & 3 &−3&−2&−1& 0\cr −2 & 2 &−1 &−2 & 2 & 2 \cr 3 & 0 & 2 &−1& 3 & 1\cr −2 &−2 & 0 & 2 & 1 &−3 \cr 2 &−1&−2&−1\cr 3 &−3 &−3 &−3 } \right ] & & & & \cr &\left [\array{ 3 & 3 &−3&−2&−1& 0\cr −2 & 2 &−1 &−2 & 2 & 2 \cr 3 & 0 & 2 &−1& 3 & 1\cr −2 &−2 & 0 & 2 & 1 &−3 \cr 2 &−1&−2&−1& 1 & 2\cr 3 &−3 &−3 &−3 & 0 } \right ] & &\left [\array{ 3 & 3 &−3&−2&−1& 0\cr −2 & 2 &−1 &−2 & 2 & 2 \cr 3 & 0 & 2 &−1& 3 & 1\cr −2 &−2 & 0 & 2 & 1 &−3 \cr 2 &−1&−2&−1& 1 & 2\cr 3 &−3 &−3 &−3 & 0 &−2 } \right ] & & & & }

\eqalignno{ L & = \left [\array{ 1 & 0 & 0 & 0 &0&0\cr −2 & 1 & 0 & 0 &0 &0 \cr 3 & 0 & 1 & 0 &0&0\cr −2 &−2 & 0 & 1 &0 &0 \cr 2 &−1&−2&−1&1&0\cr 3 &−3 &−3 &−3 &0 &1 } \right ] &U & = \left [\array{ 3&3&−3&−2&−1& 0\cr 0&2 &−1 &−2 & 2 & 2 \cr 0&0& 2 &−1& 3 & 1\cr 0&0 & 0 & 2 & 1 &−3 \cr 0&0& 0 & 0 & 1 & 2\cr 0&0 & 0 & 0 & 0 &−2 } \right ] & & & & }

The hypotheses of Theorem TD can be weakened slightly to include matrices where not every {A}_{k} is nonsingular. The introduces a rearrangement of the rows and columns of A to force as many as possible of the smaller submatrices to be nonsingular. Then permutation matrices also enter into the decomposition. We will not present the details here, but instead suggest consulting a more advanced text on matrix analysis.