Numerical Linear Algebra Background

• matrix structure and algorithm complexity

• solving linear equations with factored matrices

• LU, Cholesky, LDL^T factorization

• block elimination and the matrix inversion lemma

• solving underdetermined equations

Matrix structure and algorithm complexity

cost (execution time) of solving Ax = b with A ∈ R^n×n

• for general methods, grows as n³

• less if A is structured (banded, sparse, Toeplitz, . . . )

flop counts

• flop (floating-point operation): one addition, subtraction,
multiplication, or division of two floating-point numbers

• to estimate complexity of an algorithm: express number of flops as a
(polynomial) function of the problem dimensions, and simplify by
keeping only the leading terms

• not an accurate predictor of computation time on modern computers

• useful as a rough estimate of complexity

vector-vector operations (x, y ∈ Rⁿ)

• inner product x^T y: 2n -1 flops (or 2n if n is large)

• sum x + y, scalar multiplication x: n flops

matrix-vector product y = Ax with A ∈ R^m×n

• m(2n - 1) flops (or 2mn if n large)

• 2N if A is sparse with N nonzero elements

• 2p(n + m) if A is given as A = UV^T , U ∈ R^m×p, V ∈ R^n×p

matrix-matrix product C = AB with A ∈ R^m×n, B ∈ R^n×p

• mp(2n - 1) flops (or 2mnp if n large)

• less if A and/or B are sparse

• (1/2)m(m + 1)(2n - 1) ≈ m²n if m = p and C symmetric

Linear equations that are easy to solve

diagonal matrices

lower triangular



called forward substitution

upper triangular flops via backward substitution

orthogonal matrices: A^-1 = A^T

• 2n² flops to compute x = A^T b for general A

• less with structure, e.g., if A = I - 2uu^T with we can
compute x = A^T b = b - 2(u^T b)u in 4n flops

permutation matrices:

where is a permutation of (1, 2, . . . , n)

• interpretation:

• satisfies A^-1 = A^T , hence cost of solving Ax = b is 0 flops
example:

The factor-solve method for solving Ax = b

• factor A as a product of simple matrices (usually 2 or 3):

(A_i diagonal, upper or lower triangular, etc)
• compute by solving k ‘easy’ equations

cost of factorization step usually dominates cost of solve step

equations with multiple righthand sides

cost: one factorization plus m solves
 

LU factorization

every nonsingular matrix A can be factored as
A = PLU

with P a permutation matrix, L lower triangular, U upper triangular
cost: (2/3)n³ flops

cost: (2/3)n³ + 2n² ≈ (2/3)n³ for large n

sparse LU factorization

A = P₁LUP₂

• adding permutation matrix P₂ offers possibility of sparser L, U (hence,
cheaper factor and solve steps)

• P₁ and P₂ chosen (heuristically) to yield sparse L, U

• choice of P₁ and P2 depends on sparsity pattern and values of A

• cost is usually much less than (2/3)n³; exact value depends in a
complicated way on n, number of zeros in A, sparsity pattern

Cholesky factorization

every positive definite A can be factored as
A = LL^T

with L lower triangular
cost: (1/3)n³ flops

cost: (1/3)n³ + 2n² ≈ (1/3)n³ for large n

sparse Cholesky factorization

A = PLL^TP^T

• adding permutation matrix P offers possibility of sparser L

• P chosen (heuristically) to yield sparse L

• choice of P only depends on sparsity pattern of A (unlike sparse LU)

• cost is usually much less than (1/3)n^3; exact value depends in a
complicated way on n, number of zeros in A, sparsity pattern

LDL^T factorization

every nonsingular symmetric matrix A can be factored as

A = PLDL^TP^T

with P a permutation matrix, L lower triangular, D block diagonal with
1 × 1 or 2 × 2 diagonal blocks

cost: (1/3)n³

• cost of solving symmetric sets of linear equations by LDL^T factorization:
(1/3)n³ + 2n² ≈ (1/3)n³ for large n

• for sparse A, can choose P to yield sparse L; cost ≪ (1/3)n³

Equations with structured sub-blocks
(1)

• variables

• if A11 is nonsingular, can eliminate
to compute x₂, solve

dominant terms in flop count

• step 1:   (f is cost of factoring A₁₁; s is cost of solve step)

• step 2:   (cost dominated by product of A₂₁ and )
• step 3:

total:

examples

• general no gain over standard method
#flops

• block elimination is useful for structured
for example, diagonal

Structured matrix plus low rank term

(A + BC)x = b

• A ∈ R^n×n, B ∈ R^n×p, C ∈ R^p×n

• assume A has structure (Ax = b easy to solve)
first write as

now apply block elimination: solve

then solve Ax = b − By

this proves the matrix inversion lemma: if A and A + BC nonsingular,
(A + BC)^-1 = A^-1 - A^-1B(I + CA^-1B)^-1CA^-1

example: A diagonal, B,C dense

• method 1: form D = A + BC, then solve Dx = b
cost: (2/3)n³ + 2pn²

• method 2 (via matrix inversion lemma): solve
(I + CA−1B)y = CA^-1b, (2)
then compute x = A^-1b - A^-1By
total cost is dominated by (2): 2p²n + (2/3)p³ (i.e., linear in n)

Underdetermined linear equations

if A ∈ R^p×n with p < n, rank A = p,

• is (any) particular solution
• columns of span nullspace of A
• there exist several numerical methods for computing F
(QR factorization, rectangular LU factorization, . . . )