Systems of linear equations - continued

LU decomposition¶

To solve a system of linear equations, Gaussian elimination requires the fewest arithmetic operations.

When the coefficient matrix $\mathbf{A}$ does not change and only the constant vector $\mathbf{b}$ changes, we can avoid repeating Gaussian elimination on the coefficient matrix. By decomposing the matrix $\mathbf{A}$ we can reach the solution with fewer arithmetic operations. To this end we will look at the LU decomposition!

Any matrix can be written as a product of two matrices:

Here, there are infinitely many ways to write the matrices $\mathbf{B}$ and $\mathbf{C}$.

In the LU decomposition we require that the matrix $\mathbf{B}$ is lower triangular and the matrix $\mathbf{C}$ is upper triangular:

Each of the matrices $\mathbf{L}$ and $\mathbf{U}$ has $(n+1)\,n/2$ nonzero elements; together that is $n^2+n$ unknowns. The known matrix $\mathbf{A}$ defines $n^2$ values. So to uniquely determine the matrices $\mathbf{L}$ and $\mathbf{U}$ we are missing $n$ equations. Here we will use the LU decomposition, which obtains the additional equations through the condition $L_{ii}=1$, $i=0, 1,\dots,n-1$.

The system of linear equations:

is therefore written using the decomposition of the matrix $\mathbf{A}$:

We now reach the solution of the system $\mathbf{A}\mathbf{x}=\mathbf{b}$ by solving two triangular systems of equations.

First we compute the vector $\mathbf{y}$:

Once $\mathbf{y}$ is computed, we can determine $\mathbf{x}$ from:

.

LU decomposition of the coefficient matrix $\mathbf{A}$¶

In what follows we will show that Gaussian elimination actually represents an LU decomposition of the coefficient matrix $\mathbf{A}$. We will rely on symbolic computation, so let’s import the sympy package:

Let’s begin the demonstration with an example of symbolically written matrices $\mathbf{L}$ and $\mathbf{U}$ of dimension $3\times 3$:

The coefficient matrix $\mathbf{A}$ written with the elements of the matrices $\mathbf{L}$ and $\mathbf{U}$ is therefore:

Let’s now perform Gaussian elimination on the coefficient matrix $\mathbf{A}$.

Using the first row, we perform Gaussian elimination in the first column:

We continue in the second column:

From the elimination above we observe:

1. the matrix $\mathbf{U}$ equals the matrix we obtain if we perform Gaussian elimination on the coefficient matrix $\mathbf{A}$.

• the off-diagonal terms of $\mathbf{L}$ are the factors we used in Gaussian elimination.

Numerical implementation of the LU decomposition¶

We will look at the numerical implementation on the system defined as:

Let’s perform Gaussian elimination and store the coefficients by which we multiply the pivot row m into the matrix L at the position with the indices that the eliminated element has in the matrix $\mathbf{A}$.

Let’s fill in the diagonal of $\mathbf{L}$:

Now let’s solve the lower triangular system of equations $\mathbf{L}\,\mathbf{y}=\mathbf{b}$:

We continue by solving the upper triangular system $\mathbf{U}\,\mathbf{x}=\mathbf{y}$:

As we noted above, when the constant vector $\mathbf{b}$ changes, repeating Gaussian elimination is not necessary. We only need to perform forward and then back substitution. Let’s look at an example:

Pivoting¶

Let’s look at the following system of equations:

If we were to perform Gaussian elimination in the first column of the matrix $\mathbf{A}$:

We notice that we have a problem with division by 0 in the first row. The elementary operations we can perform on the system allow swapping arbitrary rows. We can rearrange the system so that the pivot element is not equal to 0. Even so, it can happen that the pivot element we divide by has a very small value $\varepsilon$. Because this would increase the round-off error, among all the rows we choose for the pivot row the one whose pivot has the largest absolute value.

If during Gaussian elimination we swap rows so that the pivot element is the largest, we call this row pivoting or also partial pivoting. This makes Gaussian elimination numerically stable.

It can be shown that when solving a system of equations $\mathbf{A}\,\mathbf{x}=\mathbf{b}$ in which the matrix $\mathbf{A}$ is diagonally dominant, row pivoting is not necessary. Solving without pivoting is numerically stable.

A square matrix $\mathbf{A}$ of dimension $n$ is diagonally dominant if the absolute value of the diagonal element of each row is greater than the sum of the absolute values of the other elements in the row:

Gaussian elimination with partial pivoting¶

We will look at Gaussian elimination with partial pivoting. Without partial pivoting, in the $i$-th elimination step we choose row $i$ as the pivot. With partial pivoting we first check whether the $i$-th diagonal element is, in absolute value, the largest element in column $i$ on or below the diagonal; if not, we swap row $i$ with the row below it that has the largest element in absolute value in column $i$. With partial pivoting we reduce the influence of the round-off error on the result.

If we were to perform full pivoting, in addition to swapping rows we would also use swapping the order of the variables (column swapping). Full pivoting improves stability, but it is rarely used and we will not address it here.

The algorithm for Gaussian elimination with partial pivoting is therefore:

LU decomposition with partial pivoting¶

Similarly to Gaussian elimination, we can also extend the LU decomposition with partial pivoting. Solving thus becomes numerically stable. In doing so we must store the information about the row swaps, which we then pass to the function for solving the corresponding triangular systems.

Let’s look at an example:

In the example above we used a compact way of storing the triangular matrices $\mathbf{L}$ and $\mathbf{U}$; each is in fact defined by 6 elements, and for the matrix $\mathbf{L}$ we know that the diagonal elements equal 1. The matrix lu thus contains $3\times3=9$ elements:

On the diagonal and above the diagonal are the values of the upper triangular matrix $\mathbf{U}$, while below the diagonal are the subdiagonal elements of the matrix $\mathbf{L}$. When computing the solution we will take into account that the diagonal values of $\mathbf{L}$ equal 1.

The numeric list piv tells us how the rows were swapped, which must be taken into account when computing the solution:

Let’s now also determine the solution (the code for computing the solution is in the module tools.py)

Checking the solution:

The SciPy module¶

The SciPy module is built on the numpy module and contains many different high-level programs/modules/functions. The theoretical background of the modules is of course various numerical algorithms; we will get to know some of them within this textbook as well. A good source for these numerical algorithms in connection with SciPy is the documentation. For an excellent introduction to SciPy you can watch the YouTube video: SciPy Tutorial (2022): For Physicists, Engineers, and Mathematicians.

A brief overview of the module hierarchy:

• Linear algebra (scipy.linalg)

• Integration (scipy.integrate)

• Optimization (scipy.optimize)

• Interpolation (scipy.interpolate)

• Fourier transform (scipy.fft)

• Sparse matrices (scipy.sparse)

• Statistics (scipy.stats)

• Signal processing (scipy.signal)

• Special functions (scipy.special)

• Multidimensional image processing (scipy.ndimage)

• File handling (scipy.io)

In the following lectures we will look at some of the submodules in more detail.

Let’s look at how the LU decomposition is implemented within SciPy:

The function scipy.linalg.lu_factor (documentation):

lu_factor(a, overwrite_a=False, check_finite=True)

requires the input of a coefficient matrix (or a list of coefficient matrices) a; overwrite_a, if True, overwrites the value of a with the result of the decomposition (this can be important to save memory and increase speed). The function lu_factor returns a tuple (lu, piv):

• lu - the $\mathbf{L}$ \ $\mathbf{U}$ matrix (a matrix of the same dimension as a, but below the diagonal it contains the elements of $\mathbf{L}$, while the remaining elements define $\mathbf{U}$; the diagonal elements of $\mathbf{L}$ have values 1.

• piv - the pivot indices, representing a permutation matrix P: row i of the matrix a was swapped with row piv[i] (at each step the previous state is taken into account, slightly different logic than in our function LU_decomposition_pivoting).

lu_factor is used together with the function scipy.linalg.lu_solve (documentation), which gives us the solution of the system:

lu_solve(lu_and_piv, b, trans=0, overwrite_b=False, check_finite=True)

lu_and_piv is the tuple of the result (lu, piv) from lu_factor, b is the constant vector (or a list of vectors). The other parameters are optional.

Let’s look at how to use it:

Let’s obtain the solution:

Let’s check the validity of the solution:

Computing the inverse matrix¶

The inverse matrix of a square matrix $\textbf{A}$ of order $n\times n$ is denoted by $\textbf{A}^{-1}$. It is a matrix of order $n\times n$, such that

where $\mathbf{I}$ is the identity matrix.

The most efficient way to compute the inverse of a matrix $\mathbf{A}$ is to solve the matrix equation:

The matrix $\mathbf{X}$ is the inverse of the matrix $\mathbf{A}$: $\mathbf{A}^{-1}=\mathbf{X}$.

Computing the inverse matrix is therefore equal to solving $n$ systems of $n$ linear equations:

where $\mathbf{b}_i$ is the $i$-th column of the matrix $\mathbf{B}=\mathbf{I}$.

Computational cost: we perform the LU decomposition on the matrix $\mathbf{A}$ (computational cost of order $n^3$) and then find the solution for each $\mathbf{x}_i$ ($2\,n^2$ arithmetic operations for each $i$). The total computational cost is therefore still of order $n^3$. (If we were to perform Gaussian elimination $n$ times, the cost would be of order $n^4$.)

Here the point of the LU decomposition now becomes apparent; we compute the decomposition only once and then, for each constant vector, find the solution by solving two triangular systems. If we were to solve the systems by Gaussian elimination, we would need $(n^3+n^2)\,n$ arithmetic operations.

Let’s compute the inverse of the matrix $\textbf{A}$. First we prepare the identity matrix using np.identity():

Then we solve the system of equations for each column of the identity matrix:

The solution is therefore:

Let’s check the solution:

Solution using NumPy:

With appropriate use of functions from the SciPy module we can reach the solution even faster. Into the function lu_solve we can insert a vector $\mathbf{b}$ or a matrix $\mathbf{B}$, each column $i$ of which represents a new constant vector $\mathbf{b}_i$.

Solving overdetermined systems¶

When we solve a system of $m$ linear equations with $n$ unknowns and $m>n$ holds and the rank of the augmented matrix $[\mathbf{A}\,|\,\mathbf{b}]$ equals $n+1$, we have an overdetermined (inconsistent) system.

An overdetermined (also inconsistent system):

has no solution. However, we can find the best approximate solution using the least-squares method.

The sum of squares of the residuals is defined by the dot product:

which we reshape into:

where we have used the fact that, due to the scalar value, $\mathbf{b}^T\mathbf{A}\mathbf{x}=(\mathbf{b}^T(\mathbf{A}\mathbf{x}))^T=(\mathbf{A}\mathbf{x})^T\mathbf{b}$.

The solution of the equation, the gradient of the sum of squares, determines its minimum:

Thus, from the normal equation:

we determine the best approximate solution:

By introducing the pseudo-inverse matrix:

of the matrix $\mathbf{A}$, the solution of the overdetermined system is written:

The procedure presented above is relatively simple, but it is numerically demanding and in some cases can be ill-conditioned; in practice it is recommended that the pseudo-inverse matrix be computed using the functions:

• numpy.linalg.pinv from the numpy module (documentation),

• The functions pinv or pinvh from the scipy.linalg module (the choice depends on the problem being addressed; see the documentation),

which are based on better numerical methods.

Example of a system with a unique solution:

Now let’s make an overdetermined system (the function numpy.vstack (documentation) stacks lists by columns, numpy.random.seed (documentation) resets the random number generator to the seed value seed, while numpy.random.normal (documentation) generates a normally distributed list of length size and standard deviation scale):

Now let’s solve the overdetermined system:

We see that the overdetermined system with random values (simulation of measurement noise) gives a similar result to the solution without noise. If we were to increase the noise level, the deviation from the unique solution would increase.

We can also determine the pseudo-inverse matrix ourselves and check the difference with the built-in function:

Iterative methods¶

We often encounter large systems of linear equations whose coefficient matrix has few nonzero elements (such matrices are called sparse).

Iterative methods perform very well when solving such systems of linear equations; their advantages compared to direct methods are:

• arithmetic operations are performed only on the nonzero elements (despite the iterative solving, there can be fewer of them)

• the required memory space can be incomparably smaller.

Gauss-Seidel method¶

In what follows we will look at the idea of the Gauss-Seidel iterative method. First we write the system of equations $\mathbf{A}\,\mathbf{x}=\mathbf{b}$ as:

Suppose we are in step $k-1$ of the iterative method and the approximations $x_{j}^{(k-1)}$ ($j=0,1,\dots, n-1$) are known. From the sum above we isolate the term $i$:

Since the approximations $x_{j}^{(k-1)}$ do not exactly satisfy the linear problem, we can determine a new approximation $x_{i}^{(k)}$ from the equation above:

We have split the sum into two parts and, for computing the $i$-th term, taken into account the terms with index smaller than $i$ already determined in the $k$-th iteration.

We stop the iterative approach when we reach the desired accuracy of the solution $\epsilon$:

Example¶

Initial approximation:

Let’s prepare the matrix $\mathbf{A}$ without the diagonal elements (Why? Try to answer below, once we have performed the iterations.)

Let’s perform the iterations:

(Since we will write the values into x right away, there is no need to split the sum into two parts)

Let’s check the solution

The method works well if the matrix is diagonally dominant (there are, however, methods that work even when the matrix is not diagonally dominant; see e.g.: J. Petrišič, Reševanje enačb, 1996, p. 149: Conjugate gradient method).