Systems of Equations

2.1 Gaussian Elimination

1) Linear System

Linear System: A collection of linear equations written in matrix form as

$$Ax=b,$$

where:

• $A\in\mathbb{R}^{n\times n}$ is the coefficient matrix.
• $x\in\mathbb{R}^n$ is the unknown vector.
• $b\in\mathbb{R}^n$ is the right-hand-side vector.

2) Augmented Matrix

Augmented Matrix: A matrix containing both the coefficient matrix and right-hand-side vector.

$$[A\mid b].$$

For example,

$$\begin{bmatrix} a_{11} & a_{12} & \mid & b_1\\ a_{21} & a_{22} & \mid & b_2 \end{bmatrix}.$$

3) Elementary Row Operations

The following operations preserve the solution set of a linear system:

1. Exchange two rows.
2. Add a multiple of one row to another row.
3. Multiply a row by a nonzero scalar.

4) Gaussian Elimination

Gaussian Elimination: A direct method that transforms a linear system into an upper-triangular system.

To eliminate $a_{ij}$ below the pivot $a_{jj}$, define the multiplier

$$m_{ij} = \frac{a_{ij}}{a_{jj}}.$$

Then perform

$$\operatorname{row}_i \leftarrow \operatorname{row}_i - m_{ij}\operatorname{row}_j.$$

After elimination, the system has the form

$$Ux=c,$$

where $U$ is upper triangular.

5) Pivot

Pivot: The diagonal entry used to eliminate entries below it.

At elimination step $j$, the pivot is

$$a_{jj}.$$

Naive Gaussian elimination fails if a zero pivot occurs and may be numerically unstable if a pivot is very small.

6) Upper-Triangular Matrix

Upper-Triangular Matrix: A square matrix $U$ satisfying

$$u_{ij}=0 \qquad\text{for }i>j.$$

All entries below the main diagonal are zero.

7) Back Substitution

Back Substitution: The process of solving an upper-triangular system from the last equation upward.

For

$$Ux=c,$$

the solution components are

$$x_i = \frac{ c_i-\displaystyle\sum_{j=i+1}^{n}u_{ij}x_j }{ u_{ii} }, \qquad i=n,n-1,\ldots,1.$$

Back substitution requires nonzero diagonal entries.

8) Operation Count

The elimination stage requires approximately

$$\frac{2}{3}n^3$$

arithmetic operations.

Back substitution requires approximately

$$n^2$$

operations.

Therefore, the total cost is

$$\frac{2}{3}n^3+n^2 \approx \frac{2}{3}n^3.$$

Gaussian elimination has time complexity

$$O(n^3).$$

2.2 The LU Factorization

1) Lower-Triangular Matrix

Lower-Triangular Matrix: A square matrix $L$ satisfying

$$l_{ij}=0 \qquad\text{for }i<j.$$

All entries above the main diagonal are zero.

2) Unit Lower-Triangular Matrix

Unit Lower-Triangular Matrix: A lower-triangular matrix whose diagonal entries are all $1$.

$$l_{ii}=1.$$

3) LU Factorization

LU Factorization: A decomposition of a square matrix into a lower-triangular matrix and an upper-triangular matrix.

$$A=LU.$$

The matrix $U$ contains the result of Gaussian elimination.

The elimination multipliers are stored below the diagonal of $L$.

4) Solving with LU Factorization

Given

$$Ax=b$$

and

$$A=LU,$$

the system becomes

$$LUx=b.$$

Introduce an intermediate vector $y$:

$$Ly=b,$$

followed by

$$Ux=y.$$

The solution is obtained in two stages:

1. Solve $Ly=b$ using forward substitution.
2. Solve $Ux=y$ using back substitution.

5) Forward Substitution

Forward Substitution: The process of solving a lower-triangular system from the first equation downward.

For

$$Ly=b,$$

the components are

$$y_i = \frac{ b_i-\displaystyle\sum_{j=1}^{i-1}l_{ij}y_j }{ l_{ii} }, \qquad i=1,2,\ldots,n.$$

If $L$ is unit lower triangular, then $l_{ii}=1$.

6) Complexity of LU Factorization

Computing the factorization costs approximately

$$\frac{2}{3}n^3$$

operations.

Solving each triangular system costs approximately

$$n^2$$

operations.

For $k$ different right-hand sides,

$$Ax=b_1,\ldots,Ax=b_k,$$

the total cost is approximately

$$\frac{2}{3}n^3+2kn^2.$$

LU factorization is efficient when several systems share the same coefficient matrix.

7) Existence of LU Factorization

An LU factorization without row exchanges may fail if Gaussian elimination encounters a zero pivot.

Therefore, not every nonsingular matrix has a factorization of the form

$$A=LU$$

without row permutations.

2.3 Sources of Error

1) Vector Norm

Vector Norm: A function that measures the size of a vector and satisfies:

$$\|x\|\geq0,$$ $$\|x\|=0 \iff x=0,$$ $$\|\alpha x\| = |\alpha|\|x\|,$$

and

$$\|x+y\| \leq \|x\|+\|y\|.$$

2) Infinity Norm of a Vector

Infinity Norm: The largest absolute component of a vector.

For

$$x=(x_1,\ldots,x_n)^T,$$

the infinity norm is

$$\|x\|_\infty = \max_{1\leq i\leq n}|x_i|.$$

3) Matrix Infinity Norm

Matrix Infinity Norm: The largest absolute row sum of a matrix.

$$\|A\|_\infty = \max_{1\leq i\leq n} \sum_{j=1}^{n}|a_{ij}|.$$

4) Residual

Residual: The difference between the required right-hand side and the right-hand side produced by an approximate solution $x_a$.

$$r=b-Ax_a.$$

The residual measures how accurately the computed solution satisfies the original system.

5) Backward Error

Backward Error: The norm of the residual.

$$E_{\mathrm{backward}} = \|b-Ax_a\|_\infty.$$

The relative backward error is

$$\frac{\|b-Ax_a\|_\infty}{\|b\|_\infty}.$$

6) Forward Error

Forward Error: The distance between the exact solution $x$ and the approximate solution $x_a$.

$$E_{\mathrm{forward}} = \|x-x_a\|_\infty.$$

The relative forward error is

$$\frac{\|x-x_a\|_\infty}{\|x\|_\infty}.$$

A small residual does not necessarily imply a small forward error.

7) Error Magnification Factor

Error Magnification Factor: The ratio of relative forward error to relative backward error.

$$\text{Error Magnification Factor} = \frac{ \|x-x_a\|_\infty/\|x\|_\infty }{ \|b-Ax_a\|_\infty/\|b\|_\infty }.$$

A large value means that small input errors can produce large solution errors.

8) Condition Number

Condition Number: The maximum possible error magnification when solving $Ax=b$ over all right-hand sides.

For a nonsingular matrix,

$$\operatorname{cond}(A) = \|A\|\|A^{-1}\|.$$

Using the infinity norm,

$$\operatorname{cond}_\infty(A) = \|A\|_\infty \|A^{-1}\|_\infty.$$

9) Condition Number Properties

For a nonsingular matrix,

$$\operatorname{cond}(A)\geq1.$$

If

$$\operatorname{cond}(A)\approx10^k,$$

up to approximately $k$ decimal digits may be lost through error magnification.

The expected relative forward error can be as large as

$$\operatorname{cond}(A) \times \text{relative backward error}.$$

10) Well-Conditioned System

Well-Conditioned System: A system whose solution changes only slightly when the input data changes slightly.

A well-conditioned matrix has a relatively small condition number.

11) Ill-Conditioned System

Ill-Conditioned System: A system in which small input errors can produce large changes in the solution.

An ill-conditioned matrix has a large condition number.

Ill-conditioning is a property of the mathematical problem, not the solution algorithm.

12) Swamping

Swamping: A loss of significant digits caused when a large floating-point number dominates a much smaller number during arithmetic.

In Gaussian elimination, swamping may occur when a very small pivot produces a very large multiplier.

$$m_{ij} = \frac{a_{ij}}{a_{jj}}.$$

Partial pivoting is used to reduce this problem.

2.4 The $PA=LU$ Factorization

1) Row Exchange

Row Exchange: Interchanging two rows so that a more suitable pivot is placed in the pivot position.

Row exchanges avoid zero pivots and reduce the use of dangerously small pivots.

2) Partial Pivoting

Partial Pivoting: At elimination step $j$, select the entry with the largest absolute value in column $j$ at or below the diagonal.

Choose $p$ satisfying

$$|a_{pj}| = \max_{j\leq i\leq n}|a_{ij}|.$$

Then exchange rows $p$ and $j$.

Partial pivoting keeps elimination multipliers bounded by

$$|m_{ij}|\leq1.$$

3) Permutation Matrix

Permutation Matrix: A matrix obtained by rearranging the rows of the identity matrix.

Multiplying a matrix on the left by $P$ permutes its rows.

$$PA = \text{row-permuted form of }A.$$

Permutation matrices satisfy

$$P^{-1}=P^T.$$

4) $PA=LU$ Factorization

$PA=LU$ Factorization: An LU factorization that records the row exchanges performed during partial pivoting.

$$PA=LU.$$

Here:

• $P$ is a permutation matrix.
• $L$ is unit lower triangular.
• $U$ is upper triangular.

5) Solving with $PA=LU$

Starting from

$$Ax=b,$$

multiply both sides by $P$:

$$PAx=Pb.$$

Using $PA=LU$,

$$LUx=Pb.$$

The solution is found by solving

$$Ly=Pb$$

and then

$$Ux=y.$$

6) Numerical Stability

Gaussian elimination with partial pivoting is numerically stable for most practical systems.

It is the standard general-purpose direct method for solving dense linear systems.

2.5 Iterative Methods

1) Direct Method

Direct Method: A method that obtains the solution after a finite sequence of arithmetic operations in exact arithmetic.

Examples include:

• Gaussian elimination
• LU factorization
• Cholesky factorization

2) Iterative Method

Iterative Method: A method that generates a sequence of approximations

$$x^{(0)},x^{(1)},x^{(2)},\ldots$$

that is expected to converge to the exact solution.

Iterative methods are especially useful for large, sparse systems.

3) Matrix Splitting

Write the coefficient matrix as

$$A=D+L+U,$$

where:

• $D$ contains the diagonal entries.
• $L$ contains the entries below the diagonal.
• $U$ contains the entries above the diagonal.

4) Jacobi Method

Jacobi Method: An iterative method that calculates every component of the new approximation using only values from the previous iteration.

For $Ax=b$,

$$x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{\substack{j=1\\j\neq i}}^n a_{ij}x_j^{(k)} \right).$$

In matrix form,

$$x^{(k+1)} = D^{-1} \left( b-(L+U)x^{(k)} \right).$$

5) Strict Diagonal Dominance

Strictly Diagonally Dominant Matrix: A matrix satisfying

$$|a_{ii}| > \sum_{\substack{j=1\\j\neq i}}^n |a_{ij}|$$

for every row $i$.

If $A$ is strictly diagonally dominant, then:

• $A$ is nonsingular.
• The Jacobi Method converges for every initial guess.
• The Gauss-Seidel Method converges for every initial guess.

Strict diagonal dominance is sufficient but not necessary for convergence.

6) Gauss-Seidel Method

Gauss-Seidel Method: An iterative method that immediately uses newly computed component values.

$$x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j=1}^{i-1} a_{ij}x_j^{(k+1)} - \sum_{j=i+1}^{n} a_{ij}x_j^{(k)} \right).$$

In matrix form,

$$(D+L)x^{(k+1)} = b-Ux^{(k)}.$$

Therefore,

$$x^{(k+1)} = (D+L)^{-1} \left( b-Ux^{(k)} \right).$$

Gauss-Seidel often converges faster than Jacobi when both methods converge.

7) Successive Over-Relaxation

Successive Over-Relaxation: A weighted version of the Gauss-Seidel Method.

First compute the Gauss-Seidel update $\widehat{x}_i^{(k+1)}$, then define

$$x_i^{(k+1)} = (1-\omega)x_i^{(k)} + \omega\widehat{x}_i^{(k+1)}.$$

The parameter $\omega$ is the relaxation factor.

• $\omega=1$: Gauss-Seidel Method
• $0<\omega<1$: under-relaxation
• $1<\omega<2$: over-relaxation

A suitable value of $\omega$ may accelerate convergence.

8) Iteration Matrix

A stationary iterative method can be written as

$$x^{(k+1)} = Tx^{(k)}+c.$$

The matrix $T$ is the iteration matrix.

The error satisfies

$$e^{(k+1)} = Te^{(k)}.$$

Therefore,

$$e^{(k)} = T^k e^{(0)}.$$

9) Spectral Radius

Spectral Radius: The largest absolute value among the eigenvalues of a matrix.

$$\rho(T) = \max_i|\lambda_i(T)|.$$

The iteration converges for every initial guess if and only if

$$\rho(T)<1.$$

10) Sparse Matrix

Sparse Matrix: A matrix in which most entries are zero.

Sparse matrices should be stored using formats that record only nonzero entries.

For a sparse matrix with $N$ nonzero entries, a matrix-vector multiplication can require approximately

$$O(N)$$

operations instead of $O(n^2)$.

11) Fill-In

Fill-In: The creation of new nonzero entries during Gaussian elimination at positions that were originally zero.

Fill-in can greatly increase memory use and computation time for sparse direct methods.

Iterative methods avoid much of this fill-in.

2.6 Symmetric Positive-Definite Matrices

1) Symmetric Matrix

Symmetric Matrix: A square matrix satisfying

$$A^T=A.$$

Equivalently,

$$a_{ij}=a_{ji}.$$

2) Positive-Definite Matrix

Positive-Definite Matrix: A symmetric matrix satisfying

$$x^TAx>0$$

for every nonzero vector $x$.

3) Symmetric Positive-Definite Matrix

Symmetric Positive-Definite Matrix: A matrix that is both symmetric and positive definite.

Such a matrix is nonsingular and has positive eigenvalues.

4) Principal Submatrix

Principal Submatrix: A square submatrix formed by selecting the same index set of rows and columns.

A leading principal submatrix uses the first $k$ rows and first $k$ columns.

5) Positive-Definite Test

A symmetric matrix is positive definite if and only if all its leading principal determinants are positive.

$$\det(A_k)>0, \qquad k=1,\ldots,n,$$

where $A_k$ is the leading $k\times k$ principal submatrix.

6) Cholesky Factorization

Cholesky Factorization: A factorization of a symmetric positive-definite matrix as

$$A=R^TR,$$

where $R$ is upper triangular with positive diagonal entries.

Equivalently,

$$A=LL^T,$$

where $L$ is lower triangular.

7) Cholesky Factorization Theorem

Every symmetric positive-definite matrix has a unique Cholesky factorization with positive diagonal entries.

The entries of $R$ can be calculated by

$$r_{ii} = \sqrt{ a_{ii} - \sum_{k=1}^{i-1}r_{ki}^2 },$$

and, for $j>i$,

$$r_{ij} = \frac{ a_{ij} - \displaystyle\sum_{k=1}^{i-1}r_{ki}r_{kj} }{ r_{ii} }.$$

8) Solving with Cholesky Factorization

For

$$Ax=b$$

and

$$A=R^TR,$$

solve

$$R^Ty=b$$

using forward substitution, followed by

$$Rx=y$$

using back substitution.

Cholesky factorization requires approximately half the operations and storage of a general LU factorization.

9) $A$-Inner Product

For a symmetric positive-definite matrix $A$, define

$$\langle v,w\rangle_A = v^TAw.$$

This is an inner product induced by $A$.

10) $A$-Conjugate Vectors

$A$-Conjugate Vectors: Two vectors $v$ and $w$ satisfying

$$v^TAw=0.$$

They are also called $A$-orthogonal vectors.

11) Conjugate Gradient Method

Conjugate Gradient Method: A method for solving

$$Ax=b$$

when $A$ is symmetric positive definite.

Start with an initial guess $x_0$ and define

$$r_0=b-Ax_0,$$ $$d_0=r_0.$$

For $k=0,1,2,\ldots$, calculate

$$\alpha_k = \frac{r_k^Tr_k} {d_k^TAd_k},$$ $$x_{k+1} = x_k+\alpha_kd_k,$$ $$r_{k+1} = r_k-\alpha_kAd_k,$$ $$\beta_k = \frac{r_{k+1}^Tr_{k+1}} {r_k^Tr_k},$$

and

$$d_{k+1} = r_{k+1}+\beta_kd_k.$$

12) Conjugate Gradient Properties

In exact arithmetic, the Conjugate Gradient Method finds the exact solution in at most $n$ steps.

The search directions satisfy

$$d_i^TAd_j=0 \qquad\text{for }i\neq j.$$

The residuals satisfy

$$r_i^Tr_j=0 \qquad\text{for }i\neq j.$$

In practice, it is used as an iterative method and stopped when the residual becomes sufficiently small.

It is particularly effective for large, sparse, symmetric positive-definite systems.

13) Preconditioning

Preconditioning: Transforming a linear system into an equivalent system with better convergence properties.

Choose a matrix $M$ that approximates $A$ but is easier to solve with.

The preconditioned system is

$$M^{-1}Ax=M^{-1}b.$$

A useful preconditioner reduces the condition number while remaining inexpensive to apply.

14) Jacobi Preconditioner

Jacobi Preconditioner: A preconditioner using the diagonal of $A$.

$$M=D.$$

It is simple and inexpensive but may provide limited improvement.

15) Preconditioned Conjugate Gradient Method

For the residual

$$r_k=b-Ax_k,$$

solve

$$Mz_k=r_k.$$

Then use $z_k$ in place of $r_k$ when constructing the search direction.

Preconditioning can greatly reduce the number of iterations required for convergence.

2.7 Nonlinear Systems of Equations

1) Nonlinear System

A system of $n$ nonlinear equations in $n$ variables can be written as

$$F(x)=0,$$

where

$$F(x) = \begin{bmatrix} f_1(x)\\ f_2(x)\\ \vdots\\ f_n(x) \end{bmatrix}.$$

A solution is a vector $r$ satisfying

$$F(r)=0.$$

2) Jacobian Matrix

Jacobian Matrix: The matrix of first-order partial derivatives of a vector-valued function.

$$DF(x) = J_F(x) = \begin{bmatrix} \dfrac{\partial f_1}{\partial x_1} & \cdots & \dfrac{\partial f_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \dfrac{\partial f_n}{\partial x_1} & \cdots & \dfrac{\partial f_n}{\partial x_n} \end{bmatrix}.$$

The Jacobian is the multivariable analogue of the derivative.

3) Multivariable Taylor Expansion

Near $x_k$,

$$F(x) \approx F(x_k) + DF(x_k)(x-x_k).$$

Setting the approximation equal to zero produces the multivariable Newton step.

4) Multivariate Newton’s Method

Multivariate Newton’s Method: An iterative method for solving nonlinear systems.

The formal iteration is

$$x_{k+1} = x_k - DF(x_k)^{-1}F(x_k).$$

The inverse should not be computed explicitly.

Instead, solve

$$DF(x_k)s_k = -F(x_k)$$

and update

$$x_{k+1} = x_k+s_k.$$

5) Convergence of Multivariate Newton’s Method

If:

• The initial guess is sufficiently close to a solution $r$.
• $DF(r)$ is nonsingular.
• The necessary derivatives are continuous.

Then Multivariate Newton’s Method is locally quadratically convergent.

The method may fail if the Jacobian is singular or the initial guess is unsuitable.

6) Broyden’s Method

Broyden’s Method: A quasi-Newton method that approximates the Jacobian instead of recomputing it at every iteration.

It is useful when the Jacobian is unavailable or expensive to calculate.

7) Secant Condition

Let

$$\delta_{k+1} = x_{k+1}-x_k$$

and

$$\gamma_{k+1} = F(x_{k+1})-F(x_k).$$

The updated Jacobian approximation $A_{k+1}$ is required to satisfy

$$A_{k+1}\delta_{k+1} = \gamma_{k+1}.$$

This is the multidimensional analogue of the secant condition.

8) Broyden’s Method I

Let $A_k$ approximate the Jacobian.

The iteration is

$$x_{k+1} = x_k-A_k^{-1}F(x_k).$$

The Jacobian approximation is updated by

$$A_{k+1} = A_k + \frac{ \left( \gamma_{k+1}-A_k\delta_{k+1} \right) \delta_{k+1}^T }{ \delta_{k+1}^T\delta_{k+1} }.$$

In implementation, solve

$$A_ks_k=-F(x_k)$$

and set

$$x_{k+1}=x_k+s_k.$$

9) Broyden’s Method II

Let $B_k$ approximate the inverse Jacobian.

The iteration is

$$x_{k+1} = x_k-B_kF(x_k).$$

The inverse approximation is updated by

$$B_{k+1} = B_k + \frac{ \left( \delta_{k+1}-B_k\gamma_{k+1} \right) \delta_{k+1}^TB_k }{ \delta_{k+1}^TB_k\gamma_{k+1} }.$$

Broyden’s Method II avoids solving a new linear system at every iteration.

10) Broyden Convergence

At a simple root, Broyden’s Method generally converges superlinearly.

It is usually slower than Newton’s quadratic convergence but does not require repeated Jacobian evaluations.

Neither Newton’s Method nor Broyden’s Method guarantees convergence from an arbitrary initial guess.

Comparison of Linear System Methods

Method	Matrix Requirement	Type	Main Cost	Best Use
Gaussian Elimination	General nonsingular	Direct	$O(n^3)$	Dense systems
$PA=LU$	General nonsingular	Direct	$O(n^3)$	Stable general solver
Cholesky	Symmetric positive definite	Direct	$O(n^3)$	Dense SPD systems
Jacobi	Often diagonally dominant	Iterative	One matrix-vector step	Parallel sparse solving
Gauss-Seidel	Often diagonally dominant	Iterative	One sequential sweep	Sparse systems
SOR	Suitable convergent splitting	Iterative	One weighted sweep	Accelerated stationary iteration
Conjugate Gradient	Symmetric positive definite	Iterative	Matrix-vector products	Large sparse SPD systems

Essential Concepts

1. Linear System: A system written as $Ax=b$.
2. Gaussian Elimination: Converts a system into upper-triangular form.
3. Pivot: The diagonal value used during elimination.
4. Back Substitution: Solves an upper-triangular system from bottom to top.
5. LU Factorization: Expresses a matrix as $A=LU$.
6. Forward Substitution: Solves a lower-triangular system from top to bottom.
7. Residual: The vector $r=b-Ax_a$.
8. Forward Error: The difference between the approximate and exact solutions.
9. Backward Error: The amount by which the computed solution fails to satisfy the system.
10. Condition Number: Measures the sensitivity of a linear system.
11. Partial Pivoting: Selects the largest available pivot in the current column.
12. Permutation Matrix: Records row exchanges.
13. $PA=LU$ Factorization: Stable LU factorization with row permutations.
14. Jacobi Method: Uses only values from the previous iteration.
15. Gauss-Seidel Method: Immediately uses newly calculated values.
16. SOR: Applies a relaxation factor to Gauss-Seidel updates.
17. Strict Diagonal Dominance: A sufficient condition for Jacobi and Gauss-Seidel convergence.
18. Spectral Radius: Determines convergence of stationary iterative methods.
19. Sparse Matrix: Contains mostly zero entries.
20. Fill-In: New nonzero entries created during elimination.
21. Symmetric Positive-Definite Matrix: Satisfies $A^T=A$ and $x^TAx>0$.
22. Cholesky Factorization: Expresses an SPD matrix as $A=R^TR$.
23. Conjugate Gradient Method: Solves large sparse SPD systems.
24. Preconditioning: Improves convergence by transforming the system.
25. Jacobian Matrix: Contains the partial derivatives of a vector-valued function.
26. Multivariate Newton’s Method: Solves a Jacobian system at every iteration.
27. Broyden’s Method: Approximates the Jacobian or its inverse using rank-one updates.