Linear Systems

Gaussian Elimination

Solving linear systems is really important, how do we do it quickly?

The most trivial way to solve $Ax=b$ is with our explicit formula for $A^{-1}$ (calculating determinants recursively), but this is uselessly slow ($O(n!)$ if done completely naively)

While the inverse method is nice in theory, Gaussian elimination (the manual method, named after the same dude but unrelated to the distribution), is faster

Let’s be precise about the properties of this algorithm

1. It works for systems of equations with non-singular coefficient matrices
2. We can achieve good accuracy by picking “stronger” pivots (larger)

This algorithm runs in $O(n^3)$ and gives an exact solution (barring floating-point errors)

Iterative Methods

Iterative methods are promising because they allow one to improve approximations bit by bit, which allows for flexible trade-offs in time and accuracy

Assume we pick $S,T$ such that $S$ is invertible and one can write $A=S-T$
We can write $x=Mx+c$ where $M=S^{-1}T$ and $c=S^{-1}b$

This looks like a fixed-point problem, so we attempt the scheme $x_{k+1}=M{x}_k+c$ and hope that the error decreases

1. This scheme relies only on matrix-vector multiplication and addition
2. If $M$ is sparse then we can do this faster
3. This is only helpful if $S$ is cheap to invert

We can quickly derive $e_k=M^k{e}_0$ so we need to show $\Vert M^k{e}_0\Vert \to 0$ as $k\to \infty$

To do this, we need to first define a matrix norm

We note that all norms are equivalent,
Theorem: Let $\Vert \cdot {\Vert }_{\alpha }$, $\Vert \cdot {\Vert }_{\beta }$ be norms for vectors in ${\mathbb{R}}^n$, then there exists $M,m\ge 0$ such that $m\Vert x{\Vert }_{\alpha }\le \Vert x{\Vert }_{\beta }\le M\Vert x{\Vert }_{\alpha }$ for all $x\in {\mathbb{R}}^n$

Corollary: $\Vert x_k{\Vert }_{\alpha }\to 0⟺\Vert x_k{\Vert }_{\beta }\to 0$ as $k\to \infty$

Example:
$\Vert A\Vert ={\sup }_x\frac{\Vert Ax\Vert }{\Vert x\Vert }$

We can prove each property relatively easily, for example,
$\Vert A+B\Vert ={\sup }_x\frac{\Vert (A+B)x\Vert }{\Vert x\Vert }\le {\sup }_x\frac{\Vert Ax\Vert +\Vert Bx\Vert }{\Vert x\Vert }\le {\sup }_x\frac{\Vert Ax\Vert }{\Vert x\Vert }+{\sup }_x\frac{\Vert Bx\Vert }{\Vert x\Vert }=\Vert A\Vert +\Vert B\Vert$

We can further derive some useful identities for this norm,
$\Vert A{\Vert }_{\infty }={\max }_{i=1}^n{\sum }_{j=1}^n|A_{ij}|$
$\Vert A{\Vert }_1=\Vert A^T{\Vert }_{\infty }$
$\Vert A{\Vert }_2={\lambda }^{1/2}$ where $\lambda$ is the largest eigenvalue of $A^{T}A$

Now let’s build up what we need
Does $M^k\to 0$? If it is diagonalizable, we can write $M^k=R{D}^k{R}^{-1}$ where $D$ is diagonalizable and contains the eigenvalues ${\lambda }_i$, so $M^k\to 0⟺\rho (M)={\max }_{i=1}^n|{\lambda }_i|<1$

$\Vert e_k\Vert =\Vert M^k{e}_0\Vert$
$\le \Vert R{\Vert }_{\infty }\Vert R^{-1}{\Vert }_{\infty }\Vert D^k{\Vert }_{\infty }\Vert e_0{\Vert }_{\infty }$
$\le C\cdot \rho (M{)}^k$

${\lim }_{k\to \infty }\frac{\Vert e_{k+1}{\Vert }_{\infty }}{\Vert e_k{\Vert }_{\infty }}\le \rho (M)$
I.e. at least linear convergence if $\rho (M)<1$

To obtain our convergence theorem’s error bound, we use a few important properties,

• $\rho (A)\le \Vert A\Vert$ for any norm
• $\rho (A)<1⟹I-A$ is invertible
• If $\rho (A)<1$, then $\Vert (I-A{)}^{-1}\Vert <\frac{1}{1-\Vert A\Vert }$

Deriving it from here is pretty simple, and I’m going to omit it

Theorem: Let $A=S-T$ where $S$ is non-singular and define $x_k$ by $S{x}_k=T{x}_{k-1}+b$. Then $x_k\to x_{\ast }$ where $A{x}_{\ast }=b$ if and only if $\rho (S^{-1}T)<1$. Furthermore, if $\Vert S^{-1}T\Vert <1$ then $\Vert x_k-x_{\ast }\Vert \le \frac{\Vert S^{-1}T\Vert }{1-\Vert S^{-1}T\Vert }\Vert x_k-x_{k-1}\Vert$.

As we’ve said,
$x_{k+1}=M{x}_k+c$ where $M=S^{-1}T$

This is great because it gives us an error bound, however we still require an easy to compute $S^{-1}$ for the scheme to be practical

One way begins by splitting $A=L+D+U$ (lower, diagonal, upper) and defines $S=D$, $T=-L-U$
This is the most trivial solution, since diagonal matrices are very easy to invert, and is called the Jacobi scheme

However, we note that the order on the rows of $A$ effects the convergence of our method under the Jacobi scheme! This is bad

We can find that convergence of $A$ under the Jacobi scheme rests on strict diagonal dominance of $A$, where the magnitude of each diagonal is greater than the sum of the magnitudes of the remaining values in its row

Theorem: Let $A$ be strictly diagonally dominant, then the Jacobi method for solving $Ax=b$ converges

An alternative to the Jacobi method is the Gauss-Seidel method, which splits $A$ into $S=L+D$ and $T=-U$

Lower triangular matrices are easy to invert and our same theorem applies