Chapter 7 M-File
Script Files
ShowLSFit | Fits a line to the square root function. |
ShowLSq | Sensitivity of LS solutions. |
ShowSparseLS | Examines "\" LS solving with sparse arrays. |
ShowQR | Illustrates QRrot. |
ShowBVP | Solves a two-point boundary value problem. |
CholBench | Benchmarks various Cholesky implementations. |
CholErr | Sensitivity of symmetric positive definite systems. |
ShowCholBlock | Explores block size selection in CholBlock. |
Function Files
Rotate | Computes a rotation to zero bottom of 2-vector. |
QRrot | QR factorization via rotations. |
LSq | Least squares solution via QR factorization. |
CholTrid | Cholesky factorization (tridiagonal version). |
CholTridSol | Solves a factored tridiagonal system. |
CholScalar | Cholesky factorization (scalar version). |
CholDot | Cholesky factorization (dot product version). |
CholSax | Cholesky factorization (saxpy version). |
CholGax | Cholesky factorization (gaxpy version). |
CholRecur | Cholesky factorization (recursive version). |
CholBlock | Cholesky factorization (block version). |
MakeScalar | Makes a block matrix a scalar matrix. |
% Script File: ShowLSFit % Displays two LS fits to the function f(x) = sqrt(x) on [.25,1]. close all z = linspace(.25,1); fz = sqrt(z); for m = [2 100 ] x = linspace(.25,1,m)'; A = [ones(m,1) x]; b = sqrt(x); xLS = A\b; alpha = xLS(1); beta = xLS(2); figure plot(z,fz,z,alpha+beta*z,'--') title(sprintf('m = %2.0f, alpha = %10.6f, beta = %10.6f',m,alpha,beta)) end
% Script File: Show LSq % Illustrates LSq on problems with user-specified % dimension and condition. m = input('Enter m: '); n = input('Enter n: '); logcond = input('Enter the log of the desired condition number: '); condA = 10^logcond; [Q1,R1] = qr(randn(m,m)); [Q2,R2] = qr(randn(n,n)); A = Q1(:,1:n)*diag(linspace(1,1/condA,n))*Q2; clc disp(sprintf('m = %2.0f, n = %2.0f, cond(A) = %5.3e',m,n,condA)) b = A*ones(n,1); [x,res] = LSq(A,b); disp(' ') disp('Nonzero residual problem.') disp(' ') disp(' Exact Solution Computed Solution') disp('-----------------------------------------') for i=1:n disp(sprintf(' %20.16f %20.16f',1,x(i))) end b = b+Q1(:,m); [x,res] = LSq(A,b); disp(' ') disp(' ') disp(' ') disp('Zero residual problem.') disp(' ') disp(' Exact Solution Computed Solution') disp('-----------------------------------------') for i=1:n disp(sprintf(' %20.16f %20.16f',1,x(i))) end
% Script File: ShowSparseLS % Illustrate sparse backslash solving for LS problems clc n = 50; disp(' m n full A flops sparse A flops') disp('----------------------------------------') for m = 100:100:1000 A = tril(triu(rand(m,m),-5),5); p = m/n; A = A(:,1:p:m); A_sparse = sparse(A); b = rand(m,1); % Solve an m-by-n LS problem where the A matrix has about % 10 nonzeros per column. In column j these nonzero entries % are more or less A(j*m/n+k,j), k=-5:5. flops(0) x = A\b; f1 = flops; flops(0) x = A_sparse\b; f2 = flops; disp(sprintf('%4d %4d %10d %10d ',m,n,f1,f2)) end
% Script File: ShowQR % Illustrates how the QR factorization is found via Givens Rotations. m = 5; n = 3; disp(' ') disp('Show upper triangularization of ') A = rand(m,n); for j=1:n for i=m:-1:j+1 input('Strike Any Key to Continue'); clc %Zero A(i,j) Before = A [c,s] = Rotate(A(i-1,j),A(i,j)); A(i-1:i,j:n) = [c s; -s c] *A(i-1:i,j:n); disp(sprintf('Zero A(%g,%g)',i,j)) After = A end end disp('Upper Triangularization Complete.')
% Script File: ShowBVP % Illustrates the numerical solution to % % -D^2 u + u(x) = 2xsin(x) - 2cos(x) u(0) = u(pi) = 0. % % Exact solution = u(x) = x*sin(x). close all n = 100; x = linspace(0,pi,n)'; hx = pi/(n-1); d = 2*ones(n-2,1) + hx^2; e = -ones(n-2,1); [g,h] = CholTrid(d,e); b = hx^2*( 2*x(2:n-1).*sin(x(2:n-1)) - 2*cos(x(2:n-1))); umid = CholTridSol(g,h,b); u = [0;umid;0]; plot(x,u,x,x.*sin(x)) err = norm(u - x.*sin(x),'inf'); title('Solution to -D^2 u + u = 2xsin(x) - 2cos(x), u(0)=u(pi)=0') xlabel(sprintf(' n = %3.0f norm(u - xsin(x),''inf'') = %10.6f',n,err))
% Script File: CholBench % Benchmarks five different Cholesky implemntations. clc n = input('Enter n: '); X = rand(2*n,n); A = X'*X; disp(' '); disp(sprintf('n = %2.0f',n)) disp(' '); disp('Algorithm Relative Time Flops') disp('-----------------------------------') flops(0); tic; G = CholScalar(A); t1 = toc; f = flops; disp(sprintf('CholScalar %6.3f %5.0f',t1/t1,f)); flops(0); tic; G = CholDot(A); t2 = toc; f = flops; disp(sprintf('CholDot %6.3f %5.0f',t2/t1,f)); flops(0); tic; G = CholSax(A); t3 = toc; f = flops; disp(sprintf('CholSax %6.3f %5.0f',t3/t1,f)); flops(0); tic; G = CholGax(A); t4 = toc; f = flops; disp(sprintf('CholGax %6.3f %5.0f',t4/t1,f)); flops(0); tic; G = CholRecur(A); t5 = toc; f = flops; disp(sprintf('CholRecur %6.3f %5.0f',t5/t1,f)); disp(' ') disp('Relative time = Time/CholScalar Time')
% Script File: CholErr % Error when solving a Hilbert Matrix System. clc disp(' n cond(A) relerr ') disp('----------------------------------') for n = 2:12 A = hilb(n); b = randn(n,1); G = CholScalar(A); y = LTriSol(G,b); x = UTriSol(G',y); condA = cond(A); xTrue = invhilb(n)*b; relerr = norm(x-xTrue)/norm(xTrue); disp(sprintf(' %2.0f %10.3e %10.3e ',n,condA,relerr)) end
% Script File: ShowCholBlock % Effect of block size on performance of block Cholesky. n = 192; A = randn(n,n); A = A'*A; G0 = chol(A)'; clc disp(sprintf('n = %2.0f',n)) disp(' ') disp(' Block Size Time / Unblocked Time ') disp('---------------------------------------') k =0; for p=[ 8 12 16 24 32 48 96] tic; G = CholBlock(A,p); k=k+1; t(k) = toc; disp(sprintf(' %2.0f %6.3f ',p,t(k)/t(1))) end
function [c,s] = Rotate(x1,x2) % [c,s] = Rotate(x1,x2) % x1 and x2 are real scalars and c and s is a cosine-sine pair so % -s*x1 + c*x2 = 0. if x2==0 c = 1; s = 0; else if abs(x2)>=abs(x1) cotangent = x1/x2; s = 1/sqrt(1+cotangent^2); c = s*cotangent; else tangent = x2/x1; c = 1/sqrt(1+tangent^2); s = c*tangent; end end
function [Q,R] = QRRot(A) % [Q,R] = QRRot(A) % The QR factorization of an m-by-n matrix A. (m>=n). % Q is m-by-m orthogonal and R is m-by-n upper triangular. [m,n] = size(A); Q = eye(m,m); for j=1:n for i=m:-1:j+1 %Zero A(i,j) [c,s] = Rotate(A(i-1,j),A(i,j)); A(i-1:i,j:n) = [c s; -s c]*A(i-1:i,j:n); Q(:,i-1:i) = Q(:,i-1:i)*[c s; -s c]'; end end R = triu(A);
function [xLS,res] = LSq(A,b) % [xLS,res] = LSq(A,b) % Solution to the LS problem min norm(Ax-b) where A is a full % rank m-by-n matrix with m>=n and b is a column m-vector. % xLS is the n-by-1 vector that minimizes the norm(Ax-b) and % res = norm(A*xLS-b). [m,n] = size(A); for j=1:n for i=m:-1:j+1 %Zero A(i,j) [c,s] = Rotate(A(i-1,j),A(i,j)); A(i-1:i,j:n) = [c s; -s c]*A(i-1:i,j:n); b(i-1:i) = [c s; -s c]*b(i-1:i); end end xLS = UTriSol(A(1:n,1:n),b(1:n)); if m==n res = 0; else res = norm(b(n+1:m)); end
function [g,h] = CholTrid(d,e) % G = CholTrid(d,e) % Cholesky factorization of a symmetric, tridiagonal positive definite matrix A. % d and e are column n-vectors with the property that % A = diag(d) + diag(e(2:n),-1) + diag(e(2:n),1) % % g and h are column n-vectors with the property that the lower bidiagonal % G = diag(g) + diag(h(2:n),-1) satisfies A = GG^T. n = length(d); g = zeros(n,1); h = zeros(n,1); g(1) = sqrt(d(1)); for i=2:n h(i) = e(i)/g(i-1); g(i) = sqrt(d(i) - h(i)^2); end
function x = CholTridSol(g,h,b) % x = CholTridSol(g,h,b) % Solves the linear system G*G'x = b where b is a column n-vector and % G is a nonsingular lower bidiagonal matrix. g and h are column n-vectors % with G = diag(g) + diag(h(2:n),-1). n = length(g); y = zeros(n,1); % Solve Gy = b y(1) = b(1)/g(1); for k=2:n y(k) = (b(k) - h(k)*y(k-1))/g(k); end % Solve G'x = y x = zeros(n,1); x(n) = y(n)/g(n); for k=n-1:-1:1 x(k) = (y(k) - h(k+1)*x(k+1))/g(k); end
function G = CholScalar(A) % G = CholScalar(A) % Cholesky factorization of a symmetric and positive definite matrix A. % G is lower triangular so A = G*G'. [n,n] = size(A); G = zeros(n,n); for i=1:n % Compute G(i,1:i) for j=1:i s = A(j,i); for k=1:j-1 s = s - G(j,k)*G(i,k); end if j < i G(i,j) = s/G(j,j); else G(i,i) = sqrt(s); end end end
function G = CholDot(A) % G = CholDot(A) % Cholesky factorization of a symmetric and positive definite matrix A. % G is lower triangular so A = G*G'. [n,n] = size(A); G = zeros(n,n); for i=1:n % Compute G(i,1:i) for j=1:i if j==1 s = A(j,i); else s = A(j,i) - G(j,1:j-1)*G(i,1:j-1)'; end if j < i G(i,j) = s/G(j,j); else G(i,i) = sqrt(s); end end end
function G = CholSax(A) % G = CholSax(A) % Cholesky factorization of a symmetric and positive definite matrix A. % G is lower triangular so A = G*G'. [n,n] = size(A); G = zeros(n,n); s = zeros(n,1); for j=1:n s(j:n) = A(j:n,j); for k=1:j-1 s(j:n) = s(j:n) - G(j:n,k)*G(j,k); end G(j:n,j) = s(j:n)/sqrt(s(j)); end
function G = CholGax(A) % G = CholGax(A) % Cholesky factorization of a symmetric and positive definite matrix A. % G is lower triangular so A = G*G'. [n,n] = size(A); G = zeros(n,n); s = zeros(n,1); for j=1:n if j==1 s(j:n) = A(j:n,j); else s(j:n) = A(j:n,j) - G(j:n,1:j-1)*G(j,1:j-1)'; end G(j:n,j) = s(j:n)/sqrt(s(j)); end
function G = CholRecur(A) % G = CholRecur(A) % Cholesky factorization of a symmetric and positive definite matrix A. % G is lower triangular so A = G*G'. [n,n] = size(A); G = zeros(n,n); if n==1 G = sqrt(A); else G(1:n-1,1:n-1) = CholRecur(A(1:n-1,1:n-1)); G(n,1:n-1) = LTriSol(G(1:n-1,1:n-1),A(1:n-1,n))'; G(n,n) = sqrt(A(n,n) - G(n,1:n-1)*G(n,1:n-1)'); end
function G = CholBlock(A,p) % G = CholBlock(A,p) % Cholesky factorization of a symmetric and positive definite n-by-n matrix A. % G is lower triangular so A = G*G'. p is the block size and must divide n. % Represent A and G as m-by-m block matrices where m = n/p. [n,n] = size(A); m = n/p; A = MakeBlock(A,p); G = cell(m,m); for i=1:m for j=1:i S = A{i,j}; for k=1:j-1 S = S - G{i,k}*G{j,k}'; end if j < i G{i,j} = (G{j,j}\S')'; else G{i,i} = CholScalar(S); end end end % Convert G to a scalar matrix. G = MakeScalar(G);
function A = MakeScalar(A_block) % A = MakeScalar(A_block) % Represents the m-by-m block matrix A_block as an n-by-n matrix of scalrs with % where each block is p-by-p and n=mp. [m,m] = size(A_block); [p,p] = size(A_block{1,1}); for i=1:m for j=1:m if ~isempty(A_block{i,j}) A(1+(i-1)*p:i*p,1+(j-1)*p:j*p) = A_block{i,j}; end end end