Find Lu Factorizationof a and Complete Solution to Ax B

Linear Equations and Eigensystems

George Lindfield , John Penny , in Numerical Methods (Fourth Edition), 2019

2.7 LU Decomposition

LU decomposition (or factorization) is a similar process to Gaussian elimination and is equivalent in terms of elementary row operations. The matrix A can be decomposed so that

(2.14) $\mathbf{A}=\mathbf{LU}$

where L is a lower triangular matrix with a leading diagonal of ones and U is an upper triangular matrix. Matrix A may be real or complex. Compared with Gaussian elimination, LU decomposition has a particular advantage when the equation system we wish to solve, $\mathbf{Ax}=\mathbf{b}$ , has more than one right side or when the right sides are not known in advance. This is because the factors L and U are obtained explicitly and they can be used for any right sides as they arise without recalculating L and U. Gaussian elimination does not determine L explicitly but rather forms ${\mathbf{L}}^{-1}\mathbf{b}$ so that all right sides must be known when the equation is solved.

The major steps required to solve an equation system by LU decomposition are as follows. Since $\mathbf{A}=\mathbf{LU}$ , then $\mathbf{Ax}=\mathbf{b}$ becomes

$\mathbf{LUx}=\mathbf{b}$

where b is not restricted to a single column. Letting $\mathbf{y}=\mathbf{Ux}$ leads to

$\mathbf{Ly}=\mathbf{b}$

Because L is a lower triangular matrix this equation is solved efficiently by forward substitution. To find x we then solve

$\mathbf{Ux}=\mathbf{y}$

Because U is an upper triangular matrix, this equation can also be solved efficiently by back substitution.

We now illustrate the LU decomposition process by solving (2.10) with $p=1$ . We are not concerned with b and we do not form an augmented matrix. We proceed exactly as with Gaussian elimination, see Table 2.1, except that we keep a record of the elementary row operations performed at the ith stage in ${\mathbf{T}}^{(i)}$ and place the results of these operations in a matrix ${\mathbf{U}}^{(i)}$ rather than over-writing A.

We begin with the matrix

(2.15) $\mathbf{A}=\left[\begin{array}{ccc}\hfill 3& \hfill 6& \hfill 9\\ \hfill 2& \hfill 5& \hfill 2\\ \hfill -3& \hfill -4& \hfill -11\end{array}\right]$

Following the same operations as used in Table 2.1, we will create a matrix ${\mathbf{U}}^{(1)}$ with zeros below the leading diagonal in the first column using the following elementary row operations:

(2.16) $\text{row 2 of}{\mathbf{U}}^{\left(1\right)}=\text{ row 2 of}\mathbf{A}-2(\text{row 1 of}\mathbf{A})/3$

and

(2.17) $\text{row 3 of}{\mathbf{U}}^{\left(1\right)}=\text{ row 3 of}\mathbf{A}+3(\text{row 1 of}\mathbf{A})/3$

Now A can be expressed as the product ${\mathbf{T}}^{(1)}$ ${\mathbf{U}}^{(1)}$ as follows:

$\left[\begin{array}{ccc}\hfill 3& \hfill 6& \hfill 9\\ \hfill 2& \hfill 5& \hfill 2\\ \hfill -3& \hfill -4& \hfill -11\end{array}\right]=\left[\begin{array}{ccc}\hfill 1& \hfill 0& \hfill 0\\ \hfill 2/3& \hfill 1& \hfill 0\\ \hfill -1& \hfill 0& \hfill 1\end{array}\right]\left[\begin{array}{ccc}\hfill 3& \hfill 6& \hfill 9\\ \hfill 0& \hfill 1& \hfill -4\\ \hfill 0& \hfill 2& \hfill -2\end{array}\right]$

Note that row 1 of A and row 1 of ${\mathbf{U}}^{(1)}$ are identical. Thus row 1 of ${\mathbf{T}}^{(1)}$ has a unit entry in column 1 and zero elsewhere. The remaining rows of ${\mathbf{T}}^{(1)}$ are determined from (2.16) and (2.17). For example, row 2 of ${\mathbf{T}}^{(1)}$ is derived by rearranging (2.16); thus:

(2.18) $\text{row 2 of}\mathbf{A}=\text{ row 2 of}{\mathbf{U}}^{(\mathbf{1})}-2(\text{row 1 of}\mathbf{A})/3$

or

(2.19) $\text{row 2 of}\mathbf{A}=2(\text{row 1 of}{\mathbf{U}}^{\left(1\right)})/3+\text{row 2 of}{\mathbf{U}}^{\left(1\right)}$

since row 1 of ${\mathbf{U}}^{(1)}$ is identical to row 1 of A. Hence row 2 of ${\mathbf{T}}^{(1)}$ is $[2/310]$ .

We now move to the next stage of the decomposition process. In order to bring the largest element of column 2 in ${\mathbf{U}}^{(1)}$ onto the leading diagonal we must interchange rows 2 and 3. Thus ${\mathbf{U}}^{(1)}$ becomes the product ${\mathbf{T}}^{(2)}{\mathbf{U}}^{(2)}$ as follows:

$\left[\begin{array}{ccc}\hfill 3& \hfill 6& \hfill 9\\ \hfill 0& \hfill 1& \hfill -4\\ \hfill 0& \hfill 2& \hfill -2\end{array}\right]=\left[\begin{array}{ccc}\hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 1\\ \hfill 0& \hfill 1& \hfill 0\end{array}\right]\left[\begin{array}{ccc}\hfill 3& \hfill 6& \hfill 9\\ \hfill 0& \hfill 2& \hfill -2\\ \hfill 0& \hfill 1& \hfill -4\end{array}\right]$

Finally, to complete the process of obtaining an upper triangular matrix we make

$\text{row 3 of}\mathbf{U}=\text{ row 3 of}{\mathbf{U}}^{\left(2\right)}-(\text{row 2 of}{\mathbf{U}}^{\left(2\right)})/2$

Hence, ${\mathbf{U}}^{(2)}$ becomes the product ${\mathbf{T}}^{(3)}\mathbf{U}$ as follows:

$\left[\begin{array}{ccc}\hfill 3& \hfill 6& \hfill 9\\ \hfill 0& \hfill 2& \hfill -2\\ \hfill 0& \hfill 1& \hfill -4\end{array}\right]=\left[\begin{array}{ccc}\hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 1/2& \hfill 1\end{array}\right]\left[\begin{array}{ccc}\hfill 3& \hfill 6& \hfill 9\\ \hfill 0& \hfill 2& \hfill -2\\ \hfill 0& \hfill 0& \hfill -3\end{array}\right]$

Thus $\mathbf{A}={\mathbf{T}}^{(1)}{\mathbf{T}}^{(2)}{\mathbf{T}}^{(3)}\mathbf{U}$ , implying that $\mathbf{L}={\mathbf{T}}^{(1)}{\mathbf{T}}^{(2)}{\mathbf{T}}^{(3)}$ as follows:

$\left[\begin{array}{ccc}\hfill 1& \hfill 0& \hfill 0\\ \hfill 2/3& \hfill 1& \hfill 0\\ \hfill -1& \hfill 0& \hfill 1\end{array}\right]\left[\begin{array}{ccc}\hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 1\\ \hfill 0& \hfill 1& \hfill 0\end{array}\right]\left[\begin{array}{ccc}\hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 1/2& \hfill 1\end{array}\right]=\left[\begin{array}{ccc}\hfill 1& \hfill 0& \hfill 0\\ \hfill 2/3& \hfill 1/2& \hfill 1\\ \hfill -1& \hfill 1& \hfill 0\end{array}\right]$

Note that owing to the row interchanges L is not strictly a lower triangular matrix but it can be made so by interchanging rows.

Matlab implements LU factorization by using the function lu and may produce a matrix that is not strictly a lower triangular matrix. However, a permutation matrix P may be produced, if required, such that $\mathbf{LU}=\mathbf{PA}$ with L lower triangular.

We now show how the Matlab function lu solves the example based on the matrix given in (2.15):

>> A = [3 6 9;2 5 2;-3 -4 -11]

A =

3   6   9

2   5   2

-3   -4   -11

To obtain the L and U matrices, we must use that Matlab facility of assigning two parameters simultaneously as follows:

>> [L1 U] = lu(A)

L1 =

1.0000   0   0

0.6667   0.5000   1.0000

-1.0000   1.0000   0

U =

3.0000   6.0000   9.0000

0   2.0000   -2.0000

0   0   -3.0000

Note that the L1 matrix is not in lower triangular form, although its true form can easily be deduced by interchanging rows 2 and 3 to form a triangle. To obtain a true lower triangular matrix we must assign three parameters as follows:

>> [L U P] = lu(A)

L =

1.0000   0   0

-1.0000   1.0000   0

0.6667   0.5000   1.0000

U =

3.0000   6.0000   9.0000

0   2.0000   -2.0000

0   0   -3.0000

P =

1   0   0

0   0   1

0   1   0

In the preceding output, P is the permutation matrix such that L*U = P*A or P'*L*U = A. Thus P'*L is equal to L1.

The Matlab operator \ determines the solution of $\mathbf{Ax}=\mathbf{b}$ using LU factorization. As an example of an equation system with multiple right sides we solve $\mathbf{AX}=\mathbf{B}$ where

$\mathbf{A}=\left[\begin{array}{ccc}\hfill 3& \hfill 4& \hfill -5\\ \hfill 6& \hfill -3& \hfill 4\\ \hfill 8& \hfill 9& \hfill -2\end{array}\right]\text{ and}\mathbf{B}=\left[\begin{array}{cc}\hfill 1& \hfill 3\\ \hfill 9& \hfill 5\\ \hfill 9& \hfill 4\end{array}\right]$

Performing LU decomposition such that $\mathbf{LU}=\mathbf{A}$ gives

$\mathbf{L}=\left[\begin{array}{ccc}\hfill 0.375& \hfill -0.064& \hfill 1\\ \hfill 0.750& \hfill 1& \hfill 0\\ \hfill 1& \hfill 0& \hfill 0\end{array}\right]\text{ and}\mathbf{U}=\left[\begin{array}{ccc}\hfill 8& \hfill 9& \hfill -2\\ \hfill 0& \hfill -9.75& \hfill 5.5\\ \hfill 0& \hfill 0& \hfill -3.897\end{array}\right]$

Thus $\mathbf{LY}=\mathbf{B}$ is given by

$\left[\begin{array}{ccc}\hfill 0.375& \hfill -0.064& \hfill 1\\ \hfill 0.750& \hfill 1& \hfill 0\\ \hfill 1& \hfill 0& \hfill 0\end{array}\right]\left[\begin{array}{cc}\hfill y_{11}& \hfill y_{12}\\ \hfill y_{21}& \hfill y_{22}\\ \hfill y_{31}& \hfill y_{32}\end{array}\right]=\left[\begin{array}{cc}\hfill 1& \hfill 3\\ \hfill 9& \hfill 5\\ \hfill 9& \hfill 4\end{array}\right]$

We note that implicitly we have two systems of equations which when separated can be written:

$\mathbf{L}\left[\begin{array}{c}\hfill y_{11}\\ \hfill y_{21}\\ \hfill y_{31}\end{array}\right]=\left[\begin{array}{c}\hfill 1\\ \hfill 9\\ \hfill 9\end{array}\right]\text{ and}\mathbf{L}\left[\begin{array}{c}\hfill y_{12}\\ \hfill y_{22}\\ \hfill y_{32}\end{array}\right]=\left[\begin{array}{c}\hfill 3\\ \hfill 5\\ \hfill 4\end{array}\right]$

In this example, L is not strictly a lower triangular matrix owing to the reordering of the rows. However, the solution of this equation is still found by forward substitution. For example, $1y_{11}=b_{31}=9$ , so that $y_{11}=9$ . Then $0.75y_{11}+1y_{21}=b_{21}=9$ . Hence $y_{21}=2.25$ , etc. The complete Y matrix is

$\mathbf{Y}=\left[\begin{array}{cc}\hfill 9.000& \hfill 4.000\\ \hfill 2.250& \hfill 2.000\\ \hfill -2.231& \hfill 1.628\end{array}\right]$

Finally solving $\mathbf{UX}=\mathbf{Y}$ by back substitution gives

$\mathbf{X}=\left[\begin{array}{cc}\hfill 1.165& \hfill 0.891\\ \hfill 0.092& \hfill -0.441\\ \hfill 0.572& \hfill -0.418\end{array}\right]$

The Matlab function det determines the determinant of a matrix using LU factorization as follows. Since $\mathbf{A}=\mathbf{LU}$ then $|\mathbf{A}|=|\mathbf{L}||\mathbf{U}|$ . The elements of the leading diagonal of L are all ones so that $|\mathbf{L}|=1$ . Since U is upper triangular, its determinant is the product of the elements of its leading diagonal. Thus, taking account of row interchanges the appropriately signed product of the diagonal elements of U gives the determinant.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128122563000117

Computer Solution of Large Linear Systems

In Studies in Mathematics and Its Applications, 1999

Proof

Suppose that there exists an LU factorization. From the proof of proposition 2.6, we have

$A_{k+1}=L_k\dots L_{1}A\cdot$

Therefore,

$A=L_1^{-1}\dots L_k^{-1}{A}_{k+1},$

and we have also

$A=LU\cdot$

In block form, this is written as:

$\begin{array}{ccc}A=\left(\begin{array}{cc}{A}_{1,1}& A_{1,2}\\ A_{2,1}& A_{2,2}\end{array}\right),& L=\left(\begin{array}{cc}{L}_{1,1}& 0\\ L_{2,1}& L_{2,2}\end{array}\right),& U=\left(\begin{array}{cc}{U}_{1,1}& U_{1,2}\\ 0& U_{2,2}\end{array}\right),\end{array}$

and from proposition 2.6,

$\begin{array}{cc}{L}_1^{-1}\dots L_k^{-1}=\left(\begin{array}{cc}{L}_{1,1}& 0\\ L_{2,1}& I\end{array}\right),& A_{k+1}=\left(\begin{array}{cc}{U}_{1,1}& U_{1,2}\\ 0& W_{2,2}\end{array}\right)\end{array},$

where all the matrices in block position (1, 1) are square of order k. By equating blocks, we have

$\begin{array}{l}{A}_{1,1}=L_{1,1}{U}_{1,1}\\ A_{2,2}=L_{2,1}{U}_{1,2}+L_{2,2}{U}_{2,2}\\ A_{2,2}=L_{2,1}{U}_{1,2}+W_{2,2}\end{array}$

Therefore, L _1,1 U _1,1 is the LU factorization of the leading principal submatrix of order k of A and L _2,2 U _2,2 is the factorization of the matrix W _2,2 in the bottom right hand corner of A _k+1 as we have W _2,2 = L _2,2 U _2,2

Note that det(A) = det(A _k+1). We have det(L _1,1) = 1 and det(A _1,1) = det(U _1,1). Since the matrix U _1,1 is upper triangular, its determinant is equal to the product of the diagonal elements. Therefore, for all k,

$det\left(A_{1,1}\right)=a_{1,1}^{\left(1\right)}\cdots a_{k,k}^{\left(k\right)}\cdot$

This shows that the principal minors are non–zero. The converse of the proof is easily derived by induction.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/S0168202499800034

Factorization of the Discretized Problem

Jacques Henry , Angel M. Ramos , in Factorization of Boundary Value Problems Using the Invariant Embedding Method, 2016

5.5 Row permutation

We now wish to consider the effect of block row permutation and show the connection with the continuous calculation in section 3.5. We adopt the conditions of the calculation in section 5.2. Consider the matrix A_h of problem $\left({\mathcal{P}}_h\right)$ . It is a well-known fact that the Gaussian method with row permutation is theoretically equivalent to a method that first performs all permutations and then performs LU factorization on the permuted matrix. In general, this is only relevant in a theoretical setting because at the k-th step the decision to perform a permutation is made based on the state of the factorization calculation at that step. However, in our case, we have block tridiagonal matrices and there is therefore only one single possible block row permutation at each step: the permutation that moves the sub-diagonal block into the diagonal position. Suppose we decide to perform this permutation systematically. In this case, the first row block ends up in last position, whereas all other blocks "move up" by one step:

[5.29] $A_h^{\prime }=\frac{1}{h^2}\left(\begin{array}{cccccc}\hfill -I\hfill & \hfill B_2\hfill & \hfill -I\hfill & \hfill \hfill & \hfill 0\hfill & \hfill \hfill \\ \hfill \hfill & \hfill \cdot \hfill & \hfill \cdot \hfill & \hfill \cdot \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \cdot \hfill & \hfill \cdot \hfill & \hfill \cdot \hfill & \hfill \hfill \\ \hfill \hfill & \hfill 0\hfill & \hfill \hfill & \hfill -I\hfill & \hfill \cdot \hfill & \hfill -I\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill -I\hfill & \hfill B_{N-1}\hfill \\ \hfill B_1\hfill & \hfill -I\hfill & \hfill \hfill & \hfill \hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right).$

The problem that we must solve then becomes:

[5.30] $A_h^{\prime }{u}_h=\left(\begin{array}{c}\hfill f_2\hfill \\ \hfill \cdot \hfill \\ \hfill f_{N-2}\hfill \\ \hfill f_{N-1}+\frac{u_1\left(a_N\right)}{h^2}\hfill \\ \hfill f_1+\frac{u_0\left(a_{1/2}\right)}{h}\hfill \end{array}\right).$

We apply the method of invariant embedding as in section 5.2. At step i ₀ ∈ {2, …, N}, we have:

$\frac{1}{h^2}\left(\begin{array}{cccccc}\hfill -I\hfill & \hfill B_2\hfill & \hfill -I\hfill & \hfill \hfill & \hfill 0\hfill & \hfill \hfill \\ \hfill \hfill & \hfill \cdot \hfill & \hfill \cdot \hfill & \hfill \cdot \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \cdot \hfill & \hfill \cdot \hfill & \hfill \cdot \hfill & \hfill \hfill \\ \hfill \hfill & \hfill 0\hfill & \hfill \hfill & \hfill -I\hfill & \hfill \cdot \hfill & \hfill -I\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill -I\hfill & \hfill B_{i_0-1}\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill -I\hfill \end{array}\right)\left(\begin{array}{c}\hfill u_1\hfill \\ \hfill \cdot \hfill \\ \hfill \cdot \hfill \\ \hfill \cdot \hfill \\ \hfill \cdot \hfill \\ \hfill u_{i_0-1}\hfill \end{array}\right)=\left(\begin{array}{c}\hfill f_2\hfill \\ \hfill \cdot \hfill \\ \hfill \cdot \hfill \\ \hfill f_{i_0-2}\hfill \\ \hfill f_{i_0-1}+\frac{{\gamma }_{i_0}}{h^2}\hfill \\ \hfill f_{i_0}+\frac{{\gamma }_{i_0+1}}{h^2}-\frac{1}{h^2}{B}_{i_0}{\gamma }_{i_0}\hfill \end{array}\right).$

The solution has an affine dependency on the data ${\gamma }_{i_0}$ and ${\gamma }_{i_0+1}$ . As the system is upper triangular with   − I on the diagonal, this solution can be obtained by simple successive substitution. However, as we want the solution u_h to satisfy the condition $\frac{u_1-u_0}{h}=-u_0\left(a_{1/2}\right)$ , the data cannot be independent but must satisfy an affine relation. We write this relation in the form:

[5.31] $\frac{{\gamma }_{i_0+1}-{\gamma }_{i_0}}{h}=P_{i_0+1}{\gamma }_{i_0+1}+w_{i_0+1/2}.$

We check that P_i and w _{i  +   1/2} are, respectively, solutions of [5.11] and [5.9]. For i ₀  =   0, regardless of the value of the data γ ₁, we want $\frac{{\gamma }_1-{\gamma }_0}{h}=-u_0\left(a_{1/2}\right)$ . We therefore deduce that:

[5.32] $P_1=0,w_{1/2}=-u_0\left(a_{1/2}\right).$

Consider the solution of [5.30] ${\left\{u_j\right\}}_{j=1}^{i_0-1}$ corresponding to the data ${\gamma }_{i_0},{\gamma }_{i_0+1}$ . We have:

[5.33] $\frac{{\gamma }_{i_0}-u_{i_0-1}\left({\gamma }_{i_0},{\gamma }_{i_0+1}\right)}{h}=P_{i_0}{\gamma }_{i_0}+w_{i_0-1/2},$

because ${\left\{u_j\right\}}_{j=1}^{i_0-2}$ is a solution of the same problem of dimension (i ₀  −   2) (p  −   1) for the data $u_{i_0-1}\left({\gamma }_{i_0},{\gamma }_{i_0+1}\right)$ and ${\gamma }_{i_0}$ . By subtracting [5.33] from [5.31]:

[5.34] $\begin{array}{lll}\frac{{\gamma }_{i_0+1}-2{\gamma }_{i_0}+u_{i_0-1}}{h^2}\hfill & =\hfill & -f_{i_0}-{\nabla }_{h,2}^2{\gamma }_{i_0}\hfill \\ \hfill & =\hfill & P_{i_0}\frac{{\gamma }_{i_0+1}-{\gamma }_{i_0}}{h}+\frac{P_{i_0+1}-P_{i_0}}{h}{\gamma }_{i_0}\hfill \\ \hfill & \hfill & +\frac{\left(w_{i_0+1/2}-w_{i_0-1/2}\right)}{h}.\hfill \end{array}$

By applying [5.31] to reduce to just ${\gamma }_{i_0}$ :

[5.35] $\frac{{\gamma }_{i_0+1}-{\gamma }_{i_0}}{h}={\left(I-h{P}_{i_0+1}\right)}^{-1}\left(P_{i_0+1}{\gamma }_{i_0}+w_{i_0+1/2}\right).$

By substituting this into [5.34] and using the fact that ${\gamma }_{i_0}$ is arbitrary, we have that:

[5.36] $-\frac{P_{i_0+1}-P_{i_0}}{h}=P_{i_0+1}{\left(I-h{P}_{i_0+1}\right)}^{-1}{P}_{i_0+1}+{\nabla }_{h,2}^2,$

and for the independent term ${\gamma }_{i_0}$ , we have:

[5.37] $\frac{w_{i_0+1/2}-w_{i_0-1/2}}{h}=-P_{i_0+1}{\left(I-h{P}_{i_0+1}\right)}^{-1}{w}_{i_0+1/2}-f_{i_0}.$

These equations are satisfied for i ₀ ranging from 2 to N  −   1. From the initial conditions[5.32], we also recover [5.11] and [5.9]. However, the "upward" phase for calculating u_i is now accomplished by solving [5.30] by substitution, which we can write as:

[5.38] $-u_{i-1}+2u_i-h^2{\nabla }_{h,2}^2{u}_i-u_{i+1}=h^2{f}_i\forall i\in \left\{2,\dots ,N-1\right\},$

adding the following equation to find u _{N −   1}:

[5.39] $\frac{u_1\left(a_N\right)-u_{N-1}}{h}=P_N{u}_1\left(a_N\right)+w_{N-1/2}.$

This calculation shows that the block LU factorization with systematic block row permutation may be interpreted as a discretization of the method of invariant embedding in the form of a family of Cauchy problems that was presented in section 3.5.

Remark 5.4

It is a well-known fact that the stability of the Gaussian method is improved by performing the row permutation that brings the sub-diagonal element with the highest absolute value into the pivotal position. When discretizing the second derivatives, as is the case here, the blocks with the highest spectral radius naturally end up in the diagonal position: indeed, in our case, if λ _h,i are the eigenvalues of   −   ∇ _h,2 ², the diagonal term of A_h presented in [5.1] has the eigenvalues 2(1   + hλ _h,i ), whereas the sub-diagonal blocks have the eigenvalue −   1.

There is therefore no benefit in performing block row permutations from the point of view of stability. The above calculation even shows that, in our case, the problem with row permutations leads to an ill-posed problem in the limit h  →   0: the Cauchy problem for the Laplacian. This shows the connection between the destabilization of the Gaussian factorization with pivoting and the instability of the Cauchy problem. However, as commented in section 3.5, in the formulation corresponding to the Cauchy problem, storing the matrices P_i is not necessary to calculate u_i .

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9781785481437500050

Computational and Numerical Challenges in Environmental Modeling

Giorgio Semenza , in Studies in Computational Mathematics, 2006

6.5.3 Using preconditioning

As mentioned in the previous subsection, dropping small non-zero elements can be useful in the attempts to preserve better the sparsity of the original matrix and, thus, to improve the performance of the selected sparse algorithm. While the performance is often improved when dropping of small non-zero elements is applied during the Gaussian elimination, the accuracy of the calculated triangular matrices L and U will often become very poor. This is why a preconditioning is necessary when this approach is used.

Ten preconditioned methods

Ten well-known numerical methods for solving systems of linear algebraic equations by using a direct method based on the Gaussian elimination and preconditioned methods for sparse matrices are discussed below. The methods are listed in Table 6.1. Some information about the abbreviations, as well as references, where more details about the ten methods can be found, are given in Table 6.1.

Table 6.1. Solvers for systems of linear algebraic equations with general sparse coefficient matrices, which are used in this section.

No.	Method	Abbreviation	Reference
1	Direct solution	DIR	[272], [281]
2	Iterative Refinement	IR	[274], [281]
3	Pure Modified Orthomin	PURE	[262], [281]
4	Preconditioned Modified Orthomin	ORTH	[262], [281]
5	Conjugate Gradients Squared	CGS	[230]
6	Bi-Conjugate Gradients STAB	BiCGSTAB	[255]
7	Transpose-Free Quasi Min. Residual	TFQMR	[107]
8	Generalized Minimum Residual	GMRES	[207]
9	Eirola-Nevanlinna Method	EN	[86], [263]
10	Block Eirola-Nevanlinna Method	BEN	[270]

The method for calculating directly the solution (DIR) is based on the use of a sparse matrix technique for calculation of an LU factorization of the coefficient matrix of the system of linear algebraic equations without dropping small non-zero elements, as in (Zlatev [281]).

The same method is also used to calculate approximate LU factorizations, which are used as preconditioners in the iterative methods. Let us reiterate here that a special parameter RELTOL (relative drop-tolerance), 0 ≤ RELTOL < 1, is used to calculate an approximate LU factorization. The use of this parameter to drop small non-zero elements during the Gaussian elimination was described in one of the previous subsections.

It should be added here that the relative drop-tolerance can also be used to drop all small non-zero elements in the original matrix before the beginning of the Gaussian elimination. The situation is very similar to that described in the previous subsection. If the absolute value of an element a_ij of the original coefficient matrix is smaller than the product of RELTOL and the largest in absolute value in row i, then a_ij is dropped (replaced by zero). This means that a_ij dropped when the following inequality

(6.7) $\left|a_{i,j}\right|\le \mathit{RELTOL}\underset{1\le j\le N}{max}\left|a_{i,j}\right|$

is satisfied. Note that (6.7) is very similar to (6.4), but while (6.7) is applied to the elements of the original matrix, (6.4) is applied to the elements calculated at some stage s + 1 of the Gaussian elimination.

If RELTOL is set equal to zero, then the direct solution can be calculated (see also the previous paragraph). However, if 0 < RELTOL < 1 then the LU factorization obtained by using the selected positive value of the relative drop-tolerance RELTOL must be used as a preconditioner in one of the iterative methods.

If RELTOL becomes larger, then in general the obtained LU factorization is becoming less accurate, but more sparse. This means that the number of iterations needed to obtain the required accuracy of the solution of the system of linear algebraic equations will tend to become greater for larger values of RELTOL, but the computational work per iteration will, in general, become smaller. If RELTOL is very large, then the rate of convergence can be very slow or the preconditioned iterative method may even not converge (the calculations must be repeated by using a smaller value of the relative drop-tolerance RELTOL when this happens; it might be necessary sometimes to reduce the value of RELTOL several times). This shows that it is very important to select robust and reliable stopping criteria when preconditioned iterative methods are used. The proper implementation of such stopping criteria will allow us to decide

•

whether the required accuracy is achieved or not (and to stop the iterative process when the required accuracy is achieved), and

•

whether the iterative process is slowly convergent or does not convergent at all and,

−

to reduce the relative drop-tolerance, and

−

to recalculate the LU factorization with the reduced drop-tolerance

when there are problems with the convergence.

It should be emphasized here that robust and reliable stopping criteria are also needed when pure iterative methods are applied.

Preconditioning is always used when the system of linear algebraic equations is solved by the last seven methods in Table 6.1. The Modified Orthomin Method (Zlatev [281]) is used both as a pure iterative method (PURE) and as a preconditioned iterative method (ORTH). If an approximate LU factorization is used in IR, then this method can be considered as a preconditioned iterative method.

It must be emphasized here that if the preconditioned iterative method used either is not convergent or the convergence is not sufficiently fast, then the preconditioner can automatically be improved, as described above, by reducing the relative drop-tolerance RELTOL (see more details in Zlatev [281]).

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/S1570579X06800212

Numerical Methods for Roots of Polynomials - Part II

J.M. McNamee , V.Y. Pan , in Studies in Computational Mathematics, 2013

10.4.3 Generalizations of the Quotient-Difference Method

Fox and Hayes (1968) associate a band matrix with the polynomial and perform an LU factorization, where L and U are triangular band matrices. They restrict their attention mainly to a quartic polynomial, but Kershaw (1987) extends their work to the general polynomial of degree N. Kershaw associates with $p(x)$ a semi-infinite band matrix

(10.117) $\mathbf{A}=\left[\begin{array}{cccccccccc}\hfill d_0^{(0)}\hfill & \hfill d_1^{(0)}\hfill & \hfill \dots \hfill & \hfill d_m^{(0)}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill \\ \hfill d_0^{(1)}\hfill & \hfill d_1^{(1)}\hfill & \hfill \dots \hfill & \hfill d_m^{(1)}\hfill & \hfill d_{m+1}^{(1)}\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill \\ \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \\ \hfill d_0^{(n-1)}\hfill & \hfill d_1^{(n-1)}\hfill & \hfill \dots \hfill & \hfill d_m^{(n-1)}\hfill & \hfill d_{m+1}^{(n-1)}\hfill & \hfill \dots \hfill & \hfill d_{N-1}^{(n-1)}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill \\ \hfill c_0\hfill & \hfill c_1\hfill & \hfill \dots \hfill & \hfill c_m\hfill & \hfill c_{m+1}\hfill & \hfill \dots \hfill & \hfill c_{N-1}\hfill & \hfill c_N\hfill & \hfill 0\hfill & \hfill \dots \hfill \\ \hfill 0\hfill & \hfill c_0\hfill & \hfill \dots \hfill & \hfill c_{m-1}\hfill & \hfill c_m\hfill & \hfill \dots \hfill & \hfill c_{N-2}\hfill & \hfill c_{N-1}\hfill & \hfill c_N\hfill & \hfill \dots \hfill \\ \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \end{array}\right]$

The first n rows are arbitrary, and from row n  +   1 onwards the element in the main diagonal is $c_n$ . We factorize A into LU where L is lower triangular with unit diagonals and band-width n  +   1, while U is upper triangular with band-width m   +   1. That is,

(10.118) $\mathbf{L}=\left[\begin{array}{ccccccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill \\ \hfill {\ell }_1^{(1)}\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill \\ \hfill {\ell }_2^{(2)}\hfill & \hfill {\ell }_1^{2)}\hfill & \hfill 1\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill \\ \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \\ \hfill {\ell }_n^{(n)}\hfill & \hfill {\ell }_{n-1}^{(n)}\hfill & \hfill {\ell }_{n-2}^{(n)}\hfill & \hfill \dots \hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill \dots \hfill \\ \hfill 0\hfill & \hfill {\ell }_n^{(n+1)}\hfill & \hfill {\ell }_{n-1}^{(n+1)}\hfill & \hfill \dots \hfill & \hfill {\ell }_1^{n+1)}\hfill & \hfill 1\hfill & \hfill \dots \hfill \\ \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \end{array}\right]$

while

(10.119) $\mathbf{U}=\left[\begin{array}{cccccc}\hfill u_0^{(0)}\hfill & \hfill u_1^{(0)}\hfill & \hfill \dots \hfill & \hfill u_m^{(0)}\hfill & \hfill 0\hfill & \hfill \dots \hfill \\ \hfill 0\hfill & \hfill u_0^{(1)}\hfill & \hfill \dots \hfill & \hfill u_{m-1}^{(1)}\hfill & \hfill u_m^{(1)}\hfill & \hfill \dots \hfill \\ \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \end{array}\right]$

We associate polynomials with the matrices A (top or arbitrary part), L, U as follows:

(10.120) $d_s(x)={\displaystyle \sum _{t=0}^{m+s}}{d}_t^{(s)}{x}^t(s=0\text{,}1\text{,}\dots \text{,}n-1)$

(10.121) ${\ell }_r(x)={\displaystyle \sum _{t=0}^n}{\ell }_t^{(r)}{x}^{n-t}(r=n\text{,}n+1\text{,}\dots )$

(10.122) $u_r(x)={\displaystyle \sum _{t=0}^m}{u}_t^{(r)}{x}^t(r=0\text{,}1\text{,}\dots )$

The part of A below the nth row is associated with $p(x)$ . Kershaw shows that usually ${\ell }_r(x)$ and $u_r(x)$ tend to fixed polynomials $\ell (x)$ and $u(x)$ in such a way that

(10.123) $p(x)=\ell (x)u(x)$

Of course the case n  =   1 gives a linear factor, and n  =   2 a quadratic for $\ell (x)$ . Or (although Kershaw does not mention this) we may use $n\approx \lfloor N/2\rfloor$ and repeat the process recursively until we have only linear and/or quadratic factors.

Jones and Magnus (1980) discuss what they call the F–G relations, which are very close to the Q–D scheme. Assuming that the zeros ${\zeta }_i$ of $p(x)$ satisfy

(10.124) $0<|{\zeta }_1|⩽|{\zeta }_2|⩽\dots ⩽|{\zeta }_n|$

Let

(10.125) $F_1^{(m)}=F_{n+1}^{(m)}=G_{n+1}^{(m)}=0(m=0\text{,}1\text{,}2\text{,}\dots )$

(10.126) $F_k^{(0)}=\frac{c_{n-k}}{c_{n-k+1}}(k=2\text{,}3\text{,}\dots \text{,}n)$

(10.127) $G_k^{(0)}=-\frac{c_{n-k}}{c_{n-k+1}}(k=1\text{,}2\text{,}\dots \text{,}n)$

Then for m  =   0,1,… let

(10.128) $G_1^{(m+1)}=G_1^{(m)}+F_2^{(m)}$

and

(10.129) $F_k^{(m+1)}=\frac{F_k^{(m)}(G_k^{(m)}+F_{k+1}^{(m)})}{G_{k-1}^{(m)}+F_k^{(m)}}(k=2\text{,}3\text{,}\dots \text{,}n)$

(10.130) $G_k^{(m+1)}=G_k^{(m)}+F_{k+1}^{(m)}-F_k^{(m+1)}(k=2\text{,}3\text{,}\dots \text{,}n)$

(Note: we assume that

(10.131) $G_{k-1}^{(m)}+F_k^{(m)}\ne 0\text{;}(k=2\text{,}3\text{,}\dots \text{,}n))$

The authors prove that if

(10.132) $|{\zeta }_{k-1}|<|{\zeta }_k|<|{\zeta }_{k+1}|$

Then

(10.133) ${\displaystyle \underset{m\to \infty }{\lim }}{G}_k^{(m)}={\zeta }_{n-k+1}$

(but we suspect that the above is a misprint for ${\zeta }_k$ on the right of (10.133)). In some numerical tests this method converged somewhat faster than the Q–D method, at least in some cases.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780444527301000047

Factorization by the QR Method

Jacques Henry , Angel M. Ramos , in Factorization of Boundary Value Problems Using the Invariant Embedding Method, 2016

Abstract:

In Chapter 5, we saw that in the context of a method of finite differences, one particular discretization of the factorization of problem $\left({\mathcal{P}}_0\right)$ corresponds exactly to the block Gaussian LU factorization of the matrix of the discretized problem $\left({\mathcal{P}}_0\right)$ . In other words, this factorization may be viewed as an infinite-dimensional generalization of the Gaussian factorization. It is natural to ask ourselves if there also exists a generalization of the QR method, where Q is an orthogonal matrix and R is an upper triangular matrix. We know that the QR method applies not only to square matrices, but also to rectangular matrices with more rows than columns and so corresponds to overdetermined problems that we can solve in the sense of least squares. This generalization to overdetermined problems is also possible for the factorization of the boundary value problems studied here (this is discussed in the doctoral thesis of Maria Orey and an additional boundary condition is discussed), but we will restrict our attention in this chapter to well-posed problems.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9781785481437500086

The Passive Dendritic Tree

FABRIZIO GABBIANI , STEVEN J. COX , in Mathematics for Neuroscientists, 2010

8.1 THE DISCRETE PASSIVE TREE

We work in the concrete context of Figure 8.2 on the way to a more general understanding. We have indexed the compartments, following an observation of Hines, in a manner that leads to minimal fill-in in the LU factorization, Exercise 5.2, of the resulting linear system associated with the backward Euler and trapezoid schemes. The physical lengths and radii of the three fibers are

FIGURE 8.2. The compartmentalization of a branched cell with soma, and its associated circuit diagram.

$\begin{array}{l}\begin{array}{lll}{\ell }_1,{\ell }_2,{\ell }_3,\hfill & \text{and}& a_1,\end{array}\hfill \end{array}{a}_2,a_3$

respectively, while the length of each compartment, except the soma, is dx. The soma is presumed to have surface area A _s and is not typically further compartmentalized.

If we inject I _stim at the soma then Kirchhoff's Current Law, at the node with potential v _3,4, requires

(8.1) $I_{stim}=C_m{A}_s{v}_{3,4}^{\prime }+g_{Cl}{A}_s{v}_{3,4}+a_3^2\pi \left(v_{3,4}-v_{3,3}\right)/\left(\text{d}x{R}_a\right)$

where, as in §6.1, v_i,ijj = V_i,ijj – V _Cl , while at the branch point (v _3,1) we find

$\frac{\pi a_3^2\left(v_{3,2}-v_{3,1}\right)}{R_a\text{dx}}=\frac{\pi a_2^2\left(v_{3,1}-v_{2,4}\right)}{R_a\text{dx}}+\frac{\pi a_1^2\left(v_{3,1}-v_{1,4}\right)}{R_a\text{dx}}+2\pi a_3\text{dx}\left(C_m{v}_{3,1}^{\prime }+g_{Cl}{v}_{3,1}\right).$

Current balance at the remaining nodes proceeds exactly as before, recall Eq. (6.5). In particular, with ${\lambda }_j^2\equiv a_j/\left(2R_a{g}_{Cl}\right)$ the squared space constant of branch j, we find

(8.2) $\begin{array}{rrr}\hfill \tau v_{1,1}^{\prime }+v_{1,1}-{\lambda }_1^2\left(v_{1,2}-v_{1,1}\right)/{\text{dx}}^2& \hfill =& \hfill \text{0}\\ \hfill \tau v_{1,1}^{\prime }+v_{1,2}-{\lambda }_1^2\left(v_{1,3}-2v_{1,2}+v_{1,1}\right)/{\text{dx}}^2& \hfill =& \hfill \text{0}\\ \hfill \tau v_{1,3}^{\prime }+v_{1,3}-{\lambda }_1^2\left(v_{1,4}-2v_{1,3}+v_{1,2}\right)/{\text{dx}}^2& \hfill =& \hfill \text{0}\\ \hfill \tau v_{1,4}^{\prime }+v_{1,4}-{\lambda }_1^2\left(v_{3,1}-2v_{1,4}+v_{1,3}\right)/{\text{dx}}^2& \hfill =& \hfill \text{0}\\ \hfill \tau v_{2,1}^{\prime }+v_{2,1}-{\lambda }_2^2\left(v_{2,2}-v_{2,1}\right)/{\text{dx}}^2& \hfill =& \hfill \text{0}\\ \hfill \tau v_{2,2}^{\prime }+v_{2,2}-{\lambda }_2^2\left(v_{2,3}-2v_{2,2}+v_{2,1}\right)/{\text{dx}}^2& \hfill =& \hfill \text{0}\\ \hfill \tau v_{2,3}^{\prime }+v_{2,3}-{\lambda }_2^2\left(v_{2,4}-2v_{2,3}+v_{2,2}\right)/{\text{dx}}^2& \hfill =& \hfill \text{0}\\ \hfill \tau v_{2,4}^{\prime }+v_{2,4}-{\lambda }_2^2\left(v_{3,1}-2v_{2,4}-v_{2,3}\right)/{\text{dx}}^2& \hfill =& \hfill \text{0}\\ \hfill \tau v_{3,1}^{\prime }+v_{3,1}+\frac{a_2{\lambda }_2^2(v_{3,1}-v_{2,4})-a_3{\lambda }_3^2\left(v_{3,2}-v_{3,1}\right)+a_1{\lambda }_1^2\left(v_{3,1}-v_{1,4}\right)}{a_3{\text{dx}}^2}& \hfill =& \hfill \text{0}\\ \hfill \tau v_{3,2}^{\prime }+v_{3,2}-{\lambda }_2^2\left(v_{3,3}-2v_{3,2}+v_{3,1}\right)/{\text{dx}}^2& \hfill =& \hfill \text{0}\\ \hfill \tau v_{3,3}^{\prime }+v_{3,3}-{\lambda }_3^2\left(v_{3,4}-2v_{3,3}+v_{3,2}\right)/{\text{dx}}^2& \hfill =& \hfill \text{0}\\ \hfill \tau v_{3,4}^{\prime }+v_{3,4}(A_3/A_s){\lambda }_3^2\left(v_{3,3}-v_{3,4}\right)/{\text{dx}}^2-I_{stim}/(g_{Cl}{A}_s)& \hfill =& \hfill 0\end{array}$

where A ₃ = 2πa ₃dx. We write this collection of equations as the linear system

(8.3) $\begin{array}{rr}\hfill v^{\prime }(t)=Bv(t)+f(t),& \hfill B=(H-I)/\tau ,f(t)=I_{stim}(t){\text{e}}_{12}/(C_m{A}_s)\end{array}$

and H is the Hines matrix

$\text{H}=\frac{1}{\text{d}{x}^2}\left(\begin{array}{llllllllllll}-{\lambda }_1^2\hfill & {\lambda }_1^2& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ {\lambda }_1^2& -2{\lambda }_1^2& {\lambda }_1^2& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& -{\lambda }_1^2& -2{\lambda }_1^2& -{\lambda }_1^2& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& -{\lambda }_1^2& -2{\lambda }_1^2& 0& 0& 0& 0& {\lambda }_1^2& 0& 0& 0\\ 0& 0& 0& 0& -{\lambda }_2^2& {\lambda }_2^2& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& {\lambda }_2^2& -2{\lambda }_2^2& {\lambda }_2^2& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& {\lambda }_2^2& -2{\lambda }_2^2& {\lambda }_2^2& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& {\lambda }_2^2& -2{\lambda }_2^2& {\lambda }_2^2& 0& 0& 0\\ 0& 0& 0& r_1{\lambda }_1^2& 0& 0& 0& r_2{\lambda }_2^2& -c& {\lambda }_3^2& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& {\lambda }_2^2& -2{\lambda }_3^2& {\lambda }_3^2& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& -2{\lambda }_3^2& {\lambda }_3^2\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& \rho {\lambda }_3^2& -\rho {\lambda }_3^2\end{array}\right)$

where

$\begin{array}{lllll}{r}_1=a_1/a_3,\hfill & r_2=a_2/a_3,& c=r_1{\lambda }_1^2+r_2{\lambda }_2^2+{\lambda }_3^2,& \text{and}& \rho =A_3/A_s.\end{array}$

The genius of H is at least double – it factors easily and is similar to a symmetric matrix. Regarding the former, we note that, as in the tridiagonal S of Chapter 6, Gaussian Elimination applied to this matrix requires only one elimination per column.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780123748829000083

Further Topics in Linear Programming

Bernard Kolman , Robert E. Beck , in Elementary Linear Programming with Applications (Second Edition), 1995

Algorithm

Virtually all LP codes designed for production, rather than teaching, use the revised simplex method. This method has several desirable features, including the ability to handle a large number of variables. The real limit on the size of a problem is the number of constraints (see Section 3.5). Other features will be described when we discuss error detection and correction.

Most of the large LP codes provide an option for computing B ⁻¹ that is based upon a procedure from numerical linear algebra called LU factorization . That is, B is written as LU, the product of a lower triangular matrix L and an upper triangular matrix U. The inverses of upper and lower triangular matrices are easily calculated. Then B ⁻¹ = U ⁻¹ L ⁻¹. Many large linear programming models have sparse matrices (ones with few nonzero entries). The matrix representations can then be highly compressed and L ⁻¹ and U ⁻¹ can be calculated in RAM, with special routines for sparse matrices, resulting in significant time savings. For this reason, more and more codes will provide an LU-factorization option.

The revised simplex algorithm with iterative B ⁻¹ calculation is usually programmed to check itself at specified intervals. Between checks it follows the description we gave in Section 3.4. The check involves computing the next B ⁻¹ in a manner different from the one we described. The matrix B can be constructed from the list of basic variables and the original problem as it was read in and stored. Then a very good method of numerically inverting B, such as the LU-factorization method described above, is used. This procedure of occasionally recomputing B ⁻¹ from the given problem serves to produce a more accurate basic feasible solution. However, in general the procedure is expensive in terms of computation time and must be used sparingly.

As was indicated in Section 2.2 most LP codes provide several options for handling degeneracy when it occurs.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780124179103500061

Inherently Parallel Algorithms in Feasibility and Optimization and their Applications

Y. Saad , in Studies in Computational Mathematics, 2001

4 PRECONDITIONERS

We now turn our attention to preconditioning techniques for distributed sparse matrices. Many of the ideas used in preconditioning distributed sparse matrices are borrowed from the Domain Decomposition literature; see, for example[15,36,57]. The simplest of these use block preconditionings based on the domains, and are termed Schwarz alternating procedures in the Partial Differential Equations literature. Let R_i be a restriction operator which projects a global vector onto its restriction in subdomain i. Algebraically, this is represented by an n_i   × n matrix, consisting of zeros and ones, where n_i is the size of subdomain i. The transpose of R_i represents a prolongation operator from subdomain i to the whole space. Let A_i be the local submatrix associated with subdomain i. In the Block Jacobi (or additive Schwarz) procedure, the global solution is updated as,

where r is the residual vector b  − Ax. Each of the solves with the matrices A_i is done independently. The method does not specify which technique is used for solving these local systems. It clear that a direct solution method could be used in case each subproblem is small enough.

However, it is common to use iterative solvers for the local systems since these are often fairly large. An iterative solver consists of an accelerator such as the conjugate gradient algorithm, and a preconditioner such as an Incomplete LU (ILU) factorization. For details on such techniques see, for example, [2,5,18,49]. In terms of software, a publically available package called SPARSKIT [46] includes a fair number of standard accelerators and preconditioners. Among the preconditioners available in SPARSKIT, the ILU with threshold (ILUT) reaches a good compromise between simplicity and robustness and is widely used.

In the of domain decomposition methods, however, the accelerator for the global iteration must take into account the fact that the preconditioner may not a constant operator (i.e., that it may change at each step of the outer iteration). Many variations are possible. In general, we found that iterating to solve each of the sub problems accurately is not cost-effective. Instead often a simple forward-backward sweep, with ILU factors obtained from an ILUT preconditioner, often yields the fastest combination. Subdomain partitions may be allowed to overlap. This technique works reasonably well for a small number of overlapping subdomains as was shown in experiments using a purely algebraic form in [44].

This can be extended to block Multiplicative Schwarz procedure, which is a block- Gauss-Seidel iteration in which, likewise, a block is associated with a domain. Here, the algorithm can be described by a simple loop which updates the local solutions for subdomain 1, 2, …, s in turn, taking into account the latest value reached for the current solution each time. The only difference with the block Jacobi iteration is that the solution and residual are updated immediately after each local correction δ_i . In the block Jacobi case, these local corrections are all computed using the same residual r and then they are added at the same time to the current solution x.

A technique of labeling the domains, known as "multicoloring", has been be exploited to extract parallelism from the Gauss-Seidel procedure, see e.g., [57]. In this procedure, the subdomains are colored such that subdomains that are coupled by an equation are always assigned a different color. One such coloring is illustrated in Figure 5 where the numbers refer to a color. With this coloring, the Gauss-Seidel loop can now proceed differently. Since couplings between vertices of different domains are only possible when these domains have different colors, the Gauss-Seidel sweep can proceed by color. Thus, in Figure 5 all 4 domains of color number 1 can update their variables first at the same time. Then the two domains with color number 2 can simultaneously do their update, and so on. For this example this color-based sweep will require 4 steps instead of s  =   12 steps for the domain-based sweep. The color-based loop is as follows.

Figure 5. Subdomain coloring for the multicolor Multiplicative Schwarz procedure

For color  =   1,…, num_colors Do:

If mycolor   = color Do:

x :=x  + δ

r :=r-Aδ

EndIf

EndDo

Here, mycolor represents the color of the node to which the subdomain is assigned. Though this procedure represents an improvement over the standard Gauss-Seidel procedure, it is still sequential with respect to the colors since when one color is active (one of the steps of the color loop), all processors of a different color will be idle. This may lead to substantial loss of efficiency. Several remedies have been proposed, the most common one of which is to further partition each subdomain in each processor and then find a global coloring of the new, finer, partitioning. The goal here is to have all colors of this new partitioning represented in each processor. An example of a second level partitioning of the domain in Figure 5 is shown in Figure 6.

Figure 6. Two-level multicoloring of the subdomains for the multicolor Multiplicative Schwarz procedure

This two-level coloring now presents a different disadvantage, namely that the rate of convergence of the new procedure is usually smaller than that of the original scheme, due to the fact that the subdomains are smaller. Other remedies have been proposed, see e.g., [52] for a discussion. Our own experience is that many of the remedies do not come for free in that they tend to increase the cost of the scheme either due to the increased number of iterations or other overhead.

Distributed versions of the Incomplete LU factorization ILU(0) have been developed and tested on the Connection Machine 5 as early as 1993 [38]. More recently, the idea was generalized to higher order ILU(k) factorizations by Karypis and Kumar [32] and Hysom and Pothen [26]. While the original ILU factorization is sequential in nature, it is fairly easy to devise parallel versions by what amounts to reordering on the global system. The effect of such reorderings can be detrimental to convergence and this has been the subject of recent research, see, e.g., [6].

The strategy for unraveling parallelism in ILU is illustrated in Figure 7 in the simplest case of ILU(0). The unknowns associated with the interior nodes do not depend on the other variables. As a result they can be eliminated independently and in parallel in the ILU(0) process. The rows associated with the nodes on the interface of the subdomain will require more attention. Indeed, now the order in which to eliminate these rows is important, since this order will essentially determine the amount of parallelism. At one extreme we can have the rows eliminated one by one in the natural order 1,2,…,n which results in no parallelism. We can also try to reorder the equations in such a way as to maximize parallelism as is done in sparse Gaussian elimination - but this would be too costly for incomplete LU factorizations. A somewhat intermediate solution, is to eliminate the equations using a simple "schedule", i.e., a global order in which to eliminate the rows, set by a priority rule among processors. As an example, we could impose a global order based on the labels of the processors, meaning that between two adjacent processors, the one with a lower number has a higher priority in processing rows than the other. Thus, in our illustration of Figure 7 Processor 10 processes the rows associated with its interior nodes first, then it will wait for those interface rows belonging to processors 2, 4 and 6. Once they are received, by Processor 10, it will eliminate all its rows associated with (all) its interface nodes. Finally, once this is done, the interface rows in Processor 10 that are awaited by Processors 13 and 14 will be sent to them. This sets a global ordering for the elimination process. The local ordering of the rows, i.e., the order in which we process the interface rows in the same processor (e.g. Processor 10) may not be as important. The global order based on PE-number, defines a natural priority graph and parallelism can easily be exploited in a data-driven implementation.

Figure 7. Distributed ILU(0) with a priority rule based on PE numbers

However, it seems rather artificial to use the label of a processor to determine the global order of the elimination. In fact, we can also define a proper order in which to perform the elimination by replacing the PE-numbers by any labels provided any two neighboring processors have a different label. The most natural way of doing this is to perform a multicoloring of the subdomains, and use the colors in exactly the same way as above to define an order of the tasks. The colors can be defined to be positive integers as in the illustration of Figure 5.

It is possible to eliminate the interior variables u_i from each of the local systems (3). Doing so will leave us with a system which involves only interface variables (denoted by y_i (3)), and referred to as a Schur complement system. This system can be solved and then the other variables can be computed. One of the main difficulties with this approach is to precondition the Schur complement system. Indeed, depending on the formulation, this system may be dense or have diagonal blocks that are dense [49]. It is possible to define approximations for this Schur complement from which preconditioners can then be derived [8,7,13]. An alternative is to define preconditioners for the global system based on some approximate solution of the Schur complement system [53,51]. This approach does not require that the Schur complement system be solved exactly since the solution is only used to define a preconditioner. Generally speaking, Schur complement preconditioners are often harder to implement but of better qualily than Schwarz-based preconditioners.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/S1570579X01800252

Mesh Generation

Marshall Bern , Paul Plassmann , in Handbook of Computational Geometry, 2000

2.2 Solution methods

The solution of the sparse linear system is usually the most computationally demanding of the three steps. Solution methods include direct factorization and preconditioned iterative methods, methods which can vary dramatically in required storage and computational cost for different problems. Moreover, the discrete formulation and mesh generation steps greatly influence the efficacy of a solution method. For example, if one chooses to use higher-order basis functions in the finite element method, one can use a coarser mesh, and a smaller but denser linear system.

Direct factorization methods, such as sparse Cholesky or LU factorization, tend be slower but more robust than iterative methods. Luckily the extra computational cost can be amortized when the factorization is reused to solve for more than one right-hand side. An important issue in direct methods is ordering the vertices to minimize the "fill", the number of intermediate nonzeros. Nested dissection [ 72] uses graph separators to save an asymptotic factor of O(n ^3/d ) in the fill. Any planar graph admits separators of size O(n ^1/2); reasonable three-dimensional meshes admit separators of size O(n ^2/3) [89].

Iterative methods [17] have proved effective in solving the linear systems arising in physical modeling. Most large problems cannot be effectively solved without the use of preconditioning. A popular preconditioner involves an incomplete factorization. Rather than computing the exact factors for the matrix A  = LU, approximate factors are computed such that $A\approx \tilde{L}\tilde{U}$ and the preconditioned system

${\tilde{L}}^{-1}A{\tilde{U}}^{-1}\left(\tilde{U}u\right)={\tilde{L}}^{-1}f$

is solved iteratively. Ideally, the incomplete factors should be easy to compute and require a modest amount of storage, and the condition number of the preconditioned system should be much better than the original system.

Multigrid methods [140] can achieve the ultimate goal of iterative methods, convergence in O(1) iterations, for certain classes of problems. These methods use a sequence of meshes, graded from fine to coarse. The iterative solution of the fine linear system is accelerated by solving the coarser systems. The cycle repeatedly projects the current residual from a finer mesh onto the next coarser mesh and interpolates the solution from the coarser mesh onto the finer. One drawback of multigrid is that computing a sequence of meshes may be difficult for complicated geometries.

Domain decomposition [122] represents something of a hybrid between iterative and direct approaches. This approach divides the domain into possibly overlapping small domains; it typically solves subproblems on the small domains directly, but iterates to the global solution in which neighboring subproblem solutions agree. This approach enjoys some of the superior convergence properties of multigrid methods, while imposing less of a burden on the mesh generator. In fact, the domain is often partitioned so that subproblems admit structured meshes.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780444825377500073

Find Lu Factorizationof a and Complete Solution to Ax B

Source: https://www.sciencedirect.com/topics/mathematics/lu-factorization

Share :