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)

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

**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

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