Three Ways to Solve the Same System
For a system \(\mathbf{A}\mathbf{x}=\mathbf{b}\), LU factorization, Jacobi iteration, and Gauss-Seidel iteration organize the work in different ways. LU is a direct method: it factors the coefficient matrix and then solves triangular systems. Jacobi and Gauss-Seidel are iterative methods: they generate a sequence of approximate solutions, which may converge to the exact solution under suitable conditions. No single method is best for every problem. The choice depends especially on the matrix, how many right-hand sides must be solved, whether an iteration converges quickly enough, and the accuracy required.
Here, operation counts are rough counts of arithmetic operations for dense matrices; they describe how work grows with the system size, not the precise time on every computer. A right-hand side means a particular vector \(\mathbf{b}\) in \(\mathbf{A}\mathbf{x}=\mathbf{b}\). When several systems have the same \(\mathbf{A}\) but different right-hand sides, the cost of factoring \(\mathbf{A}\) can be shared.
What Each Method Costs
For a dense \(n\times n\) matrix, LU factorization takes about \(\tfrac23 n^3\) arithmetic operations in the usual leading-order count. Once the factors are available, forward and backward substitution for one right-hand side take about \(2n^2\) operations together. Thus, a new right-hand side requires triangular solves, not a new factorization. If the matrix changes, in general the factors must be recomputed.
A Jacobi or Gauss-Seidel sweep updates every component using the entries in a row of the system. For a dense matrix, one sweep takes about \(2n^2\) arithmetic operations. The number of sweeps needed is not fixed: it depends on the matrix, the starting vector, and the desired accuracy. Consequently, a low cost per sweep is not by itself evidence that an iterative method is cheaper overall.
| Method | Main work | Useful when |
|---|---|---|
| LU | Factor once, then solve triangular systems | A reliable direct solve is wanted, especially for several right-hand sides with the same matrix |
| Jacobi | Repeated simultaneous updates | The iteration converges and its simple updates suit the problem or implementation |
| Gauss-Seidel | Repeated updates using new values within each sweep | The iteration converges and using newly updated components is advantageous |
Proof. The triangular systems have unique solutions: forward substitution determines the entries of \(\mathbf{y}\) in order because each diagonal entry of \(\mathbf{L}\) is nonzero, and backward substitution determines the entries of \(\mathbf{x}\) in reverse order because each diagonal entry of \(\mathbf{U}\) is nonzero. From \(\mathbf{L}\mathbf{y}=\mathbf{P}\mathbf{b}\) and \(\mathbf{U}\mathbf{x}=\mathbf{y}\), we get \(\mathbf{L}\mathbf{U}\mathbf{x}=\mathbf{P}\mathbf{b}\). Since \(\mathbf{L}\mathbf{U}=\mathbf{P}\mathbf{A}\), this is \(\mathbf{P}\mathbf{A}\mathbf{x}=\mathbf{P}\mathbf{b}\). A permutation matrix is invertible, so \(\mathbf{A}\mathbf{x}=\mathbf{b}\). Nonsingularity of \(\mathbf{A}\) makes this solution unique. The equations for new right-hand sides use the same \(\mathbf{L}\) and \(\mathbf{U}\); only the triangular solves change. \(\square\)
Worked Example: One Factorization, Two Right-Hand Sides
Factor \(\mathbf{A}=\begin{bmatrix}4&1\\2&3\end{bmatrix}\) without row exchanges, and use the factors to solve for \(\mathbf{b}=\begin{bmatrix}11\\13\end{bmatrix}\) and \(\mathbf{c}=\begin{bmatrix}1\\2\end{bmatrix}\).
Eliminating the entry below the first pivot uses the multiplier \(\tfrac12\). Therefore
For \(\mathbf{b}\), forward substitution gives \(y_1=11\) and \(\tfrac12(11)+y_2=13\), so \(y_2=\tfrac{15}{2}\). Back substitution gives \(\tfrac52 x_2=\tfrac{15}{2}\), hence \(x_2=3\), and \(4x_1+3=11\), hence \(x_1=2\). Check: \(\mathbf{A}\begin{bmatrix}2\\3\end{bmatrix}=\begin{bmatrix}8+3\\4+9\end{bmatrix}=\begin{bmatrix}11\\13\end{bmatrix}\).
For \(\mathbf{c}\), forward substitution gives \(y_1=1\) and \(\tfrac12+y_2=2\), so \(y_2=\tfrac32\). Back substitution gives \(x_2=\tfrac35\), then \(4x_1+\tfrac35=1\), so \(x_1=\tfrac1{10}\). Check: \(\mathbf{A}\begin{bmatrix}\tfrac1{10}\\\tfrac35\end{bmatrix}=\begin{bmatrix}\tfrac25+\tfrac35\\\tfrac15+\tfrac95\end{bmatrix}=\begin{bmatrix}1\\2\end{bmatrix}\). The factorization was used for both systems; only the solves were repeated.
Iterations, Progress, and Residuals
Jacobi and Gauss-Seidel split a system into component equations and update the components repeatedly. In Jacobi, every new component uses values from the previous iterate. In Gauss-Seidel, each new component uses the newest available values within the same sweep. The convergence tests established in the earlier material on Jacobi and Gauss-Seidel give conditions under which these sequences converge. Even when convergence is known, the number of sweeps needed to reach a useful residual is part of the practical cost.
Proof. Subtract the update from the fixed-point equation: \(\mathbf{x}^*-\mathbf{x}^{(k+1)}=(\mathbf{B}\mathbf{x}^*+\mathbf{c})-(\mathbf{B}\mathbf{x}^{(k)}+\mathbf{c})=\mathbf{B}(\mathbf{x}^*-\mathbf{x}^{(k)})\). Thus \(\mathbf{e}^{(k+1)}=\mathbf{B}\mathbf{e}^{(k)}\). For the residual, the exact solution satisfies \(\mathbf{b}=\mathbf{A}\mathbf{x}^*\), so \(\mathbf{r}^{(k)}=\mathbf{b}-\mathbf{A}\mathbf{x}^{(k)}=\mathbf{A}(\mathbf{x}^*-\mathbf{x}^{(k)})=\mathbf{A}\mathbf{e}^{(k)}\). \(\square\)
Worked Example: Comparing One Jacobi and One Gauss-Seidel Sweep
For \(\mathbf{A}=\begin{bmatrix}4&1\\2&3\end{bmatrix}\), solve \(\mathbf{A}\mathbf{x}=\begin{bmatrix}11\\13\end{bmatrix}\) by one sweep of each method from \(\mathbf{x}^{(0)}=\begin{bmatrix}0\\0\end{bmatrix}\). The component equations are \(x_1=(11-x_2)/4\) and \(x_2=(13-2x_1)/3\).
Jacobi uses only old values. Thus \(x_1^{(1)}=11/4\), followed by \(x_2^{(1)}=13/3\), and \(\mathbf{x}_{J}^{(1)}=\begin{bmatrix}11/4\\13/3\end{bmatrix}\). Substitution into the update equations checks these values: \((11-0)/4=11/4\) and \((13-2(0))/3=13/3\). Its residual is
Gauss-Seidel first obtains \(x_1^{(1)}=11/4\), then uses that new value to obtain \(x_2^{(1)}=(13-2(11/4))/3=5/2\). Its iterate is \(\mathbf{x}_{GS}^{(1)}=\begin{bmatrix}11/4\\5/2\end{bmatrix}\), and the update checks are \((11-0)/4=11/4\) and \((13-2(11/4))/3=5/2\). Its residual is
In this sweep, Gauss-Seidel has a smaller maximum absolute residual. That observation is specific to this example and this starting vector; it is not a general guarantee that every Gauss-Seidel problem will converge faster. Both methods are convergent here by the strict diagonal-dominance convergence test, since \(4>1\) and \(3>2\). The exact solution, already checked in the LU example, is \(\begin{bmatrix}2\\3\end{bmatrix}\).
Use Residuals Carefully
Residuals are useful for checking the equations actually solved, and are often monitored during iteration. The same residual definition applies to a direct-method result or an iterative approximation. To assess solution error, however, the residual must be interpreted together with information about the matrix, such as its conditioning; a small residual by itself is not a certificate of small forward error.
Worked Example: Checking Each Component of a Residual
Let \(\mathbf{A}=\begin{bmatrix}1&9/10\\4/5&1\end{bmatrix}\), \(\mathbf{b}=\begin{bmatrix}19/10\\17271/10000\end{bmatrix}\), and \(\widetilde{\mathbf{x}}=\begin{bmatrix}1729/1000\\3439/10000\end{bmatrix}\). Calculate \(\mathbf{r}=\mathbf{b}-\mathbf{A}\widetilde{\mathbf{x}}\) exactly.
The first component of \(\mathbf{A}\widetilde{\mathbf{x}}\) is \(\tfrac{1729}{1000}+\tfrac9{10}\tfrac{3439}{10000} =\tfrac{172900}{100000}+\tfrac{30951}{100000} =\tfrac{203851}{100000}\). Thus the first residual component is
The second component of \(\mathbf{A}\widetilde{\mathbf{x}}\) is \(\tfrac45\tfrac{1729}{1000}+\tfrac{3439}{10000} =\tfrac{13832}{10000}+\tfrac{3439}{10000} =\tfrac{17271}{10000}\), so the second residual component is \(0\). Therefore \(\mathbf{r}=\begin{bmatrix}-13851/100000\\0\end{bmatrix}\). The first equation is not satisfied exactly, even though the second one is. In particular, replacing \(203851\) by \(187231\) in the first-component calculation would give an incorrect residual.
A Practical Choice
A useful decision begins with the structure and reuse pattern of the problem, not just the name of a method.
- Choose LU when a dependable direct solve is wanted and the matrix is small or moderate in size. It is especially attractive when many right-hand sides share the same matrix, since the factorization is performed once.
- Consider Jacobi or Gauss-Seidel when the convergence tests support the iteration and the cost of repeated sweeps is reasonable. Check the residual during the process rather than assuming a fixed number of sweeps is enough.
- Account for the matrix structure. The dense operation counts above can overstate the work for matrices with many zero entries; in such cases the actual arithmetic depends on how that structure is stored and used.
- Check the result. Substitute a direct solution back into the system, or calculate the residual of an iterate. If the requested accuracy is not reached, continue the iteration only when convergence behavior supports doing so.
Common Pitfall: Counting Sweeps Instead of Total Work
A common mistake is to call an iterative method cheaper because one sweep is simpler than an LU factorization. For a dense system, if \(k\) sweeps are needed, the iterative work is roughly \(2kn^2\), while LU followed by one solve is roughly \(\tfrac23n^3+2n^2\). These are approximate leading-order counts, but they show why the number of sweeps matters. If the same matrix is used for \(s\) right-hand sides, LU costs roughly \(\tfrac23n^3+2sn^2\), whereas separately iterating on all right-hand sides costs roughly \(2skn^2\). The comparison can change substantially as \(s\) or \(k\) changes. These estimates do not replace convergence or accuracy checks; they help identify which method is worth considering.
Check Your Understanding
Use the definitions and comparisons above to decide what information is relevant before selecting a solver.
- Why can one LU factorization serve several systems with the same coefficient matrix but different right-hand sides?
- For the system in the worked iteration example, compute the second Jacobi iterate from its first iterate.
- What does a zero component of a residual say about the corresponding equation?
- Why does one inexpensive sweep not establish that Jacobi or Gauss-Seidel will be cheaper than LU overall?
- For what change to a sequence of systems must the LU factors generally be recomputed?
