== Assignment 2 == //{{nla:assign2.pdf|Numerical Linear Algebra February 7, 2008 1}}// --- //[[honeypot@handbasket.org|David Wagner]] 2008/02/14 11:14// //Work any 7 of the following problems. You make work more for additional credit.// ====== ✔Problem 1 ====== //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₂)||², ﬁnd a formula for the projection of b onto the subspace L. Answer: if A = [a₁ a₂], then P = A(A*A)^-1A*.// - Assume given are real scalars or vectors in ℜⁿ; * A=[a₁ a₂]; Pb = L = A[s t]*; * A=[matrix{2}{2}{a_11 a_21 a_12 a_22}]; A*=[matrix{2}{2}{a_11 a_12 a_21 a_22}]. //Note the row and column indices.// - phi(s,t) = {1/2}||b- (s a_1 + t a_2)||^2 = {1/2}sqrt{Sigma(b- (s a_1 + t a_2))^2}^2 = {1/2}Sigma(b- (s a_1 + t a_2))^2 = {1/2}((b_1- (s a_11 + t a_21))^2 + (b_2- (s a_12 + t a_22))^2)={1/2}( {b_1}^2 -2b_1(s a_11 + t a_21) + (s a_11 + t a_21)^2 + {b_2}^2 -2b_2(s a_12 + t a_22) + (s a_12 + t a_22)^2)= {1/2}({b_1}^2 -2b_1 s a_11 -2b_1 t a_21 +s^2 {a_11}^2 +2st a_11 a_21 +t^2 {a_21}^2 +{b_2}^2 -2b_2 s a_12 -2b_2 t a_22 +s^2 {a_12}^2 +2st a_12 a_22 +t^2 {a_22}^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. - {{partial phi}/{partial s}}={1/2}(-2 b_1 a_11 + 2{a_11}^2 s+2t a_11 a_21 -2b_2 a_12 +2{a_12}^2s +2 a_12 a_22 t )=-b_1 a_11 + {a_11}^2 s+t a_11 a_21 -b_2 a_12 +{a_12}^2s + a_12 a_22 t=0 - {{partial phi}/{partial t}}={1/2}(-2 b_1 a_21 +2s a_11 a_21 +2{a_21}^2 t -2b_2 a_22 +2s a_12 a_22 +2{a_22}^2 t)=-b_1 a_21 +s a_11 a_21 +{a_21}^2 t -b_2 a_22 +s a_12 a_22 +{a_22}^2 t=0 - {{partial^2 phi}/{partial s^2}}={a_11}^2 +{a_12}^2 >0✔ - {{partial^2 phi}/{partial t^2}}={a_21}^2 +{a_22}^2 >0✔ - ({a_11}^2+{a_12}^2)s +(a_11 a_21 +a_12 a_22)t=b_1 a_11+b_2 a_12 - ({a_21}^2+{a_22}^2)t +(a_11 a_21 +a_12 a_22)s =b_1 a_21 +b_2 a_22 - [matrix{2}{2}{({a_11}^2+{a_12}^2) (a_11 a_21 +a_12 a_22) (a_11 a_21 +a_12 a_22) ({a_21}^2+{a_22}^2)}][matrix{2}{1}{s t}]=[matrix{2}{2}{a_11 a_12 a_21 a_22}][matrix{2}{1}{b_1 b_2}] - ([matrix{2}{2}{a_11 a_12 a_21 a_22 }][matrix{2}{2}{a_11 a_21 a_12 a_22}])[matrix{2}{1}{s t}]=[matrix{2}{2}{a_11 a_12 a_21 a_22}][matrix{2}{1}{b_1 b_2}] - A*A[s t]* = A*b. //Recall the row and column indices are as noted above.// - (A*A)^-1A*b = [s t]* - A(A*A)^-1A*b = A[s t]* = Pb - -> P = A(A*A)^-1A*✔ ====== Problem 2 ====== //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.)// - i,j,k in [1,n] wedge bbN. Assume real numbers. - Ax = [matrix{3}{1}{{sum{j=1}{n}{A_{1j} x_j}} vdots {sum{j=1}{n}{A_{nj} x_j}} }]; (Ax)^2 = [matrix{3}{1}{{sum{j=1}{n}{a_{1j}x_j} sum{i=1}{n}{sum{j=1}{n}{a_{ij}x_j}}} vdots {sum{j=1}{n}{a_{nj}x_j} sum{i=1}{n}{sum{j=1}{n}{a_{ij}x_j}}} } ] - phi(x) = {1/2}sqrt{sum{i=1}{n}{{(Ax-b)_i}^2}}^2 = {1/2}sum{i=1}{n}{{(Ax-b)_i}^2}= {1/2}sum{i=1}{n}((Ax)^2 - 2bAx + b^2)_i= {1/2}sum{i=1}{n}(Ax(Ax - 2b) + b^2)_i - phi(x_k) = {1/2}sum{i=1}{n}(sum{j=1}{n}{a_{kj}x_j} sum{h=1}{n}{sum{j=1}{n}{a_{hj}x_j}} - 2b sum{j=1}{n}{A_{kj} x_j} + b^2)_i = {1/2}sum{i=1}{n}( sum{j=1}{n}{A_{kj} x_j}(sum{h=1}{n}{sum{j=1}{n}{a_{hj}x_j}} - 2b) + b^2)_i * 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. - phi(x) = {1/2}sqrt{(Ax-b) star (Ax-b)}^2 = ½(Ax-b)*(Ax-b) - {{partial phi}/{partial x_k}} ={1/2} 2x_k sum{j=1}{n}{Ajk x_k FIXME } -2b A_{jk} = 0FIXME - - phi(x_1,...,x_n) ====== ✔Problem 3 ====== //3. Consider the matrix// A = (matrix{3}{2}{ 1 2 0 1 1 0}) //By hand calculation ﬁnd the orthogonal projection P onto the range(A) and the image under P of the vector (1, 2, 3)^T. (See problem 1.)// - P = P² = range(A) = Ax; v=[1 2 3]^T. Find P, Pv. - A* = [matrix{2}{3}{1 0 1 2 1 0}] - P = A(A*A)^-1A* = [matrix{3}{2}{1 2 0 1 1 0}] ([matrix{2}{3}{1 0 1 2 1 0}] [matrix{3}{2}{1 2 0 1 1 0}])^{-1} [matrix{2}{3}{1 0 1 2 1 0}] = [matrix{3}{2}{1 2 0 1 1 0}] [matrix{2}{2}{2 2 2 5}]^{-1} [matrix{2}{3}{1 0 1 2 1 0}] = [matrix{3}{2}{1 2 0 1 1 0}] [matrix{2}{3}{{1/6} {-{1/3}} {{5/6}} {1/3} {1/3} {-{1/3}} }] =[matrix{3}{3}{{5/6} {1/3} {1/6} {1/3} {1/3} {-{1/3}} {1/6} {-{1/3}} {5/6}}] * P = {1/6}[matrix{3}{3}{5 2 1 2 2 {-2} 1 {-2} 5}] - Pv = {1/6}[matrix{3}{3}{5 2 1 2 2 {-2} 1 {-2} 5}][matrix{3}{1}{1 2 3}] = [matrix{3}{1}{2 0 2}] ====== Problem 4 ====== //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 =i, v_j>. (Use the fact that = 0 for all v in L.) How is this related to the stiﬀness matrix used in ﬁnite element methods?// ====== Problem 5 ====== //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)* - (I - 2P)*(I - 2P) - = (I* - 2P*)(I - 2P) - = I*I - 2P*I - 2I*P +4P*P - = I -2P* -2P +4PP - = I -2P -2P +4P - = I✔ The geometric interpretation has something to do with the fundamental theorem of linear algebra and projections into the fundamental subspaces. //FIXME Needs diagram// ====== ✔Problem 6 ====== //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 a_{ij} = a_{ji} = overline{a_{ij}}. - i,j = 1,2,...,n - Ax = λx = A*x - Ax_k = λ_kx_k = A*x_k. - lambda_k x_k = sum{j=1}{n}{a_{jk}} x_{jk} = a_{1k}x_{1k} +a_{2k}x_{2k} + cdots + a_{nk}x_{nk} = a_k*x_k = a_k x_k. * Since a_k is real, lambda_k must be real. - (λ_k x_k)*x_k = (a_k x_k)*x_k. //[Dot both sides.]// - λ_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. - Ax = λ_xx = A*x; Ay = λ_yy = A*y; λ_x≠λ_y - Ax = λ_xx - (Ax)*y = (λ_xx)*y //[Dot both sides again]// - x*(A*y) = λ_x(x*y) //[More dancing, this time a cha-cha.]// - x*(λ_yy) = λ_x(x*y) - λ_y*(x*y) = λ_x(x*y) * The complex conjugate of a real number is itself: λ_y* = λ_y. - λ_y(x*y) = λ_x(x*y) * Since λ_x and λ_y are distinct and (assumed) nonzero, (x*y) must be zero. - ∴ x is orthogonal to y. ====== ✔Problem 7 ====== //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).// - Ux = λx - U⁻¹Ux = U⁻¹λx = U*x - Ix = U*λx - x_k=U*λ_kx_k - [matrix{4}{1}{x_{1k} x_{2k} vdots x_{jk}}] = [matrix{4}{1}{ sum{j=1}{n}{u_{j1} lambda_k x_{jk}} sum{j=1}{n}{u_{j2} lambda_k x_{jk}} vdots sum{j=1}{n}{u_{jn} lambda_k x_{jk}} }] = lambda_k [matrix{4}{1}{ {u_{11}x_{1k} + u_{21}x_{2k} + cdots + u_{j1}x_{jk}} {u_{12}x_{1k} + u_{22}x_{2k} + cdots + u_{j2}x_{jk}} vdots {u_{1n}x_{1k} + u_{2n}x_{2k} + cdots + u_{jn}x_{jk}} }] = lambda_k [matrix{4}{1}{ {u_{11} + u_{21}+ cdots + u_{j1}} {u_{12} + u_{22} + cdots + u_{j2}} vdots {u_{1n} + u_{2n} + cdots + u_{jn}} }] x_k - x_k = lambda_k sum{j=1}{n}(u_{jk} x_k) = lambda_k sum{j=1}{n}(u_{jk}) x_k - lambda_k = 1/{sum{j=1}{n}{u_{jk}} } ^Each eigenvalue of a unitary matrix U is the multiplicative inverse of the sum of the corresponding matrix column.| ====== Problem 8 ====== //8. Let A be the operator deﬁned 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 ﬁnite diﬀerence 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 ﬁnding an eigenvector-eigenvalue pair amounts to solving a __diﬀerence__ 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.)// ====== ✔Problem 9 ====== //9. Let A be a symmetric matrix of order n, with eigenvalues λ₁ ≤ λ₂ ≤ · · · ≤ λ_n and deﬁne the Rayleigh quotient R(x) = /, x ≠ 0. Show that max_x R(x) = λ_n and min_x R(x) = λ₁. (Show that optimizing ˜R(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.)// - A ∈ ℜ^nxn; i ∈ [1,2,3,...,n]. - Ax = λx = (Ax)*x = (λx)*x = λ(x*x). //[Dot both sides by x, then apply scalar linearity.]// - λ = (Ax)*x/(x*x) = / = R(x). //[x ≠ 0 => (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?]// - R(x) = /, x ≠ 0. x ≠ 0 => ≠ 0. - R(x) = {1/{}} = {1/{}} hat{R}(x). * 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. - R(x) = hat{R}(hat{x}) when hat{x} = x : {1/{}} = 1 = = sqrt{} = vert x vert_2 = 1. //[// and // are positive real numbers.]// - ∴hat{R}(hat{x}) = , vert hat{x} vert_2 = 1. - Now, A has a complete set of real orthonormal eigenvectors hat{x} = [hat{x_1} | hat{x_2} |...|hat{x_n}]; hat{x} is unitary. * vert hat{x} vert_2 = vert hat{x_i} vert_2 = 1. * = 0:i<>j * A hat{x} = I hat{lambda} hat{x}. //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).// - The projection of a unit vector e onto fashionable x is P hat{e} = (hat{x} hat{x}^T) hat{e} = hat{x}.† * 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. ====== ✔Problem 10 ====== //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||₂.// - Note that for all real x, |x| = √x² ≥ 0. - 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. * max |x_i| = |x_{max}| = sqrt{{|x_{max}|}^2} ≤ sqrt{|{x_{max}}|^2 + sum{j=1}{m-1}{|y_j|^2}} = sqrt{sum{i=1}{m}{|x_i|^2}} - max |x_i| ≤ sqrt{sum{i=1}{m}{|x_i|^2}} - vert x vert_{infty} ≤ vert x vert_2✔ 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.]// - vert {a_i} vert_{infty} ≤ vert {a_i} vert_2 - {max}under{1<=i<=n} {|a_i|} ≤ sqrt{sum{i=1}{n}{|a_i|^2}} - {max}under{1<=i<=n} {|a_i|^2} ≤ sum{i=1}{n}{|a_i|^2} - {max}under{1<=i<=n} {|x_i|^2|a_i|^2} ≤ |x_i|^2 sum{i=1}{n}{|a_i|^2} - {max}under{1<=i<=n} sqrt{|x_i|^2|a_i|^2} ≤ sqrt{n sum{i=1}{n}{|x_i|^2|a_i|^2}} - {max}under{1<=i<=n} |x_i||a_i| ≤ sqrt{n} sum{i=1}{n}{|x_i|^2|a_i|^2} - vert A vert_infty ≤ sqrt{n} vert A vert_2✔ ====== ✔Problem 11 ====== //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; τ )⁻¹.// - E(u, v; σ) = E(u, v; τ )⁻¹ - E(u, v; σ) E(u, v; τ ) = I - E(u, v; σ + τ − στv∗u) = I - I − (σ + τ − στv∗u)uv∗ = I - (σ + τ − στv∗u)uv∗ = 0 - σuv∗ + τuv∗ − (στv∗u)uv∗ = 0 - σuv∗ + τuv∗ = (στv∗u)uv∗ - (σ + τ)uv∗ = (στv∗u)uv∗ - v∗u(σ + τ)uv∗ = (σ⁻¹ + τ⁻¹)(στv∗u)uv∗ - v∗u(σ + τ)uv∗ = (τv∗u)uv∗ + (σv∗u)uv∗ - (σ + τ)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'."]// ====== Problem 12 ====== //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.//