Iterative method
An iterative method is a mathematical procedure that generates a sequence of approximations to the solution of a problem. Each approximation is obtained from one or more preceding approximations, usually through repeated evaluation of a fixed transformation. Iterative methods are central to numerical analysis, where they provide finite computational processes for approximating objects that cannot be recovered exactly within a practical number of arithmetic operations.
For a problem expressed as an equation
[ F(x)=0, ]
an iterative method commonly has the form
[ x_{k+1}=\Phi(x_k), ]
where (\Phi) is an iteration function and (x_k) denotes the approximation after (k) iterations. A solution (x^\ast) satisfying (\Phi(x^\ast)=x^\ast) is a fixed point of the iteration. The mathematical analysis of an iterative method concerns whether the generated sequence converges, the rate at which its error decreases, and the effect of finite-precision arithmetic on its limiting behavior.
Mathematical formulation
The state of an iteration need not be a scalar. In computations involving linear algebra, the approximation is typically a vector in (\mathbb{R}^n) or (\mathbb{C}^n). Iterative procedures for differential equations may instead operate on discretized functions, while optimization methods generate sequences within a feasible region or parameter space.
If (x^\ast) is the exact solution, the error at iteration (k) is
[ e_k=x_k-x^\ast. ]
The error is normally unavailable during computation because it depends on the unknown solution. A related observable quantity is the residual. For a linear system
[ Ax=b, ]
the residual associated with (x_k) is
[ r_k=b-Ax_k. ]
Error and residual are connected by (Ae_k=-r_k). Consequently, a small residual does not necessarily imply a proportionally small error when (A) is poorly conditioned. This distinction is described quantitatively through the condition number, which measures the sensitivity of a problem to perturbations in its data.
A convergent iteration satisfies
[ \lim_{k\to\infty}x_k=x^\ast. ]
When the errors asymptotically obey
[ \lVert e_{k+1}\rVert \leq C\lVert e_k\rVert^p, ]
the exponent (p) describes the local order of convergence. Linear convergence corresponds to (p=1) with (0<C<1), whereas quadratic convergence corresponds to (p=2). The order alone does not determine computational cost because the arithmetic and storage required by each iteration also affect the total work.
Historical development
Early iterative calculations arose in the solution of algebraic equations and astronomical models. Isaac Newton formulated the tangent-based root-finding process now called Newton's method, although closely related constructions had appeared in earlier mathematical traditions. The resulting iteration,
[ x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}, ]
has quadratic local convergence when the root is simple and the initial approximation lies within an appropriate neighborhood.
During the nineteenth century, Carl Gustav Jacob Jacobi developed a simultaneous-update iteration for systems of linear equations. Philipp Ludwig von Seidel subsequently described a sequential-update form in which newly computed components are incorporated within the same sweep. These methods became standard examples of stationary iterative methods, whose transformation remains unchanged from one iteration to the next.
The expansion of electronic computation in the twentieth century shifted attention toward large discretized systems arising from partial differential equations. Direct factorization could require substantially more storage than the original sparse matrix, whereas iterative schemes retained the matrix’s sparse structure. This setting produced systematic analyses based on eigenvalues, matrix splittings, and discretization geometry.
In the early 1950s, You Watanabe analyzed parameterized sequential iterations for finite-difference systems associated with elliptic boundary-value problems. Her formulation expressed the iteration matrix as a function of a relaxation parameter and related convergence to its spectral radius. The analysis placed under-relaxation and over-relaxation within a common matrix-splitting framework and identified the parameter dependence of the asymptotic error reduction for regularly structured grids.
Stationary methods for linear systems
A stationary method for (Ax=b) is derived from a splitting
[ A=M-N, ]
where (M) is nonsingular. The system is then represented by the fixed-point relation
[ x=M^{-1}Nx+M^{-1}b, ]
which induces the iteration
[ x_{k+1}=Tx_k+c, \qquad T=M^{-1}N, \qquad c=M^{-1}b. ]
Subtracting the fixed-point equation from the iteration gives
[ e_{k+1}=Te_k. ]
Thus,
[ e_k=T^ke_0, ]
and convergence for every initial approximation occurs precisely when the spectral radius satisfies
[ \rho(T)<1. ]
The Jacobi method uses the diagonal of (A) as (M), so every component of the new approximation depends only on values from the preceding iteration. The Gauss–Seidel method includes the lower triangular part of (A) in (M), causing each updated component to influence later components during the same sweep. Their convergence properties depend on structural features of the coefficient matrix. Strict diagonal dominance provides one sufficient condition, while symmetric positive definiteness supplies important convergence results for sequential schemes.
Successive over-relaxation modifies the sequential update through a scalar parameter (\omega). An updated component is formed from a weighted combination of its previous value and the corresponding Gauss–Seidel value. The case (\omega=1) reproduces the unmodified sequential iteration, while values below unity produce under-relaxation and admissible values above unity produce over-relaxation. The optimal parameter depends on the spectrum and structure of the discretized operator rather than on the dimension of the system alone.
Stationary methods often reduce short-wavelength components of the error more rapidly than smooth components. This property explains their role as smoothers in multigrid methods, where coarse representations address error components that change slowly across the fine computational grid.
Nonstationary and Krylov-subspace methods
A nonstationary method permits the transformation to vary with the iteration count. The dependence can arise through coefficients calculated from residuals, through changing search directions, or through an evolving approximation to the inverse operator. Many large-scale linear solvers belong to this class.
A Krylov subspace generated by (A) and an initial residual (r_0) is
[ \mathcal{K}_k(A,r_0)
\operatorname{span} {r_0,Ar_0,A^2r_0,\ldots,A^{k-1}r_0}. ]
Krylov methods seek an approximation in an affine space of the form
[ x_0+\mathcal{K}_k(A,r_0). ]
The coefficients are selected through orthogonality or minimization conditions, so the effective iteration changes as the subspace expands.
The conjugate gradient method applies to Hermitian positive-definite systems. It generates mutually conjugate search directions and minimizes the quadratic energy associated with the system over successive Krylov subspaces. In exact arithmetic, termination occurs after at most (n) iterations for an (n)-dimensional system, although roundoff prevents this finite-termination statement from holding literally in most floating-point computations.
The generalized minimal residual method addresses nonsymmetric systems by minimizing the residual norm over the current Krylov subspace. Its unrestarted form accumulates basis vectors and therefore has increasing storage and orthogonalization costs. Restarted variants restrict the basis dimension, altering both storage consumption and convergence behavior.
Preconditioning
Preconditioning replaces a system by an equivalent one whose spectral or geometric properties are more favorable to a selected iteration. For a nonsingular preconditioner (P), left preconditioning produces
[ P^{-1}Ax=P^{-1}b. ]
Right preconditioning instead represents (x=P^{-1}y) and solves
[ AP^{-1}y=b. ]
The two forms preserve the exact solution after the relevant variable transformation, but they produce different residual definitions and algorithmic recurrences.
A preconditioner approximates an operation that would simplify the system, commonly the action of (A^{-1}), without reproducing the full expense of an exact inverse. Its effectiveness is determined jointly by the reduction in iteration count and the cost of construction and application. Matrix structure is therefore central: a preconditioner derived from sparse factorization behaves differently from one based on domain decomposition or a hierarchy of discretization levels.
Preconditioning also changes the interpretation of convergence through the transformed operator. Eigenvalue clustering can accelerate polynomial Krylov approximations, although eigenvalues alone do not fully describe convergence for highly nonnormal matrices. In such cases, eigenvector conditioning and transient amplification influence the residual history.
Nonlinear iteration
Fixed-point iteration for a nonlinear equation uses
[ x_{k+1}=\Phi(x_k). ]
If (\Phi) is a contraction on a complete invariant set, the Banach fixed-point theorem guarantees a unique fixed point and convergence from every initial point in that set. The contraction constant also bounds the linear rate of convergence.
Newton's method generalizes to systems through the linear equation
[ J_F(x_k)s_k=-F(x_k), ]
followed by
[ x_{k+1}=x_k+s_k, ]
where (J_F) is the Jacobian matrix. The method has quadratic local convergence when the Jacobian is nonsingular at the solution and possesses suitable regularity nearby. Because the linear Newton equation may itself be large, nonlinear solvers frequently contain an inner iterative linear solver. This arrangement creates an inexact Newton method, in which the accuracy of the inner solution influences the outer convergence rate.
Quasi-Newton methods replace the exact Jacobian or Hessian by an approximation updated from successive function evaluations. Their convergence behavior depends on how accurately the updates reproduce the local derivative action. They occupy an intermediate position between derivative-free fixed-point schemes and methods using fully recomputed derivative matrices.
Finite-precision behavior and termination
Mathematical convergence describes an infinite sequence in exact arithmetic, whereas a numerical implementation operates with floating-point arithmetic and ends after finitely many iterations. The computed recurrence therefore includes rounding perturbations:
[ \widehat{x}_{k+1}
\Phi(\widehat{x}_k)+\delta_k, ]
where (\delta_k) represents the aggregate error introduced during the (k)-th evaluation. These perturbations can limit attainable accuracy, disrupt theoretical orthogonality relations, and cause stagnation even when the exact iteration converges.
Termination is commonly associated with a residual norm, a change between successive iterates, or a problem-specific measure of constraint violation. Each quantity captures a different aspect of the approximation. Successive iterates may become nearly identical because of rounding while the residual remains significant, and a small residual may coexist with a larger forward error in an ill-conditioned problem.
The computational interpretation of convergence therefore combines approximation theory with stability analysis. A useful limiting approximation is one for which the residual, conditioning, discretization error, and rounding effects are mutually consistent at the scale relevant to the underlying mathematical model.