Rating: 7.8/10.
Textbook about numerical linear algebra cover many algorithms, such as solving linear systems and finding eigenvalues. This book is relatively advanced, as it rarely provides examples and focuses on detailed proofs of convergence properties of algorithms. It also has no examples of computing any algorithms or experimental results, only proofs, so it’s probably less accessible to casual readers and practitioners. I went through about half of the book (the first 5 chapters) which are on more common topics, and skipped chapters 6-9 that cover more advanced material.
Chapter 1. Two basic examples of needing to solve systems of linear equations are: interpolation problem of finding a linear combination of basis functions to approximate a set of points, and solving a differential equation where the boundary values are given. The next few sections provide a review of Big-O analysis of algorithm complexity, properties of norms, inner products, symmetrical and orthogonal matrices, eigenvectors, properties of SVD and the pseudo-inverse.
Chapter 2: Discusses error, stability, and conditioning. The IEEE 754 double-precision floating-point format is a b-adic floating-point with a base of 2, 52 bits of precision, 11 bits of exponent, and one sign bit. The machine epsilon is 2^-51, and it guarantees that operations, including arithmetic operations, sin, log, etc, are within a factor of 1+eps of the true value.
Vector norms have several choices, such as the 1-norm, 2-norm, and inf-norm. They are topologically equivalent in that they are bounded within a constant factor of each other. Matrix norms are defined by how they affect vector norms after multiplication, and an induced matrix norm is the supremum of the effect it can have on any vector norm. Given a vector, the condition number of a matrix is |A| |A^-1|, which is always at least 1. If it is an induced matrix norm, then the norm of the identity matrix is 1.
To compute the matrix norm, you have to show that it is always less than some number C after multiplying by any vector, and there is some vector that achieves the value C. The matrix norm can be calculated for the 1-norm, 2-norm, and inf-norm. The 1-norm and inf-norm are the maximum sum of any row or column, respectively, while the 2-norm is the square root of the largest eigenvalue of A^T A, also the largest singular value of A for non-square matrices. The condition number can be defined as the ratio of the largest to the smallest non-zero singular value. It is possible to define matrix norms with different vector norms p and q in the numerator and denominator. In these cases, the matrix norms can be different if defined over the reals and the complex numbers: the norm over the reals is always less than the norm over complexes.
A condition number of a function measures whether a small change in x can lead to a large change in f(x); if the condition number is large, then it is said to be ill-conditioned. This can be calculated using limits or gradient information. Some examples: multiplication and square root are always well-conditioned, but addition can sometimes be ill-conditioned if the inputs are large and opposite; the formula for calculating the root of a quadratic polynomial is also ill-conditioned if using the standard formula, but it is possible to use a modified formula that is well-conditioned. The condition number of a matrix also measures how errors will be amplified when performing matrix operations like solving Ax=b. When analyzing algorithms, the condition number is often hard to compute so instead, we use a different definition of backward stability, which means for any input, there is a nearby input that produces the correct result.
Chapter 3: Direct methods of solving linear systems. When solving Ax = b: if A is triangular, then this is straightforward with back substitution; otherwise, we use Gaussian elimination. This is a sequence of row operations to convert Ax = b into an upper triangular form, does not require additional space and, assuming it does not end prematurely, it takes O(n^2) time.
LU factorization is A = LU, where L is lower triangular with diagonals 1 and U is upper triangular. If every principal submatrix has nonzero determinant, then the matrix has an LU factorization. It basically captures Gaussian elimination where L is all the row reductions and U is the resulting matrix. LU factorization is unique if it exists, and Gaussian elimination directly translates into an algorithm for LU factorization. It is also possible to compute L and U component-wise using some formulas. The LU factorization can be computed in O(n^3) time, but then once it is computed, it can be used to solve Ax equals b in O(n^2) time by two steps of forward and backward substitution. Tridiagonal matrices have a special formula of LU factorization where L and U are also nonzero only near the diagonal.
Pivoting will be required if the diagonal entries are zero, and even if nonzero, Gaussian elimination leads to numerical instability, so it is common to swap the row with a row below it that has the highest magnitude diagonal value. Pivoting is represented theoretically by permutation matrices, so every invertible matrix has LU factorization that is represented as PA = LU, where P is a permutation matrix. In an implementation, though, P is never actually constructed; we just exchange the rows directly. If a matrix is symmetric and positive definite, then we can turn LU decomposition into L L^T or Cholesky decomposition by decomposing U into a diagonal matrix and normalized upper triangular matrix, and this can be computed faster than standard Gaussian elimination.
QR factorization transforms a matrix, which can be non-square, into an orthogonal Q and upper triangular R. This is useful because LU decomposition does not preserve the condition number, but QR transforms do. We construct the Q from a series of Householder transform matrices, where each one represents geometrically a reflection that is also orthogonal. Multiplying a sequence of Householder matrices gives us an orthogonal Q matrix. This takes O(mn^2) time, does not require pivoting, and is unique up to signs of the diagonal entries.
When Ax = b is overdetermined with m > n, we usually solve least squares by minimizing the |Ax – b|, and the pseudo-inverse of A gives the optimal solution for this. Solving least squares by inverting A^T A is ill-conditioned, and the alternative is using QR factorization of A or using the SVD-based pseudo-inverse. There is several pages of perturbation analysis of this method of solving least squares using the SVD-based pseudo-inverse, when the condition number is analyzed, we realize that it is problematic if the matrix is rank-deficient, ie, when the matrix has singular values that are very close to zero, the condition number can blow up, so Tikhonov regularization (basically L2 regularization) adds a norm penalty when minimizing least squares so that the condition number doesn’t explode for rank-deficient matrices.
Chapter 4: Iterative methods of solving Ax = b: by evaluating a linear map F(x) that has the correct fixed point as the solution and is a contracting mapping (meaning it decreases the norm when applied), this will always converge to the fixed point. The condition for the iteration to always converge is if the spectral radius of C, where F(x) = Cx + c, is less than one, and then it doesn’t matter what c is.
Jacobi iteration uses the diagonal entries to obtain a fast and approximate inverse, and we iterate until the difference between the previous and current x is below a certain threshold, taking O(n^2) per iteration. The Jacobi iteration will converge if the matrix is irreducible (meaning it’s not made up of separate block matrices) and is diagonally dominant, meaning that the diagonal entries are bigger than the sum of all other entries for all rows.
The Gauss-Seidel iteration is a slight modification of the Jacobi method to use the current step entries of the vector instead of the previous step, and this converges faster; if the matrix is symmetric and positive definite, then Gauss-Seidel will converge. In Jacobi iteration, relaxation is choosing the new x to be a linear interpolation of the old x and the new iteration, and this sometimes converges faster. The optimal value for the interpolation parameters can be determined from the eigenvalues.
A similar relaxation method may be applied to Gauss-Seidel iteration, and this is called the Successive Over-Relaxation (SOR) method. As long as the Jacobi method converges, SOR always converges faster. For symmetric matrices, there is a symmetric variant of the Gauss-Seidel iteration called SSOR that maintains the matrix being symmetric at every step.
Chapter 5: First some preliminary theorems about eigenvalues. Gershgorin’s circle theorem localizes eigenvalues to a union of disks that depend on the matrix entries. The Rayleigh and Courant-Fischer theorems characterize eigenvalues of symmetric matrices as minimum and maximum ratios. Sylvester’s theorem states that any similarity transformation preserves the number of positive, zero, and negative eigenvalues. The Bauer-Fike theorem gives bounds on eigenvalues after perturbation of the matrix, and this depends on the condition number of the transformation to a diagonal matrix.
The power method (or von Mises iteration), finds the largest eigenvalue and corresponding eigenvector through a process that starts with a vector and iteratively multiplies it by the matrix while normalizing it, which eventually converges to the largest eigenvector from which you can calculate the corresponding eigenvalue. You must be careful to ensure the sign of the eigenvector doesn’t flip for it to converge, and the convergence speed is determined by the ratio of the largest eigenvalue to the second largest one.
The inverse iteration method (or von Wielandt iteration), performs the same method by applying it to the inverse of the matrix, which finds the smallest eigenvalue and eigenvector. A common use is to apply it to the inverse of a shifted matrix that is shifted by an initial guess of the eigenvalue; applying the inverse iteration provides a refinement of the eigenvalue and also the corresponding eigenvector, which is a method of getting an eigenvector when the approximate eigenvalue is known. Since this involves repeatedly solving a linear system, it can be pre-computed with LU factorization.
The Rayleigh quotient is a function that, if x is close to the eigenvector, then r(x) is close to the eigenvalue. Rayleigh iteration is like using the inverse power iteration but using the Rayleigh quotient to approximate the eigenvalue since it is not known, and if it converges, then it will converge in cubic time. Jacobi’s method is for computing all eigenvectors and eigenvalues of a symmetric matrix; it works by applying Givens rotations, which operate on two dimensions at a time, to zero out all of the diagonal entries until the matrix is diagonal. After convergence, which is guaranteed, the diagonals are the eigenvalues and the transformation matrix gives the eigenvectors.
The Householder reduction can also be applied similar to QR factorization, but producing QAQ^T, in the non-symmetric case it is not possible to transform to an upper triangular matrix; instead, we must reduce to a Hessenberg form where the band below the diagonal is nonzero and otherwise is upper triangular. The QR method is to compute A = QR, and then swapping the RQ gives Q^T AQ, which has the same eigenvalues, and it can be proven that when this is done iteratively, it converges to an upper triangular matrix, and the diagonal converges to the eigenvalues.
Naively applying the QR iteration process requires doing a QR factorization on each step, taking O(n^3) steps, but if we first reduce to Hessenberg form (which takes O(n^3) time), then computing a QR factorization on the Hessenberg matrix only takes O(n^2) time. Once the eigenvalue is known, inverse iteration can compute the associated eigenvector, which can be further improved using the shift method to reduce the size of the matrix by one in each iteration.
The computation of the SVD: naively you can apply eigenvalue methods on A^T A but this is not recommended to compute A^T A since it will square the condition number. Instead, we can first apply the Householder transform to convert the matrix into a bidiagonal matrix, which preserves singular values. If the bidiagonal matrix is reducible, we can split it into independent subproblems. Assuming it is irreducible (all diagonals and superdiagonals are nonzero), the Golub-Kahan SVD method then finds the eigenvalues of a diagonal matrix upper bidiagonal matrix using the shifted QR method.
Chapter 6. A simple iterative method for solving Ax = b by gradient descent – if the matrix is symmetric positive definite, then you can differentiate to obtain a function to minimize, and then each step changes the direction of steepest gradient and finds the step size to minimize along that direction. Convergence bounds for this method can be derived using Kantorovich inequality and depend on the condition number of the matrix, which can be slow if the condition number is large. In the naive descent algorithm, each direction of descent will be orthogonal to the previous one, but this can zigzag many times inefficiently if the matrix is ill-conditioned.
Instead, the conjugate gradient method defines two directions as A-conjugate if P^T AQ = 0. If a set of such directions are all mutually A-conjugate, then the descent algorithm always terminates after n steps for an n×n matrix. It’s possible to choose a new search direction that is A-conjugate to all previous directions from a formula given the previous direction and the residue. Krylov space is a sequence of subspaces generated by iterating a matrix from multiplying from a starting vector. The conjugate gradient method at each step produces the best approximation to the solution from successive Krylov spaces, and the number of iterations to reach a given position depends on the square root of the condition number rather than the condition number for the naive steepest descent method, so it converges much faster.
