Rational Matrix Functions and Rank-One Updates

(Matlab Codes)

Daniel S. Bernstein  and Charles F. Van Loan

 

BooksPapers  |  Research   |  Biographical   | MiscellaneousHome |

 

 

Contents

Kry Computes Krylov Matrix
Rational_Update Computes R(A+uv') from R(A)
Pade Computes specified Pade Approximate of exp(A)
TestA Checks Rational_Update on Pade Examples
phi "Exact" phi vector using expM
phi1 phi vector using powers of a Taylor approximation to exp
phi2 phi-vector using a single Taylor approximatyion to exp
TestB Compares Efficiency and Accuracy of phi1 and phi2
Jac1 Finite Difference Approximation of phi1's Jacobian
Jac2 Finite Difference Approximation of phi2's Jacobian
FastJac2 Fast Computation  of phi2's Jacobian
TestC Compares Jac2 and FastJac2

 

Krylov Matrix Computation

  function K = Kry(A,X,j)
% A is an n-by-n matrix, X is an n-by-p matrix, and j is a positive integer.
% K = [ X , A*X , A^2*X , ... , A^j-1*X ]
 
  [n,p] = size(X);
  K = X;
  for i=2:j
     K = [ K A*K(:,(i-2)*p+1:(i-1)*p)];
  end

|top|


Update of a Rational Matrix Function

  function [XN,YN,XD,YD,X,Y] = Rational_Update(A,u,v,alfa,N,beta,Q,R)
%
% Assume that A is an n-by-n matrix and that u and v are column n-vectors.
% The n-by-n matrices Q and R define a QR factorization of a nonsingular
% matrix D. This matrix and the n-by-n matrix N are polynomials in A:      
%     
%             N = alfa(0)I + alfa(1)A + ... + alfa(r)A^r
%	      D = beta(0)I + beta(1)A + ... + beta(q)A^q
%
%  The column vectors alfa(1:r) and beta(1:q) contain the coefficients,
%  with the exception of alfa(0) and beta(0), which are not required. Note
%  that alfa and beta can be empty vectors when the underlying degree is 0.
% 
%  The output matrices XN,YN,XD,YD,X, and Y relate to
%
%             N_ = alfa(0)I + alfa(1)(A+uv') + ... + alfa(p)(A+uv')^r
%	      D_ = beta(0)I + beta(1)(A+uv') + ... + beta(q)(A+uv')^q
%
%  in the following way:
%
%               N_ = N + XN*YN'           XN and YN are n-by-r
%               D_ = N + XD*YD'           XD and YD are n-by-q
%              (D_\N_) = (D\N) + X*Y'     X  and Y  are n-by-d, d = max([r q])
%
%  Exceptions occur if r or q is zero:
%
%       If r=0, then XN=YN=zeros(n,1). 
%       If q=0, then XD=YD=zeros(n,1). 
%       If r=0 and q=0, then X=Y=zeros(n,1).
%
%  Abbreviated calls are allowed for the polynomial case (alfa ~= [], beta = []):
%
%             [XN,YN] = Rational_Update(A,u,v,alfa) 
%

% 4-argument call:

   if nargin==4
      d = length(alfa);
      XN = Kry(A,u,d); 
      YN = Kry(A'+v*u',v,d)*hankel(alfa); 
   else

% General call:

% The special case of a constant numerator polynomial:

      if alfa==[]
         alfa=0; 
      end
   
% The special case of a constant denominator polynomial:

      if beta==[], 
         beta=0; 
      end
   
      [n,n]=size(A); 
      r=length(alfa); 
      q=length(beta); 
      d=max([r q]);

% Compute the Krylov and Hankel matrices.
   
      Kd = Kry(A,u,d); 
      Ld = Kry(A'+v*u',v,d);   
      H_alfa = hankel(alfa); 
      H_beta = hankel(beta);
 
% Compute the numerator and denominator perturbations:
     
      XN = Kd(:,1:r); YN = Ld(:,1:r)*H_alfa; 
      XD = Kd(:,1:q); YD = Ld(:,1:q)*H_beta; 
   
% Now get the rational's perturbations:
     
      YN_tilde = [YN zeros(n,d-r)]; 
      YD_tilde = [YD zeros(n,d-q)];
      X = R\(Q'*Kd);   
      M = eye(d,d) + YD_tilde'*X;
      Y = (YN_tilde' - (M\YD_tilde') * (R\(Q'*(N + XN*YN'))))';
   end
   

|top|


Pade Approximation

  function [alfa,beta,N,D] = Pade(p,q,A)
%
% Computes the Pade approximation D^{-1}N where
%
%     N = N(p,q) = I + alfa(1)A + ... + alfa(p)A^p
%
%     D = D(p,q) = I + beta(1)A + ... + beta(q)A^q
%
% A call of the form Pade(p,q) returns only alfa(1:p) and beta(1:q).
%
% 
 
  alfa = zeros(p,1);
  a = 1;
  for k=1:p
     a = a*((p-k+1)/k)/(p+q-k+1); 
     alfa(k) = a;
  end

  beta = zeros(q,1);
  b = 1;
  for k=1:q
     b = -b*((q-k+1)/k)/(p+q-k+1);
     beta(k) = b;
  end
   
  if nargin==3
     [n,n] = size(A);
     N = eye(n,n);
     D = eye(n,n);
     Apower = A;
     for k=1:max([p q])
	if k<=p 
           N = N + alfa(k)*Apower;
	end
	if k<=q
           D = D + beta(k)*Apower;
	end
	if k<max([p q])
	   Apower = A*Apower;
        end
     end
  end

|top|


Test A

% TestA.m
% Confirms the correctness of Rational_Update.
% Uses (p,q) Pade approximations to  the matrix exponential.
% The A-matrix is randomly chosen and of modest norm so as to
% avoid conditioning.

% Set the size of the problem.

n = 10; 

% Set up A and A1 = A + uv':


A = randn(n,n);  A = A/(n*norm(A));
u = randn(n,1);  u = u/(n*norm(u));
v = randn(n,1);  v = v/(n*norm(v));
A1 = A+u*v';


clc
home
disp('Let F  = D\N   be the (p,q) Pade approximation of exp(A).')
disp('Let F1 = D1\N1 be the (p,q) Pade approximation of 
exp(A+uv'').')
disp('N-Error = norm(N1 - (N+XN*YN''))')
disp('D-Error = norm(D1 - (D+XD*YD''))')
disp('F-Error = norm(F1 - (F+X*Y''))')
disp(' ')


disp(' p   q    N1-Error     D1-Error     F1-error')
disp('---------------------------------------------')
for p=0:4
   for q=0:4

% Generate the (p,q) Pade approximate of A and A1.

      [alfa,beta,N,D]   = Pade(p,q,A);
	[alfa,beta,N1,D1] = Pade(p,q,A1);
	[Q,R] = qr(D);
	  
	[XN,YN,XD,YD,X,Y] = Rational_Update(A,u,v,alfa,N,beta,Q,R);

% Compute the error in the updated numerator, denominator, and F

      eN = norm(N1-(N+XN*YN'));
	eD = norm(D1-(D+XD*YD'));
	eF = norm((D1\N1)-((D\N)+X*Y'));

	disp(sprintf('%2d  %2d  %10.3e   %10.3e   %10.3e',p,q,eN,eD,eF))
   end
end
disp(' ')
disp('        (Now test 4-argument calls.)')
disp(' ')

q=0;
for p=1:8

% Generate the (p,0) Pade approximate of A and A1.

   [alfa,beta,N,D]   = Pade(p,q,A);
   [alfa,beta,N1,D1] = Pade(p,q,A1);	  
   [XN,YN] = Rational_Update(A,u,v,alfa);
	  
% Compute the error in the updated numerator.

   eN = norm(N1-(N+XN*YN'));
   eD = 0;
   eF = eN;
   disp(sprintf('%2d  %2d  %10.3e   %10.3e   %10.3e',p,q,eN,eD,eF))
 
end

Test A Output:

Let F  = D\N   be the (p,q) Pade approximation of exp(A).
Let F1 = D1\N1 be the (p,q) Pade approximation of exp(A+uv').
N-Error = norm(N1 - (N+XN*YN'))
D-Error = norm(D1 - (D+XD*YD'))
F-Error = norm(F1 - (F+X*Y'))
 
 p   q    N1-Error     D1-Error     F1-error
---------------------------------------------
 0   0  -0.000e+00   -0.000e+00   -0.000e+00
 0   1  -0.000e+00    2.220e-16    6.706e-16
 0   2  -0.000e+00    2.231e-16    4.447e-16
 0   3  -0.000e+00    2.244e-16    4.465e-16
 0   4  -0.000e+00    2.237e-16    5.571e-16
 1   0   2.220e-16   -0.000e+00    2.220e-16
 1   1   2.220e-16    2.220e-16    6.668e-16
 1   2   2.220e-16    2.221e-16    6.667e-16
 1   3   2.220e-16    2.231e-16    4.479e-16
 1   4   2.220e-16    4.443e-16    5.573e-16
 2   0   2.240e-16   -0.000e+00    2.240e-16
 2   1   1.120e-16    2.224e-16    6.668e-16
 2   2   2.242e-16    2.221e-16    4.461e-16
 2   3   1.111e-16    2.221e-16    4.525e-16
 2   4   2.223e-16    2.222e-16    6.680e-16
 3   0   4.443e-16   -0.000e+00    4.443e-16
 3   1   2.222e-16    1.110e-16    4.478e-16
 3   2   2.230e-16    2.225e-16    4.540e-16
 3   3   2.222e-16    2.232e-16    4.584e-16
 3   4   2.221e-16    2.223e-16    6.728e-16
 4   0   4.444e-16   -0.000e+00    4.444e-16
 4   1   1.479e-17    2.221e-16    6.672e-16
 4   2   2.223e-16    2.221e-16    4.528e-16
 4   3   2.255e-16    4.441e-16    6.665e-16
 4   4   4.441e-16    2.244e-16    2.257e-16
 
        (Now test 4-argument calls.)
 
 1   0   2.220e-16    0.000e+00    2.220e-16
 2   0   2.240e-16    0.000e+00    2.240e-16
 3   0   4.443e-16    0.000e+00    4.443e-16
 4   0   4.444e-16    0.000e+00    4.444e-16
 5   0   4.444e-16    0.000e+00    4.444e-16
 6   0   6.664e-16    0.000e+00    6.664e-16
 7   0   6.664e-16    0.000e+00    6.664e-16
 8   0   6.666e-16    0.000e+00    6.666e-16

|top|


The Exact Vector phi

  function z = phi(A,t,c,b)
%
% Assume that A is n-by-n and that c and b are column n-vectors.
% Assume that t is a length-m vector with 0 < t(1) < ... < t(m).
% Define X(:,0) = b and delk = t(k)-t(k-1) with t(0)=0;
% Computes a column m-vector z 
%
%                          c'*X(:,1)
%                    z  =     :
%                          c'*X(:,m)
%
% where X(:,k) = exp(delk*A)*X(:,k-1).
%

  n = length(b);  m = length(t);
%
% Generate X.
%
  X = zeros(n,m); 
  for k=1:m
     if k==1;
        X(:,k) = expm(A*t(1))*b;
     else
	X(:,k) = expm(A*(t(k)-t(k-1)))*X(:,k-1);
     end	
  end
%
% Form the z-vector from X.
%
  Z = c'*X;
  z = zeros(m,1);
  z = Z(:);

 

|top|


A Finite Difference Approximation to phi

 function z = phi1(A,t,c,b,r)
%
% Assume that A is n-by-n and that c and b are column n-vectors.
% Assume that t is a length-m vector with 0 < t(1) < ... < t(m).
% Define X(:,0) = b and delk = t(k)-t(k-1) with t(0)=0;
% Computes a column m-vector z 
%
%                          c'*X(:,1)
%                    z  =     :
%                          c'*X(:,m)
%
% where X(:,k) = (I + delk*A + ... + (delk*A)^r/r!)*X(:,k-1).
%

  n = length(b);  m = length(t);
%
% Generate X.
%
  X = zeros(n,m); xprev = b; tprev = 0;
  for k=1:m
     delk = t(k) - tprev;
     X(:,k) = xprev;
     term = xprev;
     for j=1:r
        term = (delk/j)*(A*term);
        X(:,k) = X(:,k) + term;
     end
     tprev = t(k);
     xprev = X(:,k);
  end
%
% Form the z-vector from X.
%
  Z = c'*X;
  z = zeros(m,1);
  z = Z(:);

|top|


Another Finite Difference Approximation to phi
  function z = phi2(A,t,c,b,r)
%
% Assume that A is n-by-n and that c and b are column n-vectors.
% Assume that t is a length-m vector with 0 < t(1) < ... < t(m).
% Computes a column m-vector z 
%
%                          c'*X(:,1)
%                    z  =     :
%                          c'*X(:,m)
%
% where X(:,k) = (I + t(k)*A + ... + (t(k)*A)^r/r!)*b.
%

  n = length(b);  
  m = length(t);
%
% Generate X.
%
  X = zeros(n,m); 
  for k=1:m
     X(:,k) = b;
     term = b;
     for j=1:r
        term = (t(k)/j)*(A*term);
	X(:,k) = X(:,k) + term;
     end
  end
%
% Form the z-vector from X.
%
  Z = c'*X;
  z = zeros(m,1);
  z = Z(:);

|top|


Test B
% TestB.m
% Compares the Taylor approximations phi1 and phi2 to
% the "true" phi (obtained with "exact" exponentials)

% Set the problem size. (m >= n^2)
m = 10;
n = 3;

A = randn(n,n);
t = sort(rand(m,1));
b = randn(n,1);
c = randn(n,1);

clc
home
disp(sprintf('n = %2d',n))
disp(sprintf('m = %2d',m))
disp(sprintf('norm(A*t(m)) = %6.1f\n\n',norm(A*t(m))));
disp('  r       phi1      phi2        phi1        phi2')
disp('         flops     flops       error       error')
disp('---------------------------------------------------')
flops      norm(phi-phi2)
for r = 1:10
   z0 = phi(A,t,c,b);
   flops(0);
   z1 = phi1(A,t,c,b,r);
   f1 = flops;
   flops(0);
   z2 = phi2(A,t,c,b,r);
   f2 = flops;
   e1 = norm(z1-z0);
   e2 = norm(z2-z0);
   disp(sprintf('%3d  %8d  %8d     %10.3e  %10.3e',r,f1,f2,e1,e2))
end


Test B Output:
n =  3
m = 10
norm(A*t(m)) =    1.8


  r       phi1      phi2        phi1        phi2
         flops     flops       error       error
---------------------------------------------------
  1       320       310      5.029e-02   1.168e-01
  2       570       560      6.387e-03   5.799e-02
  3       820       810      5.351e-04   2.087e-02
  4      1070      1060      2.839e-05   4.368e-03
  5      1320      1310      1.048e-06   6.069e-04
  6      1570      1560      2.960e-08   6.231e-05
  7      1820      1810      7.210e-10   5.322e-06
  8      2070      2060      1.722e-11   4.369e-07
  9      2320      2310      4.106e-13   3.744e-08
 10      2570      2560      4.480e-15   3.177e-09

|top|


Finite Difference Jacobian of phi1
  function J = Jac1(Ac,t,c,b,r,zc)
%
% Finite difference Jacobian of z(A) = phi1(A,t,c,b,r) at A = Ac. 
% Ac is n-by-n, c and b are column n-vectors, t is an m-vector, 
% r is a positive integer, and zc = z(Ac).
%
% J is an m-by-n^2 matrix. If k = (j-1)n+i  with 1<=i<=n and 1<=j<=n, 
% then J(:,k) is an approximation of the partial of z(A) with
% respect to A(i,j) at A = Ac.

   m = length(t);
   n = length(c);
   
% There are more sophisticated ways to determine the dividied
% difference parameter, but the following will do for now:
   
   delta = sqrt(eps);  
   
   
   J = zeros(m,n^2);
   col = 0;
   for j=1:n
      for i=1:n
	 col = col+1;
	 Aij = Ac; 
         Aij(i,j) = Ac(i,j)+delta;
         J(:,col) = (phi1(Aij,t,c,b,r) - zc)/delta;
      end
   end


|top|


Finite Difference Jacobian of phi2

  function J = Jac2(Ac,t,c,b,r,zc)
%
% Finite difference Jacobian of z(A) = phi2(A,t,c,b,r) at A = Ac. 
% Ac is n-by-n, c and b are column n-vectors, t is an m-vector, 
% r is a positive integer, and zc = z(Ac).
%
% J is an m-by-n^2 matrix. If k = (j-1)n+i  with 1<=i<=n and 1<=j<=n, 
% then J(:,k) is an approximation of the partial of z(A) with
% respect to A(i,j) at A = Ac.

   m = length(t);
   n = length(c);
   
% There are more sophisticated ways to determine the divided
% difference parameter, but the following will do for now:
 
   delta = sqrt(eps);
   J = zeros(m,n^2);
   col = 0;
   for j=1:n
      for i=1:n
	 col = col+1;
	 Aij = Ac; 
	 Aij(i,j) = Ac(i,j)+delta;
	 J(:,col) = (phi2(Aij,t,c,b,r) - zc)/delta;
      end
   end

|top|


Fast Jacobian of phi2

  function J = FastJac2(A,t,c,b,r)
%
% The Jacobian of z(A) = phi2(A,t,c,b,r) at A. 
% A is n-by-n, c and b are column n-vectors, t is an m-vector, 
% r is a positive integer.
%
% J is an m-by-n^2 matrix. If k = (j-1)n+i  with 1<=i<=n and 1<=j<=n, 
% then J(:,k) is the partial of z(A) with
% respect to A(i,j).

   m = length(t);
   n = length(c);
   
   K = [eye(n,n) A];
   L = [eye(n,n) A'];
   for j=2:r-1
      P = A*K(:,(j-1)*n+1:j*n);
      K = [K P];
      L = [L P'];
   end
   c_vec = K'*c;
   b_vec = L'*b;
   
   
   J = zeros(m,n^2);
   [alfa,beta] = Pade(r,0);
   
   
  
   for k=1:m
      % Set up row k
      H = hankel(alfa.*(t(k).^((1:r)')));	
      col = 0;
      for j=1:n
	 z = H*b_vec(j:n:n*r);
         for i=1:n
	    col = col+1;
	    J(k,col) = (c_vec(i:n:n*r)'*z);
	 end
      end
   end

|top|


Test C
% TestC.m
% Jac2 and Jac2Fast

% Set the problem size. (m >= n^2)
n = 6;
m = n^2;

A = randn(n,n);
t = sort(rand(m,1));
b = randn(n,1);
c = randn(n,1);

clc
home
disp(sprintf('n = %2d',n))
disp(sprintf('m = %2d',m))
disp(sprintf('norm(A*t(m)) = %6.1f\n',norm(A*t(m))));
disp('J is computed via Jac2, Jfast is computed vis FastJac2, and')
disp('rho = max(max(abs(J-Jfast)))')
disp(' ')
disp('  r       Jac2    FastJac2      rho')
disp('         flops      flops       ')
disp('---------------------------------------------------')
flops      norm(phi-phi2)
for r = 1:n
   flops(0)
   zc = phi2(A,t,c,b,r);
   J =  Jac2(A,t,c,b,r,zc);
   f = flops;
   flops(0);
   Jfast = FastJac2(A,t,c,b,r);
   f1 = flops;
   err = max(max(abs(J-Jfast)));
   disp(sprintf('%3d   %8d  %8d     %10.3e  %11.3e',r,f,f1,err))
end


Test C Output


n =  6
m = 36
norm(A*t(m)) =    5.7

J is computed via Jac2, Jfast is computed vis FastJac2, and
rho = max(max(abs(J-Jfast)))
 
  r       Jac2    FastJac2      rho
         flops      flops       
---------------------------------------------------
  1     131870      4762      6.182e-08  
  2     245090     10314      1.851e-07  
  3     358310     15798      2.518e-07  
  4     471530     22146      3.225e-07  
  5     584750     29358      3.666e-07  
  6     697970     37434      3.511e-07  
 

|top|