Numerical Linear Algebra February 7, 2008 1

— David Wagner 2008/02/14 11:14

Work any 7 of the following problems. You make work more for additional credit.

1. Suppose a₁ and a₂ are linearly independent vectors, L = span{a₁, a₂}, and b is a vector not in L. By minimizing the function φ(s, t) = ½||b − (sa₁ + ta₂)||², find a formula for the projection of b onto the subspace L. Answer: if A = [a₁ a₂], then P = A(A*A)^-1A*.

1. Assume given are real scalars or vectors in ℜⁿ;
   • A=[a₁ a₂]; Pb = L = A[s t]*;
   • A=; A*=. Note the row and column indices.
2. ( )
   • The minimum of this function occurs at the smallest distance between b and a linear combination of a₁ and a₂. This will occur when both first partials are zero and second partials are positive.
3. ==0
4. ==0
   1. >0✔
   2. >0✔
5. 
6. 
7. =
8. =
9. A*A[s t]* = A*b. Recall the row and column indices are as noted above.
10. (A*A)^-1A*b = [s t]*
11. A(A*A)^-1A*b = A[s t]* = Pb
12. → P = A(A*A)^-1A*✔

2. Suppose that A is an n × n matrix. Consider solving Ax = b, by minimizing the quadratic functional φ(x) = ½||Ax − b||². Show by taking partial derivatives that ∇φ(x) is just A∗(Ax − b). (Consider φ = φ(x₁, x₂, …, x_n), use the definition of Euclidean norm and Ax, then find ∂φ/∂x_k, k = 1, …, n.)

1. . Assume real numbers.
2. ;
3. = = =
4. =
   • The outermost sum (index i) simply adds the rows of the resultant column vector.
   • The minimum of this sum of nonnegative values occurs when each term is minimized.
   • Each term is minimal when partials x_k are all zero. (And second partials at this point are positive; see above.)
   • Only a few terms are retained in each partial.
   • But looking at this monster makes me cross-eyed. If accounting was my thing, I'd be certified already.
5. = ½(Ax-b)*(Ax-b)

3. Consider the matrix

By hand calculation find the orthogonal projection P onto the range(A) and the image under P of the vector (1, 2, 3)^T. (See problem 1.)

1. P = P² = range(A) = Ax; v=[1 2 3]^T. Find P, Pv.
2. A* =
3. P = A(A*A)^-1A* = = = =
   • P =
4. Pv =

4. Show that if {v_j} is a basis (not necessarily orthonormal) of a subspace L, then the closest point w to y in span({v_j}) is given in terms of the Gramian matrix, M_ij =<v_i, v_j>. (Use the fact that <y − w, v> = 0 for all v in L.) How is this related to the stiffness matrix used in finite element methods?

5. If P is an orthogonal projection, then I − 2P is unitary. Prove this algebraically and give a geometric interpretation.

• Projection is Idempotent: P = P² = PP
• Orthogonal: P = P*

Unitary

• (I - 2P)*(I - 2P) = (I - 2P)(I - 2P)* = I
• (I - 2P)^-1 = (I - 2P)*

1. (I - 2P)*(I - 2P)
2. = (I* - 2P*)(I - 2P)
3. = I*I - 2P*I - 2I*P +4P*P
4. = I -2P* -2P +4PP
5. = I -2P -2P +4P
6. = I✔

The geometric interpretation has something to do with the fundamental theorem of linear algebra and projections into the fundamental subspaces.

Needs diagram

6. Let A ∈ C^n×n be hermitian symmetric (A = A∗). (a) Prove that all eigenvalues of A are real, and (b) Prove that if x and y are eigenvectors corresponding to distinct eigenvalues, then x and y are orthogonal.

• Hermitian: A=A*.
• Symmetric: a_ij=a_ji. a_k (and a_k*) is the column (or row) corresponding to lambda_k. Note a_k = a_k*; whether it is a row or column is implied by context.
• Hermitian and symmetric implies a real matrix since .

1. i,j = 1,2,…,n
2. Ax = λx = A*x
3. Ax_k = λ_kx_k = A*x_k.
4. = = a_k*x_k = a_k x_k.
   • Since a_k is real, lambda_k must be real.
     1. (λ_k x_k)*x_k = (a_k x_k)*x_k. [Dot both sides.]
     2. λ_k (x_k*x_k) = x_k*(a_k* x_k) = a_k(x_k*x_k). [Do a little conjugate symmetry do-si-do.]
        • (x_k*x_k) is real, as is a_k ∴ lambda_k must also be real.

Now part b.

1. Ax = λ_xx = A*x; Ay = λ_yy = A*y; λ_x≠λ_y
2. Ax = λ_xx
3. (Ax)*y = (λ_xx)*y [Dot both sides again]
4. x*(A*y) = λ_x(x*y) [More dancing, this time a cha-cha.]
5. x*(λ_yy) = λ_x(x*y)
6. λ_y*(x*y) = λ_x(x*y)
   • The complex conjugate of a real number is itself: λ_y* = λ_y.
7. λ_y(x*y) = λ_x(x*y)
   • Since λ_x and λ_y are distinct and (assumed) nonzero, (x*y) must be zero.
8. ∴ x is orthogonal to y.

7. What can be said about the eigenvalues of a unitary matrix U? (A unitary matrix is one whose adjoint is also its inverse, i.e., U⁻¹ = U∗, equivalently U*U=I).

1. Ux = λx
2. U⁻¹Ux = U⁻¹λx = U*x
3. Ix = U*λx
4. x_k=U*λ_kx_k
5. = = =
6. =
7. 

8. Let A be the operator defined by Au = −u'' acting on the space C²₀[0, 1] of functions satisfying zero Dirichlet boundary conditions. If K is a positive integer we showed in class that λ_k = k²π² is an eigenvalue with associated eigenfunction sin(kπx). Fix n and consider the corresponding discrete problem on a finite difference grid of n − 1 interior points. Show that the eigenvalues of the matrix −n²A, where A is the matrix with ones on the sub- and super-diagonal, and -2 along the diagonal, are given by λ_k = 2(1 − cos(πk/n )) for k = 1, …, n − 1 with associated eigenvector (u_k)_j = sin(πkj/n) for j = 1, …, n − 1. What happens in the limit as n → ∞? (Hint: consider the analogy with the continuous problem above, now finding an eigenvector-eigenvalue pair amounts to solving a difference equation with boundary conditions. The role of the exponential functions is now played by vectors of the form u_j = ξ^j, for instance. Challenging, but it can be done by following the ODE boundary value problem blueprint. If your ODE skills are a little rusty, you could just use the definition of eigenvalue-eigenvector pair and verify the assertion. Note that this approach will require adroit use of various trig identities.)

9. Let A be a symmetric matrix of order n, with eigenvalues λ₁ ≤ λ₂ ≤ · · · ≤ λ_n and define the Rayleigh quotient R(x) = <Ax, x>/<x, x>, x ≠ 0. Show that max_x R(x) = λ_n and min_x R(x) = λ₁. (Show that optimizing ˜R(x) = <Ax, x>, for ||x|| = 1 gives an equivalent problem, and use the fact that A has a complete set of orthonormal eigenvectors to expand an arbitrary unit vector in terms of this basis.)

1. A ∈ ℜ^nxn; i ∈ [1,2,3,…,n].
2. Ax = λx = (Ax)*x = (λx)*x = λ(x*x). [Dot both sides by x, then apply scalar linearity.]
3. λ = (Ax)*x/(x*x) = <Ax, x>/<x, x> = R(x). [x ≠ 0 ⇒ (x*x) = <x, x> ≠ 0.]
   • Clearly, the maximum R(x) is the maximum λ, λ_n, and the minimum R(x) is the minimum λ, λ₁.

[Why does the problem statement suggest a more difficult solution?]

1. R(x) = <Ax, x>/<x,x>, x ≠ 0. x ≠ 0 ⇒ <x,x> ≠ 0.
2. .
   • There are an infinite number of solutions to the original eigenvalue problem. One set of these (hatted x) is the one for which x forms an orthonormal basis: ||x|| = 1.
3. when . [ are positive real numbers.]
4. ∴.
5. Now, A has a complete set of real orthonormal eigenvectors is unitary.
   • .
   • 
   • . Hatted x is a matrix and an orthonormal basis for A, and hatted lambda is a vector. A unit vector e wears a hat because that seems to be the fashion for norm 1 variables (though the lambda vector is not necessarily norm 1).
6. The projection of a unit vector e onto fashionable x is .†
   • This rotates or reflects the unit vector into x, without changing its length.
   • The provided hint, ”…to expand an arbitrary unit vector in terms of this [orthonormal] basis,” is misleading; an orthonormal projection does not 'expand' any vector since it does not change its length.
   • If 'expand' is meant to indicate term-by-term expansion, this really seems like a very long road to travel to duplicate the three lines before starting down this track.

10. Show that if x is a vector in ℜ^m, then ||x||_∞ ≤ ||x||₂. Show that if A is a n × n matrix, then ||A||_∞ ≤ √n||A||₂.

1. Note that for all real x, |x| = √x² ≥ 0.
2. Now suppose x_max be an element of x such that max|x_max| = max|x_i|, and let y be a vector in ℜ^m-1 containing all x_i but (one) x_max.
   • = = ≤ =
3. ≤
4. ≤ ✔

Now, part 2. [Note, because this formula layout software does not have an asterisk (of all thing to omit) I have to use a circle or superscript T instead of a star in some places.]

1. ≤
2. ≤
3. ≤
4. ≤
5. ≤
6. ≤
7. ≤ ✔

11. An elementary matrix is one of the form E(u, v; σ) = I − σuv∗. Show that since E(u, v; σ)E(u, v; τ) = E(u, v; σ + τ − στv∗u) that for σ⁻¹ + τ⁻¹ = v∗u, E(u, v; σ) = E(u, v; τ )⁻¹.

1. E(u, v; σ) = E(u, v; τ )⁻¹
2. E(u, v; σ) E(u, v; τ ) = I
3. E(u, v; σ + τ − στv∗u) = I
4. I − (σ + τ − στv∗u)uv∗ = I
5. (σ + τ − στv∗u)uv∗ = 0
6. σuv∗ + τuv∗ − (στv∗u)uv∗ = 0
7. σuv∗ + τuv∗ = (στv∗u)uv∗
8. (σ + τ)uv∗ = (στv∗u)uv∗
9. v∗u(σ + τ)uv∗ = (σ⁻¹ + τ⁻¹)(στv∗u)uv∗
10. v∗u(σ + τ)uv∗ = (τv∗u)uv∗ + (σv∗u)uv∗
11. (σ + τ)v∗uuv∗ = (σ + τ)v∗uuv∗✔
    • These steps are all reversible, so just follow through this backwards.

[“No, Holmes, it is not elementary. And please stop calling me 'dear'.”]

12. A matrix, or the transpose of a matrix, L_i(l_i) = E(l_i, e_i, ; 1), where e_j* l_i = 0, j ≤ i, is called an elementary triangular matrix. Show that L_i⁻¹ (l_i) = L_i(−l_i) and det(L_i) = 1. Explain how these matrices arise in the LU-factorization algorithm.