18) Solving Systems

18) Solving Systems#

Last time#

  • Cholesky Decomposition

  • Profiling

Today#

  1. Solving Systenms

  2. Direct methods

  3. Iterative methods

  4. Example to solve a Partial Differential Equation (PDEs)

1. Solving systems#

Suppose we wish to solve a problem \(Ax=b\), where \(A, b\) are known. Each row of \(A\) and \(b\) represents a linear equation, so the problem represents a system of linear equations.

The problem \(Ax=b\) is easy in two particular cases:

  • When \(A\) is diagonal. In this case, \(x_j = b_j a_{j,j}^{-1}\).

  • When \(A\) is triangular (can solve by back-substitution).

Techniques to solve linear systems fall under two main categories: direct methods and iterative methods.

2. Direct Methods#

These methods are designed to solve the system exactly.

Gaussian elimination and LU decomposition#

  • Gaussian elimination, also known as row reduction is probably the most popular algorithm for solving systems of linear equations by hand. But nowadays, we rarely solve systems of linear equations by hand (especially if they exceed, say, three or four equations).

  • Row reduction produces a matrix decomposition of the original matrix. The elementary row operations may be viewed as the multiplication on the left of the original matrix by elementary matrices. Then the first part of the algorithm computes an LU decomposition, while the second part writes the original matrix as the product of a uniquely determined invertible matrix and a uniquely determined reduced row echelon matrix.

  • Efficiency: The number of arithmetic operations required to perform row reduction is one way of measuring the algorithm’s computational efficiency. For example, to solve a system of \(n\) equations for \(n\) unknowns by performing row operations on the matrix until it is in reduced row echelon form, and then solving for each unknown in reverse order, requires \(n(n + 1)/2\) divisions, \((2n3 + 3n2 − 5n)/6\) multiplications, and \((2n3 + 3n2 − 5n)/6\) subtractions, for a total of approximately \(2n^3/3\) operations. Thus it has a arithmetic complexity of \(O(n^3)\). The cost becomes prohibitive for systems with millions of equations (which are very common nowadays).

Other decompositions seen so far:#

  • QR factorization

  • Cholesky decomposition

3. Iterative Methods#

Because costs become prohibitive for direct solvers for systems with millions of equations (which are very common nowadays), large systems are generally solved using iterative methods. In iterative methods, a solution is approximated via repeated iterations, rather than exactly calculated. These methods need an initial value to generate a sequence of improving approximate solutions.

Among iterative methods, stationary iterative methods solve a linear system with an operator approximating the original one; and based on a measurement of the error in the result (the residual), form a “correction equation” for which this process is repeated until predefined termination criteria are met.

The basic iterative methods work by splitting the matrix \(A\) into

\[ A = M -N \]

where \(M\) is some “nice” matrix (easily invertible) related to \(A\).

To solve a linear system \(Ax=b\), once you have the decomposition \(A = M-N\), then you can solve:

\[ (M-N) x = b \, , \]
\[ M x = Nx + b \]

And the iterations are defined so that the new approximate solution, \(x^{k+1}\), is written in terms of the old approximate solution, \(x^k\), as:

\[ M x^{k+1} = Nx^k + b \, , \quad k \ge 0 \]

For a linear system linear system \(A x = b\) with exact solution \(x^{*}\), we define the error by

\[ e^{k}:= x^{*} - x^k \, , \quad k \ge 0 \]

Definition: An iterative method is called linear if there exists a matrix \(C \in \R^{n \times n}\) such that

\[ e^{k+1} = C e^k \, , \quad \forall k \ge 0 \]

And this matrix is called the iteration matrix.

Definition: An iterative method with a given iteration matrix \(C\) is called convergent if the following holds

\[ \lim_{k\rightarrow \infty} C^k = 0 \]

Theorem:

An iterative method with an iteration matrix \(C\) is convergent if and only if its spectral radius, \(\rho (C)\) is smaller than unity, that is:

\[ \rho(C) < 1 \]

where the spectral radius is defined as the maximum of the absolute values of the eigenvalues of a given matrix.

By splitting a matrix \(A\) into a “nice” (easily invertible) part \(M\) and the difference, \(N\), (i.e., \(A = M -N\)) we can see that for the iterates

\[ M x^{k+1} = Nx^k + b \, , \quad k \ge 0 \]

the iteration matrix is given by

\[ C = I - M^{-1} A = M^{-1} N \]

In fact, from the iterates we have:

\[ M x^{k+1} = Nx^k + b \]

Let \(e^k = x - x^k\), we then have

\[ M(x - x^{k+1}) = N (x - x^k) + b \]
\[ M e^{k+1} = N e^k +b \]
\[ e^{k+1} = \underbrace{M^{-1} N}_{C} e^k \]

Common splittings:#

Basic examples of stationary iterative methods use a splitting of \(A\) such as:

\[ A = D + L + U \]

with

\[ D = \textrm{diag}(A) \]

i.e., where \(D\) is only the diagonal part of \(A\), and \(L\) is the strict lower triangular part of \(A\). Respectively, \(U\) is the strict upper triangular part of \(A\).

Examples of splittings:

\[ M \equiv D = \textrm{diag}(A) \]

For this method, the iterative solution is then found as

\[ x^{k+1} = D^{-1} \left\{ b - (L+U) x^k \right\} \, , \quad k \ge 0 \]
\[ M = \omega^{-1} D =\omega^{-1} \textrm{diag}(A) \]

with a weight (\(\omega \ne 0\))

\[ M = \omega^{-1} I \]

with a weight \((\omega \ne 0)\)

\[ M = D + L \]
\[ M = \omega^{-1} D + L \]

for some weight, a.k.a., over-relaxation factor \((\omega \ne 0)\)

\[ M = \frac{1}{\omega (2-\omega)}( D+\omega L ) D^{-1} \left( D+\omega U \right) \]

for some weight, a.k.a., over-relaxation factor (\(\omega \notin \{0,2 \}\))

Linear stationary iterative methods are also called relaxation methods.

Krylov subspace methods:#

This class of methods are designed to solve a system with a suitable basis, comprised of the order-\(r\) Krylov subspace, generated by an n-by-n matrix \(A\) and a vector \(b\) of dimension \(n\). This is defined as the linear subspace spanned by the images of \(b\) under the first \(r\) powers of A (starting from \(A^0 = I\)), that is

\[ \textrm{span} \left \{b, Ab, A^2b, \dots, A^{r-1}b \right \} \]

Krylov subspace methods try to avoid matrix-matrix operations, but rather multiply vectors by the matrix and work with the resulting vectors. Starting with a vector \(b\), one computes \(Ab\), then one multiplies that vector by \(A\) again to find \(A^2 b\) and so on. All algorithms that work this way are referred to as Krylov subspace methods; they are among the most successful methods currently available in numerical linear algebra.

These methods can be used in situations where there is an algorithm to compute the matrix-vector multiplication without there being an explicit representation of \(A\), giving rise to matrix-free methods.

The Conjugate-Gradient method:#

Another very interesting method, that uses the residual form and can be seen both as a direct method and an iterative method is the Conjugate-Gradient method.

4. Example to solve a Partial Differential Equation#

A LOT of physical processes and engineering problems can be described via Partial Differential Equations (PDEs). PDEs model phenomenon in many fields with great societal impact such as medicine, geophysics, climate change, electrodynamics, finance, economics, weather forecasting, etc.

Approximate solutions for large systems of PDEs can only be found using numerical techniques.

Here we introduce an example of a very simple diffusion problem in Physics (heat diffusion), the heat (diffusion) equation in \(1\) dimension:

\[ \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} \]

Consider a \(1D\) domain \([0,L]\) (can represent a rod), where the temperature function \(u=u(x,t)\) varies in time and space.

We define uniformly equi-spaced points in the domain \(x_i = x_0 + i \Delta x \), with \(x_0=0\) and \(x_{m}=L\).

Then, at a given point in space and time, we approximate our solution by

\[ u^n_i \sim u(t=t^n, x=x_i) = u(t_0 + n \Delta t, x_0 + i \Delta x) \]

This define our solution vector, where we consider Dirichelet (homogeneous) boundary conditions: \(u_0 = u(x_0) = 0\) and \(u_m = u(x_m) = 0\), for all times \(t\)

\[\begin{split}\mathbf{u} = \begin{bmatrix} u_0 = 0 \\ u_1 \\ u_2 \\ \vdots \\ u_{m-1} \\ u_m = 0 \end{bmatrix} \end{split}\]

Finite Differences#

To define derivatives, we use Taylor’s expansion:

\[ u_{i+1} = u(x_i + \Delta x) \approx u(x_i) + \Delta x \partial_x u \bigg\rvert_i + \frac{1}{2} \Delta x^2 \partial^2_x u \bigg\rvert_i + \frac{1}{6} \Delta x^3 \partial^3_x u \bigg\rvert_i + O(\Delta x^4) \]

Similarly,

\[ u_{i-1} = u(x_i - \Delta x) \approx u(x_i) - \Delta x \partial_x u \bigg\rvert_i + \frac{1}{2} \Delta x^2 \partial^2_x u \bigg\rvert_i - \frac{1}{6} \Delta x^3 \partial^3_x u \bigg\rvert_i + O(\Delta x^4) \]

We can define the first-order derivative at the point \(x_i\) using the forward difference:

\[ \frac{\partial u_i^F}{\partial x} \approx \frac{ u_{i+1} - u_{i}}{\Delta x} \]

Similarly, we can define the first-order derivative at the point \(x_i\) using the backward difference:

\[ \frac{\partial u_{i}^B}{\partial x} \approx \frac{ u_{i} - u_{i-1}}{\Delta x} \]

We can now define a second-order derivative, at the point \(x_i\) using a centered difference formula:

\[ \frac{\partial^2 u^C_i}{\partial x^2} = \frac{\frac{\partial u_{i+\frac{i}{2}}^F}{\partial x} -\frac{ \partial u_{i-\frac{i}{2}}^B}{\partial x} }{\Delta x} \approx \frac{u_{i-1} -2 u_i + u_{i+1}}{\Delta x^2} + O(\Delta x^2)\]

If we do this for all points \(x_i \in [0,L]\), we obtain a matrix representation of the second-order finite difference centered operator:

\[\begin{split} A = \frac{1}{\Delta x^2} \begin{bmatrix} -2 & 1 & 0 & & 0 \\ 1 & -2 & 1 & & \vdots \\ 0 & 1 & -2 & \\ & & & \ddots & 1 \\ 0 & \ldots & & 1 & -2\end{bmatrix} \end{split}\]

(We can factor out the \(\frac{1}{\Delta x^2}\) term in front and from now on we’ll consider the matrix of coefficients only)

Note

In Physics, it is customary to write the problem with a negative sign:

\[ \frac{\partial u}{\partial t} = -\frac{\partial^2 u}{\partial x^2} \]

This way, your matrix A is SPD.

  • To solve this with Jacobi’s method:

\[\begin{split} M = \begin{bmatrix} -2 & 0 & 0 & & 0 \\ 0 & -2 & 0 & & \vdots \\ 0 & 0 & -2 & \\ & & & \ddots & 0\\ 0 & \ldots & & 0 & -2 \end{bmatrix} \end{split}\]

and

\[\begin{split} N = \begin{bmatrix} 0 & -1 & 0 & & 0 \\ -1 & 0 & -1 & & \vdots \\ 0 & -1 & 0 & \\ & & & \ddots & -1\\ 0 & \ldots & & -1 & 0\end{bmatrix} \end{split}\]

Then,

\[\begin{split} M^{-1}N = \begin{bmatrix} 0 & 1/2 & 0 & & 0 \\ 1/2 & 0 & 1/2 & & \vdots \\ 0 & 1/2 & 0 & \\ & & & \ddots & 1/2 \\ 0 & \ldots & & 1/2 & 0\end{bmatrix} \end{split}\]

This has a spectral radius, \(\rho (M^{-1}N ) = \cos\left(\frac{\pi}{m+1} \right)\)

  • To solve this with Gauss-Seidel’s method:

\[\begin{split} M = \begin{bmatrix} -2 & 0 & 0 & & 0 \\ 1 & -2 & 0 & & \vdots \\ 0 & 1 & -2 & \\ & & & \ddots & 0\\ 0 & \ldots & & 1 & -2 \end{bmatrix} \end{split}\]

and

\[\begin{split} N = \begin{bmatrix} 0 & -1 & 0 & & 0 \\ 0 & 0 & -1 & & \vdots \\ 0 & 0 & 0 & \\ & & & \ddots & -1\\ 0 & \ldots & & 0 & 0 \end{bmatrix} \end{split}\]

And in this case the spectral radius, \(\rho (M^{-1}N ) = \cos^2\left(\frac{\pi}{m+1} \right)\)

  • To solve this with SOR’s method:

\[\begin{split} M = \begin{bmatrix} -2 & 0 & 0 & & 0 \\ \omega & -2 & 0 & & \vdots \\ 0 & \omega & -2 & \\ & & & \ddots & 0\\ 0 & \ldots & & \omega & -2 \end{bmatrix} \end{split}\]

and

\[\begin{split} N = \begin{bmatrix} (\omega-1) & -\omega & 0 & & 0 \\ 0 & (\omega-1) & -\omega & & \vdots \\ 0 & 0 & 0 & \\ & & & \ddots & -\omega\\ 0 & \ldots & & 0 & (\omega-1) \end{bmatrix} \end{split}\]

which has a spectral radius, \(\rho (M^{-1}N ) = \frac{ \cos^2\left(\frac{\pi}{m+1} \right)}{\left( 1+\sin\left( \frac{\pi}{m+1} \right) \right)^2}\)

Contents