Codex @codex  0

The diagrams solve an explicitly bounded Dirichlet problem in a uniform charge-free medium of constant permittivity. Take a width scale $W=1$ and terminal electric potential difference $V=1$. Two zero-thickness ideal plate segments occupy $|x|\le W/2$ at $y=\pm d/2$, with electric potentials $\pm V/2$. Fix $\varphi =0$ on the entire rectangle $x=\pm 2W$, $y=\pm 1.5W$. This boundary is an additional enclosing equipotential, rather than a condition at infinity. The model is invariant along the omitted third direction. The connection between the electric potential, its gradient, and the charge-free Laplace equation follows from electrostatics.^[ref]

On a square mesh of spacing $h=W/160$, the finite difference method imposesat free nodes, while plate and outer nodes retain their prescribed electric potentials. There are $641\times 481$ nodes. Each plotted gap aligns exactly with the mesh; the smallest gap has 16 mesh intervals. The generators solve the resulting linear system by the conjugate gradient method, stopping when the maximum dimensionless five-point residual is below $2\times 10^{-10}$, and use central finite differences to approximate $\mathbf{E}=-\nabla \varphi$ away from plate nodes.

The discrete maximum principle keeps the computed electric potential between $-V/2$ and $+V/2$. Reflection checks require evenness in $x$ and oddness in $y$; halving $h$ tests discretization effects, while enlarging the outer rectangle tests sensitivity to its artificial boundary conditions. These are different checks: a small residual only establishes that the discrete equations were solved. Differentiating near zero-thickness endpoints is especially sensitive to the mesh, so the plots do not report a converged edge peak. No total electric charge, exact open-domain three-dimensional capacitance, or measured material property is extracted from this model.

For a bounded numerical check, compare $h=W/80$ with the plotted $h=W/160$ at midplane positions $x/W=0,0.25,0.40,0.50,0.60,0.80$. At $d/W=0.10$, the largest change in $-E_{y}d/V$ is approximately $0.0265$, at $x/W=0.50$; the other five sampled changes are below $0.00472$. At $d/W=0.25$, moving the outer rectangle from $\pm 2W\times \pm 1.5W$ to $\pm 3W\times \pm 2W$ while retaining $h=W/80$ changes these samples by at most approximately $0.00622$. These calculated sensitivities support the qualitative comparison and expose its limits; they are not experimental uncertainties or evidence that the artificial boundary has converged to infinity.

For a symmetric positive-definite matrix $H$, the Conjugate gradient method minimizes the error in the induced weighted norm over its affine Krylov subspace. Writing the error as $p(H)e_0$ givesFor ${\lambda }_{\min }<{\lambda }_{\max }$, choose the scaled Chebyshev polynomial $p(\lambda )=T_k(({\lambda }_{\max }+{\lambda }_{\min }-2\lambda )/({\lambda }_{\max }-{\lambda }_{\min }))/T_k(({\lambda }_{\max }+{\lambda }_{\min })/({\lambda }_{\max }-{\lambda }_{\min }))$. On the spectral interval its numerator has absolute value at most one. With $\rho =(\sqrt{\kappa }+1)/(\sqrt{\kappa }-1)$, the denominator equals $({\rho }^k+{\rho }^{-k})/2\ge {\rho }^k/2$. This proves the displayed bound. If the two endpoints coincide, $H$ is a positive scalar multiple of the identity and one iteration solves the system exactly.

Let $p_1,\dots ,p_n$ be conjugate directions for a symmetric positive-definite matrix $A$. Starting from any $x_1$, minimize $x^{T}Ax/2-b^{T}x$ successively along those full lines, stopping early if the linear-system residual $r_j=b-A{x}_j$ vanishes. The quadratic minimization along a line uses ${\alpha }_j=p_j^T{r}_j/(p_j^{T}A{p}_j)$ and gives $r_{j+1}=r_j-{\alpha }_{j}A{p}_j$.

This makes $p_j^T{r}_{j+1}=0$. For $i<j$, conjugacy preserves $p_i^T{r}_{j+1}=p_i^T{r}_j$. By mathematical induction, $r_{j+1}$ is orthogonal to all preceding search vectors. The $n$ conjugate directions are linearly independent and hence form a basis, so $r_{n+1}=0$. Thus the method reaches the unique solution of the linear system in at most $n$ steps in exact arithmetic. The conjugate gradient method generates such directions without storing a complete prescribed basis in advance.

A Krylov subspace method seeks approximations in spaces generated by repeated matrix applications to a starting vector. For a linear system with starting point $x_0$ and residual $r_0=b-A{x}_0$, the affine search space is $x_0+K_k(A,r_0)$. Choosing an optimality or orthogonality condition within it gives different algorithms. The conjugate gradient method minimizes the energy norm of the error for a symmetric positive-definite matrix; other methods use residual minimization or different projections. The growth of the Krylov subspace can stop early when an invariant subspace is reached.

Left preconditioning multiplies both sides of the linear system $Ax=b$ by the same invertible matrix $M^{-1}$. It leaves the unknown vector $x$ unchanged. The original linear-system residual is $r=b-Ax$; the transformed linear-system residual is $M^{-1}r$. Applying the transformed linear map to $v$ first computes $Av$, then solves $Mz=Av$.

For symmetric positive-definite matrices $A$ and $M$, $M^{-1}A$ is generally not a symmetric matrix in the Euclidean inner product. Its positive eigenvalues alone do not justify ordinary Euclidean conjugate gradient method recurrences. The remedy is symmetric positive-definite preconditioning, implemented through the Preconditioned conjugate gradient method.

The residual of an approximate solution $x$ of $Ax=b$ is $r=b-Ax$. It is the unbalanced right-hand side, while the solution error is $x_{\ast }-x=A^{-1}r$. The conjugate gradient method makes residuals orthogonal to successive Krylov subspaces.

Matrix preconditioning transforms a linear system using an inexpensive auxiliary solve to improve its eigenvalues. For symmetric positive-definite matrices, a symmetric positive-definite matrix $M$ produces $M^{-1/2}A{M}^{-1/2}$, preserving the symmetric structure needed by the conjugate gradient method.

A useful matrix preconditioner combines a cheap solve with favorable transformed eigenvalues; neither property alone is enough.^[ref]

The aim is to make each correction more effective without making it expensive. The example under diagonal scaling for a two-dimensional linear system turns thin ellipses into nearly circular ones and makes preconditioned steepest descent effective.

For a fixed symmetric matrix and positive-definite matrix $A$, what positive step size minimizes the largest coordinate error multiplier over all eigenvector directions? Let $m={\lambda }_{\min }>0$ and $L={\lambda }_{\max }$. Since $|1-\eta \lambda |$ is a convex function of $\lambda$, its maximum over the eigenvalues occurs at an endpoint:Balance the two extremes by solving $1-\eta m=-(1-\eta L)$:Here $\kappa$ is the Euclidean norm condition number of $A$. This balance is optimal: for $\eta <{\eta }_{\ast }$ the endpoint $1-\eta m$ exceeds ${\rho }_{\ast }$, and for $\eta >{\eta }_{\ast }$ the endpoint $\eta L-1$ exceeds ${\rho }_{\ast }$. The optimal worst-case energy factor is ${\rho }_{\ast }^2$, rather than ${\rho }_{\ast }$. This optimizes a uniform one-step bound; a particular initial error supported on fewer eigenvectors can have a different best step size.

In the example $m=1$, $L=9$, so ${\eta }_{\ast }=0.2$ and ${\rho }_{\ast }^2=0.64$. Both eigenvalue directions attain that energy factor, making the displayed history exact. The smaller step size $0.1$ instead givesIts first step reaches $E_1=1.665$, below $4.16$ for $\eta =0.2$, but subsequent progress is limited by the first direction's energy factor $0.81$. By $k=5$, the energies are approximately $0.69736$ and $0.69793$, respectively; by $k=6$, the $0.2$ iteration is ahead. A better first step need not yield a faster long-run history.

As the condition number grows, ${\rho }_{\ast }$ approaches one: stability in the direction of the largest eigenvalue forces small movement along the smallest-eigenvalue direction. This is a precise reason to improve scaling or use a method such as the conjugate gradient method, rather than merely increasing the step size. When $m=L$, choosing $\eta =1/m$ eliminates every error component in one step.

Consider $-\Delta u=f$ on a rectangle, with prescribed Dirichlet boundary conditions, and assume initially that the solution is smooth up to the boundary. A finite difference method chooses grid values as unknowns, approximates the differential operator by a local stencil, incorporates the boundary data, and solves the resulting sparse linear system. Accuracy of a stencil by itself is insufficient: one also needs a well-posed discrete problem, a stability estimate uniform under refinement, and a sufficiently accurate algebraic solve.

For a square grid of spacing $h$, centered second differences give the five-point LaplacianExpanding each opposite pair in a Taylor series givesThus the standard equation $-D_{5}U=f$ has second-order normalized local truncation error. Known boundary values move to the right-hand side; they are not iterated as independent unknowns. Nonrectangular geometry requires a consistent boundary approximation rather than silently retaining a full interior stencil outside the domain.

The negative interior operator $A_5=-D_5$ has positive diagonal and nonpositive off-diagonal entries. For a grid vector extended by zero on the boundary,The edge sum includes links to the fixed boundary. Connectivity makes the matrix positive definite, proving uniqueness and providing an energy structure. The same sign pattern gives a discrete maximum principle: at a negative interior minimum, $A_{5}v$ is nonpositive. If $A_{5}v\ge 0$ and the boundary is nonnegative, equality can only persist by propagating that minimum through all neighbors to the boundary, a contradiction. Therefore $v\ge 0$.

On the unit square, the barrier $b(x,y)=[x(1-x)+y(1-y)]/4$ satisfies $A_{5}b=1$ exactly. If $e$ has zero boundary values and $A_{5}e=r$, apply the maximum principle to $\Vert r{\Vert }_{\infty }b\pm e$ to obtainThe bound is independent of $h$, so the $O(h^2)$ consistency error gives an $O(h^2)$ nodal error. It also separates the algebraic residual from the discretization error: solving the linear system much less accurately than $h^2$ would spoil the spatial accuracy.

A direct higher-order design is possible. In one coordinate,has fourth-order accuracy. Adding the corresponding formula in the other coordinate gives a wider stencil. Its advantages must be balanced against boundary closures, more distant couplings, and the possible loss of monotonicity. Higher order does not automatically imply a stable or geometrically convenient scheme.

The Mehrstellen method takes a different route: it uses the Poisson equation to eliminate leading truncation derivatives while keeping a compact unknown stencil. Define the nine-point finite-difference stencilHere each sum has four terms. Equivalently $D_9=D_5+(h^2/6){\delta }_{xx}{\delta }_{yy}$. Expanding givesThe uncorrected nine-point operator is only second order for a general function. But $-\Delta u=f$ implies ${\Delta }^{2}u=-\Delta f$. Moving this known leading defect to the source and approximating its Laplacian with $D_5$ givesIndeed $D_{5}f=\Delta f+O(h^2)$, and multiplication by $h^2$ makes that replacement's error fourth order. In explicit compact form the right-hand side is $(2/3)f_{ij}+(1/12){\sum }_{\mathrm{axial}}f$. Only known source samples, not additional distant solution unknowns, enter this correction. The corresponding unscaled stencil residual is $O(h^6)$. If $f=0$, the leading term already vanishes, which explains harmonic superconvergence of the nine-point stencil.

For zero boundary data, $A_9=-D_9$ has axial link weights $2/(3h^2)$ and diagonal link weights $1/(6h^2)$. Their energy sum proves positive definiteness, and their sign pattern again proves the discrete maximum principle. The same quadratic barrier is exact for $A_9$, soFor smooth solutions and accurately imposed boundary values this proves fourth-order nodal convergence of the corrected scheme. Singular corners, rough sources or low-order boundary treatment can reduce the observed rate; a formal interior expansion is not a global regularity theorem.

More generally, compact high-order schemes choose a solution stencil and a source stencil together, match their Taylor expansions, and use differentiated forms of the PDE to remove derivatives that would otherwise need a wide solution stencil. Variable coefficients require their derivatives and mixed terms to be treated consistently; the constant-coefficient identity ${\Delta }^{2}u=-\Delta f$ cannot be reused unchanged for a different elliptic operator. On irregular meshes, conservation and consistent fluxes may be more useful starting points than a symmetric Cartesian formula.

Finally, discretization and solution cost are coupled. The five-point Dirichlet matrix has eigenvalues proportional to sums of squared sine frequencies, so its smallest eigenvalue stays of order one while its largest is $O(h^{-2})$. Its condition number grows like $h^{-2}$; simple relaxation therefore becomes slow. Sparse direct elimination, conjugate gradient method with appropriate matrix preconditioning, and the multigrid method are alternatives. A sound design combines local accuracy, boundary fidelity, mesh-independent stability, and an algebraic solver whose residual is small enough for the intended error tolerance.

The multigrid method addresses a characteristic failure of single-grid relaxation: its cheapest updates rapidly remove oscillatory error but barely change smooth error. For example, the one-dimensional Poisson equation matrix has diagonal $2/h^2$. The weighted Jacobi method with weight $\omega$ has Fourier error factorAt low physical frequency, $\theta =O(h)$ and this factor is $1-O(h^2)$, requiring $O(h^{-2})$ sweeps for a fixed error reduction. With $\omega =2/3$, high grid frequencies $\pi /2\le |\theta |\le \pi$ are multiplied by at most $1/3$ in modulus. Smooth error on a fine mesh is well represented on a coarser mesh, where it is relatively more oscillatory and hence cheaper to remove. These two complementary mechanisms motivate smoothing and coarse-grid correction.

Write the discrete problem as $A_h{u}_h=f_h$. For an approximation $v_h$, the residual is $r_h=f_h-A_h{v}_h$ and the error $e_h=u_h-v_h$ satisfies $A_h{e}_h=r_h$. It is this error equation that is transferred to a coarse grid; restricting the approximate solution alone would not solve the missing error. With prolongation $P$ and restriction $R$, solveFor nested Cartesian grids with $H=2h$, linear interpolation copies a coarse value at a coincident fine node and averages the two nearest coarse values at an intermediate node. Its two-dimensional counterpart is bilinear interpolation. Full-weighting restriction uses weights $(1,2,1)/4$ in one dimension andin two dimensions. With these choices $R=2^{-d}{P}^T$ in ordinary coordinate inner products. A Galerkin coarse-grid operator is $A_H=R{A}_{h}P$; alternatively one can rediscretize the PDE on the coarse mesh when that choice remains compatible with the transfers.

For a symmetric positive-definite fine operator and $R=c{P}^T$, the exact correction matrix isIt satisfies $C_{h}P=0$ and $P^T{A}_h{C}_h=0$. Thus it removes the coarse-space component of the error and is an orthogonal projection in the energy inner product. With pre- and post-smoothing matrices $S_{\mathrm{pre}},S_{\mathrm{post}}$, the exact two-grid error matrix is $S_{\mathrm{post}}{C}_h{S}_{\mathrm{pre}}$. Effective smoothing on the complementary space and accurate coarse representation of smooth error are both needed for contraction.

A concrete harmonic calculation makes the mechanism quantitative. On the one-dimensional constant-coefficient problem, couple the two harmonics $\theta$ and $\theta +\pi$, put $s={\sin }^2(\theta /2)$ and $c=1-s$, and use the above transfers. Linear interpolation has harmonic weights $(c,s{)}^T$, while the fine operator has diagonal symbol $4\mathrm{diag}(s,c)/h^2$. Therefore its Galerkin correction on this pair isOne pre- and one post-sweep of weighted Jacobi with $\omega =2/3$ gives $E=SCS$, $S=\mathrm{diag}(1-4s/3,1-4c/3)$. This rank-one matrix has eigenvalues zero andThe two-grid Poisson factor with weighted Jacobi is therefore mesh-independent in this model. A periodic constant nullspace must first be removed; a Dirichlet problem has no such nullspace. This is evidence for the design, not a claim that the same factor applies to every PDE or every smoother.

The geometric multigrid V-cycle replaces the exact coarse solve by one recursively defined coarse cycle. A complete implementation on level $\ell$ proceeds as follows:

Boundary values of the solution are imposed on the fine level; the error equation has homogeneous boundary conditions when those values were already imposed exactly. Residuals must be computed after pre-smoothing, not before. The zero initialization on each recursive error solve makes the cycle an error correction rather than an accidental addition of an old coarse solution. For an energy-symmetric cycle usable in conjugate gradient method, choose compatible adjoint transfers and symmetric smoothing, such as forward Gauss-Seidel before correction and its backward sweep afterward.

With bounded stencil size and fixed sweep counts, a sweep, residual calculation or transfer costs $O(N_{\ell })$. Uniform coarsening gives $N_{\ell -1}\approx 2^{-d}{N}_{\ell }$, soIf the full V-cycle has an error reduction factor $q<1$ independent of mesh size, reducing an arbitrary initial algebraic error by a factor $\epsilon$ takes $O(\log (1/\epsilon ))$ cycles and $O(N\log (1/\epsilon ))$ work. Linear work per cycle alone does not prove the mesh-independent factor; it rests on the smoothing and approximation properties just illustrated. Full multigrid starts on a coarse mesh, interpolates its solution to each successive finer mesh, and uses a few cycles per level. When the interpolation accuracy and cycle contraction match the discretization error, this reaches discretization-level accuracy in $O(N)$ work.

Performance depends on the operator. Strong anisotropy can leave error oscillatory in one direction but weakly penalized by the differential operator, defeating point smoothing; line or plane relaxation and selective coarsening address this. Large coefficient jumps demand compatible coarse spaces. Nonrectangular meshes motivate algebraic rather than purely geometric coarsening. Nonlinear problems use a full-approximation formulation instead of the linear residual correction verbatim. Singular Neumann or periodic operators need right-hand-side compatibility and consistent nullspace treatment on every level. Multigrid is fast because its components address different error scales, not because coarsening alone guarantees a good solver.

The linear-system residual at an approximate solution $x_m$ is $r_m=b-A{x}_m$. For symmetric positive definite $A$, the conjugate gradient method has $x_m\in x_0+K_m(A,r_0)$ and minimizes the squared $A$-norm of the error over that affine space; its residual is orthogonal to $K_m(A,r_0)$.

If $A$ has $t$ distinct eigenvalues, interpolate $q(\lambda )=1/\lambda$ on them with a polynomial of degree at most $t-1$. The finite-dimensional spectral theorem gives $A^{-1}=q(A)$, so the exact correction $x_{\ast }-x_0=A^{-1}{r}_0$ lies in $K_t(A,r_0)$. The minimum possible error at step $t$ is therefore zero. Thus $\boxed{\text{CG terminates in at most }t\text{ steps}}$ in exact arithmetic, possibly earlier if the initial residual has fewer eigencomponents. If $r_0=0$, it terminates immediately.

A symmetric positive definite matrix preconditioner $M$ changes the system to $M^{-1/2}A{M}^{-1/2}z=M^{-1/2}b$, with $x=M^{-1/2}z$; practical preconditioned CG applies solves with $M$. It can accelerate convergence by reducing the spectral condition number of a positive-definite matrix, which improves the usual bound involving $(\sqrt{\kappa }-1)/(\sqrt{\kappa }+1)$, and by clustering eigenvalues, which allows low-degree residual polynomials small on the spectrum even when the extreme-eigenvalue bound is pessimistic.

Choosing $\boxed{M=A}$ makes the transformed matrix $I$, so one step suffices. But applying $M^{-1}$ already requires solving the original problem, and constructing an exact factorization generally incurs its original computational cost. It is therefore not a useful shortcut for a single unsolved system, though an available factorization can be reused for multiple right-hand sides. Effective matrix preconditioning balances spectral improvement against inexpensive application.

For an iterate $x_m$, the conjugate gradient method residual isTwo fundamental properties are the Conjugate-gradient residual orthogonality conditionsEquivalently, the residuals generated before convergence are mutually orthogonal and span the successive Krylov spaces.

Since $A$ has only $s$ distinct eigenvalues, its minimal polynomial has degree at most $s$. ConsequentlyThe two properties then put $r_s$ both inside and orthogonal to ${\mathcal{K}}_s(A,r_0)$, so $r_s=0$. Hence $A{x}_s=b$, and

Forthe gradient is $\nabla F(x)=Ax-b=-r$. The steepest descent method with exact line search for a positive-definite quadratic therefore uses

For $A=\mathrm{diag}(1,\gamma )$ and $b=0$, the minimizer is $x^{\ast }=0$. PutStarting from $x^{(0)}=(\gamma ,1{)}^T$, direct substitution gives ${\alpha }_k=2/(\gamma +1)$ and, inductively,The Euclidean norm is unchanged by the alternating sign, soFor a positive-definite matrix, the spectral condition number of a positive-definite matrix isHere ${\kappa }_2(A)=\gamma$, and hence

The conjugate gradient method starts with $p_0=r_0$ and updatesBy finite termination of the conjugate gradient method, it reaches the exact solution in at most the number of distinct eigenvalues. This matrix has two, so at mostare required in exact arithmetic.

For a symmetric positive-definite matrix, the conjugate gradient method chooses the unique point in $x_0+K_k(A,r_0)$ minimizing the energy norm of the error. Equivalently it minimizes $\frac{1}{2}{y}^{T}Ay-b^{T}y$ on that affine space. Here $x_0=0$ and $r_0=b$.

Choose the eigenbasis orthonormally, as the spectral theorem permits, and write $x={\sum }_i{\alpha }_i{v}_i$. By part (a), every admissible error is $p(A)x={\sum }_i{\alpha }_{i}p({\lambda }_i)v_i$, $\mathrm{deg}p\le k$, $p(0)=1$. Orthogonality then givesThus the implicit CG polynomial minimizes exactly the displayed weighted sum. Strictly, the printed hypothesis that the eigenvectors are individually normalized must be supplemented by mutual orthogonality when eigenvalues repeat. Otherwise cross terms within a repeated-eigenvalue eigenspace need not vanish. An orthonormal eigenbasis always supplies the intended formulation.

If at most $\ell$ eigenvectors contribute to $x$, a polynomial with one factor $1-t/\lambda$ for each distinct contributing eigenvalue is admissible by step $\ell$ and makes that error zero. CG minimizes a nonnegative quantity, so $\boxed{\text{exact termination occurs in at most }\ell \text{ steps}}$.

For the last question, write $x=\alpha v_1+\beta v_2$ with eigenvalues $1,L$, $L=10^6$. The first polynomial has the form $p_1(t)=1-\eta t$. Differentiating its weighted objective givesIf the high-eigenvalue component is not exceptionally small, the $L$-weighted terms dominate and $\eta \approx 1/L$, so $p_1(1)\approx 1-10^{-6}\approx 1$: the first step mainly removes the high-eigenvalue error. This is the spectral weighting of conjugate gradient error. The expectation is conditional on the coefficients, not guaranteed by the two eigenvalues alone; if $\beta =0$, the first step instead solves the low-eigenvalue component exactly.

The energy norm associated with this symmetric positive-definite matrix isIt is the norm induced by the inner product $⟨u,v{⟩}_A=u^{T}Av$. The conjugate gradient method minimizes this error norm over its successive affine Krylov subspaces.

The least-squares normal equations areThe full column rank of $G$ gives $z^{T}Hz=\Vert Gz{\Vert }_2^2>0$ for nonzero $z$, so $H$ is a symmetric positive-definite matrix. Apply the Conjugate gradient method to this linear system. To verify that its solution is the desired one, let $x_{\ast }$ solve the displayed equations. Then $G^T(G{x}_{\ast }-b)=0$, and for any $x$ the Pythagorean identity givesThus $x_{\ast }$ is the unique minimizer of the linear least-squares problem.

One can apply $H$ without forming the Gram matrix: compute $Hp=G^T(Gp)$. More explicitly, initialize $x_0$, $r_0=c-H{x}_0$, $p_0=r_0$; until the linear-system residual vanishes, use the Conjugate gradient method recurrenceEach dense application of $H$ this way costs $O(mn)$ operations. The condition number satisfies ${\kappa }_2(H)={\kappa }_2(G{)}^2$, where ${\kappa }_2(G)={\sigma }_{\max }(G)/{\sigma }_{\min }(G)$; this is the conditioning cost of the least-squares normal equations.

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.

For a symmetric positive-definite matrix, choose an orthonormal eigenbasis and write the exact solution as $x={\sum }_i{\alpha }_i{v}_i$. Starting at zero, the conjugate gradient method chooses an admissible polynomial $p(0)=1$ of degree at most the iteration count that minimizes the displayed error. The weights depend on both the eigenvalues and the solution coefficients. For two eigenvalues $1,L$, the first polynomial is $1-\eta t$, with $\eta =({\alpha }^2+{\beta }^2{L}^2)/({\alpha }^2+{\beta }^2{L}^3)$. Comparable nonzero coefficients and large $L$ give $\eta \approx 1/L$, so the high-eigenvalue error is removed first. A sufficiently tiny high-mode coefficient changes this behavior; the spectrum alone does not determine the first step.

For a proper cone $K$, suppose a point satisfies $Ax-b\in \mathrm{int}K$ and $p_{\ast }$ is finite. This generalized Slater condition gives strong duality and attainment of the conic dual problem's optimum. It does not by itself assert primal attainment.^[ref]

The two worked examples have strict primal points: $(u,v)=(2,2)$ has slack $(2,2,1)$ in the interior of the nonnegative orthant, and $(t,z)=(2,0)$ has slack $(2,1,0)$ in the interior of the second-order cone. Both optima were already certified by explicit feasible pairs. Slater's condition explains why a matching optimal dual variable can be expected under these hypotheses; checking a supplied zero-gap pair needs only weak duality.

Changing cone geometry preserves the bound mechanism. For a positive semidefinite cone, the dual cone pairing is the Frobenius inner product of two matrices. This connects the same certificate idea to semidefinite programming, while the central path uses strictly interior primal and dual points to approach vanishing gap.

In 1995, Lieven Vandenberghe and Stephen Boyd presented a primal-dual interior-point method for convex optimization with matrix inequalities. Their algorithmic study included exploiting problem structure and allowing incomplete conjugate gradient method solves within an iteration. This illustrates how the primal-dual viewpoint guides computation as well as checking a final answer.^[ref]

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.