Numerical linear algebra

Numerical linear algebra develops stable finite algorithms for matrix computations.

Matrix preconditioning transforms a system using an inexpensive auxiliary solve to improve its spectrum. For symmetric positive definite matrices, a positive definite $M$ produces $M^{-1/2}A{M}^{-1/2}$, preserving the symmetric structure needed by Conjugate gradient method.

A matrix preconditioner is an auxiliary invertible matrix whose inverse is cheap to apply and whose use improves the spectrum of a linear system. In symmetric positive definite problems, a positive definite $M$ preserves that structure under the symmetric transformation $M^{-1/2}A{M}^{-1/2}$.

Gaussian elimination applies elementary row operations to reduce a linear system to triangular form. Dense elimination uses $O(n^3)$ arithmetic operations.

An LU decomposition factors a matrix as a lower-triangular matrix times an upper-triangular matrix. It records Gaussian elimination so that multiple right-hand sides can be solved by forward and backward substitution.

An LDL decomposition of a real symmetric matrix iswhere $L$ is unit lower triangular and $D$ is diagonal. Symmetric elimination computes it one Schur complement at a time. The matrix is positive definite exactly when every diagonal pivot in $D$ is positive.

Machine precision is the characteristic relative rounding scale of a floating-point system.

For an approximate eigenpair $(\tilde{\lambda },\tilde{v})$, the residual is

If $\Vert \tilde{v}{\Vert }_2=1$, the smallest operator-norm perturbation that makes $(\tilde{\lambda },\tilde{v})$ an exact eigenpair has norm $\Vert r{\Vert }_2$. One such rank-one perturbation is

Successive Householder reflections zero subdiagonal column entries and factor a matrix into an orthogonal factor and an upper-triangular factor.

For full column rank $A\in {\mathbb{R}}^{m\times n}$, a full QR factorization gives $Q^{T}A=(R_1,0{)}^T$ with invertible upper triangular $R_1$. The equation $Ax=b$ is consistent exactly when the final $m-n$ entries of $Q^{T}b$ vanish. If consistent, back substitution gives its unique solution. This separates consistency from uniqueness in overdetermined systems.

Orthogonal Householder similarities reduce a real symmetric matrix to symmetric tridiagonal form.

The unshifted QR algorithm factors $A_k=Q_k{R}_k$ and forms $A_{k+1}=R_k{Q}_k=Q_k^T{A}_k{Q}_k$, preserving eigenvalues while driving the matrix toward block diagonal form.

LetThenThus ${\overline{Q}}_k{\overline{R}}_k$ is a QR decomposition of the matrix power $A^{k+1}$. In particular, the first $r$ columns of ${\overline{Q}}_k$ span the same subspace as the first $r$ columns of $A^{k+1}$ whenever those columns are independent.

An unshifted QR step preserves the bandwidth of a real symmetric banded matrix. A Givens-rotation implementation exposes this as bulges created during triangularization and removed when the factors are multiplied in reverse order.

The accumulated orthogonal factor in QR iteration is the $Q$ factor of a matrix power, so its leading columns span the same spaces as simultaneous power iteration.

Let a real symmetric matrix have an orthonormal eigenbasis $w_1,\dots ,w_n$, withWrite two starting vectors asFor two-column subspace iteration to recover
$\mathrm{span}(w_{n-1},w_n)$, the two projections onto that dominant subspace must be independent. The exact condition isRequiring each of the four coefficients to be nonzero does not imply this determinant condition.

When an off-block subdiagonal entry tends to zero, QR iteration asymptotically separates the corresponding invariant spectral blocks even if eigenvalues within one block have equal modulus.

A stationary iteration has the formwith a fixed iteration matrix $H$. It converges for every initial vector exactly when the spectral radius satisfies $\rho (H)<1$.

With $A=D+L+U$, the Gauss-Seidel method solves $(D+L)x^{(k+1)}=b-U{x}^{(k)}$. It uses each newly computed coordinate immediately in the remaining triangular solve. It is successive over-relaxation with parameter one, and converges from every start when $A$ is symmetric positive definite.

For a splitting $A=D+L+U$, successive over-relaxation updatesIt reduces to the Gauss-Seidel method at $\omega =1$. If $A$ is symmetric positive definite, the Householder-John theorem applies to $M=D/\omega +L$, $N=(1/\omega -1)D-U$, since $M^T+N=(2/\omega -1)D$ is positive definite for $0<\omega <2$.

For a stationary iteration $x^{(k+1)}=H{x}^{(k)}+v$, the matrix $H$ propagates the error by $e^{(k+1)}=H{e}^{(k)}$. Convergence from every initial error is equivalent to $\rho (H)<1$.

A splitting $A=M-N$ with invertible $M$ produces

Let $A=M-N$ be Hermitian positive definite. If $M^{\ast }+N$ is also Hermitian positive definite, then $M$ is invertible andThus the stationary iteration associated with the splitting converges from every initial vector.

Writing $A=D+L+U$, where $D$ is diagonal and $L,U$ are strictly triangular, the Jacobi iteration is

With relaxation parameter $\omega$, weighted Jacobi updates byIts iteration matrix is $H_{\omega }=I-\omega D^{-1}A$.

For the five-point discretization of $-\Delta +c$ on the unit square, the weighted iteration eigenvalues are $1-\omega [4{\sin }^2(\pi kh/2)+4{\sin }^2(\pi lh/2)+c{h}^2]/(4+c{h}^2)$. A fixed parameter converges for every mesh when $0<\omega \le 1$, but $\omega =1$ poorly damps the highest frequencies.

For the centered second-difference matrix with diagonal $-2$ and off-diagonals $1$, the weighted-Jacobi eigenvalues areThe method converges for every grid size exactly when $0<\omega \le 1$.

For modes with $k/(n+1)\in [1/2,1]$ and large $n$, the worst attenuation factor isIt is minimized by $\omega =2/3$, for which ${\mu }_{\omega }=1/3$.

For a three-by-three matrix with diagonal entries one and every off-diagonal entry $\alpha$, the Jacobi method iteration matrix has eigenvalue $-2\alpha$ on the all-ones vector and eigenvalue $\alpha$ with multiplicity two on its orthogonal complement. Its spectral radius is $2|\alpha |$, so Jacobi iteration converges exactly when $|\alpha |<1/2$.

If $A=D+L+L^T$ is symmetric positive definite and tridiagonal, set $S=\mathrm{diag}(1,-1,1,-1,\dots )$. Thenis positive definite. Applying the Householder-John theorem to $M=D$ and $N=-(L+L^T)$ proves that the Jacobi iteration converges.