Math 660 - Fall 1999

Summary

29 August

§1.1

Matrices

• Representing linear transformations.
• Functions on an integer grid.

Matrix and vector operations

• Matrix addition <-> addition of functions
• Multiplying a matrix by a scalar <-> multiplying a function by a scalar
• Matrix multiplication <-> composition of linear transformations
• Hadamard (Schur) product <-> (pointwise) multiplication of functions
• Partitioning

Algorithms

• SAXPY
• GAXPY
• for loop and MATLAB notation; j=1:n means increment j in steps of size 1
• j=a:h:b <-> increment j in steps of size h, starting at j=a and ending at or just below b

HW

1. Problem 1.1.1. Solution

---

2 September

§1.2 & 1.3

Partitioning/block multiplication of matrices A, B

• Column partition of A must match row partition in B
• Can multiply blocks in the same way as ordinary elements, except that the order of blocks must be preserved - A then B.

MATLAB notation

• row j of A is A(j,:)
• column k of A is A(:,k)
• rows j through l from columns k through m of A is A(j:l,k:m)

Banded matrices

• p = lower bandwidth; q = upper bandwidth
• Classification scheme
  • Diagonal: p = q = 0
  • Upper trangular: p=0 q = # of columns - 1 (or less)
  • Upper bidiagonal; p = 0 q = 1
  • Tridiagonal; p = q = 1
  • Upper Hessenberg: p = 1 and q = # of columns - 1 (or less)
  • similar classification for ``lower'' diagonal, etc.

HW

1. Problem 1.2.5 Solution
2. Problem 1.2.10 Solution
3. Problem 1.3.4
4. MATLAB matrix exercise. Let A the matrix below. Find the submatrices A(2:4,2:5), A([3,2,4],1:5), and A([1,4],[6,3]).

   -2	1	5	1	2	3
   0	4	9	-1	-6	-4
   1	-3	7	0	1	-8
   5	1	-1	1	5	3

5. Let A be the matrix above, and let B be a matrix with 6 rows, partitioned 1 2 | 3 4 5 | 6 . Give two block partitions for A compatible with B's partitioning scheme. How many partitions are there for A that are compatible with B?

---

7 September

§2.1

Vector spaces

Definition - Vector space. A vector space is a set V together with two operations, + and × . If u & v are in V, then u + v is in V; if r is a scalar, then r×v is in V. The operations satisfy the following rules.

Addition		Scalar multiplication
u + (v + w) = (u + v) + w		a×(b×u) = (ab)×u
Identitiy: u + 0 = 0 + u = u		(a + b)×u = a× u + b×u
Inverse: u + (-u) = (-u) + u = 0		a×(u + v) = a×u + a×v
u + v = v + u		1×u = u

Subspaces

Definition - Subspace. A subset S of V is a subspace if, under + & × from V, S is a vector space in its own right.

Theorem. S is a subspace of V if and only if these hold:

1. 0 is in S.
2. S is closed under + .
3. S is closed under × .

Ways of getting subspaces.

• Span of a set
• Sum and direct sum of subspaces
• Intersection of subspaces

Linear independence and linear dependence

---

9 September

§2.1 and class notes

Bases

Definition - basis. A subset B = {v₁ ... v_n} of V is a basis for v if B spans V and is linearly independent. Equivalently, B is a basis if it is maximally linearly independent; that is, B is not a proper subset of some other linearly independent set. Unless we specifically state otherwise, we will assume that B is ordered.

Theorem. Every basis for V has the same number of vectors as every other basis. This common number is defined to be the dimension of V, dim(V).

Coordinates and bases. The properties defining a basis are exactly the ones needed to define ``good'' coordinates for a vector space. Given an order basis B = {v₁ ... v_n} and a vector v, we can uniquely write the vector as v = x₁v₁ +...+ x_nv_n, and thus represent it by the column vector [v]_B = [x₁, ..., x_n]^T

Theorem A. Let V be a vector space, and let dim(V) = n.

1. V is isomorphic to Rⁿ (or Cⁿ, for the complex case).
2. Let C = {v₁ ... v_k}, with k < n. If C is linearly independent, we may add vectors to it to turn it into a basis.
3. Let D = {v₁ ... v_j}, with j > n. If D spans V, we may remove vectors from it to turn it into a basis.
4. If S is a subspace of V, then dim(S) is at most n. If dim(S) = n, then S = V.

Linear transformations

Definition - linear transformation. Let V and W be vector spaces, and let L: V -> W. L is called a linear transformation if for all scalars a, b and vectors u, v, we have L[au+bv]=aL[u]+bL[v].

Subspaces associated with L.

1. The null space of L - null(L) = {v in V | L[v] = 0} - is a subspace of V. The dimension of the null space is called the nullity.
2. The range of L - ran(L) = {w in W | w = L[v] for some v in V} - is a subspace of W. The dimension of the range is called the rank of L.

Matrix for L. Let B = {v₁ ... v_n} and C = {w₁ ... w_m} be bases for V and W, respectively. The m×n matrix A that represents L has its jth column given by A( : , j) = [ L[v_j] ]_C. To say that A represents L means that the equation w = L[v] is equivalent to the matrix equation y = Ax, where y = [w]_C and x = [v]_B.

Theorem B. Let K: U -> V and L: V -> W be linear transformations. Also, let A and B be the matrices representing L and K, respectively. (Appropriate bases are assumed given.)

1. The composition L o K is a linear transformation, and AB represents L o K. (Matrix multiplication is thus compositon for functions).
2. If both transformations take V to W, then the matrix representing L+K is A+B.
3. If c is a scalar, then the matrix for cL is cA.
4. If L: V -> V is one-to-one and onto, then the inverse mapping L^-1 is linear. The matrix for L^-1 is A^-1, the matrix inverse of A.

HW

1. Prove Theorem A, Part 3.
2. Prove Theorem A, Part 4.
3. Show that ran(L) is a subspace of W.
4. Find the matrix A of L: P₃ -> P₃, where L[p]=p' and B = {1, x, x², x³}.
   1. Find a basis for the column space of A from among its columns. (See Theorem A, Part 3.) Find the rank of A (and of L, too.) Use the correspondence between A and L to read off ran(L) in terms of polynomials.
   2. Find the null space of A, and read off the null space of L in terms of polynomials. Find the nullity of A (and, again of L). Verify that rank(A) + nullity(A) = 4.
5. Show that if A is an m×n matrix, then rank(A) + nullity(A) = n, the number of columns.

---

14 September

§2.1, §2.2, and §2.3

Determinants

1. Recursive definition and permuation definition
2. If B is A with two rows interchanged, then det(B) = - det(A).
3. If B is A with a multiple of one row added to another, det(B) = det(A).
4. If A has two rows identical, then det(A) = 0.
5. det(A^T) = det(A)
6. If A, B are both n×n matrices, then det(AB) = det(A)det(B)
7. A is invertible if and only if det(A) is not 0.

Norms

Definition -- Norm. Let V be a vector space, and let f : V --> R. f is a norm for V if these hold.

1. positivity - f(x) is nonnegative, with f(x) = 0 ==> x=0.
2. triangle inequality - f(x+y) <= f(x) + f(y)
3. positive homogeneity - f(a x) = |a|f(x)

p-norms -- ||x||_p: = (|x₁|^p + ... + |x_n|^p)^1/p, p >= 1; ||x||_infinity := max{ |x_j| },

Holder's Inequality. |x^Ty| <= || x ||_p|| y ||_q, 1/p + 1/q = 1

Minkowski's Inequality -- || x + y ||_p <= || x ||_p + || y ||_p

Equivalence of norms. c₁ || x ||_a <= || x ||_b <= c₂|| x ||_a

Convergence. x^(k) converges to x if and only if lim_{k - > infinity}|| x^(k) - x || = 0.

Matrix Norms. Frobenius norm and operator norms.

HW (Due 9/23)

1. Prove Holder's inequality. Two hints:
   1. If 0< s < t, and if a,b are positive, a + b = 1, then ln(a s + b t) >= a ln(s) + b ln(t); consequently, s^a t^b <= a s + b t.
   2. There is no loss of generality in assuming that || x ||_p = || y ||_q = 1

---

16 September

§2.3

Matrix Norms

Definition - Operator Norm. Let V and W be vector spaces having norms || v ||_V and || w ||_W, respectively. If L:V -> W, then we define the operator norm of L by this:
|| L ||_V,W := sup{ || L[v] ||_W ||v||_V^-1 } (v is not 0.)
There are two other equivalent ways to formulate this.

1. || L ||_V,W := sup_{|| v || = 1}{ || L[v] ||_W }
2. || L ||_V,W is the smallest C for which || L[v] ||_W <= C||v||_V

Matrix Operator Norms We apply the definition above to the case in which we have

• A is in R^m×n
• V = Rⁿ and || v ||_V = || v ||_r (p-norm, p=r)
• W = R^m and || w ||_W = || w ||_p (p-norm, p=p)

and we set || A ||_r,p = operator norm || A ||_V,W. Also, if r=p, || A ||_p := || A ||_p,p

Properties

• There is a vector u in Rⁿ, with || u ||_r = 1, such that || A ||_r,p = || A u ||_p
• || A ||₁ = || A ||_1,1 = max_j=1..n || A( : , j ) ||₁
• || A ||₂ is the square root of the largest eigenvalue of A^TA
• Submultiplicative property. || A B ||_p <= || A ||_p|| B ||_p (Holds in other situations.)

Application - Perturbations

• If F is in R^n×n, || F ||_p < 1, then I - F is nonsingular, and its inverse is given by the convergent Neumann series (geometric series)
  (I - F)^-1 = I + F + F² + F³ + ...
  Moreover, || (I - F)^-1 ||_p <= (1 -|| F ||_p )^-1

HW (Due 9/23)

1. 2.2.1
2. 2.2.5
3. 2.2.8
4. 2.3.3
5. 2.5.1 Solution

---

21 September

§2.4, §2.5.1

Floating Point Numbers

• Base = b, Precision = t, Lower bound on exponent = L, Upper bound = U
• Form of number: ±.d₁d₂...d_t × b^e, where e is any integer, staring at L and running up to U. the d_j's are integers starting at 0 and running up to b-1. d₁ is never 0.
• A floating point system is finite.
• Catastrophic loss of precision is possible.

Orthogonal Sets

``Dot'' product: < x, y > = y^T x (for the real case; < x, y > = y^H x in the complex case)

Angle between x and y: cos^-1(< x, y >/|| x || || y ||)

x and y are orthogonal if < x, y > = 0

An orthogonal set of vectors is a set for which each pair of vectors is orthogonal.

A set is orthonormal if it is orthogonal and if each vector in the set has norm 1.

The orthogonal complement of a set S is the space of all vectors orthogonal to every vector in S.

Orthogonal matrix. An n×n matrix U is orthogonal if U^TU = I; That is, U^T = U^-1.

Proposition A.

1. An orthogonal set of nonzero vectors is linearly independent.
2. An orthonormal set is linealy independent.

Proposition B.

1. An orthonormal set having n vectors is a basis for Rⁿ.
2. A set {u₁, ..., u_n} is orthonormal in Rⁿ if and only if the matrix U = [u₁ ... u_n] is orthogonal.
3. If the n×n matrix U is orthogonal, so its columns are a basis, then the formula for changing coordinates relative to the standard basis (SB) to the new, orthonormal basis, is
   [x]_NB = U^T[x]_SB

HW (Due 9/30)

1. Find all floating point numbers in the case that b = 3, t = 3, L=-1, and U = 2. Give an example illustrating ``catastrophic cancellation'' for this system.
2. Prove formulas 2.5.1 and 2.5.2 in the text. Solution
3. Show that the singular values of an m×n matirx A are square roots of the eigenvalues of A^TA.
4. 2.5.8
5. 2.5.9
6. Find the singular value decomposition for the matrices given below.

   1	-1	1
   3	-1	2

   1	3
   3	-1
   4	2

---

23 September

§§2.5.1-2.5.3

Gram-Schmidt process

Singular Value Decomposition

HW See the HW for 21 September.

---

28 September

§§2.5.4-2.5.5, 2.7.1

Thin SVD

Let A be an m×n matrix with SVD given by A = U^T S V. If m >= n, then A = U(1:n,:)^T S(1:n,1:n) V.

Numerical Rank

Definition.   rank(A,t) := min _{||A - B|| <= t} rank(B)

Distances to Lower Rank Matrices

Theorem.   min _{rank(B) = k} ||A - B||₂ = s_k+1 ,   k < rank(A)

Corollary.   If r_t = rank(A,t), then s_rt+1 <= t <= s_rt

Inverse of a square matrix via SVD expansion

A = s₁ u₁ v₁^T + ... + s_n u_n v_n^T   ==>   A^-1 = s₁^-1 v₁ u₁^T + ... + s_n^-1 v_n u_n^T

HW (Due 10/7)

1. This exercise requires you to use MATLAB or software capable of doing singular value decompositions. (You may, of course, write your own C or FORTRAN programs.) Let A be the matrix below. Plot rank(A,t) vs. log₁₀ t for t = 10^-7 to t = 10. Generate 3 random 6×6 matrices with coefficients < 0.05. Add each matrix to A, and repeat the exercise with it.

   1	1/2	1/3	1/4	1/5	1/6
   1/2	1/3	1/4	1/5	1/6	1/7
   1/3	1/4	1/5	1/6	1/7	1/8
   1/4	1/5	1/6	1/7	1/8	1/9
   1/5	1/6	1/7	1/8	1/9	1/10
   1/6	1/7	1/8	1/9	1/10	1/11

2. For the four matrices generated in the previous exercise, find the condition number relative to the 2-norm.
3. Let b be the column vector below. First, solve A x = b, with A as in question 1. Then Generate 3 random 6×1 column vectors with coefficients < 0.05. Add each vector to b, and resolve the system with the new value of b. Calculate the relative change in x divided by the relative change in b. (Use 2-norms here.)

   -1

   2

   3

   4

   -1

   3

---

30 September

§§2.7.2-2.7.4

Condition number

Definition. K(A) = ||A|| ||A^-1||

In solving A x = b, the relative error in x behaves like K(A) times the sum of the relative error in A plus the relative error in b.

2-norm condition number. K₂(A) = s₁ (s_n)^-1

MATLAB demo - diary file

---

5 October

§§2.7.3 & 3.1-3.2.5.

Condition number - revisited

Reviewed definition.

Discussed why det(A) is a poor measure of the conditioning of A.

Upper and lower triangular matrices

Backward substitution. Solving systems of the form Ux=b, where U is upper triangular.

Forward substitution. Solving systems of the form Lx=b, where L is lower triangular.

Properties.

1. Products of upper traingular matrices are upper triangular.
2. Inverses of upper triangular are upper triangular.
3. These two properties hold for unit upper triangular matrices.
4. All the properties listed hold for lower triangular matrices, too.

Gauss elimination

Matrix form of row operations

• Row(j) = c*Row(j)   <==>   I - (1-c)e_j e_j^T
• Row(j) = Row(j) - c*Row(k), k ~ = j,   <==>   I - c*e_je_k^T
• Zero out all entries in a column x below x_k (assume x_k ~ = 0)   <==>   M_k = I - t_ke_k^T, where
  t_k(j) = 0 if j <= k and t_k(j) = x_j x_k^-1 if j > k
• M_k^-1 = I + t_ke_k^T

LU decomposition

A can be factored in the product of a unit lower traingular matrix L and an upper triangular matrix U, provided all pivots, u(1,1), u(2,2), ... , u(n-1,n-1), are all nonzero.

det(A) = u(1,1)u(2,2) · · · u(n,n)

HW (Due 10/19)

1. Let A be a real, 2×2 matrix having the singular value decomposition A = USV^T. Using the Frobenius norm, find the distance from A to the rank 1 matrices. (Hint: Use the invariance of the Frobenius norm under orthogonal transfromations.)

---

7 October

§§3.2.5-8, 3.2.12.

An example

Let A be the 3×3 matrix below.

1	-2	3
2	-5	4
-1	3	4

The Gauss transformations used to triangularize A all have the form M_k=I - t_k e_k^T, where
t₁= [0 2 -1]^T and t₂= [0 0 -1]^T.
As we showed in class,

L = M₁^-1 M₂^-1 = I + t₁ e₁^T + t₂ e₂^T =

1	0	0
2	1	0
-1	-1	1

and the upper triangular matrix U = M₂M₁A is
U =

1	-2	3
0	-1	-2
0	0	5

When does a matrix have an LU factorization?

Proposition. The pivot u_kk for an n×n matrix A is
u_kk = det(A(1: k,1: k)) / det(A(1: (k-1),1: (k-1))) .

Theorem. An n×n matrix A has an LU factorization if and only if all submatrices A(1: k,1: k),
k=1: (n-1), are invertible.

Outer product algorithm

HW (continued from 5 October)

2. Using the matrix A in the example, go through each step in the algorithm on p. 98 of the text.
3. Determine which of these 3×3 matrices has an LU factorization. For each that does, find it.

   4	-1	1
   1	4	6
   -2	3	-7

   2	-2	1
   -5	5	-9
   -4	3	1

   2	-2	-1
   -1	-1	3
   -4	-6	2

4. For each of the matrices above, use MATLAB's lu command to find an appropriate factorization. In each case, look at L*U given by the command. Is it the same as the A you started with?
5. 3.2.5

---

12 and 14 October

Review and mid-term test.

---

19 October

§§3.4.1-3.4.4

Partial Pivoting

Permutation matrices and properties

• P is a matrix with rows or columns of I permuted.
• P is orthogonal.
• P*A permutes the rows of A
• A*P permutes the columns of A.

Interchange matrices

• E_{j <-> k} interchanges rows j and k.
• E_{j <-> k} M_k = M'_kE_{j <-> k}, where M and M' are Gauss transformations

Strategy - pivot off the largest entry on or below the diagonal in a column.

Complete pivoting selects the largest entry from the whole submatrix remaining, not just the part of the column on or below the diagonal.

---

21 October

§§4.1.1-4.2.1

LDM^T Decompositon

Theorem.  If an n×n matrix A and all its principal submatrices are invertible, then there are unit lower triangular matrices L and M, and an invertible diagonal matrix D such that A = LDM^T. Moreover, the decomposition is unique.

Proof: This is the LU decomposition in disguise. Here, D = diag(pivots) and M^T = D^-1U. We proved uniqueness in class. Also, there is a version in case pivoting is required.

Algorithm. Derived in class using partitioning.

LDL^T Decompositon

Theorem.  If A is an invertible, symmetric n×n matrix, and all its principal submatrices are invertible, then A = LDL^T, where D is an invertible diagonal matrix and L is a unit lower triangular matrix. Moreover, the decomposition is unique.

Proof: A = A^T implies that LDM^T = MDL^T. Uniqueness of the LDM^T factorization implies that M = L.

For symmetric A, the algorithm mentioned above takes roughly half the number of flops as the usual Gaussian elimination algorithm.

Positive Definiteness

Definition. A quadratic forms Q(x) = x^T A x is positive definite if Q(x)>0, unless x = 0. (The matrix A is not unique. The same quadratic form can be generated by A + B, where B is antisymmetric.)

Definition. If A is a symmetric matrix, then A is said to be positive definite if the associated quadratic form Q(x) = x^T A x is positive definite. Here are a few properties; throughout, we assume A is positive definite.

• A is invertible.
• If X is a k×n matrix with rank(X) = k, then X^T A X is a positive definite k×k matrix.
• Every principal submatrix of A is positive definite.

HW (Due 10/28.)

1. Find the permutation matrix P corresponding to p = [2 3 5 4 1]. Express this matrix as a product of interchange matrices. Determine whether the permutation is even or odd.
2. Problem 3.4.1.
3. Problem 3.4.3. (Write an m-file or C or FORTRAN program to do this.)
4. Write MATLAB m-files (or C or fORTRAN) programs to carry out the steps in Algorithm 3.4.1 on pg. 112, modified by the change mentioned on the top of pg. 114. Test this m-file (or program) on the matrix in Example 3.4.2.
5. Show that the matrices A and A+B generate the same quadratic form if and only if B is antisymmetric.

---

26 October

§§4.2.3-4.2.5

Cholesky Factorization

• A=A^T and positive definite ==> A = GG^T, where G is a lower triangular matrix with positive diagonal elements.
• G is called the Cholesky triangle
• The factorization follows from A = LDL^T; in fact, one has that G = LD^½

Algorithms for finding G

Update columns - §4.2.4

• G(j,j)*G(:,j) = A(:,j) - G(j,1)*G(:,1) - G(j,2)*G(:,2) - ... - G(j,j-1)*G(:,j-1)
• G(j,j) = (A(j,j) - G(j,1)*G(j,1) - G(j,2)*G(j,2) - ... - G(j,j-1)*G(j,j-1))^½

Outer product rank-1 update - §4.2.5

Note: Both algorithms provide a test for a symmetric matrix to be positive definite.

---

28 October

§§4.2.8-4.2.10, 4.3.1-4.3.2.

Symmetric positive semidefinite (spsd) matrices

Definition. A=A^T and x^T A x >= 0 for all x in Rⁿ

• |A(j,k)|<=(A(j,j)*A(k,k))^½
• A(j,j) = 0 ==> A(j,:) = 0 and A(:,j) = 0
• Algorithms for getting the Cholesky factorization in the positive definite case work for spsd matrices, with minor modifications in them.
• Small diagonal elements make pivoting useful.

Symmetric pivoting - PAP^T

Square root and polar decomposition

Banded matrices and the LU factorization

Theorem. Let A = LU. If A has upper bandwidth q and lower bandwidth p, then the lower triangular matrix L has lower bandwidth p, and the upper triangular matrix U has upper bandwidth q.

Forward substitution to solve Lx = b requires ~2pn flops, rather than n². Similarly, backward substitution in Ux=b requires only ~2qn flops.

HW (Due 11/9.)

1. Consider the symmetric 4×4 matrix A below. Use the Cholesky Algorithm 4.2.1 (p. 144) to find G or to show that it is not positive definite.

   A =

   12	6	4	3
   6	4	3	2
   4	3	3	2
   3	2	2	3

   Write a MATLAB (or C or FORTRAN) program to implement Algorithm 4.2.2 (p. 145).

   Problem 4.2.1(a)

   Problem 4.2.5 (Assume A=A^T.)

   Problem 4.2.11

---

2 November

§§4.3.1-4.3.3, 4.3.5, 5.2.7

Banded matrices - Cholesky factorization (finished).

Standard (``classical'') Gram-Schmidt algorithm for ``thin'' QR factorization

---

4 November

§§5.2.6-5.2.8

Modified Gram-Schmidt algorithm

QR Factorization - properties

Thin QR factorization is unique, and G = R^T in the Cholesky factorization of A^TA

ran(A(1:k)) = ran(Q(:,1:k))

null(A^T) = ran(Q(:, (n+1):m))

Connection with least squares data fitting

HW (Due 11/17.)

1. Problem 4.3.1
2. Let f(x) be a smooth function that has period 1. That is, f(x+1) = f(x).
   1. Use a Taylor series argument to show that f ''(x) = (f(x+h) + f(x-h) -2f(x)) h^-2 + O(h), for h small.
   2. Let h = 1/n. Suppose we are given y(k+1) = f(kh), k=0:(n-1). Show that

      f ''(jh) = -h^-2 A(j+1,:)*y + O(h),  j = 0:(n-1),

      where A is the matrix given in Problem 4.3.11.

3. Problem 4.3.11(a,b)
4. Problem 5.2.6 (Note: stated formula is incorrect, unless v is a column vector, which requires k=n-1. You may make this assumption.)
5. Problem 5.2.7

---

9 November

§§5.1.1-5.1.6

Householder Matrices

• Q = I - c*v*v^T, c = (2/v^Tv)
• Q is orthogonal.
• If v = x - || x ||₂ e₁ , then Q*x = || x ||₂ e₁
• Q represents a reflection about an axis perpendicular to v.
• Householder matrices play the same role in QR factorization that Gauss transformations do LU factorizations.
• Algorithm 5.1.1 computes v and c to machine accuracy; v is normalized so that v(1)=1.

Working with Householder matrices

• If A is an m×n matrix and we want the product Q*A, where Q is a Householder matrix, then use the ``update'' form of multiplication: Q*A = A - b*v*w^T, where w = A^T*v.
• Direct multiplication requires many more FLOPS.

---

11 November

§§5.2.1, 5.3.1-5.3.5

Householder QR algorithm

Least squares problems

• Arise in data-fitting problems
• Normal equations, A^TA*x - A^T*b = 0, are derived by minimizing f(x) = ½|| A*x - b ||² over all x in Rⁿ. Solve via Cholesky approach:
  1. Factor: A^TA = G*G^T, where G is the Cholesky triangle.
  2. Solve via forward substitution: G*y = A^T*b
  3. Solve via back substitution: G^T*x = y
• QR approach: A = Q*R, f(x) = ½|| R*x - Q^Tb ||². Solve R(1:n,:)*x = Q^Tb(1:n) via back substitution.
  • Householder QR - form Q^Tb using factored form of Q.
  • Modified Gram-Schmidt - Use on the augmented matrix [A b].

HW (Due 11/24.)

1. P5.1.2
2. P5.1.4
3. P5.2.2
4. P5.2.5
5. P5.3.5

---

17 November

§§5.3.5 (finished), 5.3.7, 5.4.1-5.4.3

Corrrected Problem 5.2.6 and gave hints on 5.2.7.

Sensitivity of LS - worked Problem 5.3.2

Rank deficiency and qr with column pivoting - an example

Complete orthogonal decompositions

Bidiagonalization

---

19 November

§§5.5

Pseudoinverse

• Any matrix A maps ran(A^T) one-to-one and onto ran(A). The inverse map A⁺ is the pseudoinverse.
• If the svd for A is A = USV^T, then the pseudoinverse is given by A⁺ = VS⁺U^T, where
  S⁺ = diag(s₁^-1,...,s_r^-1)
• Moore-Penrose conditions

Rank deficient LS problems

• Problem: Minimize || A x - b ||, where A is an m×n matrix, and rank(A) is less than both m and n.
• The minimizer x is not unique. However, x = A⁺ b solves the problem, and, among all minimizers, has the smallest 2-norm.

---

23 November

§10.2.1

The method of steepest descent

HW (Due 12/2.)

1. Write a MATLAB (or C or FORTRAN) program to upper bidiagonalize a matrix, and use it to upper bidiagonalize the 4×4 matrix B below.

   B =

   12	6	4	3
   6	4	3	2
   4	3	3	2
   3	2	2	3

   Consider the LS problem: minimize || A x - b ||, where b = [ 1 0 -1 2 0]^T and A is given below. Show that A is rank deficient. Find all x that solve the LS problem. (You may use MATLAB etc.)

   A =

   3	1	-4	-3
   -1	-4	5	1
   -1	0	1	1
   1	-2	1	-1
   1	3	-4	-1

   P5.5.2

   P10.2.1

   Let D = diag(1,1,10^-2), b = [ 1 1 0 ]^T, F(x) = || D x - b ||², and x₀ = [ 2 2 1 ]^T. Minimize F using the method of steepest descent with an initial guess x₀. Illustrate convergence by plotting iterates x_k. Use MATLAB to do the graphics and computations. Solution

---

30 November

§10.2.2-10.2.4, 7.1.1, 8.1.1

Conjugate gradient method

• If we have moved along directions p₁, p₂, ..., p_k-1, then we choose p_k to satisfy these conditions:
  1. p_k^T A p_j = 0, j=1:(k-1). (A-conjugacy)
  2. Let r_k-1 be the k-1 residual. p_k minimizes || p - r_k-1 || over all p satisfying the A-conjugacy conditions above.
  3. If the initial direction is p₁=r₀, then the updates are computed via
     
     p_k = r_k-1 + b_k p_k-1,   where b_k = (p_k^T A r_k-1)(p_k^T A p_k-1)^-1
• This avoids the slow convergence sometimes encountered in steepest descent. See the m-file conjgrad.m

---

2 December

§8.1.1

Schur decomposition for a symmetric matrix A.

• A Q = Q D, Q = [q₁ q₂ ... q_n] is orthogonal, and its columns are eigenvectors of A. D is diagonal and its entries are (in order) eigenvalues corresponding to the q_j's.
• Proof via Housholder transformations and induction.
• || A ||₂ is max_j=1:n { |d_j| }
• Connection with the svd.
  1. A = A^T and A = U S V^T ==> A² = U S² U^T
  2. Schur decomposition: A = Q D Q^T ==> A² = Q D² Q^T
  3. Result: D² = Q^TU S² U^TQ or, equivalently, D²W = W S², where W = Q^TU is orthogonal.
  4. Compare k^th columns: D²W(: , k) = s_k² W(: , k), so (D² - s_k² I) W(: , k) = 0. Hence, we have these equations:
     
     (d₁² - s_k²) W(1, k) = 0
     (d₂² - s_k²) W(2, k) = 0
     ...
     (d_n² - s_k²) W(n, k) = 0
     Since at least one of the W(j,k)'s is not 0, we have for that j, d_j² = s_k². Thus the singular values of A coincide with the magnitudes of the eigenvalues of A. One can also show that if all singular values are distinct, then W is a permutation matrix.

Courant-Fischer Theorem

---

7 December

Review