4.5 Iterative methods for linear systems

Iterative methods are powerful tools for solving large, sparse linear systems in numerical analysis. They start with an initial guess and refine the solution iteratively, making them efficient for problems where direct methods are impractical.

These methods, including Jacobi, Gauss-Seidel, and conjugate gradient, offer advantages in memory usage and computational cost. They're particularly useful in data science and statistics when dealing with massive datasets or complex models.

Iterative methods overview

• Iterative methods are a class of algorithms used to solve large, sparse linear systems of equations in numerical analysis, data science, and statistics
• Unlike direct methods that aim to solve the system in a finite number of operations, iterative methods start with an initial guess and iteratively improve the solution until convergence is reached
• Iterative methods are particularly useful when the linear system is too large to be solved efficiently using direct methods like Gaussian elimination or LU decomposition

Comparison to direct methods

Top images from around the web for Comparison to direct methods

• 

  Expected Ideal Solution and Verified By Numerical Analysis for a Triangular Unit Cell Model View original

  Is this image relevant?

• 

  An Effective Numerical Calculation Method for Multi-Time-Scale Mathematical Models in Systems ... View original

  Is this image relevant?

• 

  Gauss-Legendre Iterative Methods and Their Applications on Nonlinear Systems and BVP-ODEs View original

  Is this image relevant?

• 

  Expected Ideal Solution and Verified By Numerical Analysis for a Triangular Unit Cell Model View original

  Is this image relevant?

• 

  An Effective Numerical Calculation Method for Multi-Time-Scale Mathematical Models in Systems ... View original

  Is this image relevant?

1 of 3

Top images from around the web for Comparison to direct methods

• 

  Expected Ideal Solution and Verified By Numerical Analysis for a Triangular Unit Cell Model View original

  Is this image relevant?

• 

  An Effective Numerical Calculation Method for Multi-Time-Scale Mathematical Models in Systems ... View original

  Is this image relevant?

• 

  Gauss-Legendre Iterative Methods and Their Applications on Nonlinear Systems and BVP-ODEs View original

  Is this image relevant?

• 

  Expected Ideal Solution and Verified By Numerical Analysis for a Triangular Unit Cell Model View original

  Is this image relevant?

• 

  An Effective Numerical Calculation Method for Multi-Time-Scale Mathematical Models in Systems ... View original

  Is this image relevant?

1 of 3

• Direct methods, such as Gaussian elimination and LU decomposition, solve linear systems exactly in a finite number of steps
• Iterative methods, on the other hand, start with an initial approximation and iteratively refine the solution until a desired level of accuracy is achieved
• Direct methods are generally more accurate but can be computationally expensive for large, sparse systems, while iterative methods are often more efficient in terms of memory and computational resources

Advantages of iterative methods

• Iterative methods are well-suited for large, sparse linear systems where direct methods become impractical due to high computational complexity and memory requirements
• They can exploit the sparsity of the matrix, performing operations only on non-zero elements, thus reducing storage and computational costs
• Iterative methods can be easily parallelized, making them suitable for distributed computing environments and high-performance computing

Disadvantages of iterative methods

• Convergence of iterative methods is not always guaranteed and depends on the properties of the linear system, such as the matrix condition number and spectral radius
• The rate of convergence can be slow for ill-conditioned systems, requiring a large number of iterations to achieve the desired accuracy
• The choice of initial guess and other parameters, such as relaxation factors or preconditioners, can significantly impact the convergence behavior and efficiency of iterative methods

Jacobi method

• The Jacobi method is one of the simplest iterative methods for solving linear systems of equations
• It is named after Carl Gustav Jacob Jacobi, a German mathematician who introduced the method in the 19th century
• The Jacobi method is based on the idea of updating each variable independently using the values of other variables from the previous iteration

Jacobi method algorithm

• Given a linear system $Ax = b$, where $A$ is an $n \times n$ matrix and $x$ and $b$ are $n$-dimensional vectors, the Jacobi method starts with an initial guess $x^{(0)}$

• At each iteration $k$, the method updates each component of the solution vector $x^{(k+1)}_i$ using the formula: $x^{(k+1)}_i = \frac{1}{a_{ii}} \left(b_i - \sum_{j \neq i} a_{ij} x^{(k)}_j \right)$

• The process is repeated until convergence is reached, i.e., when the difference between consecutive iterations falls below a specified tolerance

Jacobi method convergence

• The convergence of the Jacobi method depends on the properties of the matrix $A$
• The method is guaranteed to converge if the matrix $A$ is strictly diagonally dominant, meaning that for each row, the absolute value of the diagonal element is greater than the sum of the absolute values of the off-diagonal elements
• The rate of convergence is determined by the spectral radius of the iteration matrix, which is derived from the matrix $A$

Jacobi method example

• Consider the linear system: $3x_1 + x_2 = 5$ $x_1 + 2x_2 = 4$

• Starting with an initial guess $x^{(0)} = (0, 0)$, the Jacobi method updates the solution as follows: $x^{(1)}_1 = \frac{1}{3}(5 - x^{(0)}_2) = \frac{5}{3}$ $x^{(1)}_2 = \frac{1}{2}(4 - x^{(0)}_1) = 2$

  $x^{(2)}_1 = \frac{1}{3}(5 - x^{(1)}_2) = 1$ $x^{(2)}_2 = \frac{1}{2}(4 - x^{(1)}_1) = \frac{7}{6}$

• The process continues until the desired convergence is achieved, with the solution approaching $x = (1, 1.5)$

Gauss-Seidel method

• The Gauss-Seidel method is an iterative method for solving linear systems that improves upon the Jacobi method
• It is named after Carl Friedrich Gauss and Philipp Ludwig von Seidel, who developed the method in the 19th century
• The Gauss-Seidel method updates each variable using the most recently computed values of other variables, as opposed to the Jacobi method, which uses values from the previous iteration

Gauss-Seidel method algorithm

• Given a linear system $Ax = b$, where $A$ is an $n \times n$ matrix and $x$ and $b$ are $n$-dimensional vectors, the Gauss-Seidel method starts with an initial guess $x^{(0)}$

• At each iteration $k$, the method updates each component of the solution vector $x^{(k+1)}_i$ using the formula: $x^{(k+1)}_i = \frac{1}{a_{ii}} \left(b_i - \sum_{j < i} a_{ij} x^{(k+1)}_j - \sum_{j > i} a_{ij} x^{(k)}_j \right)$

• The process is repeated until convergence is reached, i.e., when the difference between consecutive iterations falls below a specified tolerance

Gauss-Seidel vs Jacobi method

• The Gauss-Seidel method generally converges faster than the Jacobi method because it uses the most recently computed values of variables in each iteration
• However, the Gauss-Seidel method is more sensitive to the ordering of the equations and variables, as the convergence rate can be affected by the order in which the variables are updated
• The Gauss-Seidel method requires sequential updates of variables, making it less suitable for parallel implementation compared to the Jacobi method

Gauss-Seidel method convergence

• Like the Jacobi method, the convergence of the Gauss-Seidel method depends on the properties of the matrix $A$
• The method is guaranteed to converge if the matrix $A$ is symmetric positive definite or strictly diagonally dominant
• The rate of convergence is determined by the spectral radius of the iteration matrix, which is derived from the matrix $A$

Gauss-Seidel method example

• Consider the same linear system as in the Jacobi method example: $3x_1 + x_2 = 5$ $x_1 + 2x_2 = 4$

• Starting with an initial guess $x^{(0)} = (0, 0)$, the Gauss-Seidel method updates the solution as follows: $x^{(1)}_1 = \frac{1}{3}(5 - x^{(0)}_2) = \frac{5}{3}$ $x^{(1)}_2 = \frac{1}{2}(4 - x^{(1)}_1) = \frac{7}{6}$

  $x^{(2)}_1 = \frac{1}{3}(5 - x^{(1)}_2) = \frac{17}{18}$ $x^{(2)}_2 = \frac{1}{2}(4 - x^{(2)}_1) = \frac{35}{36}$

• The process continues until the desired convergence is achieved, with the solution approaching $x = (1, 1.5)$ faster than the Jacobi method

Successive over-relaxation (SOR)

• Successive over-relaxation (SOR) is an iterative method that improves upon the Gauss-Seidel method by introducing a relaxation parameter $\omega$
• The SOR method was developed by David M. Young and Stanley P. Frankel in the 1950s
• The relaxation parameter $\omega$ is used to control the amount of correction applied to each variable during the iterative process, allowing for faster convergence when chosen appropriately

SOR method algorithm

• Given a linear system $Ax = b$, where $A$ is an $n \times n$ matrix and $x$ and $b$ are $n$-dimensional vectors, the SOR method starts with an initial guess $x^{(0)}$ and a relaxation parameter $\omega$

• At each iteration $k$, the method updates each component of the solution vector $x^{(k+1)}_i$ using the formula: $x^{(k+1)}_i = (1 - \omega) x^{(k)}_i + \frac{\omega}{a_{ii}} \left(b_i - \sum_{j < i} a_{ij} x^{(k+1)}_j - \sum_{j > i} a_{ij} x^{(k)}_j \right)$

• The process is repeated until convergence is reached, i.e., when the difference between consecutive iterations falls below a specified tolerance

SOR vs Gauss-Seidel method

• The SOR method is a generalization of the Gauss-Seidel method, with the Gauss-Seidel method being a special case of SOR when $\omega = 1$
• When the relaxation parameter $\omega$ is chosen optimally, the SOR method can converge significantly faster than the Gauss-Seidel method
• However, finding the optimal value of $\omega$ can be challenging and may require additional computational effort or prior knowledge of the matrix properties

SOR convergence and optimal parameter

• The convergence of the SOR method depends on the properties of the matrix $A$ and the choice of the relaxation parameter $\omega$

• For the SOR method to converge, the matrix $A$ must be symmetric positive definite or strictly diagonally dominant

• The optimal value of $\omega$ that minimizes the spectral radius of the iteration matrix and maximizes the convergence rate is given by: $\omega_{opt} = \frac{2}{1 + \sqrt{1 - \rho(J)^2}}$

  where $\rho(J)$ is the spectral radius of the Jacobi iteration matrix

SOR method example

• Consider the same linear system as in the previous examples: $3x_1 + x_2 = 5$ $x_1 + 2x_2 = 4$

• Starting with an initial guess $x^{(0)} = (0, 0)$ and a relaxation parameter $\omega = 1.2$, the SOR method updates the solution as follows: $x^{(1)}_1 = (1 - 1.2) \cdot 0 + \frac{1.2}{3}(5 - 0) = 2$ $x^{(1)}_2 = (1 - 1.2) \cdot 0 + \frac{1.2}{2}(4 - 2) = 1.2$

  $x^{(2)}_1 = (1 - 1.2) \cdot 2 + \frac{1.2}{3}(5 - 1.2) = 1.04$ $x^{(2)}_2 = (1 - 1.2) \cdot 1.2 + \frac{1.2}{2}(4 - 1.04) = 1.488$

• The process continues until the desired convergence is achieved, with the solution approaching $x = (1, 1.5)$ faster than the Gauss-Seidel method when $\omega$ is chosen appropriately

Conjugate gradient method

• The conjugate gradient (CG) method is an iterative method for solving large, sparse, symmetric positive definite linear systems
• It was developed by Magnus Hestenes and Eduard Stiefel in the 1950s
• The CG method is based on the idea of minimizing a quadratic function associated with the linear system, using conjugate directions to update the solution and residual vectors

Conjugate gradient algorithm

• Given a linear system $Ax = b$, where $A$ is an $n \times n$ symmetric positive definite matrix and $x$ and $b$ are $n$-dimensional vectors, the CG method starts with an initial guess $x^{(0)}$ and initial residual $r^{(0)} = b - Ax^{(0)}$
• At each iteration $k$, the method performs the following steps:
  1. Compute the search direction $p^{(k)}$:
     • If $k = 0$, set $p^{(0)} = r^{(0)}$
     • If $k > 0$, set $\beta^{(k)} = \frac{\langle r^{(k)}, r^{(k)} \rangle}{\langle r^{(k-1)}, r^{(k-1)} \rangle}$ and $p^{(k)} = r^{(k)} + \beta^{(k)} p^{(k-1)}$
  2. Compute the step size $\alpha^{(k)} = \frac{\langle r^{(k)}, r^{(k)} \rangle}{\langle p^{(k)}, Ap^{(k)} \rangle}$
  3. Update the solution $x^{(k+1)} = x^{(k)} + \alpha^{(k)} p^{(k)}$
  4. Update the residual $r^{(k+1)} = r^{(k)} - \alpha^{(k)} Ap^{(k)}$
• The process is repeated until convergence is reached, i.e., when the norm of the residual falls below a specified tolerance

Conjugate gradient vs Jacobi and Gauss-Seidel

• The CG method is generally more efficient than the Jacobi and Gauss-Seidel methods for solving large, sparse, symmetric positive definite linear systems
• Unlike the Jacobi and Gauss-Seidel methods, which update the solution vector component-wise, the CG method updates the entire solution vector at each iteration using conjugate directions
• The CG method has a better convergence rate than the Jacobi and Gauss-Seidel methods, especially for ill-conditioned systems, as it minimizes the error in the $A$-norm rather than the Euclidean norm

Conjugate gradient convergence

• The convergence of the CG method is guaranteed for symmetric positive definite matrices
• The rate of convergence depends on the condition number of the matrix $A$, with faster convergence for well-conditioned matrices
• In exact arithmetic, the CG method converges in at most $n$ iterations, where $n$ is the size of the linear system
• However, in practice, roundoff errors may cause the method to lose conjugacy, requiring additional iterations or restarting strategies

Conjugate gradient example

• Consider the linear system: $2x_1 - x_2 = 1$ $-x_1 + 2x_2 = 0$

• The matrix $A = \begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix}$ is symmetric positive definite

• Starting with an initial guess $x^{(0)} = (0, 0)$, the CG method proceeds as follows:

  $r^{(0)} = b - Ax^{(0)} = (1, 0)$ $p^{(0)} = r^{(0)} = (1, 0)$

  $\alpha^{(0)} = \frac{\langle r^{(0)}, r^{(0)} \rangle}{\langle p^{(0)}, Ap^{(0)} \rangle} = \frac{1}{2}$ $x^{(1)} = x^{(0)} + \alpha^{(0)} p^{(0)} = (0.5, 0)$ $r^{(1)} = r^{(0)} - \alpha^{(0)} Ap^{(0)} = (0, 0.5)$

  $\beta^{(1)} = \frac{\langle r^{(1)}, r^{(1)} \rangle}{\langle r^{(0)}, r^{(0)} \rangle} = \frac{1}{4}$ $p^{(1)} = r