Chapter 7 M-File

M-File Home

Script Files

Function Files

ShowLSFit

% 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

ShowLSq

% 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

ShowSparseLS

% 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

ShowQR

% 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.')

ShowBVP

% 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))

CholBench

% 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')

CholErr

% 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

ShowCholBlock

% 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

Rotate

  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

QRrot

  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);

LSq

  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

CholTrid

  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

CholTridSol

  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

CholScalar

  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

CholDot

  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

CholSax

  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

CholGax

  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

CholRecur

  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

CholBlock

  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);

MakeScalar

  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