Problem of the Day

Oct 22 Suppose A is n-by-n, symmetric, and pentadiagonal. Assume that a(1:n) is the diagonal, b(1:n-1) is the super diagonal, and c(1:n-2) is the super-super diagonal. Let p_k(x) = det(A(1:k,1:k) - x*I). Can you develop a connection between p_k and its "predecessors" p_{k-1}, p_{k-2}, and p_{k-3}?

Oct 20Suppose Q'*A*Q = diag(d1,...,dn) is a schur decomposition. Assume that x is a unit vector and that lambda  is not an eigenvalue of A and that it is closest to dk. If z solves (A - lambda*I)z = x and y = z/norm(z), does it follow that |Q(:,k)'*y| > |Q(:,k)*x| ? That is, the angle between Q(:,k) and y is smaller than the angle between Q(:,k) and k?

Oct 18 (a) Suppose A is symmetric and positive definite. How would you use the Schur decomposition to compute a symmetric and positive definite X so A = X*X? (b) How can the SVD of A be used to compute the Schur decomposition of the matrix B = [tau*eye(n,n) A'; A zeros(m,m)] where A is m-by-n. Assume m>=n for clarity.

Oct 15 Suppose A is a 2-by-2 matrix. How would you compute the cosine-sine pair (c,s) so that if Q = [c s;-s;c] then Q'*A is symmetric. Given that you can do that, how would you compute a 2-by-2 SVD?

Oct 13  Suppose x(lambda) solves  (A'*A - lambda*I)x = A'*b  where A is m-by-n with rank n.  If lambda_1 < lambda_2, what can you say about norm(A*x(lambda_1) - b)  versus norm(A*x(lambda_2_ - b) ?

Oct 11  Suppose A and B are given 2-by-2 matrices. How would you minimize norm(A-BQ,'fro') subject to the constraint that Q = [c s ; -s c]?  ie., Q is a rotatio

Oct 8 Suppose A = USV' is the svd of an m-by-n matrix A. Assume that |V(n,n)| = max(|V(:,n)|. Show that if A = QR then |R(n,n)| <= sqrt(n)*S(n,n).

Oct 6  Show that x_TLS solves  (A'*A - sigma*I)x = A'*b where sigma is the smallest singular value of [A b] and we assume that sigma < smallest singular value of A.

Oct 4 If x = [x1;x2]. How would you determine (c,s) so if y = [c s;-s;c]*x, then y(1) = y(2)? 

Oct 1 How would you use the Matlab QR function to compute the QR factorization of A = Y*Z' where Y and Z are m-by-n matrices with m>n.

Sept 29   Without using the gradient, show that norm(Ax-b,2) is minimized by setting x = (A'*A)\(A'*b)

Sept 27 Prove or disprove. If A is a 3-by-3 symmetric and positive definite then its tridiagonal part B = tril(triu(A,-1),1)) is also symmetric and positive definit.

Sept 20. Suppose L is n-by-n with lower bandwidth p. Thus, the solution to Lx = b involves O(np). Can you deduce via benchmarking whether or not y = L\b runs in O(np) time in Matlab?

Sept 17  (a) If Ax = b+r then (A+E)x = b where E = -r*x'/(x'*x). Thus, if r is "small" then so is E. If A is symmetric and Ax = b+r, can you find a small symmetric F so (A+F)x = b? (b) A 3-by-3 symmetric matrix A with positive diagonal, has the property that det(A(k,k)) > 0 for k = [1 2], [1 3], and [2 3]. Is A positive definite?

Sept 15  If we apply Gaussian Elimination with partial pivoting to a symmetric positive definite matrix A and produce PA = LU, does it follow that P = I?

Sept 13

(a) Suppose y and z are column n-vectors and that e_k  is the kth column of the n-by-n identity I. Under what conditions can you find a column n-vector v so that (I - v*e_k')y = z?

(b) Suppose pVec is a permutation of the integer vector 1:n. Let P be the n-by-n permutation matrix whose i-th row is row pVec(i) of I, the n-by-n identity. Suppose x is a column n-vector. How would you compute y = P'*x?

Sept 10  Generalize this so that it handles general n...
  function [L,U] = MyBlockLU(A)
% A n-by-n with n a power of two.
% Assume that A has an LU factorization
% A = L*U

[n,n] = size(A);
if n==2
    L = [ 1 0 ; A(2,1)/A(1,1) 1];
    U = [ A(1,1) A(1,2) ; 0 A(2,2) - A(1,2)*L(2,1)];
else
    m = n/2;
    % Compare blocks in ...
   
    %    L11   0    U11   U12         A11  A12
    %             *              =                  (all blocks m-by-m)
    %    L21  L22    0    U22         A21  A22
  
    [L11,U11] = MyBlockLU(A(1:m,1:m));              %  L11*U11 = A11
    U12 = L11\A(1:m,m+1:n);                         %  L11*U12 = A12
    L21 = (U11'\A(m+1:n,1:m)')';                    %  L21*U11 = A21
    [L22,U22] = MyBlockLU(A(m+1:n,m+1:n)-L21*U12);  %  L21*U12 + L22*U22 = A22
    L = [L11 zeros(m,m) ; L21 L22];
    U = [U11 U12; zeros(m,m) U22];

Sept 8

If u = I(:,i) and v = I(:,j) then A = alfa*u*v' is the same as A except that the (i,j) entry is now A(i,j)+alfa. The Sherman-Morrison theorem says that A + alfa*u*v' is singular iff 1 + alfa*v'*inv(A)*u = 0.

Now suppose E is k-by-k and that A is nonsingular but n-by-n with n>>k. How can we chose alfa so that Atilde is singular where Atilde is the same as A except Atilde(1:k,1:k) = A(1:k,1:k) + alfa*E? Use the Sherman-Morrison-Woodbury formula on p. 50. Hint: express Atilde as a rank-k modification of A.

Sept 3. (a) What is the 2-norm of the matrix u*v' where u and v are n-vectors? (b) Suppose Ax = y, norm(x,2) = 1. Specify a matrix E in terms of x and y with the property that A+E is singular and norm(E,2) = norm(y,2).

Sept 1 (a) What is the exact value of |fl(1/10) - (1/10)| assuming IEEE floating point arithmetic? Show that

                                                      M<= norm(A,2) <= nM

where M = max(max(abs(A)) and A is n-by-n.

August 30

  m = 100; p = 80; n = 70; 
  A = randn(m,p); A = tril(A);
  B = randn(p,n); B = triu(B);
  C = randn(m,n);
  D = MyProd(A,B,C);
  error = norm( (C+A*B)- D,'fro')

  function D = MyProd(A,B,C)
% A,B, and C are matrices and D = C + A*B;
% A gaxpy-rich implementation...
  [m,p] = size(A); [p,n] = size(B);
  D = C;
  for j=1:n
    D(:,j) = D(:,j) + A*B(:,j);
  end

% Modify this function so that it efficiently handles (in the flop sense) the case
% when A is lower triangular and B is upper triangular. In addition,
% rearrange the computations so that it is rich in outer products.

August 27. How many flops are required to compute B = u*v'*A*v*u' where u and v are n-by-1 and A is n-by-n?

August 25. Write a Matlab script that inputs a real number M and displays what the 2-by-2 matrix A = [1 M; 0 1] does to unit 2-norm vectors. In particular, for a suffiicient number of theta, 0 <=theta<=2pi, display the point (y(1),y(2)) where y = A*[cos(theta);sin(theta)]. You should observe an ellipse. How does the ratio of big-to-little semiaxis grow as M grows?