Problem of the Day

|4220 Home|

May 7  Suppose F =

                                  1         e'            1                  e'

                                  e     (C-iS)         v           (C + iS)*E

                                  e        v'         (-1)^m          v'*E

                                  1    E(C + iS)     Ev          E(C - iS)E

where e = ones(m-1,1),  v = (-1,1,-1,...)', E  is eye(m-1,m-1) with columns in reverse order, and C and S are (m-1)-by-(m-1). Verify that if y  = F*[0; x; 0 ; -E*x] where x is (m-1)-by-1, then  S*x = (1/(2*i))*y(1:m-1). What can you say about

                                            y = F*[x(1); x(2:m); x(m+1);E*x(2:m)]   ?

May 5  Suppose y = Fn*x where n = 4m and Fn is the n-by-n DFT matrix. It can be shown that

             y =  [ I   I   I   I ;  I  -iI   -I   iI ;  I  -I   I  -I;  I   iI   -I   -iI] * [ z1; D * z2; D^2 * z3; D^3 * z4]

where I = eye(m,m), i^2 = -1, z1 = x(1:4:n), z2 = Fm*x(2:4:n), z3 = Fm*x(3:4:n), z4 = Fm*x(4:4:n), and D = diag(d) where d = exp(-2*pi*i/n).^(0:m-1). Develop a flop-efficient, vectorized, recursive implementation of  the following function:

May 3   Assume that A is n-by-n, b is n-by-1 and n = n1*n2. We wish to determine y (n1-by-1) and z (n2-by-1) so that

                                                      phi(y,z) = norm(A*kron(y,z) - b,2)

is minimum. How would you choose y so that phi(y,zc) is minimum where zc (n2-by-1) is fixed.    How would you minimize phi(yc,z) where yc (n1-by-1) is fixed? Formulate an alternating LS procedure whereby y and z are alternately updated.

April 30  Suppose a is 8-by-1. We wish to determine u (2-by-1), v (2-by-1), and w(2,1) so that

                                   phi(u,v,w) = norm(a - kron(u,kron(v,w)),2)

is minimized. How would you choose u so that phi(u,vc,wc) is minimized? How would you choose v so that phi(uc,v,wc) is minimized? How would you choose w so that phi(uc,vc,w) is minimized? Formulate an alternating LS problem whereby u, v, and w are alternately updated.

April 28  How would you use fminbnd to find the area of the largest triangle that can fit inside the ellipse (x/a)^2 + (y/b)^2 = 1.? (By largest we mean in terms of perimeter.)

April 26  Assume that  A(lambda) is an n-by-n matrix whose entries are differentiable functions of the scalar lambda. Define the function F from R^{n+1} to R^{n+1} by F([lambda;x]) = [A(lambda)x ; x'x-1. What is the Jacobian of the function F?

April 23    Suppose A is an n1-by-n1 block matrix with n2-by-n2 blocks. Assume that block (i,j) is zero if |i-j|>1. Assume that the superdiagonal and sub-diagonal blocks are each equal to C (n2-by-n2 and symmetric) . Assume that the diagonal blocks are each equal to eye(n2,n2).  (a) Write a function newX = JacobiStep(x,n1,C,b) that carries out a single Jacobi step for the linear system Ax = b where b is n1*n2-by-1.  (b) Write a function newX = GaussSeidelStep(x,n1,C,b) that carries out a single Gauss-Seidel step for the linear system Ax = b where b is n1*n2-by-1.

April 21 The nonlinear function F = [F1(x) ; ...; Fn(x)] has the property that
the kth component function Fk(x) depend only x(k) and x(n) EXCEPT Fn, which
depends on x(1),...,x(n). Describe a Newton method step.

April 19

(1) Suppose A(theta) is a matrix that depends on a real parameter theta.
Assume that A(theta) is symmetric positive definite for all theta and that
G(theta) is its cholsky factor (assumed to be differentiable).
How would you compute the derivative of the (n,n) entry of G?

(2) Complete:

   function  F = AppxExp(A,j)
 % Suppose R(z) = (1+z/2)/(1 - z/2).
 % F = M^p where p = 2^j and M = R(A/p)

April 16. Suppose t and y are column m-vectors. We wish to fit data points (t(i),y(i)), i=1:m with a function

                                                                          f(t) = a*exp(bt). 

Taking the least squares approach, the idea is to choose a and b so that norm(a*exp(bt) - y,2) is minimized. Show how this problem could be solved using the 1-D minimizer  fminbnd. Hint: for a fixed b, what is the best a? Suggest a heuristic for picking the initial bracketing interval.

April 14. Suppose H is positive definite and  Q'*A*Q = D its schur decomposition. Let s(mu) be the solution to the linear system (H + mu*I)s = -g. How would you compute mu so  that norm(s(mu) = delta where delta >0. Assume that norm(s(0)) > delta. Are there multiple solutions? What solution minimizes .5*s'*H*s + s'*g .

Apr 12  Suppose F(x) is a function from Rn to Rn with the property that its Jacobian J(x) is pentadiagonal for all x. Explain how one can generate a divided difference approximation of J with just five F-evaluations. Hint: for a cleverly chosen vector v, (F(x+hv)-F(x))/h will provide an approximation to J(:,1), J(:,5), J(:,9), etc.

Apr 9. F:Rn-->Rn has the property that F(x) is linear in x(1:m) where 1<=m<n. How could this structure be exploited in an implementation of the finite difference Newton method?

Apr 7 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?

Apr 5   Suppose p(x)= a(1)+a(2)x + a(3)x^2 + a(4)x^3 + x^4 has 4 real roots a<b<c<d. We apply golden section search with initial interval [L,R] . Under what conditions  will the process find a local minimum?

Apr 2  Consider applying the secant method to the function f(x) = (x^2 - A) where A>0. Show that

                        e(k+1)  =  e(k) * e(k -1) / (x(k) + x(k-1) )            where e(j) = x(j) - sqrt(A)

Mar 31   Write a function r = cRoots(a,b,c,d) that returns in the column 3-vector r all the roots of the cubic equation ax^3 + bx^2 + cx + d = 0. Make effective use of fzero

From this we conclude that  |e(k+1) / ( e(k) * e(k-1)) |  converges to a constant C>0.

Mar 29  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 cube" more cubical!

Mar 19 Suppose A has distinct real eigenvalues. (a)  If  [X,D] = eig(A), then how would you evaluate

                           F = I  +  A  +  A^2 / 2!   +  A^3 / 3! + ....

(b) Suppose T is upper triangular with distinct diagonal entries and define

                            F = I  +  T +  T^2 / 2!    + T^3  / 3! + ...

Note that F is upper triangular and that F(i,i) = exp(T(i,i)), i=1:n. From the equation FT = TF, derive recipes for F(i,i+1), i=1:n-1. Then derive recipes for F(i,i+2), i= 1:n-2. Then derive an algorithm for computing all of F.

Mar 17  (a) A is a real matrix and mu is a complex eigenvalue. Suppose Az = mu*z. Show that the real and imaginary parts of the (complex) eigenvector z define a real, 2-dimensional invariant subspace. Suppose the dominant eigenvalue lambda_1 of a real matrix is complex. Note that the conjugate of lambda_1 is also an eigenvalue so the dominant eigenvalue is not unique. What happens if we apply the power method with a real starting vector? Assume that all the other eigenvalues are smaller in maginitude. Suppose A = [w x; y z] has complex conjugate eigenvalues. Show how to compute Q = [c s ;-s c] so that if B = Q'*A*Q then B(1,1) = B(2,2).

Mar 15 (a) If A-sI = QR, why is B = RQ +sI  similar to A? (b) Suppose A and B are m-by-n and have full column rank. What can you say about a vector x that zeros the gradient of f(x) = norm(Ax,2)^2  / norm(Bx,2)^2 >?

Mar 12  Suppose

                               p_{k}(mu) = (a(k) - mu)*p_{k-1}(mu) - c(k)^2 * p_{k-2}(mu)

with p_{0} = 1, p_{1}(mu) = a(1) - mu. Show how to compute the derivative with respect to mu of p_{n}(mu).

Mar 5  Suppose Ax = lambda*x is the dominant eigenpair for a symmetric A. What happens if we apply the power method with the matrix A - lambda*x*x' ?

Mar  3  Show that after a Jacobi update that zeros A(p,q), then off(Anew)^2 = off(A)^2 - 2A(p,q)^2. Here, off(A)^2 is the sum of the squares of the off-diagonal entries.

Mar 1 If you have a function [c,s] = Schur2(A) that computes a rotation that diagonalizes a symmetric 2-by-2 A, show how it can be used to compute a 2-by-2 svd. Hint: How would you determine c and s so that [c s;-s c]'*[w x ; y z] is symmetric.?

Feb 26 Set up a constrained LS problem for determining the best cubic polynomial that fits the data (x(i),y(i)), i=1:m subject to the constraint that the slope of the cubic at x(1) and x(m) are equal.

Feb 24 Use the SVD to characterize

                                     norm(A*x(lambda) - b,2)^2

where x(lambda) solves the linear system  (A'*A + lambda^2*I )x = A'*b.

Feb 22 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.

Feb 19  (a) Suppose u is a unit m-vector and u(1:k-1) = 0. Suppose A is an m-by-n matrix. Carefully show how you would overwrite A with H*A where H = I - 2*u*u'. (b) Given c and s so c^2 + s^2 = 1, find c1 and s1 so that c1^2 + s1^2 = 1 and (I-2vv') = [c s;s -c] where v = [c1;s1].

Feb 17  (a) 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)

(b) Using very little pencil and paper, what is the SVD of A = [10 20 ; 25 -8]?

Feb 15   (a) For an n = 2 SVD, is it always possible to make  both U and V rotations? Recall that a 2-by-2 orthogonal matrix is a rotation if it has the form [c  s ; -s  c]. (b) 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).

Feb 12  A is symmetric positive definite and we compute A = LU. How would you compute lower triangular G so A = GG'?

Feb 10  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)?

Feb 8. What is the exact 1-norm condition of the matrix U = eye(n,n) - triu(ones(n,n),1))?

Feb 5. Suppose x and y are two positive, normalized floating numbers. (Assume IEEE double precision.)  Exactly how many floating point numbers are there  in between x and y given that

                                                           log2(x) = e1 <  e2 = log2(y)

and that neither x or y is an integral power of two.

Feb 3. Suppose A is n-by-n and b, c, and d are column n-vectors. How would you compute x, y, and z so Ax = b, A'y=d and A(n:-1:1,,n:-1:1)z = d? Hint: live off of a single LU factorization.

Feb 1.  Modify this implementation so that it efficiently performs as advertised.

Jan 29. Suppose C, D, E, F, G, H, and Z are initialized n-by-n matrices. Assume that b is a 3n-by-1 vector. How would you solve [C D E; Z F G; Z Z Z]*x = b assuming that Z = zeros(n,n)? You may use the backslash operator.

Jan 27 Suppose A is m-by-n, B is n-by-p, C is p-by-q, and D is q-by-r. How would you compute F = ABCD if minimizing flops is the goal? How does your answer change if the four matrices are all square and B and C have zeros below their diagonal?

Jan 25. Suppose A is an initialized n-by-n matrix and y and z are initialized n-by-1 vectors. What is the dimension of B = z*y'*A and how would you compute it? How many flops would be required by your method?