<< numerical_analysis

Linear Equations Solving

1. Linear Equation Solving by Elimination

We’ve kown that the equations with $n$-variables $$\{\begin{array}{l}{a}_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1n}{x}_n=b_1\\ a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2n}{x}_n=b_2\\ \cdots \\ a_{n1}{x}_1+a_{n2}{x}_2+\cdots +a_{nn}{x}_n=b_n\end{array}$$

can be written in matrix form $A\mathit{x}=\mathit{b}$ $$\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right]$$

where $A$ is coefficient matrix, $\mathit{b}$ is constant vector, and $\mathit{x}$ is solution vector.

Cramer’s rule guarantees that for coefficient matrix $A$ with $detA\ne 0$, the solution exists and is unique: $$x_i=\frac{D_i}{D},i=1,2,\dots ,n$$

where $D_i$ can be generated by replacing the $i$-th column of coefficient matrix $A$ with constant vector $\mathit{b}$.

Although Cramer’s rule plays an important role in thoery of linear equations, the time complexity $O(n!\times n)$ of computing determinant makes it not acceptable to apply Cramer’s rule in large scale of linear systems.

# Gaussian Elimination

Gaussian elimination, instead, applies elementary row operation to augmented matrix ${\left[\begin{array}{cc}A& b\end{array}\right]}^T$. It converts the matrix into an upper triangular matrix, after these operations, the equation can be easily solved by substituting (bottom-up). The time complexity of this approach is $O(n^3)$.

The Gaussian elimination requires

$${\mathrm{\Delta }}_k=\left|\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1k}\\ a_{21}& a_{22}& \cdots & a_{2k}\\ ⋮& ⋮& \ddots & ⋮\\ a_{k1}& a_{k2}& \cdots & a_{kk}\end{array}\right|\ne 0$$

for any $k<n$. While the solution exists only requires $detA\ne 0$. This indicates that Gaussian elimination cannot be applied to some solvable functions with $a_{kk}^{(k-1)}=0$ during the elimination. Besides, the round-off error will increase for small $a_{kk}^{(k-1)}$.

# Gaussian Elimination with Column Pivoting

To overcome the flaw of Gaussian elimination, the column pivoting operation for each iteration in elimination is proposed. Formally, for each $k=1,2,\dots ,n-1$, we search for the element with maximal absolute value $|a_{m,k}^{(k-1)}|$ within $|a_{k,k}^{k-1}|,|a_{k+1,k}^{k-1}|,\dots |a_{n,k}^{k-1}|$, then we swap row $k,m$ before elimination.

We can solve any linear equations with $detA\ne 0$ with column pivoting. Since the total time complexity of comparison is $O(n^2)$, the time complexity of Gaussian elimination with column pivoting is still $O(n^3)$.

# Gauss-Jordan Elimination

Gaussian elimination execute substitution process right after the matrix is converted to upper triangular matrix, while we can do further elimination to convert it into a diagonal matrix, which is called Gauss-Jordan elimination.

2. Linear Equation Solving by Matrix Decomposition

# Relation between Line Operation and Matrix Decomposition

An elementary line operation is equivalent to multiply a elementary matrix. The process of Gaussian elimination is equivalent to multiply a matrix $T$ to $A$, results in a upper triangular matrix: $$TA=U$$

Multiply $L=T^{-1}$ on both sides: $$A=LU$$ where $$L=T^{-1}=\left[\begin{array}{c}1\\ l_{21}& 1\\ ⋮& ⋮& \ddots \\ l_{n1}& l_{n2}& \cdots & 1\end{array}\right]$$ can be proved to be a lower triangular matrix. The form $A=LU$ indicates that we devide the matrix $A$ into the multiplication of two matrixes $L$, $U$. where $L$ and $U$ are lower and upper triangular matrix respectively. Equivalently, we can also decompose $A$ as $A=UL$.

Since you can move the coefficients from one matrix to another, there are infinite decompositions of form $A=LU$. If $L$ is contrainted to be the unit lower triangular matrix, that unique decomposition is called Doolittle decomposition. If $U$ is contrainted to be the unit upper triangle matrix, that unique decomposition is called Crout decomposition.

Besides, there are other forms of decomposition such as $A=LD{M}^T$, where the $L,M$ is lower triangular matrix (hence $M^T$ is upper triangular matrix), and $D$ is diagonal matrix. Specially, if $A$ is a positive definite matrix, we have $M=L$, $A=LD{L}^T$, this is what we called $LD{L}^T$ decomposition or refined Cholesky decomposition. Let $P=L\sqrt{D}$ we have $A=P{P}^T$, which is called Cholesky decomposition.

The matrix is decomposed into two or three triangular/diagonal matrixes after decomposition, which can be easily solved by substitution. Note that decomposition process is independent of $\mathit{b}$. Hence if there are multiple equations with same $A$ but various $\mathit{b}$, this approach reduces the time complexity by cancelling the repetative processing of $A$.

# LU Decomposition

Compare elements to determine the matrix after Doolittle decomposition:

$$A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right]=\left[\begin{array}{c}1\\ l_{21}& 1\\ ⋮& ⋮& \ddots \\ l_{n1}& l_{n2}& \cdots & 1\end{array}\right]\left[\begin{array}{cccc}{u}_{11}& u_{12}& \cdots & u_{1n}\\ & u_{22}& \cdots & u_{2n}\\ & & \ddots & ⋮\\ & & & u_{nn}\end{array}\right]$$

$$a_{11}=\sum _{r=1}^n{l}_{1r}{u}_{r1}=\left[\begin{array}{cccc}1& 0& \cdots & 0\end{array}\right]\left[\begin{array}{c}{u}_{11}\\ 0\\ ⋮\\ 0\end{array}\right]=\sum _{r=1}^1{l}_{1r}{u}_{r1}=u_{11}$$

$$a_{1j}=\sum _{r=1}^n{l}_{1r}{u}_{rj}=\left[\begin{array}{cccc}1& 0& \cdots & 0\end{array}\right]\left[\begin{array}{c}{u}_{1j}\\ u_{2j}\\ ⋮\\ 0\end{array}\right]=\sum _{r=1}^1{l}_{1r}{u}_{rj}=u_{1j}$$ Hence $a_{ij}=u_{1j}$, $j=1,2,\dots ,n$.

$$\begin{array}{rl}{a}_{kj}& =\sum _{r=1}^n{l}_{kr}{u}_{rj}=\left[\begin{array}{cccccccc}{l}_{k1}& l_{k2}& \cdots & l_{k,k-1}& 1& 0& \cdots & 0\end{array}\right]\left[\begin{array}{c}{u}_{1j}\\ ⋮\\ u_{jj}\\ 0\\ ⋮\\ 0\end{array}\right]\\ & =\sum _{r=1}^n{l}_{kr}{u}_{rj}=\sum _{r=1}^{k-1}{l}_{kr}{u}_{rj}+u_{kj}\end{array}$$

Now we have the algorithm for Doolittle decomposition:

Refer to https://github.com/revectores/revector for the C++ implementation of Doolittle decomposition.

Refer to https://github.com/revectores/LU-decomposition-impls for some parallelized Doolittle decomposition implemented in C.

Similarly, the algorithm of Crout decomposition:

Refer to https://github.com/revectores/revector for the C++ implementation of Crout decomposition.

# Decomposition of Positive Definite Matrix

As mentioned above, we have Cholesky decomposition $A=L{L}^T$ for positive definitive matrix

In practice we often decompose as form $A=LD{L}^T$ instaed to remove square root computation.

Proof. To prove the existence of $LD{L}^T$ decomposition, we first apply Doolittle decomposition to matrix and extract the diagonal line of $U$: $$\begin{array}{rl}A=LU& =\left[\begin{array}{c}1\\ l_{21}& 1\\ ⋮& ⋮& \ddots \\ l_{n1}& l_{n2}& \cdots & 1\end{array}\right]\left[\begin{array}{cccc}{u}_{11}& u_{12}& \cdots & u_{1n}\\ & u_{22}& \cdots & u_{2n}\\ & & \ddots & ⋮\\ & & & u_{nn}\end{array}\right]\\ & =\left[\begin{array}{c}1\\ l_{21}& 1\\ ⋮& ⋮& \ddots \\ l_{n1}& l_{n2}& \cdots & 1\end{array}\right]\left[\begin{array}{cccc}{u}_{11}& & & \\ & u_{22}& & \\ & & \ddots & \\ & & & u_{nn}\end{array}\right]\left[\begin{array}{cccc}{\bar{u}}_{11}& {\bar{u}}_{12}& \cdots & {\bar{u}}_{1n}\\ & {\bar{u}}_{22}& \cdots & {\bar{u}}_{2n}\\ & & \ddots & ⋮\\ & & & {\bar{u}}_{nn}\end{array}\right]\end{array}$$

Since $A$ is symmetric positive definite matrix, we have $u_{ii}>0$. We can prove $L=\bar{U}$ since $$A=LU=LD{\bar{U}}^T=A^T=\bar{U}(D{L}^T)$$ That is, $$\begin{array}{r}A=LD{L}^T=\left[\begin{array}{c}1\\ l_{21}& 1\\ ⋮& ⋮& \ddots \\ l_{n1}& l_{n2}& \cdots & 1\end{array}\right]\left[\begin{array}{cccc}{u}_{11}& & & \\ & u_{22}& & \\ & & \ddots & \\ & & & u_{nn}\end{array}\right]\left[\begin{array}{cccc}1& l_{21}& \cdots & l_{n1}\\ & 1& \cdots & l_{n2}\\ & & \ddots & ⋮\\ & & & 1\end{array}\right]\end{array}$$

# Condition Number of Matrix

For the non-singular matrix $A$, we define $${\kappa }_p(A)={\Vert A\Vert }_p{\Vert A^{-1}\Vert }_p$$ as the condition number of matrix $A$. Due to the equivalence of different norms, the condition numbers with different norms are also equivalent.

How the error $\delta x$ is influenced by small disturbance $\delta A$ in coefficient matrix $A$ is represented as $$\frac{\Vert \delta x\Vert }{\Vert x\Vert }\le \frac{{\kappa }_A{\displaystyle \frac{\Vert \delta A\Vert }{\Vert A\Vert }}}{1-{\kappa }_A{\displaystyle \frac{\Vert \delta A\Vert }{\Vert A\Vert }}}$$ How the error $\delta x$ is influenced by small disturbance $\delta \mathit{b}$ in constant vector $\mathit{b}$ is represented as $$\frac{\Vert \delta x\Vert }{\Vert x\Vert }\le {\kappa }_A\frac{\Vert \delta \mathit{b}\Vert }{\Vert \mathit{b}\Vert }$$ We call the matrix with large ${\kappa }_A$ ill-conditioned, which states “the small disturbance causes a huge difference in solution”.

3. Linear Equation Solving by Iteration

Transform the equations $AX=\mathit{y}$ into $X=MX+\mathit{g}$, for any $X^{(0)}\in R^n$, construct the iteration $$X^{(k+1)}=M{X}^{(k)}+\mathit{g}$$ if the iteration series $\{X^{(k)}\}$ converges, the limit of iteration series $X^{\ast }$ is the solution of equations $AX=\mathit{y}$.

In practice, the iteration is terminated once we have ${\Vert X^{(k+1)}-X^{(k)}\Vert }_p<\epsilon$ for some pre-chosen error bound $\epsilon$, and pick $X^{(k+1)}$ as the approximation of solution.

The convergence of iteration is determined by spectral radius of iteration matrix $M$. Specifically, the iteration converges if and only if the spectral radius $\rho (M)<1$.

Proof. If $X^{\ast }$ is the solution of equation $AX=\mathit{y}$, then $$X^{\ast }=M{X}^{\ast }+\mathit{g}$$

$$\begin{array}{rl}{X}^{\ast }-X^{(k+1)}& =M(X^{\ast }-X^k)\\ & =M^2(X^{\ast }-X^{k-1})\\ & =\cdots \\ & =M^{k+1}(X^{\ast }-X^{(0)})\end{array}$$

${\displaystyle \underset{k\to \mathrm{\infty }}{lim}{M}^k=0}$ if and only if $\rho (M)<1$. Hence we define the matrix with spectral radius less than 1 convergent matrix.

That is, whether the linear equations converges depends on the property of iteration matrix, regardless of the solution $\alpha$ and $X^{(0)}$.

By definition, we have to compute all the eigenvalues to get the spectral radius of matrix. We may simplify this work by computing the norm ${\Vert A\Vert }_p$. Note that ${\Vert A\Vert }_p\ge \rho (A)$, if ${\Vert A\Vert }_p<1$, the iteration matrix must be convergent.

# Jacobi Iteration

$$\{\begin{array}{l}{a}_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1n}{x}_n=y_1\\ a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2n}{x}_n=y_2\\ \cdots \\ a_{n1}{x}_1+a_{n2}{x}_2+\cdots +a_{nn}{x}_n=y_n\end{array}$$

$$\{\begin{array}{l}{x}_1={\displaystyle \frac{1}{a_{11}}}\left(-a_{12}{x}_2-\cdots -a_{1n}{x}_n+y_1\right)\\ x_2={\displaystyle \frac{1}{a_{22}}}\left(-a_{21}{x}_2-\cdots -a_{2n}{x}_n+y_2\right)\\ \cdots \\ x_n={\displaystyle \frac{1}{a_{nn}}}\left(-a_{n1}{x}_2-\cdots -a_{nn}{x}_n+y_n\right)\end{array}$$

The iteration form $$\{\begin{array}{l}{x}_1^{(k+1)}={\displaystyle \frac{1}{a_{11}}}\left(-a_{12}{x}_2^{(k)}-\cdots -a_{1n}{x}_n^{(k)}+y_1\right)\\ x_2^{(k+1)}={\displaystyle \frac{1}{a_{22}}}\left(-a_{21}{x}_2^{(k)}-\cdots -a_{2n}{x}_n^{(k)}+y_2\right)\\ \cdots \\ x_n^{(k+1)}={\displaystyle \frac{1}{a_{nn}}}\left(-a_{n1}{x}_1^k-\cdots -a_{n,n-1}{x}_n^{(k)}+y_n\right)\end{array}$$ The matrix form of Jacobi iteration is $$X^{(k+1)}=B{X}^k+\mathit{g}$$

where $B=I-D^{-1}A,\mathit{g}=D^{-1}y$.

There is a shortcut to judge whether Jacobi iteration converges: if the coefficient matrix $A$ meets one of following condition:

• $$\begin{gathered} |a_{ii}|>{\displaystyle \sum _{j=1 \\ j\ne i}^n|a_{ij}|,i=1,2,\dots ,n} \end{gathered}$$.
• $$\begin{gathered} |a_{jj}|>{\displaystyle \sum _{i=1 \\ i\ne j}^n|a_{ij}|,j=1,2,\dots ,n} \end{gathered}$$

The Jacobi iteration must converge. This can be proved by showing the spectral radius of iteration matrix is less than 1.

==TODO: Add the proofs.==

# Gauss-Seidel Iteration

The Gauss-Seidel iteration use those new values computed in current iteration instead of last one: $$\{\begin{array}{l}{x}_1^{(k+1)}={\displaystyle \frac{1}{a_{11}}}\left(-a_{12}{x}_2^{(k)}-\cdots -a_{1n}{x}_n^{(k)}+y_1\right)\\ x_2^{(k+1)}={\displaystyle \frac{1}{a_{22}}}\left(-a_{21}{x}_2^{(k+1)}-\cdots -a_{2n}{x}_n^{(k)}+y_2\right)\\ \cdots \\ x_n^{(k+1)}={\displaystyle \frac{1}{a_{nn}}}\left(-a_{n1}{x}_1^{(k+1)}-\cdots -a_{n,n-1}{x}_{n-1}^{(k+1)}+y_n\right)\end{array}$$ Denote $$D=\left[\begin{array}{c}{a}_{11}\\ & a_{22}\\ & & \ddots \\ & & & a_{nn}\end{array}\right],L=\left[\begin{array}{c}0\\ a_{21}& 0\\ ⋮& ⋮& \ddots \\ a_{n1}& a_{n2}& \cdots & 0\end{array}\right],U=\left[\begin{array}{cccc}0& a_{12}& \cdots & a_{1n}\\ & 0& \cdots & a_{2n}\\ & & \ddots & ⋮\\ & & & 0\end{array}\right]$$ We have $$AX=(D+L+U)X=(D+L)X+UX=\mathit{y}$$ That is, $$(D+L)X=-UX+\mathit{y}$$ Hence $$X^{(k+1)}=-(D+L{)}^{-1}U{X}^{(k)}+(D+L{)}^{-1}\mathit{y}$$ Denote $S=-(D+L{)}^{-1}U,\mathit{f}=(D+L{)}^{-1}\mathit{y}$, the Gauss-Seidel iteration can be expressed by $$X^{k+1}=S{X}^k+\mathit{f}$$

# Successive Over-Relaxation

The SOR method is a variant of Gauss-Seidel method resulting in faster convergence: $$\{\begin{array}{l}{x}_1^{(k+1)}=(1-\omega )x_1^{(k)}+\omega (b_{12}{x}_2^k+\cdots +b_{1n}{x}_n^{(k)}+g_1)\\ x_2^{(k+1)}=(1-\omega )x_2^{(k)}+\omega (b_{21}{x}_2^k+\cdots +b_{2n}{x}_n^{(k)}+g_2)\\ \cdots \\ x_2^{(k+1)}=(1-\omega )x_n^{(k)}+\omega (b_{n1}{x}_2^{(k+1)}+\cdots +b_{n,n-1}{x}_{n-1}^{(k+1)}+g_n)\end{array}$$ where $\omega$ is relaxation factor.

The matrix form of SOR is $X^{(k+1)}=S_{\omega }{X}^{(k)}+\mathit{f}$, where $$\begin{array}{rl}{S}_w& =(I+\omega D^{-1}L{)}^{-1}[(1-\omega )I-\omega D^{-1}U]\\ f& =\omega (I+\omega D^{-1}L)Y\end{array}$$

Proof. $$\begin{array}{rl}{X}^{(k+1)}& =(1-\omega )X^{(k)}+\omega (\tilde{L}{X}^{(k+1)}+\tilde{U}{X}^{(k)}+\mathit{g})\\ (I-\omega \tilde{L})X^{(k+1)}& =((1-\omega )I+\omega \tilde{U})X^{(k)}+\omega \mathit{g}\\ X^{(k+1)}& =(I-\omega \tilde{L}{)}^{-1}((1-\omega )I+\omega \tilde{U})X^{(k)}+(I-\omega \tilde{L}{)}^{-1}\omega \mathit{g}\end{array}$$ where $\tilde{L}=-D^{-1}L,\tilde{U}=-D^{-1}U$.

The necessary condition that SOR converges is $0<\omega <2$. Specially, if $A$ is the positive determined matrix, $0<\omega <2$ is also the sufficient condition.

The iteration is called

• Under-relaxation iteration if $\omega \in (0,1)$.
• Gauss-Seidel iteration if $\omega =1$.
• Over-relaxation iteration if $\omega \in (1,2)$.

The best $\omega$ (with fastest convergence speed) is hard to determine for specific equation.