Product Triangular Systems With Shift

Matlab Codes

TriProd

  function x = TriProd(R,lambda,b);
%
% x = TriProd(R,lambda,b)
% Computes the solution to 
%
%               (R1*...*Rp - lambda*I)x = b
%
% where R is a px1 cell array containing
% n-by-n upper triangular matrices R{1},..,R{p},
% b is a n-by-1 column vector and lambda is a real scalar.
%
% Uses the implicit method.


p = length(R);
n = length(R{1});

% Calculate the first xc, i.e., x(n)

Rprod = R{1}(n,n);
for i=2:p
   Rprod = Rprod*R{i}(n,n);
end
xc = b(n)/(Rprod-lambda);

% Create the first WC vectors.

WC(:,p) = xc;
for j=p-1:-1:1
   WC(:,j) = R{j+1}(n,n)*WC(:,j+1);
end
% WC(:,j) holds the jth WC vector

for k=1:n-1
   % Compute x(n-k)
   gamma = b(n-k);
   sigProd = 1;
   for j=1:p
      gamma   = gamma - sigProd*(R{j}(n-k,n-k+1:n)*WC(:,j));
      sigProd = sigProd*R{j}(n-k,n-k);
   end 
   gamma = gamma/(sigProd - lambda);
   
   % Update xc
   xc = [gamma;xc];
   
   % Update the WC vectors
   if k<n-1
      omega(p) = gamma;
      for j=p-1:-1:1
         omega(j) = R{j+1}(n-k,n-k)*omega(j+1)+R{j+1}(n-k,n-k+1:n)*WC(:,j+1);
      end
      WC = [omega;WC];
   end
end
x = xc;

 
TriProdExp

  function x = TriProdExp(R,lambda,b);
%
% x = TriProdExp(R,lambda,b)
% Computes the solution to 
%
%               (R1*...*Rp - lambda*I)x = b
%
% where R is a px1 cell array containing
% n-by-n upper triangular matrices R{1},..,R{p},
% b is a n-by-1 column vector and lambda is a real scalar.
%
% Uses the explicit method.

p = length(R);
n = length(R{1});

% Explicitly form the matrix of coefficients
A = R{1};
for j=2:p
   A = A*R{j};
end
A = A - lambda*eye(n,n);

% Apply back-substitution
x = A\b;

QTriProd

  function x = QTriProd(R,lambda,b)
%
% x = QTriProd(R,lambda,b)
% Computes the solution to 
%
%               (R1*...*Rp - lambda*I)x = b
%
% where R is a px1 cell array containing
% n-by-n quasi-upper triangular matrix R{1}
% and n-by-n upper triangular matrices R{2},..,R{p}.
% b is a n-by-1 column vector and lambda is a real scalar.
%
% Uses the implicit method.


p = length(R);
n = length(R{1});

% Do we start with a bump?
if abs(R{1}(n,n-1)) > 0
   m = n-1;
else
   m = n;
end

% Calculate the first xc, i.e., x(m:n)

Rprod = R{1}(m:n,m:n);
for i=2:p
   Rprod = Rprod*R{i}(m:n,m:n);
end
xc = (Rprod-lambda*eye(n-m+1,n-m+1))\b(m:n);

% Create the first WC vectors.

WC(:,p) = xc;
for j=p-1:-1:1
   WC(:,j) = R{j+1}(m:n,m:n)*WC(:,j+1);
end
% WC(:,j) holds the jth WC vector

k=n-m+1;
while k<=n-1
   
   % Have the bottom k components of x, i.e., x(n-k+1:n).
   % Do we go after the next component x(n-k) or the next two 
   % components x(n-k-1) and x(n-k).
   
   m = n-k;
   if k<n-1
      if abs(R{1}(n-k,n-k-1))>0
         m = n-k-1;
      end
   end
  
   % Compute x(m:n-k)
   
   gamma = b(m:n-k);
   sigProd = eye(length(gamma),length(gamma));
   for j=1:p
      gamma = gamma - sigProd*(R{j}(m:n-k,n-k+1:n)*WC(:,j));
      sigProd = sigProd*R{j}(m:n-k,m:n-k);
   end
   gamma = (sigProd - lambda*eye(n-k-m+1,n-k-m+1))\gamma;
   
   % Update xc
   xc = [gamma;xc];
   
   % Update the WC vectors
   if k<n-1
      omega(1:length(gamma),p) = gamma;
      for j=p-1:-1:1
         omega(:,j) = R{j+1}(m:n-k,m:n-k)*omega(:,j+1)+R{j+1}(m:n-k,n-k+1:n)*WC(:,j+1);
      end
      WC = [omega;WC];
   end
   k = k + length(gamma);
end
x = xc;

QTriProdExp

   function x = QTriProdExp(R,lambda,b)
%
% x = QTriProdExp(R,lambda,b)
% Computes the solution to 
%
%               (R1*...*Rp - lambda*I)x = b
%
% where R is a px1 cell array containing
% n-by-n upper quasi-triangular matrix R{1}
% and n-by-n upper triangular matrices R{2},..,R{p}.
% b is a n-by-1 column vector and lambda is a real scalar.
%
% Uses the explicit method.

p = length(R);
n = length(R{1});

% Explicitly form the matrix of coefficients
A = R{p};
for j=p-1:-1:1
   A = TriMult(R{j},A);
end
A = A - lambda*eye(n,n);

% Apply back-substitution

x = zeros(n,1);
i = n;
while (i>=1)
   j=i;
   if i>1
      if abs(A(i-1,i))>0
         j=i-1;
      end
   end
   
   % We solve for x(j:i).
   % Note, j=i-1 means that we have encountered a 2-by-2 bump.
   
   x(j:i) = A(j:i,j:i)\b(j:i);
   if j-1>=1
      b(1:j-1) = b(1:j-1) - A(1:j-1,j:i)*x(j:i);
   end
   i = j-1;
end

   

Test

   % Tests Error Bounds for QTriProdExp and QTriProd
home
randn('seed',0)
n = 20:20:100;
p = [2 5 10 15];
smallnn = [.01 .0001 .000001 .00000001];
m = 100;
for k = 0:length(smallnn)
   disp(' ')
   if k==0
      disp('Well-Conditioned')
   else
      disp(sprintf('smallnn = %8.2e',smallnn(k)))
   end
  
   Ratio12 = zeros(length(n),length(p));
   Ratio21 = zeros(length(n),length(p));
   for i = 1:length(n)
      for j = 1:length(p)
         for sample = 1:m
            if k==0
               [R,lambda,b] = GenProb(n(i),p(j));
            else
               [R,lambda,b] = GenProb(n(i),p(j),smallnn(k));
            end
            x1 = QTriProd(R,lambda,b);    E1 = ErrBound(R,lambda,b,x1);
            x2 = QTriProdExp(R,lambda,b); E2 = ErrBound(R,lambda,b,x2);
            
            Ratio12(i,j) = max(Ratio12(i,j),E1/E2);
            Ratio21(i,j) = max(Ratio21(i,j),E2/E1);
         end
      end
   end 

   disp(' ')
   Ratio12 = Ratio12
   disp(' ')
  
   disp(' ')
   Ratio21 = 1./Ratio21
   disp(' ')

end

  function [R,lambda,b] = GenProb(n,p,smallnn)
%
% n and p are positive integers.
% R is a p-by-1 cell array with R{1} is an n-ny-n
% random upper quasi-triangular matrix and R{2},...,R{p}
% is an n-by-n random upper triangular matrix.
% b is an n-by-1 random vector.
% If nargin == 2 then lambda is chosen so that
%
%      A = R{1}*...*R{n} - lambda*I
%
% is well-conditioned. Otherwise, it is chosen so that
% the (n,n) entry of A equals smallnn.

b = randn(n,1);
R = cell(p,1);
for j = 1:p
   R{j} = triu(randn(n,n));
end

% Make R{1} quasi-triangular.
for i=1:2:n-2
   R{1}(i+1,i) = 1;
end

% Determine lambda
if nargin==2
   lambda = 3*n;
else
   prodnn = 1;
   for i=1:p
      prodnn = prodnn*R{i}(n,n);
   end
   lambda = prodnn - smallnn;
end

   
   

  function e = ErrBound(R,lambda,b,x)
% Computes the relative residual for a computed solution to 
%
%               (R1*...*Rp - lambda*I)x = b
%
% where R is a px1 cell array containing
% n-by-n quasi-upper triangular matrix R{1}
% and n-by-n upper triangular matrices R{2},..,R{p}.
% b is a n-by-1 column vector and lambda is a real scalar.

p = length(R);
n = length(R{1});

% Form z =  R{1}* ... *R{p}x
z = x;
%d = 1;

for j=p:-1:1
   z = R{j}*z;
   
end
% Compute residual vector
res = b - (z - lambda*x);
e = (norm(res)/norm(x));