7 Direct Solution of Linear Systems

7.1 Linear Systems

We consider for \(n\in\mathbb{N},\) \(n\geq 1,\) the linear system of equations \[\tag{7.1} \begin{matrix} a_{11}x_1 &+ \cdots + & a_{1n}x_n & = b_1\\ \vdots & &\vdots & \vdots\\ a_{n1}x_1 &+ \cdots + & a_{nn}x_n & = b_n\\ \end{matrix}\] with the coefficients \(a_{ij}\in\mathbb{R},\) \(i,j=1,\ldots,n,\) the unknowns \(x_i\in\mathbb{R},\) \(i=1,\ldots,n,\) and the values of the right hand side \(b_i\in\mathbb{R},\) \(i=1,\ldots,n.\) We would like to determine the values of the \(x_i\) for given \(b_i\) and \(a_{ij}.\)

As an example for a linear system, we can take the Vandermonde Matrix, which we have encountered during polynomial interpolation, cf. Equation (6.11).

For the case \(n=2,\) we may also consider the following system \[\begin{aligned} 2 x_1 & - \phantom{2}x_2 = 0\\ -x_1 & +2 x_2 =3 \end{aligned}\] whose solution is given by \(x_1=1, x_2=2.\) For this small \(2\)-dimensional system, we can compute the solution of the system readily by, e.g. solving in the second equation for the first variable and inserting the result \(x_1=2x_2 - 3\) into the first equation, this procedure becomes tedious and unpractical for larger \(n.\)

In general, we will answer the following questions connected to the linear system (7.1):

• Does the system (7.1) have a solution?

• If the system (7.1) has a solution, is it unique or can there be several solutions?

• If a solution exists, how can we compute it in an efficient and stable way?

• How do the answers to these questions depend on the coefficients and the right hand side?

In order to answer these questions, we would like to derive more systematic ways of analysing and solving linear systems.

7.2 Triangular Systems and Forward/Backward Substitution

As a first step, we consider a linear system of the following form, which we call lower triangular: \[\tag{7.2} \begin{matrix} l_{11}x_1 & & & & & = &b_1\\ l_{21}x_1 & +& l_{22}x_2 & & & = &b_2\\ \vdots & & & & \ddots & &\vdots\\ l_{n1}x_1 & +& l_{n2}x_2 &+&l_{nn}x_n & = &b_n\, .\\ \end{matrix}\] Here, for the coefficients \(l_{ij}\in\mathbb{R},\) \(i,j=1,\ldots,n,\) we have that \(l_{ij} = 0\) for \(j>i.\) Such a lower triangular system can be easily solved by starting with the first row and computing \[x_1=l_{11}^{-1}b_1\,,\] which can be done for \(l_{11}\not = 0.\) The resulting value for \(x_1\) can then be inserted into the second row, leading to \[x_2=l_{22}^{-1}(b_2-l_{21}x_1)\,,\] which can be done for \(l_{22}\not = 0.\) Continuing with this procedure, we see that we can compute \(x_i\) as \[\tag{7.3} x_i = l_{ii}^{-1}\left (b_i- \sum_{j=1}^{i-1} l_{ij}x_j\right) \quad \text{ for } i=1,\ldots,n\, .\] Clearly, this requires \(l_{ii}\not = 0\) for \(i=1,\ldots,n.\) We call the procedure (7.3) forward substitution. Counting, we see that forward substitution requires \(\mathcal O(n^2)\) floating point operations.

The according procedure for an upper triangular system, i.e. a system with coefficients \(r_{ij} = 0\) for \(j<i,\) is called backward substitution. Backward substitution starts in the last row with computing \(x_n\) and computes the values of \(x_1\) last.

For lower (upper) triangular systems, we can therefore already partially answer our questions from above: If \(\Pi_{i=1}^n l_{ii} \not = 0\) (\(\Pi_{i=1}^n r_{ii} \not = 0\)), the lower (upper) triangular system has a unique solution \(x_1,\ldots,x_n,\) which can be computed by means of forward (backward) substitution.

7.3 Gaussian Elimination

We write the coefficients and the right hand side of the linear system (7.1) in its extended form \[\tag{7.4} \begin{matrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} & b_1\\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n} & b_2\\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n} & b_3\\ \vdots & \vdots & &\vdots & \vdots\\ a_{n1} & a_{n2} & a_{n3} & \cdots & a_{nn} & b_n\\ \end{matrix}\, .\] We will denote the rows in the above extended form by \(R_1,\ldots,R_n,\) e.g. \[R_i = [a_{i1} \cdots a_{in}|b_i]\, .\] We observe that, due to linearity, the following operations will not change the solution \(x\) of our linear system (7.1):

• swapping rows \(R_i,R_j,\) \(i,j=1,\ldots,n\)

• multiplying a row with a scalar, i.e. \(R_i\mapsto \alpha R_i,\) \(0\not=\alpha\in\mathbb{R},\) \(i=1,\ldots,n.\)

• adding the multiple of one row to another row, i.e. \(R_i\mapsto R_i+\alpha R_j,\) \(0\not=\alpha\in\mathbb{R},\) \(i,j=1,\ldots,n.\)

We exploit this in order to transform our system into upper triangular form, so that we can apply backward substitution.

We start by eliminating all non-zero entries in the tableau in the first column starting from the second row. To this end, let us assume that \(a_{11}\not = 0.\) Otherwise, we would select a row \(i\) with \(a_{i1}\not = 0\) and would swap row \(i\) and row \(1.\) We call \(a_{11}\) our pivot element.

We start with eliminating \(a_{21}.\) By setting \[\alpha_{21} = - \frac{a_{21}}{a_{11}}\, ,\] we see that with \[\tag{7.5} R_2\mapsto R_2+\alpha_{21} \cdot R_1 = R_2 - \frac{a_{21}}{a_{11}} \cdot R_1\] the second row is transformed into \[R_2 +\alpha_{21} R_1 = [a_{21} - \frac{a_{21}}{a_{11}} a_{11} \,\, *\cdots *|*] = [0 \,\, *\cdots *|*]\] and that our extended form (7.4) has been transformed into \[\begin{matrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} & b_1\\ 0 & a_{22}^1 & a_{23}^1 & \cdots & a_{2n}^1 & b_2^1\\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n} & b_3\\ \vdots & \vdots & &\vdots & \vdots\\ a_{n1} & a_{n2} & a_{n3} & \cdots & a_{nn} & b_n\\ \end{matrix}\] Here, the superscript \({}^1\) indicates that the corresponding values of the coefficients and the right hand side have been changed in this first elimination step (7.5). We now continue with the second row. Setting \[\alpha_{31} = - \frac{a_{31}}{a_{11}}\, ,\] we see that with \[\tag{7.6} R_3\mapsto R_3+\alpha_{31} \cdot R_1 = R_3 - \frac{a_{31}}{a_{11}} \cdot R_1\] we arrive at \[\begin{matrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} & b_1\\ 0 & a_{22}^1 & a_{23}^1 & \cdots & a_{2n}^1 & b_2^1\\ 0 & a_{32}^2 & a_{33}^2 & \cdots & a_{3n}^2 & b_3^2\\ \vdots & \vdots & &\vdots & \vdots\\ a_{n1} & a_{n2} & a_{n3} & \cdots & a_{nn} & b_n\\ \end{matrix}\, .\] We continue this process, until we have eliminated all entries of the first column below the first row, leading to \[\begin{matrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} & b_1\\ 0 & a_{22}^1 & a_{23}^1 & \cdots & a_{2n}^1 & b_2^1\\ 0 & a_{32}^2 & a_{33}^2 & \cdots & a_{3n}^2 & b_3^2\\ \vdots & \vdots & &\vdots & \vdots\\ 0 & a_{n2}^{n-1} & a_{n3}^{n-1} & \cdots & a_{nn}^{n-1} & b_n^{n-1}\\ \end{matrix}\, .\] We continue in the second row. Again, if \(a_{22}^1 = 0,\) we swap the second row with a row with a leading non-zero element in the second column. If such a non-zero pivot element does not exist, we proceed with the next column, if there is one.

Assuming that \(a_{22}^1\not = 0,\) we proceed as above. By setting \[\alpha_{32} = - \frac{a_{32}}{a_{22}^1}\, ,\] we see that with \[\tag{7.7} R_3\mapsto R_3+\alpha_{32} \cdot R_2 = R_3 - \frac{a_{32}}{a_{22}^1} \cdot R_2\] we arrive at \[\begin{matrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} & b_1\\ 0 & a_{22}^1 & a_{23}^1 & \cdots & a_{2n}^1 & b_2^1\\ 0 & 0 & a_{33}^n & \cdots & a_{3n}^n & b_3^n\\ \vdots & \vdots & &\vdots & \vdots\\ 0 & a_{n2}^{n-1} & a_{n3}^{n-1} & \cdots & a_{nn}^{n-1} & b_n^{n-1}\\ \end{matrix}\] We now continue this process until the tableau has the desired upper triangular form \[\begin{matrix}[cccccc|l] * & * & * & \cdots & * & * \\ 0 & * & * & \cdots & * & * \\ 0 & 0 & * & \cdots & * & * \\ \vdots & \vdots & &\ddots&\vdots & \vdots\\ 0 & 0 & 0 & \cdots & * & * \\ \end{matrix}\] If we cannot find a non-zero leading element —or pivot element— in column \(i,\) we proceed with the next column.

7.4 Matrices and Vectors

As a first step towards a systematic treatment of linear systems, we introduce the following notations.

7.4.1 Linear Combinations of Matrices and Vectors

We now introduce the multiplication of matrices with vectors and the multiplication of matrices with matrices. These operations will help us to write our linear system (7.1) in a compact form using matrices and vectors. They will furthermore be very helpful for analysing and solving linear systems.

7.4.2 Matrix-Vector Product

7.4.3 Matrix-Matrix Product

As handling the the different indices and sizes can be confusing in the beginning, we write down the matrix-matrix product for the simpler case of quadratic matrices, i.e. \(m=n=l.\) This case will be of primary interest for us in the following \[\begin{aligned} A\cdot B & = \begin{pmatrix} \sum\limits_{j=1}^n a_{1j}\cdot b_{j1} & \cdots & \sum\limits_{j=1}^n a_{1j}\cdot b_{jn}\\ \vdots & &\vdots \\ \sum\limits_{j=1}^n a_{nj}\cdot b_{j1} & \cdots & \sum\limits_{j=1}^n a_{nj}\cdot b_{jn}\\ \end{pmatrix}\\ & = \begin{pmatrix} a_{11}b_{11} + \cdots + a_{1n}b_{n1}& \cdots &a_{11}b_{1n} + \cdots + a_{1n}b_{nn} \\ \vdots & &\vdots \\ a_{n1}b_{11} + \cdots + a_{nn}b_{n1} & \cdots & a_{n1}b_{1n} + \cdots + a_{nn}b_{nn}\\ \end{pmatrix} \, . \end{aligned}\] We note that the matrix-vectorx multiplication (7.15) is the special case of matrix-matrix multiplication (7.17) with \(l=1.\)

7.4.4 Linearity

An important property of the matrix-vector product is the fact that it is linear. More precisely, for \(x,y\in\mathbb{R}^n,\) \(A\in \mathbb{R}^{m\times n}\) and \(\alpha,\beta\in\mathbb{R},\) we have \[\begin{aligned} A(\alpha x+ \beta y) = \alpha Ax + \beta Ay\, . \end{aligned}\] This can be seen from the definition of the matrix-vector (7.16) product and the definition of vector addition (7.14) for \(i=1,\ldots,m\) as follows \[\begin{aligned} \big(A(\alpha x+\beta y) \big)_i & = \sum\limits_{j=1}^n a_{ij}\cdot (\alpha x+\beta y)_j\\ & = \sum\limits_{j=1}^n a_{ij}\cdot (\alpha x_j+\beta y_j)\\ & = \alpha \sum\limits_{j=1}^n a_{ij}\cdot x_j+\beta \sum\limits_{j=1}^n a_{ij}\cdot y_j\\ &= \alpha Ax + \beta Ay\, . \end{aligned}\]

7.4.5 Identity and Inverse Matrix

The matrix product (7.17) has a neutral element, the identity matrix

The identity matrix is called identity matrix because it is the neutral element of the matrix multiplication (7.17). More precisely, let \(A\in\mathbb{R}^{n\times n}.\) Then we have \[\tag{7.21} I_n\cdot A = A\cdot I_n = A\, .\] Furthermore, if \(x\in \mathbb{R}^n,\) we also see from the definition that \[\tag{7.22} I_n\cdot x = x \, .\] The fact that we have identified a neutral element of the matrix multiplication leads us to the following definition

Using Definition 7.4, we can now write our linear system (7.1) in the compact and elegant form \[\tag{7.24} \begin{aligned} Ax &= b \end{aligned}\] with the (coefficient) matrix \(A\in \mathbb{R}^{n\times n},\) the right hand side \(b\in \mathbb{R}^n\) and the vector of unknowns \(x\in\mathbb{R}^n.\)

We also see that if \(A\) is invertible, we can formally solve the linear system (7.24) by multiplying from the left with \(A^{-1}\) \[A^{-1}\cdot A x = I_n x = x = A^{-1}b\, .\]

7.4.6 Properties of Triangular Matrices

For \(A\) being upper (lower) triangular — and we write \(A=T\) in order to emphasize that we are talking about a triangular matrix — we already have seen that backward (forward) substitution can be carried out for solving the upper (lower) triangular system \[Tx=b\, ,\] if all diagonal elements of \(T\) are different from zero. This condition gives rise to the following definition:

We emphasize that definition (7.25) of the determinant applies only to triangular matrices.

As we can see from Equation (7.3), forward or backward substitution is applicable to a triangular matrix \(T\) if and only if \(\det (T)\not = 0.\) For triangular matrices, we therefore can formulate

We now collect some properties of triangular matrices.

Let us now go back to the solution of the linear system (7.24). Although it can be formally solved by computing the solution as \(x=A^{-1}b,\) in practice, we will never compute the inverse, \(A^{-1}\) , as this is numerically inconvenient and in general not needed. We will instead transform the linear system (7.1) into triangular form and then apply forward or backward substitution. This is done by means of Gaussian elimination.

7.5 Elementary Matrix-Transformations

Gaussian elimination uses elementary operations on matrices, such as adding and multiplying rows. These elementary operations we would like to collect now for later use and reference.

We define some elementary operations on matrices, which we will use as building blocks for our algorithms.

7.5.0.0.1 Scaling

Let \(D=\mathop{\mathrm{diag}}(d_1, \ldots ,d_n)\) with \(d_i\neq0\) for \(i=1, \ldots ,n.\) We see \[\begin{aligned} (Dx)_i&=&d_ix_i\\ (DA)_{ij}&=&a_{ij} \cdot d_i\\ (AD)_{ij}&=&a_{ij}\cdot d_j \end{aligned}\] The inverse of \(D\) is given by \[D^{-1}=\mathop{\mathrm{diag}}(d_1^{-1}, \ldots ,d_n^{-1}).\] For the determinant we have \[\det{D}=d_1\cdot\ldots\cdot d_n.\] For similarity transformations with \(D\) we get \[A\mapsto DAD^{-1}\quad\text{ with }a_{ij}\mapsto\frac{d_i}{d_j}a_{ij}.\]

7.5.0.0.2 Permutation

We call \[P_{ij} = \begin{pmatrix} I\\ &0& &1\\ & &I\\ &1& &0\\ & & & &I\end{pmatrix}\] a Permutation Matrix. We can see that

• \(P_{ij}x\) swaps \(x_i\) and \(x_j\);

• \(P_{ij}A\) swaps row \(i\) and row \(j\);

• \(AP_{ij}\) swaps columns \(i\) and \(j.\)

The inverse of \(P\) is \[P_{ij}^{-1}=P_{ij}.\] The determinant can be computed as \[\det{P_{ij}}=-1.\] For similarity transformations we have \[A\mapsto P_{ij}AP_{ij}\] with \[\begin{pmatrix} a_{ii}&a_{ij}\\ a_{ji}&a_{jj}\end{pmatrix} \mapsto \begin{pmatrix} a_{jj}&a_{ji}\\ a_{ij}&a_{ii}\end{pmatrix}.\] Furthermore, we have \[P_{ij}^T=P_{ij}.\]

7.5.0.0.3 Elimination Operations

Consider the following operation \[\tag{7.27} N_{ij}(\alpha)=\begin{pmatrix} 1\\ &\ddots\\ & &\ddots\\ &\alpha& &\ddots\\ & & & &1\end{pmatrix}\, ,\] where \(\alpha\) is in the \(i\)-th row and the \(j\)-th column. We see that

• \(N_{ij}(\alpha)x\) adds \(\alpha\)-times component \(j\) to component \(i\) of \(x.\)

• \(N_{ij} A\) adds \(\alpha\) times row \(j\) to row \(i\) of \(A.\)

• \(A N_{ij}\) add \(\alpha\)-times column \(j\) to column \(i\) of \(A.\)

For the inverse, we get \[N_{ij}(\alpha)^{-1}=N_{ij}(-\alpha).\] The determinant can be computed as \[\det{N_{ij}(\alpha)}=1\] and for similarity transformations we have \[A\mapsto N_{ij}(\alpha)AN_{ij}(-\alpha).\]

In order to understand the action of the transformation \(N_{ij}(\alpha)\) better, we can write \[N_{ij}(\alpha) = I_n + W\] with \(W=(w_{kl})_{k,l=1,\ldots,n}\) and \[w_{kl} = \begin{cases} \alpha & \text{if } k = i, l = j, \\ 0 & \text{if } k \neq i, l \neq j. \end{cases}\]

This leads to \[(WA)_{ml} = \sum_{k=1}^{n} a_{mk}w_{kl}\, ,\] which can be different from zero only for \(l=j\) by the definition of \(W.\) In this case, we have, again by definition of \(W,\) \[(WA)_{mj} = \sum_{k=1}^{n} a_{mk}w_{kj} = \alpha a_{mi}\, .\] We have seen that the matrix \(WA\) contains non-zero elements in row \(i,\) more precisely: the \(i\)-th row of \(WA\) is \(\alpha\)-times the row \(j\) of \(A\) and all other rows are zero.

The following theorem will turn out to be useful. As the proof is technical and for the moment being not necessary for going ahead with our arguments, we postpone it.

7.6 LU-Decomposition

Let \(A\in\mathbb{R}^{n\times n}.\) Our goal is to write \(A=LU\) with a lower triangular matrix \(L\) with \(l_{ii}=1\) for \(i=1, \ldots ,n,\) and an upper triangular matrix \(U.\)

Our goal is to transform \(A\) into an upper-triangular matrix. We notice that the elementary transformations we have listed above, in particular the swapping of two rows and adding the multiple of row \(i\) to row \(j\) does not change the solution of the linear system, if the right hand side is transformed accordingly.

We therefore assume that \(a_{11}\not= 0.\) If this is not the case, we choose a row \(i\) with \(a_{i1}\not = 0\) and swap the rows \(1\) and \(i.\) Then we add multiples of the first row to all rows below the first row, i.e. to the rows \(2,\ldots, n.\) Our goal is to eliminate all entries in the first column, starting from the second row. This means that we would like to achieve \((a_{i1})_{i=2,\ldots,n}=0.\)

As we can see, for \(a_{11}\neq0\) we can set \(l_{21}=\frac{a_{21}}{a_{11}}.\) Now, we multiply the first row with \(l_{21}\) and subtract the result from the second row. We get in particular for \(a_{21}\) that \[a_{21} \mapsto a_{21} - l_{21} \cdot a_{11} = a_{21} - \frac{a_{21}}{a_{11}} a_{11} = 0\, .\]

Using our elementary operation (7.27) \(N_{21}(\alpha),\) which adds \(\alpha\)-times row \(1\) to row \(2,\) we see that we can write the above equation equivalently as \[\big(N_{21}(-l_{21})A\big)_{21}=0\, .\] We therefore decide as a first operation to apply \(N_{21}(-l_{21}\)) from the left to \(A,\) i.e. we add \(\alpha\)-times row \(1\) to row \(2\) and call this elimination step \(1.\) We denote the resulting matrix by adding the superscript \(A^1\): \[A^1 = N_{21}(-l_{21})\cdot A^0\, .\] Here, we have set \(A=A^0.\) We now continue with row \(3,\) leading to \[A^2 = N_{31}(-l_{31})\cdot A^1 = N_{31}(-l_{31})\cdot N_{21}(-l_{21})\cdot A^0\, .\] where we have set \(l_{31}=\frac{a^1_{31}}{a_{11}}.\)

We then repeat this procedure until we have eliminated all entries in the first column below the first row, i.e. until \((a_{i1})_{i=2,\ldots,n}=0.\) In general, this will affect all matrix entries below the first row.

Then, we continue with the second column and eliminate all entries in the second column below the second row, i.e. we eliminate \((a_{i2})_{i=3,\ldots,n}=0.\) For simplicity, we drop the superscript and simply replace the respective old values of \(A\) with the computed new ones.

The resulting algorithm can be written in pseudo-code as follows:

• for \((j=1, \ldots ,n-1) \{\)

  • if \(a_{jj} = 0\)

    • stop (no success)

  • for \(i=(j+1, \ldots ,n) \{\)

    • set \(l_{ij}:=\frac{a_{ij}}{a_{jj}};\)

    • \(\qquad~ A:=N_{ij}(-l_{ij})A;\)

  • \(\}\)

• \(\}\)

Before step \((i,j),\) which eliminates the entry at position \((i,j),\) the matrix \(A\) has the form \[\begin{pmatrix} \ast&\ast&\ast&\ast&\ast&\ast&\ast&\ast\\ &\ast&\ast&\ast&\ast&\ast&\ast&\ast\\ & &\ast&\ast&\ast&\ast&\ast&\ast\\ & & &\ast&\ast&\ast&\ast&\ast\\ & & & &\ast&\ast&\ast&\ast\\ & & &a_{ij}&\ast&\ast&\ast&\ast\\ & & &\ast&\ast&\ast&\ast&\ast\\ & & &\ast&\ast&\ast&\ast&\ast\end{pmatrix}.\] The result of a successful run of the above algorithm is the upper triangular matrix \[\tag{7.28} U = N_{n,n-1}(-l_{n,n-1})\cdots N_{21}(-l_{21})A \, .\] We now recall that \(N_{ij}^{-1}(\alpha)=N_{ij}(-\alpha).\) We therefore set \[\tag{7.29} L = N_{21}(l_{21})\cdots N_{n,n-1}(l_{n,n-1})\] and see that multiplying Equation (7.28) from the left with \(L\) means multiplying successively from the left by \(N_{n,n-1}(l_{n,n-1}),\ldots,N_{21}(l_{21}).\) Thus, we arrive at \[\begin{aligned} LU & = \\ &= N_{21}(l_{21})\cdots N_{n,n-1}(l_{n,n-1}) N_{n,n-1}(-l_{n,n-1})\cdots N_{21}(-l_{21})A\\ &= N_{21}(l_{21})\cdots I_n \cdots N_{21}(-l_{21})A\\ & \,\,\,\vdots\\ &= A \end{aligned}\] or in short \[\tag{7.30} LU = A\, .\] Theorem 7.5 ensures that \(L\) is of lower triangular form. More precisely, we exactly know \(L,\) as it is given by \[L=\begin{pmatrix} 1\\ l_{21}&\ddots\\ \vdots& &\ddots\\ l_{n1}&\cdots&l_{n,n-1}&1\end{pmatrix} \,.\] We therefore can now solve our original system (7.1) by using the above algorithm.

We emphasize that in our above considerations we have assumed that in step \(j\) for the element \(a_{jj}\neq0\) holds. If this is not the case, the above algorithm stops. Consequently, a decomposition of the form \(A=LU\) does not necessarily always exist. We will see in the next section how to resolve this problem.

In terms of computational complexity, we count approximately \(\frac13n^3\) point operations (multiplication or division) and \(\frac23n^3\) operations with additions.

7.7 Practical Considerations

If the \(LU\)-decomposition of a matrix \(A\) is known, it can be used for computing the solutions for different right hand sides. Reusing the \(LU\)-decomposition requires only \({\mathcal O}(n^2)\) operations instead of \({\mathcal O}(n^3)\)

We can compute the determinant of \(A\) using the \(LU\)-decomposition as \[\det(A)=\det(LU)=\det(L)\det(U)=\det(U)=\prod_{i=1}^n r_{ii}\, .\] We note that Cramer’s rule for the determinant is numerically not stable and has complexity \({\mathcal O}(n!)\) and therefore should be avoided.

By design, Gaussian elimination creates so-called “fill-in” during the elimination. The original sparse structure of the matrix is destroyed during elimination, as non-zero values are created through the algebraic operations on the rows. Fill-in can be reduced by

• Swapping rows and columns for optimizing the sparsity pattern

• Dropping small entries which show up during the elimination. This leads to an incomplete \(LU\) decomposition.

7.8 LU-decomposition in Detail

A decomposition of the form \(A=LU\) for invertible \(A\) in general will not always exist. Moroever, instabilities can occur during the computation. We therefore consider

• Pivoting for the \(LU\)-decomposition

• Stability of the \(LU\)-decomposition

• Scaling and iterative refinement.

The above considerations motivate the following pivoting strategy:

Column pivoting or swapping of rows with pivot-search in the columns.

Let \(A^{(0)} = A.\) Choose in elimination step \(j\) an index \(p\in\{j, \ldots ,n\}\) with \[|a_{pj}^{(j-1)}|\ge |a_{kj}^{(j-1)}|\quad\text{ for } k=j, \ldots ,n\, .\] Let \(\tilde{A}^{(j-1)}\) be the matrix obtained by swapping row \(i\) and \(p\) in the matrix \(A^{(j-1)}.\) \[A^{(j)}= \left( \prod\limits_{i>j}N_{ij}(-l_{ij})\right) \tilde{A}^{(j-1)}\] Our strategy has led to \[|l_{ij}|= |\frac{\tilde{a}_{ij}^{(j-1)}}{\tilde{a}_{jj}^{(j-1)}}|\le1\, ,\] since \[|\frac{a_{ij}^{(j-1)}}{a_{pj}^{(j-1)}}| \le 1 \quad\text{ for all } i=j, \ldots ,n\, .\] In addition to column pivoting, also row pivoting is possible. Both strategies require \(\mbox{\large$\mathcal{O}$}(n^2)\) operations. A complete pivot search (rows and columns) would require \(\mbox{\large$\mathcal{O}$}(n^3)\) operations and therefore is rarely used in practice.

Algorithmics

Home

Chapters

Contents

---

✕