Linear Systems and LU

Matrix Facts

$Ae_j = a_j$ is the $j$-th column, $e_i^T A$ the $i$-th row, $e_i^T A e_j = A_{ij}$. Column operations act from the right, e.g. $A[e_1 + 2e_3 \;\; e_2 \;\; e_3 \;\; e_4 \;\; e_5]$ adds twice the third column to the first in place, and row operations act from the left. $(AB)^T = B^T A^T$, $(AB)^{-1} = B^{-1}A^{-1}$.

Triangular Systems

Forward substitution for lower triangular $Lx = b$,

$$x_1 = \frac{b_1}{L_{11}}, \qquad x_2 = \frac{b_2 - L_{21}x_1}{L_{22}}, \qquad x_i = \frac{b_i - \sum_{j < i} L_{ij}x_j}{L_{ii}}$$

and back substitution for upper triangular $Ux = b$ from the last row up. A triangular matrix is singular iff some diagonal entry is 0. Cost $\sum_{i=1}^n (2i - 1) \sim n^2$ flops.

LU Factorization

Adding $a$ times row $j$ to row $i$ is left multiplication by the unit lower triangular elementary matrix

$$L_{ij} = I + a e_i e_j^T, \qquad L_{ij}^{-1} = I - a e_i e_j^T$$

Gaussian elimination (for $4 \times 4$) gives $U = L_{43}L_{42}L_{32}L_{41}L_{31}L_{21}A$, so $A = LU$ with $L = L_{21}^{-1}L_{31}^{-1}L_{41}^{-1}L_{32}^{-1} \cdots$ unit lower triangular, its entries the row multipliers. To solve $Ax = b$:

1. factor $LU = A$,
2. forward substitution $Lz = b$,
3. back substitution $Ux = z$.

Cost $f = O(g)$ if $\lim_{n \to \infty} f(n)/g(n) < \infty$, $f \sim g$ if the limit is 1. A flop is one of $+, -, \times, \div, \sqrt{\cdot}$. Using $\sum_{k=1}^n k^p \sim \frac{n^{p+1}}{p+1}$,

$$\sum_{j=1}^{n-1} (n - j)[2(n - j) + 3] = \sum_{k=1}^{n-1} k(2k + 3) \sim \frac{2}{3}n^3$$

for LU, against $\sim n^2$ per triangular solve. If the run time is $t = Cn^p$, then $\log t = p \log n + \log C$, so $p$ is the slope on a log-log plot.

Pivoting

Not every matrix has an LU factorization (a 0 pivot), and small pivots give large multipliers and subtractive cancellation. Partial pivoting swaps the row with the largest entry in absolute value into the pivot position at each step, giving

$$PA = LU, \qquad A = P^T LU, \qquad Ax = b \iff LUx = Pb$$

with all entries of $L$ bounded by 1 in absolute value. $P = e_{i_1}e_1^T + \cdots + e_{i_n}e_n^T$ is a permutation matrix and $P^{-1} = P^T$. The pivot search is $O(n - k)$ at step $k$, $O(n^2)$ in total.

Matrix Structure

• Strictly diagonally dominant $|A_{ii}| > \sum_{j \neq i} |A_{ij}|$ for all $i$: nonsingular, and LU exists without pivoting.
• Banded Upper bandwidth $b_u$ if $A_{ij} = 0$ for $j - i > b_u$, lower bandwidth $b_l$ if $A_{ij} = 0$ for $i - j > b_l$, total bandwidth $b_u + b_l + 1$ ($b_u = b_l = 1$ is tridiagonal). Without pivoting $L$ and $U$ keep the bandwidths, the zeros are predictable, and LU costs $O(b_u b_l n)$.
• Symmetric $A^T = A$. $L$ and $U$ are not transposes, but $A = LDL^T$ with $D$ diagonal, and using the symmetry it costs $\sim \frac{1}{3}n^3$ flops. In general it still needs (symmetric) pivoting.
• Symmetric positive definite $x^T A x = \sum_i \sum_j A_{ij}x_i x_j > 0$ for all $x \neq 0$, all eigenvalues positive.

For SPD $A = LDL^T$, taking $x = L^{-T}e_k$ gives $0 < x^T A x = e_k^T D e_k = D_{kk}$, so $D = (D^{1/2})^2$ with $D^{1/2}$ the elementwise square root, and

$$A = LD^{1/2}D^{1/2}L^T = R^T R, \qquad R = D^{1/2}L^T \text{ upper triangular}$$

(Cholesky). No pivoting is needed, $\sim \frac{1}{3}n^3$ flops.

Conditioning of Linear Systems

For vector and matrix norms, with induced norms, $\|A\|_\infty$ the max row sum and $\|A\|_1$ the max column sum. Perturb $Ax = b$ to $A(x + h) = b + d$. Then $\|h\| \le \|A^{-1}\|\|d\|$ and

$$\frac{\|h\| / \|x\|}{\|d\| / \|b\|} \le \|A^{-1}\| \cdot \|A\| = \kappa(A)$$

The same $\kappa(A)$ bounds perturbations of $A$, and $\kappa(A) \ge \|AA^{-1}\| = \|I\| = 1$. If $\kappa(A) = 10^t$, about $t$ digits are lost. The residual $r = b - A\tilde{x}$ is a backward error, and

$$\frac{\|x - \tilde{x}\|}{\|x\|} \le \kappa(A) \frac{\|r\|}{\|b\|}$$