Intuition The big picture (WHY iterate at all?)To solve A x = b A\mathbf{x}=\mathbf{b} A x = b 2 3 n 3 \frac{2}{3}n^3 3 2 n 3 large sparse systems (e.g. n = 10 6 n=10^6 n = 1 0 6 Iterative methods start from a guess x ( 0 ) \mathbf{x}^{(0)} x ( 0 ) repeatedly refine it: each step is cheap (just a matrix–vector product), and for nice matrices the guess marches toward the true answer. We trade exactness in one shot for cheap approximate steps that converge .
Intuition WHAT we're really doingEvery classical iteration is a matrix splitting . Write A = M − N A = M - N A = M − N M M M easy to invert . Then
A x = b ⟺ ( M − N ) x = b ⟺ M x = N x + b . A\mathbf{x}=\mathbf{b} \iff (M-N)\mathbf{x}=\mathbf{b} \iff M\mathbf{x} = N\mathbf{x} + \mathbf{b}. A x = b ⟺ ( M − N ) x = b ⟺ M x = N x + b . M x ( k + 1 ) = N x ( k ) + b . M\mathbf{x}^{(k+1)} = N\mathbf{x}^{(k)} + \mathbf{b}. M x ( k + 1 ) = N x ( k ) + b .
Definition Iteration matrixSolving for the new iterate: x ( k + 1 ) = M − 1 N ⏟ = : T x ( k ) + M − 1 b \mathbf{x}^{(k+1)} = \underbrace{M^{-1}N}_{=:T}\,\mathbf{x}^{(k)} + M^{-1}\mathbf{b} x ( k + 1 ) = =: T M − 1 N x ( k ) + M − 1 b T = = = M − 1 N = = T = ==M^{-1}N== T === M − 1 N == iteration matrix . It alone governs convergence.
Derivation of the error recurrence (from scratch). Let x ∗ \mathbf{x}^* x ∗ M x ∗ = N x ∗ + b M\mathbf{x}^* = N\mathbf{x}^* + \mathbf{b} M x ∗ = N x ∗ + b M x ( k + 1 ) = N x ( k ) + b M\mathbf{x}^{(k+1)} = N\mathbf{x}^{(k)} + \mathbf{b} M x ( k + 1 ) = N x ( k ) + b M ( x ( k + 1 ) − x ∗ ) = N ( x ( k ) − x ∗ ) . M(\mathbf{x}^{(k+1)} - \mathbf{x}^*) = N(\mathbf{x}^{(k)} - \mathbf{x}^*). M ( x ( k + 1 ) − x ∗ ) = N ( x ( k ) − x ∗ ) . e ( k ) = x ( k ) − x ∗ \mathbf{e}^{(k)} = \mathbf{x}^{(k)} - \mathbf{x}^* e ( k ) = x ( k ) − x ∗ e ( k + 1 ) = M − 1 N e ( k ) = T e ( k ) ⇒ e ( k ) = T k e ( 0 ) . \boxed{\;\mathbf{e}^{(k+1)} = M^{-1}N\,\mathbf{e}^{(k)} = T\,\mathbf{e}^{(k)} \;\Rightarrow\; \mathbf{e}^{(k)} = T^{k}\mathbf{e}^{(0)}.\;} e ( k + 1 ) = M − 1 N e ( k ) = T e ( k ) ⇒ e ( k ) = T k e ( 0 ) . T T T Why this matters: convergence (e ( k ) → 0 \mathbf{e}^{(k)}\to 0 e ( k ) → 0 any starting guess happens iff T k → 0 T^k\to 0 T k → 0
Choose M = D M = D M = D N = L + U N = L+U N = L + U M x ( k + 1 ) = N x ( k ) + b M\mathbf{x}^{(k+1)} = N\mathbf{x}^{(k)}+\mathbf{b} M x ( k + 1 ) = N x ( k ) + b x i ( k + 1 ) = 1 a i i ( b i − ∑ j ≠ i a i j x j ( k ) ) . x_i^{(k+1)} = \frac{1}{a_{ii}}\Big(b_i - \sum_{j\neq i} a_{ij}\,x_j^{(k)}\Big). x i ( k + 1 ) = a ii 1 ( b i − ∑ j = i a ij x j ( k ) ) .
Intuition HOW Jacobi thinks"Solve equation i i i x i x_i x i assuming all the other variables are still their old values. " Every component uses the whole old vector — so you can compute all x i ( k + 1 ) x_i^{(k+1)} x i ( k + 1 ) T J = = = D − 1 ( L + U ) = = T_J = ==D^{-1}(L+U)== T J === D − 1 ( L + U ) ==
Choose M = D − L M = D - L M = D − L N = U N = U N = U x i ( k + 1 ) = 1 a i i ( b i − ∑ j < i a i j x j ( k + 1 ) − ∑ j > i a i j x j ( k ) ) . x_i^{(k+1)} = \frac{1}{a_{ii}}\Big(b_i - \sum_{j<i} a_{ij}\,x_j^{(k+1)} - \sum_{j>i} a_{ij}\,x_j^{(k)}\Big). x i ( k + 1 ) = a ii 1 ( b i − ∑ j < i a ij x j ( k + 1 ) − ∑ j > i a ij x j ( k ) ) .
Intuition HOW Gauss-Seidel improves on JacobiThe only change: use the freshly-computed components immediately . By the time you reach x i x_i x i x 1 , … , x i − 1 x_1,\dots,x_{i-1} x 1 , … , x i − 1 this sweep — so use them! This usually halves the iterations but kills parallelism (it's sequential). Iteration matrix T G S = = = ( D − L ) − 1 U = = T_{GS} = ==(D-L)^{-1}U== T GS === ( D − L ) − 1 U ==
Definition Spectral radiusρ ( T ) = max i ∣ λ i ∣ \rho(T) = \max_i |\lambda_i| ρ ( T ) = max i ∣ λ i ∣ absolute eigenvalue of T T T
WHY ρ < 1 \rho<1 ρ < 1 Diagonalize T = P Λ P − 1 T = P\Lambda P^{-1} T = P Λ P − 1 λ i \lambda_i λ i T k = P Λ k P − 1 T^k = P\Lambda^k P^{-1} T k = P Λ k P − 1 Λ k = diag ( λ i k ) \Lambda^k = \operatorname{diag}(\lambda_i^k) Λ k = diag ( λ i k ) k → ∞ k\to\infty k → ∞ λ i k → 0 \lambda_i^k \to 0 λ i k → 0 all i i i ∣ λ i ∣ < 1 |\lambda_i|<1 ∣ λ i ∣ < 1 ρ ( T ) < 1 \rho(T)<1 ρ ( T ) < 1 e ( k ) = T k e ( 0 ) \mathbf{e}^{(k)}=T^k\mathbf{e}^{(0)} e ( k ) = T k e ( 0 )
sufficient easy check: diagonal dominance
Computing eigenvalues is hard. A practical sufficient (not necessary) condition: if A A A strictly diagonally dominant (SDD),
∣ a i i ∣ > ∑ j ≠ i ∣ a i j ∣ for every row i , |a_{ii}| > \sum_{j\neq i}|a_{ij}| \quad \text{for every row } i, ∣ a ii ∣ > ∑ j = i ∣ a ij ∣ for every row i , both Jacobi and Gauss-Seidel converge for any start. WHY it feels right: each x i x_i x i
Worked example Worked example 1 — solve a2 × 2 2\times2 2 × 2 { 4 x + y = 6 x + 3 y = 7 , exact: x = 11 11 = 1 , y = 2. \begin{cases}4x + y = 6\\ x + 3y = 7\end{cases},\qquad \text{exact: } x=\tfrac{11}{11}=1,\ y=2. { 4 x + y = 6 x + 3 y = 7 , exact: x = 11 11 = 1 , y = 2. a 11 = 4 > 1 a_{11}=4>1 a 11 = 4 > 1 a 22 = 3 > 1 a_{22}=3>1 a 22 = 3 > 1 x ( 0 ) = ( 0 , 0 ) \mathbf{x}^{(0)}=(0,0) x ( 0 ) = ( 0 , 0 )
Jacobi x = 6 − y 4 , y = 7 − x 3 x = \frac{6-y}{4},\ y=\frac{7-x}{3} x = 4 6 − y , y = 3 7 − x
k = 1 k{=}1 k = 1 x = 6 − 0 4 = 1.5 , y = 7 − 0 3 = 2.333 x=\frac{6-0}{4}=1.5,\ y=\frac{7-0}{3}=2.333 x = 4 6 − 0 = 1.5 , y = 3 7 − 0 = 2.333 Why? both use old ( 0 , 0 ) (0,0) ( 0 , 0 ) k = 2 k{=}2 k = 2 x = 6 − 2.333 4 = 0.917 , y = 7 − 1.5 3 = 1.833 x=\frac{6-2.333}{4}=0.917,\ y=\frac{7-1.5}{3}=1.833 x = 4 6 − 2.333 = 0.917 , y = 3 7 − 1.5 = 1.833 Why? use k = 1 k{=}1 k = 1 k = 3 k{=}3 k = 3 x = 6 − 1.833 4 = 1.042 , y = 7 − 0.917 3 = 2.028 x=\frac{6-1.833}{4}=1.042,\ y=\frac{7-0.917}{3}=2.028 x = 4 6 − 1.833 = 1.042 , y = 3 7 − 0.917 = 2.028 ( 1 , 2 ) (1,2) ( 1 , 2 )
Gauss-Seidel (use new x x x y y y
k = 1 k{=}1 k = 1 x = 6 − 0 4 = 1.5 x=\frac{6-0}{4}=1.5 x = 4 6 − 0 = 1.5 y = 7 − 1.5 3 = 1.833 y=\frac{7-1.5}{3}=1.833 y = 3 7 − 1.5 = 1.833 Why this step? y y y just-computed x = 1.5 x=1.5 x = 1.5 0 0 0 k = 2 k{=}2 k = 2 x = 6 − 1.833 4 = 1.042 x=\frac{6-1.833}{4}=1.042 x = 4 6 − 1.833 = 1.042 y = 7 − 1.042 3 = 1.986 y=\frac{7-1.042}{3}=1.986 y = 3 7 − 1.042 = 1.986 k = 3 k{=}3 k = 3 x = 6 − 1.986 4 = 1.004 x=\frac{6-1.986}{4}=1.004 x = 4 6 − 1.986 = 1.004 y = 7 − 1.004 3 = 1.999 y=\frac{7-1.004}{3}=1.999 y = 3 7 − 1.004 = 1.999 Faster — fewer steps to reach ( 1 , 2 ) (1,2) ( 1 , 2 )
Worked example Worked example 2 — convergence test by eigenvalueFor A = ( 2 − 1 − 1 2 ) A=\begin{pmatrix}2&-1\\-1&2\end{pmatrix} A = ( 2 − 1 − 1 2 ) D = 2 I D=2I D = 2 I L + U = ( 0 1 1 0 ) L+U=\begin{pmatrix}0&1\\1&0\end{pmatrix} L + U = ( 0 1 1 0 ) T J = D − 1 ( L + U ) = 1 2 ( 0 1 1 0 ) T_J = D^{-1}(L+U)=\tfrac12\begin{pmatrix}0&1\\1&0\end{pmatrix} T J = D − 1 ( L + U ) = 2 1 ( 0 1 1 0 ) ( 0 1 1 0 ) \begin{pmatrix}0&1\\1&0\end{pmatrix} ( 0 1 1 0 ) ± 1 \pm1 ± 1 λ ( T J ) = ± 1 2 \lambda(T_J)=\pm\tfrac12 λ ( T J ) = ± 2 1 ρ ( T J ) = 1 2 < 1 \rho(T_J)=\tfrac12 < 1 ρ ( T J ) = 2 1 < 1 Why this matters: even without iterating we predicted the speed.
Worked example Worked example 3 — Forecast-then-Verify (does it diverge?)A = ( 1 2 3 1 ) A=\begin{pmatrix}1&2\\3&1\end{pmatrix} A = ( 1 3 2 1 ) ∣ 1 ∣ < ∣ 2 ∣ |1|<|2| ∣1∣ < ∣2∣ not diagonally dominant. Forecast: likely diverges.
T J = D − 1 ( L + U ) = ( 0 − 2 − 3 0 ) T_J=D^{-1}(L+U)=\begin{pmatrix}0&-2\\-3&0\end{pmatrix} T J = D − 1 ( L + U ) = ( 0 − 3 − 2 0 ) a 11 = a 22 = 1 a_{11}=a_{22}=1 a 11 = a 22 = 1 λ 2 − 6 = 0 ⇒ λ = ± 6 ≈ ± 2.449 \lambda^2 - 6 = 0 \Rightarrow \lambda=\pm\sqrt6\approx\pm2.449 λ 2 − 6 = 0 ⇒ λ = ± 6 ≈ ± 2.449 ρ = 2.449 > 1 \rho=2.449>1 ρ = 2.449 > 1 diverges . ✓ forecast confirmed. Fix: swap rows to ( 3 1 1 2 ) \begin{pmatrix}3&1\\1&2\end{pmatrix} ( 3 1 1 2 )
Common mistake "Diagonal dominance fails ⇒ it must diverge."Why it feels right: SDD is the rule we're taught for convergence. The flaw: SDD is sufficient, not necessary . The truth is ρ ( T ) < 1 \rho(T)<1 ρ ( T ) < 1 ρ ( T ) < 1 \rho(T)<1 ρ ( T ) < 1 Fix: if SDD fails, don't conclude — actually examine ρ ( T ) \rho(T) ρ ( T )
Common mistake "Gauss-Seidel always beats Jacobi."Why it feels right: GS uses fresher info, and for SDD / SPD matrices it indeed converges faster (ρ ( T G S ) = ρ ( T J ) 2 \rho(T_{GS})=\rho(T_J)^2 ρ ( T GS ) = ρ ( T J ) 2 The flaw: there exist matrices where Jacobi converges but Gauss-Seidel diverges , and vice-versa. Fix: "usually faster," not "always converges."
Common mistake Updating Gauss-Seidel like Jacobi (storing old vector).Why it feels right: the formulas look almost identical. The flaw: if you keep using the old vector for the lower-index terms, you've just coded Jacobi. Fix: in GS, overwrite x i x_i x i in place so later components in the same sweep see it.
Common mistake Forgettinga i i ≠ 0 a_{ii}\neq0 a ii = 0 a i i a_{ii} a ii 0 0 0 permute rows first (a row swap) to put nonzeros on the diagonal.
Intuition When do we stop?Iterate until ∥ x ( k + 1 ) − x ( k ) ∥ < ε \|\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}\| < \varepsilon ∥ x ( k + 1 ) − x ( k ) ∥ < ε ∥ b − A x ( k ) ∥ < ε \|\mathbf{b}-A\mathbf{x}^{(k)}\| < \varepsilon ∥ b − A x ( k ) ∥ < ε O ( nnz ) O(\text{nnz}) O ( nnz ) A A A O ( n ) O(n) O ( n ) O ( n 3 ) O(n^3) O ( n 3 ) That's the 80/20 win: cheap repeated steps + good ρ \rho ρ
What is the iteration matrix T T T A = M − N A=M-N A = M − N T = M − 1 N T = M^{-1}N T = M − 1 N ; the iteration is
x ( k + 1 ) = T x ( k ) + M − 1 b \mathbf{x}^{(k+1)}=T\mathbf{x}^{(k)}+M^{-1}\mathbf{b} x ( k + 1 ) = T x ( k ) + M − 1 b .
Error recurrence of a stationary iteration? e ( k + 1 ) = T e ( k ) \mathbf{e}^{(k+1)}=T\mathbf{e}^{(k)} e ( k + 1 ) = T e ( k ) , so
e ( k ) = T k e ( 0 ) \mathbf{e}^{(k)}=T^k\mathbf{e}^{(0)} e ( k ) = T k e ( 0 ) .
Exact necessary-and-sufficient condition for convergence from any start? Spectral radius
ρ ( T ) < 1 \rho(T)<1 ρ ( T ) < 1 .
Splitting M , N M,N M , N M = D , N = L + U M=D,\ N=L+U M = D , N = L + U , giving
T J = D − 1 ( L + U ) T_J=D^{-1}(L+U) T J = D − 1 ( L + U ) .
Splitting M , N M,N M , N M = D − L , N = U M=D-L,\ N=U M = D − L , N = U , giving
T G S = ( D − L ) − 1 U T_{GS}=(D-L)^{-1}U T GS = ( D − L ) − 1 U .
Component formula for Jacobi? x i ( k + 1 ) = 1 a i i ( b i − ∑ j ≠ i a i j x j ( k ) ) x_i^{(k+1)}=\frac{1}{a_{ii}}\big(b_i-\sum_{j\neq i}a_{ij}x_j^{(k)}\big) x i ( k + 1 ) = a ii 1 ( b i − ∑ j = i a ij x j ( k ) ) .
Component formula for Gauss-Seidel? x i ( k + 1 ) = 1 a i i ( b i − ∑ j < i a i j x j ( k + 1 ) − ∑ j > i a i j x j ( k ) ) x_i^{(k+1)}=\frac{1}{a_{ii}}\big(b_i-\sum_{j<i}a_{ij}x_j^{(k+1)}-\sum_{j>i}a_{ij}x_j^{(k)}\big) x i ( k + 1 ) = a ii 1 ( b i − ∑ j < i a ij x j ( k + 1 ) − ∑ j > i a ij x j ( k ) ) .
Key practical difference Jacobi vs Gauss-Seidel? GS uses freshly updated components within the same sweep (in-place); Jacobi uses only the old vector (parallelisable).
A simple sufficient condition guaranteeing both converge? Strict diagonal dominance:
∣ a i i ∣ > ∑ j ≠ i ∣ a i j ∣ |a_{ii}|>\sum_{j\neq i}|a_{ij}| ∣ a ii ∣ > ∑ j = i ∣ a ij ∣ for all
i i i .
Is diagonal dominance necessary for convergence? No — it is sufficient only; the true criterion is
ρ ( T ) < 1 \rho(T)<1 ρ ( T ) < 1 .
Asymptotic error reduction factor per step? Approximately
ρ ( T ) \rho(T) ρ ( T ) .
Why iterative over direct methods for large sparse A A A Each sweep costs
O ( nnz ) ≈ O ( n ) O(\text{nnz})\approx O(n) O ( nnz ) ≈ O ( n ) vs
O ( n 3 ) O(n^3) O ( n 3 ) for elimination, and memory stays sparse.
Recall Feynman: explain it to a 12-year-oldImagine guessing everyone's height in a class, then fixing each guess using a rule that ties heights together. Jacobi: everybody updates their guess at the same time using last round's numbers. Gauss-Seidel: you go one by one, and each person uses the just-corrected numbers of the people before them — so corrections spread faster. If the rule "pulls" guesses gently toward the truth (numbers don't blow up — that's ρ < 1 \rho<1 ρ < 1 ρ > 1 \rho>1 ρ > 1
"Jacobi waits, Seidel updates." And for convergence: "Spectral radius Sub-one Settles" (ρ < 1 \rho<1 ρ < 1 D iagonally D ominant ⇒ D efinitely converges.
Gaussian Elimination — the direct O ( n 3 ) O(n^3) O ( n 3 ) Eigenvalues and Eigenvectors — convergence rests entirely on ρ ( T ) \rho(T) ρ ( T ) Spectral Radius — the master quantity here.Fixed-Point Iteration — these are the matrix version of x = g ( x ) x=g(x) x = g ( x ) SOR — Successive Over-Relaxation — accelerates Gauss-Seidel with a relaxation factor ω \omega ω Sparse Matrices — why per-sweep cost is cheap.Finite Difference Methods for PDEs — main source of the huge sparse systems.
direct too costly for large sparse
converges iff T^k to zero
uses all old values, parallel
uses fresh components immediately
Matrix splitting A = M - N
Iteration matrix T = inv M times N
Intuition Hinglish mein samjhoDekho, jab humein bahut bada system A x = b A\mathbf{x}=\mathbf{b} A x = b O ( n 3 ) O(n^3) O ( n 3 ) iterative methods : ek guess se start karo, aur baar-baar thoda-thoda sudhaarte jao jab tak answer paas na aa jaaye. Har step sasta hota hai, aur agar matrix "achhi" hai to guess slowly sahi answer ki taraf march karta hai.
Core trick hai splitting : A = M − N A = M - N A = M − N M M M x ( k + 1 ) = T x ( k ) + M − 1 b \mathbf{x}^{(k+1)} = T\mathbf{x}^{(k)} + M^{-1}\mathbf{b} x ( k + 1 ) = T x ( k ) + M − 1 b T = M − 1 N T=M^{-1}N T = M − 1 N iteration matrix hai. Jacobi mein M = D M=D M = D x i x_i x i Gauss-Seidel mein M = D − L M=D-L M = D − L x i x_i x i
Sabse important baat — convergence . Error har step pe T T T e ( k ) = T k e ( 0 ) \mathbf{e}^{(k)}=T^k\mathbf{e}^{(0)} e ( k ) = T k e ( 0 ) spectral radius ρ ( T ) < 1 \rho(T)<1 ρ ( T ) < 1 ∣ λ ∣ |\lambda| ∣ λ ∣ strictly diagonally dominant hai (∣ a i i ∣ > ∑ j ≠ i ∣ a i j ∣ |a_{ii}|>\sum_{j\neq i}|a_{ij}| ∣ a ii ∣ > ∑ j = i ∣ a ij ∣ sufficient condition hai, necessary nahi. Asli rule hamesha ρ ( T ) < 1 \rho(T)<1 ρ ( T ) < 1
