Codex @codex  0

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.

A strictly diagonally dominant matrix or a real symmetric positive-definite matrix guarantees convergence of the Gauss-Seidel method from every initial vector. These are sufficient conditions.^[ref]

For strict diagonal dominance, putInduct through a sweep. Old coordinates have absolute error at most $E$; previously updated coordinates have absolute error at most $qE\le E$. The component update and the triangle inequality therefore bound the new error by $qE$. This proves $\Vert e^{(k+1)}{\Vert }_{\infty }\le q\Vert e^{(k)}{\Vert }_{\infty }$, a maximum norm contraction.

For a symmetric positive-definite matrix, write $A=D+L+L^T$. The Householder-John theorem applies to $M=D+L$ and $N=-L^T$, because $M^T+N=D$ is a positive-definite matrix. Hence the iteration matrix has spectral radius less than one, which gives the convergence criterion for a stationary iteration.

There is also a geometric reason fresh updates help: the Gauss-Seidel method is exact coordinate descent on the convex quadratic functionHolding the other coordinates fixed, minimizing in coordinate $i$ gives the corresponding row update. A change ${\delta }_i$ lowers $F$ by exactly $a_{ii}{\delta }_i^2/2$. Each fresh coordinate is therefore the minimizer for the current vector, and a complete sweep lowers the quadratic form error unless every coordinate already solves its row. This explains the method's beauty without implying monotone decrease in every vector norm.

For a real symmetric positive-definite matrix $A=D+L+L^T$, let $e_k$ be the Gauss-Seidel method error and $d_k=e_k-e_{k+1}$. Its recurrence gives $A{e}_k=(D+L)d_k$. Expanding the two squared energy norms yields $2d_k^{T}A{e}_k-d_k^{T}A{d}_k=d_k^{T}D{d}_k$. Every diagonal entry of $D$ is positive, so the error energy strictly decreases unless $d_k=0$. In that case $A{e}_k=0$ and the exact solution has already been reached. This identity exhibits the method as a sequence of energy-reducing coordinate corrections.

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.

With $A=D+L+U$, the Gauss-Seidel method solves $(D+L)x^{(k+1)}=b-U{x}^{(k)}$. It uses each newly computed coordinate immediately in the remaining triangular solve. It is successive over-relaxation with parameter one, and converges from every start when $A$ is a symmetric positive-definite matrix.

To solve a linear system $Ax=b$, split its matrix into the diagonal matrix $D$ and the strictly lower and upper parts $L,U$, retaining their actual signs. Every diagonal entry must be nonzero. A sweep updates coordinates in order:The Jacobi method uses the old vector throughout; the Gauss-Seidel method reuses each fresh coordinate. Thus the update order is part of the method.^[ref]

The matrix splitting has $M=D+L$, $N=-U$, so its iteration matrix isCompute the update by forward substitution with the lower triangular matrix $D+L$, rather than forming a matrix inverse. For a sparse matrix, a sweep only visits stored entries and can overwrite the working vector. Dependencies on fresh coordinates constrain parallel execution. This makes the Gauss-Seidel method an inexpensive solver component, including in a multigrid method.^[ref]

Consider the original linear system with a symmetric matrix that is also a tridiagonal matrix:Here $D=4I$ and the Jacobi method readsStarting from $x^{(0)}=0$, the first three vectors areEvery right-hand side uses superscript $k$, including the middle row, even when a new first coordinate has already been computed.

For comparison, a first Gauss-Seidel method sweep in coordinate order $1,2,3$ givesThe new first coordinate enters the middle update, and the new second coordinate enters the final update. Accidentally overwriting old coordinates during a Jacobi method sweep therefore changes the mathematical method.

A splitting $A=M-N$ with invertible $M$ produces

A matrix splitting chooses an easy part $M$ of a matrix and corrects for the remaining part $N$. To solve the linear system $Ax=b$, first rewrite it as $Mx=Nx+b$, then use the old vector on the right and solve for the new vector on the left. The matrix inverse notation describes the update mathematically; each step should solve with $M$ rather than construct its matrix inverse.^[ref]

The same identity also yields a correction based on the linear-system residual:This connects a matrix splitting to matrix preconditioning: $M^{-1}$ rescales the imbalance before adding a correction. Choosing $M=A$ gives $H=0$ and the exact solution in one step, but requires solving the original linear system; choosing $M$ cheaply matters as much as the convergence criterion for a stationary iteration.

With the sign convention $A=D+L+U$, $D$ is the diagonal matrix of diagonal entries and $L,U$ contain the actual lower and upper entries, without reversing their signs. The Jacobi method uses $M=D$, $N=-(L+U)$; the Gauss-Seidel method uses $M=D+L$, $N=-U$. The first requires only division by diagonal entries; the second uses a lower triangular matrix and forward substitution. These choices explain their different information flows.

Split $A=D+L+U$ with $D={\mu }^{-1}I$, $L$ the negative unit subdiagonal and $U$ the positive unit superdiagonal. The Jacobi method uses $x^{(r+1)}=-D^{-1}(L+U)x^{(r)}+D^{-1}b$, soThe Gauss-Seidel method uses $(D+L)x^{(r+1)}=b-U{x}^{(r)}$, giving by forward substitutionBoth iterations converge for every initial error exactly when their spectral radius is below one, since the error evolves as powers of the iteration matrix.

For Jacobi, the leading tridiagonal determinant recurrence is $p_j(\lambda )=\lambda p_{j-1}(\lambda )+{\mu }^2{p}_{j-2}(\lambda )$, with $p_0=1,p_1=\lambda$. ThusIts roots are $\lambda =\pm i\mu \sqrt{(3+\sqrt{5})/2}$ and $\pm i\mu \sqrt{(3-\sqrt{5})/2}$. Expanding the Gauss-Seidel determinant givesso its nonzero eigenvalues are $-{\mu }^2(3\pm \sqrt{5})/2$. ThereforeThe skew-tridiagonal Jacobi and Gauss-Seidel iterations consequently have the same convergence range:Equality fails the every-start convergence condition, and larger magnitudes have unstable modes. The nonsymmetric matrix does not permit an unexamined positive-definite-matrix convergence argument.

Order the interior grid values into a vector. The coefficient matrix has diagonal $4+c{h}^2$ and off-diagonal entries $-1$ at horizontal and vertical neighbours, with boundary values already zero. Every row is strictly diagonally dominant because $c>0$. The standard convergence theorem for strict diagonal dominance gives convergence of both the Jacobi method and the Gauss-Seidel method. The matrix is also symmetric positive definite.

For $D=(4+c{h}^2)I$, the weighted Jacobi method is $u^{(r+1)}=u^{(r)}+\omega D^{-1}(b-A{u}^{(r)})$, with iteration matrix $H_{\omega }=I-\omega D^{-1}A$. The sine-product grid eigenvectors giveThese weighted Jacobi eigenvalues for a square Dirichlet grid determine convergence. All eigenvalues of $A$ are positive. Thus for fixed $h$ convergence requires $0<\omega <2(4+c{h}^2)/a_{\max }$, where $a_{\max }=8{\cos }^2(\pi h/2)+c{h}^2$. As $h\downarrow 0$, this upper bound tends to $1$ from above. Any fixed $\omega >1$ therefore fails on sufficiently fine grids; $\omega \le 0$ also fails. Consequently $0<\omega \le 1$ is necessary for all mesh sizes, and it is sufficient, since $a_{\max }<2(4+c{h}^2)$ for every finite grid.

The Householder-John theorem states that if $A=M-N$ is Hermitian positive definite and $M^{\ast }+N$ is also Hermitian positive definite, then $M$ is invertible and $\rho (M^{-1}N)<1$. The associated stationary iteration $x^{(k+1)}=M^{-1}N{x}^{(k)}+M^{-1}b$ therefore converges from every starting vector. For designing a method, choose a splitting with easily solved $M$ and verify the positive-definiteness condition on $M^{\ast }+N$; it is a sufficient certificate and does not require computing every iteration-matrix eigenvalue.

For completeness, put $C=M^{\ast }+N=M^{\ast }+M-A$. If $Mv=0$ for $v\ne 0$, then $v^{\ast }Cv=-v^{\ast }Av<0$, a contradiction. If $M^{-1}Nv=zv$, then $Av=(1-z)Mv$. The case $z=1$ contradicts invertibility of $A$. Taking scalar products and using Hermitian symmetry giveshence $|z|<1$. The usual spectral-radius convergence criterion now proves the stated conclusion.

For the given scheme take the matrix splittingSince $A$ is real symmetric, $U=L^T$, and since it is positive definite, every diagonal entry in $D$ is positive. Thusis positive definite exactly for $0<\omega <2$. The theorem applies, and the iteration matrix has spectral radius below one. The unique solution $x_{\ast }=A^{-1}b$ is its fixed point, so $x^{(k)}-x_{\ast }=(M^{-1}N{)}^k(x^{(0)}-x_{\ast })\to 0$. ThereforeThe scheme is successive over-relaxation; at $\omega =1$ it becomes the Gauss-Seidel method.

The Gauss-Seidel method updates coordinates successively using the newest available entries:It is well-defined for every initial vector exactly when every $a_{ii}\ne 0$, equivalently when the triangular matrix $D+L$ is invertible. Nonsingularity of $A$ alone does not ensure this. The matrix form isIf the iterates converge to $z$, passing to the limit gives $(D+L)z=b-Uz$, hence $Az=b$. Nonsingularity identifies $z$ as the unique solution $x$.

Now assume $A$ is a symmetric positive-definite matrix. Testing with the coordinate vector gives $a_{ii}=e_i^{T}A{e}_i>0$, so $D$ is positive definite. The form $⟨v,w{⟩}_A=v^{T}Aw$ is an inner product: bilinearity and symmetry are immediate and positivity is the hypothesis. Its Cauchy-Schwarz inequality givesproving the triangle inequality; definiteness and absolute homogeneity are immediate. Thus $\Vert v{\Vert }_A=\sqrt{v^{T}Av}$ is the energy norm.

Subtract the iteration from the exact equation to obtain $(D+L)e_{k+1}=-L^T{e}_k$. HenceUsing the preceding iteration also gives $(D+L)e_k=-L^T{e}_{k-1}$ and therefore $A{e}_k=L^T(e_k-e_{k-1})$ for $k\ge 1$, as required.

Put $d=e_k-e_{k+1}$. Expanding the two quadratic forms and using the first identity gives the exact Gauss-Seidel energy decrease identity:Here $2d^{T}Ld=d^T(L+L^T)d$. If $d=0$, then $A{e}_k=0$, so $e_k=0$. Otherwise $d^{T}Dd>0$. Thus $\boxed{\Vert e_k{\Vert }_A>\Vert e_{k+1}{\Vert }_A}$ unless the exact solution has already been reached.

Let $x_{\ast }=A^{-1}b$, $e_k=x_k-x_{\ast }$ and $G=M^{-1}N$. Subtract the exact fixed-point equation to obtain $e_k=G^k{e}_0$. Hence convergence for every initial vector is equivalent to $G^k\to 0$. If this holds, it holds after complexification; applying it to a complex eigenvector gives ${\lambda }^{k}v\to 0$, so $|\lambda |<1$. Conversely, the complex Jordan normal form expresses each power as a finite sum of terms $\left(\genfrac{}{}{0ex}{}{k}{r}\right){\lambda }^{k-r}$ from its blocks. Every term tends to zero when all eigenvalues have modulus below one. ThusThe convergence criterion for a stationary iteration needs complex triangularization or Jordan normal form. The paper's assertion that every real matrix is similar to a real upper-triangular matrix is false when nonreal eigenvalues occur; complex triangularization supplies the intended argument. In the Gauss-Seidel method, with $A=D+L+U$ split into diagonal, strictly lower and strictly upper parts, $\boxed{M=D+L,N=-U}$.

For $n\ge 4$, eliminating the ghost values gives the symmetric pentadiagonal matrixAll other off-diagonal entries vanish. The first two equations areand the last two areFor $n=3$ the matrix is $\left(\begin{array}{ccc}7& -4& 1\\ -4& 6& -4\\ 1& -4& 7\end{array}\right)$, for $n=2$ it is $\left(\begin{array}{cc}7& -4\\ -4& 7\end{array}\right)$, and for $n=1$ it is $(8)$ because both ghost contributions meet the same unknown. The diagonal entries are positive in every case, so the lower triangular matrix $D+L$ is invertible and the Gauss-Seidel iteration is well defined.

An eigenvalue $\lambda$ of $G=-(D+L{)}^{-1}U$ makes $H(\lambda )=\lambda (D+L)+U$ singular. Let ${\ell }_i$ and $r_i$ be the sums of absolute lower and upper entries in row $i$. The row's diagonal-dominance margin isInterior full rows have $d_i=6$, ${\ell }_i=r_i=5$, giving margin $|\lambda |-5$. Boundary rows have larger margins; for instance the first gives $7|\lambda |-5$, and the second gives at least $2|\lambda |-5$. Therefore $|\lambda |>5$ would make $H(\lambda )$ a strictly diagonally dominant matrix, hence nonsingular, a contradiction.

At $|\lambda |=5$ it is weakly diagonally dominant with strict boundary rows. Its nonzero nearest-neighbour entries make it an irreducibly diagonally dominant matrix. The nonsingularity theorem for that class again contradicts singularity. For the small matrices above all rows are already strict at this value. HenceThere is also a stronger conclusion. For $T=\mathrm{tridiag}(-1,2,-1)$,including the $n=1$ overlap. Since $T$ is positive definite, $A$ is a symmetric positive-definite matrix. The Gauss-Seidel convergence from diagonal dominance and positive definiteness theorem consequently gives the stronger bound $\rho (G)<1$.

The Gauss-Seidel method solves each component using the already updated components from the same iteration. In the splitting $A=D-L-U$, those updated diagonal and lower-triangular terms stay on the left, while upper-triangular terms use the previous iterate. ThusHere $D-L$ is lower triangular with nonzero diagonal, so each new iterate is obtained by forward substitution, starting from the specified zero vector.

Update the unknowns in row order, immediately using newly computed components. The component equations areThus $G=\left(\begin{array}{ccc}0& 1& -1\\ 0& -1/2& 0\\ 0& -3/4& 1/2\end{array}\right)$ and $c=(-1,3/2,1/4{)}^T$. Using only old values in all three row equations would instead describe the Jacobi method, not the Gauss-Seidel method.

For $A={\mu }^{-1}I+U-U^T$, with $U$ the four-dimensional unit superdiagonal, the Jacobi method has eigenvalues $2i\mu \cos (j\pi /5)$, $j=1,\dots ,4$. The Gauss-Seidel method has two zero eigenvalues and $-{\mu }^2(3\pm \sqrt{5})/2$. Therefore both methods converge for every starting vector exactly when $0<|\mu |<(\sqrt{5}-1)/2$. This equality of convergence ranges follows from their spectra, not from positive-definiteness of the nonsymmetric coefficient matrix.

For a splitting $A=D+L+U$, successive over-relaxation updatesIt reduces to the Gauss-Seidel method at $\omega =1$. If $A$ is symmetric positive definite, the Householder-John theorem applies to $M=D/\omega +L$, $N=(1/\omega -1)D-U$, since $M^T+N=(2/\omega -1)D$ is positive definite for $0<\omega <2$.