Why do we need approximation methods when we have algorithms to find exact roots? While I was studying numerical methods and optimizations recently, I observed that whenever we find a root to an equation or a system of linear equations, we always find approximate roots. However, we already have algorithms for finding exact roots.

For example, the Gauss-Seidel or Gauss-Jacobi methods of iterative solution tend to find the "most correct" approximate root to the system of linear equations, although we can find exact roots to a system of linear equations through algorithms such as Gauss elimination or Gauss-Jordan method.


  What is the reason to use approximation methods when we can find exact roots?
There are at least two important reasons.

First, iterative methods can be significantly faster, and as there are a lot of problems with big systems (hundreds of thousands or more variables), speed can be important. Also, precise methods usually can't take advantage of sparsity, which is common case too.

Second, we usually use finite-precision floating point numbers, so we have to deal with rounding errors, which tend to accumulate in precise methods more. For example, in Gauss elimination, if you get a row slightly wrong, it will affect all subsequent rows. As in iterative methods the "main part" doesn't change, they are (usually) less susceptible to it.
For polynomials: almost no polynomial's roots can be expressed exactly.  All linear, quadratic, cubic and quartic polynomials can be factored exactly.  The generic case for quintics is similar to that for $x^5 + x + 5$, whose five roots are not expressible using any finite expression of sums, differences, products, quotients, or rational powers of rational numbers.  (They can be expressed using Jacobi elliptic functions, which can only be evaluated exactly for very special arguments.)  For more on this see solvability by radicals in Galois theory.  Generically, for degree above four, nearly all polynomials are not solvable by radicals.  As a consequence, the only way forward is to find approximate roots.

Now to linear systems of equations.

There are two resources that any algorithm requires: time and space.  For matrix algorithms, it is typical to express the complexities as functions of the size of the matrix.  For $n \times n$ matrices, the size is $n$.  Also, real implementations of algorithms can increase or decrease the resource usage by nearly constant multiplicative factors.  (Implementing on this hardware may be twice as fast as that hardware.  Using this representation of the numbers may use ten times as much space as using that representation.)  Consequently, constant multipliers are typically ignored when calculating and reporting time complexity -- the "$2$" in "$2 n^2$" is meaningless when we can get a factor of $2$ just by moving to a different computing device.  We indicate that we are ignoring constant multipliers by using big O notation.  Also, efficiency is largely irrelevant for small problem sizes, we typically think of $n$ as being large.  Consequently, although a detailed accounting of resource usage may show, for instance, usage of $\frac{2}{3}n^3 + \frac{3}{2}n^2 - \frac{7}{6}n$ units of a resource, this is reported as $O(n^3)$ because, for large $n$, the lower order terms are irrelevant and as already noted, the leading coefficient is irrelevant.

The following are widely reported, so I don't provide specific citations.  We report space complexity assuming floating point representations.  (Use of exact representations in the worst case lead to exponential space consumption.)

Gaussian elimination and Gauss-Jordan elimination have time complexity $O(n^3)$ and space complexity $O(n^2)$.

(Gauss-)Jacobi iteration has time complexity $O(n^2)$ and space complexity $O(n^2)$.

Gauss-Seidel has time complexity $O(n^2)$ and space complexity $O(n^2)$.

For these algorithms, space complexity is not a concern.  However, time complexity does vary.  So, suppose we have a problem instance with $n = 10^6$.  (This is not actually so large, but is illustrative.)  Then we can expect the ratio of time required for Gauss or Gauss-Jordan to the time required for Gauss-Jacobi or Gauss-Seidel to be
$$  \frac{(10^6)^3}{(10^6)^2} = 10^6  \text{.}  $$
That is, the approximate algorithms are about a million times faster for this problem.

There are additional considerations.  For instance, we are discussing calculations using floating point numbers.  Such representations have finite precision, so long calculations risk accumulation of error, sometimes to the point that the "solution" returned by the algorithm is garbage.  We capture this by describing for which inputs the algorithm is "stable".  If the algorithm is stable for a particular input, accumulated errors do not dominate the output.  Gaussian and Gauss-Jordan are stable for diagonally-dominant matrices.  If your matrix is not diagonally dominant, then the large multiplications to cancel off-diagonal terms amplify the precision errors in those terms (and the subsequent subtractions can lead to catastrophic cancellation).  Gauss-Jacobi requires the spectral radius of the matrix to be less than $1$ to even converge.  One way to ensure the spectral radius is less than $1$ is to be diagonally dominant.  There are symmetric positive-definite matrices for which Gauss-Jacobi does not even converge, e.g. $(\begin{smallmatrix} 1 &amp; 2 \\ 2 &amp; 1 \end{smallmatrix})$.  Gauss-Seidel converges for positive-definite and for diagonally dominant matrices.  Notice that none of the methods you listed are stable/converge if your matrix is far from diagonally dominant.  (So what do you do if you have such a matrix?  Precondition.)

To summarize: We have these many algorithms because (1) the running time can vary substantially, depending on which algorithm you use and (2) different algorithms have different stability requirements.