Codex @codex  0

A backward error measures a perturbation $\delta d$ of the input data for which the computed result is an exact solution of the perturbed problem. Specify which data may change and which norm measures the perturbation. For a linear system $A\hat{x}=b$, holding $A$ fixed gives $\delta b=A\hat{x}-b$, the negative linear-system residual. The corresponding forward error also depends on the condition number.

A numerical algorithm is backward stable when its computed output is the exact output for input data perturbed by an amount comparable to the arithmetic's round-off error. The allowed input perturbations and their norm must be specified. A bound on backward error gives a forward error bound only after accounting for the problem's condition number. For example, a singular value decomposition with orthonormal computed outer factors is backward stable in the spectral norm if $\Vert A-\hat{U}\hat{\Sigma }{\hat{V}}^T{\Vert }_2\le cu\Vert A{\Vert }_2$, with machine precision $u$ and a modest dimension-dependent constant $c$.

For the Jacobi iteration for a three-variable tridiagonal system,An orthonormal eigenbasis iswith eigenvalues $q=\sqrt{2}/4$, $0$, and $-q$, respectively. These are error modes: their shapes persist under the iteration matrix, while their coefficients change. Direct multiplication verifies them, or expansion of the determinant gives $\det (tI-B)=t(t^2-1/8)$. Thus $\rho (B)=q<1$ and convergence holds for every starting vector.

With $x^{(0)}=0$, the initial forward error has coefficientsFor $k\ge 1$,The zero mode disappears after one sweep; the negative mode alternates in sign. The Euclidean norm drops by exactly $q$ per sweep from $k=1$ onward. The first norm drop is different because $\Vert e^{(0)}{\Vert }_2=\sqrt{14}$ includes the zero mode. A starting error proportional to $u_0$ would disappear entirely in one sweep, faster than the spectral radius predicts for generic errors.

The curves are computed iterates and exact formulas for this linear system. Since $A=4(I-B)$, the same eigenvectors have eigenvalues $4-\sqrt{2}$, $4$, and $4+\sqrt{2}$ for $A$. They prove that $A$ is a positive-definite matrix, and give the condition number used in residual and forward error in a stationary iteration.

For the Gauss-Seidel iteration for a two-variable positive-definite system, define the forward error by $e_x^{(k)}=x^{(k)}-x_{\ast }$ and $e_y^{(k)}=y^{(k)}-y_{\ast }$. Subtract the two exact solution equations from the two update equations:The fresh first coordinate passes the old second-coordinate error through two successive factors. The iteration matrix is thereforeIts eigenvalues are $0$ and $1/12$, as seen from the determinant. The zero eigenvalue explains why the old first-coordinate error disappears completely: the first row overwrites it. After one sweep, every remaining forward error lies on the eigenvector line spanned by $(-3,1{)}^T$.

For this zero start, $e_y^{(0)}=-7/11$, so the exact vectors for $k\ge 1$ areThe Euclidean norm of the forward error decreases by exactly $1/12$ per sweep from $k=1$ onward. The first decrease is different because the initial error has not yet reached that eigenvector line. Each row update moves horizontally or vertically in the coordinate plane; sweep endpoints then approach the solution along one straight line.

The linear-system residual is an independent way to check the calculation. Completing a sweep satisfies the second row exactly, but changing $y$ disturbs the first row:Thus a zero residual in the last updated row alone does not certify a solution; inspect the complete linear-system residual.

For the same matrix, the stale schedule of the Jacobi method gives the iteration matrixHere its eigenvalues are $\pm 1/\sqrt{12}$, whereas the fresh schedule gives $\rho (H)=\rho (B{)}^2$. This comparison belongs to this particular linear system; it is not a universal relation between the two methods' spectral radii.

The forward error is the difference between a computed result $\hat{x}$ and the exact result $x_{\ast }$ for the original data. Its size can be measured by a norm, or relative to $\Vert x_{\ast }\Vert$ when $x_{\ast }\ne 0$. A linear-system residual measures the equation imbalance instead; the condition number relates that imbalance to the forward error.

An invertible matrix and nonzero diagonal entries make the update well defined, but do not guarantee convergence. For the original linear systemthe Gauss-Seidel method updates $x^{(k+1)}=1-2y^{(k)}$ and $y^{(k+1)}=1-2x^{(k+1)}=-1+4y^{(k)}$. Starting at zero gives $(1,-1)$, then $(3,-5)$. The second-coordinate forward error is $e_y^{(k)}=-4^k/3$, and the iteration matrix $H=\left(\begin{array}{cc}0& -2\\ 0& 4\end{array}\right)$ has spectral radius $4$. The linear system has a unique solution while this particular stationary iterative method for a linear system diverges.

For a stationary iteration $x^{(k+1)}=H{x}^{(k)}+v$, the matrix $H$ propagates the error by $e^{(k+1)}=H{e}^{(k)}$. Convergence from every initial error is equivalent to $\rho (H)<1$.

For a matrix splitting $A=M-N$, the iteration matrix is $H=M^{-1}N=I-M^{-1}A$, where $I$ is the identity matrix. If $A{x}_{\ast }=b$, subtracting its fixed point equation from the update givesThus $b$ determines the fixed point, while $H$ determines how the forward error moves. An eigenvector $v$ with eigenvalue $\lambda$ is an error direction whose coefficient is multiplied by $\lambda$ on every step. A negative real eigenvalue reverses its sign; a zero eigenvalue removes that component in one step.

A linear solver is an algorithm that computes a solution or controlled approximation to a linear system $Ax=b$. For an invertible matrix, the mathematical solution is unique. A direct linear solver usually computes a matrix factorization, such as Pivoted LU decomposition, followed by triangular solves. An iterative method instead updates an approximate vector, often using matrix-vector products and a matrix preconditioner.

A stopping rule should reflect the required accuracy. The linear-system residual $r=b-A\hat{x}$ can be scaled aswhen the denominator is positive, using compatible vector norms and an induced matrix norm. This compares equation error with the input scale; a condition number is still needed to interpret it as a forward error guarantee. In a differential equation computation, the linear solver tolerance must be small enough relative to the intended discretization accuracy. A stable time step alone does not control errors left by an inaccurate inner solve.

Numerical linear algebra develops algorithms for linear systems, least squares, eigenvalue problems and other matrix computations, with attention to roundoff error, cost and storage. A mathematically equivalent formula can behave very differently in floating-point arithmetic: solving through an LU decomposition avoids forming an explicit matrix inverse, while a QR decomposition gives an orthogonal route to least squares.

The central distinction is between accuracy of the computed answer and accuracy of the problem it solves. For $Ax=b$, an approximate vector $\hat{x}$ has linear-system residual $r=b-A\hat{x}$. A small linear-system residual need not imply a small forward error when the condition number is large. For the exact examplethe solution is $x=(1,1{)}^T$, whereas $\hat{x}=(1,0{)}^T$ has Euclidean norm residual $10^{-12}$ but Euclidean norm error $1$. Here the condition number in the spectral norm is $10^{12}$. The lost component barely changes the equations because its coefficient is tiny.

Direct linear solvers such as Gaussian elimination factor a matrix; iterative methods such as conjugate gradient successively improve a vector. For a large sparse matrix, fill-in and matrix preconditioning can matter as much as arithmetic count. The task is to preserve the structure that makes a computation accurate and affordable.

For $0<\omega <2$, the coefficient $\gamma =2/\omega -1$ is positive. Moreoveris a positive-definite matrix: for a nonzero vector $z$, put $y=B^{-1}{A}^{1/2}z\ne 0$; then $z^{T}Gz=y^{T}Dy>0$. Hence every eigenvalue of $M=I-\gamma G$ is strictly below one. Part (d) also gives $M=T^{T}T$, so these eigenvalues are nonnegative. ThereforeThe energy-norm proof of convergence of SOR thus makes the iteration a strict contraction mapping in this norm. Its forward error satisfies $e^{(j+1)}=C(\omega )e^{(j)}$, so it converges geometrically to zero.

Let $\tilde{A}=\hat{U}\hat{\Sigma }{\hat{V}}^T$. Interpret the diagonal of $\hat{\Sigma }$ as nonnegative computed singular values and sort them consistently. If signed diagonal entries are used instead, take their absolute values and absorb signs into the orthonormal columns. Because the two outer factors have the stated orthonormality, these sorted numbers are exactly ${\sigma }_j(\tilde{A})$.

For completeness, the singular value perturbation bound follows from the Eckart-Young-Mirsky theorem. If $B_{j-1}$ is the truncated singular value decomposition of a matrix $B$, the triangle inequality givesExchange $A$ and $B$ to get $|{\sigma }_j(A)-{\sigma }_j(B)|\le \Vert A-B{\Vert }_2$. Applying this to $B=\tilde{A}$ and the given backward error provesThus backward stability of the singular value decomposition gives an absolute forward error bound for every singular value, even a very small one.

No uniform absolute error bound of order $ϵ$ follows. Take the Jordan blockHere $\Vert B{\Vert }_2=1$ and $\Vert E{\Vert }_2=ϵ$. The exact eigenvalues of $B$ are both zero, whereasThese are admissible outputs under a backward error guarantee of size $ϵ$, yet their forward errors are $\sqrt{ϵ}$. This is the square-root splitting of a defective double eigenvalue. More generally, eigenvalue sensitivity depends on the matrix; backward stability alone does not control it. The singular value perturbation bound used in part (ii) has no equally uniform counterpart for the eigenvalues of arbitrary nonsymmetric matrices.

The linear-system residual $r^{(k)}=b-A{x}^{(k)}$ is computable without knowing $x_{\ast }$. The forward error $e^{(k)}=x^{(k)}-x_{\ast }$ is the actual difference from the desired solution. For an invertible matrix,Keeping $A$ fixed, the approximate vector exactly solves $A{x}^{(k)}=b+\delta b$ with $\delta b=-r^{(k)}$. Thus, for $b\ne 0$, $\Vert r^{(k)}\Vert /\Vert b\Vert$ is the relative backward error when only the right-hand side may change. For a vector norm and its induced operator norm,A small relative linear-system residual ensures a small relative forward error when the condition number is moderate; an ill-conditioned matrix can amplify it. The bound follows directly from $e=-A^{-1}r$ and $\Vert b\Vert =\Vert A{x}_{\ast }\Vert \le \Vert A\Vert \Vert x_{\ast }\Vert$.^[ref]

For the Jacobi iteration for a three-variable tridiagonal system, the first linear-system residual is $(1,3,1{)}^T$, while the first forward error is $(-1/2,-1,-1/2{)}^T$. In the Euclidean norm, the relative forward error is $\sqrt{3/28}\approx 0.3273$ and the relative linear-system residual is $\sqrt{11/120}\approx 0.3028$. Its spectral condition number of a positive-definite matrix is $(4+\sqrt{2})/(4-\sqrt{2})\approx 2.094$, so the bound is about $0.6340$. The bound is a guarantee, not an equality. In floating-point arithmetic, recompute the linear-system residual from the current vector when checking termination; a small update alone need not certify a small forward error.

A stationary iterative method for a linear system uses the same iteration matrix at every step:For a specified fixed point $x_{\ast }$ satisfying $x_{\ast }=H{x}_{\ast }+v$, the forward error obeys $e^{(k)}=H^k{e}^{(0)}$. Convergence to that same $x_{\ast }$ for every starting vector is equivalent to $\rho (H)<1$, as the convergence criterion for a stationary iteration proves. Merely having a convergent sequence for every start is weaker: $H=I$, $v=0$ leaves each starting value unchanged.

For a matrix splitting $A=M-N$ with $M$ invertible, choose $H=M^{-1}N$ and $v=M^{-1}b$. Then every solution of $Ax=b$ is a fixed point. The choices of $M$ produce the Jacobi method, Gauss-Seidel method and successive over-relaxation. Fixed coefficients make the error behavior analyzable once, while the vector $b$ determines the desired solution.