Problem of the Day

Jan 22  From the syllabus page, download and unzip IntroGUI. Making sure that all 8 files are in the current working directory, type EllipseFit in the command window. Play with the GUI. Do you think it is possible to always find an Ellipse E that is "better: than the conic ellipse C? In other words, is it always possible to play with the sliders so that F(E) < F(C) where F is the objective function encoded in Fit.m.

--------------------

Jan 24  This problem gives you experience with Matlab's colon notation and Boolean capabilities. Write a Matlab function Sudoku(A) that returns "1" if A is a valid Sudoku solution and zero otherwise. Thus, A must be a 9x9 matrix with the property that every column, row, and 3x3 block is made up of distinct positive digits. Hints on checking. If sort(A(5,:)) is the same as 1:9, then row 5 is "OK".  If B= A(4:6,7:9) and sort(B(:)) is the same as (1:9)', then the  (2,3)  block is OK. Some built-in functions that you may find useful: sort, sum, any, all. Solution

--------------------

Jan 27 Suppose A, F, and G are nxn and that F'= F and G' = G. Develop a flop-efficient method for computing N = M*M where M is the 2n-by-2n matrix [A F; G -A']. Don't worry about vectorization. Hint: Let n = 3, generate a sample M explicitly, and then look at N = M*M. E.g. A = rand(3,3), F = rand(3,3); F = F+F'; G = rand(3,3); G = G+G'; M = [A F;G -A']. Solution

Jan 29 Assume that A is an nxn matrix in sparse format. Set up B = [A zeros(n,n) ; zeros(n,n) A] is sparse format.Note that if [i,j,s] = find(A) then s is the vector of A's nonzero entry and i and j are the associated index vectors. Thus, if i(2) = 100 and j(2)=200 then A(100,200) = s(2). Solution

Jan 31 Note that if ek is the k-th column of eye(n,n) and a = SolveHaar(ek), then a is the kth column of inv(Hn). Use this fact to produce a spy plot of H64, the 64-by-64 Haar matrix. Solution

Riddle: What is the most widely studied rank-1 matrix in the history of the world? Hint: You have spent hours yourself studying its properties! Answer: A = (1:9)' * (1:9) !

--------------------

Feb 3   If A = USV' is the svd of A, then what are the singular values of B = A *inv(A' * A)A' assuming that A is nxn and nonsingular? Solution

Feb 5  Using GE, LTriSol, and UTriSol, Show how to compute n-vectors x and y so that solve Ax + Cy = b and A'y = d where A (nxn), C (nxn), b (nx1), and d (nx1) are given. Solution

Feb 7  How would you solve (A^k)x = b using lu. Assume that A is nxn and nonsingular, k is a positive integer, and b is nx1. Solution

--------------------

Feb 10 Without using Matlab, what is the 1-norm condition of A = [.981 .726; .529 .384]? Using Matlab, compute its 2-norm condition from its SVD. Solution

Feb 12 Download ShowHessLU and add a subfunction SolveHessTranspose(m,U,piv,b)  that can be used to solve H' * x = b. Solution

Feb 14  If  A is n-by-n, symmetric, and positive definite, explain why A(1:n-1,1:n-1) is symmetric and positive definite. Suppose b is a given column n-vector. How would you go about computing x and y so that Ax = b and A(1:n-1,1:n-1)y = b(1:n-1)? Hint: Show that if A = GG' ,then G(1:n-1:1:n-1) is the lower triangular Cholesky factor of A(1:n-1,1:n-1). Solution

--------------------

Feb 19 Modify ShowLDLIndef (available on the assignments page) so that it prints  the number of positive, negative, and zero eigenvalues of A. Hint. The Sylvester Law of inertia states that if two matrices are congruent (like A and D) then NumPosEigs(A) = NumPosEigs(D), and NumNegEigs(A) = NumNegEigs(D). Note that the eigenvalues of D are easy to obtain. If D(i,i) is a 1x1 block, then it is an eigenvalue of D. If D(i:i+1,i:i+1) is a 2x2 block, then its eigenvalues are also eigenvalues of D. Solution

Feb 21 We say that A is structurally symmetric if A(i,j) and A(j,i) are either both zero of both non zero for all (i,j) pairs. Suppose A = LU. Does it follow that  L(i,j) and U(j,i) are either both zero or both nonzero? Solution

--------------------

Feb 24 The act of finding sqrt(A) is the act of "building a square" with area A. If we have a rectangle with base x and height A/x, then we can make it "more square" by replacing x with the average of the "current" base and current height, i.e., xNew = (x + A/x)/2. This is precisely the update recipe obtained when we apply Newton's method to get a zero of the  function f(x) = x^2 - A.

To find cube roots we can apply Newton's method to   f(x) = x^3 -A   giving

xNew = x - (x^3-A)/(3x^2) =   ( 2x^3  +  A )/( 3x^2 )

Can you derive this recipe using an intuitive geometric argument similar to the one above for square roots? Start with the fact that we want to build a cube with volume A. Develop an update  that makes the current "box" more cubical!  Solution

Feb 26 Consider the function Bisection.m. In particular, consider the code fragment where the current bracketing interval is about to be halved. Would it be OK to change the comparison  fa*fmid <= 0  to  fa*fmid < 0 ? Solution

Feb 28   How would you use fminbnd to find the area of the largest triangle whose vertices are on the ellipse (x/a)^2 + (y/b)^2 = 1? Solution

-------------------   MIDTERM CUT-OFF -------------------------------------------------------------------------

Mar 3  Under what conditions will Gauss-Seidel converge when A = [a b ; c d]? Solution

Mar 5    Suppose phi(x) = x'Ax/2 - x'*b  where A is symmetric and positive definite. Determine the scalar mu that minimizes phi( x0 - mu*g0)where g0 = Ax0 - b. Note that g0 is the gradient of phi at x0. Solution

Mar 7  Suppose A is symmetric positive definite and that its cholesky factorization A = GG' is known. If u is a vector, how would you compute the 1-norm condition of inv(G)(A + uu')inv(G')? Make use of the Sherman-Morrison formula

inv(C + yz') = inv(C) - alfa inv(C)yz'inv(C) alfa = 1 + y'i*inv(C)*z

-------------------

Mar 10   Review

Mar 12   Exam

Mar 14   For the 3-by-2 case, verify that the method of normal equations is equivalent to zeroing the gradient of

f(x) = (A*x - b)'*(A*x - b)

Solution

-------------------

Mar 17  Suppose [Q,R] = qr(A) where A is mxn and m>n. Show how to minimize || q - b|| where b (mx1) is given and q is orthogonal to every vector in the range of A. Solution

Mar 19  Given matrices A1 (m1xn), A2 (m2xn) and A3 (m3xn) and vectors b1 (m1x1), b2 (m2x1) and b3 (m3x1), how would you compute a vector x that minimizes || A1*x - b1 ||^2  + || A2*x - b2 ||^2  + ||A3*x - b3 ||^2? Solution

Mar 21  How would you determine c and s so that [c s;-s c]'*[w x ; y z] is symmetric.? Hint: Use the subfunction function Givens that is in GivensQR.m. Solution

-------------------

Mar 24  Suppose f and g are column n-vectors. Show how to compute an orthogonal Q as a product of Givens rotations so that Q' * (f * g') * Q  is zero everywhere except in the (1,1) and (1,2) entry. Hint. Put a lot of zeros into f and then put a lot of zeros in g. Hint: If (Q'*f)*(Q'*g)'  has all those zeros then Q'* had better be zero in all but component 1 and Q'*g had better be zero everywhere except components 1 and 2. Solution

Mar 26 What can you say about Q, R, and P if [Q,R,P] = qr(A) and A has more columns than rows. Assuming that A has full row rank, how could you use Q, R, and P to compute a solution to  Ax = b? Solution

Mar 28  Visit JacobiCyclic and consider the update sequence A <-- Vpq'*A;    A<-- A*Vpq. Since the new A is symmetric, show how to reduce the second update from O(n) to O(1). Solution

-------------------

Apr 7  Suppose Q is an orthogonal matrix with the property that Q'SQ = diag(d1,...,dn) where S is the symmetric tridiagonal matrix with two's on the diagonal and minus ones on the sub and superdiagonal. Assume that the eigenvalues are known and that solving linear systems with Q and Q' requires n*log n flops. Let A(alfa,beta) be the symmetric tridiagonal matrix with alfa's down the diagonal and beta's down the sub and superdiagonal. How would you solve A(alfa,beta)x = b for many different choices of alfa and beta? Solution

Apr 9  Suppose A is symmetric, tridiagonal, and positive definite. How would you use the power method to compute the smallest eigenvalue of A?Solution

Apr 11 Look at ShowTriDiagEigs where you will see how to compute pn(x) = det(T-xI) where T is an nxn symmetric tridiagonal matrix.  By differentiating the 3-term recurrence, develop an algorithm that returns both pn(x) and its derivative. Solution

------------------

Apr 14  Suppose M, N, and P are nxn matrices and that A(lambda) = M + lambda*N + lambda^2*P. What is the Jacobian of the function

F:R^{n+1}---> R^{n+1}

defined by F([lambda ; v] ) = [A(lambda)*v ; v'*v - 1 ]. Here, v is an n-vector. Solution

Apr 16   The Jacobian of a function F :R^{n} --> R^{n} is always tridiagonal. Assuming that n is a multiple of 3, what can you say about the entries in the vector

[F(x + delta*I(:,1:3:n)*ones(n/3,1) ) - F(x) ] /delta

where delta is very small? (Display the n=9 case to see what is going on.) Explain how a finite difference approximation of F's Jacobian can be obtained with just four F-evaluations. Solution

Apr 18 How would you use fsolve to find abscissas x1 and x2 and weights w1 and w2 so that if  quad(p) = w1*p(x1) + w2*p(x2), then quad(p) = the integral from -1 to 1 of p(x) where p(x) is a cubic polynomial? Hint. If quad(p) gives the exact integral if p(x) = 1, p(x) = x, p(x) = x^2 and p(x) = x^3, then it will be exact for all cubics by linearity. Solution

------------------

Apr 21 Assume that it costs \$R to paint a unit area red, \$G to paint a unit area green, \$B dollars to paint a unit area blue, and \$W to paint a unit area white. A tetrahedron has its vertices on the unit sphere and each of its faces is painted a different color. What is the most expensive tetrahedron that can be produced in this fashion? Set up an objective function F with the property that if it is maximized, then its value is the cost of the most expensive tetrahedron.Assume the availability of A = Area(phi1,theta1,phi2,theta2,phi3,theta3) that returns the area of the triangle obtained by connecting the sphere points points (phi1,theta1), (phi2,theta2), (phi3,theta3). Convention: the phi's are longitudes and thetas are latitudes. Solution

Apr 23  Assume: F=[F_{1}(x) F_{2}(x) F_{m}(x)]' where F_{k}:R^{n}-->R. Define: phi(x) = F(x)'*F(x)/2. Fact: gradient of phi = J'*F where J is the Jacobian of F. Fact: Hessian of phi = J'*J + F_{1}(x)H_{1}(x) + ... + F_{m}(x)H_{m}(x) where H_{k}(x) is the Hessian of H_{k}(x) Suppose Newton's method is used to find a zero of the gradient and that the Hessian is positive definite at the current point and that sc is the Newton step. Explain why it is possible to choose a step length mu so that phi(xc+mu*sc) < phi(xc). Hint. The inverse of a positive definite matrix is positive definite. Solution

Apr 25 In Levenberg-Marquardt algorithm for nonlinear least squares, the step sc is determined by solving the equation (Jc'*Jc + lambda*I)s = -Jc'*Fc where Jc (mxn) is the current Jacobian and Fc (mx1) . Suppose you have the SVD Jc = USV'. How could fzero be used to determine lambda so that the 2-norm of the solution s = s(lambda) has a prescribed value? Do not worry about starting values. Solution

------------------

Apr 28 A and B are n-by-3 matrices. How could lsqnonlin be used to determine a 3x3 orthogonal matrix Z of the form Z = G1*G2*G3 where

G1 = [1 0 0; 0 cos(theta1) sin(theta1); 0 -sin(theta1) cos(theta1)]

G2 = [cos(theta2) sin(theta2) 0; -sin(theta2) cos(theta2) 0; 0 0 1]

G3 = [1 0 0; 0 cos(theta3) sin(theta3); 0 -sin(theta3) cos(theta3)]

so that norm(AZ - B,'fro') is minimized? Solution

Apr 30 How could the function PageRank (available off syllabus page) be made more efficient by improving the function MyProd that eigs uses? Solution

May 2 Assume the availability of a function Amat = A(lx) that returns an mxn matrix Amat whose entries are a function of the column p-vector x. Assume that b is a column m-vector. Write a function y = F(x,b) that could be handed over to lsqnonlin in order to minimize the 2-norm of

b - A(x) * inv( A(x)'*A(x) ) * A(x)' *b.

Make effective use of the the thin QR factorization. Solution

------------------