\section{LINEAR SYSTEMS OF EQUATIONS} Consider the linear system of equations $$\eqalign{ 1x_1 + 2x_2 + 3x_3 + 4x_4 & = 5 \cr 2x_1 + 6x_2 + 9x_3 + 8x_4 & = 10 \cr -1x_1\phantom{{}+6x_2} + 4x_3 - 5x_4 & = -6 \cr 1x_1\phantom{{}+6x_2} + 12x_3 + 2x_4 & = 3 \cr} \hbox{\quad or\quad} Ax = \pmatrix{1 & 2 & 3 & 4\cr 2 & 6 & 9 & 8\cr -1 & 0 & 4 & -5\cr 1 & 0 & 12 & 2\cr} x =\pmatrix{5\cr 10\cr -6\cr 3\cr}=b $$ in matrix form. Gaussian elimination, a technique for solving $Ax=b$, systematically eliminates variables from equations, producing an equivalent linear system that is trivial to solve. The elimination operations can be expressed in terms of matrix operations: $$\eqalign{ \pmatrix{1 & 0 & 0 & 0 \cr -2 & 1 & 0 & 0 \cr 1 & 0 & 1 & 0 \cr -1 & 0 & 0 & 1 \cr} \pmatrix{1 & 2 & 3 & 4 \cr 2 & 6 & 9 & 8 \cr -1 & 0 & 4 & -5 \cr 1 & 0 & 12 & 2 \cr} x & = L_1Ax=A^{(2)}x = \pmatrix{1 & 2 & 3 & 4 \cr 0 & 2 & 3 & 0 \cr 0 & 2 & 7 & -1 \cr 0 & -2 & 9 & -2 \cr} x \cr & = L_1b = \pmatrix{5\cr 0\cr -1\cr -2\cr}, \cr \pmatrix{1 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 \cr 0 & -1 & 1 & 0 \cr 0 & 1 & 0 & 1 \cr} \pmatrix{1 & 2 & 3 & 4 \cr 0 & 2 & 3 & 0 \cr 0 & 2 & 7 & -1 \cr 0 & -2 & 9 & -2 \cr} x & = L_2A^{(2)}x=A^{(3)}x = \pmatrix{1 & 2 & 3 & 4 \cr 0 & 2 & 3 & 0 \cr 0 & 0 & 4 & -1 \cr 0 & 0 & 12 & -2 \cr} x \cr & = L_2L_1b = \pmatrix{5\cr 0\cr -1\cr -2\cr}, \cr \pmatrix{1 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 \cr 0 & 0 & 1 & 0 \cr 0 & 0 & -3 & 1 \cr} \pmatrix{1 & 2 & 3 & 4 \cr 0 & 2 & 3 & 0 \cr 0 & 0 & 4 & -1 \cr 0 & 0 & 12 & -2 \cr} x & = L_3A^{(3)}x=Ux = \pmatrix{1 & 2 & 3 & 4 \cr 0 & 2 & 3 & 0 \cr 0 & 0 & 4 & -1 \cr 0 & 0 & 0 & 1 \cr} x \cr & = L_3L_2L_1b = \pmatrix{5\cr 0\cr -1\cr 1\cr}={\tilde b}. \cr} $$ Now $Ux={\tilde b}$ can be easily solved for $x$ using backward substitution. First $x_4=1$. Then substituting $x_4$ in the third equation gives $x_3=0$. Using $x_3=0$, $x_4=1$ in the second equation gives $x_2=0$. Finally $x_1=1$ from the first equation. \noindent {\bf Observation}: $U$ is upper triangular, and the $L_i$ are unit lower triangular (ones on the diagonal). \vfil $L_3L_2L_1A=U \Rightarrow A = (L_1^{-1}L_2^{-1}L_3^{-1})U=LU$, where $L$ is unit lower triangular. $A=LU$ is called the $LU$ factorization (or decomposition) of $A$. Here, $$ A= \pmatrix{1 & 2 & 3 & 4 \cr 2 & 6 & 9 & 8 \cr -1 & 0 & 4 & -5 \cr 1 & 0 & 12 & 2 \cr} = \pmatrix{1 & 0 & 0 & 0 \cr 2 & 1 & 0 & 0 \cr -1 & 1 & 1 & 0 \cr 1 & -1 & 3 & 1 \cr} \pmatrix{1 & 2 & 3 & 4 \cr 0 & 2 & 3 & 0 \cr 0 & 0 & 4 & -1 \cr 0 & 0 & 0 & 2 \cr}=LU. $$ Note that the elements in $L$ are simply the multipliers used to reduce $A$ to $U$. Thus when reducing $A$ to $U$, $L$ is obtained for free. Note that both $L$ and $U$ can be stored on top of the original matrix $A$, so no extra storage is required to compute the $LU$ factorization. Given $A=LU$, to solve $Ax=b$: \item{1)} solve $Ly=b$ (forward substitution), \item{2)} solve $Ux=y$ (backward substitution). \bigskip \centerline{\bf Floating point operation counts} {\baselineskip=18pt \tabskip=20pt \halign{ #\hfil&#\hfil\cr To reduce $A$ to $U$ and simultaneously & \cr apply $L^{-1}$ to $p$ right hand sides: & ${2 \over 3}n(n^2-1)+n(n-1)(p-{1\over 2})$ \cr To compute $LU$ factorization: & ${2 \over 3}n(n^2-1)+{n(n-1)\over 2}$ \cr Forward substitution ($Ly=b$): & $n(n-1)$ \cr Backward substitution ($Ux=y$): & $n^2$ \cr To compute $A^{-1}$: & $2n^3-{3\over 2}n^2+{n\over 2}$ \cr Product $A^{-1}b$: & $(2n-1)n$ \cr }} {\bf Conclusion}: the most efficient way to solve a linear system $Ax=b$ for one or many right hand sides $b$ is to compute the $LU$ decomposition of $A$, and then use forward or backward substitution for each right hand side. \medskip\hrule\smallskip ``There is nothing you can do with the inverse that you can't do better without it.'' (Cleve Moler) \smallskip\hrule\bigskip \centerline{\bf PIVOTING and SCALING}\bigskip If $A_{kk}^{(k)}=0$ for some $k$, then $k$ must be interchanged with some row $m>k$. This interchange, and subsequent introduction of 0s into the $k$th column below the diagonal, is called pivoting. With row interchanges, an $LU$ decomposition is always possible. \noindent {\bf Theorem}: For any $A \in \Enn$, $\exists$ a permutation matrix $P$ (matrix of 0s and 1s with exactly one 1 in each row and column), a unit lower triangular $L \in \Enn$, and an upper triangular $U \in \Enn$, such that $PA = LU$. \noindent Example. Consider solving the linear system $$\eqalign{ 1.00\cdot 10^{-4}x_1 + 1.00x_2 & = 1.00\cr 1.00x_1 + 1.00x_2 & = 2.00\cr},$$ whose exact solution is $\pmatrix{1.00010\cr .99990}$, using 3-digit arithmetic. Pivoting on $1.00\cdot 10^{-4}$ gives $$\pmatrix{1.00\cdot 10^{-4} & 1.00 & \vdots & 1.00 \cr 0 & -1.00\cdot 10^4 & \vdots & -1.00\cdot 10^4} \Longrightarrow x_2=1.00,\ x_1=0.00.$$ Pivoting on $1.00$ gives $$\pmatrix{1.00\cdot & 1.00 & \vdots & 2.00 \cr 0 & 1.00 & \vdots & 1.00} \Longrightarrow x_2=x_1=1.00.$$ Now scale the original system to get $$ \eqalign{ 1.00\cdot 10^1x_1+1.00\cdot 10^5x_2 & = 1.00\cdot 10^5\cr 1.00x_1+1.00x_2 & = 2.00\cr}.$$ \leftline{Pivoting on $1.00\cdot 10^1$ gives $x_2=1.00,\ x_1=0.00.$} \leftline{Pivoting on $1.00$ gives $x_2=x_1=1.00$.} \leftline{Potential rule? Pivot on row which maximizes $\displaystyle \bigl | A_{mk}^{(k)} \bigr | \big / \bigl \| A_m^{(k)} \bigr \|_\infty.$} Consider these examples (R.W. Hamming) where $fl(1+\epsilon)>1, fl(1+10\epsilon^2)=1$: Scaling $$A=\pmatrix{3 & 1/\epsilon & 1/\epsilon\cr 1 & 2\epsilon & \epsilon\cr 2 & \epsilon & \epsilon}$$ by rows first and then pivoting produces $$\pmatrix{3\epsilon & 1 & 1\cr 1 & 2\epsilon & \epsilon\cr 2 & \epsilon & \epsilon} \longrightarrow 2 \pmatrix{0 & 1 & 1 \cr 0 & 3\epsilon/2 & \epsilon/2 \cr 1 & \epsilon/2 & \epsilon/2 \cr} \longrightarrow 2 \pmatrix{0 & 1 & 1 \cr 0 & 0 & -\epsilon \cr 1 & \epsilon/2 & \epsilon/2 \cr},$$ from which back substitution will be accurate. Scaling $A$ by columns first and then pivoting produces $$\pmatrix{3 & 1 & 1\cr 1 & 2\epsilon^2 & \epsilon^2\cr 2 & \epsilon^2 & \epsilon^2} \longrightarrow \pmatrix{3 & 1 & 1 \cr 0 & -1/3 & -1/3 \cr 1 & -2/3 & -2/3 \cr} \longrightarrow 2 \pmatrix{3 & 1 & 1 \cr 0 & -1/3 & -1/3 \cr 0 & 0 & 0 \cr},$$ from which back substitution is impossible! \noindent {\bf Bottom line}: reliable pivoting and scaling strategies are not known. \noindent {\bf What to do}? Obtain a posteriori error estimates (backward error analysis). \bigskip \centerline{\bf BACKWARD ERROR ANALYSIS}\bigskip Effect of perturbation in the data: Let $A \in \Enn$ be invertible, $b$, $\delta b \in \En$, $Ax=b$, \hfil\break $A(x + \delta x) = b + \delta b$, where $\| \delta b \| \ll \| b \|$. $\delta x = A^{-1} \delta b \Rightarrow \| \delta x \| = \| A^{-1} \delta b \| \leq \| A^{-1} \| \| \delta b \|.$ $\| b \| = \| Ax \| \leq \| A \| \| x \| \Rightarrow \| x \| \geq {\| b \| / \| A\|}.$ Combining these two inequalities gives $${\| \delta x \| \over \| x \|} \leq { \| A^{-1} \| \| \delta b \| \over \| b \| / \| A \|} = \| A \| \| A^{-1} \| {\| \delta b \| \over \| b \|} = \hbox{ cond } A \;{\| \delta b \| \over \| b \|}.$$ \noindent {\bf Definition}: $\hbox{cond } A=\| A \| \| A^{-1} \|$ is called the condition number of $A$ (in the norm $\|\cdot\|$). \noindent {\bf Theorem}: Let $A \in \Enn$ be invertible, $\delta A \in \Enn$, $b$, $\delta b \in \En$, $Ax=b$, $(A + \delta A)(x+\delta x) = b + \delta b$, $\| \delta A \| < {1 / \| A^{-1} \|}$. Then, $$ {\| \delta x \| \over \| x \|} \leq {\hbox{cond } A \over 1 - \hbox{cond } A\;{\| \delta A \| \over \| A \|}} \left( {\| \delta b \| \over \| b \|} + {\| \delta A \| \over \| A \|} \right).$$ If $\hat{x}$ is an estimate of the exact solution $x$ to $Ax=b$, the residual $r=b-A \hat{x}$. \par $r=b-A \hat{x}=A(x- \hat{x}) \Rightarrow \| \hat{x} - x \| = \| -A^{-1} r \| \leq \| A^{-1} \| \| r \| \Rightarrow$ $${\| \hat{x} - x \| \over \| x \|} \leq {\| A^{-1} \| \| r \| \over \| b \| / \| A \|} = \hbox{cond } A\;{\| r \| \over \| b \|}. \qquad \hbox{(This inequality is sharp.)}$$ \noindent {\bf Moral}: A small residual does not imply an accurate estimate! \noindent {\bf Theorem}: The solution $\hat{x}$ computed by Gaussian elimination with partial pivoting satisfies \hfil\break $(A + \delta A)\hat{x}=b$, where $\| \delta A \|_\infty \leq 1.01 (n^3 + 3n^2)\rho\| A \|_\infty \mu,$ and $$\rho = \max_{i,j,k} {|A^{(k)}_{ij}| \over \| A \|_\infty}.$$ \noindent {\bf Folklore}: $\| \delta A \|_\infty \leq n\mu \| A \|_\infty.$ \vfil\eject \section{MATRIX DECOMPOSITION THEOREMS} \noindent {\it Theorem 1.} Let A be a real $n\times n$ matrix and $A_k$ the submatrix of A formed from the first $k$ rows and columns. Let $\det A_k\ne 0$ for $k=1,\ldots ,n$. Then $\exists$ a unique unit (ones on the diagonal) lower triangular matrix $L$ and a unique upper triangular matrix $U$ such that $A=LU$. \noindent {\it Theorem 2.} Under the hypothesis of Theorem 1, $\exists$ a unique lower triangular matrix $L$, a unique diagonal matrix $D$, and a unique upper triangular matrix $U$ such that $A=LDU$. \noindent {\it Theorem 3.} Let $A$ be symmetric $\left(A=A^t\right)$ and satisfy the hypothesis of Theorem 1. Then $\exists$ a unique unit lower triangular matrix $L$ and a unique diagonal matrix $D$ such that $A=LDL^t$. \noindent {\it Theorem 4.} Let $A$ be symmetric and positive definite $\left(x^tAx>0\ \hbox {for all}\ x\ne 0\right)$. Then $\exists$ a unique lower triangular matrix $G$ with positive diagonal elements such that $A=GG^t$. $(A=GG^t$ is the Cholesky factorization of $A.)$ \noindent {\it Definition.} A permutation matrix $P$ is an $n\times n$ matrix consisting of 0's and 1's, such that every row and column contains exactly one 1. (Intuitively, $Py$ simply permutes the elements of the vector $y$.) \noindent {\it Theorem 5.} Let $A$ be an $n\times n$ matrix. Then $\exists$ a permutation matrix $P$, a unit lower triangular matrix $L$, and an upper triangular matrix $U$ such that $PA=LU$. \medskip\hrule\medskip \centerline {DIRECT ALGORITHMS} \noindent Suppose it is known beforehand that $A$ has an $LU$ factorization without interchanging rows. Then, carefully examining the equation $$\def\a#1#2{a_{#1#2}} \def\l#1#2{\ell_{#1#2}} \def\u#1#2{u_{#1#2}} \pmatrix {\a11&\a12&\a13&\cdots&\a1n\cr \a21&\a22&\a23&\cdots&\a2n\cr \a31&\a32&\a33&\cdots&\a3n\cr \vdots&\vdots&\vdots&\ddots&\vdots\cr \a n1&\a n2&\a n3&\cdots&\a nn\cr} = \pmatrix {1&0&0&\cdots&0\cr \l21&1&0&\cdots&0\cr \l31&\l32&1&\cdots&0\cr \vdots&\vdots&\vdots&\ddots&\vdots\cr \l n1&\l n2&\l n3&\cdots&1\cr} \pmatrix {\u11&\u12&\u13&\cdots&\u1n\cr 0&\u22&\u23&\cdots&\u2n\cr 0&0&\u33&\cdots&\u3n\cr \vdots&\vdots&\vdots&\ddots&\vdots\cr 0&0&0&\cdots&\u nn}$$ shows that $L$ and $U$ can be determined recursively. Precisely, $$\displaylines {\hbox {for } k = 1,\ldots, n\cr u_{km} = a_{km} - \sum^{k-1}_{j=1}\ \ell_{kj} u_{jm}, \quad m = k,\ldots, n\cr \ell_{pk} = \left(a_{pk} - \sum^{k-1}_{j=1} \ell_{pj} u_{jk}\right)\Big/ u_{kk},\quad p = k + 1, \ldots, n\cr}$$ These formulas (known as the Crout or Doolittle algorithm) require exactly the same computational effort as Gaussian elimination, but are inappropriate unless it is known a priori that pivoting is not necessary. \vfil\eject \noindent{\bf Definition}. Let $u\in C^n$, $u\ne0$. A Householder reflection is a matrix of the form $$P=I-{1\over\beta}uu^\ast,\qquad \hbox{where }\beta={\|u\|^2\over2}.$$ $A^\ast$ denotes the conjugate transpose $(\bar A)^t$ of a complex matrix $A\in C^{m\times n}$. \noindent{\bf Theorem}. Householder reflections $P$ are Hermitian ($P^\ast=P$), unitary ($P^\ast P=I$), and self-inverse ($P^{-1}=P$). \noindent{\bf Theorem}. Let $a\in C^n$, $a\ne0$, $\alpha=e^{i\;{\rm arg}(a_1)}\|a\|$, $u=a+\alpha e_1$, $\beta=u_1\bar\alpha$. Then $P = I - {1\over\beta}uu^\ast$ is a Householder reflection and $Pa=-\alpha e_1$. For any $b\in C^n$, $Pb=b-(u^\ast b/\beta)u$. Proof. $u^\ast u=\bigl(a^\ast+\bar\alpha e_1^t\bigr) (a+\alpha e_1) = \|a\|^2 + 2\Re(a_1)\bar\alpha + \|a\|^2 =2(\|a\|^2 +\Re(a_1)\bar\alpha )= 2(\|a\|^2 + a_1\bar\alpha) = 2(a_1+\alpha)\bar\alpha = 2u_1\bar \alpha = 2\beta$ so $P=I-{1\over\beta}uu^\ast$ is a Householder reflection. Now $\displaystyle Pa=a-(u^\ast a/\beta)u = a-{\bigl(a^\ast+\bar\alpha e_1^t\bigr)a\over\beta} (a+\alpha e_1) = a - {\bigl(\|a\|^2 + \bar \alpha a_1\bigr)\over\beta} (a+\alpha e_1) = a - (a+\alpha e_1) = -\alpha e_1$. \QED \noindent{\bf Corollary}. For $a\in C^n$, $a\ne0$, and $1\le rj+1$. \noindent{\bf Theorem}. Let $A\in C^{n\times n}$. Then $\exists$ a unitary matrix $Q$ such that $QAQ^\ast=H$ is upper Hessenberg. $Q$ may be computed as a product of Householder reflections, and $Q$ is orthogonal if $A$ is real. \noindent{\bf Corollary}. If $A\in C^{n\times n}$ is Hermitian or real symmetric, then $\exists$ a unitary matrix $Q$ such that $QAQ^\ast=T$ is tridiagonal ($T_{ij}=0$ for $|i-j|>1$). \vfil\eject