Skip to main content

Subsection 7.3.2 Gauss-Seidel iteration

A variation on the Jacobi iteration is the Gauss-Seidel iteration. It recognizes that since values at points are updated in some order, if a neighboring value has already been updated earlier in the current step, then you might as well use that updated value. For our example from Subsection 7.1.1 this is captured by the algorithm

\begin{equation*} \begin{array}{l} {\bf for~} k = 0, \ldots, \mbox{ convergence} \\ ~~~ {\bf for~} i = 0, \ldots, N \times N-1 \\ ~~~ ~~~ \upsilon_{i}^{(k+1)} = {(h^2 \phi_{i} +\upsilon_{i-N}^{(k+1)} + \upsilon_{i-1}^{(k+1)} + \upsilon_{i+1}^{(k)} + \upsilon_{i+N}^{(k)})}/{4} \\ ~~~ {\bf endfor} \\ {\bf endfor} \end{array} \end{equation*}

modified appropriately for points adjacent to the boundary. This algorithm exploits the fact that \(\upsilon_{i-N}^{(k+1)} \) and \(\upsilon_{i-1}^{(k+1)} \) have already been computed by the time \(\upsilon_{i}^{(k+1)} \) is updated. Once again, the superscripts are there to emphasize the iteration during which a value is updated. In practice, the superscripts can be dropped because of the order in which the computation happens.

Homework 7.3.2.2.

Here we repeat (7.3.1) for Jacobi's iteration applied to the example in Subsection 7.1.1:

\begin{equation} \begin{array}{l} \left( \begin{array} {r r r r | r r r r | r r } 4 \amp \amp \amp \amp \amp \amp \amp \amp \amp \\ \amp 4 \amp \amp \amp \amp \amp \amp \amp \amp \\ \amp \amp 4 \amp \amp \amp \amp \amp \amp \amp \\ \amp \amp \amp 4 \amp \amp \amp \amp \amp \amp \\ \hline \amp \amp \amp \amp 4 \amp \amp \amp \amp \amp \\ \amp \amp \amp \amp \amp 4 \amp \amp \amp \amp \\ \amp \amp \amp \amp \amp \amp 4 \amp \amp \amp \\ \amp \amp \amp \amp \amp \amp \amp 4 \amp \amp \\ \hline \amp \amp \amp \amp \amp \amp \amp \amp 4 \amp \\ \amp \amp \amp \amp \amp \amp \amp \amp \amp \ddots \end{array} \right) \left( \begin{array}{c} \upsilon_0^{(k+1)} \\ \upsilon_1^{(k+1)}\\ \upsilon_2^{(k+1)}\\ \upsilon_3^{(k+1)}\\ \hline \upsilon_4^{(k+1)}\\ \upsilon_5^{(k+1)}\\ \upsilon_6^{(k+1)}\\ \upsilon_7^{(k+1)}\\ \hline \upsilon_8^{(k+1)} \\ \vdots \end{array} \right) \\ ~~~~= \left( \begin{array} {r r r r | r r r r | r r } 0 \amp 1 \amp \amp \amp 1 \amp \amp \amp \amp \amp \\ 1 \amp 0 \amp 1 \amp \amp \amp 1 \amp \amp \amp \amp \\ \amp 1 \amp 0 \amp 1 \amp \amp \amp 1 \amp \amp \amp \\ \amp \amp 1 \amp 0 \amp \amp \amp \amp 1 \amp \amp \\ \hline 1 \amp \amp \amp \amp 0 \amp 1 \amp \amp \amp 1 \amp \\ \amp 1 \amp \amp \amp 1 \amp 0 \amp 1 \amp \amp \amp \ddots \\ \amp \amp 1 \amp \amp \amp 1 \amp 0 \amp 1 \amp \amp \\ \amp \amp \amp 1 \amp \amp \amp 1 \amp 0 \amp \amp \\ \hline \amp \amp \amp \amp 1 \amp \amp \amp \amp 0 \amp \ddots \\ \amp \amp \amp \amp \amp \amp \ddots \amp \amp \ddots \amp \ddots \end{array} \right) \left( \begin{array}{c} \upsilon_0^{(k)} \\ \upsilon_1^{(k)}\\ \upsilon_2^{(k)}\\ \upsilon_3^{(k)}\\ \hline \upsilon_4^{(k)}\\ \upsilon_5^{(k)}\\ \upsilon_6^{(k)}\\ \upsilon_7^{(k)}\\ \hline \upsilon_8^{(k)} \\ \vdots \end{array} \right) + \left( \begin{array}{c} h^2 \phi_0 \\ h^2 \phi_1 \\ h^2 \phi_2 \\ h^2 \phi_3 \\ \hline h^2 \phi_4 \\ h^2 \phi_5 \\ h^2 \phi_6\\ h^2 \phi_7 \\ \hline h^2 \phi_8 \\ \vdots \end{array} \right). \end{array}\label{chapter07-GS-matrix-form}\tag{7.3.2} \end{equation}

Modify this to reflect the Gauss-Seidel iteration.

Solution
\begin{equation*} \begin{array}{l} \left( \begin{array} {r r r r | r r r r | r r } 4 \amp \amp \amp \amp \amp \amp \amp \amp \amp \\ -1 \amp 4 \amp \amp \amp \amp \amp \amp \amp \amp \\ \amp -1 \amp 4 \amp \amp \amp \amp \amp \amp \amp \\ \amp \amp -1 \amp 4 \amp \amp \amp \amp \amp \amp \\ \hline -1 \amp \amp \amp \amp 4 \amp \amp \amp \amp \amp \\ \amp -1 \amp \amp \amp -1 \amp 4 \amp \amp \amp \amp \\ \amp \amp -1 \amp \amp \amp -1 \amp 4 \amp \amp \amp \\ \amp \amp \amp -1 \amp \amp \amp -1 \amp 4 \amp \amp \\ \hline \amp \amp \amp \amp -1 \amp \amp \amp \amp 4 \amp \\ \amp \amp \amp \amp \amp \amp \ddots \amp \amp \ddots \amp \ddots \end{array} \right) \left( \begin{array}{c} \upsilon_0^{(k+1)} \\ \upsilon_1^{(k+1)}\\ \upsilon_2^{(k+1)}\\ \upsilon_3^{(k+1)}\\ \hline \upsilon_4^{(k+1)}\\ \upsilon_5^{(k+1)}\\ \upsilon_6^{(k+1)}\\ \upsilon_7^{(k+1)}\\ \hline \upsilon_8^{(k+1)} \\ \vdots \end{array} \right) \\ ~~~~:= \left( \begin{array} {r r r r | r r r r | r r } 0 \amp 1 \amp \amp \amp 1 \amp \amp \amp \amp \amp \\ \amp 0 \amp 1 \amp \amp \amp 1 \amp \amp \amp \amp \\ \amp \amp 0 \amp 1 \amp \amp \amp 1 \amp \amp \amp \\ \amp \amp \amp 0 \amp \amp \amp \amp 1 \amp \amp \\ \hline \amp \amp \amp \amp 0 \amp 1 \amp \amp \amp 1 \amp \\ \amp \amp \amp \amp \amp 0 \amp 1 \amp \amp \amp \ddots \\ \amp \amp \amp \amp \amp \amp 0 \amp 1 \amp \amp \\ \amp \amp \amp \amp \amp \amp \amp 0 \amp \amp \\ \hline \amp \amp \amp \amp \amp \amp \amp \amp 0 \amp \ddots \\ \amp \amp \amp \amp \amp \amp \amp \amp \amp \ddots \end{array} \right) \left( \begin{array}{c} \upsilon_0^{(k)} \\ \upsilon_1^{(k)}\\ \upsilon_2^{(k)}\\ \upsilon_3^{(k)}\\ \hline \upsilon_4^{(k)}\\ \upsilon_5^{(k)}\\ \upsilon_6^{(k)}\\ \upsilon_7^{(k)}\\ \hline \upsilon_8^{(k)} \\ \vdots \end{array} \right) + \left( \begin{array}{c} h^2 \phi_0 \\ h^2 \phi_1 \\ h^2 \phi_2 \\ h^2 \phi_3 \\ \hline h^2 \phi_4 \\ h^2 \phi_5 \\ h^2 \phi_6\\ h^2 \phi_7 \\ \hline h^2 \phi_8 \\ \vdots \end{array} \right). \end{array} \end{equation*}

This homework suggests the following:

  • We wish to solve \(A x = y \text{.}\)

    We write symmetric \(A \) as

    \begin{equation*} A = \begin{array}[t]{c} \underbrace{ (D-L) } \\ M \end{array} - \begin{array}[t]{c} \underbrace{ (L^T) } \\ N \end{array}, \end{equation*}

    where \(-L \) equals the strictly lower triangular part of \(A \) and \(D \) is its diagonal.

  • We then notice that

    \begin{equation*} A x = y \end{equation*}

    can be rewritten as

    \begin{equation*} ( D - L - L^T )x = y \end{equation*}

    or, equivalently,

    \begin{equation*} ( D - L ) x = L^T x + y. \end{equation*}

    If you think about it carefully, this captures (7.3.2) for our example. Finally,

    \begin{equation*} x = (D - L)^{-1} ( L^T x + y ). \end{equation*}

  • If we now let \(x^{(k)} \) be the values of our vector \(x \) in the current step. Then the values after all elements have been updated are given by the vector

    \begin{equation*} x^{(k+1)} = (D - L)^{-1} ( L^T x^{(k)} + y ). \end{equation*}
Homework 7.3.2.3.

When the Gauss-Seidel iteration is used to solve \(A x = y \text{,}\) where \(A \in \Rnxn \text{,}\) it computes entries of \(x^{(k+1)} \) in the forward order \(\chi_0^{(k+1)}, \chi_1^{(k+1)}, \ldots \text{.}\) If \(A = D - L - L^T \text{,}\) this is captured by

\begin{equation} ( D - L ) x^{(k+1)} = L^T x^{(k)} + y . \label{chapter07-GS-forward}\tag{7.3.3} \end{equation}

Modify (7.3.3) to yield a "reverse" Gauss-Seidel method that computes the entries of vector \(x^{(k+1)} \) in the order \(\chi_{n-1}^{(k+1)}, \chi_{n-2}^{(k+1)}, \ldots \text{.}\)

Solution

The reverse order is given by \(\chi_{n-1}^{(k+1)}, \chi_{n-2}^{(k+1)}, \ldots \text{.}\) This corresponds to the splitting \(M = D - L^T \) and \(N = L \) so that

\begin{equation*} ( D - L^T ) x^{(k+1)} = L x^{(k)} + y . \end{equation*}
Homework 7.3.2.4.

A "symmetric" Gauss-Seidel iteration to solve symmetric \(A x = y \text{,}\) where \(A \in \Rnxn \text{,}\) alternates between computing entries in forward and reverse order. In other words, if \(A = M_F - N_F \) for the forward Gauss-Seidel method and \(A = M_R - N_R \) for the reverse Gauss-Seidel method, then

\begin{equation*} \begin{array}{rcl} M_F x^{(k+\frac{1}{2})} \amp = \amp N_F x^{(k)} + y \\ M_R x^{(k+1)} \amp = \amp N_R x^{(k+\frac{1}{2})} + y \end{array} \end{equation*}

constitutes one iteration of this symmetric Gauss-Seidel iteration. Determine \(M \) and \(N \) such that

\begin{equation*} M x^{(k+1)} = N x^{(k)} + y \end{equation*}

equals one iteration of the symmetric Gauss-Seidel iteration.

(You may want to follow the hint...)

Hint

  • From this unit and the last homework, we know that \(M_F = ( D - L ) \text{,}\) \(N_F = L^T \text{,}\) \(M_R = ( D - L^T ) \text{,}\) and \(N_R = L \text{.}\)

  • Show that

    \begin{equation*} ( D - L^T ) x^{(k+1)} = L ( D - L )^{-1} L^T x^{(k)} + ( I + L ( D - L )^{-1} ) y. \end{equation*}

  • Show that \(I + L ( D - L )^{-1} = D ( D - L )^{-1} \text{.}\)

  • Use these insights to determine \(M \) and \(N \text{.}\)

Solution

  • From this unit and the last homework, we know that \(M_F = ( D - L ) \text{,}\) \(N_F = L^T \text{,}\) \(M_R = ( D - L^T ) \text{,}\) and \(N_R = L \text{.}\)

  • Show that

    \begin{equation*} ( D - L^T ) x^{(k+1)} = L ( D - L )^{-1} L^T x^{(k)} + ( I + L ( D - L )^{-1} ) y. \end{equation*}

    We show this by substituting \(M_R \) and \(N_R \text{:}\)

    \begin{equation*} ( D - L^T ) x^{(k+1)} = L x^{(k+ \frac{1}{2})}+ y \end{equation*}

    and then substituting in for \(x^{(k+ \frac{1}{2})} \text{,}\) \(M_F \) and \(N_F \text{:}\)

    \begin{equation*} ( D - L^T ) x^{(k+1)} = L ( ( D - L )^{-1} L^T x^{(k)} + y ) + y. \end{equation*}

    Multiplying out the right-hand side and factoring out \(y \) yields the desired result.

  • Show that \(I + L ( D - L )^{-1} = D ( D - L )^{-1} \text{.}\)

    We show this by noting that

    \begin{equation*} \begin{array}{l} I + L ( D - L )^{-1} \\ ~~~=~~~~ \\ ( D - L )( D - L )^{-1} + L ( D - L )^{-1} ~~~=~~~~ \\ ( D - L + L ) ( D - L )^{-1} ~~~=~~~~ \\ D ( D - L )^{-1}. \end{array} \end{equation*}

  • Use these insights to determine \(M \) and \(N \text{.}\)

    We now notice that

    \begin{equation*} ( D - L^T ) x^{(k+1)} = L ( D - L )^{-1} L^T x^{(k)} + ( I + L ( D - L )^{-1} ) y \end{equation*}

    can be rewritten as (Someone check this... My brain hurts.)

    \begin{equation*} ( D - L^T ) x^{(k+1)} = L ( D - L )^{-1} L^T x^{(k)} + D ( D - L )^{-1} y \end{equation*}
    \begin{equation*} \begin{array}[t]{c} \underbrace{ ( D - L )D^{-1} ( D - L^T ) } \\ M \end{array} x^{(k+1)} = \begin{array}[t]{c} \underbrace{ ( D - L )D^{-1} L ( D - L )^{-1} L^T }\\ N \end{array} x^{(k)} + y \end{equation*}
login