Problem of the Day
% Monday, August 29
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.
% Wednesday, Aug 31
n = 100;
T = triu(randn(n,n));
clc
disp(' p Error')
disp('----------------------')
for p=50:50:1000
X = randn(n,p);
Y = fSept1(T,X);
error = norm( (T*T*X - Y),'fro');
disp(sprintf('%4d %10.3e',p,error))
end
function Y = fSept1(T,X)
% T is n-by-n upper triangular and X is n-by-p
% Y = T*T*X
Y = T*T*X;
% This involves 2n^3 + 2pn^2 flops. Produce a better implementation by fully exploiting
% T's structure and examining when Y = (T*T)*X is more flop efficient than
% Y = T*(T*X). You might want to vectorize Alg 1.2.1 on p 17 of GVL.
Solve problem P2.3.8 in GVL.
Given a square matrix A, find the closest possible matrix B that has a repeated singular value. To be specific, minimize norm(A-B,2) subject to the constraint that B has a repeated singular value. I can find a B = B_star such that
norm(A-B_star,2) = d/sqrt(2)
where d = min sigma(i)-sigma(i+1), i=1:n-1. But I'm not sure if it is a minimzer.
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.
% Generallize 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];
end
%Sept 11
n =32;
A = randn(n,n);
[L,U] = MyBlockLU(A);
errL = norm(L-tril(L))
errU = norm(U-triu(U))
err = norm(A - L*U)
% Show L and U for 4-by-4 example...
[L,U] = MyBlockLU(randn(4,4))
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?
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?
Show that if Ax = b + r and norm(r,2) = delta, then there exists a matrix E with norm(E,2) = delta so (A+E)x = b.
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?
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?
Suppose H is lower Hessenberg. How would you solve Hx = b?
Suppose A is an n-by-n symmetric matrix. Let mu = max|A(i,j)| and assume that the max occurs off the diagonal. Overwrite A with B = P*A*P' where P is a permutation and |B(2,1)| = mu.
Suppose Tx = b + e where T is symmetric positive definite Toeplitz and x and b are given n-vectors. Can we find a symmetric Toeplitz E so (T + E)x = b? Hint: set up a linear system for the n-scalars that define E.
Without using the gradient, show that norm(Ax-b,2) is minimized by setting x = (A'*A)\(A'*b)
Suppose A is n-by-n. Making effective use of the Matlab QR function, i.e., {Q,R] = qr(A), show how to compute an orthogonal U and a lower triangular L so A = U*L. Hint: If E = I(:,n:-1:1), what can you say about EUE where U is upper triangular?
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.
Suppose A is m-by-n with rank(A) = r < n. How would you find the minimzer of norm(A*x-b,2) which is closest to a given vector d?
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).
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.
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) ?
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]? iIe., Q is a rotation
Suppose W and Z are n-by-p with orthonormal columns. Show that range(W) and range(Z) have anontrivial intersection if and only if W'*Z has a unit singular value.
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?
Suppose A is skew-symmetric. How could you compute Householder matrices P_1,...P_n-1 so that if U = P_1 *...* P_n-2, then U'*A*U is tridiagonal.
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}?
Suppose 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?
If B is n-by-n and bidiagonal. By considering the tridiagonal matrix B'*B, can we conclude that the singular values of B(1:k-1,1:k-1) separate the singular values of B(1:k,1:k) for k=2:n?
Suppose M = [0 A' ; A 0] where A is m-by-n with m>= n. Let p = m+n. If Z'*M*Z = D = diag(d1,...dp) is a Schur decomposition, is it possible to extract U, V and S from Z and D so that U'*AV = S is the SVD of A?
Suppose A is 2-by-2 with complex eigenvalues. Suppose Q = [c s; -s c] withc^2 + s^2 = 1. (a) Is it possible to choose c and s so that if Q'*A*Q = T, then T(1,1) = T(2,2)? (b) Is it possible to choose c and s so that T(1,2) = -T(2,1)?
Suppose T = lambda*I + N is k-by-k with N strictly upper triangular. This means that T has a single eigenvalue lambda and that it has algebraic multiplicity k. Show that if the dimension of null(T - lambda*I ) = r then there exists an orthogonal V such that V'*T*V = [lambda*eye(r,r) T12; 0 T22] where T12 is r-by(n-r) and T22 is upper triangular and (n-r)-by-(n-r).
Nov 11 We have three n-by-n orthogonal matrices partitioned conformably:
Q = [Q1 Q2] U = [U1 U2] V = [V1 V2]
Recall the definition of distance between subspaces (Theorem 2.6.1). Show that
dist(ran(Q1),ran(V1)) <= dist(ran(Q1),ran(U1)) + dist(ran(U1),ran(V1))
Hint: I = U1*U1' + U2*U2'.
Nov 14 Suppose A has all its eigenvalues in the open left half plan. Let f(mu) be the minium singular value of A + mu*i*I.. Show that there is an n-by-n matrix E with norm(E,2) = min f(mu) so that A+E has a purely imaginary eigenvalue.
Nov 16 Suppose A's eigenvalues satisfy |lambda_1| = |lambda_2| > lambda_3| >= ...>= |lambda_n| and that it has a full basis of eigenvectors. Is it possible to deduce lambda_1 and lambda_2 from the vectors produced by the power method?
Nov 18 Suppose A is n-by-n, symmetric, and positive definite. How would you compute Q (n-by-r) with Q'*Q = I so that trace(A*Q*Q') is maximal?
Nov 21 If S is skew-symmetric, then its real Schur decomposition has the form Q'*A*Q = diag(D11,...,Dpp) where each Dii is either 1-by-1 or 2-by-2. Given that this is available, how would you compute M = (I + S)-1(I-S)? What about F = eS?
Nov 23 Suppose A, B, and C are n-by-n. Show that if [0 -A ; I -B] - mu*[I 0 ; 0 C] is singular then so is the matrix M = A + mu*B + mu^2*C.
Nov 28 Suppose after k steps of the Lanczos process we have AQk = Qk*Tk + r*ekT and that r is an eigenvector of A. Does anything special happen during the next step?
Nov 30 Suppose A is symmetric and has a repeated eigenvalue. Show that the Krylov matrix K(A,q,r) = [q Aq An-1q] is singular.
Dec 2 The Rayleigh quotient r(x) = x'*A*x/x'*x is maximized by setting x to be an eigenvector that corresponds to the most positive eigenvalue. Suppose B is p-by-n with p<n. How could you compute the maximizer of r(x) subject to the constraint that x is in the nullspace of B?