When we study an iterative method from textbooks, for example, see the Gauss-Seidel Method, the given matrix is decomposed with suitable splittings. In the example, $A = L+U$. So we can proceed with that decomposition if we solve the exercise through pen and paper.

I don't understand, why when we implement that method on a computer, that decomposition is not considered anymore?

For example, on this page- implementation with Matlab, it is not used the decomposition introduced before. Why? Why we have to struggle to find a suitable decomposition that will not be used within the code?


1 Answer 1


In Gauss-Seidel, you are using this $A=L+U$ splitting implicitly. So, you never form $L$, $U$, or $L^{-1}$ explicitly. Which is extremely good, since forming an additional matrix (not even talking about a calculation of an explicit inverse) is a huge burden.

Instead, since $L$ is lower triangular, you can change the explicit inverse of $L$, by performing a row-by-row forward subsitution.

Notice, that solving for (via the forward substitution when $L$ is lower-triangular):

$$ \begin{equation} L\mathbf{y}=\mathbf{b} \label{1} \tag{1} \end{equation} $$

is theoretically the same as performing a matrix-vector product: $$ \begin{equation} \mathbf{y}=L^{-1}\mathbf{b} \label{2} \tag{2} \end{equation} $$ However, numerically you always want $\eqref{1}$, not $\eqref{2}$. The reasons are simple: computation of the matrix inverse is numerically unstable and has a huge cost. (see Q1 (especially this answer), Q2 for some additional details).

With that in mind, take a look at the expression from Wikipedia article on Gauss-Seidel. The iteration: $$ \mathbf{x}^{(k+1)}=L^{-1}(\mathbf{b}-U\mathbf{x}^{(k)}) $$ is totally equivalent to a for-loop for $i$ $$ \begin{equation} x^{(k+1)}_i=\frac{1}{a_{ii}}\left(b_i-\sum\limits_{j=1}^{i-1}a_{ij}x_{j}^{(k+1)}\right)-\frac{1}{a_{ii}}\left(\sum\limits_{j=i+1}^{n}a_{ij}x_j^{(k)}\right) \label{3} \tag{3} \end{equation} $$

Here, $a_{ij}$ are the entries of the original matrix $A=L+U$, $k$ is the iteration number, $n$ is the dimension of the square matrix $A$. By looking at $\eqref{3}$, you still see those $L$ and $U$; however, they are present only implicitly, and all operations are done with the entries of the original matrix $A$ stored exactly as it was before.


Matrix decomposition usually implies factorization, so $A=L+U$ is not really (strictly) one according to the definition.

$A=L+U$ splitting (used in this question) has nothing to do with a more famous $A=LU$ LU decomposition.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.