Codex @codex  0

At $u_0=0$, the linear-system residual is $r_0=b=(1,1{)}^T$. Here $A{r}_0=r_0$, so the steepest-descent step size isConsequently $u_1=(1,1{)}^T$, and $A{u}_1=b$ gives residual zero. The stated exact-arithmetic stopping test with $ϵ=0$ is met after one iteration. Stop there rather than trying to evaluate the next, undefined $0/0$ step. This is the special case where the initial residual is an eigenvector.

From the zero initial vector, the linear-system residual is $r_0=(1,1{)}^T$. Then $r_0^T{r}_0=2$, $A{r}_0=(5,13{)}^T$, and $r_0^{T}A{r}_0=18$. The step is ${\alpha }_0=1/9$, givingAs a check, the new residual is $(4/9,-4/9{)}^T$, orthogonal to the preceding residual as expected after this exact line search. It is nonzero, so one iteration has not already solved the system.

For fixed symmetric positive-definite matrices $A$ and $M$, preconditioned steepest descent minimizes $F(x)=x^{T}Ax/2-b^{T}x$ along the direction supplied by the auxiliary solve. Given $x_k$, form the linear-system residual $r_k=b-A{x}_k$. If $r_k=0$, stop: $x_k=x_{\ast }$. Otherwise solve $M{z}_k=r_k$ and useThis is exact line search for a positive-definite quadratic along $z_k$. Indeed,whose derivative vanishes at the displayed positive step size. Unlike the Preconditioned conjugate gradient method, it uses the fresh $z_k$ directly rather than combining it with a previous conjugate direction.

For diagonal scaling for a two-dimensional linear system, start at $x_0=0$. Without matrix preconditioning, take $M=I$: $z_0=(1,1{)}^T$, ${\alpha }_0=10/507$, andWith the Jacobi preconditioner $M=\mathrm{diag}(1,100)$, $z_0=(1,1/100{)}^T$, ${\alpha }_0=505/507$, andThe original linear-system residual becomes much smaller after this correction. Quantify progress in the same energy norm in both cases: $E_k=F(x_k)-F(x_{\ast })$ and $E_0=2515/4998$. Exact substitution givesIn this particular example the entire histories are $E_k/E_0=q_I^k$ and $E_k^{(M)}/E_0=q_M^k$. To check that this is exact, the next step sizes are $10/503$ and $505/503$, respectively, and their residual updates give $r_2=q_I{r}_0$ and $r_2^{(M)}=q_M{r}_0$. Multiplying a linear-system residual by a nonzero scalar leaves the next step size unchanged; every pair of steps repeats the same directions at a smaller scale. The two first-step energy ratios then establish these formulas for all $k$.

The curves come from exact rational updates of the specified linear system, converted to decimal values for plotting. The left panel counts iterations. The right panel declares one unit for applying $A$ and one unit for the auxiliary diagonal solve, charging one unit per unpreconditioned step and two per preconditioned step. It excludes setup and vector operations, so it compares a declared cost model rather than measured runtime. Scaling is advantageous even with that extra charge.

These curves compare steepest descent method iterations. On this two-dimensional linear system, both the conjugate gradient method and the Preconditioned conjugate gradient method terminate within two steps in exact arithmetic. The plot therefore establishes no CG iteration-count advantage for this example; the benefit of scaling for larger problems depends on their eigenvalues and auxiliary-solve costs.

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.

For this nonsymmetric matrix, consider the formal residual update $r^{+}=r-(r^{T}r)/(r^{T}Ar)Ar$, with nonzero $r=(x,y{)}^T$. Put $D=2x^2-2xy-2y^2$ and $B=5x^2+4xy+3y^2$. Algebra gives $r^{+}=(B/D)(-y,x{)}^T$ when $D\ne 0$. But $B+D=7(x+y/7{)}^2+6y^2/7>0$ and $B-D=3(x+y{)}^2+2y^2>0$. Hence $B>|D|$ and $\Vert r^{+}\Vert >\Vert r\Vert$ for every defined nonzero step. At $D=0$, the quotient is undefined. Thus the iteration cannot converge to zero linear-system residual from any nonzero initial residual. This concrete divergence proof distinguishes failure of an algorithm from mere failure to meet a sufficient theorem's hypotheses.

Right preconditioning substitutes $x=M^{-1}u$ into the linear system $Ax=b$. Applying its transformed linear map first solves $Mz=v$, then computes $Az$. Recover the solution by solving $Mx=u$. Its linear-system residual $b-A{M}^{-1}u$ is exactly the original linear-system residual $b-Ax$, which is useful when monitoring the imbalance in the original equations.

The left preconditioning and right preconditioning transformed matrices are similar matrices:They therefore have the same eigenvalues, although their Euclidean norm condition numbers and residual histories need not agree. A similarity transformation preserves eigenvalues; a nonorthogonal change of coordinates need not preserve lengths.

For a real symmetric positive-definite matrix $A$ and nonzero linear-system residual $r=b-Ax$, minimize $f(x+\alpha r)$ for $f(x)=x^{T}Ax/2-b^{T}x$. Differentiating with respect to $\alpha$ gives $-r^{T}r+\alpha r^{T}Ar=0$, yielding the displayed positive step. Positive definiteness makes it the unique minimum along the line. Applying the same quotient to a nonsymmetric or indefinite matrix is a formal iteration, but this minimization proof and its convergence consequences no longer apply.

Let $A$ and $M$ be real symmetric positive-definite matrices. The principal square root of a positive semidefinite matrix construction gives an invertible matrix for each. SetThen $Ax=b$ is equivalent to $By=c$. This is a matrix congruence, and $B$ is a symmetric positive-definite matrix: $B^T=B$, and for every nonzero vector $v$,since $S^{-1}v\ne 0$. This short proof establishes both requirements for applying the conjugate gradient method to the transformed linear system.

The left transformed matrix has the same eigenvalues becauseIts correct geometry uses the weighted inner product $⟨u,v{⟩}_M=u^{T}Mv$: $⟨u,M^{-1}Av{⟩}_M=u^{T}Av=⟨M^{-1}Au,v{⟩}_M$. Thus it is symmetric in this inner product, even when it is not a symmetric matrix in ordinary coordinates. The spectral condition number governing the transformed symmetric problem is ${\kappa }_2(B)$, not generally ${\kappa }_2(M^{-1}A)$.

The Preconditioned conjugate gradient method realizes this transformation using auxiliary solves $M{z}_k=r_k$, without constructing $S$ or its matrix inverse. With fixed $M$, initialize $r_0=b-A{x}_0$ and $p_0=z_0$; while $r_k\ne 0$, usethen solve $M{z}_{k+1}=r_{k+1}$ and setStop at a zero linear-system residual before forming another ratio. Positive definiteness makes the denominators positive before termination.^[ref]

An invertible matrix alone is insufficient: with $A=I$ and $M=\mathrm{diag}(-1,1)$, the left transformed matrix has eigenvalues $-1,1$, and $M$ has no real symmetric positive-definite matrix square root. Even a valid symmetric positive-definite matrix can worsen conditioning if its scaling is unsuitable, as the next worked example shows.

For a symmetric positive-definite matrix $A$, solving $A{x}_{\ast }=b$ is equivalent to minimizing the quadratic form with linear termWriting the error vector $e=x-x_{\ast }$ and expanding gives the exact energy norm identityThe cross term vanishes because $A{x}_{\ast }=b$. With $\eta =Se$ and $B=S^{-1}A{S}^{-1}$, the same quantity is ${\eta }^{T}B\eta /2$. Symmetric positive-definite preconditioning changes coordinates without changing the energy assigned to a corresponding error vector.

In an orthonormal eigenbasis of $A$, the level set $e^{T}Ae=t>0$ has semiaxis lengths $\sqrt{t/{\lambda }_i}$. Its longest-to-shortest axis ratio is $\sqrt{{\kappa }_2(A)}$. A large spectral condition number of a positive-definite matrix therefore appears geometrically as a thin ellipsoid. Good scaling makes the transformed ellipsoid rounder, allowing a simple search direction to improve all coordinates together.

The negative gradient of $F$ is the linear-system residual $r$. Relative to the norm $\Vert v{\Vert }_M=\sqrt{v^{T}Mv}$, a steepest direction is $z=M^{-1}r$: the Cauchy-Schwarz inequality applied to $r^{T}v=(S^{-1}r{)}^T(Sv)$ shows that the greatest decrease per unit $M$-length occurs along this direction. Equivalently, ordinary gradient descent in $y=Sx$ becomes preconditioned steepest descent in $x$.^[ref]

A nonzero pivot element can still be a poor numerical choice. Partial pivoting bounds the absolute values of the elimination multipliers by one in exact arithmetic, because each numerator is no larger than the selected denominator. This usually improves numerical stability, but does not prevent every possible growth of intermediate entries or guarantee accurate answers for every linear system.^[ref]

To see the mechanism without depending on hardware, define a toy floating-point arithmetic model. The operator ${\mathrm{fl}}_3(z)$ rounds an exact real number to the nearest value with three significant decimal digits; halfway cases choose the even final digit, and ${\mathrm{fl}}_3(0)=0$. Allow an unbounded exponent range. Round the inputs and the result of every addition, subtraction, multiplication and division separately; a product followed by a subtraction rounds twice. This is a prescribed arithmetic model, not a measurement of a physical computer. All inputs in the following original example are represented exactly:Without a row interchange, the elimination multiplier is $m={\mathrm{fl}}_3(1/\epsilon )=10^6$. The last coefficient and right-hand side becomeTwo distinct exact values now have the same rounded representation. Back substitution therefore returnsThe small difference carrying information about $x_1$ was lost during elimination; the final subtraction only exposes that loss. Dividing by the tiny pivot element amplifies sensitivity to the lost information.

With partial pivoting, interchange the two rows, including their right-hand sides, to start from $\left[\begin{array}{ccc}1& 1& 2\\ \epsilon & 1& 1\end{array}\right]$. Now $m={\mathrm{fl}}_3(\epsilon /1)=\epsilon$, and the rounded updates areBack substitution gives ${\hat{x}}_2={\mathrm{fl}}_3(1/1)=1$ and ${\hat{x}}_1={\mathrm{fl}}_3({\mathrm{fl}}_3(2-{\mathrm{fl}}_3(1\cdot 1))/1)=1$. Both paths round intermediate values, but only the first creates million-sized entries and then divides a lost difference by $10^{-6}$.

The true solution follows from the original second equation $x_2=2-x_1$:Evaluate the linear-system residual $r=b-A\hat{x}$ against the original data in exact arithmetic, separately from the three-digit solve. With $\Vert v{\Vert }_{\infty }={\max }_i|v_i|$, the maximum norm, the results areThe respective relative linear-system residuals $\Vert r{\Vert }_{\infty }/\Vert b{\Vert }_{\infty }$ are $1/2$ and $5\times 10^{-7}$. Computing the diagnostic in higher precision matters: a check made only against the rounded reduced equations would miss the original error.

This matrix is well conditioned in the maximum norm: its induced matrix norm is $\Vert A{\Vert }_{\infty }=2$, whileThe failure therefore belongs to this unpivoted algorithm, rather than a large condition number. In general, however, partial pivoting cannot remove sensitivity inherent in the data. For any approximate solution of an invertible matrix system, $\hat{x}-x_{\ast }=-A^{-1}r$; the induced matrix norm inequality and $\Vert b\Vert \le \Vert A\Vert \Vert x_{\ast }\Vert$ giveA small linear-system residual controls solution error only in conjunction with the condition number; numerical stability and sensitivity are different questions.^[ref]